ReinSehgalRESXSec.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  @ Sep 07, 2009 - CA
14  Integrated with GNU Numerical Library (GSL) via ROOT's MathMore library.
15  @ Jan 29, 2013 - CA
16  Don't look-up depreciated $GDISABLECACHING environmental variable.
17  Use the RunOpt singleton instead.
18 */
19 //____________________________________________________________________________
20 
21 #include <TMath.h>
22 #include <Math/IFunction.h>
23 #include <Math/IntegratorMultiDim.h>
24 
35 #include "Framework/Utils/RunOpt.h"
38 #include "Framework/Utils/Cache.h"
44 
45 using namespace genie;
46 using namespace genie::constants;
47 //using namespace genie::units;
48 
49 //____________________________________________________________________________
51 ReinSehgalRESXSecWithCache("genie::ReinSehgalRESXSec")
52 {
53 
54 }
55 //____________________________________________________________________________
57 ReinSehgalRESXSecWithCache("genie::ReinSehgalRESXSec", config)
58 {
59 
60 }
61 //____________________________________________________________________________
63 {
64 
65 }
66 //____________________________________________________________________________
68  const XSecAlgorithmI * model, const Interaction * interaction) const
69 {
70  if(! model->ValidProcess(interaction) ) return 0.;
72 
73  const KPhaseSpace & kps = interaction->PhaseSpace();
74  if(!kps.IsAboveThreshold()) {
75  LOG("ReinSehgalRESXSec", pDEBUG) << "*** Below energy threshold";
76  return 0;
77  }
78 
79  //-- Get init state and process information
80  const InitialState & init_state = interaction->InitState();
81  const ProcessInfo & proc_info = interaction->ProcInfo();
82  const Target & target = init_state.Tgt();
83 
85  int nucleon_pdgc = target.HitNucPdg();
86  int nu_pdgc = init_state.ProbePdg();
87 
88  //-- Get neutrino energy in the struck nucleon rest frame
89  double Ev = init_state.ProbeE(kRfHitNucRest);
90 
91  //-- Get the requested resonance
92  Resonance_t res = interaction->ExclTag().Resonance();
93 
94  // If the input interaction is off a nuclear target, then chek whether
95  // the corresponding free nucleon cross section already exists at the
96  // cross section spline list.
97  // If yes, calculate the nuclear cross section based on that value.
98  //
100  if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() ) {
101  Interaction * in = new Interaction(*interaction);
102  if(pdg::IsProton(nucleon_pdgc)) {
104  } else {
106  }
107  if(xsl->SplineExists(model,in)) {
108  const Spline * spl = xsl->GetSpline(model, in);
109  double xsec = spl->Evaluate(Ev);
110  SLOG("ReinSehgalResT", pNOTICE)
111  << "XSec[RES/" << utils::res::AsString(res)<< "/free] (Ev = "
112  << Ev << " GeV) = " << xsec/(1E-38 *genie::units::cm2) << " x 1E-38 cm^2";
113  if(! interaction->TestBit(kIAssumeFreeNucleon) ) {
114  int NNucl = (pdg::IsProton(nucleon_pdgc)) ? target.Z() : target.N();
115  xsec *= NNucl;
116  }
117  delete in;
118  return xsec;
119  }
120  delete in;
121  }
122 
123  // There was no corresponding free nucleon spline saved in XSecSplineList that
124  // could be used to speed up this calculation.
125  // Check whether local caching of free nucleon cross sections is allowed.
126  // If yes, store free nucleon cross sections at a cache branch and use those
127  // at any subsequent call.
128  //
129  bool bare_xsec_pre_calc = RunOpt::Instance()->BareXSecPreCalc();
130  if(bare_xsec_pre_calc && !fUsePauliBlocking) {
131  Cache * cache = Cache::Instance();
132  string key = this->CacheBranchName(res, it, nu_pdgc, nucleon_pdgc);
133  LOG("ReinSehgalResT", pINFO)
134  << "Finding cache branch with key: " << key;
135  CacheBranchFx * cache_branch =
136  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
137  if(!cache_branch) {
138  LOG("ReinSehgalResT", pWARN)
139  << "No cached RES v-production data for input neutrino"
140  << " (pdgc: " << nu_pdgc << ")";
141  LOG("ReinSehgalResT", pWARN)
142  << "Wait while computing/caching RES production xsec first...";
143 
144  this->CacheResExcitationXSec(interaction);
145 
146  LOG("ReinSehgalResT", pINFO) << "Done caching resonance xsec data";
147  LOG("ReinSehgalResT", pINFO)
148  << "Finding newly created cache branch with key: " << key;
149  cache_branch =
150  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
151  assert(cache_branch);
152  }
153  const CacheBranchFx & cbranch = (*cache_branch);
154 
155  //-- Get cached resonance neutrinoproduction xsec
156  // (If E>Emax, assume xsec = xsec(Emax) - but do not evaluate the
157  // cross section spline at the end of its energy range-)
158  double rxsec = (Ev<fEMax-1) ? cbranch(Ev) : cbranch(fEMax-1);
159 
160  SLOG("ReinSehgalResT", pNOTICE)
161  << "XSec[RES/" << utils::res::AsString(res)<< "/free] (Ev = "
162  << Ev << " GeV) = " << rxsec/(1E-38 *genie::units::cm2) << " x 1E-38 cm^2";
163 
164  if( interaction->TestBit(kIAssumeFreeNucleon) ) return rxsec;
165 
166  int NNucl = (pdg::IsProton(nucleon_pdgc)) ? target.Z() : target.N();
167  rxsec*=NNucl; // nuclear xsec
168  return rxsec;
169  } // disable local caching
170 
171  // Just go ahead and integrate the input differential cross section for the
172  // specified interaction.
173  else {
174 
175  Range1D_t rW = kps.Limits(kKVW);
176  Range1D_t rQ2 = kps.Limits(kKVQ2);
177 
178  LOG("ReinSehgalResC", pINFO)
179  << "*** Integrating d^2 XSec/dWdQ^2 for R: "
180  << utils::res::AsString(res) << " at Ev = " << Ev;
181  LOG("ReinSehgalResC", pINFO)
182  << "{W} = " << rW.min << ", " << rW.max;
183  LOG("ReinSehgalResC", pINFO)
184  << "{Q^2} = " << rQ2.min << ", " << rQ2.max;
185 
186  ROOT::Math::IBaseFunctionMultiDim * func =
187  new utils::gsl::d2XSec_dWdQ2_E(model, interaction);
190  ROOT::Math::IntegratorMultiDim ig(ig_type,0,fGSLRelTol,fGSLMaxEval);
191  ig.SetFunction(*func);
192  double kine_min[2] = { rW.min, rQ2.min };
193  double kine_max[2] = { rW.max, rQ2.max };
194  double xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
195 
196  delete func;
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, 0.01 ) ;
219  GetParamDef( "gsl-max-eval", fGSLMaxEval, 100000 ) ;
220  GetParam("UsePauliBlockingForRES", fUsePauliBlocking);
221  // Get upper E limit on res xsec spline (=f(E)) before assuming xsec=const
222  GetParamDef( "ESplineMax", fEMax, 100. ) ;
223  fEMax = TMath::Max(fEMax, 20.); // don't accept user Emax if less than 20 GeV
224 
225  // Create the baryon resonance list specified in the config.
226  fResList.Clear();
227  string resonances ;
228  GetParam( "ResonanceNameList", resonances ) ;
229  fResList.DecodeFromNameList(resonances);
230 
231 }
232 //____________________________________________________________________________
Cross Section Calculation Interface.
bool fUsePauliBlocking
account for Pauli blocking?
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
InteractionType_t InteractionTypeId(void) const
string fGSLIntgType
name of GSL numerical integrator
const XML_Char * target
Definition: expat.h:268
set< int >::iterator it
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
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
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 DecodeFromNameList(string list, string delimiter=",")
Definition: config.py:1
static XSecSplineList * Instance()
double Evaluate(double x) const
Definition: Spline.cxx:362
static const double cm2
Definition: Units.h:77
enum genie::EResonance Resonance_t
void SetId(int pdgc)
Definition: Target.cxx:166
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
Summary information for an interaction.
Definition: Interaction.h:56
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
bool BareXSecPreCalc(void) const
Definition: RunOpt.h:52
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
bool IsEmpty(void) const
Float_t E
Definition: plot.C:20
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
string CacheBranchName(Resonance_t r, InteractionType_t it, int nu, int nuc) const
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
int ProbePdg(void) const
Definition: InitialState.h:65
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
#define pINFO
Definition: Messenger.h:63
Resonance_t Resonance(void) const
Definition: XclsTag.h:62
void CacheResExcitationXSec(const Interaction *interaction) const
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
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
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
void Configure(const Registry &config)
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)
An ABC that caches resonance neutrinoproduction cross sections on free nucleons according to the Rein...
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
const InitialState & InitState(void) const
Definition: Interaction.h:69
const char * AsString(Resonance_t res)
resonance id -> string
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
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
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
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...
enum genie::EInteractionType InteractionType_t
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
double fGSLRelTol
required relative tolerance (error)