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::DMELEventGenerator Class Reference

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

#include <DMELEventGenerator.h>

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

Public Member Functions

 DMELEventGenerator ()
 
 DMELEventGenerator (string config)
 
 ~DMELEventGenerator ()
 
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
 
DMELEvGen_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 DMEL 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 34 of file DMELEventGenerator.h.

Constructor & Destructor Documentation

DMELEventGenerator::DMELEventGenerator ( )

Definition at line 50 of file DMELEventGenerator.cxx.

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

Definition at line 56 of file DMELEventGenerator.cxx.

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

Definition at line 62 of file DMELEventGenerator.cxx.

63 {
64 
65 }

Member Function Documentation

void DMELEventGenerator::AddTargetNucleusRemnant ( GHepRecord evrec) const
private

add a recoiled nucleus remnant

Definition at line 284 of file DMELEventGenerator.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().

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 }
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 DMELEventGenerator::ComputeMaxXSec ( const Interaction in) const
privatevirtual

Implements genie::KineGeneratorWithCache.

Definition at line 404 of file DMELEventGenerator.cxx.

References genie::Interaction::AsString(), genie::utils::BindHitNucleon(), genie::utils::ComputeFullDMELPXSec(), 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().

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 }
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
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
void BindHitNucleon(Interaction &interaction, const NuclearModelI &nucl_model, double &Eb, DMELEvGen_BindingMode_t hitNucleonBindingMode)
Definition: DMELUtils.cxx:259
string AsString(void) const
Summary information for an interaction.
Definition: Interaction.h:56
const NuclearModelI * fNuclModel
nuclear model
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 DMELEventGenerator::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 353 of file DMELEventGenerator.cxx.

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

354 {
355  Algorithm::Configure(config);
356  this->LoadConfig();
357 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void DMELEventGenerator::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 359 of file DMELEventGenerator.cxx.

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

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

Definition at line 365 of file DMELEventGenerator.cxx.

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

Referenced by Configure().

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 }
bool fGenerateUniformly
uniform over allowed phase space + event weight?
DMELEvGen_BindingMode_t fHitNucleonBindingMode
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
const NuclearModelI * fNuclModel
nuclear model
string RgKey
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
bool GetParamDef(const RgKey &name, T &p, const T &def) const
DMELEvGen_BindingMode_t StringToDMELBindingMode(const std::string &mode_str)
Definition: DMELUtils.cxx:195
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
void DMELEventGenerator::ProcessEventRecord ( GHepRecord event_rec) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 67 of file DMELEventGenerator.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::ComputeFullDMELPXSec(), 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::kPSDMELEvGen, 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::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("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 }
int Z(void) const
bool fGenerateUniformly
uniform over allowed phase space + event weight?
DMELEvGen_BindingMode_t fHitNucleonBindingMode
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
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
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
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 NuclearModelI * fNuclModel
nuclear model
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
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 AddTargetNucleusRemnant(GHepRecord *evrec) const
add a recoiled nucleus remnant
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
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

Member Data Documentation

double genie::DMELEventGenerator::fEb
mutableprivate

Definition at line 51 of file DMELEventGenerator.h.

Referenced by ProcessEventRecord().

DMELEvGen_BindingMode_t genie::DMELEventGenerator::fHitNucleonBindingMode
private

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

Definition at line 64 of file DMELEventGenerator.h.

Referenced by LoadConfig(), and ProcessEventRecord().

int genie::DMELEventGenerator::fMaxXSecNucleonThrows
private

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

Definition at line 68 of file DMELEventGenerator.h.

Referenced by LoadConfig().

double genie::DMELEventGenerator::fMinAngleEM
mutableprivate

Definition at line 60 of file DMELEventGenerator.h.

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

const NuclearModelI* genie::DMELEventGenerator::fNuclModel
private

nuclear model

Definition at line 58 of file DMELEventGenerator.h.

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


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