GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CEvNSXSec.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 
22 
23 using namespace genie;
24 using namespace genie::constants;
25 using namespace genie::utils;
26 
27 //____________________________________________________________________________
29 XSecIntegratorI("genie::CEvNSXSec")
30 {
31 
32 }
33 //____________________________________________________________________________
34 CEvNSXSec::CEvNSXSec(string config) :
35 XSecIntegratorI("genie::CEvNSXSec", 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("CEvNS", pDEBUG) << "*** Below energy threshold";
53  return 0;
54  }
55 
56  Range1D_t Q2 = kps.Q2Lim();
57  LOG("CEvNS", pNOTICE)
58  << "Q2 integration range = [" << Q2.min << ", " << Q2.max << "] GeV^2";
59  assert(Q2.min > 0. && Q2.min < Q2.max);
60 
61  Interaction * interaction = new Interaction(*in);
62  interaction->SetBit(kISkipProcessChk);
63  interaction->SetBit(kISkipKinematicChk);
64 
65  ROOT::Math::IntegrationOneDim::Type ig_type =
67  ROOT::Math::IBaseFunctionOneDim * func =
68  new utils::gsl::dXSec_dQ2_E(model, interaction);
69  double abstol = 1; // We mostly care about relative tolerance
70  ROOT::Math::Integrator ig(*func,ig_type,abstol,fGSLRelTol,fGSLMaxEval);
71  double xsec = ig.Integral(Q2.min, Q2.max) * (1E-38 * units::cm2); // units: GeV^-2
72 
73  delete func;
74 
75  const InitialState & init_state = in->InitState();
76  double Ev = init_state.ProbeE(kRfLab);
77  LOG("CEvNS", pINFO)
78  << "XSec[CEvNS] (E = " << Ev << " GeV) = " << xsec/(units::cm2) << " cm^2";
79 
80  delete interaction;
81 
82  return xsec;
83 }
84 //____________________________________________________________________________
85 void CEvNSXSec::Configure(const Registry & config)
86 {
87  Algorithm::Configure(config);
88  this->LoadConfig();
89 }
90 //____________________________________________________________________________
91 void CEvNSXSec::Configure(string config)
92 {
93  Algorithm::Configure(config);
94  this->LoadConfig();
95 }
96 //____________________________________________________________________________
98 {
99  // Get GSL integration type & relative tolerance
100  this->GetParamDef("gsl-integration-type", fGSLIntgType, string("adaptive"));
101  this->GetParamDef("gsl-relative-tolerance", fGSLRelTol, 1E-2);
102 
103  // max and minimum number of integrand evaluations
104  int max_eval, min_eval ;
105  this->GetParamDef("gsl-max-eval", max_eval, 500000);
106  this->GetParamDef("gsl-min-eval", min_eval, 5000 );
107  fGSLMaxEval = (unsigned int) max_eval ;
108  fGSLMinEval = (unsigned int) min_eval ;
109 
110  // LOG("CEvNS", pDEBUG) << *this;
111 }
112 //____________________________________________________________________________
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
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
Definition: CEvNSXSec.cxx:45
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Range1D_t Q2Lim(void) const
Q2 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 ~CEvNSXSec()
Definition: CEvNSXSec.cxx:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
Kinematical phase space.
Definition: KPhaseSpace.h:33
void LoadConfig(void)
Definition: CEvNSXSec.cxx:97
double func(double x, double y)
#define pINFO
Definition: Messenger.h:62
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?
double max
Definition: Range1.h:53
int fGSLMaxEval
GSL max evaluations.
void Configure(const Registry &config)
Definition: CEvNSXSec.cxx:85
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.
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:52
#define pNOTICE
Definition: Messenger.h:61
bool GetParamDef(const RgKey &name, T &p, const T &def) const
double ProbeE(RefFrame_t rf) const
int fGSLMinEval
GSL min evaluations. Ignored by some integrators.
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
double fGSLRelTol
required relative tolerance (error)