GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PattonCEvNSPXSec.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 #include <Math/Integrator.h>
13 
28 
29 using namespace genie;
30 using namespace genie::utils;
31 using namespace genie::constants;
32 
33 //____________________________________________________________________________
35 XSecAlgorithmI("genie::PattonCEvNSPXSec")
36 {
37 
38 }
39 //____________________________________________________________________________
41 XSecAlgorithmI("genie::PattonCEvNSPXSec", config)
42 {
43 
44 }
45 //____________________________________________________________________________
47 {
48 
49 }
50 //____________________________________________________________________________
52  const Interaction * interaction, KinePhaseSpace_t kps) const
53 {
54  if(! this -> ValidProcess (interaction) ) return 0.;
55  if(! this -> ValidKinematics (interaction) ) return 0.;
56 
57  const InitialState & init_state = interaction -> InitState();
58  const Kinematics & kinematics = interaction -> Kine();
59  const Target & target = init_state.Tgt();
60 
61  // User inputs to the calculation
62  double E = init_state.ProbeE(kRfLab); // neutrino energy, units: GeV
63  double Q2 = kinematics.Q2(); // momentum transfer, units: GeV^2
64  int Z = target.Z(); // number of protons
65  int N = target.N(); // number of nucleons
66 
67  // Target atomic mass number and mass calculated from inputs
68  int A = Z+N;
69  int target_nucleus_pdgc = pdg::IonPdgCode(A,Z);
70  double M = PDGLibrary::Instance()->Find(target_nucleus_pdgc)->Mass(); // units: GeV
71  LOG("CEvNS", pDEBUG) << "M = " << M << " GeV";
72 
73  // Calculation of nuclear recoil kinetic energy computed from input Q2
74  double TA = Q2*E / (2*E*M+Q2); // nuclear recoil kinetic energy
75 
76  LOG("CEvNS", pDEBUG)
77  << "Q2 = " << Q2 << " GeV^2, E = " << E << " GeV "
78  << "--> TA = " << TA << " GeV";
79 
80  // auxiliary variables
81  double E2 = E*E;
82  double TA2 = TA*TA;
83  double Q4 = Q2*Q2;
84  double Q6 = Q2*Q4;
85 
86  // Calculation of weak charge
87  // double Qw = N - Z*(1-fSin2thw);
88  // Qw^2/4 x-section factor in arXiv:1207.0693v1 not needed here.
89  // 1/4 was absorbed in the constant front factor (below) and Qw^2 factor would
90  // have cancelled with ignored 1/Qw factor in the form factor F.
91 
92  // Calculation of nuclear density moments used for the evaluation
93  // of the neutron form factor
94  double avg_density = this->NuclearDensityMoment(A, 0); // units:: fm^-3
95  double Rn2 = this->NuclearDensityMoment(A, 2) / avg_density; // units: fm^2
96  double Rn4 = this->NuclearDensityMoment(A, 4) / avg_density; // units: fm^4
97  double Rn6 = this->NuclearDensityMoment(A, 6) / avg_density; // units: fm^6
98 
99  LOG("CEvNS", pDEBUG)
100  << "Nuclear density moments:"
101  << " <Rn^2> = " << Rn2 << " fm^2,"
102  << " <Rn^4> = " << Rn4 << " fm^4,"
103  << " <Rn^6> = " << Rn6 << " fm^6";
104 
105  Rn2 *= TMath::Power(units::fm, 2.); // units: GeV^-2
106  Rn4 *= TMath::Power(units::fm, 4.); // units: GeV^-4
107  Rn6 *= TMath::Power(units::fm, 6.); // units: GeV^-6
108 
109  // Calculation of proton form factor
110  // Form factor is neglected since it is multiplied with a small factor 1-4sin^2(\theta_{w})
111  double Fp = 0; // units: -
112  // Calculation of neutron form factor
113  // Using a Taylor expansion of sin(Qr) and keeping the first three terms (shown to be
114  // sufficient for approximating the full Fn calculation, even for heavy nuclei)
115  double Fn = N * (1 - Q2*Rn2/6. + Q4*Rn4/120. - Q6*Rn6/5040.); // units: -
116  // Overall form factor
117  double F = (Fn - (1-4*fSin2thw)*Fp); // units: -
118  F = TMath::Max(0.,F);
119  double F2 = F*F; // units: -
120 
121  LOG("CEvNS", pDEBUG)
122  << "Form factors: Fp = " << Fp << ", Fn = " << Fn << ", F = " << F;
123 
124  // dsig/dTA calculation
125  double const_factor = 0.125*kGF2/kPi; // units: GeV^-4
126  double kinematic_term = M * (2 - 2*TA/E + TA2/E2 - M*TA/E2); // units: GeV
127  kinematic_term = TMath::Max(0., kinematic_term);
128 
129  LOG("CEvNS", pDEBUG)
130  << "kinematic term: " << kinematic_term;
131 
132  double xsec = const_factor * kinematic_term * F2; // units: GeV^-3 (area/GeV)
133 
134  LOG("CEvNS", pINFO)
135  << "dsig[vA,CEvNS]/dTA (Ev = "
136  << E << " GeV, Q2 = "<< Q2 << " GeV^2; TA = " << TA << " GeV) = "
137  << xsec/(units::cm2) << " cm^2/GeV";
138 
139  // The algorithm computes dxsec/dTA
140  // Check whether variable tranformation is needed
141  if(kps!=kPSTAfE) {
142  // TODO: Move the calculation in utils::kinematics
143  // double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
144  double J = 0;
145  if(kps==kPSQ2fE) {
146  J = 2*E2*M / TMath::Power(2*E*M+Q2, 2.); // units: GeV^-1
147  }
148  xsec *= J; // units: GeV^-4 (area/GeV^2)
149  }
150 
151  return xsec;
152 }
153 //____________________________________________________________________________
155 {
156  // Calculate moments of the nuclear density
157  // Inputs:
158  // - atomic mass number, A
159  // - integer k specifying required nuclear density moment
160  // Output:
161  // - nuclear density moment in units of fm^k
162  //
163  // THINGS TO DO:
164  // 1) The calculation can be stored, as it is required only once per nucleus.
165  // The calculation is very fast so it doesn't matter.
166 
167  ROOT::Math::IBaseFunctionOneDim * integrand = new
169 
170  ROOT::Math::IntegrationOneDim::Type ig_type =
172 
173  double R0 = utils::nuclear::Radius(A); // units: fm
174  double rmin = 0; // units: fm
175  double rmax = fNuclDensMomentCalc_UpperIntegrationLimit * R0; // units: fm
176 
177  ROOT::Math::Integrator ig(
178  *integrand,ig_type,
182  double moment = 2 * constants::kPi * ig.Integral(rmin, rmax); // units: fm^k
183 
184  delete integrand;
185 
186  return moment;
187 }
188 //____________________________________________________________________________
189 double PattonCEvNSPXSec::Integral(const Interaction * interaction) const
190 {
191  double xsec = fXSecIntegrator->Integrate(this,interaction);
192  return xsec;
193 }
194 //____________________________________________________________________________
195 bool PattonCEvNSPXSec::ValidProcess(const Interaction * interaction) const
196 {
197  if(interaction->TestBit(kISkipProcessChk)) return true;
198 
199  const ProcessInfo & proc_info = interaction->ProcInfo();
200  if(!proc_info.IsCoherentElastic()) return false;
201 
202  const InitialState & init_state = interaction->InitState();
203  const Target & target = init_state.Tgt();
204  if(!target.IsNucleus()) return false;
205 
206  return true;
207 }
208 //____________________________________________________________________________
210 {
211  Algorithm::Configure(config);
212  this->LoadConfig();
213 }
214 //____________________________________________________________________________
215 void PattonCEvNSPXSec::Configure(string config)
216 {
217  Algorithm::Configure(config);
218  this->LoadConfig();
219 }
220 //____________________________________________________________________________
222 {
223  double thw = 0.;
224  this->GetParam("WeinbergAngle", thw);
225  fSin2thw = TMath::Power(TMath::Sin(thw), 2.);
226 
227  this->GetParamDef(
228  "nuclear-density-moment-gsl-upper-limit",
230  10.); // in nuclear radii
231  this->GetParamDef(
232  "nuclear-density-moment-gsl-rel-tol",
234  1E-3);
235  this->GetParamDef(
236  "nuclear-density-moment-gsl-abs-tol",
238  1.);
239  this->GetParamDef(
240  "nuclear-density-moment-gsl-max-eval",
242  10000);
243 
245  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
246  assert(fXSecIntegrator);
247 }
248 //____________________________________________________________________________
Cross Section Calculation Interface.
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition: GSLUtils.cxx:23
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double fNuclDensMomentCalc_AbsoluteTolerance
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
double fNuclDensMomentCalc_RelativeTolerance
bool IsNucleus(void) const
Definition: Target.cxx:272
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
const XSecIntegratorI * fXSecIntegrator
cross section integrator
enum genie::EKinePhaseSpace KinePhaseSpace_t
double Integral(const Interaction *i) const
double Radius(int A, double Ro=constants::kNucRo)
void Configure(const Registry &config)
double NuclearDensityMoment(int A, int k) const
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?
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool IsCoherentElastic(void) const
static constexpr double A
Definition: Units.h:74
static constexpr double cm2
Definition: Units.h:69
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
double fSin2thw
sin^2(weinberg angle)
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 Z(void) const
Definition: Target.h:68
#define pINFO
Definition: Messenger.h:62
double fNuclDensMomentCalc_UpperIntegrationLimit
int N(void) const
Definition: Target.h:69
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
Integrand for the calculation of the k^th nuclear density moment: {0}^{} {A}(r) r^k d^{3}r where {A}(...
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
static constexpr double fm
Definition: Units.h:75
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
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
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
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:48
#define pDEBUG
Definition: Messenger.h:63
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345