GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ReinSehgalRESXSec.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 
17 #include "Framework/Conventions/GBuild.h"
25 #include "Framework/Utils/RunOpt.h"
28 #include "Framework/Utils/Cache.h"
34 
35 using namespace genie;
36 using namespace genie::constants;
37 //using namespace genie::units;
38 
39 //____________________________________________________________________________
41 ReinSehgalRESXSecWithCache("genie::ReinSehgalRESXSec")
42 {
43 
44 }
45 //____________________________________________________________________________
47 ReinSehgalRESXSecWithCache("genie::ReinSehgalRESXSec", config)
48 {
49 
50 }
51 //____________________________________________________________________________
53 {
54 
55 }
56 //____________________________________________________________________________
58  const XSecAlgorithmI * model, const Interaction * interaction) const
59 {
60  if(! model->ValidProcess(interaction) ) return 0.;
61  fSingleResXSecModel = model;
62 
63  const KPhaseSpace & kps = interaction->PhaseSpace();
64  if(!kps.IsAboveThreshold()) {
65  LOG("ReinSehgalRESXSec", pDEBUG) << "*** Below energy threshold";
66  return 0;
67  }
68 
69  //-- Get init state and process information
70  const InitialState & init_state = interaction->InitState();
71  const ProcessInfo & proc_info = interaction->ProcInfo();
72  const Target & target = init_state.Tgt();
73 
74  InteractionType_t it = proc_info.InteractionTypeId();
75  int nucleon_pdgc = target.HitNucPdg();
76  int nu_pdgc = init_state.ProbePdg();
77 
78  //-- Get neutrino energy in the struck nucleon rest frame
79  double Ev = init_state.ProbeE(kRfHitNucRest);
80 
81  //-- Get the requested resonance
82  Resonance_t res = interaction->ExclTag().Resonance();
83 
84  // If the input interaction is off a nuclear target, then chek whether
85  // the corresponding free nucleon cross section already exists at the
86  // cross section spline list.
87  // If yes, calculate the nuclear cross section based on that value.
88  //
90  if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() ) {
91  Interaction * in = new Interaction(*interaction);
92  if(pdg::IsProton(nucleon_pdgc)) {
94  } else {
96  }
97  if(xsl->SplineExists(model,in)) {
98  const Spline * spl = xsl->GetSpline(model, in);
99  double xsec = spl->Evaluate(Ev);
100  SLOG("ReinSehgalResT", pNOTICE)
101  << "XSec[RES/" << utils::res::AsString(res)<< "/free] (Ev = "
102  << Ev << " GeV) = " << xsec/(1E-38 *genie::units::cm2) << " x 1E-38 cm^2";
103  if(! interaction->TestBit(kIAssumeFreeNucleon) ) {
104  int NNucl = (pdg::IsProton(nucleon_pdgc)) ? target.Z() : target.N();
105  xsec *= NNucl;
106  }
107  delete in;
108  return xsec;
109  }
110  delete in;
111  }
112 
113  // There was no corresponding free nucleon spline saved in XSecSplineList that
114  // could be used to speed up this calculation.
115  // Check whether local caching of free nucleon cross sections is allowed.
116  // If yes, store free nucleon cross sections at a cache branch and use those
117  // at any subsequent call.
118  //
119  bool bare_xsec_pre_calc = RunOpt::Instance()->BareXSecPreCalc();
120  if(bare_xsec_pre_calc && !fUsePauliBlocking) {
121  Cache * cache = Cache::Instance();
122  string key = this->CacheBranchName(res, it, nu_pdgc, nucleon_pdgc);
123  LOG("ReinSehgalResT", pINFO)
124  << "Finding cache branch with key: " << key;
125  CacheBranchFx * cache_branch =
126  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
127  if(!cache_branch) {
128  LOG("ReinSehgalResT", pWARN)
129  << "No cached RES v-production data for input neutrino"
130  << " (pdgc: " << nu_pdgc << ")";
131  LOG("ReinSehgalResT", pWARN)
132  << "Wait while computing/caching RES production xsec first...";
133 
134  this->CacheResExcitationXSec(interaction);
135 
136  LOG("ReinSehgalResT", pINFO) << "Done caching resonance xsec data";
137  LOG("ReinSehgalResT", pINFO)
138  << "Finding newly created cache branch with key: " << key;
139  cache_branch =
140  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
141  assert(cache_branch);
142  }
143  const CacheBranchFx & cbranch = (*cache_branch);
144 
145  //-- Get cached resonance neutrinoproduction xsec
146  // (If E>Emax, assume xsec = xsec(Emax) - but do not evaluate the
147  // cross section spline at the end of its energy range-)
148  double rxsec = (Ev<fEMax-1) ? cbranch(Ev) : cbranch(fEMax-1);
149 
150  SLOG("ReinSehgalResT", pNOTICE)
151  << "XSec[RES/" << utils::res::AsString(res)<< "/free] (Ev = "
152  << Ev << " GeV) = " << rxsec/(1E-38 *genie::units::cm2) << " x 1E-38 cm^2";
153 
154  if( interaction->TestBit(kIAssumeFreeNucleon) ) return rxsec;
155 
156  int NNucl = (pdg::IsProton(nucleon_pdgc)) ? target.Z() : target.N();
157  rxsec*=NNucl; // nuclear xsec
158  return rxsec;
159  } // disable local caching
160 
161  // Just go ahead and integrate the input differential cross section for the
162  // specified interaction.
163  else {
164 
165  Range1D_t rW = kps.Limits(kKVW);
166  Range1D_t rQ2 = kps.Limits(kKVQ2);
167 
168  LOG("ReinSehgalResC", pINFO)
169  << "*** Integrating d^2 XSec/dWdQ^2 for R: "
170  << utils::res::AsString(res) << " at Ev = " << Ev;
171  LOG("ReinSehgalResC", pINFO)
172  << "{W} = " << rW.min << ", " << rW.max;
173  LOG("ReinSehgalResC", pINFO)
174  << "{Q^2} = " << rQ2.min << ", " << rQ2.max;
175 
176  ROOT::Math::IBaseFunctionMultiDim * func =
177  new utils::gsl::d2XSec_dWdQ2_E(model, interaction);
178  ROOT::Math::IntegrationMultiDim::Type ig_type =
180  ROOT::Math::IntegratorMultiDim ig(ig_type,0,fGSLRelTol,fGSLMaxEval);
181  ig.SetFunction(*func);
182  double kine_min[2] = { rW.min, rQ2.min };
183  double kine_max[2] = { rW.max, rQ2.max };
184  double xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
185 
186  delete func;
187  return xsec;
188  }
189  return 0;
190 }
191 //____________________________________________________________________________
193 {
194  Algorithm::Configure(config);
195  this->LoadConfig();
196 }
197 //____________________________________________________________________________
198 void ReinSehgalRESXSec::Configure(string config)
199 {
200  Algorithm::Configure(config);
201  this->LoadConfig();
202 }
203 //____________________________________________________________________________
205 {
206  // Get GSL integration type & relative tolerance
207  GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ;
208  GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 0.01 ) ;
209  GetParamDef( "gsl-max-eval", fGSLMaxEval, 100000 ) ;
210  GetParam("UsePauliBlockingForRES", fUsePauliBlocking);
211  // Get upper E limit on res xsec spline (=f(E)) before assuming xsec=const
212  GetParamDef( "ESplineMax", fEMax, 100. ) ;
213  fEMax = TMath::Max(fEMax, 20.); // don't accept user Emax if less than 20 GeV
214 
215  // Create the baryon resonance list specified in the config.
216  fResList.Clear();
217  string resonances ;
218  GetParam( "ResonanceNameList", resonances ) ;
219  fResList.DecodeFromNameList(resonances);
220 
221 }
222 //____________________________________________________________________________
Cross Section Calculation Interface.
bool fUsePauliBlocking
account for Pauli blocking?
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
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=",")
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
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: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
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: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
Resonance_t Resonance(void) const
Definition: XclsTag.h:69
void CacheResExcitationXSec(const Interaction *interaction) const
#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
double Integrate(const XSecAlgorithmI *model, const Interaction *i) const
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
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:67
bool IsAboveThreshold(void) const
Checks whether the interaction is above the energy threshold.
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 -&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)