GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
QELXSec.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/Integrator.h>
14 
17 #include "Framework/Conventions/GBuild.h"
24 
32 #include "Framework/Utils/Range1.h"
34 
35 using namespace genie;
36 using namespace genie::constants;
37 
38 //____________________________________________________________________________
40 XSecIntegratorI("genie::QELXSec")
41 {
42 
43 }
44 //____________________________________________________________________________
45 QELXSec::QELXSec(string config) :
46 XSecIntegratorI("genie::QELXSec", config)
47 {
48 
49 }
50 //____________________________________________________________________________
52 {
53 
54 }
55 //____________________________________________________________________________
57  const XSecAlgorithmI * model, const Interaction * in) const
58 {
59  LOG("QELXSec",pDEBUG) << "Beginning integrate";
60  if(! model->ValidProcess(in)) return 0.;
61 
62  const KPhaseSpace & kps = in->PhaseSpace();
63  if(!kps.IsAboveThreshold()) {
64  LOG("QELXSec", pDEBUG) << "*** Below energy threshold";
65  return 0;
66  }
67  Range1D_t rQ2 = kps.Limits(kKVQ2);
68  if(rQ2.min<0 || rQ2.max<0) return 0;
69  LOG("QELXSec", pDEBUG)
70  << "Q2 integration range = (" << rQ2.min << ", " << rQ2.max << ")";
71 
72  Interaction * interaction = new Interaction(*in);
73  interaction->SetBit(kISkipProcessChk);
74  interaction->SetBit(kISkipKinematicChk);
75 
76  ROOT::Math::IBaseFunctionOneDim * func = new
77  utils::gsl::dXSec_dQ2_E(model, interaction);
78  ROOT::Math::IntegrationOneDim::Type ig_type =
80 
81  double abstol = 0; //We mostly care about relative tolerance
82  ROOT::Math::Integrator ig(*func,ig_type,abstol,fGSLRelTol,fGSLMaxSizeOfSubintervals, fGSLRule);
83  double xsec = ig.Integral(rQ2.min, rQ2.max) * (1E-38 * units::cm2);
84 
85  //LOG("QELXSec", pDEBUG) << "XSec[QEL] (E = " << E << ") = " << xsec;
86 
87  delete func;
88  delete interaction;
89 
90  return xsec;
91 }
92 //____________________________________________________________________________
93 void QELXSec::Configure(const Registry & config)
94 {
95  Algorithm::Configure(config);
96  this->LoadConfig();
97 }
98 //____________________________________________________________________________
99 void QELXSec::Configure(string config)
100 {
101  Algorithm::Configure(config);
102  this->LoadConfig();
103 }
104 //____________________________________________________________________________
106 {
107  // Get GSL integration type & relative tolerance
108  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive"));
109  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 0.001 ) ;
110  int max_size_of_subintervals;
111  GetParamDef( "gsl-max-size-of-subintervals", max_size_of_subintervals, 40000);
112  fGSLMaxSizeOfSubintervals = (unsigned int) max_size_of_subintervals;
113  int rule;
114  GetParamDef( "gsl-rule", rule, 3);
115  fGSLRule = (unsigned int) rule;
116  if (fGSLRule>6) fGSLRule=3;
117 }
118 //____________________________________________________________________________
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition: GSLUtils.cxx:23
string fGSLIntgType
name of GSL numerical integrator
Cross Section Integrator Interface.
A simple [min,max] interval for doubles.
Definition: Range1.h:42
unsigned int fGSLRule
GSL Gauss-Kronrod integration rule (only for GSL 1D adaptive type)
void LoadConfig(void)
Definition: QELXSec.cxx:105
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 ~QELXSec()
Definition: QELXSec.cxx:51
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
Kinematical phase space.
Definition: KPhaseSpace.h:33
double func(double x, double y)
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual bool ValidProcess(const Interaction *i) const =0
Can this cross section algorithm handle the input process?
void Configure(const Registry &config)
Definition: QELXSec.cxx:93
double max
Definition: Range1.h:53
unsigned int fGSLMaxSizeOfSubintervals
GSL maximum number of sub-intervals for 1D integrator.
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 Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
Definition: QELXSec.cxx:56
double min
Definition: Range1.h:52
bool GetParamDef(const RgKey &name, T &p, const T &def) const
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)