GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
QELEventGenerator.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  New QE event generator written by:
7  Andy Furmanski
8  Manchester
9 
10  Using a skeleton and existing QE event generator templates by:
11  Costas Andreopoulos <c.andreopoulos \at cern.ch>
12  University of Liverpool
13 */
14 //____________________________________________________________________________
15 
16 #include <TMath.h>
17 
20 #include "Framework/Conventions/GBuild.h"
38 
43 
44 using namespace genie;
45 using namespace genie::controls;
46 using namespace genie::constants;
47 using namespace genie::utils;
48 
49 //___________________________________________________________________________
51  KineGeneratorWithCache("genie::QELEventGenerator")
52 {
53 
54 }
55 //___________________________________________________________________________
57  KineGeneratorWithCache("genie::QELEventGenerator", config)
58 {
59 
60 }
61 //___________________________________________________________________________
63 {
64 
65 }
66 //___________________________________________________________________________
68 {
69  LOG("QELEvent", pDEBUG) << "Generating QE event kinematics...";
70 
71  // Get the random number generators
73 
74  // Access cross section algorithm for running thread
76  const EventGeneratorI * evg = rtinfo->RunningThread();
77  fXSecModel = evg->CrossSectionAlg();
78 
79  // Get the interaction and check we are working with a nuclear target
80  Interaction * interaction = evrec->Summary();
81  // Skip if not a nuclear target
82  if(interaction->InitState().Tgt().IsNucleus()) {
83  // Skip if no hit nucleon is set
84  if(! evrec->HitNucleon()) {
85  LOG("QELEvent", pFATAL) << "No hit nucleon was set";
86  gAbortingInErr = true;
87  exit(1);
88  }
89  } // is nuclear target
90 
91  // set the 'trust' bits
92  interaction->SetBit(kISkipProcessChk);
93  interaction->SetBit(kISkipKinematicChk);
94 
95  // Access the hit nucleon and target nucleus entries at the GHEP record
96  GHepParticle * nucleon = evrec->HitNucleon();
97  GHepParticle * nucleus = evrec->TargetNucleus();
98  bool have_nucleus = (nucleus != 0);
99 
100  // Store the hit nucleon radius before computing the maximum differential
101  // cross section (important when using the local Fermi gas model)
102  Target* tgt = interaction->InitState().TgtPtr();
103  double hitNucPos = nucleon->X4()->Vect().Mag();
104  tgt->SetHitNucPosition( hitNucPos );
105 
106  //-- For the subsequent kinematic selection with the rejection method:
107  // Calculate the max differential cross section or retrieve it from the
108  // cache. Throw an exception and quit the evg thread if a non-positive
109  // value is found.
110  // If the kinematics are generated uniformly over the allowed phase
111  // space the max xsec is irrelevant
112  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
113 
114  // For a composite nuclear target, check to make sure that the
115  // final nucleus has a recognized PDG code
116  if ( have_nucleus ) {
117  // compute A,Z for final state nucleus & get its PDG code
118  int nucleon_pdgc = nucleon->Pdg();
119  bool is_p = pdg::IsProton(nucleon_pdgc);
120  int Z = (is_p) ? nucleus->Z()-1 : nucleus->Z();
121  int A = nucleus->A() - 1;
122  TParticlePDG * fnucleus = 0;
123  int ipdgc = pdg::IonPdgCode(A, Z);
124  fnucleus = PDGLibrary::Instance()->Find(ipdgc);
125  if (!fnucleus) {
126  LOG("QELEvent", pFATAL)
127  << "No particle with [A = " << A << ", Z = " << Z
128  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
129  exit(1);
130  }
131  }
132 
133  // In the accept/reject loop, each iteration samples a new value of
134  // - the hit nucleon 3-momentum,
135  // - its binding energy (only actually used if fHitNucleonBindingMode == kUseNuclearModel)
136  // - the final lepton scattering angles in the neutrino-and-hit-nucleon COM frame
137  // (measured with respect to the velocity of the COM frame as seen in the lab frame)
138  unsigned int iter = 0;
139  bool accept = false;
140  while (1) {
141 
142  iter++;
143  LOG("QELEvent", pINFO) << "Attempt #: " << iter;
144  if(iter > kRjMaxIterations) {
145  LOG("QELEvent", pWARN)
146  << "Couldn't select a valid (pNi, Eb, cos_theta_0, phi_0) tuple after "
147  << iter << " iterations";
148  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
150  exception.SetReason("Couldn't select kinematics");
151  exception.SwitchOnFastForward();
152  throw exception;
153  }
154 
155  // If the target is a composite nucleus, then sample an initial nucleon
156  // 3-momentum and removal energy from the nuclear model.
157  if ( tgt->IsNucleus() ) {
158  fNuclModel->GenerateNucleon(*tgt, hitNucPos);
159  }
160  else {
161  // Otherwise, just set the nucleon to be at rest in the lab frame and
162  // unbound. Use the nuclear model to make these assignments. The call
163  // to BindHitNucleon() will apply them correctly below.
164  fNuclModel->SetMomentum3( TVector3(0., 0., 0.) );
166  }
167 
168  // Put the hit nucleon off-shell (if needed) so that we can get the correct
169  // value of cos_theta0_max
172 
173  double cos_theta0_max = std::min(1., CosTheta0Max(*interaction));
174 
175  // If the allowed range of cos(theta_0) is vanishing, skip doing the
176  // full differential cross section calculation (it will be zero)
177  if ( cos_theta0_max <= -1. ) continue;
178 
179  // Pick a direction
180  // NOTE: In the kPSQELEvGen phase space used by this generator,
181  // these angles are specified with respect to the velocity of the
182  // probe + hit nucleon COM frame as measured in the lab frame. That is,
183  // costheta = 1 means that the outgoing lepton's COM frame 3-momentum
184  // points parallel to the velocity of the COM frame.
185  double costheta = rnd->RndKine().Uniform(-1., cos_theta0_max); // cosine theta
186  double phi = rnd->RndKine().Uniform( 2.*kPi ); // phi: [0, 2pi]
187 
188  // Set the "bind_nucleon" flag to false in this call to ComputeFullQELPXSec
189  // since we've already done that above
190  LOG("QELEvent", pDEBUG) << "cth0 = " << costheta << ", phi0 = " << phi;
191  double xsec = genie::utils::ComputeFullQELPXSec(interaction, fNuclModel,
192  fXSecModel, costheta, phi, fEb, fHitNucleonBindingMode, fMinAngleEM, false);
193 
194  // select/reject event
195  this->AssertXSecLimits(interaction, xsec, xsec_max);
196 
197  double t = xsec_max * rnd->RndKine().Rndm();
198 
199 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
200  LOG("QELEvent", pDEBUG)
201  << "xsec= " << xsec << ", Rnd= " << t;
202 #endif
203  accept = (t < xsec);
204 
205  // If the generated kinematics are accepted, finish-up module's job
206  if(accept) {
207  double gQ2 = interaction->KinePtr()->Q2(false);
208  LOG("QELEvent", pINFO) << "*Selected* Q^2 = " << gQ2 << " GeV^2";
209 
210  // reset bits
211  interaction->ResetBit(kISkipProcessChk);
212  interaction->ResetBit(kISkipKinematicChk);
213  interaction->ResetBit(kIAssumeFreeNucleon);
214 
215  // get neutrino energy at struck nucleon rest frame and the
216  // struck nucleon mass (can be off the mass shell)
217  const InitialState & init_state = interaction->InitState();
218  double E = init_state.ProbeE(kRfHitNucRest);
219  double M = init_state.Tgt().HitNucP4().M();
220  LOG("QELKinematics", pNOTICE) << "E = " << E << ", M = "<< M;
221 
222  // The hadronic inv. mass is equal to the recoil nucleon on-shell mass.
223  // For QEL/Charm events it is set to be equal to the on-shell mass of
224  // the generated charm baryon (Lamda_c+, Sigma_c+ or Sigma_c++)
225  // Similarly for strange baryons
226  //
227  const XclsTag & xcls = interaction->ExclTag();
228  int rpdgc = 0;
229  if (xcls.IsCharmEvent()) {
230  rpdgc = xcls.CharmHadronPdg();
231  } else if (xcls.IsStrangeEvent()) {
232  rpdgc = xcls.StrangeHadronPdg();
233  } else {
234  rpdgc = interaction->RecoilNucleonPdg();
235  }
236  assert(rpdgc);
237  double gW = PDGLibrary::Instance()->Find(rpdgc)->Mass();
238  LOG("QELEvent", pNOTICE) << "Selected: W = "<< gW;
239 
240  // (W,Q2) -> (x,y)
241  double gx=0, gy=0;
242  kinematics::WQ2toXY(E,M,gW,gQ2,gx,gy);
243 
244  // lock selected kinematics & clear running values
245  interaction->KinePtr()->SetQ2(gQ2, true);
246  interaction->KinePtr()->SetW (gW, true);
247  interaction->KinePtr()->Setx (gx, true);
248  interaction->KinePtr()->Sety (gy, true);
249  interaction->KinePtr()->ClearRunningValues();
250 
251  // set the cross section for the selected kinematics
252  evrec->SetDiffXSec(xsec, kPSQELEvGen);
253 
254  TLorentzVector lepton(interaction->KinePtr()->FSLeptonP4());
255  TLorentzVector outNucleon(interaction->KinePtr()->HadSystP4());
256  TLorentzVector x4l(*(evrec->Probe())->X4());
257 
258  // Add the final-state lepton to the event record
259  evrec->AddParticle(interaction->FSPrimLeptonPdg(), kIStStableFinalState,
260  evrec->ProbePosition(), -1, -1, -1, interaction->KinePtr()->FSLeptonP4(), x4l);
261 
262  // Set its polarization
264 
265  // Add the final-state nucleon to the event record
267  evrec->AddParticle(interaction->RecoilNucleonPdg(), ist, evrec->HitNucleonPosition(),
268  -1, -1, -1, interaction->KinePtr()->HadSystP4(), x4l);
269 
270  // Store struck nucleon momentum and binding energy
271  TLorentzVector p4ptr = interaction->InitStatePtr()->TgtPtr()->HitNucP4();
272  LOG("QELEvent",pNOTICE) << "pn: " << p4ptr.X() << ", "
273  << p4ptr.Y() << ", " << p4ptr.Z() << ", " << p4ptr.E();
274  nucleon->SetMomentum(p4ptr);
275  nucleon->SetRemovalEnergy(fEb);
276 
277  // add a recoiled nucleus remnant
278  this->AddTargetNucleusRemnant(evrec);
279 
280  break; // done
281  } else { // accept throw
282  LOG("QELEvent", pDEBUG) << "Reject current throw...";
283  }
284 
285  } // iterations - while(1) loop
286  LOG("QELEvent", pINFO) << "Done generating QE event kinematics!";
287 }
288 //___________________________________________________________________________
290 {
291  // add the remnant nuclear target at the GHEP record
292 
293  LOG("QELEvent", pINFO) << "Adding final state nucleus";
294 
295  double Px = 0;
296  double Py = 0;
297  double Pz = 0;
298  double E = 0;
299 
300  GHepParticle * nucleus = evrec->TargetNucleus();
301  bool have_nucleus = nucleus != 0;
302  if (!have_nucleus) return;
303 
304  int A = nucleus->A();
305  int Z = nucleus->Z();
306 
307  int fd = nucleus->FirstDaughter();
308  int ld = nucleus->LastDaughter();
309 
310  for(int id = fd; id <= ld; id++) {
311 
312  // compute A,Z for final state nucleus & get its PDG code and its mass
313  GHepParticle * particle = evrec->Particle(id);
314  assert(particle);
315  int pdgc = particle->Pdg();
316  bool is_p = pdg::IsProton (pdgc);
317  bool is_n = pdg::IsNeutron(pdgc);
318 
319  if (is_p) Z--;
320  if (is_p || is_n) A--;
321 
322  Px += particle->Px();
323  Py += particle->Py();
324  Pz += particle->Pz();
325  E += particle->E();
326 
327  }//daughters
328 
329  TParticlePDG * remn = 0;
330  int ipdgc = pdg::IonPdgCode(A, Z);
331  remn = PDGLibrary::Instance()->Find(ipdgc);
332  if(!remn) {
333  LOG("HadronicVtx", pFATAL)
334  << "No particle with [A = " << A << ", Z = " << Z
335  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
336  assert(remn);
337  }
338 
339  double Mi = nucleus->Mass();
340  Px *= -1;
341  Py *= -1;
342  Pz *= -1;
343  E = Mi-E;
344 
345  // Add the nucleus to the event record
346  LOG("QELEvent", pINFO)
347  << "Adding nucleus [A = " << A << ", Z = " << Z
348  << ", pdgc = " << ipdgc << "]";
349 
350  int imom = evrec->TargetNucleusPosition();
351  evrec->AddParticle(
352  ipdgc,kIStStableFinalState, imom,-1,-1,-1, Px,Py,Pz,E, 0,0,0,0);
353 
354  LOG("QELEvent", pINFO) << "Done";
355  LOG("QELEvent", pINFO) << *evrec;
356 }
357 //___________________________________________________________________________
359 {
360  Algorithm::Configure(config);
361  this->LoadConfig();
362 }
363 //____________________________________________________________________________
364 void QELEventGenerator::Configure(string config)
365 {
366  Algorithm::Configure(config);
367  this->LoadConfig();
368 }
369 //____________________________________________________________________________
371 {
372  // Load sub-algorithms and config data to reduce the number of registry
373  // lookups
374  fNuclModel = 0;
375 
376  RgKey nuclkey = "NuclearModel";
377 
378  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
379  assert(fNuclModel);
380 
381  // Safety factor for the maximum differential cross section
382  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.6 ) ;
383 
384  // Minimum energy for which max xsec would be cached, forcing explicit
385  // calculation for lower eneries
386  GetParamDef( "Cache-MinEnergy", fEMin, 1.00 ) ;
387 
388  // Maximum allowed fractional cross section deviation from maxim cross
389  // section used in rejection method
390  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
391  assert(fMaxXSecDiffTolerance>=0);
392 
393  // Generate kinematics uniformly over allowed phase space and compute
394  // an event weight?
395  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
396 
397  GetParamDef( "SF-MinAngleEMscattering", fMinAngleEM, 0. ) ;
398 
399  // Decide how to handle the binding energy of the initial state struck
400  // nucleon
401  std::string binding_mode;
402  GetParamDef( "HitNucleonBindingMode", binding_mode, std::string("UseNuclearModel") );
403 
405 
406  GetParamDef( "MaxXSecNucleonThrows", fMaxXSecNucleonThrows, 800 );
407 }
408 //____________________________________________________________________________
410 {
411  // Computes the maximum differential cross section in the requested phase
412  // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
413  // method and the value is cached at a circular cache branch for retrieval
414  // during subsequent event generation.
415  // The computed max differential cross section does not need to be the exact
416  // maximum. The number used in the rejection method will be scaled up by a
417  // safety factor. But it needs to be fast.
418  LOG("QELEvent", pINFO) << "Computing maximum cross section to throw against";
419 
420  double xsec_max = -1;
421  double dummy_Eb = 0.;
422 
423  // Clone the input interaction so that we can modify it a bit
424  Interaction * interaction = new Interaction( *in );
425  interaction->SetBit( kISkipProcessChk );
426  interaction->SetBit( kISkipKinematicChk );
427 
428  // We'll select the max momentum and zero binding energy.
429  // That should give us the nucleon with the highest xsec
430  const Target& tgt = interaction->InitState().Tgt();
431  // Pick a really slow nucleon to start, but not one at rest,
432  // since Prob() for the Fermi gas family of models is zero
433  // for a vanishing nucleon momentum
434  double max_momentum = 0.010;
435  double search_step = 0.1;
436  const double STEP_STOP = 1e-6;
437  while ( search_step > STEP_STOP ) {
438  double pNi_next = max_momentum + search_step;
439 
440  // Set the nucleon we're using to be upstream at max energy and unbound
441  fNuclModel->SetMomentum3( TVector3(0., 0., -pNi_next) );
443 
444  // Set the nucleon total energy to be on-shell with a quick call to
445  // BindHitNucleon()
446  genie::utils::BindHitNucleon(*interaction, *fNuclModel, dummy_Eb, kOnShell);
447 
448  // TODO: document this, won't work for spectral functions
449  double dummy_w = -1.;
450  double prob = fNuclModel->Prob(pNi_next, dummy_w, tgt,
451  tgt.HitNucPosition());
452 
453  double costh0_max = genie::utils::CosTheta0Max( *interaction );
454 
455  if ( prob > 0. && costh0_max > -1. ) max_momentum = pNi_next;
456  else search_step /= 2.;
457  }
458 
459  {
460  // Set the nucleon we're using to be upstream at max energy and unbound
461  fNuclModel->SetMomentum3( TVector3(0., 0., -max_momentum) );
463 
464  // Set the nucleon total energy to be on-shell with a quick call to
465  // BindHitNucleon()
466  genie::utils::BindHitNucleon(*interaction, *fNuclModel, dummy_Eb, kOnShell);
467 
468  // Just a scoping block for now
469  // OK, we're going to scan the COM frame angles to get the point of max xsec
470  // We'll bin in solid angle, and find the maximum point
471  // Then we'll bin/scan again inside that point
472  // Rinse and repeat until the xsec stabilises to within some fraction of our safety factor
473  const double acceptable_fraction_of_safety_factor = 0.2;
474  const int max_n_layers = 100;
475  const int N_theta = 10;
476  const int N_phi = 10;
477  double phi_at_xsec_max = -1;
478  double costh_at_xsec_max = 0;
479  double this_nuc_xsec_max = -1;
480 
481  double costh_range_min = -1.;
482  double costh_range_max = genie::utils::CosTheta0Max( *interaction );
483  double phi_range_min = 0.;
484  double phi_range_max = 2*TMath::Pi();
485  for (int ilayer = 0 ; ilayer < max_n_layers ; ilayer++) {
486  double last_layer_xsec_max = this_nuc_xsec_max;
487  double costh_increment = (costh_range_max-costh_range_min) / N_theta;
488  double phi_increment = (phi_range_max-phi_range_min) / N_phi;
489  // Now scan through centre-of-mass angles coarsely
490  for (int itheta = 0; itheta < N_theta; itheta++){
491  double costh = costh_range_min + itheta * costh_increment;
492  for (int iphi = 0; iphi < N_phi; iphi++) { // Scan around phi
493  double phi = phi_range_min + iphi * phi_increment;
494  // We're after an upper limit on the cross section, so just
495  // put the nucleon on-shell and call it good. The last
496  // argument is false because we've already called
497  // BindHitNucleon() above
498  double xs = genie::utils::ComputeFullQELPXSec(interaction,
499  fNuclModel, fXSecModel, costh, phi, dummy_Eb, kOnShell, fMinAngleEM, false);
500 
501  if (xs > this_nuc_xsec_max){
502  phi_at_xsec_max = phi;
503  costh_at_xsec_max = costh;
504  this_nuc_xsec_max = xs;
505  }
506  //
507  } // Done with phi scan
508  }// Done with centre-of-mass angles coarsely
509 
510  // Calculate the range for the next layer
511  costh_range_min = costh_at_xsec_max - costh_increment;
512  costh_range_max = costh_at_xsec_max + costh_increment;
513  phi_range_min = phi_at_xsec_max - phi_increment;
514  phi_range_max = phi_at_xsec_max + phi_increment;
515 
516  double improvement_factor = this_nuc_xsec_max/last_layer_xsec_max;
517  if (ilayer && (improvement_factor-1) < acceptable_fraction_of_safety_factor * (fSafetyFactor-1)) {
518  break;
519  }
520  }
521  if (this_nuc_xsec_max > xsec_max) {
522  xsec_max = this_nuc_xsec_max;
523  LOG("QELEvent", pINFO) << "best estimate for xsec_max = " << xsec_max;
524  }
525 
526  delete interaction;
527  }
528  // Apply safety factor, since value retrieved from the cache might
529  // correspond to a slightly different energy
530  xsec_max *= fSafetyFactor;
531 
532 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
533  SLOG("QELEvent", pDEBUG) << interaction->AsString();
534  SLOG("QELEvent", pDEBUG) << "Max xsec in phase space = " << max_xsec;
535  SLOG("QELEvent", pDEBUG) << "Computed using alg = " << *fXSecModel;
536 #endif
537 
538  LOG("QELEvent", pINFO) << "Computed maximum cross section to throw against - value is " << xsec_max;
539  return xsec_max;
540 }
541 //____________________________________________________________________________
void ProcessEventRecord(GHepRecord *event_rec) const
int Z(void) const
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double E(void) const
Get energy.
Definition: GHepParticle.h:91
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
const NuclearModelI * fNuclModel
nuclear model
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
int FirstDaughter(void) const
Definition: GHepParticle.h:68
int RecoilNucleonPdg(void) const
recoil nucleon pdg
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
int CharmHadronPdg(void) const
Definition: XclsTag.h:52
#define pFATAL
Definition: Messenger.h:56
bool IsStrangeEvent(void) const
Definition: XclsTag.h:53
bool IsNucleus(void) const
Definition: Target.cxx:272
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec&gt;maxxsec
Defines the EventGeneratorI interface.
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
virtual int HitNucleonPosition(void) const
Definition: GHepRecord.cxx:418
double ComputeMaxXSec(const Interaction *in) const
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
void SetHitNucPosition(double r)
Definition: Target.cxx:210
void SetMomentum(const TLorentzVector &p4)
double Mass(void) const
Mass that corresponds to the PDG code.
void BindHitNucleon(Interaction &interaction, const NuclearModelI &nucl_model, double &Eb, DMELEvGen_BindingMode_t hitNucleonBindingMode)
Definition: DMELUtils.cxx:259
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)
Definition: QELUtils.cxx:93
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...
Definition: RandomGen.h:29
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
string AsString(void) const
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:284
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
int Pdg(void) const
Definition: GHepParticle.h:63
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:91
int LastDaughter(void) const
Definition: GHepParticle.h:69
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:65
const double e
int StrangeHadronPdg(void) const
Definition: XclsTag.h:55
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void AddTargetNucleusRemnant(GHepRecord *evrec) const
add a recoiled nucleus remnant
static constexpr double A
Definition: Units.h:74
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1132
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
virtual double Prob(double p, double w, const Target &) const =0
#define pINFO
Definition: Messenger.h:62
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
#define pWARN
Definition: Messenger.h:60
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
void SetRemovalEnergy(double Erm)
static RunningThreadInfo * Instance(void)
void SetRemovalEnergy(double E) const
Definition: NuclearModelI.h:90
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:117
void SetMomentum3(const TVector3 &mom) const
Definition: NuclearModelI.h:86
string RgKey
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
double CosTheta0Max(const genie::Interaction &interaction)
Definition: DMELUtils.cxx:217
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
double HitNucPosition(void) const
Definition: Target.h:89
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:313
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
Definition: QELUtils.cxx:194
Target * TgtPtr(void) const
Definition: InitialState.h:67
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
double gQ2
Definition: gtestDISSF.cxx:55
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
const InitialState & InitState(void) const
Definition: Interaction.h:69
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
int A(void) const
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
Definition: InitialState.h:66
virtual int ProbePosition(void) const
Definition: GHepRecord.cxx:352
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
const EventGeneratorI * RunningThread(void)
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
virtual bool GenerateNucleon(const Target &) const =0
bool gAbortingInErr
Definition: Messenger.cxx:34
void SetPrimaryLeptonPolarization(GHepRecord *ev)
void Configure(const Registry &config)
double ProbeE(RefFrame_t rf) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
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.
Definition: GHepParticle.h:39
Keep info on the event generation thread currently on charge. This is used so that event generation m...
virtual int TargetNucleusPosition(void) const
Definition: GHepRecord.cxx:370
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
virtual void SetDiffXSec(double xsec, KinePhaseSpace_t ps)
Definition: GHepRecord.h:133
enum genie::EGHepStatus GHepStatus_t
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
QELEvGen_BindingMode_t fHitNucleonBindingMode
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345