GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
KovalenkoQELCharmPXSec.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 #include <Math/Integrator.h>
13 
17 #include "Framework/Conventions/GBuild.h"
23 /////#include "Numerical/IntegratorI.h"
31 
32 using namespace genie;
33 using namespace genie::constants;
34 
35 //____________________________________________________________________________
37 XSecAlgorithmI("genie::KovalenkoQELCharmPXSec")
38 {
39 
40 }
41 //____________________________________________________________________________
43 XSecAlgorithmI("genie::KovalenkoQELCharmPXSec", config)
44 {
45 
46 }
47 //____________________________________________________________________________
49 {
50 
51 }
52 //____________________________________________________________________________
54  const Interaction * interaction, KinePhaseSpace_t kps) const
55 {
56  if(! this -> ValidProcess (interaction) ) return 0.;
57  if(! this -> ValidKinematics (interaction) ) return 0.;
58 
59  //----- get kinematics & init state - compute auxiliary vars
60  const Kinematics & kinematics = interaction->Kine();
61  const InitialState & init_state = interaction->InitState();
62  const Target & target = init_state.Tgt();
63 
64  //neutrino energy & momentum transfer
65  double E = init_state.ProbeE(kRfHitNucRest);
66  double E2 = E * E;
67  double Q2 = kinematics.Q2();
68 
69 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
70  LOG("QELCharmXSec", pDEBUG) << "E = " << E << ", Q2 = " << Q2;
71 #endif
72 
73  //resonance mass & nucleon mass
74  double MR = this->MRes (interaction);
75  double MR2 = TMath::Power(MR,2);
76  double Mnuc = target.HitNucMass();
77  double Mnuc2 = TMath::Power(Mnuc,2);
78 
79 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
80  LOG("QELCharmXSec", pDEBUG) << "M(RES) = " << MR;
81 #endif
82 
83  //----- Calculate the differential cross section dxsec/dQ^2
84  double Gf = kGF2 / (2*kPi);
85  double vR = (MR2 - Mnuc2 + Q2) / (2*Mnuc);
86  double xiR = this->xiBar(Q2, Mnuc, vR);
87  double vR2 = vR*vR;
88  double vR_E = vR/E;
89  double Q2_4E2 = Q2/(4*E2);
90  double Q2_2MExiR = Q2/(2*Mnuc*E*xiR);
91  double Z = this->ZR(interaction);
92  double D = this->DR(interaction);
93 
94 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
95  LOG("QELCharmXSec", pDEBUG)
96  << "Z = " << Z << ", D = " << D << ". xiR = " << xiR << ", vR = " << vR;
97 #endif
98 
99  double xsec = Gf*Z*D * (1 - vR_E + Q2_4E2 + Q2_2MExiR) *
100  TMath::Sqrt(vR2 + Q2) / (vR*xiR);
101 
102  //----- The algorithm computes dxsec/dQ2
103  // Check whether variable tranformation is needed
104  if(kps!=kPSQ2fE) {
105  double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
106  xsec *= J;
107  }
108 
109  //----- If requested return the free nucleon xsec even for input nuclear tgt
110  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
111 
112  //----- Nuclear cross section (simple scaling here)
113  int nuc = target.HitNucPdg();
114  int NNucl = (pdg::IsProton(nuc)) ? target.Z() : target.N();
115  xsec *= NNucl;
116 
117 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
118  LOG("QELCharmXSec", pINFO)
119  << "dsigma/dQ2(E=" << E << ", Q2=" << Q2 << ") = "
120  << xsec / (1E-40*units::cm2) << " x 1E-40 cm^2";
121 #endif
122 
123  return xsec;
124 }
125 //____________________________________________________________________________
126 double KovalenkoQELCharmPXSec::ZR(const Interaction * interaction) const
127 {
128  const XclsTag & xcls = interaction->ExclTag();
129  const InitialState & init_state = interaction->InitState();
130 
131  int pdgc = xcls.CharmHadronPdg();
132  bool isP = pdg::IsProton ( init_state.Tgt().HitNucPdg() );
133  bool isN = pdg::IsNeutron( init_state.Tgt().HitNucPdg() );
134 
135  if ( pdgc == kPdgLambdaPc && isN ) return fScLambdaP;
136  else if ( pdgc == kPdgSigmaPc && isN ) return fScSigmaP;
137  else if ( pdgc == kPdgSigmaPPc && isP ) return fScSigmaPP;
138  else abort();
139 }
140 //____________________________________________________________________________
141 double KovalenkoQELCharmPXSec::DR(const Interaction * interaction) const
142 {
143  const InitialState & init_state = interaction -> InitState();
144 
145  // Compute PDFs
146  PDF pdfs;
147  pdfs.SetModel(fPDFModel); // <-- attach algorithm
148 
149  // Compute integration area = [xi_bar_plus, xi_bar_minus]
150  const Kinematics & kinematics = interaction->Kine();
151 
152  double Q2 = kinematics.Q2();
153  double Mnuc = init_state.Tgt().HitNucMass();
154  double Mnuc2 = TMath::Power(Mnuc,2);
155  double MR = this->MRes(interaction);
156  double DeltaR = this->ResDM(interaction);
157 
158  double vR_minus = ( TMath::Power(MR-DeltaR,2) - Mnuc2 + Q2 ) / (2*Mnuc);
159  double vR_plus = ( TMath::Power(MR+DeltaR,2) - Mnuc2 + Q2 ) / (2*Mnuc);
160 
161  LOG("QELCharmXSec", pDEBUG)
162  << "vR = [plus: " << vR_plus << ", minus: " << vR_minus << "]";
163 
164  double xi_bar_minus = 0.999;
165  double xi_bar_plus = this->xiBar(Q2, Mnuc, vR_plus);
166 
167  LOG("QELCharmXSec", pDEBUG)
168  << "Integration limits = [" << xi_bar_plus << ", " << xi_bar_minus << "]";
169 
170  int pdgc = init_state.Tgt().HitNucPdg();
171 
172  ROOT::Math::IBaseFunctionOneDim * integrand = new
174  ROOT::Math::IntegrationOneDim::Type ig_type =
176 
177  double abstol = 1; // We mostly care about relative tolerance
178  double reltol = 1E-4;
179  int nmaxeval = 100000;
180  ROOT::Math::Integrator ig(*integrand,ig_type,abstol,reltol,nmaxeval);
181  double D = ig.Integral(xi_bar_plus, xi_bar_minus);
182 
183  delete integrand;
184 
185  return D;
186 }
187 //____________________________________________________________________________
188 double KovalenkoQELCharmPXSec::xiBar(double Q2, double Mnuc, double v) const
189 {
190  double Mo2 = fMo*fMo;
191  double v2 = v *v;
192  double xi = (Q2/Mnuc) / (v + TMath::Sqrt(v2+Q2));
193  double xib = xi * ( 1 + (1 + Mo2/(Q2+Mo2))*Mo2/Q2 );
194  return xib;
195 }
196 //____________________________________________________________________________
197 double KovalenkoQELCharmPXSec::ResDM(const Interaction * interaction) const
198 {
199 // Resonance Delta M obeys the constraint DM(R+/-) <= |M(R+/-) - M(R)|
200 // where M(R-) <= M(R) <= M(R+) are the masses of the neighboring
201 // resonances R+, R-.
202 // Get the values from the algorithm conf. registry, and if they do not exist
203 // set them to default values (Eq.(20) in Sov.J.Nucl.Phys.52:934 (1990)
204 
205  const XclsTag & xcls = interaction->ExclTag();
206 
207  int pdgc = xcls.CharmHadronPdg();
208 
209  bool isLambda = (pdgc == kPdgLambdaPc);
210  bool isSigma = (pdgc == kPdgSigmaPc || pdgc == kPdgSigmaPPc);
211 
212  if ( isLambda ) return fResDMLambda;
213  else if ( isSigma ) return fResDMSigma;
214  else abort();
215 
216  return 0;
217 }
218 //____________________________________________________________________________
219 double KovalenkoQELCharmPXSec::MRes(const Interaction * interaction) const
220 {
221  const XclsTag & xcls = interaction->ExclTag();
222 
223  int pdgc = xcls.CharmHadronPdg();
224  double MR = PDGLibrary::Instance()->Find(pdgc)->Mass();
225  return MR;
226 }
227 //____________________________________________________________________________
228 double KovalenkoQELCharmPXSec::Integral(const Interaction * interaction) const
229 {
230  double xsec = fXSecIntegrator->Integrate(this,interaction);
231  return xsec;
232 }
233 //____________________________________________________________________________
235  const Interaction * interaction) const
236 {
237  // Make sure we are dealing with one of the following channels:
238  // v + n --> mu- + Lambda_{c}^{+} (2285)
239  // v + n --> mu- + Sigma_{c}^{+} (2455)
240  // v + p --> mu- + Sigma_{c}^{++} (2455)
241 
242  if(interaction->TestBit(kISkipProcessChk)) return true;
243 
244  const XclsTag & xcls = interaction->ExclTag();
245  const InitialState & init_state = interaction->InitState();
246  const ProcessInfo & proc_info = interaction->ProcInfo();
247 
248  bool is_exclusive_charm = (xcls.IsCharmEvent() && !xcls.IsInclusiveCharm());
249  if(!is_exclusive_charm) return false;
250 
251  if(!proc_info.IsQuasiElastic()) return false;
252  if(!proc_info.IsWeak()) return false;
253 
254  bool isP = pdg::IsProton ( init_state.Tgt().HitNucPdg() );
255  bool isN = pdg::IsNeutron( init_state.Tgt().HitNucPdg() );
256 
257  int pdgc = xcls.CharmHadronPdg();
258 
259  bool can_handle = (
260  (pdgc == kPdgLambdaPc && isN) || /* v + n -> l + #Lambda_{c}^{+} */
261  (pdgc == kPdgSigmaPc && isN) || /* v + n -> l + #Sigma_{c}^{+} */
262  (pdgc == kPdgSigmaPPc && isP) /* v + p -> l + #Sigma_{c}^{++} */
263  );
264  return can_handle;
265 }
266 //____________________________________________________________________________
268  const Interaction * interaction) const
269 {
270  if(interaction->TestBit(kISkipKinematicChk)) return true;
271 
272  const InitialState & init_state = interaction->InitState();
273  double E = init_state.ProbeE(kRfHitNucRest);
274 
275  //resonance, final state primary lepton & nucleon mass
276  double MR = this -> MRes (interaction);
277  double ml = interaction->FSPrimLepton()->Mass();
278  double Mnuc = init_state.Tgt().HitNucP4Ptr()->M();
279  double Mnuc2 = TMath::Power(Mnuc,2);
280 
281  //resonance threshold
282  double ER = ( TMath::Power(MR+ml,2) - Mnuc2 ) / (2*Mnuc);
283 
284  if(E <= ER) return false;
285 
286  return true;
287 }
288 //____________________________________________________________________________
290 {
291  Algorithm::Configure(config);
292  this->LoadConfig();
293 }
294 //____________________________________________________________________________
295 void KovalenkoQELCharmPXSec::Configure(string param_set)
296 {
297  Algorithm::Configure(param_set);
298  this->LoadConfig();
299 }
300 //____________________________________________________________________________
302 {
303  fPDFModel = 0;
304  ///fIntegrator = 0;
305 
306  // Get config values or set defaults
307  GetParamDef( "Scale-LambdaP", fScLambdaP, 0.8 * 0.0102 ) ;
308  GetParamDef( "Scale-SigmaP", fScSigmaP , 0.8 * 0.0028 ) ;
309  GetParamDef( "Scale-SigmaPP", fScSigmaPP, 0.8 * 0.0159 ) ;
310  GetParamDef( "Res-DeltaM-Lambda", fResDMLambda, 0.56 ) ; //GeV
311  GetParamDef( "Res-DeltaM-Sigma", fResDMSigma, 0.20 ) ; //GeV
312  GetParamDef( "Mo", fMo, sqrt(0.1) ); //GeV
313 
314  // get PDF model and integrator
315 
316  fPDFModel = dynamic_cast<const PDFModelI *>(this->SubAlg("PDF-Set"));
317  assert(fPDFModel);
318 
319  // load XSec Integrator
321  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
322  assert(fXSecIntegrator);
323 
324  // load numerical integrator for integrand in diff x-section calc.
325 // fIntegrator = dynamic_cast<const IntegratorI *>(this->SubAlg("Integrator"));
326 // assert(fIntegrator);
327 }
328 //____________________________________________________________________________
329 // Auxiliary scalar function for internal integration
330 //____________________________________________________________________________
332  PDF * pdf, double Q2, int nucleon_pdgc) :
333 ROOT::Math::IBaseFunctionOneDim()
334 {
335  fPDF = pdf;
336  fQ2 = TMath::Max(0.3, Q2);
337  fPdgC = nucleon_pdgc;
338 }
339 //____________________________________________________________________________
341 {
342 
343 }
344 //____________________________________________________________________________
346 {
347  return 1;
348 }
349 //____________________________________________________________________________
351 {
352  if(xin<0 || xin>1) return 0.;
353 
354  fPDF->Calculate(xin, fQ2);
355  bool isP = pdg::IsProton(fPdgC);
356  double f = (isP) ? fPDF->DownValence() : fPDF->UpValence();
357 
358  LOG("QELCharmXSec", pDEBUG)
359  << "f(xin = " << xin << ", Q2 = " << fQ2 << ") = " << f;
360 
361  return f;
362 }
363 //____________________________________________________________________________
364 ROOT::Math::IBaseFunctionOneDim *
366 {
367  return new utils::gsl::wrap::KovQELCharmIntegrand(fPDF, fQ2, fPdgC);
368 }
369 //____________________________________________________________________________
Cross Section Calculation Interface.
bool IsWeak(void) const
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition: GSLUtils.cxx:23
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
int HitNucPdg(void) const
Definition: Target.cxx:304
void Configure(const Registry &config)
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
double HitNucMass(void) const
Definition: Target.cxx:233
A class to store PDFs.
Definition: PDF.h:37
int CharmHadronPdg(void) const
Definition: XclsTag.h:52
double MRes(const Interaction *interaction) const
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double ZR(const Interaction *interaction) const
bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
const XSecIntegratorI * fXSecIntegrator
const IntegratorI * fIntegrator;
enum genie::EKinePhaseSpace KinePhaseSpace_t
Pure abstract base class. Defines the PDFModelI interface to be implemented by wrapper classes to exi...
Definition: PDFModelI.h:28
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
const int kPdgSigmaPPc
Definition: PDGCodes.h:102
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
Summary information for an interaction.
Definition: Interaction.h:56
double DR(const Interaction *interaction) const
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double cm2
Definition: Units.h:69
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
const Kinematics & Kine(void) const
Definition: Interaction.h:71
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
const int kPdgLambdaPc
Definition: PDGCodes.h:99
void SetModel(const PDFModelI *model)
Definition: PDF.cxx:42
double ResDM(const Interaction *interaction) const
int Z(void) const
Definition: Target.h:68
#define pINFO
Definition: Messenger.h:62
Auxiliary scalar function for the internal integration in Kovalenko QEL charm production cross sectio...
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
int N(void) const
Definition: Target.h:69
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
const int kPdgSigmaPc
Definition: PDGCodes.h:101
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
bool IsInclusiveCharm(void) const
Definition: XclsTag.cxx:54
double Integral(const Interaction *i) const
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
double xiBar(double Q2, double Mnuc, double v) const
const InitialState & InitState(void) const
Definition: Interaction.h:69
KovQELCharmIntegrand(PDF *pdf, double Q2, int nucleon_pdgc)
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
Definition: InitialState.h:66
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345