Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | List of all members
genie::MECGenerator Class Reference

Simulate the primary MEC interaction. More...

#include "/cvmfs/nova.opensciencegrid.org/externals/genie/v3_00_06_p01/Linux64bit+3.10-2.17-e17-debug/GENIE-Generator/src/Physics/Multinucleon/EventGen/MECGenerator.h"

Inheritance diagram for genie::MECGenerator:
genie::EventRecordVisitorI genie::Algorithm

Public Member Functions

 MECGenerator ()
 
 MECGenerator (string config)
 
 ~MECGenerator ()
 
void ProcessEventRecord (GHepRecord *event) const
 
void Configure (const Registry &config)
 
void Configure (string config)
 
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...
 

Protected Member Functions

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 >
bool GetParamVect (const std::string &comm_name, std::vector< T > &v, unsigned int max, 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

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...
 

Private Member Functions

void LoadConfig (void)
 
void AddNucleonCluster (GHepRecord *event) const
 
void AddTargetRemnant (GHepRecord *event) const
 
void GenerateFermiMomentum (GHepRecord *event) const
 
void SelectEmpiricalKinematics (GHepRecord *event) const
 
void AddFinalStateLepton (GHepRecord *event) const
 
void RecoilNucleonCluster (GHepRecord *event) const
 
void DecayNucleonCluster (GHepRecord *event) const
 
void SelectNSVLeptonKinematics (GHepRecord *event) const
 
void GenerateNSVInitialHadrons (GHepRecord *event) const
 
PDGCodeList NucleonClusterConstituents (int pdgc) const
 

Private Attributes

const XSecAlgorithmIfXSecModel
 
TGenPhaseSpace fPhaseSpaceGenerator
 
const NuclearModelIfNuclModel
 
double fQ3Max
 

Detailed Description

Simulate the primary MEC interaction.

Author
Costas Andreopoulos <costas.andreopoulos stfc.ac.uk> University of Liverpool & STFC Rutherford Appleton Lab

Steve Dytman <dytman+ pitt.edu> Pittsburgh University

Sep. 22, 2008

Copyright (c) 2003-2019, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Definition at line 35 of file MECGenerator.h.

Constructor & Destructor Documentation

MECGenerator::MECGenerator ( )

Definition at line 56 of file MECGenerator.cxx.

56  :
57 EventRecordVisitorI("genie::MECGenerator")
58 {
59 
60 }
MECGenerator::MECGenerator ( string  config)

Definition at line 62 of file MECGenerator.cxx.

62  :
63 EventRecordVisitorI("genie::MECGenerator", config)
64 {
65 
66 }
Definition: config.py:1
MECGenerator::~MECGenerator ( )

Definition at line 68 of file MECGenerator.cxx.

69 {
70 
71 }

Member Function Documentation

void MECGenerator::AddFinalStateLepton ( GHepRecord event) const
private

Definition at line 319 of file MECGenerator.cxx.

References beta, genie::Interaction::FSPrimLepton(), genie::Target::HitNucP4(), genie::Interaction::InitState(), genie::RandomGen::Instance(), genie::Interaction::Kine(), genie::kIStStableFinalState, kPi, genie::kRfHitNucRest, LOG, pNOTICE, genie::GHepRecord::Probe(), genie::InitialState::ProbeE(), genie::Kinematics::Q2(), genie::utils::kinematics::Q2(), generate_hists::rnd, genie::RandomGen::RndLep(), ana::Sqrt(), genie::InitialState::Tgt(), genie::GHepParticle::X4(), genie::Kinematics::y(), and submit_syst::y.

Referenced by ProcessEventRecord().

320 {
321 // Add the final-state primary lepton in the event record.
322 // Compute its 4-momentum based on the selected interaction kinematics.
323 //
324  Interaction * interaction = event->Summary();
325 
326  // The boost back to the lab frame was missing, that is included now with the introduction of the beta factor
327  const InitialState & init_state = interaction->InitState();
328  const TLorentzVector & pnuc4 = init_state.Tgt().HitNucP4(); //[@LAB]
329  TVector3 beta = pnuc4.BoostVector();
330 
331  // Boosting the incoming neutrino to the NN-cluster rest frame
332  // Neutrino 4p
333  TLorentzVector * p4v = event->Probe()->GetP4(); // v 4p @ LAB
334  p4v->Boost(-1.*beta); // v 4p @ NN-cluster rest frame
335 
336  // Look-up selected kinematics
337  double Q2 = interaction->Kine().Q2(true);
338  double y = interaction->Kine().y(true);
339 
340  // Auxiliary params
341  double Ev = interaction->InitState().ProbeE(kRfHitNucRest);
342  LOG("MEC", pNOTICE) << "neutrino energy = " << Ev;
343  double ml = interaction->FSPrimLepton()->Mass();
344  double ml2 = TMath::Power(ml,2);
345 
346  // Compute the final state primary lepton energy and momentum components
347  // along and perpendicular the neutrino direction
348  double El = (1-y)*Ev;
349  double plp = El - 0.5*(Q2+ml2)/Ev; // p(//)
350  double plt = TMath::Sqrt(TMath::Max(0.,El*El-plp*plp-ml2)); // p(-|)
351 
352  LOG("MEC", pNOTICE)
353  << "fsl: E = " << El << ", |p//| = " << plp << ", |pT| = " << plt;
354 
355  // Randomize transverse components
357  double phi = 2*kPi * rnd->RndLep().Rndm();
358  double pltx = plt * TMath::Cos(phi);
359  double plty = plt * TMath::Sin(phi);
360 
361  // Take a unit vector along the neutrino direction
362  // WE NEED THE UNIT VECTOR ALONG THE NEUTRINO DIRECTION IN THE NN-CLUSTER REST FRAME, NOT IN THE LAB FRAME
363  TVector3 unit_nudir = p4v->Vect().Unit(); //We use this one, which is in the NN-cluster rest frame
364  // Rotate lepton momentum vector from the reference frame (x'y'z') where
365  // {z':(neutrino direction), z'x':(theta plane)} to the LAB
366  TVector3 p3l(pltx,plty,plp);
367  p3l.RotateUz(unit_nudir);
368 
369  // Lepton 4-momentum in LAB
370  TLorentzVector p4l(p3l,El);
371 
372  // Boost final state primary lepton to the lab frame
373  p4l.Boost(beta); // active Lorentz transform
374 
375  // Figure out the final-state primary lepton PDG code
376  int pdgc = interaction->FSPrimLepton()->PdgCode();
377 
378  // Lepton 4-position (= interacton vtx)
379  TLorentzVector v4(*event->Probe()->X4());
380 
381  int momidx = event->ProbePosition();
382  event->AddParticle(
383  pdgc, kIStStableFinalState, momidx, -1, -1, -1, p4l, v4);
384 }
const double kPi
TRandom3 & RndLep(void) const
rnd number generator used by final state primary lepton generators
Definition: RandomGen.h:63
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
Double_t beta
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
double y(bool selected=false) const
Definition: Kinematics.cxx:122
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:92
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const Kinematics & Kine(void) const
Definition: Interaction.h:71
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
const InitialState & InitState(void) const
Definition: Interaction.h:69
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
#define pNOTICE
Definition: Messenger.h:62
const Target & Tgt(void) const
Definition: InitialState.h:67
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:49
int Algorithm::AddLowRegistry ( Registry rp,
bool  owns = true 
)
protectedinherited

add registry with lowest priority, also update ownership

Definition at line 601 of file Algorithm.cxx.

Referenced by genie::EventGenerator::Configure().

601  {
602 
603  fConfVect.push_back( rp ) ;
604  fOwnerships.push_back( own ) ;
605 
606  if ( fConfig ) {
607  delete fConfig ;
608  fConfig = 0 ;
609  }
610 
611  return fConfVect.size() ;
612 
613 }
vector< Registry * > fConfVect
Definition: Algorithm.h:161
vector< bool > fOwnerships
ownership for every registry in fConfVect
Definition: Algorithm.h:164
Registry * fConfig
Summary configuration derived from fConvVect, not necessarily allocated.
Definition: Algorithm.h:194
void genie::MECGenerator::AddNucleonCluster ( GHepRecord event) const
private
void MECGenerator::AddTargetRemnant ( GHepRecord event) const
private

Definition at line 110 of file MECGenerator.cxx.

References genie::units::A, genie::GHepParticle::A(), genie::pdg::IonPdgCode(), genie::kIStStableFinalState, genie::kPdgClusterNN, genie::kPdgClusterNP, genie::kPdgClusterPP, genie::GHepParticle::P4(), make_associated_cosmic_defs::p4, genie::GHepParticle::Pdg(), Z, and genie::GHepParticle::Z().

Referenced by ProcessEventRecord().

111 {
112 // Add the remnant nucleus (= initial nucleus - nucleon cluster) in the
113 // event record.
114 
115  GHepParticle * target = event->TargetNucleus();
116  GHepParticle * cluster = event->HitNucleon();
117 
118  int Z = target->Z();
119  int A = target->A();
120 
121  if(cluster->Pdg() == kPdgClusterNN) { A-=2; ; }
122  if(cluster->Pdg() == kPdgClusterNP) { A-=2; Z-=1; }
123  if(cluster->Pdg() == kPdgClusterPP) { A-=2; Z-=2; }
124 
125  int ipdgc = pdg::IonPdgCode(A, Z);
126 
127  const TLorentzVector & p4cluster = *(cluster->P4());
128  const TLorentzVector & p4tgt = *(target ->P4());
129 
130  const TLorentzVector p4 = p4tgt - p4cluster;
131  const TLorentzVector v4(0.,0.,0., 0.);
132 
133  int momidx = event->TargetNucleusPosition();
134 
135  /*
136  if( A == 2 && Z == 2){
137  // residual nucleus was just two protons, not a nucleus we know.
138  // this might not conserve energy...
139  event->AddParticle(kPdgProton,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
140  // v4,v4 because I'm lazy, give the four momentum to one of the protons, not the other
141  event->AddParticle(kPdgProton,kIStStableFinalState, momidx,-1,-1,-1, v4,v4);
142  } else if ( A == 2 && Z == 0){
143  // residual nucleus was just two neutrons, not a nucleus we know.
144  // this might not conserve energy...
145  event->AddParticle(kPdgNeutron,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
146  // v4,v4 because I'm lazy, give the four momentum to one of the protons, not the other
147  event->AddParticle(kPdgNeutron,kIStStableFinalState, momidx,-1,-1,-1, v4,v4);
148  } else {
149  // regular nucleus, including deuterium
150  event->AddParticle(ipdgc,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
151  }
152  */
153 
154  event->AddParticle(ipdgc,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
155 
156 }
int Z(void) const
const XML_Char * target
Definition: expat.h:268
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
const int kPdgClusterNP
Definition: PDGCodes.h:192
const int kPdgClusterNN
Definition: PDGCodes.h:191
Float_t Z
Definition: plot.C:38
int Pdg(void) const
Definition: GHepParticle.h:64
static const double A
Definition: Units.h:82
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
int A(void) const
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
const int kPdgClusterPP
Definition: PDGCodes.h:193
int Algorithm::AddTopRegisties ( const vector< Registry * > &  rs,
bool  owns = false 
)
protectedinherited

Add registries with top priority, also udated Ownerships.

Definition at line 653 of file Algorithm.cxx.

653  {
654 
655  fConfVect.insert( fConfVect.begin(), rs.begin(), rs.end() ) ;
656 
657  fOwnerships.insert( fOwnerships.begin(), rs.size(), own ) ;
658 
659  if ( fConfig ) {
660  delete fConfig ;
661  fConfig = 0 ;
662  }
663 
664  return fConfVect.size() ;
665 
666 }
vector< Registry * > fConfVect
Definition: Algorithm.h:161
vector< bool > fOwnerships
ownership for every registry in fConfVect
Definition: Algorithm.h:164
Registry * fConfig
Summary configuration derived from fConvVect, not necessarily allocated.
Definition: Algorithm.h:194
int Algorithm::AddTopRegistry ( Registry rp,
bool  owns = true 
)
protectedinherited

add registry with top priority, also update ownership

Definition at line 585 of file Algorithm.cxx.

Referenced by genie::EventGeneratorListAssembler::AssembleGeneratorList().

585  {
586 
587  fConfVect.insert( fConfVect.begin(), rp ) ;
588  fOwnerships.insert( fOwnerships.begin(), own ) ;
589 
590  if ( fConfig ) {
591  delete fConfig ;
592  fConfig = 0 ;
593  }
594 
595  return fConfVect.size() ;
596 
597 }
vector< Registry * > fConfVect
Definition: Algorithm.h:161
vector< bool > fOwnerships
ownership for every registry in fConfVect
Definition: Algorithm.h:164
Registry * fConfig
Summary configuration derived from fConvVect, not necessarily allocated.
Definition: Algorithm.h:194
void Algorithm::AdoptConfig ( void  )
inherited

Clone the configuration registry looked up from the configuration pool and take its ownership

Definition at line 394 of file Algorithm.cxx.

References Configure(), GetConfig(), LOG, and pNOTICE.

Referenced by genie::Algorithm::AllowReconfig().

394  {
395 
396  LOG("Algorithm", pNOTICE)
397  << this->Id().Key() << " is taking ownership of its configuration";
398 
399  // if(fOwnsConfig) {
400  // LOG("Algorithm", pWARN)
401  // << this->Id().Key() << " already owns its configuration!";
402  // return;
403  // }
404 
405  this->Configure( GetConfig() );
406 }
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
#define pNOTICE
Definition: Messenger.h:62
string Key(void) const
Definition: AlgId.h:47
void Algorithm::AdoptSubstructure ( void  )
inherited

Take ownership of the algorithms subtructure (sub-algorithms,...) by copying them from the AlgFactory pool to the local pool Also bring all the configuration variables to the top level config Registry. This can be used to group together a series of algorithms & their configurations and extract (a clone of) this group from the shared pools. Having a series of algorithms/configurations behaving as a monolithic block, with a single point of configuration (the top level) is to be used when bits & pieces of GENIE are used in isolation for data fitting or reweighting

Definition at line 408 of file Algorithm.cxx.

References genie::AlgFactory::AdoptAlgorithm(), genie::Algorithm::AdoptSubstructure(), GetConfig(), genie::AlgFactory::Instance(), genie::kRgAlg, LOG, pDEBUG, pNOTICE, and genie::RegistryItemI::TypeInfo().

Referenced by genie::Algorithm::AdoptSubstructure(), genie::Algorithm::AllowReconfig(), main(), and testReconfigInOwnedModules().

409 {
410 // Take ownership of the algorithms subtructure (sub-algorithms,..) by copying
411 // them from the AlgFactory pool to the local pool. Also bring all the
412 // configuration variables to the top level. See the header for more details.
413 // A secial naming convention is required for configuration parameter keys
414 // for parameters belonging to sub-algorithms (or sub-algorithms of these
415 // sub-algorithms and so on...).
416 // The convention is: "sub-alg-key/sub-sub-alg-key/.../original name"
417 // This is a recursive method: The AdoptSubtructure()of all sub-algorithms is
418 // invoked.
419 //
420  LOG("Algorithm", pNOTICE)
421  << "Algorithm: " << this->Id().Key() << " is adopting its substructure";
422 
423 // Registry deep_config;
424 // deep_config.UnLock();
425 // deep_config.SetName(this->Id().Key());
426 
427  // deep_config.SetName(this->Id().Config() + ";D");
428  // fID.SetConfig(this->Id().Config() + ";D");
429 
430  if(fOwnsSubstruc) this->DeleteSubstructure();
431 
432  fOwnedSubAlgMp = new AlgMap;
433  fOwnsSubstruc = true;
434 
435  AlgFactory * algf = AlgFactory::Instance();
436 
437  const RgIMap & rgmap = GetConfig().GetItemMap();
438 
439  RgIMapConstIter iter = rgmap.begin();
440  for( ; iter != rgmap.end(); ++iter) {
441 
442  RgKey reg_key = iter->first;
443  RegistryItemI * ri = iter->second;
444 
445  if(ri->TypeInfo() == kRgAlg) {
446  LOG("Algorithm", pDEBUG)
447  << "Found sub-algorithm pointed to by " << reg_key;
448  RgAlg reg_alg = fConfig->GetAlg(reg_key);
449  AlgId id(reg_alg);
450 
451  LOG("Algorithm", pDEBUG) << "Adopting sub-algorithm = " << id.Key();
452  Algorithm * subalg = algf->AdoptAlgorithm(id.Name(),id.Config());
453  subalg->AdoptSubstructure();
454 
455  LOG("Algorithm", pDEBUG) << "Adding sub-algorithm to local pool";
456  AlgMapPair key_alg_pair(reg_key, subalg);
457  fOwnedSubAlgMp->insert(key_alg_pair);
458 
459  }
460 
461  }
462 
463 
464  if ( fConfig ) {
465  delete fConfig ;
466  fConfig = 0 ;
467  }
468 
469 }
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
void DeleteSubstructure(void)
Definition: Algorithm.cxx:496
AlgMap * fOwnedSubAlgMp
local pool for owned sub-algs (taken out of the factory pool)
Definition: Algorithm.h:167
bool fOwnsSubstruc
true if it owns its substructure (sub-algs,...)
Definition: Algorithm.h:155
Algorithm abstract base class.
Definition: Algorithm.h:54
map< string, Algorithm * > AlgMap
Definition: Algorithm.h:49
Registry item pABC.
Definition: RegistryItemI.h:30
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
virtual RgType_t TypeInfo(void) const =0
map< RgKey, RegistryItemI * >::const_iterator RgIMapConstIter
Definition: Registry.h:50
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const RgIMap & GetItemMap(void) const
Definition: Registry.h:162
void AdoptSubstructure(void)
Definition: Algorithm.cxx:408
pair< string, Algorithm * > AlgMapPair
Definition: Algorithm.h:52
Algorithm * AdoptAlgorithm(const AlgId &algid) const
Definition: AlgFactory.cxx:127
Algorithm ID (algorithm name + configuration set name)
Definition: AlgId.h:35
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
string RgKey
Registry * fConfig
Summary configuration derived from fConvVect, not necessarily allocated.
Definition: Algorithm.h:194
#define pNOTICE
Definition: Messenger.h:62
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
string Key(void) const
Definition: AlgId.h:47
RgAlg GetAlg(RgKey key) const
Definition: Registry.cxx:503
#define pDEBUG
Definition: Messenger.h:64
map< RgKey, RegistryItemI * > RgIMap
Definition: Registry.h:46
virtual bool genie::Algorithm::AllowReconfig ( void  ) const
inlinevirtualinherited
AlgCmp_t Algorithm::Compare ( const Algorithm alg) const
virtualinherited

Compare with input algorithm.

Definition at line 294 of file Algorithm.cxx.

References genie::AlgId::Config(), genie::Algorithm::Id(), genie::kAlgCmpDiffAlg, genie::kAlgCmpDiffConfig, genie::kAlgCmpIdentical, genie::kAlgCmpUnknown, and genie::AlgId::Name().

Referenced by genie::Algorithm::AllowReconfig().

295 {
296 // Compares itself with the input algorithm
297 
298  string alg1 = this->Id().Name();
299  string config1 = this->Id().Config();
300  string alg2 = algo->Id().Name();
301  string config2 = algo->Id().Config();
302 
303  if(alg1 == alg2)
304  {
305  if(config1 == config2) return kAlgCmpIdentical;
306  else return kAlgCmpDiffConfig;
307  }
308  else return kAlgCmpDiffAlg;
309 
310  return kAlgCmpUnknown;
311 }
string Name(void) const
Definition: AlgId.h:45
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
string Config(void) const
Definition: AlgId.h:46
void MECGenerator::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 1058 of file MECGenerator.cxx.

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

1059 {
1060  Algorithm::Configure(config);
1061  this->LoadConfig();
1062 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
void MECGenerator::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 1064 of file MECGenerator.cxx.

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

1065 {
1067  this->LoadConfig();
1068 }
Definition: config.py:1
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
void MECGenerator::DecayNucleonCluster ( GHepRecord event) const
private

accept_decay

Definition at line 425 of file MECGenerator.cxx.

References ana::assert(), genie::PDGLibrary::Find(), fPhaseSpaceGenerator, genie::GHepParticle::GetP4(), genie::GHepParticle::GetX4(), MECModelEnuComparisons::i, genie::PDGLibrary::Instance(), genie::RandomGen::Instance(), genie::kHadroSysGenErr, genie::kIStHadronInTheNucleus, genie::controls::kMaxUnweightDecayIterations, LOG, m, NucleonClusterConstituents(), genie::utils::print::P4AsString(), genie::GHepParticle::Pdg(), pERROR, pINFO, pNOTICE, pWARN, generate_hists::rnd, genie::RandomGen::RndDec(), genie::exceptions::EVGThreadException::SetReason(), genie::exceptions::EVGThreadException::SetReturnStep(), sum, genie::exceptions::EVGThreadException::SwitchOnStepBack(), and w.

Referenced by ProcessEventRecord().

426 {
427 // Perform a phase-space decay of the nucleon cluster and add its decay
428 // products in the event record
429 //
430  LOG("MEC", pINFO) << "Decaying nucleon cluster...";
431 
432  // get di-nucleon cluster
433  int nucleon_cluster_id = 5;
434  GHepParticle * nucleon_cluster = event->Particle(nucleon_cluster_id);
435  assert(nucleon_cluster);
436 
437  // get decay products
438  PDGCodeList pdgv = this->NucleonClusterConstituents(nucleon_cluster->Pdg());
439  LOG("MEC", pINFO) << "Decay product IDs: " << pdgv;
440 
441  // Get the decay product masses
442  vector<int>::const_iterator pdg_iter;
443  int i = 0;
444  double * mass = new double[pdgv.size()];
445  double sum = 0;
446  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
447  int pdgc = *pdg_iter;
448  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
449  mass[i++] = m;
450  sum += m;
451  }
452 
453  LOG("MEC", pINFO)
454  << "Performing a phase space decay to "
455  << pdgv.size() << " particles / total mass = " << sum;
456 
457  TLorentzVector * p4d = nucleon_cluster->GetP4();
458  TLorentzVector * v4d = nucleon_cluster->GetX4();
459 
460  LOG("MEC", pINFO)
461  << "Decaying system p4 = " << utils::print::P4AsString(p4d);
462 
463  // Set the decay
464  bool permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass);
465  if(!permitted) {
466  LOG("MEC", pERROR)
467  << " *** Phase space decay is not permitted \n"
468  << " Total particle mass = " << sum << "\n"
469  << " Decaying system p4 = " << utils::print::P4AsString(p4d);
470  // clean-up
471  delete [] mass;
472  delete p4d;
473  delete v4d;
474  // throw exception
475  event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
477  exception.SetReason("Decay not permitted kinematically");
478  exception.SwitchOnStepBack();
479  exception.SetReturnStep(0);
480  throw exception;
481  }
482 
483  // Get the maximum weight
484  double wmax = -1;
485  for(int idec=0; idec<200; idec++) {
486  double w = fPhaseSpaceGenerator.Generate();
487  wmax = TMath::Max(wmax,w);
488  }
489  assert(wmax>0);
490  wmax *= 2;
491 
492  LOG("MEC", pNOTICE)
493  << "Max phase space gen. weight = " << wmax;
494 
496  bool accept_decay=false;
497  unsigned int itry=0;
498  while(!accept_decay)
499  {
500  itry++;
501 
503  // report, clean-up and return
504  LOG("MEC", pWARN)
505  << "Couldn't generate an unweighted phase space decay after "
506  << itry << " attempts";
507  // clean up
508  delete [] mass;
509  delete p4d;
510  delete v4d;
511  // throw exception
512  event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
514  exception.SetReason("Couldn't select decay after N attempts");
515  exception.SwitchOnStepBack();
516  exception.SetReturnStep(0);
517  throw exception;
518  }
519  double w = fPhaseSpaceGenerator.Generate();
520  if(w > wmax) {
521  LOG("MEC", pWARN)
522  << "Decay weight = " << w << " > max decay weight = " << wmax;
523  }
524  double gw = wmax * rnd->RndDec().Rndm();
525  accept_decay = (gw<=w);
526 
527  LOG("MEC", pINFO)
528  << "Decay weight = " << w << " / R = " << gw
529  << " - accepted: " << accept_decay;
530 
531  } //!accept_decay
532 
533  // Insert the decay products in the event record
534  TLorentzVector v4(*v4d);
536  int idp = 0;
537  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
538  int pdgc = *pdg_iter;
539  TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(idp);
540  event->AddParticle(pdgc, ist, nucleon_cluster_id,-1,-1,-1, *p4fin, v4);
541  idp++;
542  }
543 
544  // Clean-up
545  delete [] mass;
546  delete p4d;
547  delete v4d;
548 }
TLorentzVector * GetX4(void) const
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
enum genie::EGHepStatus GHepStatus_t
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:34
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
A list of PDG codes.
Definition: PDGCodeList.h:33
int Pdg(void) const
Definition: GHepParticle.h:64
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
TLorentzVector * GetP4(void) const
#define pINFO
Definition: Messenger.h:63
PDGCodeList NucleonClusterConstituents(int pdgc) const
#define pWARN
Definition: Messenger.h:61
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
assert(nhit_max >=nhit_nbins)
#define pNOTICE
Definition: Messenger.h:62
Double_t sum
Definition: plot.C:31
Float_t w
Definition: plot.C:20
TRandom3 & RndDec(void) const
rnd number generator used by decay models
Definition: RandomGen.h:57
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
TGenPhaseSpace fPhaseSpaceGenerator
Definition: MECGenerator.h:65
void Algorithm::DeleteConfig ( void  )
protectedinherited

