KovalenkoQELCharmPXSec.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13 
14  @ Sep 13, 2007 - CA
15  Debugged the model in order to be included in the default event generation
16  threads in the next physics release (2.0.2). Rather than using Kovalenko's
17  expression for the ZR scaling factor, I apply an ad-hoc scaling factor
18  maintaining the relative strength of the QELC channels but lowering their
19  sum to be consistent with recent NOMAD measurement. The default value of
20  M0 has been changed from 0.1 to sqrt(0.1) as in M.Bischofberger's (ETHZ)
21  PhD thesis (DISS.ETH NO 16034). For more details see GENIE-PUB/2007/006.
22 */
23 //____________________________________________________________________________
24 
25 #include <TMath.h>
26 #include <Math/Integrator.h>
27 
37 /////#include "Numerical/IntegratorI.h"
45 
46 using namespace genie;
47 using namespace genie::constants;
48 
49 //____________________________________________________________________________
51 XSecAlgorithmI("genie::KovalenkoQELCharmPXSec")
52 {
53 
54 }
55 //____________________________________________________________________________
57 XSecAlgorithmI("genie::KovalenkoQELCharmPXSec", config)
58 {
59 
60 }
61 //____________________________________________________________________________
63 {
64 
65 }
66 //____________________________________________________________________________
68  const Interaction * interaction, KinePhaseSpace_t kps) const
69 {
70  if(! this -> ValidProcess (interaction) ) return 0.;
71  if(! this -> ValidKinematics (interaction) ) return 0.;
72 
73  //----- get kinematics & init state - compute auxiliary vars
74  const Kinematics & kinematics = interaction->Kine();
75  const InitialState & init_state = interaction->InitState();
76  const Target & target = init_state.Tgt();
77 
78  //neutrino energy & momentum transfer
79  double E = init_state.ProbeE(kRfHitNucRest);
80  double E2 = E * E;
81  double Q2 = kinematics.Q2();
82 
83 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
84  LOG("QELCharmXSec", pDEBUG) << "E = " << E << ", Q2 = " << Q2;
85 #endif
86 
87  //resonance mass & nucleon mass
88  double MR = this->MRes (interaction);
89  double MR2 = TMath::Power(MR,2);
90  double Mnuc = target.HitNucMass();
91  double Mnuc2 = TMath::Power(Mnuc,2);
92 
93 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
94  LOG("QELCharmXSec", pDEBUG) << "M(RES) = " << MR;
95 #endif
96 
97  //----- Calculate the differential cross section dxsec/dQ^2
98  double Gf = kGF2 / (2*kPi);
99  double vR = (MR2 - Mnuc2 + Q2) / (2*Mnuc);
100  double xiR = this->xiBar(Q2, Mnuc, vR);
101  double vR2 = vR*vR;
102  double vR_E = vR/E;
103  double Q2_4E2 = Q2/(4*E2);
104  double Q2_2MExiR = Q2/(2*Mnuc*E*xiR);
105  double Z = this->ZR(interaction);
106  double D = this->DR(interaction);
107 
108 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
109  LOG("QELCharmXSec", pDEBUG)
110  << "Z = " << Z << ", D = " << D << ". xiR = " << xiR << ", vR = " << vR;
111 #endif
112 
113  double xsec = Gf*Z*D * (1 - vR_E + Q2_4E2 + Q2_2MExiR) *
114  TMath::Sqrt(vR2 + Q2) / (vR*xiR);
115 
116  //----- The algorithm computes dxsec/dQ2
117  // Check whether variable tranformation is needed
118  if(kps!=kPSQ2fE) {
119  double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
120  xsec *= J;
121  }
122 
123  //----- If requested return the free nucleon xsec even for input nuclear tgt
124  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
125 
126  //----- Nuclear cross section (simple scaling here)
127  int nuc = target.HitNucPdg();
128  int NNucl = (pdg::IsProton(nuc)) ? target.Z() : target.N();
129  xsec *= NNucl;
130 
131 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
132  LOG("QELCharmXSec", pINFO)
133  << "dsigma/dQ2(E=" << E << ", Q2=" << Q2 << ") = "
134  << xsec / (1E-40*units::cm2) << " x 1E-40 cm^2";
135 #endif
136 
137  return xsec;
138 }
139 //____________________________________________________________________________
140 double KovalenkoQELCharmPXSec::ZR(const Interaction * interaction) const
141 {
142  const XclsTag & xcls = interaction->ExclTag();
143  const InitialState & init_state = interaction->InitState();
144 
145  int pdgc = xcls.CharmHadronPdg();
146  bool isP = pdg::IsProton ( init_state.Tgt().HitNucPdg() );
147  bool isN = pdg::IsNeutron( init_state.Tgt().HitNucPdg() );
148 
149  if ( pdgc == kPdgLambdaPc && isN ) return fScLambdaP;
150  else if ( pdgc == kPdgSigmaPc && isN ) return fScSigmaP;
151  else if ( pdgc == kPdgSigmaPPc && isP ) return fScSigmaPP;
152  else abort();
153 }
154 //____________________________________________________________________________
155 double KovalenkoQELCharmPXSec::DR(const Interaction * interaction) const
156 {
157  const InitialState & init_state = interaction -> InitState();
158 
159  // Compute PDFs
160  PDF pdfs;
161  pdfs.SetModel(fPDFModel); // <-- attach algorithm
162 
163  // Compute integration area = [xi_bar_plus, xi_bar_minus]
164  const Kinematics & kinematics = interaction->Kine();
165 
166  double Q2 = kinematics.Q2();
167  double Mnuc = init_state.Tgt().HitNucMass();
168  double Mnuc2 = TMath::Power(Mnuc,2);
169  double MR = this->MRes(interaction);
170  double DeltaR = this->ResDM(interaction);
171 
172  double vR_minus = ( TMath::Power(MR-DeltaR,2) - Mnuc2 + Q2 ) / (2*Mnuc);
173  double vR_plus = ( TMath::Power(MR+DeltaR,2) - Mnuc2 + Q2 ) / (2*Mnuc);
174 
175  LOG("QELCharmXSec", pDEBUG)
176  << "vR = [plus: " << vR_plus << ", minus: " << vR_minus << "]";
177 
178  double xi_bar_minus = 0.999;
179  double xi_bar_plus = this->xiBar(Q2, Mnuc, vR_plus);
180 
181  LOG("QELCharmXSec", pDEBUG)
182  << "Integration limits = [" << xi_bar_plus << ", " << xi_bar_minus << "]";
183 
184  int pdgc = init_state.Tgt().HitNucPdg();
185 
186  ROOT::Math::IBaseFunctionOneDim * integrand = new
190 
191  double abstol = 1; // We mostly care about relative tolerance
192  double reltol = 1E-4;
193  int nmaxeval = 100000;
194  ROOT::Math::Integrator ig(*integrand,ig_type,abstol,reltol,nmaxeval);
195  double D = ig.Integral(xi_bar_plus, xi_bar_minus);
196 
197  delete integrand;
198 
199  return D;
200 }
201 //____________________________________________________________________________
202 double KovalenkoQELCharmPXSec::xiBar(double Q2, double Mnuc, double v) const
203 {
204  double Mo2 = fMo*fMo;
205  double v2 = v *v;
206  double xi = (Q2/Mnuc) / (v + TMath::Sqrt(v2+Q2));
207  double xib = xi * ( 1 + (1 + Mo2/(Q2+Mo2))*Mo2/Q2 );
208  return xib;
209 }
210 //____________________________________________________________________________
211 double KovalenkoQELCharmPXSec::ResDM(const Interaction * interaction) const
212 {
213 // Resonance Delta M obeys the constraint DM(R+/-) <= |M(R+/-) - M(R)|
214 // where M(R-) <= M(R) <= M(R+) are the masses of the neighboring
215 // resonances R+, R-.
216 // Get the values from the algorithm conf. registry, and if they do not exist
217 // set them to default values (Eq.(20) in Sov.J.Nucl.Phys.52:934 (1990)
218 
219  const XclsTag & xcls = interaction->ExclTag();
220 
221  int pdgc = xcls.CharmHadronPdg();
222 
223  bool isLambda = (pdgc == kPdgLambdaPc);
224  bool isSigma = (pdgc == kPdgSigmaPc || pdgc == kPdgSigmaPPc);
225 
226  if ( isLambda ) return fResDMLambda;
227  else if ( isSigma ) return fResDMSigma;
228  else abort();
229 
230  return 0;
231 }
232 //____________________________________________________________________________
233 double KovalenkoQELCharmPXSec::MRes(const Interaction * interaction) const
234 {
235  const XclsTag & xcls = interaction->ExclTag();
236 
237  int pdgc = xcls.CharmHadronPdg();
238  double MR = PDGLibrary::Instance()->Find(pdgc)->Mass();
239  return MR;
240 }
241 //____________________________________________________________________________
242 double KovalenkoQELCharmPXSec::Integral(const Interaction * interaction) const
243 {
244  double xsec = fXSecIntegrator->Integrate(this,interaction);
245  return xsec;
246 }
247 //____________________________________________________________________________
249  const Interaction * interaction) const
250 {
251  // Make sure we are dealing with one of the following channels:
252  // v + n --> mu- + Lambda_{c}^{+} (2285)
253  // v + n --> mu- + Sigma_{c}^{+} (2455)
254  // v + p --> mu- + Sigma_{c}^{++} (2455)
255 
256  if(interaction->TestBit(kISkipProcessChk)) return true;
257 
258  const XclsTag & xcls = interaction->ExclTag();
259  const InitialState & init_state = interaction->InitState();
260  const ProcessInfo & proc_info = interaction->ProcInfo();
261 
262  bool is_exclusive_charm = (xcls.IsCharmEvent() && !xcls.IsInclusiveCharm());
263  if(!is_exclusive_charm) return false;
264 
265  if(!proc_info.IsQuasiElastic()) return false;
266  if(!proc_info.IsWeak()) return false;
267 
268  bool isP = pdg::IsProton ( init_state.Tgt().HitNucPdg() );
269  bool isN = pdg::IsNeutron( init_state.Tgt().HitNucPdg() );
270 
271  int pdgc = xcls.CharmHadronPdg();
272 
273  bool can_handle = (
274  (pdgc == kPdgLambdaPc && isN) || /* v + n -> l + #Lambda_{c}^{+} */
275  (pdgc == kPdgSigmaPc && isN) || /* v + n -> l + #Sigma_{c}^{+} */
276  (pdgc == kPdgSigmaPPc && isP) /* v + p -> l + #Sigma_{c}^{++} */
277  );
278  return can_handle;
279 }
280 //____________________________________________________________________________
282  const Interaction * interaction) const
283 {
284  if(interaction->TestBit(kISkipKinematicChk)) return true;
285 
286  const InitialState & init_state = interaction->InitState();
287  double E = init_state.ProbeE(kRfHitNucRest);
288 
289  //resonance, final state primary lepton & nucleon mass
290  double MR = this -> MRes (interaction);
291  double ml = interaction->FSPrimLepton()->Mass();
292  double Mnuc = init_state.Tgt().HitNucP4Ptr()->M();
293  double Mnuc2 = TMath::Power(Mnuc,2);
294 
295  //resonance threshold
296  double ER = ( TMath::Power(MR+ml,2) - Mnuc2 ) / (2*Mnuc);
297 
298  if(E <= ER) return false;
299 
300  return true;
301 }
302 //____________________________________________________________________________
304 {
305  Algorithm::Configure(config);
306  this->LoadConfig();
307 }
308 //____________________________________________________________________________
309 void KovalenkoQELCharmPXSec::Configure(string param_set)
310 {
311  Algorithm::Configure(param_set);
312  this->LoadConfig();
313 }
314 //____________________________________________________________________________
316 {
317  fPDFModel = 0;
318  ///fIntegrator = 0;
319 
320  // Get config values or set defaults
321  GetParamDef( "Scale-LambdaP", fScLambdaP, 0.8 * 0.0102 ) ;
322  GetParamDef( "Scale-SigmaP", fScSigmaP , 0.8 * 0.0028 ) ;
323  GetParamDef( "Scale-SigmaPP", fScSigmaPP, 0.8 * 0.0159 ) ;
324  GetParamDef( "Res-DeltaM-Lambda", fResDMLambda, 0.56 ) ; //GeV
325  GetParamDef( "Res-DeltaM-Sigma", fResDMSigma, 0.20 ) ; //GeV
326  GetParamDef( "Mo", fMo, sqrt(0.1) ); //GeV
327 
328  // get PDF model and integrator
329 
330  fPDFModel = dynamic_cast<const PDFModelI *>(this->SubAlg("PDF-Set"));
331  assert(fPDFModel);
332 
333  // load XSec Integrator
335  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
337 
338  // load numerical integrator for integrand in diff x-section calc.
339 // fIntegrator = dynamic_cast<const IntegratorI *>(this->SubAlg("Integrator"));
340 // assert(fIntegrator);
341 }
342 //____________________________________________________________________________
343 // Auxiliary scalar function for internal integration
344 //____________________________________________________________________________
346  PDF * pdf, double Q2, int nucleon_pdgc) :
347 ROOT::Math::IBaseFunctionOneDim()
348 {
349  fPDF = pdf;
350  fQ2 = TMath::Max(0.3, Q2);
351  fPdgC = nucleon_pdgc;
352 }
353 //____________________________________________________________________________
355 {
356 
357 }
358 //____________________________________________________________________________
360 {
361  return 1;
362 }
363 //____________________________________________________________________________
365 {
366  if(xin<0 || xin>1) return 0.;
367 
368  fPDF->Calculate(xin, fQ2);
369  bool isP = pdg::IsProton(fPdgC);
370  double f = (isP) ? fPDF->DownValence() : fPDF->UpValence();
371 
372  LOG("QELCharmXSec", pDEBUG)
373  << "f(xin = " << xin << ", Q2 = " << fQ2 << ") = " << f;
374 
375  return f;
376 }
377 //____________________________________________________________________________
378 ROOT::Math::IBaseFunctionOneDim *
380 {
382 }
383 //____________________________________________________________________________
Cross Section Calculation Interface.
const double kPi
Basic constants.
bool IsWeak(void) const
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition: GSLUtils.cxx:31
const XML_Char * target
Definition: expat.h:268
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
int HitNucPdg(void) const
Definition: Target.cxx:321
void Configure(const Registry &config)
double DownValence(void) const
Definition: PDF.h:52
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
double HitNucMass(void) const
Definition: Target.cxx:250
A class to store PDFs.
Definition: PDF.h:38
int CharmHadronPdg(void) const
Definition: XclsTag.h:50
T sqrt(T number)
Definition: d0nt_math.hpp:156
double MRes(const Interaction *interaction) const
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double ZR(const Interaction *interaction) const
Definition: config.py:1
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:29
static const double cm2
Definition: Units.h:77
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:37
const int kPdgSigmaPPc
Definition: PDGCodes.h:86
bool IsCharmEvent(void) const
Definition: XclsTag.h:48
Float_t Z
Definition: plot.C:38
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
Summary information for an interaction.
Definition: Interaction.h:56
double DR(const Interaction *interaction) const
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Float_t E
Definition: plot.C:20
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
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:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
const int kPdgLambdaPc
Definition: PDGCodes.h:83
void SetModel(const PDFModelI *model)
Definition: PDF.cxx:48
double ResDM(const Interaction *interaction) const
int Z(void) const
Definition: Target.h:69
#define pINFO
Definition: Messenger.h:63
Auxiliary scalar function for the internal integration in Kovalenko QEL charm production cross sectio...
void Calculate(double x, double q2)
Definition: PDF.cxx:55
Double_t xsec[nknots]
Definition: testXsec.C:47
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:70
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
double UpValence(void) const
Definition: PDF.h:51
const int kPdgSigmaPc
Definition: PDGCodes.h:85
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
double integrand(double *x, double *par)
bool IsInclusiveCharm(void) const
Definition: XclsTag.cxx:63
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
assert(nhit_max >=nhit_nbins)
double Integral(const Interaction *i) const
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
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
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
Definition: InitialState.h:67
static const double kGF2
Definition: Constants.h:60
void kinematics()
Definition: kinematics.C:10
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:49
#define pDEBUG
Definition: Messenger.h:64
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353