GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DMELEventGenerator.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 
7  Costas Andreopoulos <c.andreopoulos \at cern.ch>
8  University of Liverpool
9  Joshua Berger <josh.berger \at pitt.edu>
10  Univeristy of Pittsburgh
11 
12  For the class documentation see the corresponding header file.
13 
14 */
15 //____________________________________________________________________________
16 
17 #include <TMath.h>
18 
21 #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::DMELEventGenerator")
52 {
53 
54 }
55 //___________________________________________________________________________
57  KineGeneratorWithCache("genie::DMELEventGenerator", config)
58 {
59 
60 }
61 //___________________________________________________________________________
63 {
64 
65 }
66 //___________________________________________________________________________
68 {
69  LOG("DMELEvent", 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("DMELEvent", 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("DMELEvent", 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("DMELEvent", pINFO) << "Attempt #: " << iter;
144  if(iter > kRjMaxIterations) {
145  LOG("DMELEvent", 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 kPSDMELEvGen 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 ComputeFullDMELPXSec
189  // since we've already done that above
190  LOG("DMELEvent", pDEBUG) << "cth0 = " << costheta << ", phi0 = " << phi;
191  double xsec = genie::utils::ComputeFullDMELPXSec(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("DMELEvent", 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("DMELEvent", 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("DMELKinematics", pNOTICE) << "E = " << E << ", M = "<< M;
221 
222  // The hadronic inv. mass is equal to the recoil nucleon on-shell mass.
223  // For DMEL/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("DMELEvent", 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, kPSDMELEvGen);
253 
254  TLorentzVector lepton(interaction->KinePtr()->FSLeptonP4());
255  TLorentzVector outNucleon(interaction->KinePtr()->HadSystP4());
256  TLorentzVector x4l(*(evrec->Probe())->X4());
257 
258  evrec->AddParticle(interaction->FSPrimLeptonPdg(), kIStStableFinalState,
259  evrec->ProbePosition(), -1, -1, -1, interaction->KinePtr()->FSLeptonP4(), x4l);
260 
262  evrec->AddParticle(interaction->RecoilNucleonPdg(), ist, evrec->HitNucleonPosition(),
263  -1, -1, -1, interaction->KinePtr()->HadSystP4(), x4l);
264 
265  // Store struck nucleon momentum and binding energy
266  TLorentzVector p4ptr = interaction->InitStatePtr()->TgtPtr()->HitNucP4();
267  LOG("DMELEvent",pNOTICE) << "pn: " << p4ptr.X() << ", "
268  << p4ptr.Y() << ", " << p4ptr.Z() << ", " << p4ptr.E();
269  nucleon->SetMomentum(p4ptr);
270  nucleon->SetRemovalEnergy(fEb);
271 
272  // add a recoiled nucleus remnant
273  this->AddTargetNucleusRemnant(evrec);
274 
275  break; // done
276  } else { // accept throw
277  LOG("DMELEvent", pDEBUG) << "Reject current throw...";
278  }
279 
280  } // iterations - while(1) loop
281  LOG("DMELEvent", pINFO) << "Done generating QE event kinematics!";
282 }
283 //___________________________________________________________________________
285 {
286  // add the remnant nuclear target at the GHEP record
287 
288  LOG("DMELEvent", pINFO) << "Adding final state nucleus";
289 
290  double Px = 0;
291  double Py = 0;
292  double Pz = 0;
293  double E = 0;
294 
295  GHepParticle * nucleus = evrec->TargetNucleus();
296  bool have_nucleus = nucleus != 0;
297  if (!have_nucleus) return;
298 
299  int A = nucleus->A();
300  int Z = nucleus->Z();
301 
302  int fd = nucleus->FirstDaughter();
303  int ld = nucleus->LastDaughter();
304 
305  for(int id = fd; id <= ld; id++) {
306 
307  // compute A,Z for final state nucleus & get its PDG code and its mass
308  GHepParticle * particle = evrec->Particle(id);
309  assert(particle);
310  int pdgc = particle->Pdg();
311  bool is_p = pdg::IsProton (pdgc);
312  bool is_n = pdg::IsNeutron(pdgc);
313 
314  if (is_p) Z--;
315  if (is_p || is_n) A--;
316 
317  Px += particle->Px();
318  Py += particle->Py();
319  Pz += particle->Pz();
320  E += particle->E();
321 
322  }//daughters
323 
324  TParticlePDG * remn = 0;
325  int ipdgc = pdg::IonPdgCode(A, Z);
326  remn = PDGLibrary::Instance()->Find(ipdgc);
327  if(!remn) {
328  LOG("HadronicVtx", pFATAL)
329  << "No particle with [A = " << A << ", Z = " << Z
330  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
331  assert(remn);
332  }
333 
334  double Mi = nucleus->Mass();
335  Px *= -1;
336  Py *= -1;
337  Pz *= -1;
338  E = Mi-E;
339 
340  // Add the nucleus to the event record
341  LOG("DMELEvent", pINFO)
342  << "Adding nucleus [A = " << A << ", Z = " << Z
343  << ", pdgc = " << ipdgc << "]";
344 
345  int imom = evrec->TargetNucleusPosition();
346  evrec->AddParticle(
347  ipdgc,kIStStableFinalState, imom,-1,-1,-1, Px,Py,Pz,E, 0,0,0,0);
348 
349  LOG("DMELEvent", pINFO) << "Done";
350  LOG("DMELEvent", pINFO) << *evrec;
351 }
352 //___________________________________________________________________________
354 {
355  Algorithm::Configure(config);
356  this->LoadConfig();
357 }
358 //____________________________________________________________________________
359 void DMELEventGenerator::Configure(string config)
360 {
361  Algorithm::Configure(config);
362  this->LoadConfig();
363 }
364 //____________________________________________________________________________
366 {
367  // Load sub-algorithms and config data to reduce the number of registry
368  // lookups
369  fNuclModel = 0;
370 
371  RgKey nuclkey = "NuclearModel";
372 
373  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
374  assert(fNuclModel);
375 
376  // Safety factor for the maximum differential cross section
377  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.6 ) ;
378 
379  // Minimum energy for which max xsec would be cached, forcing explicit
380  // calculation for lower eneries
381  GetParamDef( "Cache-MinEnergy", fEMin, 1.00 ) ;
382 
383  // Maximum allowed fractional cross section deviation from maxim cross
384  // section used in rejection method
385  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
386  assert(fMaxXSecDiffTolerance>=0);
387 
388  // Generate kinematics uniformly over allowed phase space and compute
389  // an event weight?
390  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
391 
392  GetParamDef( "SF-MinAngleEMscattering", fMinAngleEM, 0. ) ;
393 
394  // Decide how to handle the binding energy of the initial state struck
395  // nucleon
396  std::string binding_mode;
397  GetParamDef( "HitNucleonBindingMode", binding_mode, std::string("UseNuclearModel") );
398 
400 
401  GetParamDef( "MaxXSecNucleonThrows", fMaxXSecNucleonThrows, 800 );
402 }
403 //____________________________________________________________________________
405 {
406  // Computes the maximum differential cross section in the requested phase
407  // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
408  // method and the value is cached at a circular cache branch for retrieval
409  // during subsequent event generation.
410  // The computed max differential cross section does not need to be the exact
411  // maximum. The number used in the rejection method will be scaled up by a
412  // safety factor. But it needs to be fast.
413  LOG("DMELEvent", pINFO) << "Computing maximum cross section to throw against";
414 
415  double xsec_max = -1;
416  double dummy_Eb = 0.;
417 
418  // Clone the input interaction so that we can modify it a bit
419  Interaction * interaction = new Interaction( *in );
420  interaction->SetBit( kISkipProcessChk );
421  interaction->SetBit( kISkipKinematicChk );
422 
423  // We'll select the max momentum and zero binding energy.
424  // That should give us the nucleon with the highest xsec
425  const Target& tgt = interaction->InitState().Tgt();
426  // Pick a really slow nucleon to start, but not one at rest,
427  // since Prob() for the Fermi gas family of models is zero
428  // for a vanishing nucleon momentum
429  double max_momentum = 0.010;
430  double search_step = 0.1;
431  const double STEP_STOP = 1e-6;
432  while ( search_step > STEP_STOP ) {
433  double pNi_next = max_momentum + search_step;
434 
435  // Set the nucleon we're using to be upstream at max energy and unbound
436  fNuclModel->SetMomentum3( TVector3(0., 0., -pNi_next) );
438 
439  // Set the nucleon total energy to be on-shell with a quick call to
440  // BindHitNucleon()
441  genie::utils::BindHitNucleon(*interaction, *fNuclModel, dummy_Eb, kOnShell);
442 
443  // TODO: document this, won't work for spectral functions
444  double dummy_w = -1.;
445  double prob = fNuclModel->Prob(pNi_next, dummy_w, tgt,
446  tgt.HitNucPosition());
447 
448  double costh0_max = genie::utils::CosTheta0Max( *interaction );
449 
450  if ( prob > 0. && costh0_max > -1. ) max_momentum = pNi_next;
451  else search_step /= 2.;
452  }
453 
454  {
455  // Set the nucleon we're using to be upstream at max energy and unbound
456  fNuclModel->SetMomentum3( TVector3(0., 0., -max_momentum) );
458 
459  // Set the nucleon total energy to be on-shell with a quick call to
460  // BindHitNucleon()
461  genie::utils::BindHitNucleon(*interaction, *fNuclModel, dummy_Eb, kOnShell);
462 
463  // Just a scoping block for now
464  // OK, we're going to scan the COM frame angles to get the point of max xsec
465  // We'll bin in solid angle, and find the maximum point
466  // Then we'll bin/scan again inside that point
467  // Rinse and repeat until the xsec stabilises to within some fraction of our safety factor
468  const double acceptable_fraction_of_safety_factor = 0.2;
469  const int max_n_layers = 100;
470  const int N_theta = 10;
471  const int N_phi = 10;
472  double phi_at_xsec_max = -1;
473  double costh_at_xsec_max = 0;
474  double this_nuc_xsec_max = -1;
475 
476  double costh_range_min = -1.;
477  double costh_range_max = genie::utils::CosTheta0Max( *interaction );
478  double phi_range_min = 0.;
479  double phi_range_max = 2*TMath::Pi();
480  for (int ilayer = 0 ; ilayer < max_n_layers ; ilayer++) {
481  double last_layer_xsec_max = this_nuc_xsec_max;
482  double costh_increment = (costh_range_max-costh_range_min) / N_theta;
483  double phi_increment = (phi_range_max-phi_range_min) / N_phi;
484  // Now scan through centre-of-mass angles coarsely
485  for (int itheta = 0; itheta < N_theta; itheta++){
486  double costh = costh_range_min + itheta * costh_increment;
487  for (int iphi = 0; iphi < N_phi; iphi++) { // Scan around phi
488  double phi = phi_range_min + iphi * phi_increment;
489  // We're after an upper limit on the cross section, so just
490  // put the nucleon on-shell and call it good. The last
491  // argument is false because we've already called
492  // BindHitNucleon() above
493  double xs = genie::utils::ComputeFullDMELPXSec(interaction,
494  fNuclModel, fXSecModel, costh, phi, dummy_Eb, kOnShell, fMinAngleEM, false);
495 
496  if (xs > this_nuc_xsec_max){
497  phi_at_xsec_max = phi;
498  costh_at_xsec_max = costh;
499  this_nuc_xsec_max = xs;
500  }
501  //
502  } // Done with phi scan
503  }// Done with centre-of-mass angles coarsely
504 
505  // Calculate the range for the next layer
506  costh_range_min = costh_at_xsec_max - costh_increment;
507  costh_range_max = costh_at_xsec_max + costh_increment;
508  phi_range_min = phi_at_xsec_max - phi_increment;
509  phi_range_max = phi_at_xsec_max + phi_increment;
510 
511  double improvement_factor = this_nuc_xsec_max/last_layer_xsec_max;
512  if (ilayer && (improvement_factor-1) < acceptable_fraction_of_safety_factor * (fSafetyFactor-1)) {
513  break;
514  }
515  }
516  if (this_nuc_xsec_max > xsec_max) {
517  xsec_max = this_nuc_xsec_max;
518  LOG("DMELEvent", pINFO) << "best estimate for xsec_max = " << xsec_max;
519  }
520 
521  delete interaction;
522  }
523  // Apply safety factor, since value retrieved from the cache might
524  // correspond to a slightly different energy
525  xsec_max *= fSafetyFactor;
526 
527 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
528  SLOG("DMELEvent", pDEBUG) << interaction->AsString();
529  SLOG("DMELEvent", pDEBUG) << "Max xsec in phase space = " << max_xsec;
530  SLOG("DMELEvent", pDEBUG) << "Computed using alg = " << *fXSecModel;
531 #endif
532 
533  LOG("DMELEvent", pINFO) << "Computed maximum cross section to throw against - value is " << xsec_max;
534  return xsec_max;
535 }
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
DMELEvGen_BindingMode_t fHitNucleonBindingMode
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
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
double ComputeFullDMELPXSec(Interaction *interaction, const NuclearModelI *nucl_model, const XSecAlgorithmI *xsec_model, double cos_theta_0, double phi_0, double &Eb, DMELEvGen_BindingMode_t hitNucleonBindingMode, double min_angle_EM=0., bool bind_nucleon=true)
Definition: DMELUtils.cxx:94
virtual int HitNucleonPosition(void) const
Definition: GHepRecord.cxx:418
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
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 NuclearModelI * fNuclModel
nuclear model
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
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
void ProcessEventRecord(GHepRecord *event_rec) const
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 AddTargetNucleusRemnant(GHepRecord *evrec) const
add a recoiled nucleus remnant
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
double ComputeMaxXSec(const Interaction *in) const
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 Configure(const Registry &config)
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
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
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
DMELEvGen_BindingMode_t StringToDMELBindingMode(const std::string &mode_str)
Definition: DMELUtils.cxx:195
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
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345