Definition at line 471 of file Algorithm.cxx.

References MECModelEnuComparisons::i.

Referenced by genie::Algorithm::AllowReconfig().

472 {
473  // there is nothing to delete if the configuration is not owned but is
474  // rather looked up from the configuration pool
475  //
476 
477  for ( unsigned int i = 0 ; i < fConfVect.size() ; ++i ) {
478  if ( fOwnerships[i] ) {
479  delete fConfVect[i] ;
480  }
481  }
482 
483  fConfVect.clear() ;
484  fOwnerships.clear() ;
485 
486  // delete owned configuration registry
487 
488  if(fConfig) {
489  delete fConfig;
490  fConfig=0;
491  }
492 
493 }
vector< Registry * > fConfVect
Definition: Algorithm.h:161
vector< bool > fOwnerships
ownership for every registry in fConfVect
Definition: Algorithm.h:164
Registry * fConfig
Summary configuration derived from fConvVect, not necessarily allocated.
Definition: Algorithm.h:194
void Algorithm::DeleteSubstructure ( void  )
protectedinherited

Definition at line 496 of file Algorithm.cxx.

Referenced by genie::Algorithm::AllowReconfig().

497 {
498  // there is nothing to delete if the sub-algorithms are not owned but rather
499  // taken from the AlgFactory's pool
500  //
501  if(!fOwnsSubstruc) return;
502 
503  // delete local algorithm pool
504  //
505  AlgMapIter iter = fOwnedSubAlgMp->begin();
506  for( ; iter != fOwnedSubAlgMp->end(); ++iter) {
507  Algorithm * alg = iter->second;
508  if(alg) {
509  delete alg;
510  alg=0;
511  }
512  }
513  delete fOwnedSubAlgMp;
514  fOwnedSubAlgMp = 0;
515 }
AlgMap * fOwnedSubAlgMp
local pool for owned sub-algs (taken out of the factory pool)
Definition: Algorithm.h:167
bool fOwnsSubstruc
true if it owns its substructure (sub-algs,...)
Definition: Algorithm.h:155
Algorithm abstract base class.
Definition: Algorithm.h:54
map< string, Algorithm * >::iterator AlgMapIter
Definition: Algorithm.h:50
Registry * Algorithm::ExtractLocalConfig ( const Registry in) const
protectedinherited

