GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AhrensNCELPXSec.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 
24 
25 using namespace genie;
26 using namespace genie::utils;
27 using namespace genie::constants;
28 
29 //____________________________________________________________________________
31 XSecAlgorithmI("genie::AhrensNCELPXSec")
32 {
33 
34 }
35 //____________________________________________________________________________
37 XSecAlgorithmI("genie::AhrensNCELPXSec", config)
38 {
39 
40 }
41 //____________________________________________________________________________
43 {
44 
45 }
46 //____________________________________________________________________________
48  const Interaction * interaction, KinePhaseSpace_t kps) const
49 {
50  if(! this -> ValidProcess (interaction) ) return 0.;
51  if(! this -> ValidKinematics (interaction) ) return 0.;
52 
53  const InitialState & init_state = interaction -> InitState();
54  const Kinematics & kinematics = interaction -> Kine();
55  const Target & target = init_state.Tgt();
56 
57  double E = init_state.ProbeE(kRfHitNucRest);
58  double Q2 = kinematics.Q2();
59  double M = target.HitNucMass();
60  double M2 = TMath::Power(M, 2.);
61  double E2 = TMath::Power(E, 2.);
62  double qmv2 = TMath::Power(1 + Q2/fMv2, 2);
63  double qma2 = TMath::Power(1 + Q2/fMa2, 2);
64 
65  //-- handle terms changing sign for antineutrinos and isospin rotations
66  int nusign = 1;
67  int nucsign = 1;
68  int nupdgc = init_state.ProbePdg();
69  int nucpdgc = target.HitNucPdg();
70  if( pdg::IsAntiNeutrino(nupdgc) ) nusign = -1;
71  if( pdg::IsNeutron(nucpdgc) ) nucsign = -1;
72 
73  //-- compute isoscalar form factor terms
74  double Ge0 = 1.5 * fkGamma / qmv2;
75  double Gm0 = 1.5 * fkGamma * (fMuP+fMuN) / qmv2;
76 
77  //-- compute isovector form factor terms
78  double Ge1 = 0.5 * fkAlpha / qmv2;
79  double Gm1 = 0.5 * fkAlpha * (fMuP-fMuN) / qmv2;
80  double Ga1 = -0.5 * fFa0 * (1 + (nucsign) * fEta) / qma2;
81 
82  //-- compute form factors
83  double Ge = Ge0 + (nucsign) * Ge1;
84  double Gm = Gm0 + (nucsign) * Gm1;
85  double Ga = (nucsign) * Ga1;
86  double Ge2 = TMath::Power(Ge,2);
87  double Gm2 = TMath::Power(Gm,2);
88  double Ga2 = TMath::Power(Ga,2);
89 
90  //-- compute the free nucleon cross section
91  double tau = 0.25 * Q2/M2;
92  double fa = 1-M*tau/E;
93  double fa2 = TMath::Power(fa,2);
94  double fb = tau*(tau+1)*M2/E2;
95  double A = (Ge2/(1+tau)) * (fa2-fb);
96  double B = (Ga2 + tau*Gm2/(1+tau)) * (fa2+fb);
97  double C = 4*tau*(M/E)*Gm*Ga * fa;
98  double xsec0 = 0.5*kGF2/kPi;
99  double xsec = xsec0 * (A + B + (nusign)*C);
100 
101  LOG("AhrensNCEL", pDEBUG)
102  << "dXSec[vN,El]/dQ2 [FreeN](Ev = "<< E<< ", Q2 = "<< Q2 << ") = "<< xsec;
103 
104  //-- The algorithm computes dxsec/dQ2
105  // Check whether variable tranformation is needed
106  if(kps!=kPSQ2fE) {
107  double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
108  xsec *= J;
109  }
110 
111  //-- if requested return the free nucleon xsec even for input nuclear tgt
112  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
113 
114  //-- compute nuclear suppression factor
115  // (R(Q2) is adapted from NeuGEN - see comments therein)
116  double R = nuclear::NuclQELXSecSuppression("Default", 0.5, interaction);
117 
118  //-- number of scattering centers in the target
119  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
120 
121  LOG("AhrensNCEL", pDEBUG)
122  << "Nuclear suppression factor R(Q2) = " << R << ", NNucl = " << NNucl;
123 
124  //-- compute nuclear cross section
125  xsec *= (R*NNucl);
126 
127  return xsec;
128 }
129 //____________________________________________________________________________
130 double AhrensNCELPXSec::Integral(const Interaction * interaction) const
131 {
132  double xsec = fXSecIntegrator->Integrate(this,interaction);
133  return xsec;
134 }
135 //____________________________________________________________________________
136 bool AhrensNCELPXSec::ValidProcess(const Interaction * interaction) const
137 {
138  if(interaction->TestBit(kISkipProcessChk)) return true;
139  return true;
140 }
141 //____________________________________________________________________________
143 {
144  Algorithm::Configure(config);
145  this->LoadConfig();
146 }
147 //____________________________________________________________________________
148 void AhrensNCELPXSec::Configure(string config)
149 {
150  Algorithm::Configure(config);
151  this->LoadConfig();
152 }
153 //____________________________________________________________________________
155 {
156  // alpha and gamma
157  double thw ;
158  GetParam( "WeinbergAngle", thw ) ;
159  double sin2thw = TMath::Power(TMath::Sin(thw), 2);
160  fkAlpha = 1.-2.*sin2thw;
161  fkGamma = -0.66666667*sin2thw;
162 
163  // eta and Fa(q2=0)
164  GetParam( "EL-Axial-Eta", fEta ) ;
165  GetParam( "QEL-FA0", fFa0 ) ;
166 
167  // axial and vector masses
168  double ma, mv ;
169  GetParam( "QEL-Ma", ma ) ;
170  GetParam( "QEL-Mv", mv ) ;
171  fMa2 = TMath::Power(ma,2);
172  fMv2 = TMath::Power(mv,2);
173 
174  // anomalous magnetic moments
175  GetParam( "AnomMagnMoment-P", fMuP ) ;
176  GetParam( "AnomMagnMoment-N", fMuN ) ;
177 
178  // load XSec Integrator
180  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
181  assert(fXSecIntegrator);
182 }
183 //____________________________________________________________________________
Cross Section Calculation Interface.
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
int HitNucPdg(void) const
Definition: Target.cxx:304
double HitNucMass(void) const
Definition: Target.cxx:233
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
enum genie::EKinePhaseSpace KinePhaseSpace_t
const XSecIntegratorI * fXSecIntegrator
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
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:336
#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
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
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
int N(void) const
Definition: Target.h:69
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
double Integral(const Interaction *i) const
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double ProbeE(RefFrame_t rf) const
void Configure(const Registry &config)
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
double NuclQELXSecSuppression(string kftable, double pmax, const Interaction *in)
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345