GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ReinSehgalRESXSecFast.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  Igor Kakorin <kakorin@jinr.ru>
7  Joint Institute for Nuclear Research
8 
9  based on code of
10  Costas Andreopoulos <c.andreopoulos \at cern.ch>
11  University of Liverpool
12  */
13 //____________________________________________________________________________
14 
15 #include <TMath.h>
16 #include <Math/IFunction.h>
17 #include <Math/IntegratorMultiDim.h>
18 
20 #include "Framework/Conventions/GBuild.h"
30 #include "Framework/Utils/RunOpt.h"
33 #include "Framework/Utils/Cache.h"
37 
38 using namespace genie;
39 using namespace genie::constants;
40 using namespace genie::units;
41 
42 //____________________________________________________________________________
44 ReinSehgalRESXSecWithCacheFast("genie::ReinSehgalRESXSecFast")
45 {
46 
47 }
48 //____________________________________________________________________________
50 ReinSehgalRESXSecWithCacheFast("genie::ReinSehgalRESXSecFast", config)
51 {
52 
53 }
54 //____________________________________________________________________________
56 {
57 
58 }
59 //____________________________________________________________________________
61  const XSecAlgorithmI * model, const Interaction * interaction) const
62 {
63  if(! model->ValidProcess(interaction) ) return 0.;
64  fSingleResXSecModel = model;
65 
66  const KPhaseSpace & kps = interaction->PhaseSpace();
67  if(!kps.IsAboveThreshold()) {
68  LOG("ReinSehgalRESXSecFast", pDEBUG) << "*** Below energy threshold";
69  return 0;
70  }
71 
72  //-- Get init state and process information
73  const InitialState & init_state = interaction->InitState();
74  const ProcessInfo & proc_info = interaction->ProcInfo();
75  const Target & target = init_state.Tgt();
76 
77  InteractionType_t it = proc_info.InteractionTypeId();
78  int nucleon_pdgc = target.HitNucPdg();
79  int nu_pdgc = init_state.ProbePdg();
80 
81  //-- Get neutrino energy in the struck nucleon rest frame
82  double Ev = init_state.ProbeE(kRfHitNucRest);
83 
84  //-- Get the requested resonance
85  Resonance_t res = interaction->ExclTag().Resonance();
86 
87  // If the input interaction is off a nuclear target, then chek whether
88  // the corresponding free nucleon cross section already exists at the
89  // cross section spline list.
90  // If yes, calculate the nuclear cross section based on that value.
91  //
93  if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() ) {
94  Interaction * in = new Interaction(*interaction);
95  if(pdg::IsProton(nucleon_pdgc)) {
97  } else {
99  }
100  if(xsl->SplineExists(model,in)) {
101  const Spline * spl = xsl->GetSpline(model, in);
102  double xsec = spl->Evaluate(Ev);
103  SLOG("ReinSehgalResTF", pNOTICE)
104  << "XSec[RES/" << utils::res::AsString(res)<< "/free] (Ev = "
105  << Ev << " GeV) = " << xsec/(1E-38 *cm2)<< " x 1E-38 cm^2";
106  if(! interaction->TestBit(kIAssumeFreeNucleon) ) {
107  int NNucl = (pdg::IsProton(nucleon_pdgc)) ? target.Z() : target.N();
108  xsec *= NNucl;
109  }
110  delete in;
111  return xsec;
112  }
113  delete in;
114  }
115 
116  // There was no corresponding free nucleon spline saved in XSecSplineList that
117  // could be used to speed up this calculation.
118  // Check whether local caching of free nucleon cross sections is allowed.
119  // If yes, store free nucleon cross sections at a cache branch and use those
120  // at any subsequent call.
121  //
122  bool bare_xsec_pre_calc = RunOpt::Instance()->BareXSecPreCalc();
123  if(bare_xsec_pre_calc && !fUsePauliBlocking) {
124  Cache * cache = Cache::Instance();
125  string key = this->CacheBranchName(res, it, nu_pdgc, nucleon_pdgc);
126  LOG("ReinSehgalResTF", pINFO)
127  << "Finding cache branch with key: " << key;
128  CacheBranchFx * cache_branch =
129  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
130  if(!cache_branch) {
131  LOG("ReinSehgalResTF", pWARN)
132  << "No cached RES v-production data for input neutrino"
133  << " (pdgc: " << nu_pdgc << ")";
134  LOG("ReinSehgalResTF", pWARN)
135  << "Wait while computing/caching RES production xsec first...";
136 
137  this->CacheResExcitationXSec(interaction);
138 
139  LOG("ReinSehgalResTF", pINFO) << "Done caching resonance xsec data";
140  LOG("ReinSehgalResTF", pINFO)
141  << "Finding newly created cache branch with key: " << key;
142  cache_branch =
143  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
144  assert(cache_branch);
145  }
146  const CacheBranchFx & cbranch = (*cache_branch);
147 
148  //-- Get cached resonance neutrinoproduction xsec
149  // (If E>Emax, assume xsec = xsec(Emax) - but do not evaluate the
150  // cross section spline at the end of its energy range-)
151  double rxsec = (Ev<fEMax-1) ? cbranch(Ev) : cbranch(fEMax-1);
152 
153  SLOG("ReinSehgalResTF", pNOTICE)
154  << "XSec[RES/" << utils::res::AsString(res)<< "/free] (Ev = "
155  << Ev << " GeV) = " << rxsec/(1E-38 *cm2)<< " x 1E-38 cm^2";
156 
157  if( interaction->TestBit(kIAssumeFreeNucleon) ) return rxsec;
158 
159  int NNucl = (pdg::IsProton(nucleon_pdgc)) ? target.Z() : target.N();
160  rxsec*=NNucl; // nuclear xsec
161  return rxsec;
162  } // disable local caching
163 
164  // Just go ahead and integrate the input differential cross section for the
165  // specified interaction.
166  else {
167 
168  Range1D_t rW = Range1D_t(0.0,1.0);
169  Range1D_t rQ2 = Range1D_t(0.0,1.0);
170 
171  LOG("ReinSehgalResTF", pINFO)
172  << "*** Integrating d^2 XSec/dWdQ^2 for R: "
173  << utils::res::AsString(res) << " at Ev = " << Ev;
174  LOG("ReinSehgalResTF", pINFO)
175  << "{W} = " << rW.min << ", " << rW.max;
176  LOG("ReinSehgalResTF", pINFO)
177  << "{Q^2} = " << rQ2.min << ", " << rQ2.max;
178 
179  ROOT::Math::IBaseFunctionMultiDim * func= new utils::gsl::d2XSecRESFast_dWQ2_E(model, interaction);
180  ROOT::Math::IntegrationMultiDim::Type ig_type = utils::gsl::IntegrationNDimTypeFromString(fGSLIntgType);
181  ROOT::Math::IntegratorMultiDim ig(ig_type,0,fGSLRelTol,fGSLMaxEval);
182  ig.SetFunction(*func);
183  double kine_min[2] = { rW.min, rQ2.min };
184  double kine_max[2] = { rW.max, rQ2.max };
185  double xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
186 
187  delete func;
188  return xsec;
189  }
190  return 0;
191 }
192 //____________________________________________________________________________
194 {
195  Algorithm::Configure(config);
196  this->LoadConfig();
197 }
198 //____________________________________________________________________________
200 {
201  Algorithm::Configure(config);
202  this->LoadConfig();
203 }
204 //____________________________________________________________________________
206 {
207 
208  // Get GSL integration type & relative tolerance
209  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
210  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 0.01 ) ;
211  GetParamDef( "gsl-max-eval", fGSLMaxEval, 100000 ) ;
212  GetParam("UsePauliBlockingForRES", fUsePauliBlocking);
213  // Get upper E limit on res xsec spline (=f(E)) before assuming xsec=const
214  GetParamDef( "ESplineMax", fEMax, 100. ) ;
215  fEMax = TMath::Max(fEMax, 20.); // don't accept user Emax if less than 20 GeV
216 
217  // Create the baryon resonance list specified in the config.
218  fResList.Clear();
219  string resonances ;
220  GetParam( "ResonanceNameList", resonances ) ;
221  fResList.DecodeFromNameList(resonances);
222 }
223 //____________________________________________________________________________
Cross Section Calculation Interface.
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Class that caches resonance neutrinoproduction cross sections on free nucleons according to the Rein-...
InteractionType_t InteractionTypeId(void) const
string fGSLIntgType
name of GSL numerical integrator
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
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 DecodeFromNameList(string list, string delimiter=",")
void CacheResExcitationXSec(const Interaction *interaction) const
static XSecSplineList * Instance()
double Evaluate(double x) const
Definition: Spline.cxx:363
enum genie::EResonance Resonance_t
void SetId(int pdgc)
Definition: Target.cxx:149
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:336
bool BareXSecPreCalc(void) const
Definition: RunOpt.h:53
#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
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
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
int ProbePdg(void) const
Definition: InitialState.h:64
Kinematical phase space.
Definition: KPhaseSpace.h:33
const int kPdgTgtFreeN
Definition: PDGCodes.h:200
const int kPdgTgtFreeP
Definition: PDGCodes.h:199
double func(double x, double y)
int Z(void) const
Definition: Target.h:68
#define pINFO
Definition: Messenger.h:62
bool fUsePauliBlocking
account for Pauli blocking?
Resonance_t Resonance(void) const
Definition: XclsTag.h:69
#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
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
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:67
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
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 -&gt; string
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
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
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
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...
enum genie::EInteractionType InteractionType_t
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
double fGSLRelTol
required relative tolerance (error)