GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
COHXSec.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/IntegratorMultiDim.h>
14 #include "Math/AdaptiveIntegratorMultiDim.h"
15 
17 #include "Framework/Conventions/GBuild.h"
26 #include "Framework/Utils/Range1.h"
28 
29 using namespace genie;
30 using namespace genie::constants;
31 using namespace genie::utils;
32 
33 //____________________________________________________________________________
35 XSecIntegratorI("genie::COHXSec")
36 {
37 
38 }
39 //____________________________________________________________________________
40 COHXSec::COHXSec(string config) :
41 XSecIntegratorI("genie::COHXSec", config)
42 {
43 
44 }
45 //____________________________________________________________________________
47 {
48 
49 }
50 //____________________________________________________________________________
52  const XSecAlgorithmI * model, const Interaction * in) const
53 {
54  if(! model->ValidProcess(in) ) return 0.;
55 
56  const KPhaseSpace & kps = in->PhaseSpace();
57  if(!kps.IsAboveThreshold()) {
58  LOG("COHXSec", pDEBUG) << "*** Below energy threshold";
59  return 0;
60  }
61  Range1D_t xl = kps.Limits(kKVx);
62  Range1D_t yl = kps.Limits(kKVy);
63  Range1D_t Q2l;
65  Q2l.max = fQ2Max;
66 
67  Interaction * interaction = new Interaction(*in);
68  interaction->SetBit(kISkipProcessChk);
69  //interaction->SetBit(kISkipKinematicChk);
70 
71  double xsec = 0.0;
72 
73  if (model->Id().Name() == "genie::ReinSehgalCOHPiPXSec") {
74  LOG("COHXSec", pINFO)
75  << "x integration range = [" << xl.min << ", " << xl.max << "]";
76  LOG("COHXSec", pINFO)
77  << "y integration range = [" << yl.min << ", " << yl.max << "]";
78 
79  ROOT::Math::IBaseFunctionMultiDim * func =
80  new utils::gsl::d2XSec_dxdy_E(model, interaction);
81  ROOT::Math::IntegrationMultiDim::Type ig_type =
83 
84  double abstol = 1; //We mostly care about relative tolerance.
85  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
86  if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) {
87  ROOT::Math::AdaptiveIntegratorMultiDim * cast =
88  dynamic_cast<ROOT::Math::AdaptiveIntegratorMultiDim*>( ig.GetIntegrator() );
89  assert(cast);
90  cast->SetMinPts(fGSLMinEval);
91  }
92 
93  double kine_min[2] = { xl.min, yl.min };
94  double kine_max[2] = { xl.max, yl.max };
95  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
96  delete func;
97  }
98  else if (model->Id().Name() == "genie::BergerSehgalCOHPiPXSec2015")
99  {
100  ROOT::Math::IBaseFunctionMultiDim * func =
101  new utils::gsl::d2XSec_dQ2dy_E(model, interaction);
102  ROOT::Math::IntegrationMultiDim::Type ig_type =
104  ROOT::Math::IntegratorMultiDim ig(ig_type);
105  ig.SetRelTolerance(fGSLRelTol);
106  ig.SetFunction(*func);
107  if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) {
108  ROOT::Math::AdaptiveIntegratorMultiDim * cast =
109  dynamic_cast<ROOT::Math::AdaptiveIntegratorMultiDim*>( ig.GetIntegrator() );
110  assert(cast);
111  cast->SetMinPts(fGSLMinEval);
112  }
113  double kine_min[2] = { Q2l.min, yl.min };
114  double kine_max[2] = { Q2l.max, yl.max };
115  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
116  delete func;
117  }
118  else if (model->Id().Name() == "genie::BergerSehgalFMCOHPiPXSec2015")
119  {
120  Range1D_t tl;
122  tl.max = fTMax;
123 
124  ROOT::Math::IBaseFunctionMultiDim * func =
125  new utils::gsl::d2XSec_dQ2dydt_E(model, interaction);
126  ROOT::Math::IntegrationMultiDim::Type ig_type =
128  ROOT::Math::IntegratorMultiDim ig(ig_type);
129  ig.SetRelTolerance(fGSLRelTol);
130  ig.SetFunction(*func);
131  if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) {
132  ROOT::Math::AdaptiveIntegratorMultiDim * cast =
133  dynamic_cast<ROOT::Math::AdaptiveIntegratorMultiDim*>( ig.GetIntegrator() );
134  assert(cast);
135  cast->SetMinPts(fGSLMinEval);
136  }
137  double kine_min[3] = { Q2l.min, yl.min, tl.min };
138  double kine_max[3] = { Q2l.max, yl.max, tl.max };
139  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
140  delete func;
141  }
142 
143  const InitialState & init_state = in->InitState();
144  double Ev = init_state.ProbeE(kRfLab);
145  LOG("COHXSec", pINFO) << "XSec[COH] (E = " << Ev << " GeV) = " << xsec;
146 
147  delete interaction;
148 
149  return xsec;
150 }
151 //____________________________________________________________________________
152 void COHXSec::Configure(const Registry & config)
153 {
154  Algorithm::Configure(config);
155  this->LoadConfig();
156 }
157 //____________________________________________________________________________
158 void COHXSec::Configure(string config)
159 {
160  Algorithm::Configure(config);
161  this->LoadConfig();
162 }
163 //____________________________________________________________________________
165 {
166 
167  // Get GSL integration type & relative tolerance
168  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
169  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-2 ) ;
170 
171  int max_eval, min_eval ;
172  GetParamDef( "gsl-max-eval", max_eval, 500000 ) ;
173  GetParamDef( "gsl-min-eval", min_eval, 5000 ) ;
174 
175  fGSLMaxEval = (unsigned int) max_eval ;
176  fGSLMinEval = (unsigned int) min_eval ;
177 
178  //-- COH model parameter t_max for t = (q - p_pi)^2
179  GetParam("COH-t-max", fTMax ) ;
180 
181  //-- COH model bounds of integration for Q^2
182  GetParam( "COH-Q2-min", fQ2Min ) ;
183  GetParam( "COH-Q2-max", fQ2Max ) ;
184 
185 }
186 //____________________________________________________________________________
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
string fGSLIntgType
name of GSL numerical integrator
void LoadConfig(void)
Definition: COHXSec.cxx:164
double fQ2Max
upper bound of integration for Q^2 in Berger-Sehgal Model
Definition: COHXSec.h:45
Cross Section Integrator Interface.
double fTMax
upper bound for t = (q - p_pi)^2
Definition: COHXSec.h:46
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:59
A simple [min,max] interval for doubles.
Definition: Range1.h:42
void Configure(const Registry &config)
Definition: COHXSec.cxx:152
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:96
string Name(void) const
Definition: AlgId.h:44
static constexpr double cm2
Definition: Units.h:69
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
Kinematical phase space.
Definition: KPhaseSpace.h:33
double func(double x, double y)
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
virtual bool ValidProcess(const Interaction *i) const =0
Can this cross section algorithm handle the input process?
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
Definition: COHXSec.cxx:51
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
double max
Definition: Range1.h:53
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
double min
Definition: Range1.h:52
virtual ~COHXSec()
Definition: COHXSec.cxx:46
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 fQ2Min
lower bound of integration for Q^2 in Berger-Sehgal Model
Definition: COHXSec.h:44
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)