ReinSehgalRESXSecWithCacheFast.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 #include <sstream>
18 #include <cassert>
19 
20 #include <TMath.h>
21 #include <Math/IFunction.h>
22 #include <Math/IntegratorMultiDim.h>
23 
40 #include "Framework/Utils/Cache.h"
43 #include "Framework/Utils/Range1.h"
44 
45 
46 
47 
48 
49 
50 
51 using std::ostringstream;
52 
53 using namespace genie;
54 using namespace genie::controls;
55 using namespace genie::constants;
56 //using namespace genie::units;
57 
58 //____________________________________________________________________________
61 {
62 
63 }
64 //____________________________________________________________________________
67 {
68 
69 }
70 //____________________________________________________________________________
72 XSecIntegratorI(nm,conf)
73 {
74 
75 }
76 //____________________________________________________________________________
78 {
79 
80 }
81 //____________________________________________________________________________
83  const Interaction * in) const
84 {
85 // Cache resonance neutrino production data from free nucleons
86 
87  Cache * cache = Cache::Instance();
88 
90 // assert(fIntegrator);
91 
92  // Compute the number of spline knots - use at least 10 knots per decade
93  // && at least 40 knots in the full energy range
94  const double Emin = 0.01;
95  const int nknots_min = (int) (10*(TMath::Log(fEMax)-TMath::Log(Emin)));
96  const int nknots = TMath::Max(100, nknots_min);
97  double * E = new double[nknots]; // knot 'x'
98 
99  TLorentzVector p4(0,0,0,0);
100 
101  int nu_code = in->InitState().ProbePdg();
102  int nuc_code = in->InitState().Tgt().HitNucPdg();
103  int tgt_code = (nuc_code==kPdgProton) ? kPdgTgtFreeP : kPdgTgtFreeN;
104 
105  Interaction * interaction = new Interaction(*in);
106  interaction->InitStatePtr()->SetPdgs(tgt_code, nu_code);
107  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(nuc_code);
108 
109  InteractionType_t wkcur = interaction->ProcInfo().InteractionTypeId();
110  unsigned int nres = fResList.NResonances();
111  for(unsigned int ires = 0; ires < nres; ires++) {
112 
113  // Get next resonance from the resonance list
114  Resonance_t res = fResList.ResonanceId(ires);
115 
116  interaction->ExclTagPtr()->SetResonance(res);
117 
118  // Get a unique cache branch name
119  string key = this->CacheBranchName(res, wkcur, nu_code, nuc_code);
120 
121  // Make sure the cache branch does not already exists
122  CacheBranchFx * cache_branch =
123  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
124  assert(!cache_branch);
125 
126  // Create the new cache branch
127  LOG("ReinSehgalResCF", pNOTICE)
128  << "\n ** Creating cache branch - key = " << key;
129  cache_branch = new CacheBranchFx("RES Excitation XSec");
130  cache->AddCacheBranch(key, cache_branch);
131  assert(cache_branch);
132 
133  const KPhaseSpace & kps = interaction->PhaseSpace();
134  double Ethr = kps.Threshold();
135  LOG("ReinSehgalResCF", pNOTICE) << "E threshold = " << Ethr;
136 
137  // Distribute the knots in the energy range as is being done in the
138  // XSecSplineList so that the energy threshold is treated correctly
139  // in the spline - see comments there in.
140  int nkb = (Ethr>Emin) ? 5 : 0; // number of knots < threshold
141  int nka = nknots-nkb; // number of knots >= threshold
142  // knots < energy threshold
143  double dEb = (Ethr>Emin) ? (Ethr - Emin) / nkb : 0;
144  for(int i=0; i<nkb; i++) {
145  E[i] = Emin + i*dEb;
146  }
147  // knots >= energy threshold
148  double E0 = TMath::Max(Ethr,Emin);
149  double dEa = (TMath::Log10(fEMax) - TMath::Log10(E0)) /(nka-1);
150  for(int i=0; i<nka; i++) {
151  E[i+nkb] = TMath::Power(10., TMath::Log10(E0) + i * dEa);
152  }
153 
154  // Compute cross sections at the given set of energies
155  for(int ie=0; ie<nknots; ie++) {
156  double xsec = 0.;
157  double Ev = E[ie];
158  p4.SetPxPyPzE(0,0,Ev,Ev);
159  interaction->InitStatePtr()->SetProbeP4(p4);
160 
161  if(Ev>Ethr+kASmallNum) {
162  // Get integration ranges
163  Range1D_t rW = Range1D_t(0.0,1.0);
164  Range1D_t rQ2 = Range1D_t(0.0,1.0);
165 
166  LOG("ReinSehgalResCF", pINFO)
167  << "*** Integrating d^2 XSec/dWdQ^2 for R: "
168  << utils::res::AsString(res) << " at Ev = " << Ev;
169  LOG("ReinSehgalResCF", pINFO)
170  << "{W} = " << rW.min << ", " << rW.max;
171  LOG("ReinSehgalResCF", pINFO)
172  << "{Q^2} = " << rQ2.min << ", " << rQ2.max;
173 
174  if(rW.max<rW.min || rQ2.max<rQ2.min || rW.min<0 || rQ2.min<0) {
175  LOG("ReinSehgalResCF", pINFO)
176  << "** Not allowed kinematically, xsec=0";
177  } else {
178  ROOT::Math::IBaseFunctionMultiDim * func= new utils::gsl::d2XSecRESFast_dWQ2_E(fSingleResXSecModel, interaction);
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  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
185  delete func;
186  }
187  } else {
188  LOG("ReinSehgalResCF", pINFO)
189  << "** Below threshold E = " << Ev << " <= " << Ethr;
190  }
191  cache_branch->AddValues(Ev,xsec);
192  SLOG("ReinSehgalResCF", pNOTICE)
193  << "RES XSec (R:" << utils::res::AsString(res)
194  << ", E="<< Ev << ") = "<< xsec/(1E-38 *genie::units::cm2)
195  << " x 1E-38 cm^2";
196  }//spline knots
197 
198  // Build the spline
199  cache_branch->CreateSpline();
200  }//ires
201 
202  delete [] E;
203  delete interaction;
204 }
205 //____________________________________________________________________________
207  Resonance_t res, InteractionType_t it, int nupdgc, int nucleonpdgc) const
208 {
209 // Build a unique name for the cache branch
210 
211  Cache * cache = Cache::Instance();
212  string res_name = utils::res::AsString(res);
213  string it_name = InteractionType::AsString(it);
214  string nc_nuc = ((nucleonpdgc==kPdgProton) ? "p" : "n");
215 
216  ostringstream intk;
217  intk << "ResExcitationXSec/R:" << res_name << ";nu:" << nupdgc
218  << ";int:" << it_name << nc_nuc;
219 
220  string algkey = fSingleResXSecModel->Id().Key();
221  string ikey = intk.str();
222  string key = cache->CacheBranchKey(algkey, ikey);
223 
224  return key;
225 }
226 //____________________________________________________________________________
227 // GSL wrappers
228 //____________________________________________________________________________
230  const XSecAlgorithmI * m, const Interaction * i) :
231 ROOT::Math::IBaseFunctionMultiDim(),
232 fModel(m),
233 fInteraction(i)
234 {
236  Range1D_t Wl = kps->WLim();
237  fWmin = Wl.min;
238  fWmax = Wl.max;
239  Registry fConfig = (const_cast<XSecAlgorithmI *>(fModel))->GetConfig();
240  bool fUsingDisResJoin = fConfig.GetBool("UseDRJoinScheme");
241  double fWcut = 999999;
242  if(fUsingDisResJoin)
243  {
244  fWcut = fConfig.GetDouble("Wcut");
245  }
246  fWmax=TMath::Min(fWcut, fWmax);
247  if (fWcut<fWmin)
248  isfWcutLessfWmin=true;
249  else
250  isfWcutLessfWmin=false;
251  bool fNormBW = fConfig.GetBoolDef("BreitWignerNorm", true);
252  if (fNormBW)
253  {
254  double fN2ResMaxNWidths = fConfig.GetDoubleDef("MaxNWidthForN2Res", 2.0);
255  double fN0ResMaxNWidths = fConfig.GetDoubleDef("MaxNWidthForN0Res", 6.0);
256  double fGNResMaxNWidths = fConfig.GetDoubleDef("MaxNWidthForGNRes", 4.0);
257  Resonance_t resonance = fInteraction->ExclTag().Resonance();
258  int IR = utils::res::ResonanceIndex (resonance);
259  double MR = utils::res::Mass (resonance);
260  double WR = utils::res::Width (resonance);
261  if (IR==0)
262  fWcut = MR + fN0ResMaxNWidths * WR;
263  else if (IR==2)
264  fWcut = MR + fN2ResMaxNWidths * WR;
265  else
266  fWcut = MR + fGNResMaxNWidths * WR;
267  fWmax=TMath::Min(fWcut, fWmax);
268  if (fWcut<fWmin)
269  isfWcutLessfWmin=true;
270  }
271 
272 }
274 {
275 
276 }
278 {
279  return 2;
280 }
281 double genie::utils::gsl::d2XSecRESFast_dWQ2_E::DoEval(const double * xin) const
282 {
283 // inputs:
284 // W
285 // Q2
286 // outputs:
287 // differential cross section [10^-38 cm^2/GeV^3] for Resonance production
288 //
289 
290  if (isfWcutLessfWmin)
291  return 0;
292  double W = fWmin+(fWmax-fWmin)*xin[0];
293  fInteraction->KinePtr()->SetW(W);
294  Range1D_t Q2l = kps->Q2Lim_W();
295  if (Q2l.min<0 || Q2l.max<0)
296  return 0.0;
297  double Q2 = Q2l.min+(Q2l.max-Q2l.min)*xin[1];
298  fInteraction->KinePtr()->SetQ2(Q2);
299  double xsec = fModel->XSec(fInteraction, kPSWQ2fE)*(fWmax-fWmin)*(Q2l.max-Q2l.min);
300  return xsec/(1E-38 * units::cm2);
301 }
302 ROOT::Math::IBaseFunctionMultiDim *
304 {
305  return
307 }
308 //____________________________________________________________________________
Cross Section Calculation Interface.
RgDbl GetDoubleDef(RgKey key, RgDbl def_opt, bool set_def=true)
Definition: Registry.cxx:550
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.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
Definition: GSLUtils.cxx:67
int HitNucPdg(void) const
Definition: Target.cxx:321
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:265
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
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
double Mass(Resonance_t res)
resonance mass (GeV)
void CacheResExcitationXSec(const Interaction *interaction) const
double Width(Resonance_t res)
resonance width (GeV)
RgDbl GetDouble(RgKey key) const
Definition: Registry.cxx:489
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
static const double cm2
Definition: Units.h:77
enum genie::EResonance Resonance_t
string CacheBranchName(Resonance_t r, InteractionType_t it, int nu, int nuc) const
unsigned int NResonances(void) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
void SetResonance(Resonance_t res)
Definition: XclsTag.cxx:123
void AddCacheBranch(string key, CacheBranchI *branch)
Definition: Cache.cxx:97
KPhaseSpace * PhaseSpacePtr(void) const
Definition: Interaction.h:78
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
static keras::KerasModel * fModel
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
Resonance_t Resonance(void) const
Definition: XclsTag.h:62
Misc GENIE control constants.
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
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
int fGSLMaxEval
GSL max evaluations.
ifstream in
Definition: comparison.C:7
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:188
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
Target * TgtPtr(void) const
Definition: InitialState.h:68
RgBool GetBool(RgKey key) const
Definition: Registry.cxx:475
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
#define W(x)
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
enum genie::EInteractionType InteractionType_t
RgBool GetBoolDef(RgKey key, RgBool def_opt, bool set_def=true)
Definition: Registry.cxx:540
d2XSecRESFast_dWQ2_E(const XSecAlgorithmI *m, const Interaction *i)
string Key(void) const
Definition: AlgId.h:47
Range1D_t WLim(void) const
W limits.
int ResonanceIndex(Resonance_t res)
resonance idx, quark model / SU(6)
double fGSLRelTol
required relative tolerance (error)
Resonance_t ResonanceId(unsigned int ires) const