21 #include <Math/IFunction.h> 22 #include <Math/IntegratorMultiDim.h> 51 using std::ostringstream;
53 using namespace genie;
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];
99 TLorentzVector
p4(0,0,0,0);
111 for(
unsigned int ires = 0; ires < nres; ires++) {
128 <<
"\n ** Creating cache branch - key = " <<
key;
135 LOG(
"ReinSehgalResCF",
pNOTICE) <<
"E threshold = " << Ethr;
140 int nkb = (Ethr>Emin) ? 5 : 0;
141 int nka = nknots-nkb;
143 double dEb = (Ethr>Emin) ? (Ethr - Emin) / nkb : 0;
144 for(
int i=0;
i<nkb;
i++) {
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);
155 for(
int ie=0; ie<
nknots; ie++) {
158 p4.SetPxPyPzE(0,0,Ev,Ev);
167 <<
"*** Integrating d^2 XSec/dWdQ^2 for R: " 170 <<
"{W} = " << rW.
min <<
", " << rW.
max;
172 <<
"{Q^2} = " << rQ2.
min <<
", " << rQ2.
max;
176 <<
"** Not allowed kinematically, xsec=0";
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);
189 <<
"** Below threshold E = " << Ev <<
" <= " << Ethr;
214 string nc_nuc = ((nucleonpdgc==
kPdgProton) ?
"p" :
"n");
217 intk <<
"ResExcitationXSec/R:" << res_name <<
";nu:" << nupdgc
218 <<
";int:" << it_name << nc_nuc;
221 string ikey = intk.str();
231 ROOT::Math::IBaseFunctionMultiDim(),
240 bool fUsingDisResJoin = fConfig.
GetBool(
"UseDRJoinScheme");
241 double fWcut = 999999;
251 bool fNormBW = fConfig.
GetBoolDef(
"BreitWignerNorm",
true);
254 double fN2ResMaxNWidths = fConfig.
GetDoubleDef(
"MaxNWidthForN2Res", 2.0);
255 double fN0ResMaxNWidths = fConfig.
GetDoubleDef(
"MaxNWidthForN0Res", 6.0);
256 double fGNResMaxNWidths = fConfig.
GetDoubleDef(
"MaxNWidthForGNRes", 4.0);
262 fWcut = MR + fN0ResMaxNWidths * WR;
264 fWcut = MR + fN2ResMaxNWidths * WR;
266 fWcut = MR + fGNResMaxNWidths * WR;
295 if (Q2l.
min<0 || Q2l.
max<0)
302 ROOT::Math::IBaseFunctionMultiDim *
Cross Section Calculation Interface.
RgDbl GetDoubleDef(RgKey key, RgDbl def_opt, bool set_def=true)
const KPhaseSpace & PhaseSpace(void) const
InteractionType_t InteractionTypeId(void) const
string fGSLIntgType
name of GSL numerical integrator
void SetProbeP4(const TLorentzVector &P4)
THE MAIN GENIE PROJECT NAMESPACE
Cross Section Integrator Interface.
double DoEval(const double *xin) const
double Q2(const Interaction *const i)
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
int HitNucPdg(void) const
void SetQ2(double Q2, bool selected=false)
Kinematics * KinePtr(void) const
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.
double Threshold(void) const
Energy threshold.
static constexpr Double_t nm
const XSecAlgorithmI * fSingleResXSecModel
const Interaction * fInteraction
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
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
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)
void AddCacheBranch(string key, CacheBranchI *branch)
KPhaseSpace * PhaseSpacePtr(void) const
unsigned int NDim(void) const
Summary information for an interaction.
void AddValues(double x, double y)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
static keras::KerasModel * fModel
string CacheBranchKey(string k0, string k1="", string k2="") const
double func(double x, double y)
void SetPdgs(int tgt_pdgc, int probe_pdgc)
static const double kASmallNum
virtual ~ReinSehgalRESXSecWithCacheFast()
Resonance_t Resonance(void) const
Misc GENIE control constants.
CacheBranchI * FindCacheBranch(string key)
finding/adding cache branches
ReinSehgalRESXSecWithCacheFast()
static string AsString(InteractionType_t type)
XclsTag * ExclTagPtr(void) const
virtual const AlgId & Id(void) const
Get algorithm ID.
void SetW(double W, bool selected=false)
int fGSLMaxEval
GSL max evaluations.
A registry. Provides the container for algorithm configuration parameters.
void SetHitNucPdg(int pdgc)
const XclsTag & ExclTag(void) const
Target * TgtPtr(void) const
RgBool GetBool(RgKey key) const
assert(nhit_max >=nhit_nbins)
InitialState * InitStatePtr(void) const
const InitialState & InitState(void) const
const char * AsString(Resonance_t res)
resonance id -> string
const ProcessInfo & ProcInfo(void) const
const XSecAlgorithmI * fModel
const Target & Tgt(void) const
static Cache * Instance(void)
A simple cache branch storing the cached data in a TNtuple.
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
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)
d2XSecRESFast_dWQ2_E(const XSecAlgorithmI *m, const Interaction *i)
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