GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DMElectronPXSec.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 
7  Costas Andreopoulos <c.andreopoulos \at cern.ch>
8  University of Liverpool
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 
20 #include "Framework/Conventions/GBuild.h"
31 
32 using namespace genie;
33 using namespace genie::constants;
34 using namespace genie::controls;
35 
36 //____________________________________________________________________________
38 XSecAlgorithmI("genie::DMElectronPXSec")
39 {
40 
41 }
42 //____________________________________________________________________________
44 XSecAlgorithmI("genie::DMElectronPXSec", config)
45 {
46 
47 }
48 //____________________________________________________________________________
50 {
51 
52 }
53 //____________________________________________________________________________
55  const Interaction * interaction, KinePhaseSpace_t kps) const
56 {
57  if(! this -> ValidProcess (interaction) ) return 0.;
58  if(! this -> ValidKinematics (interaction) ) return 0.;
59 
60  //----- get initial state & kinematics
61  const InitialState & init_state = interaction -> InitState();
62  const Kinematics & kinematics = interaction -> Kine();
63  const ProcessInfo & proc_info = interaction -> ProcInfo();
64 
65  double Ev = init_state.ProbeE(kRfLab);
66  double Ev2 = TMath::Power(Ev,2);
67  double me = kElectronMass;
68  double me2 = TMath::Power(me,2);
69  double ml = interaction->FSPrimLepton()->Mass();
70  double ml2 = TMath::Power(ml,2);
71  double y = kinematics.y();
72  double MZ2 = TMath::Power(fMedMass,2);
73  double A = fgZp4 * TMath::Power(Ev,3) * me / (4.0 * kPi * (Ev2 - ml2) * TMath::Power(MZ2 + 2.0*Ev*me*(1.0 - y),2));
74 
75  // y = 1 - me/Ev - y; // FSPL = electron. XSec below are expressed in Marciano's y!
76  // if(y > 1/(1+0.5*me/Ev)) return 0;
77  // if(y < 0) return 0;
78 
79  double QeV2 = TMath::Power(0.5*(fQeL+fQeR),2);
80  double QeA2 = TMath::Power(0.5*(-fQeL+fQeR),2);
81  double QeVA = 0.25*(TMath::Power(fQeR,2) - TMath::Power(fQeL,2));
82 
83  double xsec = 0; // <-- dxsec/dy
84 
85  int inu = init_state.ProbePdg();
86 
87  if( proc_info.IsDarkMatter() && fVelMode == 0 )
88  {
89  double QdmV2 = TMath::Power(0.5*(fQdmL+fQdmR),2);
90  double QdmA2 = TMath::Power(0.5*(-fQdmL+fQdmR),2);
91  double QdmVA = 0.25*(TMath::Power(fQdmR,2) - TMath::Power(fQdmL,2));
92  double T1 = 1 + TMath::Power(y,2);
93  double T2 = (1.0 - y)* me / Ev;
94  double T3 = (1.0 - y)*ml2 / Ev / me;
95  double T4 = 2.0 * ml2 / Ev2;
96  double TL = 2.0*TMath::Power(2.0*(1.0 - y) * me * ml / MZ2 + ml/Ev,2);
97  xsec = (T1 - T2 - T3)*QdmV2*QeV2;
98  xsec += (T1 - T2 + T3 - T4)*QdmA2*QeV2;
99  xsec += (T1 + T2 - T3 - T4)*QdmV2*QeA2;
100  xsec += (T1 + T2 + T3 + T4 + TL)*QdmA2*QeA2;
101  if ( pdg::IsDarkMatter(inu) ) {
102  xsec += 4.0 * (1.0 - TMath::Power(y,2))*QdmVA*QeVA;
103  }
104  else {
105  xsec -= 4.0 * (1.0 - TMath::Power(y,2))*QdmVA*QeVA;
106  }
107  xsec *= A;
108  }
109 
110  if( proc_info.IsDarkMatter() && fVelMode == 2 )
111  {
112  double QdmS2 = TMath::Power(fQdmS,2);
113  double T1 = 2.0 * y;
114  double T2 = (1.0 - y)*me2 / Ev / me;
115  double T3 = (1.0 - y)*ml2 / Ev / me;
116  double T4 = 2.0 * ml2 / Ev2;
117  xsec = (T1 - T3)*QeV2*QdmS2;
118  xsec += (T1 - T2 - T3 - T4)*QeA2*QdmS2;
119  xsec *= A;
120  }
121 
122  #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
123  LOG("Elastic", pDEBUG)
124  << "*** dxsec(ve-)/dy [free e-](Ev="<< Ev << ", y= "<< y<< ") = "<< xsec;
125  #endif
126 
127  //----- The algorithm computes dxsec/dy
128  // Check whether variable tranformation is needed
129  if(kps!=kPSyfE) {
130  double J = utils::kinematics::Jacobian(interaction,kPSyfE,kps);
131  LOG("Elastic", pDEBUG) << "Multiplying by jacobian " << J;
132  xsec *= J;
133  }
134 
135  //----- If requested return the free electron xsec even for nuclear target
136  if( interaction->TestBit(kIAssumeFreeElectron) ) return xsec;
137 
138  //----- Scale for the number of scattering centers at the target
139  int Ne = init_state.Tgt().Z(); // num of scattering centers
140  xsec *= Ne;
141  LOG("Elastic", pDEBUG) << "Multiplying by Ne " << Ne;
142 
143  return xsec;
144 }
145 //____________________________________________________________________________
146 double DMElectronPXSec::Integral(const Interaction * interaction) const
147 {
148  double xsec = fXSecIntegrator->Integrate(this,interaction);
149  return xsec;
150 }
151 //____________________________________________________________________________
152 bool DMElectronPXSec::ValidProcess(const Interaction * interaction) const
153 {
154  if(interaction->TestBit(kISkipProcessChk)) return true;
155  return true;
156 }
157 //____________________________________________________________________________
158 bool DMElectronPXSec::ValidKinematics(const Interaction* interaction) const
159 {
160  if(interaction->TestBit(kISkipKinematicChk)) return true;
161  return true;
162 }
163 //____________________________________________________________________________
165 {
166  Algorithm::Configure(config);
167  this->LoadConfig();
168 }
169 //____________________________________________________________________________
170 void DMElectronPXSec::Configure(string config)
171 {
172  Algorithm::Configure(config);
173  this->LoadConfig();
174 }
175 //____________________________________________________________________________
177 {
178  // Dark matter couplings
179  double gZp;
180  GetParam("ZpCoupling", gZp);
181  GetParam("DarkLeftCharge", fQdmL);
182  GetParam("DarkRightCharge", fQdmR);
183  GetParam("DarkScalarCharge", fQdmS);
184  GetParam("ElectronLeftCharge", fQeL);
185  GetParam("ElectronRightCharge", fQeR);
186 
187  fgZp4 = TMath::Power(gZp, 4);
188 
189  // Mediator mass
191 
192  // velocity dependence of interaction
193  GetParamDef("velocity-mode", fVelMode, 0 );
194 
195  // load XSec Integrator
197  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
198  assert(fXSecIntegrator);
199 }
200 //____________________________________________________________________________
201 
Cross Section Calculation Interface.
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
Cross Section Integrator Interface.
const int kPdgMediator
Definition: PDGCodes.h:220
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:127
enum genie::EKinePhaseSpace KinePhaseSpace_t
static const double kElectronMass
double y(bool selected=false) const
Definition: Kinematics.cxx:112
const XSecIntegratorI * fXSecIntegrator
Summary information for an interaction.
Definition: Interaction.h:56
void Configure(const Registry &config)
#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 A
Definition: Units.h:74
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int ProbePdg(void) const
Definition: InitialState.h:64
int Z(void) const
Definition: Target.h:68
double Integral(const Interaction *i) const
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
bool IsDarkMatter(void) const
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
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 UInt_t kIAssumeFreeElectron
Definition: Interaction.h:50
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345