ReinSehgalRESXSecWithCache.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 - March 09, 2006
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Jan 18, 2008 - CA
14  Simplify the way free nucleon channels (for which we cache cross sections)
15  are built from the input interaction
16  @ Jan 19, 2008 - CA
17  Modify the way knots are distributed in the cached free nucleon resonance
18  neutrino production splines so that the energy threshold is treated more
19  accurately (see also XSecSplineList.cxx).
20  @ Sep 07, 2009 - CA
21  Integrated with GNU Numerical Library (GSL) via ROOT's MathMore library.
22 
23 */
24 //____________________________________________________________________________
25 
26 #include <sstream>
27 
28 #include <TMath.h>
29 #include <Math/IFunction.h>
30 #include <Math/IntegratorMultiDim.h>
31 
43 #include "Framework/Utils/Cache.h"
48 
49 using std::ostringstream;
50 
51 using namespace genie;
52 using namespace genie::controls;
53 using namespace genie::constants;
54 //using namespace genie::units;
55 
56 //____________________________________________________________________________
59 {
60 
61 }
62 //____________________________________________________________________________
65 {
66 
67 }
68 //____________________________________________________________________________
70 XSecIntegratorI(nm,conf)
71 {
72 
73 }
74 //____________________________________________________________________________
76 {
77 
78 }
79 //____________________________________________________________________________
81  const Interaction * in) const
82 {
83 // Cache resonance neutrino production data from free nucleons
84 
85  Cache * cache = Cache::Instance();
86 
88 // assert(fIntegrator);
89 
90  // Compute the number of spline knots - use at least 10 knots per decade
91  // && at least 40 knots in the full energy range
92  const double Emin = 0.01;
93  const int nknots_min = (int) (10*(TMath::Log(fEMax)-TMath::Log(Emin)));
94  const int nknots = TMath::Max(40, nknots_min);
95  double * E = new double[nknots]; // knot 'x'
96 
97  TLorentzVector p4(0,0,0,0);
98 
99  int nu_code = in->InitState().ProbePdg();
100  int nuc_code = in->InitState().Tgt().HitNucPdg();
101  int tgt_code = (nuc_code==kPdgProton) ? kPdgTgtFreeP : kPdgTgtFreeN;
102 
103  Interaction * interaction = new Interaction(*in);
104  interaction->InitStatePtr()->SetPdgs(tgt_code, nu_code);
105  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(nuc_code);
106 
107  InteractionType_t wkcur = interaction->ProcInfo().InteractionTypeId();
108  unsigned int nres = fResList.NResonances();
109  for(unsigned int ires = 0; ires < nres; ires++) {
110 
111  // Get next resonance from the resonance list
112  Resonance_t res = fResList.ResonanceId(ires);
113 
114  interaction->ExclTagPtr()->SetResonance(res);
115 
116  // Get a unique cache branch name
117  string key = this->CacheBranchName(res, wkcur, nu_code, nuc_code);
118 
119  // Make sure the cache branch does not already exists
120  CacheBranchFx * cache_branch =
121  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
122  assert(!cache_branch);
123 
124  // Create the new cache branch
125  LOG("ReinSehgalResC", pNOTICE)
126  << "\n ** Creating cache branch - key = " << key;
127  cache_branch = new CacheBranchFx("RES Excitation XSec");
128  cache->AddCacheBranch(key, cache_branch);
129  assert(cache_branch);
130 
131  const KPhaseSpace & kps = interaction->PhaseSpace();
132  double Ethr = kps.Threshold();
133  LOG("ReinSehgalResC", pNOTICE) << "E threshold = " << Ethr;
134 
135  // Distribute the knots in the energy range as is being done in the
136  // XSecSplineList so that the energy threshold is treated correctly
137  // in the spline - see comments there in.
138  int nkb = (Ethr>Emin) ? 5 : 0; // number of knots < threshold
139  int nka = nknots-nkb; // number of knots >= threshold
140  // knots < energy threshold
141  double dEb = (Ethr>Emin) ? (Ethr - Emin) / nkb : 0;
142  for(int i=0; i<nkb; i++) {
143  E[i] = Emin + i*dEb;
144  }
145  // knots >= energy threshold
146  double E0 = TMath::Max(Ethr,Emin);
147  double dEa = (TMath::Log10(fEMax) - TMath::Log10(E0)) /(nka-1);
148  for(int i=0; i<nka; i++) {
149  E[i+nkb] = TMath::Power(10., TMath::Log10(E0) + i * dEa);
150  }
151 
152  // Compute cross sections at the given set of energies
153  for(int ie=0; ie<nknots; ie++) {
154  double xsec = 0.;
155  double Ev = E[ie];
156  p4.SetPxPyPzE(0,0,Ev,Ev);
157  interaction->InitStatePtr()->SetProbeP4(p4);
158 
159  if(Ev>Ethr+kASmallNum) {
160  // Get W integration range and the wider possible Q2 range
161  // (for all W)
162  Range1D_t rW = kps.Limits(kKVW);
163  Range1D_t rQ2 = kps.Limits(kKVQ2);
164 
165  LOG("ReinSehgalResC", pINFO)
166  << "*** Integrating d^2 XSec/dWdQ^2 for R: "
167  << utils::res::AsString(res) << " at Ev = " << Ev;
168  LOG("ReinSehgalResC", pINFO)
169  << "{W} = " << rW.min << ", " << rW.max;
170  LOG("ReinSehgalResC", pINFO)
171  << "{Q^2} = " << rQ2.min << ", " << rQ2.max;
172 
173  if(rW.max<rW.min || rQ2.max<rQ2.min || rW.min<0 || rQ2.min<0) {
174  LOG("ReinSehgalResC", pINFO)
175  << "** Not allowed kinematically, xsec=0";
176  } else {
177 
178  ROOT::Math::IBaseFunctionMultiDim * func =
182  ROOT::Math::IntegratorMultiDim ig(ig_type,0,fGSLRelTol,fGSLMaxEval);
183  ig.SetFunction(*func);
184  double kine_min[2] = { rW.min, rQ2.min };
185  double kine_max[2] = { rW.max, rQ2.max };
186  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
187  delete func;
188  }
189  } else {
190  LOG("ReinSehgalResC", pINFO)
191  << "** Below threshold E = " << Ev << " <= " << Ethr;
192  }
193  cache_branch->AddValues(Ev,xsec);
194  SLOG("ReinSehgalResC", pNOTICE)
195  << "RES XSec (R:" << utils::res::AsString(res)
196  << ", E="<< Ev << ") = "<< xsec/(1E-38 *genie::units::cm2) << " x 1E-38 cm^2";
197  }//spline knots
198 
199  // Build the spline
200  cache_branch->CreateSpline();
201  }//ires
202 
203  delete [] E;
204  delete interaction;
205 }
206 //____________________________________________________________________________
208  Resonance_t res, InteractionType_t it, int nupdgc, int nucleonpdgc) const
209 {
210 // Build a unique name for the cache branch
211 
212  Cache * cache = Cache::Instance();
213  string res_name = utils::res::AsString(res);
214  string it_name = InteractionType::AsString(it);
215  string nc_nuc = ((nucleonpdgc==kPdgProton) ? "p" : "n");
216 
217  ostringstream intk;
218  intk << "ResExcitationXSec/R:" << res_name << ";nu:" << nupdgc
219  << ";int:" << it_name << nc_nuc;
220 
221  string algkey = fSingleResXSecModel->Id().Key();
222  string ikey = intk.str();
223  string key = cache->CacheBranchKey(algkey, ikey);
224 
225  return key;
226 }
227 //____________________________________________________________________________
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
Basic constants.
InteractionType_t InteractionTypeId(void) const
string fGSLIntgType
name of GSL numerical integrator
void SetProbeP4(const TLorentzVector &P4)
set< int >::iterator it
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
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
double Threshold(void) const
Energy threshold.
Definition: KPhaseSpace.cxx:81
static constexpr Double_t nm
Definition: Munits.h:133
static const double cm2
Definition: Units.h:77
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:123
void AddCacheBranch(string key, CacheBranchI *branch)
Definition: Cache.cxx:97
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:97
const Int_t nknots
Definition: testXsec.C:24
Float_t E
Definition: plot.C:20
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:102
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)
void SetPdgs(int tgt_pdgc, int probe_pdgc)
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:63
Misc GENIE control constants.
void CacheResExcitationXSec(const Interaction *interaction) const
Double_t xsec[nknots]
Definition: testXsec.C:47
CacheBranchI * FindCacheBranch(string key)
finding/adding cache branches
Definition: Cache.cxx:89
static string AsString(InteractionType_t type)
GENIE Cache Memory.
Definition: Cache.h:39
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:54
int fGSLMaxEval
GSL max evaluations.
ifstream in
Definition: comparison.C:7
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:188
Target * TgtPtr(void) const
Definition: InitialState.h:68
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
const int kPdgProton
Definition: PDGCodes.h:65
#define pNOTICE
Definition: Messenger.h:62
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
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:47
double fGSLRelTol
required relative tolerance (error)
Resonance_t ResonanceId(unsigned int ires) const