Split an incoming configuration Registry into a block valid for this algorithm Ownership of the returned registry belongs to the algo

Definition at line 518 of file Algorithm.cxx.

References genie::RegistryItemI::Clone(), genie::Registry::GetItemMap(), genie::Registry::Name(), and confusionMatrixTree::out.

Referenced by genie::Algorithm::AllowReconfig().

518  {
519 
520  const RgIMap & rgmap = in.GetItemMap();
521  Registry * out = new Registry( in.Name(), false );
522 
523  for( RgIMapConstIter reg_iter = rgmap.begin();
524  reg_iter != rgmap.end(); ++reg_iter ) {
525 
526  RgKey reg_key = reg_iter->first;
527  if( reg_key.find( '/' ) != string::npos) continue;
528 
529  // at this point
530  // this key is referred to the local algorithm
531  // it has to be copied in out;
532 
533  RegistryItemI * ri = reg_iter->second;
534  RgIMapPair key_item_pair( reg_key, ri->Clone() );
535  out -> Set(key_item_pair);
536 
537  }
538 
539  if ( out -> NEntries() <= 0 ) {
540  delete out ;
541  out = 0 ;
542  }
543 
544  return out ;
545 }
Registry item pABC.
Definition: RegistryItemI.h:30
string Name(void) const
get the registry name
Definition: Registry.cxx:612
map< RgKey, RegistryItemI * >::const_iterator RgIMapConstIter
Definition: Registry.h:50
const RgIMap & GetItemMap(void) const
Definition: Registry.h:162
pair< RgKey, RegistryItemI * > RgIMapPair
Definition: Registry.h:47
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
virtual RegistryItemI * Clone(void) const =0
map< RgKey, RegistryItemI * > RgIMap
Definition: Registry.h:46
Registry * Algorithm::ExtractLowerConfig ( const Registry in,
const string alg_key 
) const
protectedinherited

Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key.

Definition at line 549 of file Algorithm.cxx.

References genie::RegistryItemI::Clone(), genie::Registry::GetItemMap(), genie::Registry::Name(), and confusionMatrixTree::out.

Referenced by genie::Algorithm::AllowReconfig().

549  {
550 
551  const RgIMap & rgmap = in.GetItemMap();
552  Registry * out = new Registry( in.Name(), false );
553 
554  for( RgIMapConstIter reg_iter = rgmap.begin();
555  reg_iter != rgmap.end(); ++reg_iter ) {
556 
557  RgKey reg_key = reg_iter->first;
558  if( reg_key.find(alg_key+"/") == string::npos) continue;
559 
560  // at this point
561  // this key is referred to the sub-algorithm
562  // indicated by alg_key: it has to be copied in out;
563 
564  int new_key_start = reg_key.find_first_of('/')+1;
565  RgKey new_reg_key = reg_key.substr( new_key_start, reg_key.length() );
566 
567  RegistryItemI * ri = reg_iter->second;
568  RgIMapPair key_item_pair(new_reg_key, ri->Clone());
569  out -> Set(key_item_pair);
570 
571  }
572 
573  if ( out -> NEntries() <= 0 ) {
574  delete out ;
575  out = 0 ;
576  }
577 
578  return out ;
579 
580 }
Registry item pABC.
Definition: RegistryItemI.h:30
string Name(void) const
get the registry name
Definition: Registry.cxx:612
map< RgKey, RegistryItemI * >::const_iterator RgIMapConstIter
Definition: Registry.h:50
const RgIMap & GetItemMap(void) const
Definition: Registry.h:162
pair< RgKey, RegistryItemI * > RgIMapPair
Definition: Registry.h:47
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
virtual RegistryItemI * Clone(void) const =0
map< RgKey, RegistryItemI * > RgIMap
Definition: Registry.h:46
void Algorithm::FindConfig ( void  )
virtualinherited

Lookup configuration from the config pool Similar logic from void Configure(string)

Definition at line 135 of file Algorithm.cxx.

References gen_hdf5record::config, exit(), genie::AlgConfigPool::FindRegistry(), genie::Registry::GetItemMap(), genie::Registry::GetString(), MECModelEnuComparisons::i, genie::AlgConfigPool::Instance(), it, genie::Registry::ItemIsLocal(), parse_dependency_file_t::list, LOG, pDEBUG, pFATAL, time_estimates::pool, pWARN, moon_position_table_new3::second, genie::utils::str::Split(), string, and APDHVSetting::temp.

136 {
137 // Finds its configration Registry from the ConfigPool and gets a pointer to
138 // it. If the Registry comes from the ConfigPool then the Algorithm does not
139 // own its configuration (the ConfigPool singleton keeps the ownership and the
140 // responsibility to -eventually- delete all the Registries it instantiates
141 // by parsing the XML config files).
142 
143  DeleteConfig() ;
144 
146 
147  Registry * config = 0 ;
148 
149  // load the Default config if config is not the default
150  if ( fID.Config() != "Default" ) {
151  config = pool -> FindRegistry( fID.Name(), "Default" );
152  if ( config ) {
153  if ( config -> NEntries() > 0 ) {
154  AddTopRegistry( config, false ) ;
155  LOG("Algorithm", pDEBUG) << "\n" << *config;
156  }
157  }
158  }
159 
160  // Load the right config
161  config = pool->FindRegistry(this);
162 
163  if(!config)
164  // notify & keep whatever config Registry was used before.
165  LOG("Algorithm", pWARN)
166  << "No Configuration available for "
167  << this->Id().Key() << " at the ConfigPool";
168  else {
169  if ( config -> NEntries() > 0 ) {
170  AddTopRegistry( config, false ) ;
171  LOG("Algorithm", pDEBUG) << "\n" << config;
172  }
173  }
174 
175  const string common_key_root = "Common" ;
176  std::map<string, string> common_lists;
177 
178  // Load Common Parameters if key that start with "Common" is found
179  for ( unsigned int i = 0 ; i < fConfVect.size() ; ++i ) {
180  const Registry & temp = * fConfVect[i] ;
181  for ( RgIMapConstIter it = temp.GetItemMap().begin() ; it != temp.GetItemMap().end() ; ++it ) {
182 
183  // check if it is a "Common" entry
184  if ( it -> first.find( common_key_root ) == 0 ) {
185  // retrieve the type of the common entry
186  std::string type = it -> first.substr(common_key_root.size() ) ;
187 
188  if ( temp.ItemIsLocal( it -> first ) ) {
189 
190  string temp_list = temp.GetString( it -> first ) ;
191  if ( temp_list.length() > 0 ) {
192  common_lists[type] = temp_list ;
193  }
194  }
195  }
196 
197  }
198 
199  } // loop over the local registries
200 
201 
202  for ( std::map<string, string>::const_iterator it = common_lists.begin() ;
203  it != common_lists.end() ; ++it ) {
204 
205  vector<string> list = str::Split( it -> second , "," ) ;
206 
207  for ( unsigned int i = 0; i < list.size(); ++i ) {
208 
209  config = pool -> CommonList( it -> first, list[i] ) ;
210 
211  if ( ! config ) {
212  LOG("Algorithm", pFATAL)
213  << "No Commom parameters available for " << it -> first << " list "
214  << list[i] << " at the ConfigPool";
215 
216  exit( 78 ) ;
217  }
218  else {
219  AddLowRegistry( config, false ) ;
220  LOG("Algorithm", pDEBUG) << "Loading "
221  << it -> first << " registry "
222  << list[i] << " \n" << config;
223  }
224 
225  }
226 
227  }
228 
229 
230  // Load Tunable from CommonParameters
231  // only if the option is specified in RunOpt
232  config = pool -> CommonList( "Param", "Tunable" ) ;
233  if ( config ) {
234  if ( config -> NEntries() > 0 ) {
235  AddTopRegistry( config, false ) ;
236  LOG("Algorithm", pDEBUG) << "Loading Tunable registry \n" << config;
237  }
238  }
239  else {
240  // notify & keep whatever config Registry was used before.
241  LOG("Algorithm", pWARN)
242  << "No Tunable parameter set available at the ConfigPool";
243  }
244 
245  if ( fConfig ) {
246  delete fConfig ;
247  fConfig = 0 ;
248  }
249 
250 }
set< int >::iterator it
A singleton class holding all configuration registries built while parsing all loaded XML configurati...
Definition: AlgConfigPool.h:41
#define pFATAL
Definition: Messenger.h:57
Definition: config.py:1
AlgId fID
algorithm name and configuration set
Definition: Algorithm.h:156
map< RgKey, RegistryItemI * >::const_iterator RgIMapConstIter
Definition: Registry.h:50
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const RgIMap & GetItemMap(void) const
Definition: Registry.h:162
string Name(void) const
Definition: AlgId.h:45
bool ItemIsLocal(RgKey key) const
local or global?
Definition: Registry.cxx:193
int AddTopRegistry(Registry *rp, bool owns=true)
add registry with top priority, also update ownership
Definition: Algorithm.cxx:585
#define pWARN
Definition: Messenger.h:61
void DeleteConfig(void)
Definition: Algorithm.cxx:471
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
RgStr GetString(RgKey key) const
Definition: Registry.cxx:496
vector< Registry * > fConfVect
Definition: Algorithm.h:161
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
exit(0)
Registry * fConfig
Summary configuration derived from fConvVect, not necessarily allocated.
Definition: Algorithm.h:194
Registry * FindRegistry(string key) const
int AddLowRegistry(Registry *rp, bool owns=true)
add registry with lowest priority, also update ownership
Definition: Algorithm.cxx:601
string Key(void) const
Definition: AlgId.h:47
static AlgConfigPool * Instance()
#define pDEBUG
Definition: Messenger.h:64
string Config(void) const
Definition: AlgId.h:46
enum BeamMode string
void MECGenerator::GenerateFermiMomentum ( GHepRecord event) const
private

Definition at line 158 of file MECGenerator.cxx.

References ana::assert(), genie::PDGLibrary::Find(), fNuclModel, genie::NuclearModelI::GenerateNucleon(), genie::PDGLibrary::Instance(), LOG, genie::utils::res::Mass(), genie::NuclearModelI::Momentum3(), NucleonClusterConstituents(), make_associated_cosmic_defs::p3, genie::GHepParticle::Pdg(), pINFO, genie::GHepParticle::SetMomentum(), and ana::Sqrt().

Referenced by ProcessEventRecord().

159 {
160 // Generate the initial state di-nucleon cluster 4-momentum.
161 // Draw Fermi momenta for each of the two nucleons.
162 // Sum the two Fermi momenta to calculate the di-nucleon momentum.
163 // For simplicity, keep the di-nucleon cluster on the mass shell.
164 //
165  GHepParticle * target_nucleus = event->TargetNucleus();
166  assert(target_nucleus);
167  GHepParticle * nucleon_cluster = event->HitNucleon();
168  assert(nucleon_cluster);
169  GHepParticle * remnant_nucleus = event->RemnantNucleus();
170  assert(remnant_nucleus);
171 
172  // generate a Fermi momentum for each nucleon
173 
174  Target tgt(target_nucleus->Pdg());
175  PDGCodeList pdgv = this->NucleonClusterConstituents(nucleon_cluster->Pdg());
176  assert(pdgv.size()==2);
177  tgt.SetHitNucPdg(pdgv[0]);
179  TVector3 p3a = fNuclModel->Momentum3();
180  tgt.SetHitNucPdg(pdgv[1]);
182  TVector3 p3b = fNuclModel->Momentum3();
183 
184  LOG("FermiMover", pINFO)
185  << "1st nucleon (code = " << pdgv[0] << ") generated momentum: ("
186  << p3a.Px() << ", " << p3a.Py() << ", " << p3a.Pz() << "), "
187  << "|p| = " << p3a.Mag();
188  LOG("FermiMover", pINFO)
189  << "2nd nucleon (code = " << pdgv[1] << ") generated momentum: ("
190  << p3b.Px() << ", " << p3b.Py() << ", " << p3b.Pz() << "), "
191  << "|p| = " << p3b.Mag();
192 
193  // calcute nucleon cluster momentum
194 
195  TVector3 p3 = p3a + p3b;
196 
197  LOG("FermiMover", pINFO)
198  << "di-nucleon cluster momentum: ("
199  << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), "
200  << "|p| = " << p3.Mag();
201 
202  // target (initial) nucleus and nucleon cluster mass
203 
204  double Mi = PDGLibrary::Instance()->Find(target_nucleus->Pdg() )-> Mass(); // initial nucleus mass
205  double M2n = PDGLibrary::Instance()->Find(nucleon_cluster->Pdg())-> Mass(); // nucleon cluster mass
206 
207  // nucleon cluster energy
208 
209  double EN = TMath::Sqrt(p3.Mag2() + M2n*M2n);
210 
211  // set the remnant nucleus and nucleon cluster 4-momenta at GHEP record
212 
213  TLorentzVector p4nclust ( p3.Px(), p3.Py(), p3.Pz(), EN );
214  TLorentzVector p4remnant (-1*p3.Px(), -1*p3.Py(), -1*p3.Pz(), Mi-EN);
215 
216  nucleon_cluster->SetMomentum(p4nclust);
217  remnant_nucleus->SetMomentum(p4remnant);
218 
219  // set the nucleon cluster 4-momentum at the interaction summary
220 
221  event->Summary()->InitStatePtr()->TgtPtr()->SetHitNucP4(p4nclust);
222 }
const NuclearModelI * fNuclModel
Definition: MECGenerator.h:66
double Mass(Resonance_t res)
resonance mass (GeV)
void SetMomentum(const TLorentzVector &p4)
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:67
A list of PDG codes.
Definition: PDGCodeList.h:33
int Pdg(void) const
Definition: GHepParticle.h:64
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
#define pINFO
Definition: Messenger.h:63
PDGCodeList NucleonClusterConstituents(int pdgc) const
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
assert(nhit_max >=nhit_nbins)
virtual bool GenerateNucleon(const Target &) const =0
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
void MECGenerator::GenerateNSVInitialHadrons ( GHepRecord event) const
private

