GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AlamSimoAtharVacasSKXSec.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  Chris Marshall and Martti Nirkko
7 */
8 //____________________________________________________________________________
9 
10 #include <TMath.h>
11 #include <Math/IFunction.h>
12 #include <Math/IntegratorMultiDim.h>
13 #include "Math/AdaptiveIntegratorMultiDim.h"
14 
15 #include "Framework/Conventions/GBuild.h"
28 #include "Framework/Utils/Range1.h"
33 
34 using namespace genie;
35 using namespace genie::constants;
36 using namespace genie::controls;
37 using namespace genie::utils;
38 
39 //____________________________________________________________________________
41 XSecIntegratorI("genie::AlamSimoAtharVacasSKXSec")
42 {
43 
44 }
45 //____________________________________________________________________________
47 XSecIntegratorI("genie::AlamSimoAtharVacasSKXSec", config)
48 {
49 
50 }
51 //____________________________________________________________________________
53 {
54 
55 }
56 //____________________________________________________________________________
58  const XSecAlgorithmI * model, const Interaction * in) const
59 {
60  LOG("SKXSec", pDEBUG) << "Integrating the Alam Simo Athar Vacas model";
61 
62  const InitialState & init_state = in -> InitState();
63 
64  if(! model->ValidProcess(in) ) return 0.;
65 
66  const KPhaseSpace & kps = in->PhaseSpace(); // only OK phase space for this
67  if(!kps.IsAboveThreshold()) {
68  LOG("SKXSec", pDEBUG) << "*** Below energy threshold";
69  return 0;
70  }
71 
72  // If the input interaction is off a nuclear target, then chek whether
73  // the corresponding free nucleon cross section already exists at the
74  // cross section spline list.
75  // Cross section for PP scales with number of protons, NP and NN scale
76  // with number of neutrons
77  int nucpdgc = init_state.Tgt().HitNucPdg();
78  int NNucl = (pdg::IsProton(nucpdgc)) ?
79  init_state.Tgt().Z() : init_state.Tgt().N();
80  double Ev = init_state.ProbeE(kRfHitNucRest);
81 
83  if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() ) {
84  Interaction * interaction = new Interaction(*in);
85  Target * target = interaction->InitStatePtr()->TgtPtr();
86  if(pdg::IsProton(nucpdgc)) { target->SetId(kPdgTgtFreeP); }
87  else { target->SetId(kPdgTgtFreeN); }
88  if(xsl->SplineExists(model,interaction)) {
89  const Spline * spl = xsl->GetSpline(model, interaction);
90  double xsec = spl->Evaluate(Ev);
91  LOG("SKXSec", pINFO)
92  << "From XSecSplineList: XSec[SK,free nucleon] (E = " << Ev << " GeV) = " << xsec;
93  if(! interaction->TestBit(kIAssumeFreeNucleon) ) {
94  xsec *= NNucl;
95  LOG("SKXSec", pINFO) << "XSec[SK] (E = " << Ev << " GeV) = " << xsec;
96  }
97  delete interaction;
98  return xsec;
99  }
100  delete interaction;
101  }
102 
103  // no free nucelon spline exists -- do the integration
104 
105  // Check this
106  double Enu = init_state.ProbeE(kRfLab);
107  int kpdg = in->ExclTag().StrangeHadronPdg();
108  double mk = PDGLibrary::Instance()->Find(kpdg)->Mass();
109  double ml = PDGLibrary::Instance()->Find(in->FSPrimLeptonPdg())->Mass();
110 
111  // integration bounds for T (kinetic energy)
112  double zero = 0.0;
113  double tmax = Enu - mk - ml;
114 
115  LOG("SKXSec", pDEBUG)
116  << "Lepton/Kaon KE integration range = [" << 0.0 << ", " << tmax << "]";
117 
118  Interaction * interaction = new Interaction(*in);
119  interaction->SetBit(kISkipProcessChk);
120  interaction->SetBit(kISkipKinematicChk);
121 
122  double xsec = 0;
123 
124  // do the integration over log(1-costheta) so it's not so sharply peaked
125 
126  ROOT::Math::IBaseFunctionMultiDim * func =
127  new utils::gsl::d3Xsec_dTldTkdCosThetal(model, interaction);
128  double kine_min[3] = { zero, zero, -20 }; // Tlep, Tkaon, cosine theta lep
129  double kine_max[3] = { tmax, tmax, 0.69314718056 }; // Tlep, Tkaon, cosine theta lep
130 
131  ROOT::Math::IntegrationMultiDim::Type ig_type =
133 
134  double abstol = 1; //We mostly care about relative tolerance.
135  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
136 
137  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
138  delete func;
139 
140  delete interaction;
141 
142  return xsec;
143 }
144 //____________________________________________________________________________
146 {
147  Algorithm::Configure(config);
148  this->LoadConfig();
149 }
150 //____________________________________________________________________________
152 {
153  Algorithm::Configure(config);
154  this->LoadConfig();
155 }
156 //____________________________________________________________________________
158 {
159  // Get GSL integration type & relative tolerance
160  this->GetParamDef("gsl-integration-type" , fGSLIntgType, string("vegas") );
161  this->GetParamDef("gsl-max-evals", fGSLMaxEval, 20000);
162  this->GetParamDef("gsl-relative-tolerance", fGSLRelTol, 0.01);
163  this->GetParamDef("split-integral", fSplitIntegral, true);
164 }
165 //____________________________________________________________________________
166 
167 //_____________________________________________________________________________
168 // GSL wrappers
169 //____________________________________________________________________________
170 
172  const XSecAlgorithmI * m, const Interaction * i) :
173 ROOT::Math::IBaseFunctionMultiDim(),
174 fModel(m),
175 fInteraction(i)
176 {
177 
178 }
179 //____________________________________________________________________________
181 {
182 
183 }
184 //____________________________________________________________________________
186 {
187  // phi_kq is important for kinematics generation
188  // But dependence is weak so we will not use it in the integration
189  return 3;
190 }
191 //____________________________________________________________________________
193 {
194 // inputs:
195 // Tl [GeV]
196 // Tk [GeV]
197 // cosine theta l
198 // * calculate phi_kq based on neutrino energy -- this is for the integral only
199 // outputs:
200 // differential cross section [10^-38 cm^2]
201 //
202 
203  double Enu = fInteraction->InitState().ProbeE(kRfLab);
204 
205  double phikq = 0.5*constants::kPi;
206  if( Enu > 3.0 ) phikq = 0.55*constants::kPi;
207  else if( Enu > 1.0 ) phikq = constants::kPi*(0.5 + 0.025*(Enu-1.0));
208 
209  Kinematics * kinematics = fInteraction->KinePtr();
210 
211  double T_l = xin[0];
212  double T_k = xin[1];
213  //double cos_theta_l = xin[2];
214  double log_oneminuscostheta = xin[2];
215  double cos_theta_l = 1.0 - TMath::Exp(log_oneminuscostheta);
216  double J = 1.0 - cos_theta_l; // Jacobian for transformation
217 
218  kinematics->SetKV(kKVTl, T_l);
219  kinematics->SetKV(kKVTk, T_k);
220  kinematics->SetKV(kKVctl, cos_theta_l);
221  kinematics->SetKV(kKVphikq, phikq);
222 
223  double xsec = fModel->XSec(fInteraction);
224  LOG( "GXSecFunc", pDEBUG )
225  << "t_l = " << T_l << " t_k = " << T_k
226  << " costhetal = " << cos_theta_l << " phikq = " << phikq
227  << " enu = " << Enu << " Xsec = " << xsec;
228 
229  // integrate out phi_kq by multiplying by 2pi
230 
231  xsec *= 2.0 * genie::constants::kPi * J;
232 
233  return xsec/(1E-38 * units::cm2);
234 }
235 //____________________________________________________________________________
236 ROOT::Math::IBaseFunctionMultiDim *
238 {
239  return
240  new genie::utils::gsl::d3Xsec_dTldTkdCosThetal(fModel,fInteraction);
241 }
242 //____________________________________________________________________________
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
string fGSLIntgType
name of GSL numerical integrator
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
Cross Section Integrator Interface.
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:59
int HitNucPdg(void) const
Definition: Target.cxx:304
bool SplineExists(const XSecAlgorithmI *alg, const Interaction *i) const
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:58
bool IsNucleus(void) const
Definition: Target.cxx:272
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double Mass(Resonance_t res)
resonance mass (GeV)
static XSecSplineList * Instance()
double Evaluate(double x) const
Definition: Spline.cxx:363
void SetId(int pdgc)
Definition: Target.cxx:149
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
Summary information for an interaction.
Definition: Interaction.h:56
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
int StrangeHadronPdg(void) const
Definition: XclsTag.h:55
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool IsEmpty(void) const
static constexpr double cm2
Definition: Units.h:69
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
Kinematical phase space.
Definition: KPhaseSpace.h:33
const int kPdgTgtFreeN
Definition: PDGCodes.h:200
const int kPdgTgtFreeP
Definition: PDGCodes.h:199
d3Xsec_dTldTkdCosThetal(const XSecAlgorithmI *m, const Interaction *i)
double func(double x, double y)
int Z(void) const
Definition: Target.h:68
#define pINFO
Definition: Messenger.h:62
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
void Configure(const Registry &config)
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual bool ValidProcess(const Interaction *i) const =0
Can this cross section algorithm handle the input process?
int fGSLMaxEval
GSL max evaluations.
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
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
Target * TgtPtr(void) const
Definition: InitialState.h:67
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
Definition: InitialState.h:66
const Spline * GetSpline(const XSecAlgorithmI *alg, const Interaction *i) const
List of cross section vs energy splines.
double ProbeE(RefFrame_t rf) const
static constexpr double m
Definition: Units.h:71
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
double fGSLRelTol
required relative tolerance (error)