GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ReinSehgalRESXSecWithCacheFast.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 <sstream>
16 #include <cassert>
17 
18 #include <TMath.h>
19 #include <Math/IFunction.h>
20 #include <Math/IntegratorMultiDim.h>
21 
24 #include "Framework/Conventions/GBuild.h"
38 #include "Framework/Utils/Cache.h"
41 #include "Framework/Utils/Range1.h"
42 
43 
44 
45 
46 
47 
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 
87  assert(fSingleResXSecModel);
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(100, 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("ReinSehgalResCF", 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("ReinSehgalResCF", 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 integration ranges
161  Range1D_t rW = Range1D_t(0.0,1.0);
162  Range1D_t rQ2 = Range1D_t(0.0,1.0);
163 
164  LOG("ReinSehgalResCF", pINFO)
165  << "*** Integrating d^2 XSec/dWdQ^2 for R: "
166  << utils::res::AsString(res) << " at Ev = " << Ev;
167  LOG("ReinSehgalResCF", pINFO)
168  << "{W} = " << rW.min << ", " << rW.max;
169  LOG("ReinSehgalResCF", pINFO)
170  << "{Q^2} = " << rQ2.min << ", " << rQ2.max;
171 
172  if(rW.max<rW.min || rQ2.max<rQ2.min || rW.min<0 || rQ2.min<0) {
173  LOG("ReinSehgalResCF", pINFO)
174  << "** Not allowed kinematically, xsec=0";
175  } else {
176  ROOT::Math::IBaseFunctionMultiDim * func= new utils::gsl::d2XSecRESFast_dWQ2_E(fSingleResXSecModel, interaction);
177  ROOT::Math::IntegrationMultiDim::Type ig_type = utils::gsl::IntegrationNDimTypeFromString(fGSLIntgType);
178  ROOT::Math::IntegratorMultiDim ig(ig_type,0,fGSLRelTol,fGSLMaxEval);
179  ig.SetFunction(*func);
180  double kine_min[2] = { rW.min, rQ2.min };
181  double kine_max[2] = { rW.max, rQ2.max };
182  xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2);
183  delete func;
184  }
185  } else {
186  LOG("ReinSehgalResCF", pINFO)
187  << "** Below threshold E = " << Ev << " <= " << Ethr;
188  }
189  cache_branch->AddValues(Ev,xsec);
190  SLOG("ReinSehgalResCF", pNOTICE)
191  << "RES XSec (R:" << utils::res::AsString(res)
192  << ", E="<< Ev << ") = "<< xsec/(1E-38 *genie::units::cm2)
193  << " x 1E-38 cm^2";
194  }//spline knots
195 
196  // Build the spline
197  cache_branch->CreateSpline();
198  }//ires
199 
200  delete [] E;
201  delete interaction;
202 }
203 //____________________________________________________________________________
205  Resonance_t res, InteractionType_t it, int nupdgc, int nucleonpdgc) const
206 {
207 // Build a unique name for the cache branch
208 
209  Cache * cache = Cache::Instance();
210  string res_name = utils::res::AsString(res);
211  string it_name = InteractionType::AsString(it);
212  string nc_nuc = ((nucleonpdgc==kPdgProton) ? "p" : "n");
213 
214  ostringstream intk;
215  intk << "ResExcitationXSec/R:" << res_name << ";nu:" << nupdgc
216  << ";int:" << it_name << nc_nuc;
217 
218  string algkey = fSingleResXSecModel->Id().Key();
219  string ikey = intk.str();
220  string key = cache->CacheBranchKey(algkey, ikey);
221 
222  return key;
223 }
224 //____________________________________________________________________________
225 // GSL wrappers
226 //____________________________________________________________________________
228  const XSecAlgorithmI * m, const Interaction * i) :
229 ROOT::Math::IBaseFunctionMultiDim(),
230 fModel(m),
231 fInteraction(i)
232 {
234  Range1D_t Wl = kps->WLim();
235  fWmin = Wl.min;
236  fWmax = Wl.max;
237  Registry fConfig = (const_cast<XSecAlgorithmI *>(fModel))->GetConfig();
238  bool fUsingDisResJoin = fConfig.GetBool("UseDRJoinScheme");
239  double fWcut = 999999;
240  if(fUsingDisResJoin)
241  {
242  fWcut = fConfig.GetDouble("Wcut");
243  }
244  fWmax=TMath::Min(fWcut, fWmax);
245  if (fWcut<fWmin)
246  isfWcutLessfWmin=true;
247  else
248  isfWcutLessfWmin=false;
249  bool fNormBW = fConfig.GetBoolDef("BreitWignerNorm", true);
250  if (fNormBW)
251  {
252  double fN2ResMaxNWidths = fConfig.GetDoubleDef("MaxNWidthForN2Res", 2.0);
253  double fN0ResMaxNWidths = fConfig.GetDoubleDef("MaxNWidthForN0Res", 6.0);
254  double fGNResMaxNWidths = fConfig.GetDoubleDef("MaxNWidthForGNRes", 4.0);
255  Resonance_t resonance = fInteraction->ExclTag().Resonance();
256  int IR = utils::res::ResonanceIndex (resonance);
257  double MR = utils::res::Mass (resonance);
258  double WR = utils::res::Width (resonance);
259  if (IR==0)
260  fWcut = MR + fN0ResMaxNWidths * WR;
261  else if (IR==2)
262  fWcut = MR + fN2ResMaxNWidths * WR;
263  else
264  fWcut = MR + fGNResMaxNWidths * WR;
265  fWmax=TMath::Min(fWcut, fWmax);
266  if (fWcut<fWmin)
267  isfWcutLessfWmin=true;
268  }
269 
270 }
272 {
273 
274 }
276 {
277  return 2;
278 }
279 double genie::utils::gsl::d2XSecRESFast_dWQ2_E::DoEval(const double * xin) const
280 {
281 // inputs:
282 // W
283 // Q2
284 // outputs:
285 // differential cross section [10^-38 cm^2/GeV^3] for Resonance production
286 //
287 
288  if (isfWcutLessfWmin)
289  return 0;
290  double W = fWmin+(fWmax-fWmin)*xin[0];
291  fInteraction->KinePtr()->SetW(W);
292  Range1D_t Q2l = kps->Q2Lim_W();
293  if (Q2l.min<0 || Q2l.max<0)
294  return 0.0;
295  double Q2 = Q2l.min+(Q2l.max-Q2l.min)*xin[1];
296  fInteraction->KinePtr()->SetQ2(Q2);
297  double xsec = fModel->XSec(fInteraction, kPSWQ2fE)*(fWmax-fWmin)*(Q2l.max-Q2l.min);
298  return xsec/(1E-38 * units::cm2);
299 }
300 ROOT::Math::IBaseFunctionMultiDim *
302 {
303  return
304  new genie::utils::gsl::d2XSecRESFast_dWQ2_E(fModel,fInteraction);
305 }
306 //____________________________________________________________________________
Cross Section Calculation Interface.
RgDbl GetDoubleDef(RgKey key, RgDbl def_opt, bool set_def=true)
Definition: Registry.cxx:535
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.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
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")
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:474
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:128
void AddCacheBranch(string key, CacheBranchI *branch)
Definition: Cache.cxx:88
KPhaseSpace * PhaseSpacePtr(void) const
Definition: Interaction.h:78
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
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 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
Resonance_t Resonance(void) const
Definition: XclsTag.h:69
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.
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:171
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
Target * TgtPtr(void) const
Definition: InitialState.h:67
RgBool GetBool(RgKey key) const
Definition: Registry.cxx:460
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
static constexpr double m
Definition: Units.h:71
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:525
d2XSecRESFast_dWQ2_E(const XSecAlgorithmI *m, const Interaction *i)
string Key(void) const
Definition: AlgId.h:46
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