DISXSec.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  @ Jan 19, 2008 - CA
14  Modify the way knots are distributed in the cached free nucleon DIS cross
15  section splines so that the energy threshold is treated more accurately
16  (see also XSecSplineList.cxx).
17  @ Sep 07, 2009 - CA
18  Integrated with GNU Numerical Library (GSL) via ROOT's MathMore library.
19  @ Oct 30, 2009 - CA
20  Fix problem reported by Hyupwoo Lee (Rochester) using GENIE in electron
21  scattering mode. Check kinematical limits before integration to avoid
22  problems when users override the physical limits raising the minimum Q2
23  (for computational efficiency in certain cases; depending on the detector
24  acceptance).
25  @ Jan 29, 2013 - CA
26  Don't look-up depreciated $GDISABLECACHING environmental variable.
27  Use the RunOpt singleton instead.
28 
29 */
30 //____________________________________________________________________________
31 
32 #include <TMath.h>
33 #include <Math/IFunction.h>
34 #include <Math/IntegratorMultiDim.h>
35 #include "Math/AdaptiveIntegratorMultiDim.h"
36 
47 #include "Framework/Utils/RunOpt.h"
49 #include "Framework/Utils/Range1.h"
50 #include "Framework/Utils/Cache.h"
54 
55 using namespace genie;
56 using namespace genie::controls;
57 using namespace genie::constants;
58 
59 //____________________________________________________________________________
61 XSecIntegratorI("genie::DISXSec")
62 {
63 
64 }
65 //____________________________________________________________________________
67 XSecIntegratorI("genie::DISXSec", config)
68 {
69 
70 }
71 //____________________________________________________________________________
73 {
74 
75 }
76 //____________________________________________________________________________
78  const XSecAlgorithmI * model, const Interaction * in) const
79 {
80  if(! model->ValidProcess(in) ) return 0.;
81 
82  const KPhaseSpace & kps = in->PhaseSpace();
83  if(!kps.IsAboveThreshold()) {
84  LOG("DISXSec", pDEBUG) << "*** Below energy threshold";
85  return 0;
86  }
87 
88  const InitialState & init_state = in->InitState();
89  double Ev = init_state.ProbeE(kRfHitNucRest);
90 
91  int nucpdgc = init_state.Tgt().HitNucPdg();
92  int NNucl = (pdg::IsProton(nucpdgc)) ?
93  init_state.Tgt().Z() : init_state.Tgt().N();
94 
95  // If the input interaction is off a nuclear target, then chek whether
96  // the corresponding free nucleon cross section already exists at the
97  // cross section spline list.
98  // If yes, calculate the nuclear cross section based on that value.
99  //
101  if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() ) {
102  Interaction * interaction = new Interaction(*in);
103  Target * target = interaction->InitStatePtr()->TgtPtr();
104  if(pdg::IsProton(nucpdgc)) { target->SetId(kPdgTgtFreeP); }
105  else { target->SetId(kPdgTgtFreeN); }
106  if(xsl->SplineExists(model,interaction)) {
107  const Spline * spl = xsl->GetSpline(model, interaction);
108  double xsec = spl->Evaluate(Ev);
109  LOG("DISXSec", pINFO)
110  << "From XSecSplineList: XSec[DIS,free nucleon] (E = " << Ev << " GeV) = " << xsec;
111  if(! interaction->TestBit(kIAssumeFreeNucleon) ) {
112  xsec *= NNucl;
113  LOG("DISXSec", pINFO) << "XSec[DIS] (E = " << Ev << " GeV) = " << xsec;
114  }
115  delete interaction;
116  return xsec;
117  }
118  delete interaction;
119  }
120 
121  // There was no corresponding free nucleon spline saved in XSecSplineList that
122  // could be used to speed up this calculation.
123  // Check whether local caching of free nucleon cross sections is allowed.
124  // If yes, store free nucleon cross sections at a cache branch and use those
125  // at any subsequent call.
126  //
127  bool precalc_bare_xsec = RunOpt::Instance()->BareXSecPreCalc();
128  if(precalc_bare_xsec) {
129  Cache * cache = Cache::Instance();
130  Interaction * interaction = new Interaction(*in);
131  string key = this->CacheBranchName(model,interaction);
132  LOG("DISXSec", pINFO) << "Finding cache branch with key: " << key;
133  CacheBranchFx * cache_branch =
134  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
135  if(!cache_branch) {
136  this->CacheFreeNucleonXSec(model,interaction);
137  cache_branch =
138  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
139  assert(cache_branch);
140  }
141  const CacheBranchFx & cb = (*cache_branch);
142  double xsec = cb(Ev);
143  if(! interaction->TestBit(kIAssumeFreeNucleon) ) { xsec *= NNucl; }
144  LOG("DISXSec", pINFO) << "XSec[DIS] (E = " << Ev << " GeV) = " << xsec;
145  delete interaction;
146  return xsec;
147  }
148  else {
149  // Just go ahead and integrate the input differential cross section for the
150  // specified interaction.
151  //
152  Interaction * interaction = new Interaction(*in);
153  interaction->SetBit(kISkipProcessChk);
154 // interaction->SetBit(kISkipKinematicChk);
155 
156  // **Important note**
157  // Based on discussions with Hugh at the GENIE mini-workshop / RAL - July '07
158  // The DIS nuclear corrections re-distribute the strength in x,y but do not
159  // affect the total cross-section They should be disabled at this step.
160  // But they should be enabled at the DIS thread's kinematical selection.
161  // Since nuclear corrections don't need to be included at this stage, all the
162  // nuclear cross sections can be trivially built from the free nucleon ones.
163  //
164  interaction->SetBit(kINoNuclearCorrection);
165 
166  Range1D_t Wl = kps.WLim();
167  Range1D_t Q2l = kps.Q2Lim();
168  LOG("DISXSec", pINFO)
169  << "W integration range = [" << Wl.min << ", " << Wl.max << "]";
170  LOG("DISXSec", pINFO)
171  << "Q2 integration range = [" << Q2l.min << ", " << Q2l.max << "]";
172 
173  bool phsp_ok =
174  (Q2l.min >= 0. && Q2l.max >= 0. && Q2l.max >= Q2l.min &&
175  Wl.min >= 0. && Wl.max >= 0. && Wl.max >= Wl.min);
176 
177  double xsec = 0.;
178 
179  if(phsp_ok) {
180  ROOT::Math::IBaseFunctionMultiDim * func =
181  new utils::gsl::d2XSec_dWdQ2_E(model, interaction);
184 
185  double abstol = 1; //We mostly care about relative tolerance.
186  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
187  double kine_min[2] = { Wl.min, Q2l.min };
188  double kine_max[2] = { Wl.max, Q2l.max };
189  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
190  delete func;
191  }//phase space ok?
192 
193  LOG("DISXSec", pINFO) << "XSec[DIS] (E = " << Ev << " GeV) = " << xsec;
194 
195  delete interaction;
196 
197  return xsec;
198  }
199  return 0;
200 }
201 //____________________________________________________________________________
203 {
204  Algorithm::Configure(config);
205  this->LoadConfig();
206 }
207 //____________________________________________________________________________
209 {
210  Algorithm::Configure(config);
211  this->LoadConfig();
212 }
213 //____________________________________________________________________________
215 {
216  // Get GSL integration type & relative tolerance
217  GetParamDef("gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
218  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-2 ) ;
219 
220  int max_eval, min_eval ;
221  GetParamDef( "gsl-max-eval", max_eval, 500000 ) ;
222  GetParamDef( "gsl-min-eval", min_eval, 10000 ) ;
223 
224  fGSLMaxEval = (unsigned int) max_eval ;
225  fGSLMinEval = (unsigned int) min_eval ;
226 
227  // Energy range for cached splines
228  GetParam( "GVLD-Emin", fVldEmin) ;
229  GetParam( "GVLD-Emax", fVldEmax) ;
230 
231 }
232 //____________________________________________________________________________
234  const XSecAlgorithmI * model, const Interaction * interaction) const
235 {
236  LOG("DISXSec", pWARN)
237  << "Wait while computing/caching free nucleon DIS xsections first...";
238 
239  // Create the cache branch
240  Cache * cache = Cache::Instance();
241  string key = this->CacheBranchName(model,interaction);
242  CacheBranchFx * cache_branch =
243  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
244  assert(!cache_branch);
245  cache_branch = new CacheBranchFx("DIS XSec");
246  cache->AddCacheBranch(key, cache_branch);
247 
248  // Tweak interaction to be on a free nucleon target
249  Target * target = interaction->InitStatePtr()->TgtPtr();
250  int nucpdgc = target->HitNucPdg();
251  if(pdg::IsProton(nucpdgc)) { target->SetId(kPdgTgtFreeP); }
252  else { target->SetId(kPdgTgtFreeN); }
253 
254  // Compute threshold
255  const KPhaseSpace & kps = interaction->PhaseSpace();
256  double Ethr = kps.Threshold();
257 
258  // Compute the number of spline knots - use at least 10 knots per decade
259  // && at least 40 knots in the full energy range
260  const double Emin = fVldEmin/3.;
261  const double Emax = fVldEmax*3.;
262  const int nknots_min = (int) (10*(TMath::Log(Emax) - TMath::Log(Emin)));
263  const int nknots = TMath::Max(40, nknots_min);
264 
265  // Distribute the knots in the energy range as is being done in the
266  // XSecSplineList so that the energy threshold is treated correctly
267  // in the spline - see comments there in.
268  double * E = new double[nknots];
269  int nkb = (Ethr>Emin) ? 5 : 0; // number of knots < threshold
270  int nka = nknots-nkb; // number of knots >= threshold
271  // knots < energy threshold
272  double dEb = (Ethr>Emin) ? (Ethr - Emin) / nkb : 0;
273  for(int i=0; i<nkb; i++) {
274  E[i] = Emin + i*dEb;
275  }
276  // knots >= energy threshold
277  double E0 = TMath::Max(Ethr,Emin);
278  double dEa = (TMath::Log10(Emax) - TMath::Log10(E0)) /(nka-1);
279  for(int i=0; i<nka; i++) {
280  E[i+nkb] = TMath::Power(10., TMath::Log10(E0) + i * dEa);
281  }
282 
283  // Create the integrand
284  ROOT::Math::IBaseFunctionMultiDim * func =
285  new utils::gsl::d2XSec_dWdQ2_E(model, interaction);
286 
287  // Compute the cross section at the given set of knots
288  for(int ie=0; ie<nknots; ie++) {
289  double Ev = E[ie];
290  TLorentzVector p4(0,0,Ev,Ev);
291  interaction->InitStatePtr()->SetProbeP4(p4);
292  double xsec = 0.;
293  if(Ev>Ethr+kASmallNum) {
294  Range1D_t Wl = kps.WLim();
295  Range1D_t Q2l = kps.Q2Lim();
296  LOG("DISXSec", pINFO)
297  << "W integration range = [" << Wl.min << ", " << Wl.max << "]";
298  LOG("DISXSec", pINFO)
299  << "Q2 integration range = [" << Q2l.min << ", " << Q2l.max << "]";
300 
301  bool phsp_ok =
302  (Q2l.min >= 0. && Q2l.max >= 0. && Q2l.max >= Q2l.min &&
303  Wl.min >= 0. && Wl.max >= 0. && Wl.max >= Wl.min);
304 
305  if(phsp_ok) {
308  double abstol = 1; //We mostly care about relative tolerance.
309  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
310 
311  if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) {
312  ROOT::Math::AdaptiveIntegratorMultiDim * cast =
313  dynamic_cast<ROOT::Math::AdaptiveIntegratorMultiDim*>( ig.GetIntegrator() );
314  assert(cast);
315  cast->SetMinPts(fGSLMinEval);
316  }
317 
318  double kine_min[2] = { Wl.min, Q2l.min };
319  double kine_max[2] = { Wl.max, Q2l.max };
320  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
321  }// phase space limits ok?
322  }//Ev>threshold
323 
324  LOG("DISXSec", pNOTICE)
325  << "Caching: XSec[DIS] (E = " << Ev << " GeV) = "
326  << xsec / (1E-38 * units::cm2) << " x 1E-38 cm^2";
327  cache_branch->AddValues(Ev,xsec);
328  }//ie
329 
330  // Create the spline
331  cache_branch->CreateSpline();
332 
333  delete [] E;
334  delete func;
335 }
336 //____________________________________________________________________________
338  const XSecAlgorithmI * model, const Interaction * interaction) const
339 {
340 // Build a unique name for the cache branch
341 
342  Cache * cache = Cache::Instance();
343 
344  string algkey = model->Id().Key();
345  string ikey = interaction->AsString();
346  string key = cache->CacheBranchKey(algkey, ikey);
347  return key;
348 }
349 //____________________________________________________________________________
350 
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
string fGSLIntgType
name of GSL numerical integrator
const XML_Char * target
Definition: expat.h:268
void SetProbeP4(const TLorentzVector &P4)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:67
int HitNucPdg(void) const
Definition: Target.cxx:321
A simple [min,max] interval for doubles.
Definition: Range1.h:43
double Threshold(void) const
Energy threshold.
Definition: KPhaseSpace.cxx:81
bool SplineExists(const XSecAlgorithmI *alg, const Interaction *i) const
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:47
bool IsNucleus(void) const
Definition: Target.cxx:289
Definition: config.py:1
void LoadConfig(void)
Definition: DISXSec.cxx:214
static XSecSplineList * Instance()
Range1D_t Q2Lim(void) const
Q2 limits.
double Evaluate(double x) const
Definition: Spline.cxx:362
static const double cm2
Definition: Units.h:77
void SetId(int pdgc)
Definition: Target.cxx:166
string AsString(void) const
void Configure(const Registry &config)
Definition: DISXSec.cxx:202
void AddCacheBranch(string key, CacheBranchI *branch)
Definition: Cache.cxx:97
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
Definition: DISXSec.cxx:77
Summary information for an interaction.
Definition: Interaction.h:56
void AddValues(double x, double y)
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
bool BareXSecPreCalc(void) const
Definition: RunOpt.h:52
const UInt_t kINoNuclearCorrection
if set, inhibit nuclear corrections
Definition: Interaction.h:51
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const Int_t nknots
Definition: testXsec.C:24
bool IsEmpty(void) const
Float_t E
Definition: plot.C:20
virtual ~DISXSec()
Definition: DISXSec.cxx:72
string CacheBranchKey(string k0, string k1="", string k2="") const
Definition: Cache.cxx:102
double fVldEmax
Definition: DISXSec.h:49
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
Kinematical phase space.
Definition: KPhaseSpace.h:34
const int kPdgTgtFreeN
Definition: PDGCodes.h:177
const int kPdgTgtFreeP
Definition: PDGCodes.h:176
void CacheFreeNucleonXSec(const XSecAlgorithmI *model, const Interaction *in) const
Definition: DISXSec.cxx:233
double func(double x, double y)
int Z(void) const
Definition: Target.h:69
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:63
Misc GENIE control constants.
Double_t xsec[nknots]
Definition: testXsec.C:47
#define pWARN
Definition: Messenger.h:61
CacheBranchI * FindCacheBranch(string key)
finding/adding cache branches
Definition: Cache.cxx:89
GENIE Cache Memory.
Definition: Cache.h:39
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
int N(void) const
Definition: Target.h:70
static RunOpt * Instance(void)
Definition: RunOpt.cxx:62
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
Target * TgtPtr(void) const
Definition: InitialState.h:68
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)
double fVldEmin
Definition: DISXSec.h:48
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
double Emax
const InitialState & InitState(void) const
Definition: Interaction.h:69
string CacheBranchName(const XSecAlgorithmI *model, const Interaction *in) const
Definition: DISXSec.cxx:337
double min
Definition: Range1.h:53
#define pNOTICE
Definition: Messenger.h:62
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:67
static Cache * Instance(void)
Definition: Cache.cxx:76
A simple cache branch storing the cached data in a TNtuple.
Definition: CacheBranchFx.h:38
const Spline * GetSpline(const XSecAlgorithmI *alg, const Interaction *i) const
List of cross section vs energy splines.
const XML_Char XML_Content * model
Definition: expat.h:151
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
int fGSLMinEval
GSL min evaluations. Ignored by some integrators.
string Key(void) const
Definition: AlgId.h:47
Range1D_t WLim(void) const
W limits.
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)