SmithMonizQELCCXSec.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2017, GENIE Neutrino MC Generator Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Igor Kakorin <kakorin@jinr.ru>, Joint Institute for Nuclear Research
8  adapted from fortran code provided by
9  Konstantin Kuzmin <kkuzmin@theor.jinr.ru>,
10  Joint Institute for Nuclear Research, Institute for Theoretical and Experimental Physics
11  Vladimir Lyubushkin,
12  Joint Institute for Nuclear Research
13  Vadim Naumov <vnaumov@theor.jinr.ru>,
14  Joint Institute for Nuclear Research
15  based on code of Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
16  University of Liverpool & STFC Rutherford Appleton Lab
17 
18  For the class documentation see the corresponding header file.
19 
20 
21 */
22 //____________________________________________________________________________
23 
24 #include <sstream>
25 
26 #include <TMath.h>
27 #include <Math/IFunction.h>
28 #include <Math/Integrator.h>
29 
38 #include "Framework/Utils/Range1.h"
41 
42 using namespace genie;
43 using std::ostringstream;
44 
45 //____________________________________________________________________________
47 XSecIntegratorI("genie::SmithMonizQELCCXSec")
48 {
49 
50 }
51 //____________________________________________________________________________
53 XSecIntegratorI("genie::SmithMonizQELCCXSec", config)
54 {
55 
56 }
57 //____________________________________________________________________________
59 {
60 
61 }
62 //____________________________________________________________________________
64  const XSecAlgorithmI * model, const Interaction * in) const
65 {
66  LOG("SMQELXSec",pDEBUG) << "Beginning integrate";
67  if(! model->ValidProcess(in)) return 0.;
68 
69  const InitialState & init_state = in -> InitState();
70  const Target & target = init_state.Tgt();
71  if (target.A()<3)
72  {
73  const KPhaseSpace & kps = in->PhaseSpace();
74  if(!kps.IsAboveThreshold()) {
75  LOG("SMQELXSec", pDEBUG) << "*** Below energy threshold";
76  return 0;
77  }
78  Range1D_t rQ2 = kps.Limits(kKVQ2);
79  if(rQ2.min<0 || rQ2.max<0) return 0;
80  Interaction * interaction = new Interaction(*in);
81  interaction->SetBit(kISkipProcessChk);
82  interaction->SetBit(kISkipKinematicChk);
83  ROOT::Math::IBaseFunctionOneDim * func = new utils::gsl::dXSec_dQ2_E(model, interaction);
85  double abstol = 0; //We mostly care about relative tolerance
86  ROOT::Math::Integrator ig(*func,ig_type,abstol,fGSLRelTol,fGSLMaxSizeOfSubintervals, fGSLRule);
87  double xsec = ig.Integral(rQ2.min, rQ2.max) * (1E-38 * units::cm2);
88  delete func;
89  delete interaction;
90  return xsec;
91  }
92  else
93  {
94  Interaction * interaction = new Interaction(*in);
96  if (interaction->InitState().ProbeE(kRfLab)<sm_utils->E_nu_thr_SM()) return 0;
97  interaction->SetBit(kISkipProcessChk);
98  interaction->SetBit(kISkipKinematicChk);
99  double xsec = 0;
100 
101 
102  ROOT::Math::IBaseFunctionMultiDim * func = new utils::gsl::d2Xsec_dQ2dv(model, interaction);
103  double kine_min[2] = { 0, 0};
104  double kine_max[2] = { 1, 1};
105 
108 
109  double abstol = 0; //We mostly care about relative tolerance.
110  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol2D, fGSLMaxEval);
111 
112  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
113  delete func;
114  delete interaction;
115 
116  return xsec;
117 
118 
119  }
120 
121 }
122 //____________________________________________________________________________
124 {
125  Algorithm::Configure(config);
126  this->LoadConfig();
127 }
128 //____________________________________________________________________________
130 {
131  Algorithm::Configure(config);
132 
133  Registry r("SmithMonizQELCCXSec_specific", false ) ;
134  r.Set("sm_utils_algo", RgAlg("genie::SmithMonizUtils","Default") ) ;
135 
137 
138  this->LoadConfig();
139 }
140 //____________________________________________________________________________
142 {
143 
144  // Get GSL integration type & relative tolerance
145  GetParamDef( "gsl-integration-type", fGSLIntgType, string("gauss") ) ;
146  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1e-3 ) ;
147  int max_size_of_subintervals;
148  GetParamDef( "gsl-max-size-of-subintervals", max_size_of_subintervals, 40000);
149  fGSLMaxSizeOfSubintervals = (unsigned int) max_size_of_subintervals;
150  int rule;
151  GetParamDef( "gsl-rule", rule, 3);
152  fGSLRule = (unsigned int) rule;
153  if (fGSLRule>6) fGSLRule=3;
154  GetParamDef( "gsl-integration-type-2D", fGSLIntgType2D, string("adaptive") ) ;
155  GetParamDef( "gsl-relative-tolerance-2D", fGSLRelTol2D, 1e-7 ) ;
156  GetParamDef( "gsl-max-eval", fGSLMaxEval, 1000000000 ) ;
157 
158  sm_utils = const_cast<genie::SmithMonizUtils *>(
159  dynamic_cast<const genie::SmithMonizUtils *>(
160  this -> SubAlg("sm_utils_algo") ) );
161 }
162 
163 
164 //_____________________________________________________________________________
165 // GSL wrappers
166 //____________________________________________________________________________
168 ROOT::Math::IBaseFunctionMultiDim(),
169 fModel(m),
170 fInteraction(interaction)
171 {
172  AlgFactory * algf = AlgFactory::Instance();
173  sm_utils = const_cast<genie::SmithMonizUtils *>(dynamic_cast<const genie::SmithMonizUtils *>(algf->GetAlgorithm("genie::SmithMonizUtils","Default")));
174  sm_utils->SetInteraction(interaction);
175 }
176 //____________________________________________________________________________
178 {
179 
180 }
181 //____________________________________________________________________________
183 {
184  return 2;
185 }
186 //____________________________________________________________________________
187 double genie::utils::gsl::d2Xsec_dQ2dv::DoEval(const double * xin) const
188 {
189 // inputs:
190 // normalized Q2 from 0 to 1
191 // normalized v from 0 to 1
192 // outputs:
193 // differential cross section [10^-38 cm^2]
194 //
195 
197  double Q2 = (rQ2.max-rQ2.min)*xin[0]+rQ2.min;
198  Range1D_t rv = sm_utils->vQES_SM_lim(Q2);
199  double v = (rv.max-rv.min)*xin[1]+rv.min;
200  double J = (rQ2.max-rQ2.min)*(rv.max-rv.min); // Jacobian for transformation
201 
203  kinematics->SetKV(kKVQ2, Q2);
204  kinematics->SetKV(kKVv, v);
205 
207 
208  xsec *= J;
209 
210  return xsec/(1E-38 * units::cm2);
211 }
212 //____________________________________________________________________________
213 ROOT::Math::IBaseFunctionMultiDim *
215 {
216  return
218 }
219 
220 
void SetInteraction(const Interaction *i)
Cross Section Calculation Interface.
Range1D_t Q2QES_SM_lim(void) const
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition: GSLUtils.cxx:31
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.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
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.
double E_nu_thr_SM(void) const
int A(void) const
Definition: Target.h:71
A simple [min,max] interval for doubles.
Definition: Range1.h:43
double fGSLRelTol2D
required relative tolerance (error) for 2D integrator
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
Definition: config.py:1
unsigned int fGSLRule
GSL Gauss-Kronrod integration rule (only for GSL 1D adaptive type)
Range1D_t vQES_SM_lim(double Q2) const
static const double cm2
Definition: Units.h:77
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
string fGSLIntgType2D
name of GSL 2D numerical integrator
void Configure(const Registry &config)
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
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 Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:86
double func(double x, double y)
Double_t xsec[nknots]
Definition: testXsec.C:47
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:345
double max
Definition: Range1.h:54
int fGSLMaxEval
GSL max evaluations.
ifstream in
Definition: comparison.C:7
Contains auxiliary functions for Smith-Moniz model. Is a concrete implementation of the XSecAlgorit...
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
unsigned int fGSLMaxSizeOfSubintervals
GSL maximum number of sub-intervals for 1D integrator.
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.
double DoEval(const double *xin) const
virtual bool ValidProcess(const Interaction *i) const =0
Can this cross section algorithm handle the input process?
TRandom3 r(0)
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:53
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
Definition: InitialState.h:67
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
Float_t e
Definition: plot.C:35
void kinematics()
Definition: kinematics.C:10
const XML_Char XML_Content * model
Definition: expat.h:151
void Set(RgIMapPair entry)
Definition: Registry.cxx:282
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:49
#define pDEBUG
Definition: Messenger.h:64
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
double fGSLRelTol
required relative tolerance (error)
d2Xsec_dQ2dv(const XSecAlgorithmI *m, const Interaction *i)
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353