COHXSecAR.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  STFC, Rutherford Appleton Laboratory
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/Integrator.h>
18 #include <Math/IntegratorMultiDim.h>
19 #include "Math/AdaptiveIntegratorMultiDim.h"
20 
30 #include "Framework/Utils/Range1.h"
32 
33 using namespace genie;
34 using namespace genie::constants;
35 using namespace genie::controls;
36 using namespace genie::utils;
37 
38 //____________________________________________________________________________
40 XSecIntegratorI("genie::COHXSecAR")
41 {
42 
43 }
44 //____________________________________________________________________________
46 XSecIntegratorI("genie::COHXSecAR", config)
47 {
48 
49 }
50 //____________________________________________________________________________
52 {
53 
54 }
55 //____________________________________________________________________________
57  const XSecAlgorithmI * model, const Interaction * in) const
58 {
59  const InitialState & init_state = in -> InitState();
60 
61  if(! model->ValidProcess(in) ) return 0.;
62 
63  const KPhaseSpace & kps = in->PhaseSpace();
64  if(!kps.IsAboveThreshold()) {
65  LOG("COHXSecAR", pDEBUG) << "*** Below energy threshold";
66  return 0;
67  }
68 
69  Range1D_t y_lim = kps.Limits(kKVy);
70 
71  // Check this
72  double Enu = init_state.ProbeE(kRfLab);
73  double Elep_min = (1.-y_lim.max) * Enu;
74  double Elep_max = (1.-y_lim.min) * Enu;
75 
76  LOG("COHXSecAR", pINFO)
77  << "Lepton energy integration range = [" << Elep_min << ", " << Elep_max << "]";
78 
79  Interaction * interaction = new Interaction(*in);
80  interaction->SetBit(kISkipProcessChk);
81  //interaction->SetBit(kISkipKinematicChk);
82 
83  double xsec = 0;
84  if (fSplitIntegral) {
87 
88  //~ ROOT::Math::IntegrationOneDim::Type ig_type = ROOT::Math::IntegrationOneDim::kNONADAPTIVE;
89  ROOT::Math::IntegrationOneDim::Type ig_type = ROOT::Math::IntegrationOneDim::kADAPTIVE;
90 
91  double abstol = 1; // Pretty sure this parameter is unused by ROOT.
92  int size = 1000; // Max number of subintervals, won't reach nearly this.
93  int rule = 2; // See https://www.gnu.org/software/gsl/manual/gsl-ref_17.html#SEC283
94  // Rule 2 is 21 points min
95  ROOT::Math::Integrator ig(*func,ig_type,abstol,fGSLRelTol,size,rule);
96 
97  xsec = ig.Integral(Elep_min, Elep_max) * (1E-38 * units::cm2);
98  delete func;
99  }
100  else {
101  double zero = kASmallNum;
102  double pi = kPi-kASmallNum ;
103  double twopi = 2*kPi-kASmallNum ;
104 
105  //~ ROOT::Math::IBaseFunctionMultiDim * func =
106  //~ new utils::gsl::wrap::d5Xsec_dEldOmegaldOmegapi(model, interaction);
107  //~ double kine_min[5] = { Elep_min, zero , zero , zero, zero };
108  //~ double kine_max[5] = { Elep_max, pi , twopi , pi , twopi};
109 
110  ROOT::Math::IBaseFunctionMultiDim * func =
111  new utils::gsl::d4Xsec_dEldThetaldOmegapi(model, interaction);
112  double kine_min[4] = { Elep_min, zero , zero , zero };
113  double kine_max[4] = { Elep_max, pi , pi , twopi };
114 
117 
118  double abstol = 1; //We mostly care about relative tolerance.
119  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
120 
121  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
122  delete func;
123  }
124 
125  delete interaction;
126 
127  return xsec;
128 }
129 //____________________________________________________________________________
131 {
132  Algorithm::Configure(config);
133  this->LoadConfig();
134 }
135 //____________________________________________________________________________
137 {
138  Algorithm::Configure(config);
139  this->LoadConfig();
140 }
141 //____________________________________________________________________________
143 {
144  // Get GSL integration type & relative tolerance
145  GetParamDef( "gsl-integration-type", fGSLIntgType, string("vegas") ) ;
146 
147  int max_eval;
148  GetParamDef( "gsl-max-eval", max_eval, 4000 ) ;
149  fGSLMaxEval = (unsigned int) max_eval ;
150 
151  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 0.01) ;
152  GetParamDef( "split-integral", fSplitIntegral, true ) ;
153 
154 }
155 //_____________________________________________________________________________
156 
157 
Cross Section Calculation Interface.
const double kPi
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
string fGSLIntgType
name of GSL numerical integrator
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
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
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
Definition: COHXSecAR.cxx:56
Definition: config.py:1
static const double cm2
Definition: Units.h:77
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
void Configure(const Registry &config)
Definition: COHXSecAR.cxx:130
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)
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:63
bool fSplitIntegral
Definition: COHXSecAR.h:42
Misc GENIE control constants.
Double_t xsec[nknots]
Definition: testXsec.C:47
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
void LoadConfig(void)
Definition: COHXSecAR.cxx:142
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
virtual ~COHXSecAR()
Definition: COHXSecAR.cxx:51
const XML_Char XML_Content * model
Definition: expat.h:151
auto zero()
Definition: PMNS.cxx:47
double ProbeE(RefFrame_t rf) const
Root of GENIE utility namespaces.
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
double fGSLRelTol
required relative tolerance (error)