COHXSec.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  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Mar 03, 2009 - CA
14  Renamed COHPiXSec -> COHXSec. Adapt to naming changes made to the coherent
15  generator for including coherent vector meson production.
16  @ Sep 07, 2009 - CA
17  Integrated with GNU Numerical Library (GSL) via ROOT's MathMore library.
18 
19 */
20 //____________________________________________________________________________
21 
22 #include <TMath.h>
23 #include <Math/IFunction.h>
24 #include <Math/IntegratorMultiDim.h>
25 #include "Math/AdaptiveIntegratorMultiDim.h"
26 
37 #include "Framework/Utils/Range1.h"
39 
40 using namespace genie;
41 using namespace genie::constants;
42 using namespace genie::utils;
43 
44 //____________________________________________________________________________
46 XSecIntegratorI("genie::COHXSec")
47 {
48 
49 }
50 //____________________________________________________________________________
52 XSecIntegratorI("genie::COHXSec", config)
53 {
54 
55 }
56 //____________________________________________________________________________
58 {
59 
60 }
61 //____________________________________________________________________________
63  const XSecAlgorithmI * model, const Interaction * in) const
64 {
65  if(! model->ValidProcess(in) ) return 0.;
66 
67  const KPhaseSpace & kps = in->PhaseSpace();
68  if(!kps.IsAboveThreshold()) {
69  LOG("COHXSec", pDEBUG) << "*** Below energy threshold";
70  return 0;
71  }
72  Range1D_t xl = kps.Limits(kKVx);
73  Range1D_t yl = kps.Limits(kKVy);
74  Range1D_t Q2l;
76  Q2l.max = fQ2Max;
77 
78  Interaction * interaction = new Interaction(*in);
79  interaction->SetBit(kISkipProcessChk);
80  //interaction->SetBit(kISkipKinematicChk);
81 
82  double xsec = 0.0;
83 
84  if (model->Id().Name() == "genie::ReinSehgalCOHPiPXSec") {
85  LOG("COHXSec", pINFO)
86  << "x integration range = [" << xl.min << ", " << xl.max << "]";
87  LOG("COHXSec", pINFO)
88  << "y integration range = [" << yl.min << ", " << yl.max << "]";
89 
90  ROOT::Math::IBaseFunctionMultiDim * func =
91  new utils::gsl::d2XSec_dxdy_E(model, interaction);
94 
95  double abstol = 1; //We mostly care about relative tolerance.
96  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
97  if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) {
98  ROOT::Math::AdaptiveIntegratorMultiDim * cast =
99  dynamic_cast<ROOT::Math::AdaptiveIntegratorMultiDim*>( ig.GetIntegrator() );
100  assert(cast);
101  cast->SetMinPts(fGSLMinEval);
102  }
103 
104  double kine_min[2] = { xl.min, yl.min };
105  double kine_max[2] = { xl.max, yl.max };
106  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
107  delete func;
108  }
109  else if (model->Id().Name() == "genie::BergerSehgalCOHPiPXSec2015")
110  {
111  ROOT::Math::IBaseFunctionMultiDim * func =
112  new utils::gsl::d2XSec_dQ2dy_E(model, interaction);
115  ROOT::Math::IntegratorMultiDim ig(ig_type);
116  ig.SetRelTolerance(fGSLRelTol);
117  ig.SetFunction(*func);
118  if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) {
119  ROOT::Math::AdaptiveIntegratorMultiDim * cast =
120  dynamic_cast<ROOT::Math::AdaptiveIntegratorMultiDim*>( ig.GetIntegrator() );
121  assert(cast);
122  cast->SetMinPts(fGSLMinEval);
123  }
124  double kine_min[2] = { Q2l.min, yl.min };
125  double kine_max[2] = { Q2l.max, yl.max };
126  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
127  delete func;
128  }
129  else if (model->Id().Name() == "genie::BergerSehgalFMCOHPiPXSec2015")
130  {
131  Range1D_t tl;
133  tl.max = fTMax;
134 
135  ROOT::Math::IBaseFunctionMultiDim * func =
136  new utils::gsl::d2XSec_dQ2dydt_E(model, interaction);
139  ROOT::Math::IntegratorMultiDim ig(ig_type);
140  ig.SetRelTolerance(fGSLRelTol);
141  ig.SetFunction(*func);
142  if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) {
143  ROOT::Math::AdaptiveIntegratorMultiDim * cast =
144  dynamic_cast<ROOT::Math::AdaptiveIntegratorMultiDim*>( ig.GetIntegrator() );
145  assert(cast);
146  cast->SetMinPts(fGSLMinEval);
147  }
148  double kine_min[3] = { Q2l.min, yl.min, tl.min };
149  double kine_max[3] = { Q2l.max, yl.max, tl.max };
150  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
151  delete func;
152  }
153 
154  const InitialState & init_state = in->InitState();
155  double Ev = init_state.ProbeE(kRfLab);
156  LOG("COHXSec", pINFO) << "XSec[COH] (E = " << Ev << " GeV) = " << xsec;
157 
158  delete interaction;
159 
160  return xsec;
161 }
162 //____________________________________________________________________________
164 {
165  Algorithm::Configure(config);
166  this->LoadConfig();
167 }
168 //____________________________________________________________________________
170 {
171  Algorithm::Configure(config);
172  this->LoadConfig();
173 }
174 //____________________________________________________________________________
176 {
177 
178  // Get GSL integration type & relative tolerance
179  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
180  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-2 ) ;
181 
182  int max_eval, min_eval ;
183  GetParamDef( "gsl-max-eval", max_eval, 500000 ) ;
184  GetParamDef( "gsl-min-eval", min_eval, 5000 ) ;
185 
186  fGSLMaxEval = (unsigned int) max_eval ;
187  fGSLMinEval = (unsigned int) min_eval ;
188 
189  //-- COH model parameter t_max for t = (q - p_pi)^2
190  GetParam("COH-t-max", fTMax ) ;
191 
192  //-- COH model bounds of integration for Q^2
193  GetParam( "COH-Q2-min", fQ2Min ) ;
194  GetParam( "COH-Q2-max", fQ2Max ) ;
195 
196 }
197 //____________________________________________________________________________
Cross Section Calculation Interface.
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
void LoadConfig(void)
Definition: COHXSec.cxx:175
double fQ2Max
upper bound of integration for Q^2 in Berger-Sehgal Model
Definition: COHXSec.h:46
Cross Section Integrator Interface.
double fTMax
upper bound for t = (q - p_pi)^2
Definition: COHXSec.h:47
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:67
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Definition: config.py:1
void Configure(const Registry &config)
Definition: COHXSec.cxx:163
static const double cm2
Definition: Units.h:77
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:97
string Name(void) const
Definition: AlgId.h:45
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
Double_t xsec[nknots]
Definition: testXsec.C:47
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
Definition: COHXSec.cxx:62
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
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
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?
assert(nhit_max >=nhit_nbins)
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:53
virtual ~COHXSec()
Definition: COHXSec.cxx:57
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:45
const XML_Char XML_Content * model
Definition: expat.h:151
double ProbeE(RefFrame_t rf) const
int fGSLMinEval
GSL min evaluations. Ignored by some integrators.
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)