GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
IMDXSec.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 
15 #include "Framework/Conventions/GBuild.h"
23 
24 using namespace genie;
25 using namespace genie::constants;
26 
27 //____________________________________________________________________________
29 XSecIntegratorI("genie::IMDXSec")
30 {
31 
32 }
33 //____________________________________________________________________________
34 IMDXSec::IMDXSec(string config) :
35 XSecIntegratorI("genie::IMDXSec", config)
36 {
37 
38 }
39 //____________________________________________________________________________
41 {
42 
43 }
44 //____________________________________________________________________________
46  const XSecAlgorithmI * model, const Interaction * in) const
47 {
48  if(! model->ValidProcess(in) ) return 0.;
49 
50  const KPhaseSpace & kps = in->PhaseSpace();
51  if(!kps.IsAboveThreshold()) {
52  LOG("IMDXSec", pDEBUG) << "*** below energy threshold";
53  return 0;
54  }
55  Range1D_t yl = kps.Limits(kKVy);
56 
57  Interaction * interaction = new Interaction(*in);
58  interaction->SetBit(kISkipProcessChk);
59  //interaction->SetBit(kISkipKinematicChk);
60 
61  double abstol = 1; // We mostly care about relative tolerance
62  ROOT::Math::IBaseFunctionOneDim * func = new utils::gsl::dXSec_dy_E(model, interaction);
63  ROOT::Math::IntegrationOneDim::Type ig_type = utils::gsl::Integration1DimTypeFromString(fGSLIntgType);
64  ROOT::Math::Integrator ig(*func,ig_type,abstol,fGSLRelTol,fGSLMaxEval);
65  double xsec = ig.Integral(yl.min, yl.max) * (1E-38 * units::cm2);
66 
67  //LOG("IMDXSec", pDEBUG) << "XSec[IMD] (E = " << E << ") = " << xsec;
68 
69  delete interaction;
70  delete func;
71 
72  return xsec;
73 }
74 //____________________________________________________________________________
75 void IMDXSec::Configure(const Registry & config)
76 {
77  Algorithm::Configure(config);
78  this->LoadConfig();
79 }
80 //____________________________________________________________________________
81 void IMDXSec::Configure(string config)
82 {
83  Algorithm::Configure(config);
84  this->LoadConfig();
85 }
86 //____________________________________________________________________________
88 {
89  // Get GSL integration type & relative tolerance
90  GetParamDef( "gsl-integration-type", fGSLIntgType, string( "adaptive" ) ) ;
91  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-4 ) ;
92  int max;
93  GetParamDef( "gsl-max-eval", max, 100000 ) ;
94  fGSLMaxEval = (unsigned int) max ;
95 
96 }
97 //____________________________________________________________________________
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
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
virtual ~IMDXSec()
Definition: IMDXSec.cxx:40
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
void Configure(const Registry &config)
Definition: IMDXSec.cxx:75
double func(double x, double y)
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
bool GetParamDef(const RgKey &name, T &p, const T &def) const
void LoadConfig(void)
Definition: IMDXSec.cxx:87
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)
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
Definition: IMDXSec.cxx:45