52 using namespace genie;
78 <<
"Generating kinematics uniformly over the allowed phase space";
104 LOG(
"RESKinematics",
pWARN) <<
"No available phase space";
107 exception.
SetReason(
"No available phase space");
112 const InitialState & init_state = interaction -> InitState();
126 double dW = W.
max - W.
min;
129 unsigned int iter = 0;
135 <<
"*** Could not select a valid (W,Q^2) pair after " 136 << iter <<
" iterations";
139 exception.
SetReason(
"Couldn't select kinematics");
154 if(Q2.
max<=0. || Q2.
min>=Q2.
max)
continue;
168 LOG(
"RESKinematics",
pINFO) <<
"Initializing the sampling envelope";
170 LOG(
"RESKinematics",
pFATAL) <<
"Null sampling envelope!";
184 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 186 <<
"Q^2: [" << Q2min <<
", " << Q2max <<
"] => " 187 <<
"QD^2: [" << QD2min <<
", " << QD2max <<
"]";
196 gR = (E>mR) ? 0.220 : 0.400;
198 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 200 <<
"(m,g) = (" << mR <<
", " << gR
201 <<
"), max(xsec,W) = (" << xsec_max <<
", " << W.
max <<
")";
217 LOG(
"RESKinematics",
pINFO) <<
"Trying: W = " << gW <<
", Q2 = " <<
gQ2;
231 double t = max * rnd->
RndKine().Rndm();
236 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 238 <<
"xsec= " << xsec <<
", J= " << J <<
", Rnd= " <<
t;
240 accept = (t < J*
xsec);
249 <<
"Selected: W = " << gW <<
", Q2 = " <<
gQ2;
268 double totxsec = evrec->
XSec();
269 double wght = (vol/totxsec)*xsec;
270 LOG(
"RESKinematics",
pNOTICE) <<
"Kinematics wght = "<< wght;
274 LOG(
"RESKinematics",
pNOTICE) <<
"Current event wght = " << wght;
329 gROOT->GetListOfFunctions()->Remove(
fEnvelope);
343 double max_xsec = 0.;
345 const InitialState & init_state = interaction -> InitState();
350 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 351 LOG(
"RESKinematics",
pDEBUG) <<
"Scanning phase space for E= " <<
E;
366 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 368 <<
"Will search for max{xsec} along W(=MRes) = " <<
md;
375 if( rQ2.
max < Q2Thres || rQ2.
min <=0 )
return 0.;
381 double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
383 double xseclast = -1;
384 bool increasing =
true;
386 for(
int iq2=0; iq2<NQ2; iq2++) {
387 double Q2 = TMath::Exp(logQ2min + iq2 * dlogQ2);
390 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 392 <<
"xsec(W= " << md <<
", Q2= " << Q2 <<
") = " <<
xsec;
394 max_xsec = TMath::Max(xsec, max_xsec);
395 increasing = xsec-xseclast>=0;
403 for(
int iq2b=0; iq2b<NQ2b; iq2b++) {
404 Q2 = TMath::Exp(TMath::Log(Q2) - dlogQ2);
405 if(Q2 < rQ2.
min)
continue;
408 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 410 <<
"xsec(W= " << md <<
", Q2= " << Q2 <<
") = " <<
xsec;
412 max_xsec = TMath::Max(xsec, max_xsec);
429 Wmax = TMath::Min(Wmax,
fWcut);
431 Wmin = TMath::Max(Wmin, md-.3);
432 Wmax = TMath::Min(Wmax, md+.3);
434 if(Wmax-Wmin<0.05) { NW=1; Wmin=
Wmax; }
436 double dW = (NW>1) ? (Wmax-Wmin)/(NW-1) : 0.;
438 for(
int iw=0; iw<NW; iw++) {
439 double W = Wmin + iw*dW;
446 if( rQ2.
max < Q2Thres || rQ2.
min <=0 )
continue;
447 if( rQ2.
max-rQ2.
min<0.02 ) {NQ2=5; NQ2b=3;}
451 double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
452 double xseclast = -1;
453 bool increasing =
true;
455 for(
int iq2=0; iq2<NQ2; iq2++) {
456 double Q2 = TMath::Exp(logQ2min + iq2 * dlogQ2);
460 <<
"xsec(W= " << W <<
", Q2= " << Q2 <<
") = " <<
xsec;
461 max_xsec = TMath::Max(xsec, max_xsec);
462 increasing = xsec-xseclast>=0;
470 for(
int iq2b=0; iq2b<NQ2b; iq2b++) {
471 Q2 = TMath::Exp(TMath::Log(Q2) - dlogQ2);
472 if(Q2 < rQ2.
min)
continue;
476 <<
"xsec(W= " << W <<
", Q2= " << Q2 <<
") = " <<
xsec;
477 max_xsec = TMath::Max(xsec, max_xsec);
490 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 492 LOG(
"RESKinematics",
pDEBUG) <<
"Max xsec in phase space = " << max_xsec;
double fWcut
Wcut parameter in DIS/RES join scheme.
virtual double MaxXSec(GHepRecord *evrec) const
const KPhaseSpace & PhaseSpace(void) const
virtual void SetWeight(double wght)
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double J(double q0, double q3, double Enu, double ml)
THE MAIN GENIE PROJECT NAMESPACE
double Q2(const Interaction *const i)
virtual Interaction * Summary(void) const
static RandomGen * Instance()
Access instance.
void SetQ2(double Q2, bool selected=false)
Kinematics * KinePtr(void) const
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
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.
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
bool KnownResonance(void) const
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec
Defines the EventGeneratorI interface.
virtual double Weight(void) const
::xsd::cxx::tree::exception< char > exception
double Mass(Resonance_t res)
resonance mass (GeV)
static const double kMinQ2Limit
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
Abstract class. Provides a data caching mechanism for for concrete implementations of the EventRecord...
enum genie::EResonance Resonance_t
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
string AsString(void) const
const XSecAlgorithmI * fXSecModel
static const double kMinQ2Limit
double PhaseSpaceVolume(const Interaction *const i, KinePhaseSpace_t ps)
Summary information for an interaction.
const TLorentzVector & HitNucP4(void) const
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
double QD2toQ2(double QD2)
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
virtual void Configure(const Registry &config)
static const double kASmallNum
~RESKinematicsGenerator()
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
TF2 * fEnvelope
2-D envelope used for importance sampling
Resonance_t Resonance(void) const
double Q2toQD2(double Q2)
Misc GENIE control constants.
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
void Setx(double x, bool selected=false)
static RunningThreadInfo * Instance(void)
virtual const AlgId & Id(void) const
Get algorithm ID.
void SetW(double W, bool selected=false)
void SwitchOnFastForward(void)
void SetReason(string reason)
virtual TBits * EventFlags(void) const
A registry. Provides the container for algorithm configuration parameters.
void Sety(double y, bool selected=false)
const XclsTag & ExclTag(void) const
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
void ProcessEventRecord(GHepRecord *event_rec) const
assert(nhit_max >=nhit_nbins)
virtual double XSec(void) const
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
const ProcessInfo & ProcInfo(void) const
void ClearRunningValues(void)
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
static const unsigned int kRjMaxIterations
void Configure(const Registry &config)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
const EventGeneratorI * RunningThread(void)
double ProbeE(RefFrame_t rf) const
GENIE's GHEP MC event record.
Keep info on the event generation thread currently on charge. This is used so that event generation m...
double ComputeMaxXSec(const Interaction *interaction) const
Root of GENIE utility namespaces.
Range1D_t WLim(void) const
W limits.
const UInt_t kISkipProcessChk
if set, skip process validity checks
virtual void SetDiffXSec(double xsec, KinePhaseSpace_t ps)
Initial State information.
double RESImportanceSamplingEnvelope(double *x, double *par)