GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::QELEventGenerator Class Reference

Generates values for the kinematic variables describing QEL neutrino interaction events. Is a concrete implementation of the EventRecordVisitorI interface. More...

#include <QELEventGenerator.h>

Inheritance diagram for genie::QELEventGenerator:
Inheritance graph
[legend]
Collaboration diagram for genie::QELEventGenerator:
Collaboration graph
[legend]

Public Member Functions

 QELEventGenerator ()
 
 QELEventGenerator (string config)
 
 ~QELEventGenerator ()
 
void ProcessEventRecord (GHepRecord *event_rec) const
 
void Configure (const Registry &config)
 
void Configure (string config)
 
- Public Member Functions inherited from genie::EventRecordVisitorI
virtual ~EventRecordVisitorI ()
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 
virtual const RegistryGetConfig (void) const
 
RegistryGetOwnedConfig (void)
 
virtual const AlgIdId (void) const
 Get algorithm ID. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
virtual void SetId (string name, string config)
 
const AlgorithmSubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Private Member Functions

void LoadConfig (void)
 
double ComputeMaxXSec (const Interaction *in) const
 
void AddTargetNucleusRemnant (GHepRecord *evrec) const
 add a recoiled nucleus remnant More...
 

Private Attributes

double fEb
 
const NuclearModelIfNuclModel
 nuclear model More...
 
double fMinAngleEM
 
QELEvGen_BindingMode_t fHitNucleonBindingMode
 
int fMaxXSecNucleonThrows
 

Additional Inherited Members

- Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
 
static string BuildParamVectSizeKey (const std::string &comm_name)
 
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
 
static string BuildParamMatRowSizeKey (const std::string &comm_name)
 
static string BuildParamMatColSizeKey (const std::string &comm_name)
 
- Protected Member Functions inherited from genie::KineGeneratorWithCache
 KineGeneratorWithCache ()
 
 KineGeneratorWithCache (string name)
 
 KineGeneratorWithCache (string name, string config)
 
 ~KineGeneratorWithCache ()
 
virtual double ComputeMaxXSec (const Interaction *in, const int nkey) const
 
virtual double MaxXSec (GHepRecord *evrec, const int nkey=0) const
 
virtual double FindMaxXSec (const Interaction *in, const int nkey=0) const
 
virtual void CacheMaxXSec (const Interaction *in, double xsec, const int nkey=0) const
 
virtual double Energy (const Interaction *in) const
 
virtual CacheBranchFxAccessCacheBranch (const Interaction *in, const int nkey=0) const
 
virtual void AssertXSecLimits (const Interaction *in, double xsec, double xsec_max) const
 
- Protected Member Functions inherited from genie::EventRecordVisitorI
 EventRecordVisitorI ()
 
 EventRecordVisitorI (string name)
 
 EventRecordVisitorI (string name, string config)
 
- Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 
 Algorithm (string name)
 
 Algorithm (string name, string config)
 
void Initialize (void)
 
void DeleteConfig (void)
 
void DeleteSubstructure (void)
 
RegistryExtractLocalConfig (const Registry &in) const
 
RegistryExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...
 
template<class T >
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
 
template<class T >
bool GetParamDef (const RgKey &name, T &p, const T &def) const
 
template<class T >
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters. More...
 
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
template<class T >
int GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
 Handle to load matrix of parameters. More...
 
template<class T >
int GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
 
int GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 
- Protected Attributes inherited from genie::KineGeneratorWithCache
const XSecAlgorithmIfXSecModel
 
double fSafetyFactor
 ComputeMaxXSec -> ComputeMaxXSec * fSafetyFactor. More...
 
std::vector< double > vSafetyFactors
 MaxXSec -> MaxXSec * fSafetyFactors[nkey]. More...
 
int fNumOfSafetyFactors
 Number of given safety factors. More...
 
std::vector< string > vInterpolatorTypes
 Type of interpolator for each key in a branch. More...
 
