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