Definition at line 869 of file MECGenerator.cxx.

References ana::assert(), energy, genie::PDGLibrary::Find(), fNuclModel, genie::NuclearModelI::GenerateNucleon(), genie::PDGLibrary::Instance(), genie::Interaction::KinePtr(), genie::kIStDecayedState, genie::kKineGenErr, genie::kPdgClusterNN, genie::kPdgClusterNP, genie::kPdgClusterPP, genie::controls::kRjMaxIterations, LOG, genie::utils::res::Mass(), genie::NuclearModelI::Momentum3(), NucleonClusterConstituents(), plot_validation_datamc::p1, genie::GHepParticle::P4(), pDEBUG, genie::GHepParticle::Pdg(), pERROR, pWARN, genie::NuclearModelI::RemovalEnergy(), genie::Kinematics::SetHadSystP4(), genie::GHepParticle::SetMomentum(), genie::exceptions::EVGThreadException::SetReason(), ana::Sqrt(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), and genie::GHepParticle::X4().

Referenced by ProcessEventRecord().

870 {
871  // We need a kinematic limits accept/reject loop here, so generating the
872  // initial hadrons is combined with generating the recoil hadrons...
873 
874  LOG("MEC",pDEBUG) << "Generate Initial Hadrons - Start";
875 
876  Interaction * interaction = event->Summary();
877  GHepParticle * neutrino = event->Probe();
878  assert(neutrino);
879  TLorentzVector p4nu(*neutrino->P4());
880 
881  // get final state primary lepton & its 4-momentum
882  GHepParticle * fsl = event->FinalStatePrimaryLepton();
883  assert(fsl);
884  TLorentzVector p4l(*fsl->P4());
885 
886  // calculate 4-momentum transfer from these
887  TLorentzVector Q4 = p4nu - p4l;
888 
889  // get the target two-nucleon cluster and nucleus.
890  // the remnant nucleus is apparently set, except for its momentum.
891  GHepParticle * target_nucleus = event->TargetNucleus();
892  assert(target_nucleus);
893  GHepParticle * initial_nucleon_cluster = event->HitNucleon();
894  assert(initial_nucleon_cluster);
895  GHepParticle * remnant_nucleus = event->RemnantNucleus();
896  assert(remnant_nucleus);
897 
898  // -- make a two-nucleon system, then give them some momenta.
899 
900  // instantiate an empty local target nucleus, so I can use existing methods
901  // to get a momentum from the prevailing Fermi-motion distribution
902  Target tgt(target_nucleus->Pdg());
903 
904  // NucleonClusterConstituents is an implementation within this class, called with this
905  // It using the nucleon cluster from the earlier tests for a pn state,
906  // the method returns a vector of pdgs, which hopefully will be of size two.
907 
908  PDGCodeList pdgv = this->NucleonClusterConstituents(initial_nucleon_cluster->Pdg());
909  assert(pdgv.size()==2);
910 
911  // These things need to be saved through to the end of the accept loop.
912  bool accept = false;
913  TVector3 p31i;
914  TVector3 p32i;
915  unsigned int iter = 0;
916 
917  int initial_nucleon_cluster_pdg = initial_nucleon_cluster->Pdg();
918  int final_nucleon_cluster_pdg = 0;
919  if (neutrino->Pdg() > 0) {
920  if (initial_nucleon_cluster->Pdg() == kPdgClusterNP) {
921  final_nucleon_cluster_pdg = kPdgClusterPP;
922  }
923  else if (initial_nucleon_cluster->Pdg() == kPdgClusterNN) {
924  final_nucleon_cluster_pdg = kPdgClusterNP;
925  }
926  else {
927  LOG("MEC", pERROR) << "Wrong pdg for a CC neutrino MEC interaction"
928  << initial_nucleon_cluster->Pdg();
929  }
930  }
931  else if (neutrino->Pdg() < 0) {
932  if (initial_nucleon_cluster->Pdg() == kPdgClusterNP) {
933  final_nucleon_cluster_pdg = kPdgClusterNN;
934  }
935  else if (initial_nucleon_cluster->Pdg() == kPdgClusterPP) {
936  final_nucleon_cluster_pdg = kPdgClusterNP;
937  }
938  else {
939  LOG("MEC", pERROR) << "Wrong pdg for a CC anti-neutrino MEC interaction"
940  << initial_nucleon_cluster->Pdg();
941  }
942  }
943 
944  TLorentzVector p4initial_cluster;
945  TLorentzVector p4final_cluster;
946  TLorentzVector p4remnant_nucleus;
947  double removalenergy1;
948  double removalenergy2;
949 
950  //===========================================================================
951  // Choose two nucleons from the prevailing fermi-motion distribution.
952  // Some produce kinematically unallowed combinations initial cluster and Q2
953  // Find out, and if so choose them again with this accept/reject loop.
954  // Some kinematics are especially tough
955  while(!accept){
956  iter++;
957  if(iter > kRjMaxIterations*1000) {
958  // error if try too many times
959  LOG("MEC", pWARN)
960  << "Couldn't select a valid W, Q^2 pair after "
961  << iter << " iterations";
962  event->EventFlags()->SetBitNumber(kKineGenErr, true);
964  exception.SetReason("Couldn't select initial hadron kinematics");
965  exception.SwitchOnFastForward();
966  throw exception;
967  }
968 
969  // generate nucleons
970  // tgt is a Target object for local use, just waiting to be filled.
971  // this sets the pdg of each nucleon and its momentum from a Fermi-gas
972  // Nieves et al. would use a local Fermi gas here, not this, but ok.
973  // so momentum from global Fermi gas, local Fermi gas, or spectral function
974  // and removal energy ~0.025 GeV, correlated with density, or from SF distribution
975  tgt.SetHitNucPdg(pdgv[0]);
977  p31i = fNuclModel->Momentum3();
978  removalenergy1 = fNuclModel->RemovalEnergy();
979  tgt.SetHitNucPdg(pdgv[1]);
981  p32i = fNuclModel->Momentum3();
982  removalenergy2 = fNuclModel->RemovalEnergy();
983 
984  // not sure -- could give option to use Nieves q-value here.
985 
986  // Now write down the initial cluster four-vector for this choice
987  TVector3 p3i = p31i + p32i;
988  double mass2 = PDGLibrary::Instance()->Find( initial_nucleon_cluster_pdg )->Mass();
989  mass2 *= mass2;
990  double energy = TMath::Sqrt(p3i.Mag2() + mass2);
991  p4initial_cluster.SetPxPyPzE(p3i.Px(),p3i.Py(),p3i.Pz(),energy);
992 
993  // cast the removal energy as the energy component of a 4-vector for later.
994  TLorentzVector tLVebind(0., 0., 0., -1.0 * (removalenergy1 + removalenergy2));
995 
996  // RIK: You might ask why is this the right place to subtract ebind?
997  // It is okay. Physically, I'm subtracting it from q. The energy
998  // transfer to the nucleon is 50 MeV less. the energy cost to put this
999  // cluster on-shell. What Jan says he does in PRC.86.015504 is this:
1000  // The nucleons are assumed to be in a potential well
1001  // V = Efermi + 8 MeV.
1002  // The Fermi energy is subtracted from each initial-state nucleon
1003  // (I guess he does it dynamically based on Ef = p2/2M or so which
1004  // is what we are doing above, on average). Then after the reaction,
1005  // another 8 MeV is subtracted at that point; a small adjustment to the
1006  // momentum is needed to keep the nucleon on shell.
1007 
1008  // calculate recoil nucleon cluster 4-momentum (tLVebind is negative)
1009  p4final_cluster = p4initial_cluster + Q4 + tLVebind;
1010 
1011  // Test if the resulting four-vector corresponds to a high-enough invariant mass.
1012  // Fail the accept if we couldn't put this thing on-shell.
1013  if (p4final_cluster.M() <
1014  PDGLibrary::Instance()->Find(final_nucleon_cluster_pdg )->Mass()) {
1015  accept = false;
1016  } else {
1017  accept = true;
1018  }
1019 
1020  } // end accept loop
1021 
1022  // we got here if we accepted the final state two-nucleon system
1023  // so now we need to write everything to ghep
1024 
1025  // First the initial state nucleons.
1026  initial_nucleon_cluster->SetMomentum(p4initial_cluster);
1027 
1028  // and the remnant nucleus
1029  double Mi = PDGLibrary::Instance()->Find(target_nucleus->Pdg() )-> Mass();
1030  remnant_nucleus->SetMomentum(-1.0*p4initial_cluster.Px(),
1031  -1.0*p4initial_cluster.Py(),
1032  -1.0*p4initial_cluster.Pz(),
1033  Mi - p4initial_cluster.E() + removalenergy1 + removalenergy2);
1034 
1035  // Now the final nucleon cluster.
1036 
1037  // Getting this v4 is a position, i.e. a position within the nucleus (?)
1038  // possibly it takes on a non-trivial value only for local Fermi gas
1039  // or for sophisticated treatments of intranuclear rescattering.
1040  TLorentzVector v4(*neutrino->X4());
1041 
1042  // Now write the (undecayed) final two-nucleon system
1043  GHepParticle p1(final_nucleon_cluster_pdg, kIStDecayedState,
1044  2, -1, -1, -1, p4final_cluster, v4);
1045 
1046  //p1.SetRemovalEnergy(removalenergy1 + removalenergy2);
1047  // The "bound particle" concept applies only to p or n.
1048  // Instead, add this directly to the remnant nucleon a few lines above.
1049 
1050  // actually, this is not an status1 particle, so it is not picked up
1051  // by the aggregator. anyway, the aggregator does not run until the very end.
1052 
1053  event->AddParticle(p1);
1054 
1055  interaction->KinePtr()->SetHadSystP4(p4final_cluster);
1056 }
double RemovalEnergy(void) const
Definition: NuclearModelI.h:57
const NuclearModelI * fNuclModel
Definition: MECGenerator.h:66
#define pERROR
Definition: Messenger.h:60
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
const int kPdgClusterNP
Definition: PDGCodes.h:192
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
double Mass(Resonance_t res)
resonance mass (GeV)
void SetMomentum(const TLorentzVector &p4)
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:67
A list of PDG codes.
Definition: PDGCodeList.h:33
const int kPdgClusterNN
Definition: PDGCodes.h:191
int Pdg(void) const
Definition: GHepParticle.h:64
Summary information for an interaction.
Definition: Interaction.h:56
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
double energy
Definition: plottest35.C:25
PDGCodeList NucleonClusterConstituents(int pdgc) const
#define pWARN
Definition: Messenger.h:61
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
assert(nhit_max >=nhit_nbins)
void SetHadSystP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:317
virtual bool GenerateNucleon(const Target &) const =0
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
#define pDEBUG
Definition: Messenger.h:64
const int kPdgClusterPP
Definition: PDGCodes.h:193
const Registry & Algorithm::GetConfig ( void  ) const
virtualinherited

Get configuration registry Evaluate the summary of the configuration and returns it The summary of a configuration is a merge of all the registries known to the algorithm (see Configure() methods) but every parameter is appearing only once and in case of repetitions, only the parameter from the registry with the highest prioriry is considered.

Definition at line 254 of file Algorithm.cxx.

References febshutoff_auto::end, genie::Algorithm::GetConfig(), MECModelEnuComparisons::i, LOG, pDEBUG, r(), and moon_position_table_new3::second.

Referenced by genie::EventGeneratorListAssembler::AssembleGeneratorList(), GetAlgorithms(), genie::Algorithm::GetConfig(), genie::GRV98LO::GRV98LO(), genie::NewQELXSec::Integrate(), genie::LHAPDF5::LHAPDF5(), genie::IBDXSecMap::LoadConfig(), genie::Decayer::LoadConfig(), genie::PythiaHadronization::LoadConfig(), genie::FGMBodekRitchie::LoadConfig(), genie::NuclearModelMap::LoadConfig(), genie::SmithMonizUtils::LoadConfig(), main(), genie::AlgFactory::Print(), TestPythiaTauDecays(), testReconfigInOwnedModules(), and genie::P33PaschosLalakulichPXSec::XSec().

254  {
255 
256  if ( fConfig ) return * fConfig ;
257 
258  const_cast<Algorithm*>( this ) -> fConfig = new Registry( fID.Key() + "_summary", false ) ;
259 
260  // loop and append
261  // understand the append mechanism
262  for ( unsigned int i = 0 ; i < fConfVect.size(); ++i ) {
263  fConfig -> Append( * fConfVect[i] ) ;
264  }
265 
266  if ( fOwnsSubstruc ) {
267 
268  for ( AlgMapConstIter iter = fOwnedSubAlgMp -> begin() ;
269  iter != fOwnedSubAlgMp -> end() ; ++iter ) {
270 
271  Algorithm * subalg = iter -> second ;
272 
273  LOG("Algorithm", pDEBUG) << "Appending config from " << iter -> first << " -> " << subalg -> Id() ;
274  const Registry & r = subalg->GetConfig();
275  RgKey prefix = iter -> first + "/";
276  fConfig -> Append(r,prefix);
277 
278  }
279 
280  } //if owned substructure
281 
282  return * fConfig ;
283 }
AlgMap * fOwnedSubAlgMp
local pool for owned sub-algs (taken out of the factory pool)
Definition: Algorithm.h:167
bool fOwnsSubstruc
true if it owns its substructure (sub-algs,...)
Definition: Algorithm.h:155
Algorithm abstract base class.
Definition: Algorithm.h:54
AlgId fID
algorithm name and configuration set
Definition: Algorithm.h:156
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const XML_Char * prefix
Definition: expat.h:380
map< string, Algorithm * >::const_iterator AlgMapConstIter
Definition: Algorithm.h:51
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
vector< Registry * > fConfVect
Definition: Algorithm.h:161
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
TRandom3 r(0)
Registry * fConfig
Summary configuration derived from fConvVect, not necessarily allocated.
Definition: Algorithm.h:194
string Key(void) const
Definition: AlgId.h:47
#define pDEBUG
Definition: Messenger.h:64
Registry * Algorithm::GetOwnedConfig ( void  )
inherited

