EmpiricalMECPXSec2015.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  Steve Dytman <dytman+ \at pitt.edu>
11  Pittsburgh University
12 
13  For the class documentation see the corresponding header file.
14 
15  Important revisions after version 2.0.0 :
16  @ Nov 28, 2011 - CA
17  Integrated cross section for CCMEC is taken to be a fraction of the
18  CCQE cross section for the given neutrino energy and nucleus.
19  @ Dec 9, 2011 - SD
20  Using a simple model now - 2N mass chosen with a Gaussian.
21  Strength is tuned to get agreement with MiniBoone and NOMAD tot xs.
22  Parameters tuned to get best agreement with MiniBoone muon angle/energy dist.
23  @ May 14, 2012 - SD
24  Changed empirical model to have correct Q2 behavior. Shape is
25  now Q2*Gaussian*dipole form factor.
26  Add (e,e') MEC, start development to make them dependent on same model.
27 */
28 //____________________________________________________________________________
29 
30 #include <TMath.h>
31 
44 
45 using namespace genie;
46 using namespace genie::constants;
47 using namespace genie::controls;
48 
49 //____________________________________________________________________________
51 XSecAlgorithmI("genie::EmpiricalMECPXSec2015")
52 {
53 
54 }
55 //____________________________________________________________________________
57 XSecAlgorithmI("genie::EmpiricalMECPXSec2015", config)
58 {
59 
60 }
61 //____________________________________________________________________________
63 {
64 
65 }
66 //____________________________________________________________________________
68  const Interaction * interaction, KinePhaseSpace_t kps) const
69 {
70 
71 // meson exchange current contribution depends a lot on QE model.
72 // This is an empirical model in development, not used in default event generation.
73 
74 // if(! this -> ValidProcess (interaction) ) return 0.;
75 // if(! this -> ValidKinematics (interaction) ) return 0.;
76  bool iscc = interaction->ProcInfo().IsWeakCC();
77  bool isnc = interaction->ProcInfo().IsWeakNC();
78  bool isem = interaction->ProcInfo().IsEM();
79 
80  const Kinematics & kinematics = interaction -> Kine();
81  double W = kinematics.W();
82  double Q2 = kinematics.Q2();
83  //LOG("MEC", pINFO) << "W, Q2 trial= " << W << " " << Q2 ;
84 
85  //
86  // Do a check whether W,Q2 is allowed. Return 0 otherwise.
87  //
88  double Ev = interaction->InitState().ProbeE(kRfHitNucRest); // kRfLab
89  int nucleon_cluster_pdg = interaction->InitState().Tgt().HitNucPdg();
90  double M2n = PDGLibrary::Instance()->Find(nucleon_cluster_pdg)-> Mass(); // nucleon cluster mass
91  double M2n2 = M2n*M2n;
92  double ml = interaction->FSPrimLepton()->Mass();
94 
95  //LOG("MEC", pINFO) << "Ev, ml, M2n = " << Ev << " " << ml << " " << M2n;
96  //LOG("MEC", pINFO) << "Wlim= " << Wlim.min << " " <<Wlim.max ;
97  if(W < Wlim.min || W > Wlim.max)
98  {double xsec = 0.;
99  return xsec;
100  }
101  //use proper Q2 limit from Controls.h
103 
104  //LOG("MEC", pINFO) << "Q2lim= " << Q2lim.min << " " <<Q2lim.max ;
105  if(Q2 < Q2lim.min || Q2 > Q2lim.max)
106  {double xsec = 0.;
107  return xsec;
108  }
109 
110  //get x and y
111  double x = 0.;
112  double y = 0.;
113  genie::utils::kinematics::WQ2toXY(Ev,M2n,W,Q2,x,y);
114  // LOG("MEC", pINFO) << "x = " << x << ", y = " << y;
115  // double Tmu = (1.-y)*Ev; // UNUSED - comment to quiet compiler warnings
116 
117  // Calculate d^2xsec/dWdQ2 - first form factor which is common to both
118  double Wdep = TMath::Gaus(W, fMass, fWidth);
119  double Q2dep = Q2*TMath::Power((1+Q2/fMq2d),-8.);
120  // double nudep = TMath::Power(Tmu,2.5);
121  // LOG("MEC", pINFO) << "Tmu = " << Tmu << ", nudep = " << nudep;
122  double FF2 = Wdep * Q2dep;// * nudep;
123  //LOG("MEC", pINFO) << "form factor = " << FF2 << ", Q2 = " << Q2 << ", W = " << W;
124 
125 // using formulas in Bodek and Budd for (e,e') inclusive cross section
126  double xsec = 1.;
127  if(isem) {
128  // Calculate scattering angle
129  //
130  // Q^2 = 4 * E^2 * sin^2 (theta/2) / ( 1 + 2 * (E/M) * sin^2(theta/2) ) =>
131  // sin^2 (theta/2) = MQ^2 / (4ME^2 - 2EQ^2)
132 
133  double E = Ev;
134  double E2 = E*E;
135  double sin2_halftheta = M2n*Q2 / (4*M2n*E2 - 2*E*Q2);
136  // double sin4_halftheta = TMath::Power(sin2_halftheta, 2.);
137  double cos2_halftheta = 1.-sin2_halftheta;
138  // double cos_halftheta = TMath::Sqrt(cos2_halftheta);
139  double tan2_halftheta = sin2_halftheta/cos2_halftheta;
140  double Q4 = Q2*Q2;
141 
142  // Calculate tau and the virtual photon polarization (epsilon)
143  double tau = Q2/(4*M2n2);
144  // double epsilon = 1. / (1. + 2.*(tau/x))*tan2_halftheta); //different than RosenbluthPXSec.cxx
145 
146  // Calculate the scattered lepton energy
147  double Ep = E / (1. + 2.*(E/M2n)*sin2_halftheta);
148  double Ep2 = Ep*Ep;
149 
150  //calculate cross section - d2sig/dOmega dE for purely transverse process
151  xsec = 4*kAem2*Ep2*cos2_halftheta/Q4 * FF2 * (tau/(1+tau) +2*tau*tan2_halftheta);
152  }
153  // use BB formula which seems to be same as Llewlyn-Smith
154  // note B term is only difference between nu and antinu, so both same here
155  else if(isnc||iscc){
156  double tau = Q2/(4*M2n2);
157  double tau2 = tau*tau;
158  double smufac = 4*M2n*Ev - Q2 - ml*ml;
159  double A = (ml*ml+Q2)/M2n2 * (tau*(1+tau) - tau2*(1-tau)+4*tau2)/TMath::Power(1+tau,2.) * FF2;
160  double C = tau/4/(1+tau) * FF2;
161  xsec = A + smufac*smufac*C; // CC or NC case - Llewelyn-Smith for transverse vector process.
162  }
163  // Check whether variable tranformation is needed
164  if(kps!=kPSWQ2fE) {
165  double J = utils::kinematics::Jacobian(interaction,kPSWQ2fE,kps);
166 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
167  LOG("MEC", pDEBUG)
168  << "Jacobian for transformation to: "
169  << KinePhaseSpace::AsString(kps) << ", J = " << J;
170 #endif
171  xsec *= J;
172  }
173 
174  return xsec;
175 }
176 //____________________________________________________________________________
177 double EmpiricalMECPXSec2015::Integral(const Interaction * interaction) const
178 {
179 // Calculate the CCMEC cross section as a fraction of the CCQE cross section
180 // for the given nuclear target at the given energy.
181 // Alternative strategy is to calculate the MEC cross section as the difference
182 // of CCQE cross section for two different M_A values (eg ~1.3 GeV and ~1.0 GeV)
183 // Include hit-object combinatorial factor? Would yield different A-dependence
184 // for MEC and QE.
185 //
186 
187  bool iscc = interaction->ProcInfo().IsWeakCC();
188  bool isnc = interaction->ProcInfo().IsWeakNC();
189  bool isem = interaction->ProcInfo().IsEM();
190 
191  int nupdg = interaction->InitState().ProbePdg();
192  int tgtpdg = interaction->InitState().Tgt().Pdg();
193  double E = interaction->InitState().ProbeE(kRfLab);
194  int nucleon_cluster_pdg = interaction->InitState().Tgt().HitNucPdg();
195  double Z=interaction->InitState().Tgt().Z();
196  double A=interaction->InitState().Tgt().A();
197  double N=A-Z;
198 
199  if(iscc) {
200 
201  int nucpdg = 0;
202  // neutrino CC: calculate the CCQE cross section resetting the
203  // hit nucleon cluster to neutron
204  if(pdg::IsNeutrino(nupdg)) {
205  nucpdg = kPdgNeutron;
206  }
207  // anti-neutrino CC: calculate the CCQE cross section resetting the
208  // hit nucleon cluster to proton
209  else
210  if(pdg::IsAntiNeutrino(nupdg)) {
211  nucpdg = kPdgProton;
212  }
213  else {
214  exit(1);
215  }
216 
217  // Create a tmp QE process
218  Interaction * in = Interaction::QELCC(tgtpdg,nucpdg,nupdg,E);
219 
220  // Calculate cross section for the QE process
221  double xsec = fXSecAlgCCQE->Integral(in);
222 
223  // Add A dependence which is not known from theory
224  double fFracADep = 1.;
225  if(A>=12) fFracADep = TMath::Power((N/6.),fMECAPower-1.);
226 
227  // Use tunable fraction
228  // FFracCCQE is fraction of QE going to MEC
229  // fFracCCQE_cluster is fraction of MEC going to each NN pair
230  // double fFracCCQE = fFracCCQElo;
231 
232  double fFracCCQE_cluster=0.;
233  if(pdg::IsNeutrino(nupdg) && nucleon_cluster_pdg==2000000201) fFracCCQE_cluster= fFracPN_CC; //n+p
234  if(pdg::IsNeutrino(nupdg) && nucleon_cluster_pdg==2000000200) fFracCCQE_cluster= 1.0-fFracPN_CC; //n+n
235  if(pdg::IsAntiNeutrino(nupdg) && nucleon_cluster_pdg==2000000201) fFracCCQE_cluster= fFracPN_CC; //n+p
236  if(pdg::IsAntiNeutrino(nupdg) && nucleon_cluster_pdg==2000000202) fFracCCQE_cluster= 1.0-fFracPN_CC; //p+p
237 
238 
239  xsec *= fFracCCQE*fFracCCQE_cluster*fFracADep;
240 
241  // Use gross combinatorial factor (number of 2-nucleon targets over number
242  // of 1-nucleon targets) : (A-1)/2
243  // double combfact = (in->InitState().Tgt().A()-1)/2.;
244  // xsec *= combfact;
245 
246  delete in;
247  return xsec;
248  }
249 
250  else if(isnc) {
251  int nucpdg = kPdgProton;
252  // Create a tmp QE process
253  Interaction * inp = Interaction::QELNC(tgtpdg,nucpdg,nupdg,E);
254  nucpdg = kPdgNeutron;
255  // Create a tmp QE process
256  Interaction * inn = Interaction::QELNC(tgtpdg,nucpdg,nupdg,E);
257 
258  // Calculate cross section for the QE process - avg of p and n - best for isoscalar nuclei
259  double xsec = (Z*fXSecAlgNCQE->Integral(inp) + N*fXSecAlgNCQE->Integral(inn))/A;
260 
261  // Add A dependence which is not known from theory
262  double fFracADep = 1.;
263  if(A>=12) fFracADep = TMath::Power((A/12.),fMECAPower-1.);
264 
265  // Use tunable fraction
266  // FFracNCQE is fraction of QE going to MEC
267  // fFracNCQE_cluster is fraction of MEC going to each NN pair
268  double fFracNCQE_cluster=0.;
269  if(nucleon_cluster_pdg==2000000200) fFracNCQE_cluster= 0.5*(1-fFracPN_NC); //n+n
270  if(nucleon_cluster_pdg==2000000201) fFracNCQE_cluster= fFracPN_NC; //n+p
271  if(nucleon_cluster_pdg==2000000202) fFracNCQE_cluster= 0.5*(1-fFracPN_NC); //p+p
272  xsec *= fFracNCQE*fFracNCQE_cluster*fFracADep;
273  delete inn;
274  delete inp;
275  return xsec;
276  }
277 
278  else if(isem) {
279  int nucpdg = kPdgProton;
280  // Create a tmp QE process
281  Interaction * inp = Interaction::QELEM(tgtpdg,nucpdg,nupdg,E);
282  nucpdg = kPdgNeutron;
283  // Create a tmp QE process
284  Interaction * inn = Interaction::QELEM(tgtpdg,nucpdg,nupdg,E);
285 
286  // Calculate cross section for the QE process - avg of p and n - best for isoscalar nuclei
287  double xsec = (Z*fXSecAlgEMQE->Integral(inp) + N*fXSecAlgEMQE->Integral(inn))/A;
288 
289  // Add A dependence which is not known from theory, data wants high A suppression
290  double fFracADep = 1.;
291  if(A>=12) fFracADep = TMath::Power((A/12.),fMECAPower-1.);
292 
293  // Use tunable fraction
294  // FFracEMQE is fraction of QE going to MEC
295  // fFracEMQE_cluster is fraction of MEC going to each NN pair
296  double fFracEMQE_cluster=0.;
297  if(nucleon_cluster_pdg==2000000200) fFracEMQE_cluster= 0.5*(1-fFracPN_EM); //n+n
298  if(nucleon_cluster_pdg==2000000201) fFracEMQE_cluster= fFracPN_EM; //n+p
299  if(nucleon_cluster_pdg==2000000202) fFracEMQE_cluster= 0.5*(1-fFracPN_EM); //p+p
300  xsec *= fFracEMQE*fFracEMQE_cluster*fFracADep;
301  delete inn;
302  delete inp;
303  return xsec;
304  }
305 
306  return 0;
307 }
308 //____________________________________________________________________________
309 bool EmpiricalMECPXSec2015::ValidProcess(const Interaction * interaction) const
310 {
311  if(interaction->TestBit(kISkipProcessChk)) return true;
312 
313  const ProcessInfo & proc_info = interaction->ProcInfo();
314  if(!proc_info.IsMEC()) return false;
315 
316  return true;
317 }
318 //____________________________________________________________________________
320 {
321  Algorithm::Configure(config);
322  this->LoadConfig();
323 }
324 //____________________________________________________________________________
326 {
327  Algorithm::Configure(config);
328 
329  Registry* algos = AlgConfigPool::Instance() -> GlobalParameterList() ;
330  string global_key_head = "XSecModel@genie::EventGenerator/" ;
331  string local_key_head = "XSecModel-" ;
332 
333  Registry r( "EmpiricalMECPXSec2015_specific", false ) ;
334  r.Set( local_key_head + "QEL-NC", algos -> GetAlg( global_key_head + "QEL-NC") ) ;
335  r.Set( local_key_head + "QEL-CC", algos -> GetAlg( global_key_head + "QEL-CC") ) ;
336  r.Set( local_key_head + "QEL-EM", algos -> GetAlg( global_key_head + "QEL-EM") ) ;
337 
339 
340  this->LoadConfig();
341 }
342 //____________________________________________________________________________
344 {
345  fXSecAlgCCQE = 0;
346  fXSecAlgNCQE = 0;
347  fXSecAlgEMQE = 0;
348 
349  GetParam( "EmpiricalMEC-Mq2d", fMq2d ) ;
350  GetParam( "EmpiricalMEC-Mass", fMass ) ;
351  GetParam( "EmpiricalMEC-Width", fWidth ) ;
352  GetParam( "EmpiricalMEC-APower", fMECAPower ) ;
353 
354  GetParam( "EmpiricalMEC-FracPN_NC", fFracPN_NC ) ;
355  GetParam( "EmpiricalMEC-FracPN_CC", fFracPN_CC ) ;
356  GetParam( "EmpiricalMEC-FracPN_EM", fFracPN_EM ) ;
357 
358  GetParam( "EmpiricalMEC-FracCCQE", fFracCCQE ) ;
359  GetParam( "EmpiricalMEC-FracNCQE", fFracNCQE ) ;
360  GetParam( "EmpiricalMEC-FracEMQE", fFracEMQE ) ;
361 
362  string key_head = "XSecModel-" ;
363 
364  fXSecAlgNCQE =
365  dynamic_cast<const XSecAlgorithmI *> ( this -> SubAlg( key_head + "QEL-NC" ) ) ;
367 
368  fXSecAlgCCQE =
369  dynamic_cast<const XSecAlgorithmI *> ( this -> SubAlg( key_head + "QEL-CC" ) ) ;
371 
372  fXSecAlgEMQE =
373  dynamic_cast<const XSecAlgorithmI *> ( this -> SubAlg( key_head + "QEL-EM" ) ) ;
375 }
376 //____________________________________________________________________________
Cross Section Calculation Interface.
double W(bool selected=false) const
Definition: Kinematics.cxx:167
Basic constants.
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Range1D_t InelWLim(double El, double ml, double M)
Definition: KineUtils.cxx:477
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
int HitNucPdg(void) const
Definition: Target.cxx:321
void Configure(const Registry &config)
Range1D_t InelWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:288
int A(void) const
Definition: Target.h:71
A simple [min,max] interval for doubles.
Definition: Range1.h:43
double fFracPN_NC
toy model param: fraction of nucleon pairs that are pn, not nn or pp
int Pdg(void) const
Definition: Target.h:72
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double Mass(Resonance_t res)
resonance mass (GeV)
Definition: config.py:1
const XSecAlgorithmI * fXSecAlgCCQE
cross section algorithm for CCQE
const XSecAlgorithmI * fXSecAlgNCQE
cross section algorithm for NCQE
enum genie::EKinePhaseSpace KinePhaseSpace_t
static Interaction * QELCC(int tgt, int nuc, int probe, double E=0)
double fMass
toy model param: peak of W distribution (GeV)
static const double kMinQ2Limit
Definition: Controls.h:41
Range1D_t InelQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:309
Float_t Z
Definition: plot.C:38
const double C
Summary information for an interaction.
Definition: Interaction.h:56
int iscc
bool IsWeakNC(void) const
double fFracPN_CC
toy model param: fraction of nucleon pairs that are pn, not nn or pp
#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
static const double kAem2
Definition: Constants.h:58
static string AsString(KinePhaseSpace_t kps)
static Interaction * QELNC(int tgt, int nuc, int probe, double E=0)
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1046
virtual double Integral(const Interaction *i) const =0
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
int ProbePdg(void) const
Definition: InitialState.h:65
double fWidth
toy model param: width of W distribution (GeV)
int Z(void) const
Definition: Target.h:69
double Gaus(TH1D *h, double &err, bool isTruth)
Definition: absCal.cxx:489
Misc GENIE control constants.
Double_t xsec[nknots]
Definition: testXsec.C:47
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsMEC(void) const
bool IsEM(void) const
double fMq2d
toy model param: `mass&#39; in dipole (Q2 - dependence) form factor (GeV)
double fFracPN_EM
toy model param: fraction of nucleon pairs that are pn, not nn or pp
static const double A
Definition: Units.h:82
double max
Definition: Range1.h:54
ifstream in
Definition: comparison.C:7
Range1D_t InelQ2Lim_W(double El, double ml, double M, double W)
Definition: KineUtils.cxx:499
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
const XSecAlgorithmI * fXSecAlgEMQE
cross section algorithm for EMQE
double Integral(const Interaction *i) const
static Interaction * QELEM(int tgt, int nuc, int probe, double E=0)
double fFracEMQE
empirical model param: MEC cross section is taken to be this fraction of Rosenbluth xs ...
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
double fFracNCQE
empirical model param: MEC cross section is taken to be this fraction of NCQE cross section ...
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double fFracCCQE
empirical model param: MEC cross section is taken to be this fraction of CCQE cross section ...
const int kPdgProton
Definition: PDGCodes.h:65
double fMECAPower
power of A relative to carbon
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
void kinematics()
Definition: kinematics.C:10
void Set(RgIMapPair entry)
Definition: Registry.cxx:282
#define W(x)
double ProbeE(RefFrame_t rf) const
const int kPdgNeutron
Definition: PDGCodes.h:67
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
static AlgConfigPool * Instance()
#define pDEBUG
Definition: Messenger.h:64
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353