DFRXSec.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 */
13 //____________________________________________________________________________
14 
15 #include <TMath.h>
16 #include <Math/IFunction.h>
17 #include <Math/IntegratorMultiDim.h>
18 #include "Math/AdaptiveIntegratorMultiDim.h"
19 
34 
36 
37 using namespace genie;
38 using namespace genie::utils;
39 
40 //____________________________________________________________________________
42  : XSecIntegratorI("genie::DFRXSec")
43 {
44 
45 }
46 
47 //____________________________________________________________________________
49  : XSecIntegratorI("genie::DFRXSec", config)
50 {
51 
52 }
53 
54 //____________________________________________________________________________
56 {
57 
58 }
59 
60 //____________________________________________________________________________
61 double DFRXSec::Integrate (const XSecAlgorithmI* model, const Interaction* in) const
62 {
63  if(! model->ValidProcess(in) ) return 0.;
64 
65  const KPhaseSpace & kps = in->PhaseSpace();
66  if(!kps.IsAboveThreshold()) {
67  LOG("DFRXSec", pDEBUG) << "*** Below energy threshold";
68  return 0;
69  }
70  Range1D_t xl = kps.Limits(kKVx);
71  Range1D_t yl = kps.Limits(kKVy);
72  // the actual t lower limit depends on Q^2 and nu, or, equivalently, x and y.
73  // defer to KPhaseSpace::IsAlllowed() on where to start the t integral.
74  Range1D_t tl;
76  tl.max = fTMax;
77 
78  Interaction * interaction = new Interaction(*in);
79  interaction->SetBit(kISkipProcessChk); // todo: was enabled in COH model. do I need it here?
80  //interaction->SetBit(kISkipKinematicChk);
81 
82  double xsec = 0.0;
83  LOG("DFRXSec", pINFO)
84  << "x integration range = [" << xl.min << ", " << xl.max << "]";
85  LOG("DFRXSec", pINFO)
86  << "y integration range = [" << yl.min << ", " << yl.max << "]";
87  LOG("DFRXSec", pINFO)
88  << "t integration range = [" << tl.min << ", " << tl.max << "]";
89 
90  ROOT::Math::IBaseFunctionMultiDim * func =
91  new utils::gsl::d3XSec_dxdydt_E(model, interaction);
94 
95  double abstol = 1; //We mostly care about relative tolerance.
96  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
97  if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) {
98  ROOT::Math::AdaptiveIntegratorMultiDim * cast =
99  dynamic_cast<ROOT::Math::AdaptiveIntegratorMultiDim*>( ig.GetIntegrator() );
100  assert(cast);
101  cast->SetMinPts(fGSLMinEval);
102  }
103 
104  double kine_min[3] = { xl.min, yl.min, tl.min };
105  double kine_max[3] = { xl.max, yl.max, tl.max };
106  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
107  delete func;
108  return xsec;
109 }
110 
111 //____________________________________________________________________________
113 {
114  Algorithm::Configure(config);
115  this->LoadConfig();
116 }
117 
118 //____________________________________________________________________________
120 {
121  Algorithm::Configure(config);
122  this->LoadConfig();
123 }
124 
125 //____________________________________________________________________________
127 {
128  // Get GSL integration type & relative tolerance
129  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
130  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-2 ) ;
131 
132  int max, min;
133  GetParamDef( "gsl-max-eval", max, 500000 ) ;
134  GetParamDef( "gsl-min-eval", min, 5000 ) ;
135  fGSLMaxEval = (unsigned int) max ;
136  fGSLMinEval = (unsigned int) min ;
137 
138  //-- DFR model parameter t_max for t = (q - p_pi)^2
139  GetParam( "DFR-t-max", fTMax ) ;
140 
141 }
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
string fGSLIntgType
name of GSL numerical integrator
void Configure(const Registry &config)
Definition: DFRXSec.cxx:112
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
double fTMax
upper bound for t = (q - p_pi)^2
Definition: DFRXSec.h:54
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
virtual ~DFRXSec()
Definition: DFRXSec.cxx:55
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 Integrate(const XSecAlgorithmI *model, const Interaction *i) const
Definition: DFRXSec.cxx:61
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:63
Double_t xsec[nknots]
Definition: testXsec.C:47
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
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?
assert(nhit_max >=nhit_nbins)
double min
Definition: Range1.h:53
void LoadConfig(void)
Definition: DFRXSec.cxx:126
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
const XML_Char XML_Content * model
Definition: expat.h:151
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
int fGSLMinEval
GSL min evaluations. Ignored by some integrators.
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)