Returns the pointer of the summary registry, see previous method Gives access to the summary so it could be changed. The usage of this method is deprecated as it is mantained only for back compatibility. If you need to add or chage a parter (or more), use the AddTopRegistry() instead

Definition at line 287 of file Algorithm.cxx.

References GetConfig().

Referenced by genie::TransverseEnhancementFFModel::LoadConfig(), and genie::EffectiveSF::LoadConfig().

288 {
289 
290  GetConfig() ;
291  return fConfig;
292 }
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
Registry * fConfig
Summary configuration derived from fConvVect, not necessarily allocated.
Definition: Algorithm.h:194
template<class T >
bool genie::Algorithm::GetParam ( const RgKey name,
T p,
bool  is_top_call = true 
) const
protectedinherited

Ideal access to a parameter value from the vector of registries Returns true if the value is found and the parameters is set

Referenced by genie::CollinsSpillerFragm::BuildFunction(), genie::PetersonFragm::BuildFunction(), genie::COHXSec::LoadConfig(), genie::DISXSec::LoadConfig(), genie::INukeDeltaPropg::LoadConfig(), genie::HadronTransporter::LoadConfig(), genie::DFRKinematicsGenerator::LoadConfig(), genie::RSHelicityAmplModelNCn::LoadConfig(), genie::RSHelicityAmplModelNCp::LoadConfig(), genie::BaryonResonanceDecayer::LoadConfig(), genie::RESKinematicsGenerator::LoadConfig(), genie::DMDISXSec::LoadConfig(), genie::DipoleAxialFormFactorModel::LoadConfig(), genie::DipoleELFormFactorsModel::LoadConfig(), genie::COHKinematicsGenerator::LoadConfig(), genie::VertexGenerator::LoadConfig(), genie::H3AMNuGammaPXSec::LoadConfig(), genie::RSPPResonanceSelector::LoadConfig(), genie::IBDXSecMap::LoadConfig(), genie::Decayer::LoadConfig(), genie::DISHadronicSystemGenerator::LoadConfig(), genie::EmpiricalMECPXSec2015::LoadConfig(), genie::COHElasticPXSec::LoadConfig(), genie::KuzminNaumov2016AxialFormFactorModel::LoadConfig(), genie::SlowRsclCharmDISPXSecLO::LoadConfig(), genie::UnstableParticleDecayer::LoadConfig(), genie::AhrensNCELPXSec::LoadConfig(), genie::AlamSimoAtharVacasSKPXSec2014::LoadConfig(), genie::PythiaHadronization::LoadConfig(), genie::ReinDFRPXSec::LoadConfig(), genie::StrumiaVissaniIBDPXSec::LoadConfig(), LoadConfig(), genie::BYPDF::LoadConfig(), genie::QPMDISPXSec::LoadConfig(), genie::DFRXSec::LoadConfig(), genie::RosenbluthPXSec::LoadConfig(), genie::KNOPythiaHadronization::LoadConfig(), genie::P33PaschosLalakulichPXSec::LoadConfig(), genie::AhrensDMELPXSec::LoadConfig(), genie::MECXSec::LoadConfig(), genie::AivazisCharmPXSecLO::LoadConfig(), genie::BergerSehgalFMCOHPiPXSec2015::LoadConfig(), genie::ZExpAxialFormFactorModel::LoadConfig(), genie::QPMDMDISPXSec::LoadConfig(), genie::BBA05ELFormFactorsModel::LoadConfig(), genie::BergerSehgalCOHPiPXSec2015::LoadConfig(), genie::BBA03ELFormFactorsModel::LoadConfig(), genie::LwlynSmithQELCCPXSec::LoadConfig(), genie::ReinSehgalRESXSec::LoadConfig(), genie::NuElectronPXSec::LoadConfig(), genie::PrimaryLeptonGenerator::LoadConfig(), genie::PaisQELLambdaPXSec::LoadConfig(), genie::FGMBodekRitchie::LoadConfig(), genie::SpectralFunc1d::LoadConfig(), genie::ReinSehgalCOHPiPXSec::LoadConfig(), genie::OutgoingDarkGenerator::LoadConfig(), genie::LHAPDF6::LoadConfig(), genie::NievesSimoVacasMECPXSec2016::LoadConfig(), genie::CharmHadronization::LoadConfig(), genie::ReinSehgalRESXSecFast::LoadConfig(), genie::EventGenerator::LoadConfig(), genie::NuclearModelMap::LoadConfig(), genie::ReinSehgalSPPXSec::LoadConfig(), genie::ReinSehgalRESPXSec::LoadConfig(), genie::LwlynSmithFF::LoadConfig(), genie::SmithMonizQELCCPXSec::LoadConfig(), genie::QPMDISStrucFuncBase::LoadConfig(), genie::BBA07ELFormFactorsModel::LoadConfig(), genie::HAIntranuke::LoadConfig(), genie::NievesQELCCPXSec::LoadConfig(), genie::HAIntranuke2018::LoadConfig(), genie::HNIntranuke2018::LoadConfig(), genie::LocalFGM::LoadConfig(), genie::BSKLNBaseRESPXSec2014::LoadConfig(), genie::EffectiveSF::LoadConfig(), genie::ReinSehgalSPPPXSec::LoadConfig(), genie::KNOHadronization::LoadConfig(), genie::SmithMonizUtils::LoadConfig(), genie::MECInteractionListGenerator::LoadConfigData(), genie::PhysInteractionSelector::LoadConfigData(), genie::RESInteractionListGenerator::LoadConfigData(), genie::PauliBlocker::LoadModelType(), genie::BYStrucFunc::ReadBYParams(), and genie::LHAPDF5::SetPDFSetFromConfig().

template<class T >
bool genie::Algorithm::GetParamDef ( const RgKey name,
T p,
const T def 
) const
protectedinherited

Ideal access to a parameter value from the vector of registries, With default value. Returns true if the value is set from the registries, false if the value is the default

Referenced by genie::IMDXSec::LoadConfig(), genie::COHXSec::LoadConfig(), genie::RESXSec::LoadConfig(), genie::DISXSec::LoadConfig(), genie::DFRKinematicsGenerator::LoadConfig(), genie::COHXSecAR::LoadConfig(), genie::DMDISXSec::LoadConfig(), genie::BaryonResonanceDecayer::LoadConfig(), genie::SKKinematicsGenerator::LoadConfig(), genie::COHElKinematicsGenerator::LoadConfig(), genie::NuEKinematicsGenerator::LoadConfig(), genie::QELXSec::LoadConfig(), genie::RESKinematicsGenerator::LoadConfig(), genie::COHKinematicsGenerator::LoadConfig(), genie::IBDKinematicsGenerator::LoadConfig(), genie::NuEInteractionListGenerator::LoadConfig(), genie::QELKinematicsGenerator::LoadConfig(), genie::DMELXSec::LoadConfig(), genie::DISHadronicSystemGenerator::LoadConfig(), genie::DISKinematicsGenerator::LoadConfig(), genie::NucBindEnergyAggregator::LoadConfig(), genie::DMELKinematicsGenerator::LoadConfig(), genie::DMDISKinematicsGenerator::LoadConfig(), genie::QPMDISPXSec::LoadConfig(), genie::DFRXSec::LoadConfig(), genie::AhrensDMELPXSec::LoadConfig(), genie::NuElectronXSec::LoadConfig(), genie::QELEventGenerator::LoadConfig(), genie::P33PaschosLalakulichPXSec::LoadConfig(), genie::MECXSec::LoadConfig(), genie::FermiMover::LoadConfig(), genie::AlamSimoAtharVacasSKXSec::LoadConfig(), genie::QPMDMDISPXSec::LoadConfig(), genie::LwlynSmithQELCCPXSec::LoadConfig(), genie::ReinSehgalRESXSec::LoadConfig(), genie::FGMBodekRitchie::LoadConfig(), genie::ReinSehgalRESXSecFast::LoadConfig(), genie::KovalenkoQELCharmPXSec::LoadConfig(), genie::SmithMonizQELCCXSec::LoadConfig(), genie::ReinSehgalSPPXSec::LoadConfig(), genie::ReinSehgalRESPXSec::LoadConfig(), genie::QELEventGeneratorSM::LoadConfig(), genie::QPMDISStrucFuncBase::LoadConfig(), genie::SmithMonizQELCCPXSec::LoadConfig(), genie::NievesQELCCPXSec::LoadConfig(), genie::HAIntranuke::LoadConfig(), genie::LocalFGM::LoadConfig(), genie::HNIntranuke2018::LoadConfig(), genie::HAIntranuke2018::LoadConfig(), genie::BSKLNBaseRESPXSec2014::LoadConfig(), genie::EffectiveSF::LoadConfig(), genie::KNOHadronization::LoadConfig(), genie::NewQELXSec::LoadConfig(), genie::QELInteractionListGenerator::LoadConfigData(), genie::MECInteractionListGenerator::LoadConfigData(), genie::DFRInteractionListGenerator::LoadConfigData(), genie::COHInteractionListGenerator::LoadConfigData(), genie::RESInteractionListGenerator::LoadConfigData(), genie::SKInteractionListGenerator::LoadConfigData(), genie::DMELInteractionListGenerator::LoadConfigData(), genie::RSPPInteractionListGenerator::LoadConfigData(), genie::DISInteractionListGenerator::LoadConfigData(), and genie::DMDISInteractionListGenerator::LoadConfigData().

template<class T >
bool genie::Algorithm::GetParamVect ( const std::string comm_name,
std::vector< T > &  v,
unsigned int  max,
bool  is_top_call = true 
) const
protectedinherited

Handle to load vectors of parameters It looks for different registry item with name comm_name0, comm_name1, etc...

virtual AlgStatus_t genie::Algorithm::GetStatus ( void  ) const
inlinevirtualinherited

Get algorithm status.

Definition at line 101 of file Algorithm.h.

References genie::Algorithm::fStatus.

101 { return fStatus; }
AlgStatus_t fStatus
algorithm execution status
Definition: Algorithm.h:166
virtual const AlgId& genie::Algorithm::Id ( void  ) const
inlinevirtualinherited

Get algorithm ID.

Definition at line 98 of file Algorithm.h.

References genie::Algorithm::fID.

Referenced by genie::KineGeneratorWithCache::AccessCacheBranch(), genie::QELEventGeneratorSM::AccessCacheBranch2(), genie::QELEventGeneratorSM::AccessCacheBranchDiffv(), genie::InteractionListAssembler::AssembleInteractionList(), genie::XSecAlgorithmMap::BuildMap(), genie::InteractionGeneratorMap::BuildMap(), genie::XSecSplineList::BuildSplineKey(), genie::DISXSec::CacheBranchName(), genie::ReinSehgalRESXSecWithCache::CacheBranchName(), genie::DMDISXSec::CacheBranchName(), genie::ReinSehgalRESXSecWithCacheFast::CacheBranchName(), genie::Algorithm::Compare(), genie::RESKinematicsGenerator::ComputeMaxXSec(), genie::COHElKinematicsGenerator::ComputeMaxXSec(), genie::SKKinematicsGenerator::ComputeMaxXSec(), genie::COHKinematicsGenerator::ComputeMaxXSec(), genie::Algorithm::Configure(), genie::GEVGDriver::CreateSplines(), genie::QPMDISPXSec::DISRESJoinSuppressionFactor(), genie::QPMDMDISPXSec::DMDISRESJoinSuppressionFactor(), genie::AlgConfigPool::FindRegistry(), genie::AlgFactory::ForceReconfiguration(), genie::GEVGDriver::GenerateEvent(), GetAlgorithms(), genie::LwlynSmithQELCCPXSec::Integral(), genie::NievesQELCCPXSec::Integral(), genie::COHXSec::Integrate(), genie::QPMDISPXSec::LoadConfig(), genie::QPMDMDISPXSec::LoadConfig(), genie::EventGenerator::LoadConfig(), genie::EventGeneratorListAssembler::LoadGenerator(), main(), genie::COHKinematicsGenerator::MaxXSec_AlvarezRuso(), genie::XSecAlgorithmMap::Print(), genie::InteractionGeneratorMap::Print(), genie::AlgFactory::Print(), genie::COHHadronicSystemGenerator::ProcessEventRecord(), genie::COHPrimaryLeptonGenerator::ProcessEventRecord(), genie::COHKinematicsGenerator::ProcessEventRecord(), ProcessEventRecord(), genie::EventGenerator::ProcessEventRecord(), genie::KNOPythiaHadronization::SelectHadronizer(), TestPythiaTauDecays(), and genie::GEVGDriver::UseSplines().

98 { return fID; }
AlgId fID
algorithm name and configuration set
Definition: Algorithm.h:156
void Algorithm::Initialize ( void  )
protectedinherited

Definition at line 343 of file Algorithm.cxx.

Referenced by genie::Algorithm::AllowReconfig().

344 {
345 // Algorithm initialization
346 //
347  fAllowReconfig = true;
348  fOwnsSubstruc = false;
349  fConfig = 0;
350  fOwnedSubAlgMp = 0;
351 }
AlgMap * fOwnedSubAlgMp
local pool for owned sub-algs (taken out of the factory pool)
Definition: Algorithm.h:167
bool fOwnsSubstruc
true if it owns its substructure (sub-algs,...)
Definition: Algorithm.h:155
Registry * fConfig
Summary configuration derived from fConvVect, not necessarily allocated.
Definition: Algorithm.h:194
void MECGenerator::LoadConfig ( void  )
private

Definition at line 1070 of file MECGenerator.cxx.

References ana::assert(), fNuclModel, fQ3Max, genie::Algorithm::GetParam(), and genie::Algorithm::SubAlg().

Referenced by Configure().

1071 {
1072  fNuclModel = 0;
1073  RgKey nuclkey = "NuclearModel";
1074  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
1075  assert(fNuclModel);
1076 
1077  GetParam( "NSV-Q3Max", fQ3Max ) ;
1078 }
const NuclearModelI * fNuclModel
Definition: MECGenerator.h:66
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
string RgKey
assert(nhit_max >=nhit_nbins)
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353
int Algorithm::MergeTopRegistry ( const Registry r)
protectedinherited

