55 using namespace genie;
80 LOG(
"QELEvent",
pDEBUG) <<
"Generating QE event kinematics...";
96 LOG(
"QELEvent",
pFATAL) <<
"No hit nucleon was set";
109 bool have_nucleus = (nucleus != 0);
114 double hitNucPos = nucleon->
X4()->Vect().Mag();
127 if ( have_nucleus ) {
129 int nucleon_pdgc = nucleon->
Pdg();
131 int Z = (is_p) ? nucleus->
Z()-1 : nucleus->
Z();
132 int A = nucleus->
A() - 1;
133 TParticlePDG * fnucleus = 0;
138 <<
"No particle with [A = " << A <<
", Z = " << Z
139 <<
", pdgc = " << ipdgc <<
"] in PDGLibrary!";
149 unsigned int iter = 0;
154 LOG(
"QELEvent",
pINFO) <<
"Attempt #: " << iter;
157 <<
"Couldn't select a valid (pNi, Eb, cos_theta_0, phi_0) tuple after " 158 << iter <<
" iterations";
161 exception.
SetReason(
"Couldn't select kinematics");
188 if ( cos_theta0_max <= -1. )
continue;
196 double costheta = rnd->
RndKine().Uniform(-1., cos_theta0_max);
197 double phi = rnd->
RndKine().Uniform( 2.*
kPi );
201 LOG(
"QELEvent",
pDEBUG) <<
"cth0 = " << costheta <<
", phi0 = " << phi;
208 double t = xsec_max * rnd->
RndKine().Rndm();
210 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 212 <<
"xsec= " << xsec <<
", Rnd= " <<
t;
219 LOG(
"QELEvent",
pINFO) <<
"*Selected* Q^2 = " << gQ2 <<
" GeV^2";
231 LOG(
"QELKinematics",
pNOTICE) <<
"E = " << E <<
", M = "<< M;
249 LOG(
"QELEvent",
pNOTICE) <<
"Selected: W = "<< gW;
267 TLorentzVector x4l(*(evrec->
Probe())->X4());
278 LOG(
"QELEvent",
pNOTICE) <<
"pn: " << p4ptr.X() <<
", " 279 << p4ptr.Y() <<
", " << p4ptr.Z() <<
", " << p4ptr.E();
288 LOG(
"QELEvent",
pDEBUG) <<
"Reject current throw...";
292 LOG(
"QELEvent",
pINFO) <<
"Done generating QE event kinematics!";
299 LOG(
"QELEvent",
pINFO) <<
"Adding final state nucleus";
307 bool have_nucleus = nucleus != 0;
308 if (!have_nucleus)
return;
310 int A = nucleus->
A();
311 int Z = nucleus->
Z();
316 for(
int id = fd;
id <= ld;
id++) {
321 int pdgc = particle->
Pdg();
326 if (is_p || is_n) A--;
328 Px += particle->
Px();
329 Py += particle->
Py();
330 Pz += particle->
Pz();
335 TParticlePDG * remn = 0;
340 <<
"No particle with [A = " << A <<
", Z = " << Z
341 <<
", pdgc = " << ipdgc <<
"] in PDGLibrary!";
345 double Mi = nucleus->
Mass();
353 <<
"Adding nucleus [A = " << A <<
", Z = " << Z
354 <<
", pdgc = " << ipdgc <<
"]";
358 ipdgc,
kIStStableFinalState, imom,-1,-1,-1, Px,Py,Pz,E, 0,0,0,0);
382 RgKey nuclkey =
"NuclearModel";
424 LOG(
"QELEvent",
pINFO) <<
"Computing maximum cross section to throw against";
426 double xsec_max = -1;
427 double dummy_Eb = 0.;
440 double max_momentum = 0.010;
441 double search_step = 0.1;
442 const double STEP_STOP = 1
e-6;
443 while ( search_step > STEP_STOP ) {
444 double pNi_next = max_momentum + search_step;
455 double dummy_w = -1.;
461 if ( prob > 0. && costh0_max > -1. ) max_momentum = pNi_next;
462 else search_step /= 2.;
479 const double acceptable_fraction_of_safety_factor = 0.2;
480 const int max_n_layers = 100;
481 const int N_theta = 10;
482 const int N_phi = 10;
483 double phi_at_xsec_max = -1;
484 double costh_at_xsec_max = 0;
485 double this_nuc_xsec_max = -1;
487 double costh_range_min = -1.;
489 double phi_range_min = 0.;
490 double phi_range_max = 2*TMath::Pi();
491 for (
int ilayer = 0 ; ilayer < max_n_layers ; ilayer++) {
492 double last_layer_xsec_max = this_nuc_xsec_max;
493 double costh_increment = (costh_range_max-costh_range_min) / N_theta;
494 double phi_increment = (phi_range_max-phi_range_min) / N_phi;
496 for (
int itheta = 0; itheta < N_theta; itheta++){
497 double costh = costh_range_min + itheta * costh_increment;
498 for (
int iphi = 0; iphi < N_phi; iphi++) {
499 double phi = phi_range_min + iphi * phi_increment;
507 if (xs > this_nuc_xsec_max){
508 phi_at_xsec_max = phi;
509 costh_at_xsec_max = costh;
510 this_nuc_xsec_max = xs;
517 costh_range_min = costh_at_xsec_max - costh_increment;
518 costh_range_max = costh_at_xsec_max + costh_increment;
519 phi_range_min = phi_at_xsec_max - phi_increment;
520 phi_range_max = phi_at_xsec_max + phi_increment;
522 double improvement_factor = this_nuc_xsec_max/last_layer_xsec_max;
523 if (ilayer && (improvement_factor-1) < acceptable_fraction_of_safety_factor * (
fSafetyFactor-1)) {
527 if (this_nuc_xsec_max > xsec_max) {
528 xsec_max = this_nuc_xsec_max;
529 LOG(
"QELEvent",
pINFO) <<
"best estimate for xsec_max = " << xsec_max;
538 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 540 SLOG(
"QELEvent",
pDEBUG) <<
"Max xsec in phase space = " << max_xsec;
544 LOG(
"QELEvent",
pINFO) <<
"Computed maximum cross section to throw against - value is " << xsec_max;
void ProcessEventRecord(GHepRecord *event_rec) const
virtual double MaxXSec(GHepRecord *evrec) const
virtual GHepParticle * Particle(int position) const
bool fGenerateUniformly
uniform over allowed phase space + event weight?
THE MAIN GENIE PROJECT NAMESPACE
double E(void) const
Get energy.
virtual Interaction * Summary(void) const
const NuclearModelI * fNuclModel
nuclear model
static RandomGen * Instance()
Access instance.
void SetQ2(double Q2, bool selected=false)
Kinematics * KinePtr(void) const
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
int FirstDaughter(void) const
int RecoilNucleonPdg(void) const
recoil nucleon pdg
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
enum genie::EGHepStatus GHepStatus_t
int CharmHadronPdg(void) const
bool IsStrangeEvent(void) const
bool IsNucleus(void) const
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec
Defines the EventGeneratorI interface.
virtual int HitNucleonPosition(void) const
double ComputeMaxXSec(const Interaction *in) const
const TLorentzVector & HadSystP4(void) const
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
::xsd::cxx::tree::exception< char > exception
void SetHitNucPosition(double r)
void SetMomentum(const TLorentzVector &p4)
double Mass(void) const
Mass that corresponds to the PDG code.
double ComputeFullQELPXSec(Interaction *interaction, const NuclearModelI *nucl_model, const XSecAlgorithmI *xsec_model, double cos_theta_0, double phi_0, double &Eb, QELEvGen_BindingMode_t hitNucleonBindingMode, double min_angle_EM=0., bool bind_nucleon=true)
void BindHitNucleon(Interaction &interaction, const NuclearModelI &nucl_model, double &Eb, QELEvGen_BindingMode_t hitNucleonBindingMode)
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...
double Pz(void) const
Get Pz.
string AsString(void) const
const XSecAlgorithmI * fXSecModel
Contains minimal information for tagging exclusive processes.
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
double Px(void) const
Get Px.
virtual GHepParticle * Probe(void) const
bool IsCharmEvent(void) const
Summary information for an interaction.
const TLorentzVector & HitNucP4(void) const
int LastDaughter(void) const
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
const TLorentzVector & FSLeptonP4(void) const
int StrangeHadronPdg(void) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
void AddTargetNucleusRemnant(GHepRecord *evrec) const
add a recoiled nucleus remnant
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
virtual void Configure(const Registry &config)
virtual GHepParticle * TargetNucleus(void) const
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
void Setx(double x, bool selected=false)
void SetRemovalEnergy(double Erm)
static RunningThreadInfo * Instance(void)
void SetRemovalEnergy(double E) const
void SetW(double W, bool selected=false)
void SwitchOnFastForward(void)
static PDGLibrary * Instance(void)
void SetReason(string reason)
virtual TBits * EventFlags(void) const
void SetMomentum3(const TVector3 &mom) const
const TLorentzVector * X4(void) const
A registry. Provides the container for algorithm configuration parameters.
int fMaxXSecNucleonThrows
void Sety(double y, bool selected=false)
double CosTheta0Max(const genie::Interaction &interaction)
const UInt_t kIAssumeFreeNucleon
const XclsTag & ExclTag(void) const
double HitNucPosition(void) const
virtual GHepParticle * HitNucleon(void) const
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
Target * TgtPtr(void) const
TParticlePDG * Find(int pdgc)
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
int IonPdgCode(int A, int Z)
assert(nhit_max >=nhit_nbins)
InitialState * InitStatePtr(void) const
virtual bool GenerateNucleon(const Target &) const =0
virtual void AddParticle(const GHepParticle &p)
const InitialState & InitState(void) const
double Q2(bool selected=false) const
void ClearRunningValues(void)
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
virtual int ProbePosition(void) const
T min(const caf::Proxy< T > &a, T b)
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...
void Configure(const Registry &config)
virtual double Prob(double p, double w, const Target &) const =0
double ProbeE(RefFrame_t rf) const
GENIE's GHEP MC event record.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Keep info on the event generation thread currently on charge. This is used so that event generation m...
virtual int TargetNucleusPosition(void) const
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.
QELEvGen_BindingMode_t fHitNucleonBindingMode
double Py(void) const
Get Py.
const Algorithm * SubAlg(const RgKey ®istry_key) const