RESXSec.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: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Sep 07, 2009 - CA
14  Integrated with GNU Numerical Library (GSL) via ROOT's MathMore library.
15 
16 */
17 //____________________________________________________________________________
18 
19 #include <TMath.h>
20 #include <Math/IFunction.h>
21 #include <Math/IntegratorMultiDim.h>
22 
34 
35 using namespace genie;
36 using namespace genie::constants;
37 
38 //____________________________________________________________________________
40 XSecIntegratorI("genie::RESXSec")
41 {
42 
43 }
44 //____________________________________________________________________________
46 XSecIntegratorI("genie::RESXSec", config)
47 {
48 
49 }
50 //____________________________________________________________________________
52 {
53 
54 }
55 //____________________________________________________________________________
57  const XSecAlgorithmI * model, const Interaction * in) const
58 {
59  if(! model->ValidProcess(in) ) return 0.;
60 
61  const KPhaseSpace & kps = in->PhaseSpace();
62  if(!kps.IsAboveThreshold()) {
63  LOG("COHXSec", pDEBUG) << "*** Below energy threshold";
64  return 0;
65  }
66  Range1D_t Q2l = kps.Limits(kKVQ2);
67  Range1D_t Wl = kps.Limits(kKVW);
68 
69  Interaction * interaction = new Interaction(*in);
70  interaction->SetBit(kISkipProcessChk);
71  //interaction->SetBit(kISkipKinematicChk);
72 
73  ROOT::Math::IBaseFunctionMultiDim * func =
74  new utils::gsl::d2XSec_dWdQ2_E(model, interaction);
75 
78  double abstol = 1E-16; //We mostly care about relative tolerance.
79 
80  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
81 
82  double kine_min[2] = { Wl.min, Q2l.min };
83  double kine_max[2] = { Wl.max, Q2l.max };
84  double xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
85 
86  LOG("RESXSec", pERROR) << "Integrator opt / Integrator = " << ig.Options().Integrator();
87 
88  if(xsec < 0) {
89  LOG("RESXSec", pERROR) << "Algorithm " << *model << " returns a negative cross-section (xsec = " << xsec << " 1E-38 * cm2)";
90  LOG("RESXSec", pERROR) << "for process" << *interaction;
91  LOG("RESXSec", pERROR) << "Integrator status code = " << ig.Status();
92  LOG("RESXSec", pERROR) << "Integrator error code = " << ig.Error();
93  }
94 
95  //LOG("RESXSec", pINFO) << "XSec[RES] (Ev = " << Ev << " GeV) = " << xsec;
96 
97  delete interaction;
98  delete func;
99  return xsec;
100 }
101 //____________________________________________________________________________
103 {
104  Algorithm::Configure(config);
105  this->LoadConfig();
106 }
107 //____________________________________________________________________________
109 {
110  Algorithm::Configure(config);
111  this->LoadConfig();
112 }
113 //____________________________________________________________________________
115 {
116  // Get GSL integration type & relative tolerance
117  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
118  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-2 ) ;
119  int max, min ;
120  GetParamDef( "gsl-max-eval", max, 500000 ) ;
121  GetParamDef( "gsl-min-eval", min, 5000 ) ;
122  fGSLMaxEval = (unsigned int) max ;
123  fGSLMinEval = (unsigned int) min ;
124 }
125 //____________________________________________________________________________
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
string fGSLIntgType
name of GSL numerical integrator
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
Cross Section Integrator Interface.
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:67
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Definition: config.py:1
static const double cm2
Definition: Units.h:77
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
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
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
Kinematical phase space.
Definition: KPhaseSpace.h:34
double func(double x, double y)
Double_t xsec[nknots]
Definition: testXsec.C:47
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
void Configure(const Registry &config)
Definition: RESXSec.cxx:102
double max
Definition: Range1.h:54
int fGSLMaxEval
GSL max evaluations.
ifstream in
Definition: comparison.C:7
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
virtual ~RESXSec()
Definition: RESXSec.cxx:51
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?
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
Definition: RESXSec.cxx:56
double min
Definition: Range1.h:53
bool GetParamDef(const RgKey &name, T &p, const T &def) const
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void LoadConfig(void)
Definition: RESXSec.cxx:114
const XML_Char XML_Content * model
Definition: expat.h:151
int fGSLMinEval
GSL min evaluations. Ignored by some integrators.
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)