49 using namespace genie;
76 <<
"Generating kinematics uniformly over the allowed phase space";
93 double hitNucPos = evrec->
HitNucleon()->
X4()->Vect().Mag();
110 LOG(
"QELKinematics",
pWARN) <<
"No available phase space";
113 exception.
SetReason(
"No available phase space");
136 unsigned int iter = 0;
142 <<
"Couldn't select a valid Q^2 after " << iter <<
" iterations";
145 exception.
SetReason(
"Couldn't select kinematics");
160 gQ2 = Q2min + (Q2max-Q2min) * rnd->
RndKine().Rndm();
162 LOG(
"QELKinematics",
pINFO) <<
"Trying: Q^2 = " <<
gQ2;
171 double t = xsec_max * rnd->
RndKine().Rndm();
175 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 177 <<
"xsec= " << xsec <<
", J= " << J <<
", Rnd= " <<
t;
179 accept = (t < J*
xsec);
186 LOG(
"QELKinematics",
pINFO) <<
"Selected: Q^2 = " <<
gQ2;
201 LOG(
"QELKinematics",
pNOTICE) <<
"E = " << E <<
", M = "<< M;
216 LOG(
"QELKinematics",
pNOTICE) <<
"Selected: W = "<< gW;
229 double totxsec = evrec->
XSec();
230 double wght = (vol/totxsec)*xsec;
231 LOG(
"QELKinematics",
pNOTICE) <<
"Kinematics wght = "<< wght;
235 LOG(
"QELKinematics",
pNOTICE) <<
"Current event wght = " << wght;
268 double hitNucPos = evrec->
HitNucleon()->
X4()->Vect().Mag();
286 double En0 =
TMath::Sqrt(pxn*pxn + pyn*pyn + pzn*pzn + Mn*Mn);
287 double Eb = En0 - En;
295 LOG(
"QELKinematics",
pWARN) <<
"No available phase space";
298 exception.
SetReason(
"No available phase space");
310 double xsec_max = this->
MaxXSec(evrec);
317 LOG(
"QELKinematics",
pNOTICE) <<
"E = " << E <<
", M = "<< Mn;
330 unsigned int iter = 0;
336 <<
"Couldn't select a valid Q^2 after " << iter <<
" iterations";
339 exception.
SetReason(
"Couldn't select kinematics");
345 gQ2 = Q2min + (Q2max-Q2min) * rnd->
RndKine().Rndm();
373 double gvtilde = gv + Mn - Eb -
TMath::Sqrt(Mn*Mn+pxn*pxn+pyn*pyn+pzn*pzn);
374 double gvtilde2 = gvtilde*gvtilde;
376 LOG(
"QELKinematics",
pNOTICE) <<
"v~ = "<< gvtilde;
379 double gQ2tilde = gQ2 - gv2 + gvtilde2;
381 LOG(
"QELKinematics",
pNOTICE) <<
"Q~^2 = "<< gQ2tilde;
393 double t = xsec_max * rnd->
RndKine().Rndm();
394 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 396 <<
"xsec= " << xsec <<
", Rnd= " <<
t;
494 double max_xsec = 0.0;
498 if(rQ2.
min <=0 || rQ2.
max <= rQ2.
min)
return 0.;
506 double dlogQ2 = (logQ2max - logQ2min) /(N-1);
507 double xseclast = -1;
510 for(
int i=0;
i<N;
i++) {
511 double Q2 = TMath::Exp(logQ2min +
i * dlogQ2);
514 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 515 LOG(
"QELKinematics",
pDEBUG) <<
"xsec(Q2= " << Q2 <<
") = " <<
xsec;
517 max_xsec = TMath::Max(xsec, max_xsec);
518 increasing = xsec-xseclast>=0;
526 for(
int ib=0; ib<
Nb; ib++) {
527 Q2 = TMath::Exp(TMath::Log(Q2) - dlogQ2);
528 if(Q2 < rQ2.
min)
continue;
531 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 532 LOG(
"QELKinematics",
pDEBUG) <<
"xsec(Q2= " << Q2 <<
") = " <<
xsec;
534 max_xsec = TMath::Max(xsec, max_xsec);
544 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 546 SLOG(
"QELKinematics",
pDEBUG) <<
"Max xsec in phase space = " << max_xsec;
virtual double MaxXSec(GHepRecord *evrec) const
void SpectralFuncExperimentalCode(GHepRecord *event_rec) 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.
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 Weight(void) const
::xsd::cxx::tree::exception< char > exception
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)
void ProcessEventRecord(GHepRecord *event_rec) const
static const double kASmallNum
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Misc GENIE control constants.
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
double ComputeMaxXSec(const Interaction *in) const
void Setx(double x, bool selected=false)
static RunningThreadInfo * Instance(void)
~QELKinematicsGenerator()
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
Var Sqrt(const Var &v)
Use to take sqrt of a var.
const TLorentzVector * X4(void) const
A registry. Provides the container for algorithm configuration parameters.
void Sety(double y, bool selected=false)
const UInt_t kIAssumeFreeNucleon
const XclsTag & ExclTag(void) const
void Configure(const Registry &config)
virtual GHepParticle * HitNucleon(void) const
Target * TgtPtr(void) const
TParticlePDG * Find(int pdgc)
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
assert(nhit_max >=nhit_nbins)
virtual double XSec(void) const
InitialState * InitStatePtr(void) const
const InitialState & InitState(void) const
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...
Root of GENIE utility namespaces.
const UInt_t kISkipProcessChk
if set, skip process validity checks
virtual void SetDiffXSec(double xsec, KinePhaseSpace_t ps)
Initial State information.