IMDXSec.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/Integrator.h>
22 
31 
32 using namespace genie;
33 using namespace genie::constants;
34 
35 //____________________________________________________________________________
37 XSecIntegratorI("genie::IMDXSec")
38 {
39 
40 }
41 //____________________________________________________________________________
43 XSecIntegratorI("genie::IMDXSec", config)
44 {
45 
46 }
47 //____________________________________________________________________________
49 {
50 
51 }
52 //____________________________________________________________________________
54  const XSecAlgorithmI * model, const Interaction * in) const
55 {
56  if(! model->ValidProcess(in) ) return 0.;
57 
58  const KPhaseSpace & kps = in->PhaseSpace();
59  if(!kps.IsAboveThreshold()) {
60  LOG("IMDXSec", pDEBUG) << "*** below energy threshold";
61  return 0;
62  }
63  Range1D_t yl = kps.Limits(kKVy);
64 
65  Interaction * interaction = new Interaction(*in);
66  interaction->SetBit(kISkipProcessChk);
67  //interaction->SetBit(kISkipKinematicChk);
68 
69  double abstol = 1; // We mostly care about relative tolerance
70  ROOT::Math::IBaseFunctionOneDim * func = new utils::gsl::dXSec_dy_E(model, interaction);
72  ROOT::Math::Integrator ig(*func,ig_type,abstol,fGSLRelTol,fGSLMaxEval);
73  double xsec = ig.Integral(yl.min, yl.max) * (1E-38 * units::cm2);
74 
75  //LOG("IMDXSec", pDEBUG) << "XSec[IMD] (E = " << E << ") = " << xsec;
76 
77  delete interaction;
78  delete func;
79 
80  return xsec;
81 }
82 //____________________________________________________________________________
84 {
85  Algorithm::Configure(config);
86  this->LoadConfig();
87 }
88 //____________________________________________________________________________
90 {
91  Algorithm::Configure(config);
92  this->LoadConfig();
93 }
94 //____________________________________________________________________________
96 {
97  // Get GSL integration type & relative tolerance
98  GetParamDef( "gsl-integration-type", fGSLIntgType, string( "adaptive" ) ) ;
99  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-4 ) ;
100  int max;
101  GetParamDef( "gsl-max-eval", max, 100000 ) ;
102  fGSLMaxEval = (unsigned int) max ;
103 
104 }
105 //____________________________________________________________________________
106 
107 
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition: GSLUtils.cxx:31
string fGSLIntgType
name of GSL numerical integrator
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
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.
virtual ~IMDXSec()
Definition: IMDXSec.cxx:48
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
void Configure(const Registry &config)
Definition: IMDXSec.cxx:83
double func(double x, double y)
Double_t xsec[nknots]
Definition: testXsec.C:47
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
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 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
const XML_Char XML_Content * model
Definition: expat.h:151
void LoadConfig(void)
Definition: IMDXSec.cxx:95
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)
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
Definition: IMDXSec.cxx:53