GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ReinSehgalRESXSecWithCache.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 <sstream>
12 
13 #include <TMath.h>
14 #include <Math/IFunction.h>
15 #include <Math/IntegratorMultiDim.h>
16 
18 #include "Framework/Conventions/GBuild.h"
28 #include "Framework/Utils/Cache.h"
33 
34 using std::ostringstream;
35 
36 using namespace genie;
37 using namespace genie::controls;
38 using namespace genie::constants;
39 //using namespace genie::units;
40 
41 //____________________________________________________________________________
44 {
45 
46 }
47 //____________________________________________________________________________
50 {
51 
52 }
53 //____________________________________________________________________________
55 XSecIntegratorI(nm,conf)
56 {
57 
58 }
59 //____________________________________________________________________________
61 {
62 
63 }
64 //____________________________________________________________________________
66  const Interaction * in) const
67 {
68 // Cache resonance neutrino production data from free nucleons
69 
70  Cache * cache = Cache::Instance();
71 
72  assert(fSingleResXSecModel);
73 // assert(fIntegrator);
74 
75  // Compute the number of spline knots - use at least 10 knots per decade
76  // && at least 40 knots in the full energy range
77  const double Emin = 0.01;
78  const int nknots_min = (int) (10*(TMath::Log(fEMax)-TMath::Log(Emin)));
79  const int nknots = TMath::Max(40, nknots_min);
80  double * E = new double[nknots]; // knot 'x'
81 
82  TLorentzVector p4(0,0,0,0);
83 
84  int nu_code = in->InitState().ProbePdg();
85  int nuc_code = in->InitState().Tgt().HitNucPdg();
86  int tgt_code = (nuc_code==kPdgProton) ? kPdgTgtFreeP : kPdgTgtFreeN;
87 
88  Interaction * interaction = new Interaction(*in);
89  interaction->InitStatePtr()->SetPdgs(tgt_code, nu_code);
90  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(nuc_code);
91 
92  InteractionType_t wkcur = interaction->ProcInfo().InteractionTypeId();
93  unsigned int nres = fResList.NResonances();
94  for(unsigned int ires = 0; ires < nres; ires++) {
95 
96  // Get next resonance from the resonance list
97  Resonance_t res = fResList.ResonanceId(ires);
98 
99  interaction->ExclTagPtr()->SetResonance(res);
100 
101  // Get a unique cache branch name
102  string key = this->CacheBranchName(res, wkcur, nu_code, nuc_code);
103 
104  // Make sure the cache branch does not already exists
105  CacheBranchFx * cache_branch =
106  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
107  assert(!cache_branch);
108 
109  // Create the new cache branch
110  LOG("ReinSehgalResC", pNOTICE)
111  << "\n ** Creating cache branch - key = " << key;
112  cache_branch = new CacheBranchFx("RES Excitation XSec");
113  cache->AddCacheBranch(key, cache_branch);
114  assert(cache_branch);
115 
116  const KPhaseSpace & kps = interaction->PhaseSpace();
117  double Ethr = kps.Threshold();
118  LOG("ReinSehgalResC", pNOTICE) << "E threshold = " << Ethr;
119 
120  // Distribute the knots in the energy range as is being done in the
121  // XSecSplineList so that the energy threshold is treated correctly
122  // in the spline - see comments there in.
123  int nkb = (Ethr>Emin) ? 5 : 0; // number of knots < threshold
124  int nka = nknots-nkb; // number of knots >= threshold
125  // knots < energy threshold
126  double dEb = (Ethr>Emin) ? (Ethr - Emin) / nkb : 0;
127  for(int i=0; i<nkb; i++) {
128  E[i] = Emin + i*dEb;
129  }
130  // knots >= energy threshold
131  double E0 = TMath::Max(Ethr,Emin);
132  double dEa = (TMath::Log10(fEMax) - TMath::Log10(E0)) /(nka-1);
133  for(int i=0; i<nka; i++) {
134  E[i+nkb] = TMath::Power(10., TMath::Log10(E0) + i * dEa);
135  }
136 
137  // Compute cross sections at the given set of energies
138  for(int ie=0; ie<nknots; ie++) {
139  double xsec = 0.;
140  double Ev = E[ie];
141  p4.SetPxPyPzE(0,0,Ev,Ev);
142  interaction->InitStatePtr()->SetProbeP4(p4);
143 
144  if(Ev>Ethr+kASmallNum) {
145  // Get W integration range and the wider possible Q2 range
146  // (for all W)
147  Range1D_t rW = kps.Limits(kKVW);
148  Range1D_t rQ2 = kps.Limits(kKVQ2);
149 
150  LOG("ReinSehgalResC", pINFO)
151  << "*** Integrating d^2 XSec/dWdQ^2 for R: "
152  << utils::res::AsString(res) << " at Ev = " << Ev;
153  LOG("ReinSehgalResC", pINFO)
154  << "{W} = " << rW.min << ", " << rW.max;
155  LOG("ReinSehgalResC", pINFO)
156  << "{Q^2} = " << rQ2.min << ", " << rQ2.max;
157 
158  if(rW.max<rW.min || rQ2.max<rQ2.min || rW.min<0 || rQ2.min<0) {
159  LOG("ReinSehgalResC", pINFO)
160  << "** Not allowed kinematically, xsec=0";
161  } else {
162 
163  ROOT::Math::IBaseFunctionMultiDim * func =
165  ROOT::Math::IntegrationMultiDim::Type ig_type =
167  ROOT::Math::IntegratorMultiDim ig(ig_type,0,fGSLRelTol,fGSLMaxEval);
168  ig.SetFunction(*func);
169  double kine_min[2] = { rW.min, rQ2.min };
170  double kine_max[2] = { rW.max, rQ2.max };
171  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
172  delete func;
173  }
174  } else {
175  LOG("ReinSehgalResC", pINFO)
176  << "** Below threshold E = " << Ev << " <= " << Ethr;
177  }
178  cache_branch->AddValues(Ev,xsec);
179  SLOG("ReinSehgalResC", pNOTICE)
180  << "RES XSec (R:" << utils::res::AsString(res)
181  << ", E="<< Ev << ") = "<< xsec/(1E-38 *genie::units::cm2) << " x 1E-38 cm^2";
182  }//spline knots
183 
184  // Build the spline
185  cache_branch->CreateSpline();
186  }//ires
187 
188  delete [] E;
189  delete interaction;
190 }
191 //____________________________________________________________________________
193  Resonance_t res, InteractionType_t it, int nupdgc, int nucleonpdgc) const
194 {
195 // Build a unique name for the cache branch
196 
197  Cache * cache = Cache::Instance();
198  string res_name = utils::res::AsString(res);
199  string it_name = InteractionType::AsString(it);
200  string nc_nuc = ((nucleonpdgc==kPdgProton) ? "p" : "n");
201 
202  ostringstream intk;
203  intk << "ResExcitationXSec/R:" << res_name << ";nu:" << nupdgc
204  << ";int:" << it_name << nc_nuc;
205 
206  string algkey = fSingleResXSecModel->Id().Key();
207  string ikey = intk.str();
208  string key = cache->CacheBranchKey(algkey, ikey);
209 
210  return key;
211 }
212 //____________________________________________________________________________
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
InteractionType_t InteractionTypeId(void) const
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
void CreateSpline(string type="TSpline3")
enum genie::EResonance Resonance_t
unsigned int NResonances(void) const
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
void SetResonance(Resonance_t res)
Definition: XclsTag.cxx:128
void AddCacheBranch(string key, CacheBranchI *branch)
Definition: Cache.cxx:88
Summary information for an interaction.
Definition: Interaction.h:56
void AddValues(double x, double y)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double cm2
Definition: Units.h:69
string CacheBranchName(Resonance_t r, InteractionType_t it, int nu, int nuc) const
string CacheBranchKey(string k0, string k1="", string k2="") const
Definition: Cache.cxx:93
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)
void SetPdgs(int tgt_pdgc, int probe_pdgc)
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
void CacheResExcitationXSec(const Interaction *interaction) const
CacheBranchI * FindCacheBranch(string key)
finding/adding cache branches
Definition: Cache.cxx:80
static string AsString(InteractionType_t type)
GENIE Cache Memory.
Definition: Cache.h:38
XclsTag * ExclTagPtr(void) const
Definition: Interaction.h:77
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
double max
Definition: Range1.h:53
int fGSLMaxEval
GSL max evaluations.
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:171
Target * TgtPtr(void) const
Definition: InitialState.h:67
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
const int kPdgProton
Definition: PDGCodes.h:81
#define pNOTICE
Definition: Messenger.h:61
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
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
enum genie::EInteractionType InteractionType_t
string Key(void) const
Definition: AlgId.h:46
double fGSLRelTol
required relative tolerance (error)
Resonance_t ResonanceId(unsigned int ires) const