GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
RESXSec.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 #include <Math/IFunction.h>
13 #include <Math/IntegratorMultiDim.h>
14 
15 #include "Framework/Conventions/GBuild.h"
26 
27 using namespace genie;
28 using namespace genie::constants;
29 
30 //____________________________________________________________________________
32 XSecIntegratorI("genie::RESXSec")
33 {
34 
35 }
36 //____________________________________________________________________________
37 RESXSec::RESXSec(string config) :
38 XSecIntegratorI("genie::RESXSec", config)
39 {
40 
41 }
42 //____________________________________________________________________________
44 {
45 
46 }
47 //____________________________________________________________________________
49  const XSecAlgorithmI * model, const Interaction * in) const
50 {
51  if(! model->ValidProcess(in) ) return 0.;
52 
53  const KPhaseSpace & kps = in->PhaseSpace();
54  if(!kps.IsAboveThreshold()) {
55  LOG("COHXSec", pDEBUG) << "*** Below energy threshold";
56  return 0;
57  }
58  Range1D_t Q2l = kps.Limits(kKVQ2);
59  Range1D_t Wl = kps.Limits(kKVW);
60 
61  Interaction * interaction = new Interaction(*in);
62  interaction->SetBit(kISkipProcessChk);
63  //interaction->SetBit(kISkipKinematicChk);
64 
65  ROOT::Math::IBaseFunctionMultiDim * func =
66  new utils::gsl::d2XSec_dWdQ2_E(model, interaction);
67 
68  ROOT::Math::IntegrationMultiDim::Type ig_type =
70  double abstol = 1E-16; //We mostly care about relative tolerance.
71 
72  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
73 
74  double kine_min[2] = { Wl.min, Q2l.min };
75  double kine_max[2] = { Wl.max, Q2l.max };
76  double xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
77 
78  LOG("RESXSec", pERROR) << "Integrator opt / Integrator = " << ig.Options().Integrator();
79 
80  if(xsec < 0) {
81  LOG("RESXSec", pERROR) << "Algorithm " << *model << " returns a negative cross-section (xsec = " << xsec << " 1E-38 * cm2)";
82  LOG("RESXSec", pERROR) << "for process" << *interaction;
83  LOG("RESXSec", pERROR) << "Integrator status code = " << ig.Status();
84  LOG("RESXSec", pERROR) << "Integrator error code = " << ig.Error();
85  }
86 
87  //LOG("RESXSec", pINFO) << "XSec[RES] (Ev = " << Ev << " GeV) = " << xsec;
88 
89  delete interaction;
90  delete func;
91  return xsec;
92 }
93 //____________________________________________________________________________
94 void RESXSec::Configure(const Registry & config)
95 {
96  Algorithm::Configure(config);
97  this->LoadConfig();
98 }
99 //____________________________________________________________________________
100 void RESXSec::Configure(string config)
101 {
102  Algorithm::Configure(config);
103  this->LoadConfig();
104 }
105 //____________________________________________________________________________
107 {
108  // Get GSL integration type & relative tolerance
109  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
110  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-2 ) ;
111  int max, min ;
112  GetParamDef( "gsl-max-eval", max, 500000 ) ;
113  GetParamDef( "gsl-min-eval", min, 5000 ) ;
114  fGSLMaxEval = (unsigned int) max ;
115  fGSLMinEval = (unsigned int) min ;
116 }
117 //____________________________________________________________________________
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
string fGSLIntgType
name of GSL numerical integrator
#define pERROR
Definition: Messenger.h:59
Cross Section Integrator Interface.
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:59
A simple [min,max] interval for doubles.
Definition: Range1.h:42
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:96
static constexpr double cm2
Definition: Units.h:69
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
Kinematical phase space.
Definition: KPhaseSpace.h:33
double func(double x, double y)
virtual bool ValidProcess(const Interaction *i) const =0
Can this cross section algorithm handle the input process?
void Configure(const Registry &config)
Definition: RESXSec.cxx:94
double max
Definition: Range1.h:53
int fGSLMaxEval
GSL max evaluations.
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
virtual ~RESXSec()
Definition: RESXSec.cxx:43
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
Definition: RESXSec.cxx:48
double min
Definition: Range1.h:52
bool GetParamDef(const RgKey &name, T &p, const T &def) const
void LoadConfig(void)
Definition: RESXSec.cxx:106
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:63
double fGSLRelTol
required relative tolerance (error)