GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MECXSec.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 */
7 //____________________________________________________________________________
8 
9 #include <TMath.h>
10 #include <Math/IFunction.h>
11 #include <Math/IntegratorMultiDim.h>
12 #include "Math/AdaptiveIntegratorMultiDim.h"
13 
15 #include "Framework/Conventions/GBuild.h"
31 #include "Framework/Utils/Range1.h"
34 
35 using namespace genie;
36 using namespace genie::constants;
37 using namespace genie::controls;
38 using namespace genie::utils;
39 
40 //____________________________________________________________________________
42 XSecIntegratorI("genie::MECXSec")
43 {
44 
45 }
46 //____________________________________________________________________________
47 MECXSec::MECXSec(string config) :
48 XSecIntegratorI("genie::MECXSec", config)
49 {
50 
51 }
52 //____________________________________________________________________________
54 {
55 
56 }
57 //____________________________________________________________________________
59  const XSecAlgorithmI * model, const Interaction * in) const
60 {
61  if(! model->ValidProcess(in) ) return 0.;
62 
63  const KPhaseSpace & kps = in->PhaseSpace(); // only OK phase space for this
64  if(!kps.IsAboveThreshold()) {
65  LOG("MECXSec", pDEBUG) << "*** Below energy threshold";
66  return 0;
67  }
68 
69  Interaction interaction(*in);
70  interaction.SetBit(kISkipProcessChk);
71  interaction.SetBit(kISkipKinematicChk);
72 
73  // T, costh limits
74  double Enu = in->InitState().ProbeE(kRfLab);
75  double LepMass = in->FSPrimLepton()->Mass();
76  double TMax = Enu - LepMass;
77  double TMin = 0.0;
78  double CosthMax = 1.0;
79  double CosthMin = -1.0;
80  if (Enu < fQ3Max) {
81  TMin = 0 ;
82  CosthMin = -1 ;
83  } else {
84  TMin = TMath::Sqrt(TMath::Power(LepMass, 2) + TMath::Power((Enu - fQ3Max), 2)) - LepMass;
85  CosthMin = TMath::Sqrt(1 - TMath::Power((fQ3Max / Enu ), 2));
86  }
87 
88  double kine_min[2] = { TMin, CosthMin };
89  double kine_max[2] = { TMax, CosthMax };
90 
91  double xsec = 0;
92 
93  double abstol = 1; //We mostly care about relative tolerance.
94  genie::utils::mec::gsl::d2Xsec_dTCosth func(model, interaction, Enu, LepMass );
95  ROOT::Math::IntegrationMultiDim::Type ig_type =
97  ROOT::Math::IntegratorMultiDim ig(func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
98 
99  xsec = ig.Integral(kine_min, kine_max);
100 
101  return xsec;
102 }
103 //____________________________________________________________________________
104 void MECXSec::Configure(const Registry & config)
105 {
106  Algorithm::Configure(config);
107  this->LoadConfig();
108 }
109 //____________________________________________________________________________
110 void MECXSec::Configure(string config)
111 {
112  Algorithm::Configure(config);
113  this->LoadConfig();
114 }
115 //____________________________________________________________________________
117 {
118  GetParam( "NSV-Q3Max", fQ3Max ) ;
119 
120  // Get GSL integration type & relative tolerance
121  GetParamDef( "gsl-integration-type", fGSLIntgType, string("vegas") ) ;
122 
123  int max ;
124  GetParamDef( "gsl-max-eval", max, 20000 ) ;
125  fGSLMaxEval = (unsigned int) max ;
126 
127  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 0.01 ) ;
128  GetParamDef( "split-integral", fSplitIntegral, true ) ;
129 
130 }
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
virtual ~MECXSec()
Definition: MECXSec.cxx:53
string fGSLIntgType
name of GSL numerical integrator
Cross Section Integrator Interface.
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:59
bool fSplitIntegral
Definition: MECXSec.h:49
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
Definition: MECXSec.cxx:58
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
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
Kinematical phase space.
Definition: KPhaseSpace.h:33
void Configure(const Registry &config)
Definition: MECXSec.cxx:104
double fQ3Max
Definition: MECXSec.h:53
void LoadConfig(void)
Definition: MECXSec.cxx:116
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?
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
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.
const InitialState & InitState(void) const
Definition: Interaction.h:69
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
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)