Merge with top level registry if first reg of the vector is owned Otherwise an owned copy is added as a top registry

Definition at line 618 of file Algorithm.cxx.

618  {
619 
620  if ( fOwnerships.empty() ) {
621 
622  // this algorithm is not configured right now, the incoming registry is the only configuration
623  Registry * p = new Registry( r ) ;
624  AddTopRegistry( p ) ;
625 
626  return 1 ;
627  }
628 
629  if ( fOwnerships[0] ) {
630  //the top registry is owned: it can be changed with no consequences for other algorithms
631  fConfVect[0] -> Merge( r ) ;
632  }
633  else {
634  // The top registry is not owned so it cannot be changed
635  // The registry will be added with top priority
636 
637  Registry * p = new Registry( r ) ;
638  AddTopRegistry( p ) ;
639  }
640 
641  // The configuration has changed so the summary is not updated anymore and must be deleted
642  if ( fConfig ) {
643  delete fConfig ;
644  fConfig = 0 ;
645  }
646 
647  return fConfVect.size() ;
648 }
const char * p
Definition: xmltok.h:285
int AddTopRegistry(Registry *rp, bool owns=true)
add registry with top priority, also update ownership
Definition: Algorithm.cxx:585
vector< Registry * > fConfVect
Definition: Algorithm.h:161
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
vector< bool > fOwnerships
ownership for every registry in fConfVect
Definition: Algorithm.h:164
Registry * fConfig
Summary configuration derived from fConvVect, not necessarily allocated.
Definition: Algorithm.h:194
PDGCodeList MECGenerator::NucleonClusterConstituents ( int  pdgc) const
private

Definition at line 550 of file MECGenerator.cxx.

References genie::kPdgClusterNN, genie::kPdgClusterNP, genie::kPdgClusterPP, genie::kPdgNeutron, genie::kPdgProton, LOG, pERROR, and genie::PDGCodeList::push_back().

Referenced by DecayNucleonCluster(), GenerateFermiMomentum(), and GenerateNSVInitialHadrons().

551 {
552  bool allowdup = true;
553  PDGCodeList pdgv(allowdup);
554 
555  if(pdgc == kPdgClusterNN) {
556  pdgv.push_back(kPdgNeutron);
557  pdgv.push_back(kPdgNeutron);
558  }
559  else
560  if(pdgc == kPdgClusterNP) {
561  pdgv.push_back(kPdgNeutron);
562  pdgv.push_back(kPdgProton);
563  }
564  else
565  if(pdgc == kPdgClusterPP) {
566  pdgv.push_back(kPdgProton);
567  pdgv.push_back(kPdgProton);
568  }
569  else
570  {
571  LOG("MEC", pERROR)
572  << "Unknown di-nucleon cluster PDG code (" << pdgc << ")";
573  }
574 
575  return pdgv;
576 }
#define pERROR
Definition: Messenger.h:60
const int kPdgClusterNP
Definition: PDGCodes.h:192
A list of PDG codes.
Definition: PDGCodeList.h:33
const int kPdgClusterNN
Definition: PDGCodes.h:191
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const int kPdgProton
Definition: PDGCodes.h:65
const int kPdgNeutron
Definition: PDGCodes.h:67
const int kPdgClusterPP
Definition: PDGCodes.h:193
void Algorithm::Print ( ostream &  stream) const
virtualinherited

Print algorithm info.

Definition at line 323 of file Algorithm.cxx.

References GetConfig(), and r().

Referenced by genie::Algorithm::AllowReconfig(), and genie::operator<<().

324 {
325  // print algorithm name & parameter-set
326  stream << "\nAlgorithm Key: " << this->fID.Key();
327  stream << " - Owns Substruc: " << ((fOwnsSubstruc) ? "[true]" : "[false]");
328 
329  // print algorithm configuration
330  const Registry & r = this->GetConfig();
331  stream << r;
332 
333  if(fOwnsSubstruc) {
334  AlgMapConstIter iter = fOwnedSubAlgMp->begin();
335  for(; iter!=fOwnedSubAlgMp->end(); ++iter) {
336  Algorithm * alg = iter->second;
337  stream << "<Next algorithm is owned by : " << this->fID.Key() << ">";
338  stream << *alg;
339  }
340  }
341 }
AlgMap * fOwnedSubAlgMp
local pool for owned sub-algs (taken out of the factory pool)
Definition: Algorithm.h:167
bool fOwnsSubstruc
true if it owns its substructure (sub-algs,...)
Definition: Algorithm.h:155
Algorithm abstract base class.
Definition: Algorithm.h:54
AlgId fID
algorithm name and configuration set
Definition: Algorithm.h:156
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
map< string, Algorithm * >::const_iterator AlgMapConstIter
Definition: Algorithm.h:51
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
TRandom3 r(0)
string Key(void) const
Definition: AlgId.h:47
void MECGenerator::ProcessEventRecord ( GHepRecord event) const
virtual

shortly, this will be handled by the InitialStateAppender module

Implements genie::EventRecordVisitorI.

Definition at line 73 of file MECGenerator.cxx.

References AddFinalStateLepton(), AddTargetRemnant(), genie::EventGeneratorI::CrossSectionAlg(), DecayNucleonCluster(), fXSecModel, GenerateFermiMomentum(), GenerateNSVInitialHadrons(), genie::Algorithm::Id(), genie::RunningThreadInfo::Instance(), LOG, genie::AlgId::Name(), pFATAL, RecoilNucleonCluster(), genie::RunningThreadInfo::RunningThread(), SelectEmpiricalKinematics(), and SelectNSVLeptonKinematics().

74 {
75  //-- Access cross section algorithm for running thread
77  const EventGeneratorI * evg = rtinfo->RunningThread();
78  fXSecModel = evg->CrossSectionAlg();
79  if (fXSecModel->Id().Name() == "genie::EmpiricalMECPXSec2015") {
80  this -> AddTargetRemnant (event); /// shortly, this will be handled by the InitialStateAppender module
81  this -> GenerateFermiMomentum(event);
82  this -> SelectEmpiricalKinematics(event);
83  // TODO: test removing `AddFinalStateLepton` and replacing it with
84  // PrimaryLeptonGenerator::ProcessEventRecord(evrec);
85  this -> AddFinalStateLepton(event);
86  // TODO: maybe put `RecoilNucleonCluster` in a `MECHadronicSystemGenerator` class,
87  // name it `GenerateEmpiricalDiNucleonCluster` or something...
88  this -> RecoilNucleonCluster(event);
89  // TODO: `DecayNucleonCluster` should probably be in `MECHadronicSystemGenerator`,
90  // if we make that...
91  this -> DecayNucleonCluster(event);
92  } else if (fXSecModel->Id().Name() == "genie::NievesSimoVacasMECPXSec2016") {
93  this -> SelectNSVLeptonKinematics(event);
94  this -> AddTargetRemnant(event);
95  this -> GenerateNSVInitialHadrons(event);
96  // Note: this method in `MECTensor/MECTensorGenerator.cxx` appeared to be a straight
97  // copy of an earlier version of the `DecayNucleonCluster` method here - but, watch
98  // for this...
99  this -> DecayNucleonCluster(event);
100  }
101  else {
102  LOG("MECGenerator",pFATAL) <<
103  "ProcessEventRecord >> Cannot calculate kinematics for " <<
104  fXSecModel->Id().Name();
105  }
106 
107 
108 }
void SelectNSVLeptonKinematics(GHepRecord *event) const
#define pFATAL
Definition: Messenger.h:57
Defines the EventGeneratorI interface.
void GenerateNSVInitialHadrons(GHepRecord *event) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
string Name(void) const
Definition: AlgId.h:45
void RecoilNucleonCluster(GHepRecord *event) const
void DecayNucleonCluster(GHepRecord *event) const
void SelectEmpiricalKinematics(GHepRecord *event) const
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
static RunningThreadInfo * Instance(void)
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
const XSecAlgorithmI * fXSecModel
Definition: MECGenerator.h:64
void AddTargetRemnant(GHepRecord *event) const
void GenerateFermiMomentum(GHepRecord *event) const
const EventGeneratorI * RunningThread(void)
Keep info on the event generation thread currently on charge. This is used so that event generation m...
void AddFinalStateLepton(GHepRecord *event) const
void MECGenerator::RecoilNucleonCluster ( GHepRecord event) const
private

Definition at line 386 of file MECGenerator.cxx.

References ana::assert(), genie::GHepParticle::GetP4(), genie::kIStDecayedState, LOG, genie::GHepParticle::P4(), pINFO, tmp, and genie::GHepParticle::X4().

Referenced by ProcessEventRecord().

387 {
388  // get di-nucleon cluster & its 4-momentum
389  GHepParticle * nucleon_cluster = event->HitNucleon();
390  assert(nucleon_cluster);
391  TLorentzVector * tmp=nucleon_cluster->GetP4();
392  TLorentzVector p4cluster(*tmp);
393  delete tmp;
394 
395  // get neutrino & its 4-momentum
396  GHepParticle * neutrino = event->Probe();
397  assert(neutrino);
398  TLorentzVector p4v(*neutrino->P4());
399 
400  // get final state primary lepton & its 4-momentum
401  GHepParticle * fsl = event->FinalStatePrimaryLepton();
402  assert(fsl);
403  TLorentzVector p4l(*fsl->P4());
404 
405  // calculate 4-momentum transfer
406  TLorentzVector q = p4v - p4l;
407 
408  // calculate recoil nucleon cluster 4-momentum
409  TLorentzVector p4cluster_recoil = p4cluster + q;
410 
411  // work-out recoil nucleon cluster code
412  LOG("MEC", pINFO) << "Interaction summary";
413  LOG("MEC", pINFO) << *event->Summary();
414  int recoil_nucleon_cluster_pdg = event->Summary()->RecoilNucleonPdg();
415 
416  // vtx position in nucleus coord system
417  TLorentzVector v4(*neutrino->X4());
418 
419  // add to the event record
420  event->AddParticle(
421  recoil_nucleon_cluster_pdg, kIStDecayedState,
422  2, -1, -1, -1, p4cluster_recoil, v4);
423 }
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
Float_t tmp
Definition: plot.C:36
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
TLorentzVector * GetP4(void) const
#define pINFO
Definition: Messenger.h:63
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
assert(nhit_max >=nhit_nbins)
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
void MECGenerator::SelectEmpiricalKinematics ( GHepRecord event) const
private

Definition at line 224 of file MECGenerator.cxx.

References genie::Kinematics::ClearRunningValues(), genie::EventGeneratorI::CrossSectionAlg(), genie::PDGLibrary::Find(), fXSecModel, gQ2, genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::RunningThreadInfo::Instance(), genie::PDGLibrary::Instance(), genie::RandomGen::Instance(), genie::ProcessInfo::IsEM(), genie::utils::mec::J(), genie::Interaction::KinePtr(), genie::kKineGenErr, genie::kPSWQ2fE, genie::kRfHitNucRest, genie::controls::kRjMaxIterations, LOG, pINFO, pNOTICE, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), pWARN, genie::utils::kinematics::Q2(), generate_hists::rnd, genie::RandomGen::RndKine(), genie::RunningThreadInfo::RunningThread(), genie::Kinematics::SetQ2(), genie::exceptions::EVGThreadException::SetReason(), genie::Kinematics::SetW(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), confusionMatrixTree::t, genie::InitialState::Tgt(), W, Wmax, Wmin, genie::utils::kinematics::WQ2toXY(), genie::XSecAlgorithmI::XSec(), and xsec.

Referenced by ProcessEventRecord().

