MECXSec.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  For the class documentation see the corresponding header file.
8 
9 */
10 //____________________________________________________________________________
11 
12 #include <TMath.h>
13 #include <Math/IFunction.h>
14 #include <Math/IntegratorMultiDim.h>
15 #include "Math/AdaptiveIntegratorMultiDim.h"
16 
33 #include "Framework/Utils/Range1.h"
36 
37 using namespace genie;
38 using namespace genie::constants;
39 using namespace genie::controls;
40 using namespace genie::utils;
41 
42 //____________________________________________________________________________
44 XSecIntegratorI("genie::MECXSec")
45 {
46 
47 }
48 //____________________________________________________________________________
50 XSecIntegratorI("genie::MECXSec", config)
51 {
52 
53 }
54 //____________________________________________________________________________
56 {
57 
58 }
59 //____________________________________________________________________________
61  const XSecAlgorithmI * model, const Interaction * in) const
62 {
63  if(! model->ValidProcess(in) ) return 0.;
64 
65  const KPhaseSpace & kps = in->PhaseSpace(); // only OK phase space for this
66  if(!kps.IsAboveThreshold()) {
67  LOG("MECXSec", pDEBUG) << "*** Below energy threshold";
68  return 0;
69  }
70 
71  Interaction * interaction = new Interaction(*in);
72  interaction->SetBit(kISkipProcessChk);
73  interaction->SetBit(kISkipKinematicChk);
74 
75  // T, costh limits
76  double Enu = in->InitState().ProbeE(kRfLab);
77  double LepMass = in->FSPrimLepton()->Mass();
78  double TMax = Enu - LepMass;
79  double TMin = 0.0;
80  double CosthMax = 1.0;
81  double CosthMin = -1.0;
82  if (Enu < fQ3Max) {
83  TMin = 0 ;
84  CosthMin = -1 ;
85  } else {
86  TMin = TMath::Sqrt(TMath::Power(LepMass, 2) + TMath::Power((Enu - fQ3Max), 2)) - LepMass;
87  CosthMin = TMath::Sqrt(1 - TMath::Power((fQ3Max / Enu ), 2));
88  }
89 
90  double kine_min[2] = { TMin, CosthMin };
91  double kine_max[2] = { TMax, CosthMax };
92 
93  double xsec = 0;
94 
95  double abstol = 1; //We mostly care about relative tolerance.
96  ROOT::Math::IBaseFunctionMultiDim * func =
97  new utils::gsl::d2Xsec_dTCosth(model, interaction);
100  ROOT::Math::IntegratorMultiDim ig(
101  *func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
102 
103  xsec = ig.Integral(kine_min, kine_max);
104 
105  delete func;
106  delete interaction;
107 
108  return xsec;
109 }
110 //____________________________________________________________________________
112 {
113  Algorithm::Configure(config);
114  this->LoadConfig();
115 }
116 //____________________________________________________________________________
118 {
119  Algorithm::Configure(config);
120  this->LoadConfig();
121 }
122 //____________________________________________________________________________
124 {
125  GetParam( "NSV-Q3Max", fQ3Max ) ;
126 
127  // Get GSL integration type & relative tolerance
128  GetParamDef( "gsl-integration-type", fGSLIntgType, string("vegas") ) ;
129 
130  int max ;
131  GetParamDef( "gsl-max-eval", max, 20000 ) ;
132  fGSLMaxEval = (unsigned int) max ;
133 
134  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 0.01 ) ;
135  GetParamDef( "split-integral", fSplitIntegral, true ) ;
136 
137 }
138 //_____________________________________________________________________________
139 // GSL wrappers
140 //____________________________________________________________________________
142  const XSecAlgorithmI * m, const Interaction * i) :
143 ROOT::Math::IBaseFunctionMultiDim(),
144 fModel(m),
145 fInteraction(i)
146 {
147 
148 }
149 //____________________________________________________________________________
151 {
152 
153 }
154 //____________________________________________________________________________
156 {
157  return 2;
158 }
159 //____________________________________________________________________________
160 double genie::utils::gsl::d2Xsec_dTCosth::DoEval(const double * xin) const
161 {
162 // inputs:
163 // T [GeV]
164 // cos(theta)
165 // outputs:
166 // differential cross section (hbar=c=1 units)
167 //
168 
169  double T = xin[0];
170  double costh = xin[1];
171 
173  kinematics->SetKV(kKVTl, T);
174  kinematics->SetKV(kKVctl, costh);
175 
176  double xsec = fModel->XSec(fInteraction, kPSTlctl);
177  return xsec;
178 }
179 //____________________________________________________________________________
180 ROOT::Math::IBaseFunctionMultiDim *
182 {
183  return
185 }
186 //____________________________________________________________________________
187 
188 
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
virtual ~MECXSec()
Definition: MECXSec.cxx:55
string fGSLIntgType
name of GSL numerical integrator
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:67
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.
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
Definition: config.py:1
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
Definition: MECXSec.cxx:181
unsigned int NDim(void) const
Definition: MECXSec.cxx:155
bool fSplitIntegral
Definition: MECXSec.h:50
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
Definition: MECXSec.cxx:60
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double DoEval(const double *xin) const
Definition: MECXSec.cxx:160
static keras::KerasModel * fModel
Definition: NusVarsTemp.cxx:27
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
Kinematical phase space.
Definition: KPhaseSpace.h:34
void Configure(const Registry &config)
Definition: MECXSec.cxx:111
double fQ3Max
Definition: MECXSec.h:54
void LoadConfig(void)
Definition: MECXSec.cxx:123
double func(double x, double y)
d2Xsec_dTCosth(const XSecAlgorithmI *m, const Interaction *i)
Definition: MECXSec.cxx:141
Misc GENIE control constants.
Double_t xsec[nknots]
Definition: testXsec.C:47
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
const Interaction * fInteraction
Definition: MECXSec.h:77
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:345
int fGSLMaxEval
GSL max evaluations.
ifstream in
Definition: comparison.C:7
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
const XSecAlgorithmI * fModel
Definition: MECXSec.h:76
virtual bool ValidProcess(const Interaction *i) const =0
Can this cross section algorithm handle the input process?
const InitialState & InitState(void) const
Definition: Interaction.h:69
double T
Definition: Xdiff_gwt.C:5
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void kinematics()
Definition: kinematics.C:10
const XML_Char XML_Content * model
Definition: expat.h:151
double ProbeE(RefFrame_t rf) const
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
#define pDEBUG
Definition: Messenger.h:64
double fGSLRelTol
required relative tolerance (error)