AhrensDMELPXSec.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: Joshua Berger <jberger \at physics.wisc.edu>
8  University of Wisconsin-Madison
9 
10  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
11  University of Liverpool & STFC Rutherford Appleton Lab
12 */
13 //____________________________________________________________________________
14 
15 #include <TMath.h>
16 
17 //#include "Framework/Algorithm/AlgConfigPool.h"
30 #include "Framework/Utils/RunOpt.h"
31 
32 using namespace genie;
33 using namespace genie::utils;
34 using namespace genie::constants;
35 
36 //____________________________________________________________________________
38 XSecAlgorithmI("genie::AhrensDMELPXSec")
39 {
40 
41 }
42 //____________________________________________________________________________
44 XSecAlgorithmI("genie::AhrensDMELPXSec", 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  const InitialState & init_state = interaction -> InitState();
61  const Kinematics & kinematics = interaction -> Kine();
62  const Target & target = init_state.Tgt();
63 
64  LOG("AhrensDMEL", pDEBUG) << "Using v^" << fVelMode << " dependence";
65 
66  double E = init_state.ProbeE(kRfHitNucRest);
67  double ml = init_state.GetProbeP4(kRfHitNucRest)->M();
68  double Q2 = kinematics.Q2();
69  double M = target.HitNucMass();
70  double M2 = TMath::Power(M, 2.);
71  double E2 = TMath::Power(E, 2.);
72  double ml2 = TMath::Power(ml,2.);
73  double qma2 = TMath::Power(1 + Q2/fMa2, 2);
74 
75  //-- handle terms changing sign for antineutrinos and isospin rotations
76  // int nusign = 1; // comment-out unused variable to eliminate warnings
77  int nucsign = 1;
78  // int nupdgc = init_state.ProbePdg(); // // comment-out unused variable to eliminate warnings
79  int nucpdgc = target.HitNucPdg();
80  if( pdg::IsNeutron(nucpdgc) ) nucsign = -1;
81 
82  //-- compute axial form factor terms
83  // For dark matter scattering, the factor of 1/2 should be removed
84  // This factor should go to ~ 1 in the limit q^2 -> 0
85  // The 1/2 is an isospin piece that is irrelevant for DM
86  double Ga1 = -fFa0 * (1 + (nucsign) * fEta) / qma2;
87 
88  //-- compute form factors
89  double Ga = (nucsign) * Ga1;
90  double Ga2 = TMath::Power(Ga,2);
91 
92  //-- compute the free nucleon cross section
93  double xsec = 0.;
94  double tau = 0.25 * Q2/M2;
95  double fa = 1-M*tau/E;
96  double fa2 = TMath::Power(fa,2);
97  double fb = tau*(tau+1)*M2/E2;
98  double fc = (tau+2)*ml2/E2;
99  LOG("AhrensDMEL", pDEBUG)
100  << "Using a mediator mass of " << fMedMass;
101  double MZ2 = TMath::Power(fMedMass,2);
102  double fd = 8*ml2*M2*tau*(2*M2*tau+MZ2) / (MZ2*MZ2*E2);
103  // double gZp = RunOpt::Instance()->ZpCoupling();
104  double gZp = fgZp;
105  double gZp4 = TMath::Power(gZp,4);
106  double prop = 1. / (Q2 + MZ2);
107  double prop2 = TMath::Power(prop,2);
108  double B = 0.;
109  double xsec0 = 0.;
110  switch (fVelMode) {
111  // v^0 fermion cross-section
112  case 0:
113  B = (Ga2) * (fa2+fb+fc+fd);
114  xsec0 = 0.25*prop2*gZp4*E2/kPi/(E2-ml2);
115  xsec = xsec0 * (B);
116  break;
117  // v^2 scalar cross-section
118  case 2:
119  fc = fc - ml2/E2;
120  B = (Ga2) * (fa2-fb-fc);
121  xsec0 = 0.25*prop2*gZp4*E2/kPi/(E2-ml2);
122  xsec = xsec0 * (B);
123  break;
124  }
125 
126  LOG("AhrensDMEL", pDEBUG)
127  << "dXSec[vN,El]/dQ2 [FreeN](Ev = "<< E<< ", Q2 = "<< Q2 << ") = "<< xsec;
128 
129  //-- The algorithm computes dxsec/dQ2
130  // Check whether variable tranformation is needed
131  if(kps!=kPSQ2fE) {
132  double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
133  xsec *= J;
134  }
135 
136  //-- if requested return the free nucleon xsec even for input nuclear tgt
137  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
138 
139  //-- compute nuclear suppression factor
140  // (R(Q2) is adapted from NeuGEN - see comments therein)
141  double R = nuclear::NuclQELXSecSuppression("Default", 0.5, interaction);
142 
143  //-- number of scattering centers in the target
144  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
145 
146  LOG("AhrensDMEL", pDEBUG)
147  << "Nuclear suppression factor R(Q2) = " << R << ", NNucl = " << NNucl;
148 
149  //-- compute nuclear cross section
150  xsec *= (R*NNucl);
151 
152  return xsec;
153 }
154 //____________________________________________________________________________
155 double AhrensDMELPXSec::Integral(const Interaction * interaction) const
156 {
157  double xsec = fXSecIntegrator->Integrate(this,interaction);
158  return xsec;
159 }
160 //____________________________________________________________________________
161 bool AhrensDMELPXSec::ValidProcess(const Interaction * interaction) const
162 {
163  if(interaction->TestBit(kISkipProcessChk)) return true;
164  return true;
165 }
166 //____________________________________________________________________________
168 {
169  Algorithm::Configure(config);
170  this->LoadConfig();
171 }
172 //____________________________________________________________________________
174 {
175  Algorithm::Configure(config);
176  this->LoadConfig();
177 }
178 //____________________________________________________________________________
180 {
181  // alpha and gamma
182  double thw ;
183  this->GetParam( "WeinbergAngle", thw ) ;
184  double sin2thw = TMath::Power(TMath::Sin(thw), 2);
185  fkAlpha = 1.-2.*sin2thw;
186  fkGamma = -0.66666667*sin2thw;
187 
188  // eta and Fa(q2=0)
189  this->GetParam( "EL-Axial-Eta", fEta ) ;
190  this->GetParam( "QEL-FA0", fFa0 ) ;
191 
192  // axial and vector masses
193  double ma, mv ;
194  this->GetParam( "QEL-Ma", ma ) ;
195  this->GetParam( "QEL-Mv", mv ) ;
196  fMa2 = TMath::Power(ma,2);
197  fMv2 = TMath::Power(mv,2);
198 
199  // anomalous magnetic moments
200  this->GetParam( "AnomMagnMoment-P", fMuP ) ;
201  this->GetParam( "AnomMagnMoment-N", fMuN ) ;
202 
203  // velocity dependence of interaction
204  this->GetParamDef("velocity-mode", fVelMode, 0 );
205 
206  // mediator coupling
207  this->GetParam("ZpCoupling", fgZp ) ;
208 
209  // mediator mass
211 
212  // load XSec Integrator
214  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
216 }
217 //____________________________________________________________________________
Cross Section Calculation Interface.
const double kPi
Basic constants.
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
int HitNucPdg(void) const
Definition: Target.cxx:321
const int kPdgMediator
Definition: PDGCodes.h:196
double HitNucMass(void) const
Definition: Target.cxx:250
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
Definition: config.py:1
enum genie::EKinePhaseSpace KinePhaseSpace_t
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
Summary information for an interaction.
Definition: Interaction.h:56
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
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
double Integral(const Interaction *i) const
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
#define R(x)
int Z(void) const
Definition: Target.h:69
Double_t xsec[nknots]
Definition: testXsec.C:47
const XSecIntegratorI * fXSecIntegrator
int N(void) const
Definition: Target.h:70
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
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 Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
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:67
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
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...
Root of GENIE utility namespaces.
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
void Configure(const Registry &config)
double NuclQELXSecSuppression(string kftable, double pmax, const Interaction *in)
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
TFile fa("Li7.root")
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353