225 {
226 // Select interaction kinematics using the rejection method.
227 //
228 
229  // Access cross section algorithm for running thread
231  const EventGeneratorI * evg = rtinfo->RunningThread();
232  fXSecModel = evg->CrossSectionAlg();
233 
234  Interaction * interaction = event->Summary();
235  double Ev = interaction->InitState().ProbeE(kRfHitNucRest);
236 
237  // **** NOTE / TODO:
238  // **** Hardcode bogus limits for the time-being
239  // **** Should be able to get limits via Interaction::KPhaseSpace
240  double Q2min = 0.01;
241  double Q2max = 8.00;
242  double Wmin = 1.88;
243  double Wmax = 3.00;
244 
245  // Scan phase-space for the maximum differential cross section
246  // at the current neutrino energy
247  const int nq=30;
248  const int nw=20;
249  double dQ2 = (Q2max-Q2min) / (nq-1);
250  double dW = (Wmax-Wmin ) / (nw-1);
251  double xsec_max = 0;
252  for(int iw=0; iw<nw; iw++) {
253  for(int iq=0; iq<nq; iq++) {
254  double Q2 = Q2min + iq*dQ2;
255  double W = Wmin + iw*dW;
256  interaction->KinePtr()->SetQ2(Q2);
257  interaction->KinePtr()->SetW (W);
258  double xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
259  xsec_max = TMath::Max(xsec, xsec_max);
260  }
261  }
262  LOG("MEC", pNOTICE) << "xsec_max (E = " << Ev << " GeV) = " << xsec_max;
263 
264  // Select kinematics
266  unsigned int iter = 0;
267  bool accept = false;
268  while(1) {
269  iter++;
270  if(iter > kRjMaxIterations) {
271  LOG("MEC", pWARN)
272  << "Couldn't select a valid W, Q^2 pair after "
273  << iter << " iterations";
274  event->EventFlags()->SetBitNumber(kKineGenErr, true);
276  exception.SetReason("Couldn't select kinematics");
277  exception.SwitchOnFastForward();
278  throw exception;
279  }
280 
281  // Generate next pair
282  double gQ2 = Q2min + (Q2max-Q2min) * rnd->RndKine().Rndm();
283  double gW = Wmin + (Wmax -Wmin ) * rnd->RndKine().Rndm();
284 
285  // Calculate d2sigma/dQ2dW
286  interaction->KinePtr()->SetQ2(gQ2);
287  interaction->KinePtr()->SetW (gW);
288  double xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
289 
290  // Decide whether to accept the current kinematics
291  double t = xsec_max * rnd->RndKine().Rndm();
292  double J = 1; // jacobean
293  accept = (t < J*xsec);
294 
295  // If the generated kinematics are accepted, finish-up module's job
296  if(accept) {
297  LOG("MEC", pINFO) << "Selected: Q^2 = " << gQ2 << ", W = " << gW;
298  double gx = 0;
299  double gy = 0;
300  // More accurate calculation of the mass of the cluster than 2*Mnucl
301  int nucleon_cluster_pdg = interaction->InitState().Tgt().HitNucPdg();
302  double M2n = PDGLibrary::Instance()->Find(nucleon_cluster_pdg)->Mass();
303  bool is_em = interaction->ProcInfo().IsEM();
304  kinematics::WQ2toXY(Ev,M2n,gW,gQ2,gx,gy);
305 
306  LOG("MEC", pINFO) << "x = " << gx << ", y = " << gy;
307  // lock selected kinematics & clear running values
308  interaction->KinePtr()->SetQ2(gQ2, true);
309  interaction->KinePtr()->SetW (gW, true);
310  interaction->KinePtr()->Setx (gx, true);
311  interaction->KinePtr()->Sety (gy, true);
312  interaction->KinePtr()->ClearRunningValues();
313 
314  return;
315  }//accepted?
316  }//iter
317 }
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
int HitNucPdg(void) const
Definition: Target.cxx:321
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:265
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
Defines the EventGeneratorI interface.
double Wmax[kNWBins]
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
Summary information for an interaction.
Definition: Interaction.h:56
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1046
#define pINFO
Definition: Messenger.h:63
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:51
Double_t xsec[nknots]
Definition: testXsec.C:47
#define pWARN
Definition: Messenger.h:61
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
bool IsEM(void) const
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:241
static RunningThreadInfo * Instance(void)
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:253
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
const XSecAlgorithmI * fXSecModel
Definition: MECGenerator.h:64
double gQ2
Definition: gtestDISSF.cxx:55
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double Wmin[kNWBins]
#define pNOTICE
Definition: Messenger.h:62
void ClearRunningValues(void)
Definition: Kinematics.cxx:357
const Target & Tgt(void) const
Definition: InitialState.h:67
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
const EventGeneratorI * RunningThread(void)
#define W(x)
double ProbeE(RefFrame_t rf) const
Keep info on the event generation thread currently on charge. This is used so that event generation m...
void MECGenerator::SelectNSVLeptonKinematics ( GHepRecord event) const
private

Definition at line 578 of file MECGenerator.cxx.

References ana::assert(), genie::Interaction::ExclTagPtr(), fQ3Max, genie::Interaction::FSPrimLepton(), fXSecModel, genie::Interaction::InitState(), genie::Interaction::InitStatePtr(), genie::RandomGen::Instance(), genie::pdg::IonPdgCodeToA(), kinematics(), genie::Interaction::KinePtr(), genie::kIStNucleonTarget, genie::kIStStableFinalState, genie::kKineGenErr, genie::kKVctl, genie::kKVTl, genie::kNoResonance, genie::constants::kNucleonMass, genie::kP33_1232, genie::kPdgClusterNN, genie::kPdgClusterNP, genie::kPdgClusterPP, genie::kPdgTgtC12, genie::kPdgTgtFreeN, kPi, genie::kPSTlctl, genie::kRfHitNucRest, genie::controls::kRjMaxIterations, LOG, cet::sqlite::max(), pDEBUG, pERROR, pINFO, genie::GHepRecord::Probe(), genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), pWARN, genie::utils::kinematics::Q2(), genie::utils::kinematics::Q2YtoX(), generate_hists::rnd, genie::RandomGen::RndKine(), genie::RandomGen::RndLep(), genie::Kinematics::SetFSLeptonP4(), genie::Target::SetHitNucPdg(), genie::Kinematics::SetKV(), genie::Kinematics::SetQ2(), genie::exceptions::EVGThreadException::SetReason(), genie::XclsTag::SetResonance(), genie::Kinematics::SetW(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), ana::Sqrt(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), T, genie::InitialState::TgtPdg(), genie::InitialState::TgtPtr(), genie::GHepParticle::X4(), genie::XSecAlgorithmI::XSec(), and genie::utils::kinematics::XYtoW().

Referenced by ProcessEventRecord().

579 {
580  // -- implementation -- //
581  // The IFIC Valencia model can provide three different hadron tensors.
582  // The user probably wants all CC QE-like 2p2h events
583  // But could in principle get the no-delta component if they want (deactivated incode)
584  int FullDeltaNodelta = 1; // 1: full, 2: only delta, 3: zero delta
585 
586  // -- limit the maximum XS for the accept/reject loop -- //
587  //
588  // MaxXSec parameters. This whole calculation could be in it's own function?
589  // these need to lead to a number that is safely large enough, or crash the run.
590  double XSecMaxPar1 = 2.2504;
591  double XSecMaxPar2 = 9.41158;
592 
593  // -- Event Properties -----------------------------//
594  Interaction * interaction = event->Summary();
595  Kinematics * kinematics = interaction->KinePtr();
596 
597  double Enu = interaction->InitState().ProbeE(kRfHitNucRest);
598 
599  int NuPDG = interaction->InitState().ProbePdg();
600  int TgtPDG = interaction->InitState().TgtPdg();
601  // interacton vtx
602  TLorentzVector v4(*event->Probe()->X4());
603  TLorentzVector tempp4(0.,0.,0.,0.);
604 
605  // -- Lepton Kinematic Limits ----------------------------------------- //
606  double Costh = 0.0; // lepton angle
607  double CosthMax = 1.0;
608  double CosthMin = -1.0;
609 
610  double T = 0.0; // lepton kinetic energy
611  double TMax = std::numeric_limits<double>::max();
612  double TMin = 0.0;
613 
614  double Plep = 0.0; // lepton 3 momentum
615  double Elep = 0.0; // lepton energy
616  double LepMass = interaction->FSPrimLepton()->Mass();
617 
618  double Q0 = 0.0; // energy component of q four vector
619  double Q3 = 0.0; // magnitude of transfered 3 momentum
620  double Q2 = 0.0; // properly Q^2 (Q squared) - transfered 4 momentum.
621 
622  // Set lepton KE TMax for for throwing rndm in the accept/reject loop.
623  // We can accidentally set it too high, because the xsec will return zero.
624  // This way if someone reuses this code, they are not tripped up by it.
625  TMax = Enu - LepMass;
626 
627  // Set Tmin for throwing rndm in the accept/reject loop
628  // the hadron tensors we expect will be limited in q3
629  // therefore also the outgoing lepton KE can't be too low or costheta too backward
630  // make the accept/reject loop more efficient by using Min values.
631  if(Enu < fQ3Max){
632  TMin = 0 ;
633  CosthMin = -1 ;
634  } else {
635  TMin = TMath::Sqrt(TMath::Power(LepMass, 2) + TMath::Power((Enu - fQ3Max), 2)) - LepMass;
636  CosthMin = TMath::Sqrt(1 - TMath::Power((fQ3Max / Enu ), 2));
637  }
638 
639  // The accept/reject loop tests a rand against a maxxsec - must scale with A.
640  int NuclearA = 12;
641  double NuclearAfactorXSecMax = 1.0;
642  if (TgtPDG != kPdgTgtC12) {
643  if (TgtPDG > kPdgTgtFreeN && TgtPDG) {
644  NuclearA = pdg::IonPdgCodeToA(TgtPDG);
645  // The QE-like portion scales as A, but the Delta portion increases faster, not simple.
646  // so this gives additional safety factor. Remember, we need a safe max, not precise max.
647  if (NuclearA < 12) NuclearAfactorXSecMax *= NuclearA / 12.0;
648  else NuclearAfactorXSecMax *= TMath::Power(NuclearA/12.0, 1.4);
649  }
650  else {
651  LOG("MEC", pERROR) << "Trying to scale XSecMax for larger nuclei, but "
652  << TgtPDG << " isn't a nucleus?";
653  assert(false);
654  }
655  }
656 
657  // -- Generate and Test the Kinematics----------------------------------//
658 
660  bool accept = false;
661  unsigned int iter = 0;
662 
663  // loop over different (randomly) selected T and Costh
664  while (!accept) {
665  iter++;
666  if(iter > kRjMaxIterations) {
667  // error if try too many times
668  LOG("MEC", pWARN)
669  << "Couldn't select a valid Tmu, CosTheta pair after "
670  << iter << " iterations";
671  event->EventFlags()->SetBitNumber(kKineGenErr, true);
673  exception.SetReason("Couldn't select lepton kinematics");
674  exception.SwitchOnFastForward();
675  throw exception;
676  }
677 
678  // generate random kinetic energy T and Costh
679  T = TMin + (TMax-TMin)*rnd->RndKine().Rndm();
680  Costh = CosthMin + (CosthMax-CosthMin)*rnd->RndKine().Rndm();
681 
682  // Calculate useful values for judging this choice
683  Plep = TMath::Sqrt( T * (T + (2.0 * LepMass))); // ok is sqrt(E2 - m2)
684  Q3 = TMath::Sqrt(Plep*Plep + Enu*Enu - 2.0 * Plep * Enu * Costh);
685 
686  // Don't bother doing hard work if the selected Q3 is greater than Q3Max
687  if (Q3 < fQ3Max){
688 
689  kinematics->SetKV(kKVTl, T);
690  kinematics->SetKV(kKVctl, Costh);
691 
692  // decide whether to accept or reject these kinematics
693  // AND set the chosen two-nucleon system
694 
695  // to save time, use a pre-calculated max cross-section XSecMax
696  // it doesn't matter what it is, as long as it is big enough.
697  // RIK asks can XSecMax can be pushed to the q0q3 part of the calculation
698  // where the XS doesn't depend much on Enu.
699  // instead, this implementation uses a rough dependence on log10(Enu).
700  // starting around 0.5 GeV, the log10(max) is linear vs. log10(Enu)
701  // 1.10*TMath::Power(10.0, 2.2504 * TMath::Log10(Enu) - 9.41158)
702 
703  if (FullDeltaNodelta == 1){
704  // this block for the user who wants all CC QE-like 2p2h events
705 
706  // extract xsecmax from the spline making process for C12 and other nuclei.
707  // plot Log10(E) on horizontal and Log10(xsecmax) vertical
708  // and fit a line. Use that plus 1.35 safety factors to limit the accept/reject loop.
709  double XSecMax = 1.35 * TMath::Power(10.0, XSecMaxPar1 * TMath::Log10(Enu) - XSecMaxPar2);
710  if (NuclearA > 12) XSecMax *= NuclearAfactorXSecMax; // Scale it by A, precomputed above.
711 
712  LOG("MEC", pDEBUG) << " T, Costh: " << T << ", " << Costh ;
713 
714 
715  // We need four different cross sections. Right now, pursue the
716  // inelegant method of calling XSec four times - there is
717  // definitely some runtime inefficiency here, but it is not awful
718 
719  // first, get delta-less all
720  if (NuPDG > 0) {
721  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNN);
722  }
723  else {
724  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterPP);
725  }
726  double XSec = fXSecModel->XSec(interaction, kPSTlctl);
727  // now get all with delta
728  interaction->ExclTagPtr()->SetResonance(genie::kP33_1232);
729  double XSecDelta = fXSecModel->XSec(interaction, kPSTlctl);
730  // get PN with delta
731  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNP);
732  double XSecDeltaPN = fXSecModel->XSec(interaction, kPSTlctl);
733  // now get delta-less PN
735  double XSecPN = fXSecModel->XSec(interaction, kPSTlctl);
736 
737  if (XSec > XSecMax) {
738  LOG("MEC", pERROR) << "XSec is > XSecMax for nucleus " << TgtPDG << " "
739  << XSec << " > " << XSecMax
740  << " don't let this happen.";
741  }
742  assert(XSec <= XSecMax);
743  accept = XSec > XSecMax*rnd->RndKine().Rndm();
744  LOG("MEC", pINFO) << "Xsec, Max, Accept: " << XSec << ", "
745  << XSecMax << ", " << accept;
746 
747  if(accept){
748  // If it passes the All cross section we still need to do two things:
749  // * Was the initial state pn or not?
750  // * Do we assign the reaction to have had a Delta on the inside?
751 
752  // PDD means from the part of the XSec with an internal Delta line
753  // that (at the diagram level) did not produce a pion in the final state.
754 
755  bool isPDD = false;
756 
757  // Find out if we should use a pn initial state
758  double myrand = rnd->RndKine().Rndm();
759  double pnFraction = XSecPN / XSec;
760  LOG("MEC", pDEBUG) << "Test for pn: xsec_pn = " << XSecPN
761  << "; xsec = " << XSec
762  << "; pn_fraction = " << pnFraction
763  << "; random number val = " << myrand;
764 
765  if (myrand <= pnFraction) {
766  // yes it is, add a PN initial state to event record
767  event->AddParticle(kPdgClusterNP, kIStNucleonTarget,
768  1, -1, -1, -1, tempp4, v4);
769  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNP);
770 
771  // Its a pn, so test for Delta by comparing DeltaPN/PN
772  if (rnd->RndKine().Rndm() <= XSecDeltaPN / XSecPN) {
773  isPDD = true;
774  }
775  }
776  else {
777  // no it is not a PN, add either NN or PP initial state to event record.
778  if (NuPDG > 0) {
779  event->AddParticle(kPdgClusterNN, kIStNucleonTarget,
780  1, -1, -1, -1, tempp4, v4);
781  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNN);
782  }
783  else {
784  event->AddParticle(kPdgClusterPP, kIStNucleonTarget,
785  1, -1, -1, -1, tempp4, v4);
786  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterPP);
787  }
788  // its not pn, so test for Delta (XSecDelta-XSecDeltaPN)/(XSec-XSecPN)
789  // right, both numerator and denominator are total not pn.
790  if (rnd->RndKine().Rndm() <=
791  (XSecDelta - XSecDeltaPN) / (XSec - XSecPN)) {
792  isPDD = true;
793  }
794  }
795 
796  // now test whether we tagged this as a pion event
797  // and assign that fact to the Exclusive State tag
798  // later, we can query const XclsTag & xcls = interaction->ExclTag()
799  if (isPDD){
800  interaction->ExclTagPtr()->SetResonance(genie::kP33_1232);
801  }
802 
803 
804  } // end if accept
805  } // end if delta == 1
806 
807  /* One can make simpler versions of the above for the
808  FullDeltaNodelta == 2 (only delta)
809  or
810  FullDeltaNodelta == 3 (set Delta FF = 1, lose interference effect).
811  but I (Rik) don't see what the use-case is for these, genratorly speaking.
812  */
813 
814  }// end if passes q3 test
815  } // end while
816 
817  // -- finish lepton kinematics
818  // If the code got here, then we accepted some kinematics
819  // and we can proceed to generate the final state.
820 
821  // define coordinate system wrt neutrino: z along neutrino, xy perp
822 
823  // Cos theta gives us z, the rest in xy:
824  double PlepZ = Plep * Costh;
825  double PlepXY = Plep * TMath::Sqrt(1. - TMath::Power(Costh,2));
826 
827  // random rotation about unit vector for phi direction
828  double phi= 2 * kPi * rnd->RndLep().Rndm();
829  // now fill x and y from PlepXY
830  double PlepX = PlepXY * TMath::Cos(phi);
831  double PlepY = PlepXY * TMath::Sin(phi);
832 
833  // Rotate lepton momentum vector from the reference frame (x'y'z') where
834  // {z':(neutrino direction), z'x':(theta plane)} to the LAB
835  TVector3 unit_nudir = event->Probe()->P4()->Vect().Unit();
836  TVector3 p3l(PlepX, PlepY, PlepZ);
837  p3l.RotateUz(unit_nudir);
838 
839  // Lepton 4-momentum in LAB
840  Elep = TMath::Sqrt(LepMass*LepMass + PlepX*PlepX + PlepY*PlepY + PlepZ*PlepZ);
841  TLorentzVector p4l(p3l,Elep);
842 
843  // Figure out the final-state primary lepton PDG code
844  int pdgc = interaction->FSPrimLepton()->PdgCode();
845  int momidx = event->ProbePosition();
846 
847  // -- Store Values ------------------------------------------//
848  // -- Interaction: Q2
849  Q0 = Enu - Elep;
850  Q2 = Q3*Q3 - Q0*Q0;
851  double gy = Q0 / Enu;
852  double gx = kinematics::Q2YtoX(Enu, 2 * kNucleonMass, Q2, gy);
853  double gW = kinematics::XYtoW(Enu, 2 * kNucleonMass, gx, gy);
854 
855  interaction->KinePtr()->SetQ2(Q2, true);
856  interaction->KinePtr()->Sety(gy, true);
857  interaction->KinePtr()->Setx(gx, true);
858  interaction->KinePtr()->SetW(gW, true);
859  interaction->KinePtr()->SetFSLeptonP4(p4l);
860  // in later methods
861  // will also set the four-momentum and W^2 of the hadron system.
862 
863  // -- Lepton
864  event->AddParticle( pdgc, kIStStableFinalState, momidx, -1, -1, -1, p4l, v4);
865 
866  LOG("MEC",pDEBUG) << "~~~ LEPTON DONE ~~~";
867 }
const double kPi
#define pERROR
Definition: Messenger.h:60
TRandom3 & RndLep(void) const
rnd number generator used by final state primary lepton generators
Definition: RandomGen.h:63
static const double kNucleonMass
Definition: Constants.h:78
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
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
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
const int kPdgClusterNP
Definition: PDGCodes.h:192
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:61
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
void SetResonance(Resonance_t res)
Definition: XclsTag.cxx:123
double XYtoW(double Ev, double M, double x, double y)
Definition: KineUtils.cxx:1087
const int kPdgClusterNN
Definition: PDGCodes.h:191
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
Summary information for an interaction.
Definition: Interaction.h:56
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:307
int ProbePdg(void) const
Definition: InitialState.h:65
const int kPdgTgtFreeN
Definition: PDGCodes.h:177
#define pINFO
Definition: Messenger.h:63
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:51
#define pWARN
Definition: Messenger.h:61
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:241
XclsTag * ExclTagPtr(void) const
Definition: Interaction.h:77
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:345
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
double Q2YtoX(double Ev, double M, double Q2, double y)
Definition: KineUtils.cxx:1117
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:188
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:253
int TgtPdg(void) const
Target * TgtPtr(void) const
Definition: InitialState.h:68
const XSecAlgorithmI * fXSecModel
Definition: MECGenerator.h:64
assert(nhit_max >=nhit_nbins)
const int kPdgTgtC12
Definition: PDGCodes.h:179
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
const InitialState & InitState(void) const
Definition: Interaction.h:69
double T
Definition: Xdiff_gwt.C:5
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void kinematics()
Definition: kinematics.C:10
double ProbeE(RefFrame_t rf) const
#define pDEBUG
Definition: Messenger.h:64
const int kPdgClusterPP
Definition: PDGCodes.h:193
void Algorithm::SetId ( const AlgId id)
virtualinherited