int fNumOfInterpolatorTypes
 Number of given interpolators types. More...
 
double fMaxXSecDiffTolerance
 max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec More...
 
double fEMin
 min E for which maxxsec is cached - forcing explicit calc. More...
 
bool fGenerateUniformly
 uniform over allowed phase space + event weight? More...
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
vector< Registry * > fConfVect
 
vector< bool > fOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

Generates values for the kinematic variables describing QEL neutrino interaction events. Is a concrete implementation of the EventRecordVisitorI interface.

Author
Andrew Furmanski
Created:
August 04, 2014
License:
Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 30 of file QELEventGenerator.h.

Constructor & Destructor Documentation

QELEventGenerator::QELEventGenerator ( )

Definition at line 50 of file QELEventGenerator.cxx.

50  :
51  KineGeneratorWithCache("genie::QELEventGenerator")
52 {
53 
54 }
QELEventGenerator::QELEventGenerator ( string  config)

Definition at line 56 of file QELEventGenerator.cxx.

56  :
57  KineGeneratorWithCache("genie::QELEventGenerator", config)
58 {
59 
60 }
QELEventGenerator::~QELEventGenerator ( )

Definition at line 62 of file QELEventGenerator.cxx.

63 {
64 
65 }

Member Function Documentation

void QELEventGenerator::AddTargetNucleusRemnant ( GHepRecord evrec) const
private

add a recoiled nucleus remnant

Definition at line 289 of file QELEventGenerator.cxx.

References genie::units::A, genie::GHepParticle::A(), genie::GHepRecord::AddParticle(), genie::GHepParticle::E(), genie::PDGLibrary::Find(), genie::GHepParticle::FirstDaughter(), genie::PDGLibrary::Instance(), genie::pdg::IonPdgCode(), genie::pdg::IsNeutron(), genie::pdg::IsProton(), genie::kIStStableFinalState, genie::GHepParticle::LastDaughter(), LOG, genie::GHepParticle::Mass(), genie::GHepRecord::Particle(), genie::GHepParticle::Pdg(), pFATAL, pINFO, genie::GHepParticle::Px(), genie::GHepParticle::Py(), genie::GHepParticle::Pz(), genie::GHepRecord::TargetNucleus(), genie::GHepRecord::TargetNucleusPosition(), and genie::GHepParticle::Z().

Referenced by ProcessEventRecord().

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 }
int Z(void) const
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
double E(void) const
Get energy.
Definition: GHepParticle.h:91
int FirstDaughter(void) const
Definition: GHepParticle.h:68
#define pFATAL
Definition: Messenger.h:56
double Mass(void) const
Mass that corresponds to the PDG code.
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
int Pdg(void) const
Definition: GHepParticle.h:63
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
int LastDaughter(void) const
Definition: GHepParticle.h:69
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#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
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
#define pINFO
Definition: Messenger.h:62
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
int A(void) const
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
virtual int TargetNucleusPosition(void) const
Definition: GHepRecord.cxx:370
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
double QELEventGenerator::ComputeMaxXSec ( const Interaction in) const
privatevirtual

Implements genie::KineGeneratorWithCache.

Definition at line 409 of file QELEventGenerator.cxx.

