GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HEDISXSec.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  or see $GENIE/LICENSE
6 
7  Author: Alfonso Garcia <alfonsog \at nikhef.nl>
8  NIKHEF
9 
10  For the class documentation see the corresponding header file.
11 
12 */
13 //____________________________________________________________________________
14 
21 #include "Framework/Utils/Range1.h"
24 #include "Framework/Utils/RunOpt.h"
27 
28 #include <TMath.h>
29 #include <Math/IFunction.h>
30 #include <Math/IntegratorMultiDim.h>
31 #include "Math/AdaptiveIntegratorMultiDim.h"
32 
33 using namespace genie;
34 using namespace genie::constants;
35 
36 //____________________________________________________________________________
38 XSecIntegratorI("genie::HEDISXSec")
39 {
40 
41 }
42 //____________________________________________________________________________
43 HEDISXSec::HEDISXSec(string config) :
44 XSecIntegratorI("genie::HEDISXSec", config)
45 {
46 
47 }
48 //____________________________________________________________________________
50 {
51 
52 }
53 //____________________________________________________________________________
55  const XSecAlgorithmI * model, const Interaction * in) const
56 {
57 
58  if(! model->ValidProcess(in) ) return 0.;
59 
60  const KPhaseSpace & kps = in->PhaseSpace();
61  if(!kps.IsAboveThreshold()) {
62  LOG("DISXSec", pDEBUG) << "*** Below energy threshold";
63  return 0;
64  }
65 
66  const InitialState & init_state = in->InitState();
67  double Ev = init_state.ProbeE(kRfLab);
68 
69  // If the input interaction is off a nuclear target, then chek whether
70  // the corresponding free nucleon cross section already exists at the
71  // cross section spline list.
72  // If yes, calculate the nuclear cross section based on that value.
73  //
75  if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() ) {
76 
77  int nucpdgc = init_state.Tgt().HitNucPdg();
78  int NNucl = (pdg::IsProton(nucpdgc)) ? init_state.Tgt().Z() : init_state.Tgt().N();
79 
80  Interaction * interaction = new Interaction(*in);
81  Target * target = interaction->InitStatePtr()->TgtPtr();
82  if(pdg::IsProton(nucpdgc)) { target->SetId(kPdgTgtFreeP); }
83  else { target->SetId(kPdgTgtFreeN); }
84  if(xsl->SplineExists(model,interaction)) {
85  const Spline * spl = xsl->GetSpline(model, interaction);
86  double xsec = spl->Evaluate(Ev);
87  LOG("HEDISXSec", pINFO) << "From XSecSplineList: XSec[HEDIS,free nucleon] (E = " << Ev << " GeV) = " << xsec;
88  if( !interaction->TestBit(kIAssumeFreeNucleon) ) {
89  xsec *= NNucl;
90  LOG("HEDISXSec", pINFO) << "XSec[HEDIS] (E = " << Ev << " GeV) = " << xsec;
91  }
92  delete interaction;
93  return xsec;
94  }
95  delete interaction;
96  }
97 
98  LOG("HEDISXSec", pINFO) << in->AsString();
99 
100  Range1D_t xl = kps.XLim();
101  Range1D_t Q2l = kps.Q2Lim();
102  LOG("HEDISXSec", pDEBUG) << "X only kinematic range = [" << xl.min << ", " << xl.max << "]";
103  LOG("HEDISXSec", pDEBUG) << "Q2 only kinematic range = [" << Q2l.min << ", " << Q2l.max << "]";
104 
105  if (xl.min < fSFXmin) xl.min=fSFXmin;
106  if (Q2l.min < fSFQ2min) Q2l.min=fSFQ2min;
107  if (Q2l.max > fSFQ2max) Q2l.max=fSFQ2max;
108 
109  LOG("HEDISXSec", pDEBUG) << "X kinematic+PDF range = [" << xl.min << ", " << xl.max << "]";
110  LOG("HEDISXSec", pDEBUG) << "Q2 kinematic+PDF range = [" << Q2l.min << ", " << Q2l.max << "]";
111 
112 
113  bool phsp_ok =
114  (Q2l.min >= 0. && Q2l.max >= 0. && Q2l.max >= Q2l.min &&
115  xl.min >= 0. && xl.max <= 1. && xl.max >= xl.min);
116 
117  if (!phsp_ok) return 0.;
118 
119 
120  // Just go ahead and integrate the input differential cross section for the
121  // specified interaction.
122  //
123  double xsec = 0.;
124 
125  Interaction * interaction = new Interaction(*in);
126 
127  // If a GSL option has been chosen, then the total xsec is recomptued
128  ROOT::Math::IBaseFunctionMultiDim * func = new utils::gsl::d2XSec_dlog10xdlog10Q2_E(model, interaction);
129  ROOT::Math::IntegrationMultiDim::Type ig_type = utils::gsl::IntegrationNDimTypeFromString(fGSLIntgType);
130  double abstol = 0; //In vegas we care about number of interactions
131  double reltol = 0; //In vegas we care about number of interactions
132  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, reltol, fGSLMaxEval);
133  double kine_min[2] = { TMath::Log10(xl.min), TMath::Log10(Q2l.min) };
134  double kine_max[2] = {TMath::Log10(xl.max), TMath::Log10(Q2l.max) };
135  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
136  delete func;
137 
138  delete interaction;
139 
140  LOG("HEDISXSec", pINFO) << "XSec[HEDIS] (E = " << Ev << " GeV) = " << xsec * (1E+38/units::cm2) << " x 1E-38 cm^2";
141 
142  return xsec;
143 
144 }
145 //____________________________________________________________________________
146 void HEDISXSec::Configure(const Registry & config)
147 {
148  Algorithm::Configure(config);
149  this->LoadConfig();
150 }
151 //____________________________________________________________________________
152 void HEDISXSec::Configure(string config)
153 {
154  Algorithm::Configure(config);
155  this->LoadConfig();
156 }
157 //____________________________________________________________________________
159 {
160 
161  // Get GSL integration type & relative tolerance
162  GetParamDef( "gsl-integration-type", fGSLIntgType, string("vegas") ) ;
163 
164  int max_eval ;
165  GetParamDef( "gsl-max-eval", max_eval, 500000 ) ;
166  fGSLMaxEval = (unsigned int) max_eval ;
167 
168  // Limits from the SF tables that are useful to reduce computation
169  // time of the total cross section
170  GetParam("XGrid-Min", fSFXmin ) ;
171  GetParam("Q2Grid-Min", fSFQ2min ) ;
172  GetParam("Q2Grid-Max", fSFQ2max ) ;
173 
174 }
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
string fGSLIntgType
name of GSL numerical integrator
Cross Section Integrator Interface.
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:59
int HitNucPdg(void) const
Definition: Target.cxx:304
A simple [min,max] interval for doubles.
Definition: Range1.h:42
bool SplineExists(const XSecAlgorithmI *alg, const Interaction *i) const
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:58
bool IsNucleus(void) const
Definition: Target.cxx:272
double fSFQ2max
maximum value of Q2 for which SF tables are computed
Definition: HEDISXSec.h:49
static XSecSplineList * Instance()
Range1D_t Q2Lim(void) const
Q2 limits.
double Evaluate(double x) const
Definition: Spline.cxx:363
void SetId(int pdgc)
Definition: Target.cxx:149
string AsString(void) const
Summary information for an interaction.
Definition: Interaction.h:56
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool IsEmpty(void) const
static constexpr double cm2
Definition: Units.h:69
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
Definition: HEDISXSec.cxx:54
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
Kinematical phase space.
Definition: KPhaseSpace.h:33
double fSFQ2min
minimum value of Q2 for which SF tables are computed
Definition: HEDISXSec.h:48
const int kPdgTgtFreeN
Definition: PDGCodes.h:200
const int kPdgTgtFreeP
Definition: PDGCodes.h:199
void LoadConfig(void)
Definition: HEDISXSec.cxx:158
double func(double x, double y)
int Z(void) const
Definition: Target.h:68
#define pINFO
Definition: Messenger.h:62
virtual bool ValidProcess(const Interaction *i) const =0
Can this cross section algorithm handle the input process?
void Configure(const Registry &config)
Definition: HEDISXSec.cxx:146
double max
Definition: Range1.h:53
int fGSLMaxEval
GSL max evaluations.
int N(void) const
Definition: Target.h:69
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
Target * TgtPtr(void) const
Definition: InitialState.h:67
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
double fSFXmin
minimum value of x for which SF tables are computed
Definition: HEDISXSec.h:47
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
Range1D_t XLim(void) const
x limits
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:52
virtual ~HEDISXSec()
Definition: HEDISXSec.cxx:49
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
const Spline * GetSpline(const XSecAlgorithmI *alg, const Interaction *i) const
List of cross section vs energy splines.
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63