ReinSehgalRESXSecFast.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2017, GENIE Neutrino MC Generator Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Igor Kakorin <kakorin@jinr.ru>
8  Joint Institute for Nuclear Research - March 01, 2017
9  based on code of Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
10  University of Liverpool & STFC Rutherford Appleton Lab
11 
12  For the class documentation see the corresponding header file.
13 
14 
15 */
16 //____________________________________________________________________________
17 
18 #include <TMath.h>
19 #include <Math/IFunction.h>
20 #include <Math/IntegratorMultiDim.h>
21 
33 #include "Framework/Utils/RunOpt.h"
36 #include "Framework/Utils/Cache.h"
40 
41 using namespace genie;
42 using namespace genie::constants;
43 using namespace genie::units;
44 
45 //____________________________________________________________________________
47 ReinSehgalRESXSecWithCacheFast("genie::ReinSehgalRESXSecFast")
48 {
49 
50 }
51 //____________________________________________________________________________
53 ReinSehgalRESXSecWithCacheFast("genie::ReinSehgalRESXSecFast", config)
54 {
55 
56 }
57 //____________________________________________________________________________
59 {
60 
61 }
62 //____________________________________________________________________________
64  const XSecAlgorithmI * model, const Interaction * interaction) const
65 {
66  if(! model->ValidProcess(interaction) ) return 0.;
68 
69  const KPhaseSpace & kps = interaction->PhaseSpace();
70  if(!kps.IsAboveThreshold()) {
71  LOG("ReinSehgalRESXSecFast", pDEBUG) << "*** Below energy threshold";
72  return 0;
73  }
74 
75  //-- Get init state and process information
76  const InitialState & init_state = interaction->InitState();
77  const ProcessInfo & proc_info = interaction->ProcInfo();
78  const Target & target = init_state.Tgt();
79 
81  int nucleon_pdgc = target.HitNucPdg();
82  int nu_pdgc = init_state.ProbePdg();
83 
84  //-- Get neutrino energy in the struck nucleon rest frame
85  double Ev = init_state.ProbeE(kRfHitNucRest);
86 
87  //-- Get the requested resonance
88  Resonance_t res = interaction->ExclTag().Resonance();
89 
90  // If the input interaction is off a nuclear target, then chek whether
91  // the corresponding free nucleon cross section already exists at the
92  // cross section spline list.
93  // If yes, calculate the nuclear cross section based on that value.
94  //
96  if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() ) {
97  Interaction * in = new Interaction(*interaction);
98  if(pdg::IsProton(nucleon_pdgc)) {
100  } else {
102  }
103  if(xsl->SplineExists(model,in)) {
104  const Spline * spl = xsl->GetSpline(model, in);
105  double xsec = spl->Evaluate(Ev);
106  SLOG("ReinSehgalResTF", pNOTICE)
107  << "XSec[RES/" << utils::res::AsString(res)<< "/free] (Ev = "
108  << Ev << " GeV) = " << xsec/(1E-38 *cm2)<< " x 1E-38 cm^2";
109  if(! interaction->TestBit(kIAssumeFreeNucleon) ) {
110  int NNucl = (pdg::IsProton(nucleon_pdgc)) ? target.Z() : target.N();
111  xsec *= NNucl;
112  }
113  delete in;
114  return xsec;
115  }
116  delete in;
117  }
118 
119  // There was no corresponding free nucleon spline saved in XSecSplineList that
120  // could be used to speed up this calculation.
121  // Check whether local caching of free nucleon cross sections is allowed.
122  // If yes, store free nucleon cross sections at a cache branch and use those
123  // at any subsequent call.
124  //
125  bool bare_xsec_pre_calc = RunOpt::Instance()->BareXSecPreCalc();
126  if(bare_xsec_pre_calc && !fUsePauliBlocking) {
127  Cache * cache = Cache::Instance();
128  string key = this->CacheBranchName(res, it, nu_pdgc, nucleon_pdgc);
129  LOG("ReinSehgalResTF", pINFO)
130  << "Finding cache branch with key: " << key;
131  CacheBranchFx * cache_branch =
132  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
133  if(!cache_branch) {
134  LOG("ReinSehgalResTF", pWARN)
135  << "No cached RES v-production data for input neutrino"
136  << " (pdgc: " << nu_pdgc << ")";
137  LOG("ReinSehgalResTF", pWARN)
138  << "Wait while computing/caching RES production xsec first...";
139 
140  this->CacheResExcitationXSec(interaction);
141 
142  LOG("ReinSehgalResTF", pINFO) << "Done caching resonance xsec data";
143  LOG("ReinSehgalResTF", pINFO)
144  << "Finding newly created cache branch with key: " << key;
145  cache_branch =
146  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
147  assert(cache_branch);
148  }
149  const CacheBranchFx & cbranch = (*cache_branch);
150 
151  //-- Get cached resonance neutrinoproduction xsec
152  // (If E>Emax, assume xsec = xsec(Emax) - but do not evaluate the
153  // cross section spline at the end of its energy range-)
154  double rxsec = (Ev<fEMax-1) ? cbranch(Ev) : cbranch(fEMax-1);
155 
156  SLOG("ReinSehgalResTF", pNOTICE)
157  << "XSec[RES/" << utils::res::AsString(res)<< "/free] (Ev = "
158  << Ev << " GeV) = " << rxsec/(1E-38 *cm2)<< " x 1E-38 cm^2";
159 
160  if( interaction->TestBit(kIAssumeFreeNucleon) ) return rxsec;
161 
162  int NNucl = (pdg::IsProton(nucleon_pdgc)) ? target.Z() : target.N();
163  rxsec*=NNucl; // nuclear xsec
164  return rxsec;
165  } // disable local caching
166 
167  // Just go ahead and integrate the input differential cross section for the
168  // specified interaction.
169  else {
170 
171  Range1D_t rW = Range1D_t(0.0,1.0);
172  Range1D_t rQ2 = Range1D_t(0.0,1.0);
173 
174  LOG("ReinSehgalResTF", pINFO)
175  << "*** Integrating d^2 XSec/dWdQ^2 for R: "
176  << utils::res::AsString(res) << " at Ev = " << Ev;
177  LOG("ReinSehgalResTF", pINFO)
178  << "{W} = " << rW.min << ", " << rW.max;
179  LOG("ReinSehgalResTF", pINFO)
180  << "{Q^2} = " << rQ2.min << ", " << rQ2.max;
181 
182  ROOT::Math::IBaseFunctionMultiDim * func= new utils::gsl::d2XSecRESFast_dWQ2_E(model, interaction);
184  ROOT::Math::IntegratorMultiDim ig(ig_type,0,fGSLRelTol,fGSLMaxEval);
185  ig.SetFunction(*func);
186  double kine_min[2] = { rW.min, rQ2.min };
187  double kine_max[2] = { rW.max, rQ2.max };
188  double xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
189 
190  delete func;
191  return xsec;
192  }
193  return 0;
194 }
195 //____________________________________________________________________________
197 {
198  Algorithm::Configure(config);
199  this->LoadConfig();
200 }
201 //____________________________________________________________________________
203 {
204  Algorithm::Configure(config);
205  this->LoadConfig();
206 }
207 //____________________________________________________________________________
209 {
210 
211  // Get GSL integration type & relative tolerance
212  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
213  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 0.01 ) ;
214  GetParamDef( "gsl-max-eval", fGSLMaxEval, 100000 ) ;
215  GetParam("UsePauliBlockingForRES", fUsePauliBlocking);
216  // Get upper E limit on res xsec spline (=f(E)) before assuming xsec=const
217  GetParamDef( "ESplineMax", fEMax, 100. ) ;
218  fEMax = TMath::Max(fEMax, 20.); // don't accept user Emax if less than 20 GeV
219 
220  // Create the baryon resonance list specified in the config.
221  fResList.Clear();
222  string resonances ;
223  GetParam( "ResonanceNameList", resonances ) ;
224  fResList.DecodeFromNameList(resonances);
225 }
226 //____________________________________________________________________________
227 
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
Class that caches resonance neutrinoproduction cross sections on free nucleons according to the Rein-...
InteractionType_t InteractionTypeId(void) const
Physical System of Units.
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
static constexpr Double_t cm2
Definition: Munits.h:141
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=",")
void CacheResExcitationXSec(const Interaction *interaction) const
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
string CacheBranchName(Resonance_t r, InteractionType_t it, int nu, int nuc) const
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
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
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
bool fUsePauliBlocking
account for Pauli blocking?
Resonance_t Resonance(void) const
Definition: XclsTag.h:62
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
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
void Configure(const Registry &config)
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
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)
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)