DMELXSec.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: Joshua Berger <jberger \at physics.wisc.edu>
8  University of Wisconsin-Madison
9 
10  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
11  University of Liverpool & STFC Rutherford Appleton Lab
12 */
13 //____________________________________________________________________________
14 
15 #include <TMath.h>
16 #include <Math/IFunction.h>
17 #include <Math/Integrator.h>
18 
28 
36 #include "Framework/Utils/Range1.h"
38 
39 using namespace genie;
40 using namespace genie::constants;
41 using namespace genie::controls;
42 
43 //____________________________________________________________________________
45 XSecIntegratorI("genie::DMELXSec")
46 {
47 
48 }
49 //____________________________________________________________________________
51 XSecIntegratorI("genie::DMELXSec", config)
52 {
53 
54 }
55 //____________________________________________________________________________
57 {
58 
59 }
60 //____________________________________________________________________________
62  const XSecAlgorithmI * model, const Interaction * in) const
63 {
64  LOG("DMELXSec",pDEBUG) << "Beginning integrate";
65  if(! model->ValidProcess(in)) return 0.;
66 
67  const KPhaseSpace & kps = in->PhaseSpace();
68  if(!kps.IsAboveThreshold()) {
69  LOG("DMELXSec", pDEBUG) << "*** Below energy threshold";
70  return 0;
71  }
72  Range1D_t rQ2 = kps.Limits(kKVQ2);
73 
74  if(rQ2.min<0 || rQ2.max<0) return 0;
75  LOG("DMELXSec", pDEBUG)
76  << "Q2 integration range = (" << rQ2.min << ", " << rQ2.max << ")";
77 
78  Interaction * interaction = new Interaction(*in);
79  interaction->SetBit(kISkipProcessChk);
80  interaction->SetBit(kISkipKinematicChk);
81 
82  ROOT::Math::IBaseFunctionOneDim * func = new
83  utils::gsl::dXSec_dQ2_E(model, interaction);
86 
87  double abstol = 1; //We mostly care about relative tolerance
88  ROOT::Math::Integrator ig(*func,ig_type,abstol,fGSLRelTol,fGSLMaxEval);
89  double xsec = ig.Integral(rQ2.min, rQ2.max) * (1E-38 * units::cm2);
90 
91  //LOG("DMELXSec", pDEBUG) << "XSec[DMEL] (E = " << E << ") = " << xsec;
92 
93  delete func;
94  delete interaction;
95 
96  return xsec;
97 }
98 //____________________________________________________________________________
100 {
101  Algorithm::Configure(config);
102  this->LoadConfig();
103 }
104 //____________________________________________________________________________
106 {
107  Algorithm::Configure(config);
108  this->LoadConfig();
109 }
110 //____________________________________________________________________________
112 {
113  // Get GSL integration type & relative tolerance
114  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive"));
115  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 0.01 ) ;
116  int max;
117  GetParamDef( "gsl-max-eval", max, 100000) ;
118  fGSLMaxEval = (unsigned int) max ;
119 }
120 //____________________________________________________________________________
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
virtual ~DMELXSec()
Definition: DMELXSec.cxx:56
void LoadConfig(void)
Definition: DMELXSec.cxx:111
Definition: config.py:1
static const double cm2
Definition: Units.h:77
void Configure(const Registry &config)
Definition: DMELXSec.cxx:99
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)
Misc GENIE control constants.
Double_t xsec[nknots]
Definition: testXsec.C:47
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
double max
Definition: Range1.h:54
int fGSLMaxEval
GSL max evaluations.
ifstream in
Definition: comparison.C:7
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
Definition: DMELXSec.cxx:61
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
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)