DMDISXSec.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: Joshua Berger <jberger \at physics.wisc.edu>
8  University of Wisconsin-Madison
9 
10  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
11  University of Liverpool & STFC Rutherford Appleton Lab
12 */
13 //____________________________________________________________________________
14 
15 #include <TMath.h>
16 #include <Math/IFunction.h>
17 #include <Math/IntegratorMultiDim.h>
18 #include "Math/AdaptiveIntegratorMultiDim.h"
19 
30 #include "Framework/Utils/RunOpt.h"
32 #include "Framework/Utils/Range1.h"
33 #include "Framework/Utils/Cache.h"
37 
38 using namespace genie;
39 using namespace genie::controls;
40 using namespace genie::constants;
41 
42 //____________________________________________________________________________
44 XSecIntegratorI("genie::DMDISXSec")
45 {
46 
47 }
48 //____________________________________________________________________________
50 XSecIntegratorI("genie::DMDISXSec", config)
51 {
52 
53 }
54 //____________________________________________________________________________
56 {
57 
58 }
59 //____________________________________________________________________________
61  const XSecAlgorithmI * model, const Interaction * in) const
62 {
63  if(! model->ValidProcess(in) ) return 0.;
64 
65  const KPhaseSpace & kps = in->PhaseSpace();
66  if(!kps.IsAboveThreshold()) {
67  LOG("DMDISXSec", pDEBUG) << "*** Below energy threshold";
68  return 0;
69  }
70 
71  const InitialState & init_state = in->InitState();
72  double Ed = init_state.ProbeE(kRfHitNucRest);
73 
74  int nucpdgc = init_state.Tgt().HitNucPdg();
75  int NNucl = (pdg::IsProton(nucpdgc)) ?
76  init_state.Tgt().Z() : init_state.Tgt().N();
77 
78  // If the input interaction is off a nuclear target, then chek whether
79  // the corresponding free nucleon cross section already exists at the
80  // cross section spline list.
81  // If yes, calculate the nuclear cross section based on that value.
82  //
84  if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() ) {
85  Interaction * interaction = new Interaction(*in);
86  Target * target = interaction->InitStatePtr()->TgtPtr();
87  if(pdg::IsProton(nucpdgc)) { target->SetId(kPdgTgtFreeP); }
88  else { target->SetId(kPdgTgtFreeN); }
89  if(xsl->SplineExists(model,interaction)) {
90  const Spline * spl = xsl->GetSpline(model, interaction);
91  double xsec = spl->Evaluate(Ed);
92  LOG("DMDISXSec", pINFO)
93  << "From XSecSplineList: XSec[DIS,free nucleon] (E = " << Ed << " GeV) = " << xsec;
94  if(! interaction->TestBit(kIAssumeFreeNucleon) ) {
95  xsec *= NNucl;
96  LOG("DMDISXSec", pINFO) << "XSec[DIS] (E = " << Ed << " GeV) = " << xsec;
97  }
98  delete interaction;
99  return xsec;
100  }
101  delete interaction;
102  }
103 
104  // There was no corresponding free nucleon spline saved in XSecSplineList that
105  // could be used to speed up this calculation.
106  // Check whether local caching of free nucleon cross sections is allowed.
107  // If yes, store free nucleon cross sections at a cache branch and use those
108  // at any subsequent call.
109  //
110  bool precalc_bare_xsec = RunOpt::Instance()->BareXSecPreCalc();
111  if(precalc_bare_xsec) {
112  Cache * cache = Cache::Instance();
113  Interaction * interaction = new Interaction(*in);
114  string key = this->CacheBranchName(model,interaction);
115  LOG("DMDISXSec", pINFO) << "Finding cache branch with key: " << key;
116  CacheBranchFx * cache_branch =
117  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
118  if(!cache_branch) {
119  this->CacheFreeNucleonXSec(model,interaction);
120  cache_branch =
121  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
122  assert(cache_branch);
123  }
124  const CacheBranchFx & cb = (*cache_branch);
125  double xsec = cb(Ed);
126  if(! interaction->TestBit(kIAssumeFreeNucleon) ) { xsec *= NNucl; }
127  LOG("DMDISXSec", pINFO) << "XSec[DIS] (E = " << Ed << " GeV) = " << xsec;
128  delete interaction;
129  return xsec;
130  }
131  else {
132  // Just go ahead and integrate the input differential cross section for the
133  // specified interaction.
134  //
135  Interaction * interaction = new Interaction(*in);
136  interaction->SetBit(kISkipProcessChk);
137 // interaction->SetBit(kISkipKinematicChk);
138 
139  // **Important note**
140  // Based on discussions with Hugh at the GENIE mini-workshop / RAL - July '07
141  // The DIS nuclear corrections re-distribute the strength in x,y but do not
142  // affect the total cross-section They should be disabled at this step.
143  // But they should be enabled at the DIS thread's kinematical selection.
144  // Since nuclear corrections don't need to be included at this stage, all the
145  // nuclear cross sections can be trivially built from the free nucleon ones.
146  //
147  interaction->SetBit(kINoNuclearCorrection);
148 
149  Range1D_t Wl = kps.WLim();
150  Range1D_t Q2l = kps.Q2Lim();
151  LOG("DMDISXSec", pINFO)
152  << "W integration range = [" << Wl.min << ", " << Wl.max << "]";
153  LOG("DMDISXSec", pINFO)
154  << "Q2 integration range = [" << Q2l.min << ", " << Q2l.max << "]";
155 
156  bool phsp_ok =
157  (Q2l.min >= 0. && Q2l.max >= 0. && Q2l.max >= Q2l.min &&
158  Wl.min >= 0. && Wl.max >= 0. && Wl.max >= Wl.min);
159 
160  double xsec = 0.;
161 
162  if(phsp_ok) {
163  ROOT::Math::IBaseFunctionMultiDim * func =
164  new utils::gsl::d2XSec_dWdQ2_E(model, interaction);
167 
168  double abstol = 1; //We mostly care about relative tolerance.
169  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
170  double kine_min[2] = { Wl.min, Q2l.min };
171  double kine_max[2] = { Wl.max, Q2l.max };
172  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
173  delete func;
174  }//phase space ok?
175 
176  LOG("DMDISXSec", pINFO) << "XSec[DIS] (E = " << Ed << " GeV) = " << xsec;
177 
178  delete interaction;
179 
180  return xsec;
181  }
182  return 0;
183 }
184 //____________________________________________________________________________
186 {
187  Algorithm::Configure(config);
188  this->LoadConfig();
189 }
190 //____________________________________________________________________________
192 {
193  Algorithm::Configure(config);
194  this->LoadConfig();
195 }
196 //____________________________________________________________________________
198 {
199  // Get GSL integration type & relative tolerance
200  GetParamDef("gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
201  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-2 ) ;
202 
203  int max_eval, min_eval ;
204  GetParamDef( "gsl-max-eval", max_eval, 500000 ) ;
205  GetParamDef( "gsl-min-eval", min_eval, 10000 ) ;
206 
207  fGSLMaxEval = (unsigned int) max_eval ;
208  fGSLMinEval = (unsigned int) min_eval ;
209 
210  // Energy range for cached splines
211  GetParam( "GVLD-Emin", fVldEmin) ;
212  GetParam( "GVLD-Emax", fVldEmax) ;
213 
214 }
215 //____________________________________________________________________________
217  const XSecAlgorithmI * model, const Interaction * interaction) const
218 {
219  LOG("DMDISXSec", pWARN)
220  << "Wait while computing/caching free nucleon DIS xsections first...";
221 
222  // Create the cache branch
223  Cache * cache = Cache::Instance();
224  string key = this->CacheBranchName(model,interaction);
225  CacheBranchFx * cache_branch =
226  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
227  assert(!cache_branch);
228  cache_branch = new CacheBranchFx("DMDIS XSec");
229  cache->AddCacheBranch(key, cache_branch);
230 
231  // Tweak interaction to be on a free nucleon target
232  Target * target = interaction->InitStatePtr()->TgtPtr();
233  int nucpdgc = target->HitNucPdg();
234  if(pdg::IsProton(nucpdgc)) { target->SetId(kPdgTgtFreeP); }
235  else { target->SetId(kPdgTgtFreeN); }
236 
237  // Compute threshold
238  const KPhaseSpace & kps = interaction->PhaseSpace();
239  double Ethr = kps.Threshold();
240 
241  // Compute the number of spline knots - use at least 10 knots per decade
242  // && at least 40 knots in the full energy range
243  const double Emin = fVldEmin/3.;
244  const double Emax = fVldEmax*3.;
245  const int nknots_min = (int) (10*(TMath::Log(Emax) - TMath::Log(Emin)));
246  const int nknots = TMath::Max(40, nknots_min);
247 
248  // Distribute the knots in the energy range as is being done in the
249  // XSecSplineList so that the energy threshold is treated correctly
250  // in the spline - see comments there in.
251  double * E = new double[nknots];
252  int nkb = (Ethr>Emin) ? 5 : 0; // number of knots < threshold
253  int nka = nknots-nkb; // number of knots >= threshold
254  // knots < energy threshold
255  double dEb = (Ethr>Emin) ? (Ethr - Emin) / nkb : 0;
256  for(int i=0; i<nkb; i++) {
257  E[i] = Emin + i*dEb;
258  }
259  // knots >= energy threshold
260  double E0 = TMath::Max(Ethr,Emin);
261  double dEa = (TMath::Log10(Emax) - TMath::Log10(E0)) /(nka-1);
262  for(int i=0; i<nka; i++) {
263  E[i+nkb] = TMath::Power(10., TMath::Log10(E0) + i * dEa);
264  }
265 
266  // Create the integrand
267  ROOT::Math::IBaseFunctionMultiDim * func =
268  new utils::gsl::d2XSec_dWdQ2_E(model, interaction);
269 
270  // Compute the cross section at the given set of knots
271  double Md = interaction->InitStatePtr()->GetProbeP4()->M();
272  double Md2 = Md*Md;
273  for(int ie=0; ie<nknots; ie++) {
274  LOG("DMDISXSec", pDEBUG) << "Dealing with knot " << ie << " out of " << nknots;
275  double Ed = E[ie];
276  double pd = TMath::Max(Ed*Ed - Md2,0.);
277  pd = TMath::Sqrt(pd);
278  TLorentzVector p4(0,0,pd,Ed);
279  interaction->InitStatePtr()->SetProbeP4(p4);
280  double xsec = 0.;
281  if(Ed>Ethr+kASmallNum) {
282  Range1D_t Wl = kps.WLim();
283  Range1D_t Q2l = kps.Q2Lim();
284  LOG("DMDISXSec", pINFO)
285  << "W integration range = [" << Wl.min << ", " << Wl.max << "]";
286  LOG("DMDISXSec", pINFO)
287  << "Q2 integration range = [" << Q2l.min << ", " << Q2l.max << "]";
288 
289  bool phsp_ok =
290  (Q2l.min >= 0. && Q2l.max >= 0. && Q2l.max >= Q2l.min &&
291  Wl.min >= 0. && Wl.max >= 0. && Wl.max >= Wl.min);
292 
293  if(phsp_ok) {
296  double abstol = 1; //We mostly care about relative tolerance.
297  ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval);
298 
299  if (ig_type == ROOT::Math::IntegrationMultiDim::kADAPTIVE) {
300  ROOT::Math::AdaptiveIntegratorMultiDim * cast =
301  dynamic_cast<ROOT::Math::AdaptiveIntegratorMultiDim*>( ig.GetIntegrator() );
302  assert(cast);
303  cast->SetMinPts(fGSLMinEval);
304  }
305  double kine_min[2] = { Wl.min, Q2l.min };
306  double kine_max[2] = { Wl.max, Q2l.max };
307  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
308  }// phase space limits ok?
309  }//Ev>threshold
310 
311  LOG("DMDISXSec", pNOTICE)
312  << "Caching: XSec[DMDIS] (E = " << Ed << " GeV) = "
313  << xsec / (1E-38 * units::cm2) << " x 1E-38 cm^2";
314  cache_branch->AddValues(Ed,xsec);
315  }//ie
316 
317  // Create the spline
318  cache_branch->CreateSpline();
319 
320  delete [] E;
321  delete func;
322 }
323 //____________________________________________________________________________
325  const XSecAlgorithmI * model, const Interaction * interaction) const
326 {
327 // Build a unique name for the cache branch
328 
329  Cache * cache = Cache::Instance();
330 
331  string algkey = model->Id().Key();
332  string ikey = interaction->AsString();
333  string key = cache->CacheBranchKey(algkey, ikey);
334  return key;
335 }
336 //____________________________________________________________________________
337 
Cross Section Calculation Interface.
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
XSecIntegratorI interface implementation.
Definition: DMDISXSec.cxx:60
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
void Configure(const Registry &config)
Definition: DMDISXSec.cxx:185
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
void CacheFreeNucleonXSec(const XSecAlgorithmI *model, const Interaction *in) const
Definition: DMDISXSec.cxx:216
Definition: config.py:1
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 AddCacheBranch(string key, CacheBranchI *branch)
Definition: Cache.cxx:97
Summary information for an interaction.
Definition: Interaction.h:56
virtual ~DMDISXSec()
Definition: DMDISXSec.cxx:55
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
void LoadConfig(void)
Definition: DMDISXSec.cxx:197
string CacheBranchKey(string k0, string k1="", string k2="") const
Definition: Cache.cxx:102
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
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
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
double fVldEmin
Definition: DMDISXSec.h:51
Target * TgtPtr(void) const
Definition: InitialState.h:68
string CacheBranchName(const XSecAlgorithmI *model, const Interaction *in) const
Definition: DMDISXSec.cxx:324
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 fVldEmax
Definition: DMDISXSec.h:52
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
double Emax
const InitialState & InitState(void) const
Definition: Interaction.h:69
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
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
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)