AhrensNCELPXSec.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  @ Sep 19, 2009 - CA
14  NuNucElasticPXSec -> AhrensNCELPXSec
15 
16 */
17 //____________________________________________________________________________
18 
19 #include <TMath.h>
20 
32 
33 using namespace genie;
34 using namespace genie::utils;
35 using namespace genie::constants;
36 
37 //____________________________________________________________________________
39 XSecAlgorithmI("genie::AhrensNCELPXSec")
40 {
41 
42 }
43 //____________________________________________________________________________
45 XSecAlgorithmI("genie::AhrensNCELPXSec", config)
46 {
47 
48 }
49 //____________________________________________________________________________
51 {
52 
53 }
54 //____________________________________________________________________________
56  const Interaction * interaction, KinePhaseSpace_t kps) const
57 {
58  if(! this -> ValidProcess (interaction) ) return 0.;
59  if(! this -> ValidKinematics (interaction) ) return 0.;
60 
61  const InitialState & init_state = interaction -> InitState();
62  const Kinematics & kinematics = interaction -> Kine();
63  const Target & target = init_state.Tgt();
64 
65  double E = init_state.ProbeE(kRfHitNucRest);
66  double Q2 = kinematics.Q2();
67  double M = target.HitNucMass();
68  double M2 = TMath::Power(M, 2.);
69  double E2 = TMath::Power(E, 2.);
70  double qmv2 = TMath::Power(1 + Q2/fMv2, 2);
71  double qma2 = TMath::Power(1 + Q2/fMa2, 2);
72 
73  //-- handle terms changing sign for antineutrinos and isospin rotations
74  int nusign = 1;
75  int nucsign = 1;
76  int nupdgc = init_state.ProbePdg();
77  int nucpdgc = target.HitNucPdg();
78  if( pdg::IsAntiNeutrino(nupdgc) ) nusign = -1;
79  if( pdg::IsNeutron(nucpdgc) ) nucsign = -1;
80 
81  //-- compute isoscalar form factor terms
82  double Ge0 = 1.5 * fkGamma / qmv2;
83  double Gm0 = 1.5 * fkGamma * (fMuP+fMuN) / qmv2;
84 
85  //-- compute isovector form factor terms
86  double Ge1 = 0.5 * fkAlpha / qmv2;
87  double Gm1 = 0.5 * fkAlpha * (fMuP-fMuN) / qmv2;
88  double Ga1 = -0.5 * fFa0 * (1 + (nucsign) * fEta) / qma2;
89 
90  //-- compute form factors
91  double Ge = Ge0 + (nucsign) * Ge1;
92  double Gm = Gm0 + (nucsign) * Gm1;
93  double Ga = (nucsign) * Ga1;
94  double Ge2 = TMath::Power(Ge,2);
95  double Gm2 = TMath::Power(Gm,2);
96  double Ga2 = TMath::Power(Ga,2);
97 
98  //-- compute the free nucleon cross section
99  double tau = 0.25 * Q2/M2;
100  double fa = 1-M*tau/E;
101  double fa2 = TMath::Power(fa,2);
102  double fb = tau*(tau+1)*M2/E2;
103  double A = (Ge2/(1+tau)) * (fa2-fb);
104  double B = (Ga2 + tau*Gm2/(1+tau)) * (fa2+fb);
105  double C = 4*tau*(M/E)*Gm*Ga * fa;
106  double xsec0 = 0.5*kGF2/kPi;
107  double xsec = xsec0 * (A + B + (nusign)*C);
108 
109  LOG("AhrensNCEL", pDEBUG)
110  << "dXSec[vN,El]/dQ2 [FreeN](Ev = "<< E<< ", Q2 = "<< Q2 << ") = "<< xsec;
111 
112  //-- The algorithm computes dxsec/dQ2
113  // Check whether variable tranformation is needed
114  if(kps!=kPSQ2fE) {
115  double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
116  xsec *= J;
117  }
118 
119  //-- if requested return the free nucleon xsec even for input nuclear tgt
120  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
121 
122  //-- compute nuclear suppression factor
123  // (R(Q2) is adapted from NeuGEN - see comments therein)
124  double R = nuclear::NuclQELXSecSuppression("Default", 0.5, interaction);
125 
126  //-- number of scattering centers in the target
127  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
128 
129  LOG("AhrensNCEL", pDEBUG)
130  << "Nuclear suppression factor R(Q2) = " << R << ", NNucl = " << NNucl;
131 
132  //-- compute nuclear cross section
133  xsec *= (R*NNucl);
134 
135  return xsec;
136 }
137 //____________________________________________________________________________
138 double AhrensNCELPXSec::Integral(const Interaction * interaction) const
139 {
140  double xsec = fXSecIntegrator->Integrate(this,interaction);
141  return xsec;
142 }
143 //____________________________________________________________________________
144 bool AhrensNCELPXSec::ValidProcess(const Interaction * interaction) const
145 {
146  if(interaction->TestBit(kISkipProcessChk)) return true;
147  return true;
148 }
149 //____________________________________________________________________________
151 {
152  Algorithm::Configure(config);
153  this->LoadConfig();
154 }
155 //____________________________________________________________________________
157 {
158  Algorithm::Configure(config);
159  this->LoadConfig();
160 }
161 //____________________________________________________________________________
163 {
164  // alpha and gamma
165  double thw ;
166  GetParam( "WeinbergAngle", thw ) ;
167  double sin2thw = TMath::Power(TMath::Sin(thw), 2);
168  fkAlpha = 1.-2.*sin2thw;
169  fkGamma = -0.66666667*sin2thw;
170 
171  // eta and Fa(q2=0)
172  GetParam( "EL-Axial-Eta", fEta ) ;
173  GetParam( "QEL-FA0", fFa0 ) ;
174 
175  // axial and vector masses
176  double ma, mv ;
177  GetParam( "QEL-Ma", ma ) ;
178  GetParam( "QEL-Mv", mv ) ;
179  fMa2 = TMath::Power(ma,2);
180  fMv2 = TMath::Power(mv,2);
181 
182  // anomalous magnetic moments
183  GetParam( "AnomMagnMoment-P", fMuP ) ;
184  GetParam( "AnomMagnMoment-N", fMuN ) ;
185 
186  // load XSec Integrator
188  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
190 }
191 //____________________________________________________________________________
Cross Section Calculation Interface.
const double kPi
Basic constants.
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?
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
double HitNucMass(void) const
Definition: Target.cxx:250
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
Definition: config.py:1
enum genie::EKinePhaseSpace KinePhaseSpace_t
const XSecIntegratorI * fXSecIntegrator
const double C
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
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
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
int ProbePdg(void) const
Definition: InitialState.h:65
#define R(x)
int Z(void) const
Definition: Target.h:69
Double_t xsec[nknots]
Definition: testXsec.C:47
static const double A
Definition: Units.h:82
int N(void) const
Definition: Target.h:70
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
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 Integral(const Interaction *i) const
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
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
void Configure(const Registry &config)
Root of GENIE utility namespaces.
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:49
#define pDEBUG
Definition: Messenger.h:64
TFile fa("Li7.root")
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353