References genie::Interaction::AsString(), genie::utils::BindHitNucleon(), genie::utils::ComputeFullQELPXSec(), genie::utils::CosTheta0Max(), e, fMinAngleEM, fNuclModel, genie::KineGeneratorWithCache::fSafetyFactor, genie::KineGeneratorWithCache::fXSecModel, genie::Target::HitNucPosition(), genie::Interaction::InitState(), genie::kISkipKinematicChk, genie::kISkipProcessChk, genie::kOnShell, LOG, pDEBUG, pINFO, genie::NuclearModelI::Prob(), genie::NuclearModelI::SetMomentum3(), genie::NuclearModelI::SetRemovalEnergy(), SLOG, and genie::InitialState::Tgt().

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 }
const NuclearModelI * fNuclModel
nuclear model
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
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
string AsString(void) const
Summary information for an interaction.
Definition: Interaction.h:56
const double e
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual double Prob(double p, double w, const Target &) const =0
#define pINFO
Definition: Messenger.h:62
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
void SetRemovalEnergy(double E) const
Definition: NuclearModelI.h:90
void SetMomentum3(const TVector3 &mom) const
Definition: NuclearModelI.h:86
double CosTheta0Max(const genie::Interaction &interaction)
Definition: DMELUtils.cxx:217
double HitNucPosition(void) const
Definition: Target.h:89
const InitialState & InitState(void) const
Definition: Interaction.h:69
const Target & Tgt(void) const
Definition: InitialState.h:66
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
#define pDEBUG
Definition: Messenger.h:63
void QELEventGenerator::Configure ( const Registry config)
virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 358 of file QELEventGenerator.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

