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