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