359 {
360  Algorithm::Configure(config);
361  this->LoadConfig();
362 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void QELEventGenerator::Configure ( string  config)
virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 364 of file QELEventGenerator.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

365 {
366  Algorithm::Configure(config);
367  this->LoadConfig();
368 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void QELEventGenerator::LoadConfig ( void  )
private

Definition at line 370 of file QELEventGenerator.cxx.

References genie::KineGeneratorWithCache::fEMin, genie::KineGeneratorWithCache::fGenerateUniformly, fHitNucleonBindingMode, genie::KineGeneratorWithCache::fMaxXSecDiffTolerance, fMaxXSecNucleonThrows, fMinAngleEM, fNuclModel, genie::KineGeneratorWithCache::fSafetyFactor, genie::Algorithm::GetParamDef(), genie::utils::StringToQELBindingMode(), and genie::Algorithm::SubAlg().

Referenced by Configure().

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 }
bool fGenerateUniformly
uniform over allowed phase space + event weight?
const NuclearModelI * fNuclModel
nuclear model
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec&gt;maxxsec
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
string RgKey
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
Definition: QELUtils.cxx:194
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
bool GetParamDef(const RgKey &name, T &p, const T &def) const
QELEvGen_BindingMode_t fHitNucleonBindingMode
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
void QELEventGenerator::ProcessEventRecord ( GHepRecord event_rec) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 67 of file QELEventGenerator.cxx.

References genie::units::A, genie::GHepParticle::A(), genie::GHepRecord::AddParticle(), AddTargetNucleusRemnant(), genie::KineGeneratorWithCache::AssertXSecLimits(), genie::utils::BindHitNucleon(), genie::XclsTag::CharmHadronPdg(), genie::Kinematics::ClearRunningValues(), genie::utils::ComputeFullQELPXSec(), genie::utils::CosTheta0Max(), genie::EventGeneratorI::CrossSectionAlg(), genie::GHepRecord::EventFlags(), genie::Interaction::ExclTag(), fEb, genie::KineGeneratorWithCache::fGenerateUniformly, fHitNucleonBindingMode, genie::PDGLibrary::Find(), fMinAngleEM, fNuclModel, genie::Kinematics::FSLeptonP4(), genie::Interaction::FSPrimLeptonPdg(), genie::KineGeneratorWithCache::fXSecModel, genie::gAbortingInErr, genie::NuclearModelI::GenerateNucleon(), gQ2, genie::Kinematics::HadSystP4(), genie::GHepRecord::HitNucleon(), genie::GHepRecord::HitNucleonPosition(), genie::Target::HitNucP4(), genie::Interaction::InitState(), genie::Interaction::InitStatePtr(), genie::RunningThreadInfo::Instance(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::pdg::IonPdgCode(), genie::XclsTag::IsCharmEvent(), genie::Target::IsNucleus(), genie::pdg::IsProton(), genie::XclsTag::IsStrangeEvent(), genie::kIAssumeFreeNucleon, genie::Interaction::KinePtr(), genie::kISkipKinematicChk, genie::kISkipProcessChk, genie::kIStHadronInTheNucleus, genie::kIStStableFinalState, genie::kKineGenErr, genie::constants::kPi, genie::kPSQELEvGen, genie::kRfHitNucRest, genie::controls::kRjMaxIterations, LOG, genie::KineGeneratorWithCache::MaxXSec(), pDEBUG, genie::GHepParticle::Pdg(), pFATAL, pINFO, pNOTICE, genie::GHepRecord::Probe(), genie::InitialState::ProbeE(), genie::GHepRecord::ProbePosition(), pWARN, genie::Kinematics::Q2(), genie::Interaction::RecoilNucleonPdg(), genie::RandomGen::RndKine(), genie::RunningThreadInfo::RunningThread(), genie::GHepRecord::SetDiffXSec(), genie::Target::SetHitNucPosition(), genie::GHepParticle::SetMomentum(), genie::NuclearModelI::SetMomentum3(), genie::utils::SetPrimaryLeptonPolarization(), genie::Kinematics::SetQ2(), genie::exceptions::EVGThreadException::SetReason(), genie::NuclearModelI::SetRemovalEnergy(), genie::GHepParticle::SetRemovalEnergy(), genie::Kinematics::SetW(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::XclsTag::StrangeHadronPdg(), genie::GHepRecord::Summary(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), genie::GHepRecord::TargetNucleus(), genie::InitialState::Tgt(), genie::InitialState::TgtPtr(), genie::utils::kinematics::WQ2toXY(), genie::GHepParticle::X4(), and genie::GHepParticle::Z().

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 }
int Z(void) const
bool fGenerateUniformly
uniform over allowed phase space + event weight?
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
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
Defines the EventGeneratorI interface.
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
void SetHitNucPosition(double r)
Definition: Target.cxx:210
void SetMomentum(const TLorentzVector &p4)
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
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
bool IsCharmEvent(void) const
Definition: XclsTag.h:50
int Pdg(void) const
Definition: GHepParticle.h:63
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:91
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
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
#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
void SetMomentum3(const TVector3 &mom) const
Definition: NuclearModelI.h:86
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
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
Target * TgtPtr(void) const
Definition: InitialState.h:67
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
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
const Target & Tgt(void) const
Definition: InitialState.h:66
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
const EventGeneratorI * RunningThread(void)
virtual bool GenerateNucleon(const Target &) const =0
bool gAbortingInErr
Definition: Messenger.cxx:34
void SetPrimaryLeptonPolarization(GHepRecord *ev)
double ProbeE(RefFrame_t rf) const
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...
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
enum genie::EGHepStatus GHepStatus_t
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
QELEvGen_BindingMode_t fHitNucleonBindingMode

Member Data Documentation

double genie::QELEventGenerator::fEb
mutableprivate

Definition at line 47 of file QELEventGenerator.h.

Referenced by ProcessEventRecord().

QELEvGen_BindingMode_t genie::QELEventGenerator::fHitNucleonBindingMode
private

Enum that indicates which approach should be used to handle the binding energy of the struck nucleon

Definition at line 60 of file QELEventGenerator.h.

Referenced by LoadConfig(), and ProcessEventRecord().

int genie::QELEventGenerator::fMaxXSecNucleonThrows
private

The number of nucleons to sample from the nuclear model when choosing a maximum momentum to use in ComputeMaxXSec()

Definition at line 64 of file QELEventGenerator.h.

Referenced by LoadConfig().

double genie::QELEventGenerator::fMinAngleEM
mutableprivate

Definition at line 56 of file QELEventGenerator.h.

Referenced by ComputeMaxXSec(), LoadConfig(), and ProcessEventRecord().

const NuclearModelI* genie::QELEventGenerator::fNuclModel
private

nuclear model

Definition at line 54 of file QELEventGenerator.h.

Referenced by ComputeMaxXSec(), LoadConfig(), and ProcessEventRecord().


The documentation for this class was generated from the following files: