GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DFRXSec.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 #include "Math/AdaptiveIntegratorMultiDim.h"
15 
18 #include "Framework/Conventions/GBuild.h"
30 
32 
33 using namespace genie;
34 using namespace genie::utils;
35 
36 //____________________________________________________________________________
38  : XSecIntegratorI("genie::DFRXSec")
39 {
40 
41 }
42 
43 //____________________________________________________________________________
44 DFRXSec::DFRXSec (std::string config)
45  : XSecIntegratorI("genie::DFRXSec", config)
46 {
47 
48 }
49 
50 //____________________________________________________________________________
52 {
53 
54 }
55 
56 //____________________________________________________________________________
57 double DFRXSec::Integrate (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("DFRXSec", pDEBUG) << "*** Below energy threshold";
64  return 0;
65  }
66  Range1D_t xl = kps.Limits(kKVx);
67  Range1D_t yl = kps.Limits(kKVy);
68  // the actual t lower limit depends on Q^2 and nu, or, equivalently, x and y.
69  // defer to KPhaseSpace::IsAlllowed() on where to start the t integral.
70  Range1D_t tl;
72  tl.max = fTMax;
73 
74  Interaction * interaction = new Interaction(*in);
75  interaction->SetBit(kISkipProcessChk); // todo: was enabled in COH model. do I need it here?
76  //interaction->SetBit(kISkipKinematicChk);
77 
78  double xsec = 0.0;
79  LOG("DFRXSec", pINFO)
80  << "x integration range = [" << xl.min << ", " << xl.max << "]";
81  LOG("DFRXSec", pINFO)
82  << "y integration range = [" << yl.min << ", " << yl.max << "]";
83  LOG("DFRXSec", pINFO)
84  << "t integration range = [" << tl.min << ", " << tl.max << "]";
85 
86  ROOT::Math::IBaseFunctionMultiDim * func =
87  new utils::gsl::d3XSec_dxdydt_E(model, interaction);
88  ROOT::Math::IntegrationMultiDim::Type ig_type =
90 
91  double abstol = 1; //We mostly care about relative tolerance.
92  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
93  if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) {
94  ROOT::Math::AdaptiveIntegratorMultiDim * cast =
95  dynamic_cast<ROOT::Math::AdaptiveIntegratorMultiDim*>( ig.GetIntegrator() );
96  assert(cast);
97  cast->SetMinPts(fGSLMinEval);
98  }
99 
100  double kine_min[3] = { xl.min, yl.min, tl.min };
101  double kine_max[3] = { xl.max, yl.max, tl.max };
102  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
103  delete func;
104  return xsec;
105 }
106 
107 //____________________________________________________________________________
108 void DFRXSec::Configure (const Registry& config)
109 {
110  Algorithm::Configure(config);
111  this->LoadConfig();
112 }
113 
114 //____________________________________________________________________________
115 void DFRXSec::Configure (std::string config)
116 {
117  Algorithm::Configure(config);
118  this->LoadConfig();
119 }
120 
121 //____________________________________________________________________________
123 {
124  // Get GSL integration type & relative tolerance
125  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
126  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-2 ) ;
127 
128  int max, min;
129  GetParamDef( "gsl-max-eval", max, 500000 ) ;
130  GetParamDef( "gsl-min-eval", min, 5000 ) ;
131  fGSLMaxEval = (unsigned int) max ;
132  fGSLMinEval = (unsigned int) min ;
133 
134  //-- DFR model parameter t_max for t = (q - p_pi)^2
135  GetParam( "DFR-t-max", fTMax ) ;
136 
137 }
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:108
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:59
A simple [min,max] interval for doubles.
Definition: Range1.h:42
virtual ~DFRXSec()
Definition: DFRXSec.cxx:51
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)
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
Definition: DFRXSec.cxx:57
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
virtual bool ValidProcess(const Interaction *i) const =0
Can this cross section algorithm handle the input process?
double max
Definition: Range1.h:53
int fGSLMaxEval
GSL max evaluations.
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
double min
Definition: Range1.h:52
void LoadConfig(void)
Definition: DFRXSec.cxx:122
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
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.
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)