Set algorithm ID.

Definition at line 313 of file Algorithm.cxx.

Referenced by genie::Algorithm::AllowReconfig().

314 {
315  fID.Copy(id);
316 }
AlgId fID
algorithm name and configuration set
Definition: Algorithm.h:156
void Copy(const AlgId &id)
Definition: AlgId.cxx:78
void Algorithm::SetId ( string  name,
string  config 
)
virtualinherited

Definition at line 318 of file Algorithm.cxx.

319 {
320  fID.SetId(name, config);
321 }
const XML_Char * name
Definition: expat.h:151
Definition: config.py:1
AlgId fID
algorithm name and configuration set
Definition: Algorithm.h:156
void SetId(string name, string config="")
Definition: AlgId.cxx:70
const Algorithm * Algorithm::SubAlg ( const RgKey registry_key) const
inherited

Access the sub-algorithm pointed to by the input key, either from the local pool or from AlgFactory's pool

Definition at line 353 of file Algorithm.cxx.

References ana::assert(), genie::AlgFactory::GetAlgorithm(), genie::AlgFactory::Instance(), LOG, pERROR, and pINFO.

Referenced by genie::Algorithm::AllowReconfig(), genie::utils::gsl::FullQELdXSec::FullQELdXSec(), genie::NewQELXSec::Integrate(), genie::HadronTransporter::LoadConfig(), genie::NucleonDecayPrimaryVtxGenerator::LoadConfig(), genie::IBDXSecMap::LoadConfig(), genie::DISHadronicSystemGenerator::LoadConfig(), genie::EmpiricalMECPXSec2015::LoadConfig(), genie::COHElasticPXSec::LoadConfig(), genie::SlowRsclCharmDISPXSecLO::LoadConfig(), genie::AhrensNCELPXSec::LoadConfig(), genie::AlamSimoAtharVacasSKPXSec2014::LoadConfig(), genie::UnstableParticleDecayer::LoadConfig(), genie::ReinDFRPXSec::LoadConfig(), genie::PythiaHadronization::LoadConfig(), genie::BYPDF::LoadConfig(), genie::QPMDISPXSec::LoadConfig(), genie::AlvarezRusoCOHPiPXSec::LoadConfig(), genie::RosenbluthPXSec::LoadConfig(), genie::StrumiaVissaniIBDPXSec::LoadConfig(), LoadConfig(), genie::NNBarOscPrimaryVtxGenerator::LoadConfig(), genie::FermiMover::LoadConfig(), genie::AhrensDMELPXSec::LoadConfig(), genie::IMDAnnihilationPXSec::LoadConfig(), genie::QELEventGenerator::LoadConfig(), genie::KNOPythiaHadronization::LoadConfig(), genie::AivazisCharmPXSecLO::LoadConfig(), genie::RESHadronicSystemGenerator::LoadConfig(), genie::P33PaschosLalakulichPXSec::LoadConfig(), genie::BergerSehgalFMCOHPiPXSec2015::LoadConfig(), genie::QPMDMDISPXSec::LoadConfig(), genie::BergerSehgalCOHPiPXSec2015::LoadConfig(), genie::LwlynSmithQELCCPXSec::LoadConfig(), genie::NuElectronPXSec::LoadConfig(), genie::PaisQELLambdaPXSec::LoadConfig(), genie::ReinSehgalCOHPiPXSec::LoadConfig(), genie::NievesSimoVacasMECPXSec2016::LoadConfig(), genie::KovalenkoQELCharmPXSec::LoadConfig(), genie::CharmHadronization::LoadConfig(), genie::NuclearModelMap::LoadConfig(), genie::EventGenerator::LoadConfig(), genie::SmithMonizQELCCXSec::LoadConfig(), genie::BardinIMDRadCorPXSec::LoadConfig(), genie::QELEventGeneratorSM::LoadConfig(), genie::LwlynSmithFF::LoadConfig(), genie::MartiniEricsonChanfrayMarteauMECPXSec2016::LoadConfig(), genie::ReinSehgalRESPXSec::LoadConfig(), genie::QPMDISStrucFuncBase::LoadConfig(), genie::SmithMonizQELCCPXSec::LoadConfig(), genie::NievesQELCCPXSec::LoadConfig(), genie::HAIntranuke::LoadConfig(), genie::HAIntranuke2018::LoadConfig(), genie::HNIntranuke2018::LoadConfig(), genie::BSKLNBaseRESPXSec2014::LoadConfig(), genie::ReinSehgalSPPPXSec::LoadConfig(), genie::KNOHadronization::LoadConfig(), and genie::EventGeneratorListAssembler::LoadGenerator().

354 {
355 // Returns the sub-algorithm pointed to this algorithm's XML config file using
356 // the the values of the key.
357 // This method asserts the existence of these keys in the XML config.
358 // Note: Since only 1 parameter is used, the key value should contain both the
359 // algorithm name and its configuration set according to the usual scheme:
360 // namespace::algorithm_name/configuration_set
361 //
362  LOG("Algorithm", pINFO)
363  << "Fetching sub-alg within alg: " << this->Id().Key()
364  << " pointed to by key: " << registry_key;
365 
366  //-- if the algorithm owns its substructure:
367  // return the sub-algorithm from the local pool
368  //
369  if(fOwnsSubstruc) {
370  AlgMapConstIter iter = fOwnedSubAlgMp->find(registry_key);
371  if(iter!=fOwnedSubAlgMp->end()) return iter->second;
372  LOG("Algorithm", pERROR)
373  << "Owned sub-alg pointed to by key: " << registry_key
374  << " was not found within alg: " << this->Id().Key();
375  return 0;
376  }
377 
378  //-- if the algorithm does not own its substructure:
379  // return the sub-algorithm from the AlgFactory's pool
380  RgAlg alg ;
381  GetParam( registry_key, alg ) ;
382 
383  LOG("Algorithm", pINFO)
384  << "Registry key: " << registry_key << " points to algorithm: " << alg;
385 
386  // retrieve the Algorithm object from the the Algorithm factory
387  AlgFactory * algf = AlgFactory::Instance();
388  const Algorithm * algbase = algf->GetAlgorithm(alg.name, alg.config);
389  assert(algbase);
390 
391  return algbase;
392 }
#define pERROR
Definition: Messenger.h:60
AlgMap * fOwnedSubAlgMp
local pool for owned sub-algs (taken out of the factory pool)
Definition: Algorithm.h:167
bool fOwnsSubstruc
true if it owns its substructure (sub-algs,...)
Definition: Algorithm.h:155
Algorithm abstract base class.
Definition: Algorithm.h:54
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:86
#define pINFO
Definition: Messenger.h:63
map< string, Algorithm * >::const_iterator AlgMapConstIter
Definition: Algorithm.h:51
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
assert(nhit_max >=nhit_nbins)
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
string Key(void) const
Definition: AlgId.h:47

Member Data Documentation

bool genie::Algorithm::fAllowReconfig
protectedinherited
vector<Registry*> genie::Algorithm::fConfVect
protectedinherited

ideally these members should go private Registry will be access only through the GetParam method configurations registries from various sources the order of the vector is the precedence in case of repeated parameters position 0 -> Highest precedence

Definition at line 161 of file Algorithm.h.

AlgId genie::Algorithm::fID
protectedinherited

algorithm name and configuration set

Definition at line 156 of file Algorithm.h.

Referenced by genie::Algorithm::Id().

const NuclearModelI* genie::MECGenerator::fNuclModel
private

Definition at line 66 of file MECGenerator.h.

Referenced by GenerateFermiMomentum(), GenerateNSVInitialHadrons(), and LoadConfig().

AlgMap* genie::Algorithm::fOwnedSubAlgMp
protectedinherited

local pool for owned sub-algs (taken out of the factory pool)

Definition at line 167 of file Algorithm.h.

vector<bool> genie::Algorithm::fOwnerships
protectedinherited

ownership for every registry in fConfVect

Definition at line 164 of file Algorithm.h.

bool genie::Algorithm::fOwnsSubstruc
protectedinherited

true if it owns its substructure (sub-algs,...)

Definition at line 155 of file Algorithm.h.

TGenPhaseSpace genie::MECGenerator::fPhaseSpaceGenerator
mutableprivate

Definition at line 65 of file MECGenerator.h.

Referenced by DecayNucleonCluster().

double genie::MECGenerator::fQ3Max
private

Definition at line 68 of file MECGenerator.h.

Referenced by LoadConfig(), and SelectNSVLeptonKinematics().

AlgStatus_t genie::Algorithm::fStatus
protectedinherited

algorithm execution status

Definition at line 166 of file Algorithm.h.

Referenced by genie::Algorithm::GetStatus().

const XSecAlgorithmI* genie::MECGenerator::fXSecModel
mutableprivate

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