NuElectronPXSec.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  @ Feb 12, 2013 - CA (code from Rosen Matev)
14  Tweak kinematics. The final state primary lepton is always the electron.
15  The kinematical variable y has the definition used in Marciano's paper.
16 */
17 //____________________________________________________________________________
18 
29 
30 using namespace genie;
31 using namespace genie::constants;
32 using namespace genie::controls;
33 
34 //____________________________________________________________________________
36 XSecAlgorithmI("genie::NuElectronPXSec")
37 {
38 
39 }
40 //____________________________________________________________________________
42 XSecAlgorithmI("genie::NuElectronPXSec", config)
43 {
44 
45 }
46 //____________________________________________________________________________
48 {
49 
50 }
51 //____________________________________________________________________________
53  const Interaction * interaction, KinePhaseSpace_t kps) const
54 {
55  if(! this -> ValidProcess (interaction) ) return 0.;
56  if(! this -> ValidKinematics (interaction) ) return 0.;
57 
58  //----- get initial state & kinematics
59  const InitialState & init_state = interaction -> InitState();
60  const Kinematics & kinematics = interaction -> Kine();
61  const ProcessInfo & proc_info = interaction -> ProcInfo();
62 
63  double Ev = init_state.ProbeE(kRfLab);
64  double me = kElectronMass;
65  double y = kinematics.y();
66  double A = kGF2*2*me*Ev/kPi;
67 
68  y = 1 - me/Ev - y; // FSPL = electron. XSec below are expressed in Marciano's y!
69  if(y > 1/(1+0.5*me/Ev)) return 0;
70  if(y < 0) return 0;
71 
72  double xsec = 0; // <-- dxsec/dy
73 
74  int inu = init_state.ProbePdg();
75 
76  // nue + e- -> nue + e- [CC + NC + interference]
77  if(pdg::IsNuE(inu))
78  {
79  double em = -0.5 - fSin28w;
80  double ep = -fSin28w;
81  xsec = TMath::Power(em,2) + TMath::Power(ep*(1-y),2) - ep*em*me*y/Ev;
82  xsec *= A;
83  }
84 
85  // nuebar + e- -> nue + e- [CC + NC + interference]
86  if(pdg::IsAntiNuE(inu))
87  {
88  double em = -0.5 - fSin28w;
89  double ep = -fSin28w;
90  xsec = TMath::Power(ep,2) + TMath::Power(em*(1-y),2) - ep*em*me*y/Ev;
91  xsec *= A;
92  }
93 
94  // numu/nutau + e- -> numu/nutau + e- [NC]
95  if( (pdg::IsNuMu(inu)||pdg::IsNuTau(inu)) && proc_info.IsWeakNC() )
96  {
97  double em = 0.5 - fSin28w;
98  double ep = -fSin28w;
99  xsec = TMath::Power(em,2) + TMath::Power(ep*(1-y),2) - ep*em*me*y/Ev;
100  xsec *= A;
101  }
102 
103  // numubar/nutaubar + e- -> numubar/nutaubar + e- [NC]
104  if( (pdg::IsAntiNuMu(inu)||pdg::IsAntiNuTau(inu)) && proc_info.IsWeakNC() )
105  {
106  double em = 0.5 - fSin28w;
107  double ep = -fSin28w;
108  xsec = TMath::Power(ep,2) + TMath::Power(em*(1-y),2) - ep*em*me*y/Ev;
109  xsec *= A;
110  }
111 
112  // numu/nutau + e- -> l- + nu_e [CC}
113  if( (pdg::IsNuMu(inu)||pdg::IsNuTau(inu)) && proc_info.IsWeakCC() ) xsec=0;
114 /*
115  double ml = (pdg::IsNuMu(inu)) ? kMuonMass : kTauMass;
116  double ml2 = TMath::Power(ml,2);
117  xsec = (kGF2*s/kPi)*(1-ml2/s);
118  xsec = TMath::Max(0.,xsec); // if s<ml2 => xsec<0 : force to xsec=0
119 */
120 
121 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
122  LOG("Elastic", pDEBUG)
123  << "*** dxsec(ve-)/dy [free e-](Ev="<< Ev << ", y= "<< y<< ") = "<< xsec;
124 #endif
125 
126  //----- The algorithm computes dxsec/dy
127  // Check whether variable tranformation is needed
128  if(kps!=kPSyfE) {
129  double J = utils::kinematics::Jacobian(interaction,kPSyfE,kps);
130  xsec *= J;
131  }
132 
133  //----- If requested return the free electron xsec even for nuclear target
134  if( interaction->TestBit(kIAssumeFreeElectron) ) return xsec;
135 
136  //----- Scale for the number of scattering centers at the target
137  int Ne = init_state.Tgt().Z(); // num of scattering centers
138  xsec *= Ne;
139 
140  return xsec;
141 }
142 //____________________________________________________________________________
143 double NuElectronPXSec::Integral(const Interaction * interaction) const
144 {
145  double xsec = fXSecIntegrator->Integrate(this,interaction);
146  return xsec;
147 }
148 //____________________________________________________________________________
149 bool NuElectronPXSec::ValidProcess(const Interaction * interaction) const
150 {
151  if(interaction->TestBit(kISkipProcessChk)) return true;
152  return true;
153 }
154 //____________________________________________________________________________
155 bool NuElectronPXSec::ValidKinematics(const Interaction* interaction) const
156 {
157  if(interaction->TestBit(kISkipKinematicChk)) return true;
158  return true;
159 }
160 //____________________________________________________________________________
162 {
163  Algorithm::Configure(config);
164  this->LoadConfig();
165 }
166 //____________________________________________________________________________
168 {
169  Algorithm::Configure(config);
170  this->LoadConfig();
171 }
172 //____________________________________________________________________________
174 {
175  // weinberg angle
176  double thw ;
177  GetParam( "WeinbergAngle", thw ) ;
178  fSin28w = TMath::Power(TMath::Sin(thw), 2);
179  fSin48w = TMath::Power(TMath::Sin(thw), 4);
180 
181  // load XSec Integrator
183  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
185 }
186 //____________________________________________________________________________
187 
Cross Section Calculation Interface.
const double kPi
bool IsNuTau(int pdgc)
Definition: PDGUtils.cxx:160
Basic constants.
bool IsWeakCC(void) const
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.
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
bool IsAntiNuTau(int pdgc)
Definition: PDGUtils.cxx:175
bool IsNuE(int pdgc)
Definition: PDGUtils.cxx:150
const XSecIntegratorI * fXSecIntegrator
Definition: config.py:1
enum genie::EKinePhaseSpace KinePhaseSpace_t
static const double kElectronMass
Definition: Constants.h:71
double Integral(const Interaction *i) const
double y(bool selected=false) const
Definition: Kinematics.cxx:122
bool IsNuMu(int pdgc)
Definition: PDGUtils.cxx:155
Summary information for an interaction.
Definition: Interaction.h:56
bool IsWeakNC(void) const
Definition: NueSkimmer.h:24
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
int ProbePdg(void) const
Definition: InitialState.h:65
int Z(void) const
Definition: Target.h:69
Misc GENIE control constants.
bool IsAntiNuMu(int pdgc)
Definition: PDGUtils.cxx:170
Double_t xsec[nknots]
Definition: testXsec.C:47
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
static const double A
Definition: Units.h:82
bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
void Configure(const Registry &config)
assert(nhit_max >=nhit_nbins)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
static const double kGF2
Definition: Constants.h:60
void kinematics()
Definition: kinematics.C:10
double ProbeE(RefFrame_t rf) const
bool IsAntiNuE(int pdgc)
Definition: PDGUtils.cxx:165
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 UInt_t kIAssumeFreeElectron
Definition: Interaction.h:50
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353