IMDAnnihilationPXSec.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: Rosen Matev (r.matev@gmail.com)
8 
9  For the class documentation see the corresponding header file.
10 
11  Important revisions after version 2.8.0 :
12 
13 */
14 //____________________________________________________________________________
15 
26 
27 using namespace genie;
28 using namespace genie::constants;
29 using namespace genie::controls;
30 
31 //____________________________________________________________________________
33 XSecAlgorithmI("genie::IMDAnnihilationPXSec")
34 {
35 
36 }
37 //____________________________________________________________________________
39 XSecAlgorithmI("genie::IMDAnnihilationPXSec", config)
40 {
41 
42 }
43 //____________________________________________________________________________
45 {
46 
47 }
48 //____________________________________________________________________________
50  const Interaction * interaction, KinePhaseSpace_t kps) const
51 {
52  if(! this -> ValidProcess (interaction) ) return 0.;
53  if(! this -> ValidKinematics (interaction) ) return 0.;
54 
55  //----- get initial state & kinematics
56  const InitialState & init_state = interaction -> InitState();
57  const Kinematics & kinematics = interaction -> Kine();
58 
59  double Ev = init_state.ProbeE(kRfLab);
60  double twoMeEv = 2*kElectronMass*Ev;
61  double ymax = 1 - kMuonMass2/(twoMeEv + kElectronMass2);
62  double y = kinematics.y();
63  double A = kGF2/kPi;
64 
65  LOG("IMDAnnihilation", pDEBUG) << "Ev = " << Ev << ", y = " << y << ", ymax = " << ymax;
66  //Note: y = (Ev-El)/Ev but in Marciano's paper y=(El-(m_l^2+m_e^2)/2m_e)/Ev.
67  y = 1 - y - (kMuonMass2 + kElectronMass2)/twoMeEv;
68 
69  if(y > ymax) return 0;
70  if(y < 0) return 0;
71 
72  double xsec = A*(twoMeEv*TMath::Power(1-y,2) - (kMuonMass2 - kElectronMass2)*(1-y)); // <-- dxsec/dy
73 
74 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
75  LOG("IMDAnnihilation", pDEBUG)
76  << "*** dxsec(ve-)/dy [free e-](Ev="<< Ev << ", y= "<< y<< ") = "<< xsec;
77 #endif
78 
79  //----- The algorithm computes dxsec/dy
80  // Check whether variable tranformation is needed
81  if(kps!=kPSyfE) {
82  double J = utils::kinematics::Jacobian(interaction,kPSyfE,kps);
83  xsec *= J;
84  }
85 
86  //----- If requested return the free electron xsec even for nuclear target
87  if( interaction->TestBit(kIAssumeFreeElectron) ) return xsec;
88 
89  //----- Scale for the number of scattering centers at the target
90  int Ne = init_state.Tgt().Z(); // num of scattering centers
91  xsec *= Ne;
92 
93  return xsec;
94 }
95 //____________________________________________________________________________
96 double IMDAnnihilationPXSec::Integral(const Interaction * interaction) const
97 {
98  double xsec = fXSecIntegrator->Integrate(this,interaction);
99  return xsec;
100 }
101 //____________________________________________________________________________
102 bool IMDAnnihilationPXSec::ValidProcess(const Interaction * interaction) const
103 {
104  if(interaction->TestBit(kISkipProcessChk)) return true;
105  return true;
106 }
107 //____________________________________________________________________________
109 {
110  if(interaction->TestBit(kISkipKinematicChk)) return true;
111  return true;
112 }
113 //____________________________________________________________________________
115 {
116  Algorithm::Configure(config);
117  this->LoadConfig();
118 }
119 //____________________________________________________________________________
121 {
122  Algorithm::Configure(config);
123  this->LoadConfig();
124 }
125 //____________________________________________________________________________
127 {
128  // load XSec Integrator
130  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
132 }
133 //____________________________________________________________________________
134 
Cross Section Calculation Interface.
const double kPi
Basic constants.
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
Definition: config.py:1
void Configure(const Registry &config)
enum genie::EKinePhaseSpace KinePhaseSpace_t
static const double kElectronMass
Definition: Constants.h:71
double y(bool selected=false) const
Definition: Kinematics.cxx:122
Double_t ymax
Definition: plot.C:25
bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
static const double kElectronMass2
Definition: Constants.h:84
static const double kMuonMass2
Definition: Constants.h:85
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
const XSecIntegratorI * fXSecIntegrator
int Z(void) const
Definition: Target.h:69
Misc GENIE control constants.
Double_t xsec[nknots]
Definition: testXsec.C:47
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
double Integral(const Interaction *i) const
static const double A
Definition: Units.h:82
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
assert(nhit_max >=nhit_nbins)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
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
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