17 #include "Framework/Conventions/GBuild.h"
35 using namespace genie;
36 using namespace genie::controls;
37 using namespace genie::constants;
38 using namespace genie::utils;
62 <<
"Generating kinematics uniformly over the allowed phase space";
79 double hitNucPos = evrec->
HitNucleon()->
X4()->Vect().Mag();
96 LOG(
"DMELKinematics",
pWARN) <<
"No available phase space";
99 exception.
SetReason(
"No available phase space");
110 LOG(
"DMELKinematics",
pNOTICE) <<
"Setting max XSec";
112 LOG(
"DMELKinematics",
pNOTICE) <<
"Set max XSec to " << xsec_max;
124 unsigned int iter = 0;
130 <<
"Couldn't select a valid Q^2 after " << iter <<
" iterations";
133 exception.
SetReason(
"Couldn't select kinematics");
148 LOG(
"DMELKinematics",
pNOTICE) <<
"Attempting to sample";
149 gQ2 = Q2min + (Q2max-Q2min) * rnd->
RndKine().Rndm();
151 LOG(
"DMELKinematics",
pINFO) <<
"Trying: Q^2 = " <<
gQ2;
160 double t = xsec_max * rnd->
RndKine().Rndm();
164 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
166 <<
"xsec= " << xsec <<
", J= " << J <<
", Rnd= " << t;
168 accept = (t < J*xsec);
175 LOG(
"DMELKinematics",
pINFO) <<
"Selected: Q^2 = " <<
gQ2;
190 LOG(
"DMELKinematics",
pNOTICE) <<
"E = " << E <<
", M = "<< M;
205 LOG(
"DMELKinematics",
pNOTICE) <<
"Selected: W = "<< gW;
218 double totxsec = evrec->
XSec();
219 double wght = (vol/totxsec)*xsec;
220 LOG(
"DMELKinematics",
pNOTICE) <<
"Kinematics wght = "<< wght;
224 LOG(
"DMELKinematics",
pNOTICE) <<
"Current event wght = " << wght;
257 double hitNucPos = evrec->
HitNucleon()->
X4()->Vect().Mag();
275 double En0 = TMath::Sqrt(pxn*pxn + pyn*pyn + pzn*pzn + Mn*Mn);
276 double Eb = En0 - En;
285 LOG(
"DMELKinematics",
pWARN) <<
"No available phase space";
288 exception.
SetReason(
"No available phase space");
300 double xsec_max = this->
MaxXSec(evrec);
307 LOG(
"DMELKinematics",
pNOTICE) <<
"E = " << E <<
", M = "<< Mn;
320 unsigned int iter = 0;
326 <<
"Couldn't select a valid Q^2 after " << iter <<
" iterations";
329 exception.
SetReason(
"Couldn't select kinematics");
335 gQ2 = Q2min + (Q2max-Q2min) * rnd->
RndKine().Rndm();
352 LOG(
"DMELKinematics",
pNOTICE) <<
"W = "<< gW;
353 LOG(
"DMELKinematics",
pNOTICE) <<
"x = "<< gx;
354 LOG(
"DMELKinematics",
pNOTICE) <<
"y = "<< gy;
360 LOG(
"DMELKinematics",
pNOTICE) <<
"v = "<< gv;
363 double gvtilde = gv + Mn - Eb - TMath::Sqrt(Mn*Mn+pxn*pxn+pyn*pyn+pzn*pzn);
364 double gvtilde2 = gvtilde*gvtilde;
366 LOG(
"DMELKinematics",
pNOTICE) <<
"v~ = "<< gvtilde;
369 double gQ2tilde = gQ2 - gv2 + gvtilde2;
371 LOG(
"DMELKinematics",
pNOTICE) <<
"Q~^2 = "<< gQ2tilde;
383 double t = xsec_max * rnd->
RndKine().Rndm();
384 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
386 <<
"xsec= " << xsec <<
", Rnd= " << t;
484 double max_xsec = 0.0;
488 LOG(
"DMELKinematics",
pDEBUG) <<
"Range of Q^2: " << rQ2.
min <<
" to " << rQ2.
max;
489 if(rQ2.
min <=0 || rQ2.
max <= rQ2.
min)
return 0.;
497 double dlogQ2 = (logQ2max - logQ2min) /(N-1);
498 double xseclast = -1;
501 for(
int i=0; i<N; i++) {
502 double Q2 = TMath::Exp(logQ2min + i * dlogQ2);
505 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
506 LOG(
"DMELKinematics",
pDEBUG) <<
"xsec(Q2= " << Q2 <<
") = " << xsec;
508 max_xsec = TMath::Max(xsec, max_xsec);
509 increasing = xsec-xseclast>=0;
517 for(
int ib=0; ib<Nb; ib++) {
518 Q2 = TMath::Exp(TMath::Log(Q2) - dlogQ2);
519 if(Q2 < rQ2.
min)
continue;
522 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
523 LOG(
"DMELKinematics",
pDEBUG) <<
"xsec(Q2= " << Q2 <<
") = " << xsec;
525 max_xsec = TMath::Max(xsec, max_xsec);
535 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
537 SLOG(
"DMELKinematics",
pDEBUG) <<
"Max xsec in phase space = " << max_xsec;
~DMELKinematicsGenerator()
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)
void ProcessEventRecord(GHepRecord *event_rec) const
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
ComputeMaxXSec -> ComputeMaxXSec * fSafetyFactor.
int RecoilNucleonPdg(void) const
recoil nucleon pdg
A simple [min,max] interval for doubles.
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
double HitNucMass(void) const
int CharmHadronPdg(void) const
bool IsStrangeEvent(void) const
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec
Defines the EventGeneratorI interface.
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
void Configure(const Registry &config)
virtual double Weight(void) const
void SetHitNucPosition(double r)
Abstract class. Provides a data caching mechanism for for concrete implementations of the EventRecord...
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
Contains minimal information for tagging exclusive processes.
double PhaseSpaceVolume(const Interaction *const i, KinePhaseSpace_t ps)
bool IsCharmEvent(void) const
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...
int StrangeHadronPdg(void) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
virtual void Configure(const Registry &config)
static const double kASmallNum
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
double ComputeMaxXSec(const Interaction *in) const
void Setx(double x, bool selected=false)
static RunningThreadInfo * Instance(void)
void SetW(double W, bool selected=false)
void SwitchOnFastForward(void)
TLorentzVector * HitNucP4Ptr(void) const
static PDGLibrary * Instance(void)
void SetReason(string reason)
virtual TBits * EventFlags(void) const
const TLorentzVector * X4(void) const
A registry. Provides the container for algorithm configuration parameters.
void Sety(double y, bool selected=false)
DMELKinematicsGenerator()
const UInt_t kIAssumeFreeNucleon
const XclsTag & ExclTag(void) const
virtual GHepParticle * HitNucleon(void) const
Target * TgtPtr(void) const
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
virtual double XSec(void) const
InitialState * InitStatePtr(void) const
const InitialState & InitState(void) const
void SpectralFuncExperimentalCode(GHepRecord *event_rec) const
TParticlePDG * Find(int pdgc, bool must_exist=true)
void ClearRunningValues(void)
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
static const unsigned int kRjMaxIterations
const EventGeneratorI * RunningThread(void)
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
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...
const UInt_t kISkipProcessChk
if set, skip process validity checks
virtual void SetDiffXSec(double xsec, KinePhaseSpace_t ps)
Initial State information.