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

#include <HAIntranuke.h>

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

Public Member Functions

 HAIntranuke ()
 
 HAIntranuke (string config)
 
 ~HAIntranuke ()
 
void ProcessEventRecord (GHepRecord *event_rec) const
 
virtual void Configure (string param_set)
 
- Public Member Functions inherited from genie::Intranuke
 Intranuke ()
 
 Intranuke (string name)
 
 Intranuke (string name, string config)
 
 ~Intranuke ()
 
void Configure (const Registry &config)
 
void Configure (string param_set)
 
- Public Member Functions inherited from genie::EventRecordVisitorI
virtual ~EventRecordVisitorI ()
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 
virtual const Registry & GetConfig (void) const
 
Registry * GetOwnedConfig (void)
 
virtual const AlgId & Id (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 Algorithm * SubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Private Member Functions

void LoadConfig (void)
 
void SimulateHadronicFinalState (GHepRecord *ev, GHepParticle *p) const
 
void SimulateHadronicFinalStateKinematics (GHepRecord *ev, GHepParticle *p) const
 
INukeFateHA_t HadronFateHA (const GHepParticle *p) const
 
void Inelastic (GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
 
void ElasHA (GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
 
void InelasticHA (GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
 
double PiBounce (void) const
 
double PnBounce (void) const
 
bool HandleCompoundNucleus (GHepRecord *ev, GHepParticle *p, int mom) const
 

Private Attributes

unsigned int fNumIterations
 

Friends

class IntranukeTester
 

Additional Inherited Members

- Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
 
static string BuildParamVectSizeKey (const std::string &comm_name)
 
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
 
static string BuildParamMatRowSizeKey (const std::string &comm_name)
 
static string BuildParamMatColSizeKey (const std::string &comm_name)
 
- Protected Member Functions inherited from genie::Intranuke
void TransportHadrons (GHepRecord *ev) const
 
void GenerateVertex (GHepRecord *ev) const
 
bool NeedsRescattering (const GHepParticle *p) const
 
bool CanRescatter (const GHepParticle *p) const
 
bool IsInNucleus (const GHepParticle *p) const
 
void SetTrackingRadius (const GHepParticle *p) const
 
double GenerateStep (GHepRecord *ev, GHepParticle *p) const
 
- Protected Member Functions inherited from genie::EventRecordVisitorI
 EventRecordVisitorI ()
 
 EventRecordVisitorI (string name)
 
 EventRecordVisitorI (string name, string config)
 
- Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 
 Algorithm (string name)
 
 Algorithm (string name, string config)
 
void Initialize (void)
 
void DeleteConfig (void)
 
void DeleteSubstructure (void)
 
Registry * ExtractLocalConfig (const Registry &in) const
 
Registry * ExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...
 
template<class T >
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
 
template<class T >
bool GetParamDef (const RgKey &name, T &p, const T &def) const
 
template<class T >
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters. More...
 
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
template<class T >
int GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
 Handle to load matrix of parameters. More...
 
template<class T >
int GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
 
int GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 
- Protected Attributes inherited from genie::Intranuke
double fTrackingRadius
 tracking radius for the nucleus in the current event More...
 
TGenPhaseSpace fGenPhaseSpace
 a phase space generator More...
 
INukeHadroData * fHadroData
 a collection of h+N,h+A data & calculations More...
 
AlgFactory * fAlgf
 algorithm factory instance More...
 
const NuclearModelI * fNuclmodel
 nuclear model used to generate fermi momentum More...
 
int fRemnA
 remnant nucleus A More...
 
int fRemnZ
 remnant nucleus Z More...
 
TLorentzVector fRemnP4
 P4 of remnant system. More...
 
GEvGenMode_t fGMode
 event generation mode (lepton+A, hadron+A, ...) More...
 
double fR0
 effective nuclear size param More...
 
double fNR
 param multiplying the nuclear radius, determining how far to track hadrons beyond the "nuclear boundary" More...
 
double fNucRmvE
 binding energy to subtract from cascade nucleons More...
 
double fDelRPion
 factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement More...
 
double fDelRNucleon
 factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement More...
 
double fHadStep
 step size for intranuclear hadron transport More...
 
double fNucAbsFac
 absorption xsec correction factor (hN Mode) More...
 
double fNucCEXFac
 charge exchange xsec correction factor (hN Mode) More...
 
double fEPreEq
 threshold for pre-equilibrium reaction More...
 
double fFermiFac
 testing parameter to modify fermi momentum More...
 
double fFermiMomentum
 whether or not particle collision is pauli blocked More...
 
bool fDoFermi
 whether or not to do fermi mom. More...
 
bool fDoMassDiff
 whether or not to do mass diff. mode More...
 
bool fDoCompoundNucleus
 whether or not to do compound nucleus considerations More...
 
double fChPionMFPScale
 tweaking factors for tuning More...
 
double fNeutralPionMFPScale
 
double fPionFracCExScale
 
double fPionFracElasScale
 
double fPionFracInelScale
 
double fPionFracAbsScale
 
double fPionFracPiProdScale
 
double fNucleonMFPScale
 
double fNucleonFracCExScale
 
double fNucleonFracElasScale
 
double fNucleonFracInelScale
 
double fNucleonFracAbsScale
 
double fNucleonFracPiProdScale
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
vector< Registry * > fConfVect
 
vector< bool > fOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMap * fOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

Definition at line 51 of file HAIntranuke.h.

Constructor & Destructor Documentation

HAIntranuke::HAIntranuke ( )

Definition at line 91 of file HAIntranuke.cxx.

91  :
92 Intranuke("genie::HAIntranuke")
93 {
94 
95 }
HAIntranuke::HAIntranuke ( string  config)

Definition at line 97 of file HAIntranuke.cxx.

97  :
98 Intranuke("genie::HAIntranuke",config)
99 {
100 
101 }
HAIntranuke::~HAIntranuke ( )

Definition at line 103 of file HAIntranuke.cxx.

104 {
105 
106 }

Member Function Documentation

void HAIntranuke::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 1379 of file HAIntranuke.cxx.

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

1379  {
1380 
1381  Algorithm::Configure(param_set) ;
1382  this -> LoadConfig() ;
1383 }
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void HAIntranuke::ElasHA ( GHepRecord *  ev,
GHepParticle *  p,
INukeFateHA_t  fate 
) const
private

Definition at line 468 of file HAIntranuke.cxx.

References genie::GHepParticle::A(), genie::GHepRecord::AddParticle(), genie::INukeHadroFates::AsString(), genie::PDGLibrary::Find(), genie::Intranuke::fRemnA, genie::Intranuke::fRemnP4, genie::Intranuke::fRemnZ, genie::PDGLibrary::Instance(), genie::kIHAFtElas, genie::kIStStableFinalState, genie::kPdgNeutron, genie::kPdgProton, LOG, genie::utils::res::Mass(), genie::GHepParticle::Mass(), genie::GHepParticle::Name(), genie::GHepParticle::P4(), pDEBUG, genie::GHepParticle::Pdg(), PiBounce(), pINFO, PnBounce(), pNOTICE, pWARN, genie::GHepParticle::SetMomentum(), genie::exceptions::INukeException::SetReason(), genie::GHepParticle::SetStatus(), genie::GHepRecord::TargetNucleus(), and genie::utils::intranuke::TwoBodyKinematics().

Referenced by SimulateHadronicFinalStateKinematics().

470 {
471  // scatters particle within nucleus, copy of hN code meant to run only once
472  // in hA mode
473 
474 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
475  LOG("HAIntranuke", pDEBUG)
476  << "ElasHA() is invoked for a : " << p->Name()
477  << " whose fate is : " << INukeHadroFates::AsString(fate);
478 #endif
479 
480  if(fate!=kIHAFtElas)
481  {
482  LOG("HAIntranuke", pWARN)
483  << "ElasHA() cannot handle fate: " << INukeHadroFates::AsString(fate);
484  return;
485  }
486 
487  // check remnants
488  if(fRemnA<0 || fRemnZ<0) // best to stop it here and not try again.
489  {
490  LOG("HAIntranuke", pWARN) << "Invalid Nucleus! : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
491  p->SetStatus(kIStStableFinalState);
492  ev->AddParticle(*p);
493  return;
494  }
495 
496  // vars for incoming particle, target, and scattered pdg codes
497  int pcode = p->Pdg();
498  double Mp = p->Mass();
499  double Mt = 0.;
500  if (ev->TargetNucleus()->A()==fRemnA)
501  { Mt = PDGLibrary::Instance()->Find(ev->TargetNucleus()->Pdg())->Mass(); }
502  else
503  {
504  Mt = fRemnP4.M();
505  }
506  TLorentzVector t4PpL = *p->P4();
507  TLorentzVector t4PtL = fRemnP4;
508  double C3CM = 0.0;
509 
510  // calculate scattering angle
511  if(pcode==kPdgNeutron||pcode==kPdgProton) C3CM = TMath::Cos(this->PnBounce());
512  else C3CM = TMath::Cos(this->PiBounce());
513 
514  // calculate final 4 momentum of probe
515  TLorentzVector t4P3L, t4P4L;
516 
517  if (!utils::intranuke::TwoBodyKinematics(Mp,Mt,t4PpL,t4PtL,t4P3L,t4P4L,C3CM,fRemnP4))
518  {
519  LOG("HAIntranuke", pNOTICE) << "ElasHA() failed";
520  exceptions::INukeException exception;
521  exception.SetReason("TwoBodyKinematics failed in ElasHA");
522  throw exception;
523  }
524 
525  // Update probe particle
526  p->SetMomentum(t4P3L);
527  p->SetStatus(kIStStableFinalState);
528 
529  // Update Remnant nucleus
530  fRemnP4 = t4P4L;
531  LOG("HAIntranuke",pINFO)
532  << "C3cm = " << C3CM;
533  LOG("HAIntranuke",pINFO)
534  << "|p3| = " << t4P3L.Vect().Mag() << ", E3 = " << t4P3L.E() << ",Mp = " << Mp;
535  LOG("HAIntranuke",pINFO)
536  << "|p4| = " << fRemnP4.Vect().Mag() << ", E4 = " << fRemnP4.E() << ",Mt = " << Mt;
537 
538  ev->AddParticle(*p);
539 
540 }
genie::exceptions::INukeException
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
Definition: INukeException.h:38
genie::Intranuke::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke.h:98
genie::GHepParticle::P4
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
genie::HAIntranuke::PnBounce
double PnBounce(void) const
Definition: HAIntranuke.cxx:415
genie::HAIntranuke::PiBounce
double PiBounce(void) const
Definition: HAIntranuke.cxx:363
genie::utils::res::Mass
double Mass(Resonance_t res)
resonance mass (GeV)
Definition: BaryonResUtils.cxx:436
genie::GHepParticle::SetMomentum
void SetMomentum(const TLorentzVector &p4)
Definition: GHepParticle.cxx:236
genie::GHepParticle::Mass
double Mass(void) const
Mass that corresponds to the PDG code.
Definition: GHepParticle.cxx:143
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:63
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:135
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::utils::intranuke::TwoBodyKinematics
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition: INukeUtils.cxx:754
genie::GHepRecord::TargetNucleus
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
pINFO
#define pINFO
Definition: Messenger.h:62
genie::Intranuke::fRemnA
int fRemnA
remnant nucleus A
Definition: Intranuke.h:96
pWARN
#define pWARN
Definition: Messenger.h:60
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::GHepParticle::SetStatus
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:126
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::Intranuke::fRemnZ
int fRemnZ
remnant nucleus Z
Definition: Intranuke.h:97
genie::GHepRecord::AddParticle
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
genie::GHepParticle::A
int A(void) const
Definition: GHepParticle.cxx:184
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
INukeFateHA_t HAIntranuke::HadronFateHA ( const GHepParticle *  p) const
private

Definition at line 212 of file HAIntranuke.cxx.

References genie::INukeHadroFates::AsString(), genie::Intranuke::fHadroData, genie::Intranuke::fNucleonFracAbsScale, genie::Intranuke::fNucleonFracCExScale, genie::Intranuke::fNucleonFracElasScale, genie::Intranuke::fNucleonFracInelScale, genie::Intranuke::fNucleonFracPiProdScale, genie::Intranuke::fPionFracAbsScale, genie::Intranuke::fPionFracCExScale, genie::Intranuke::fPionFracElasScale, genie::Intranuke::fPionFracInelScale, genie::Intranuke::fPionFracPiProdScale, genie::INukeHadroData::Frac(), genie::RandomGen::Instance(), genie::kIHAFtAbs, genie::kIHAFtCEx, genie::kIHAFtElas, genie::kIHAFtInelas, genie::kIHAFtPiProd, genie::kIHAFtUndefined, genie::GHepParticle::KinE(), genie::kPdgKM, genie::kPdgKP, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, genie::controls::kRjMaxIterations, LOG, genie::units::MeV, genie::GHepParticle::Name(), pDEBUG, genie::GHepParticle::Pdg(), pINFO, pWARN, and genie::RandomGen::RndFsi().

Referenced by SimulateHadronicFinalState().

213 {
214 // Select a hadron fate in HA mode
215 //
216  RandomGen * rnd = RandomGen::Instance();
217 
218  // get pdgc code & kinetic energy in MeV
219  int pdgc = p->Pdg();
220  double ke = p->KinE() / units::MeV;
221 
222  LOG("HAIntranuke", pINFO)
223  << "Selecting hA fate for " << p->Name() << " with KE = " << ke << " MeV";
224 
225  // try to generate a hadron fate
226  unsigned int iter = 0;
227  while(iter++ < kRjMaxIterations) {
228 
229  // handle pions
230  //
231  if (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0) {
232 
233  double frac_cex = fHadroData->Frac(pdgc, kIHAFtCEx, ke);
234  double frac_elas = fHadroData->Frac(pdgc, kIHAFtElas, ke);
235  double frac_inel = fHadroData->Frac(pdgc, kIHAFtInelas, ke);
236  double frac_abs = fHadroData->Frac(pdgc, kIHAFtAbs, ke);
237  double frac_piprod = fHadroData->Frac(pdgc, kIHAFtPiProd, ke);
238 
239  LOG("HAIntranuke", pINFO)
240  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCEx) << "} = " << frac_cex
241  << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
242  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
243  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs
244  << "\n frac{" << INukeHadroFates::AsString(kIHAFtPiProd) << "} = " << frac_piprod;
245 
246  // apply external tweaks to fractions
247  frac_cex *= fPionFracCExScale;
248  frac_elas *= fPionFracElasScale;
249  frac_inel *= fPionFracInelScale;
250  frac_abs *= fPionFracAbsScale;
251  frac_piprod *= fPionFracPiProdScale;
252  double frac_rescale = 1./(frac_cex + frac_elas + frac_inel + frac_abs + frac_piprod);
253  frac_cex *= frac_rescale;
254  frac_elas *= frac_rescale;
255  frac_inel *= frac_rescale;
256  frac_abs *= frac_rescale;
257  frac_piprod *= frac_rescale;
258 
259  LOG("HAIntranuke", pINFO)
260  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCEx) << "} = " << frac_cex
261  << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
262  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
263  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs
264  << "\n frac{" << INukeHadroFates::AsString(kIHAFtPiProd) << "} = " << frac_piprod;
265 
266  // compute total fraction (can be <1 if fates have been switched off)
267  double tf = frac_cex +
268  frac_elas +
269  frac_inel +
270  frac_abs +
271  frac_piprod;
272 
273  double r = tf * rnd->RndFsi().Rndm();
274 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
275  LOG("HAIntranuke", pDEBUG) << "r = " << r << " (max = " << tf << ")";
276 #endif
277  double cf=0; // current fraction
278  if(r < (cf += frac_cex )) return kIHAFtCEx; // cex
279  if(r < (cf += frac_elas )) return kIHAFtElas; // elas
280  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
281  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
282  if(r < (cf += frac_piprod )) return kIHAFtPiProd; // pi prod
283 
284  LOG("HAIntranuke", pWARN)
285  << "No selection after going through all fates! "
286  << "Total fraction = " << tf << " (r = " << r << ")";
287  }
288 
289  // handle nucleons
290  else if (pdgc==kPdgProton || pdgc==kPdgNeutron) {
291 
292  double frac_cex = fHadroData->Frac(pdgc, kIHAFtCEx, ke);
293  double frac_elas = fHadroData->Frac(pdgc, kIHAFtElas, ke);
294  double frac_inel = fHadroData->Frac(pdgc, kIHAFtInelas, ke);
295  double frac_abs = fHadroData->Frac(pdgc, kIHAFtAbs, ke);
296  double frac_pipro = fHadroData->Frac(pdgc, kIHAFtPiProd, ke);
297 
298  // apply external tweaks to fractions
299  frac_cex *= fNucleonFracCExScale;
300  frac_elas *= fNucleonFracElasScale;
301  frac_inel *= fNucleonFracInelScale;
302  frac_abs *= fNucleonFracAbsScale;
303  frac_pipro *= fNucleonFracPiProdScale;
304  double frac_rescale = 1./(frac_cex + frac_elas + frac_inel + frac_abs + frac_pipro);
305  frac_cex *= frac_rescale;
306  frac_elas *= frac_rescale;
307  frac_inel *= frac_rescale;
308  frac_abs *= frac_rescale;
309  frac_pipro *= frac_rescale;
310 
311  LOG("HAIntranuke", pDEBUG)
312  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCEx) << "} = " << frac_cex
313  << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
314  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
315  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs
316  << "\n frac{" << INukeHadroFates::AsString(kIHAFtPiProd) << "} = " << frac_pipro;
317 
318  // compute total fraction (can be <1 if fates have been switched off)
319  double tf = frac_cex +
320  frac_elas +
321  frac_inel +
322  frac_abs +
323  frac_pipro;
324 
325  double r = tf * rnd->RndFsi().Rndm();
326 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
327  LOG("HAIntranuke", pDEBUG) << "r = " << r << " (max = " << tf << ")";
328 #endif
329  double cf=0; // current fraction
330  if(r < (cf += frac_cex )) return kIHAFtCEx; // cex
331  if(r < (cf += frac_elas )) return kIHAFtElas; // elas
332  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
333  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
334  if(r < (cf += frac_pipro )) return kIHAFtPiProd; // pi prod
335 
336  LOG("HAIntranuke", pWARN)
337  << "No selection after going through all fates! "
338  << "Total fraction = " << tf << " (r = " << r << ")";
339  }
340  // handle kaons
341  else if (pdgc==kPdgKP || pdgc==kPdgKM) {
342  double frac_inel = fHadroData->Frac(pdgc, kIHAFtInelas, ke);
343  double frac_abs = fHadroData->Frac(pdgc, kIHAFtAbs, ke);
344  LOG("HAIntranuke", pDEBUG)
345  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
346  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs;
347  // compute total fraction (can be <1 if fates have been switched off)
348  double tf = frac_inel +
349  frac_abs;
350  double r = tf * rnd->RndFsi().Rndm();
351 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
352  LOG("HAIntranuke", pDEBUG) << "r = " << r << " (max = " << tf << ")";
353 #endif
354  double cf=0; // current fraction
355  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
356  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
357  }
358  }//iterations
359 
360  return kIHAFtUndefined;
361 }
genie::RandomGen::RndFsi
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:59
genie::Intranuke::fNucleonFracInelScale
double fNucleonFracInelScale
Definition: Intranuke.h:127
genie::Intranuke::fHadroData
INukeHadroData * fHadroData
a collection of h+N,h+A data &amp; calculations
Definition: Intranuke.h:93
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::Intranuke::fNucleonFracCExScale
double fNucleonFracCExScale
Definition: Intranuke.h:125
genie::Intranuke::fNucleonFracAbsScale
double fNucleonFracAbsScale
Definition: Intranuke.h:128
genie::units::MeV
static constexpr double MeV
Definition: Units.h:129
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
genie::INukeHadroData::Frac
double Frac(int hpdgc, INukeFateHA_t fate, double ke) const
Definition: INukeHadroData.cxx:1053
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:63
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:135
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::Intranuke::fPionFracPiProdScale
double fPionFracPiProdScale
Definition: Intranuke.h:123
genie::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:173
genie::Intranuke::fPionFracInelScale
double fPionFracInelScale
Definition: Intranuke.h:121
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
pINFO
#define pINFO
Definition: Messenger.h:62
pWARN
#define pWARN
Definition: Messenger.h:60
genie::GHepParticle::KinE
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
Definition: GHepParticle.cxx:159
genie::Intranuke::fPionFracAbsScale
double fPionFracAbsScale
Definition: Intranuke.h:122
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::Intranuke::fPionFracCExScale
double fPionFracCExScale
Definition: Intranuke.h:119
genie::Intranuke::fPionFracElasScale
double fPionFracElasScale
Definition: Intranuke.h:120
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::Intranuke::fNucleonFracPiProdScale
double fNucleonFracPiProdScale
Definition: Intranuke.h:129
genie::controls::kRjMaxIterations
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
genie::Intranuke::fNucleonFracElasScale
double fNucleonFracElasScale
Definition: Intranuke.h:126
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
bool HAIntranuke::HandleCompoundNucleus ( GHepRecord *  ev,
GHepParticle *  p,
int  mom 
) const
privatevirtual

Implements genie::Intranuke.

Definition at line 1372 of file HAIntranuke.cxx.

1374 {
1375  // only relevant for hN mode
1376  return false;
1377 }
void HAIntranuke::Inelastic ( GHepRecord *  ev,
GHepParticle *  p,
INukeFateHA_t  fate 
) const
private

Definition at line 704 of file HAIntranuke.cxx.

References genie::GHepRecord::AddParticle(), genie::INukeHadroFates::AsString(), genie::Intranuke::fDoFermi, genie::Intranuke::fFermiFac, genie::Intranuke::fFermiMomentum, genie::Intranuke::fHadroData, genie::PDGLibrary::Find(), genie::GHepParticle::FirstMother(), genie::Intranuke::fNuclmodel, genie::Intranuke::fNucRmvE, genie::Intranuke::fRemnA, genie::Intranuke::fRemnP4, genie::Intranuke::fRemnZ, genie::NuclearModelI::GenerateNucleon(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::INukeHadroData::IntBounce(), genie::pdg::IsKaon(), genie::pdg::IsNeutron(), genie::pdg::IsNeutronOrProton(), genie::pdg::IsPion(), genie::pdg::IsProton(), genie::kIHAFtAbs, genie::kIHAFtPiProd, genie::kIHNFtAbs, genie::kIMdHA, genie::GHepParticle::KinE(), genie::kIStNucleonClusterTarget, genie::kIStStableFinalState, genie::kPdgCompNuclCluster, genie::kPdgKM, genie::kPdgKP, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, genie::constants::kPi, genie::GHepParticle::LastMother(), LOG, genie::units::MeV, genie::NuclearModelI::Momentum3(), genie::GHepParticle::Name(), genie::units::ns, genie::GHepParticle::P4(), pDEBUG, genie::GHepParticle::Pdg(), genie::utils::intranuke::PhaseSpaceDecay(), pINFO, genie::utils::intranuke::PionProduction(), pNOTICE, genie::PDGCodeList::push_back(), pWARN, genie::RandomGen::RndFsi(), genie::GHepParticle::SetFirstMother(), genie::Target::SetHitNucPdg(), genie::GHepParticle::SetLastMother(), genie::GHepParticle::SetMomentum(), genie::GHepParticle::SetPdgCode(), genie::exceptions::INukeException::SetReason(), genie::GHepParticle::SetStatus(), genie::GHepRecord::TargetNucleus(), genie::utils::intranuke::TwoBodyKinematics(), and genie::GHepParticle::X4().

Referenced by SimulateHadronicFinalStateKinematics().

706 {
707 
708  // Aaron Meyer (05/25/10)
709  //
710  // Called to handle all absorption and pi production reactions
711  //
712  // Nucleons -> Reaction approximated by exponential decay in p+n (sum) space,
713  // gaussian in p-n (difference) space
714  // -fit to hN simulations p C, Fe, Pb at 200, 800 MeV
715  // -get n from isospin, np-nn smaller by 2
716  // Pions -> Reaction approximated with a modified gaussian in p+n space,
717  // normal gaussian in p-n space
718  // -based on fits to multiplicity distributions of hN model
719  // for pi+ C, Fe, Pb at 250, 500 MeV
720  // -fit sum and diff of nn, np to Gaussian
721  // -get pi0 from isospin, np-nn smaller by 2
722  // -get pi- from isospin, np-nn smaller by 4
723  // -add 2-body absorption to better match McKeown data
724  // Kaons -> no guidance, use same code as pions.
725  //
726  // Normally distributed random number generated using Box-Muller transformation
727  //
728  // Pion production reactions rescatter pions in nucleus, otherwise unchanged from
729  // older versions of GENIE
730  //
731 
732 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
733  LOG("HAIntranuke", pDEBUG)
734  << "Inelastic() is invoked for a : " << p->Name()
735  << " whose fate is : " << INukeHadroFates::AsString(fate);
736 #endif
737 
738  bool allow_dup = true;
739  PDGCodeList list(allow_dup); // list of final state particles
740 
741  // only absorption/pipro fates allowed
742  if (fate == kIHAFtPiProd) {
743 
744  GHepParticle s1(*p);
745  GHepParticle s2(*p);
746  GHepParticle s3(*p);
747 
748  bool success = utils::intranuke::PionProduction(
749  ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4, fDoFermi,fFermiFac,fFermiMomentum,fNuclmodel);
750 
751  if (success){
752  LOG ("HAIntranuke",pINFO) << " successful pion production fate";
753  // set status of particles and conserve charge/baryon number
754  s1.SetStatus(kIStStableFinalState); //should be redundant
755  // if (pdg::IsPion(s2->Pdg())) s2->SetStatus(kIStHadronInTheNucleus);
756  s2.SetStatus(kIStStableFinalState);
757  // if (pdg::IsPion(s3->Pdg())) s3->SetStatus(kIStHadronInTheNucleus);
758  s3.SetStatus(kIStStableFinalState);
759 
760  ev->AddParticle(s1);
761  ev->AddParticle(s2);
762  ev->AddParticle(s3);
763 
764  return;
765  }
766  else {
767  LOG("HAIntranuke", pNOTICE) << "Error: could not create pion production final state";
768  exceptions::INukeException exception;
769  exception.SetReason("PionProduction kinematics failed - retry kinematics");
770  throw exception;
771  }
772  }
773 
774  else if (fate==kIHAFtAbs)
775 // tuned for pions - mixture of 2-body and many-body
776 // use same for kaons as there is no guidance
777  {
778  // Instances for reference
779  PDGLibrary * pLib = PDGLibrary::Instance();
780  RandomGen * rnd = RandomGen::Instance();
781 
782  double ke = p->KinE() / units::MeV;
783  int pdgc = p->Pdg();
784 
785  if (fRemnA<2)
786  {
787  LOG("HAIntranuke", pNOTICE) << "stop propagation - could not create absorption final state: too few particles";
788  p->SetStatus(kIStStableFinalState);
789  ev->AddParticle(*p);
790  return;
791  }
792  if (fRemnZ<1 && (pdgc==kPdgPiM || pdgc==kPdgKM))
793  {
794  LOG("HAIntranuke", pNOTICE) << "stop propagation - could not create absorption final state: Pi- or K- cannot be absorbed by only neutrons";
795  p->SetStatus(kIStStableFinalState);
796  ev->AddParticle(*p);
797  return;
798  }
799  if (fRemnA-fRemnZ<1 && (pdgc==kPdgPiP || pdgc==kPdgKP))
800  {
801  LOG("HAIntranuke", pINFO) << "stop propagation - could not create absorption final state: Pi+ or K+ cannot be absorbed by only protons";
802  p->SetStatus(kIStStableFinalState);
803  ev->AddParticle(*p);
804  return;
805  }
806 
807  // for now, empirical split between multi-nucleon absorption and pi d -> N N
808  //
809  // added 03/21/11 - Aaron Meyer
810  //
811  if (pdg::IsPion(pdgc) && rnd->RndFsi().Rndm()<1.14*(.903-0.00189*fRemnA)*(1.35-0.00467*ke))
812  { // pi d -> N N, probability determined empirically with McKeown data
813 
814  INukeFateHN_t fate_hN=kIHNFtAbs;
815  int t1code,t2code,scode,s2code;
816  double ppcnt = (double) fRemnZ / (double) fRemnA; // % of protons
817 
818  // choose target nucleon
819  // -- fates weighted by values from Engel, Mosel...
820  if (pdgc==kPdgPiP) {
821  double Prob_pipd_pp=2.*ppcnt*(1.-ppcnt);
822  double Prob_pipnn_pn=.083*(1.-ppcnt)*(1.-ppcnt);
823  if (rnd->RndFsi().Rndm()*(Prob_pipd_pp+Prob_pipnn_pn)<Prob_pipd_pp){
824  t1code=kPdgNeutron; t2code=kPdgProton;
825  scode=kPdgProton; s2code=kPdgProton;}
826  else{
827  t1code=kPdgNeutron; t2code=kPdgNeutron;
828  scode=kPdgProton; s2code=kPdgNeutron;}
829  }
830  if (pdgc==kPdgPiM) {
831  double Prob_pimd_nn=2.*ppcnt*(1.-ppcnt);
832  double Prob_pimpp_pn=.083*ppcnt*ppcnt;
833  if (rnd->RndFsi().Rndm()*(Prob_pimd_nn+Prob_pimpp_pn)<Prob_pimd_nn){
834  t1code=kPdgProton; t2code=kPdgNeutron;
835  scode=kPdgNeutron; s2code=kPdgNeutron; }
836  else{
837  t1code=kPdgProton; t2code=kPdgProton;
838  scode=kPdgProton; s2code=kPdgNeutron;}
839  }
840  else { // pi0
841  double Prob_pi0d_pn=0.88*ppcnt*(1.-ppcnt); // 2 * .44
842  double Prob_pi0pp_pp=.14*ppcnt*ppcnt;
843  double Prob_pi0nn_nn=.14*(1.-ppcnt)*(1.-ppcnt);
844  if (rnd->RndFsi().Rndm()*(Prob_pi0d_pn+Prob_pi0pp_pp+Prob_pi0nn_nn)<Prob_pi0d_pn){
845  t1code=kPdgNeutron; t2code=kPdgProton;
846  scode=kPdgNeutron; s2code=kPdgProton; }
847  else if (rnd->RndFsi().Rndm()*(Prob_pi0d_pn+Prob_pi0pp_pp+Prob_pi0nn_nn)<(Prob_pi0d_pn+Prob_pi0pp_pp)){
848  t1code=kPdgProton; t2code=kPdgProton;
849  scode=kPdgProton; s2code=kPdgProton; }
850  else {
851  t1code=kPdgNeutron; t2code=kPdgNeutron;
852  scode=kPdgNeutron; s2code=kPdgNeutron; }
853  }
854  LOG("HAIntranuke",pINFO) << "choose 2 body absorption, probe, fs = " << pdgc <<" "<< scode <<" "<<s2code;
855  // assign proper masses
856  //double M1 = pLib->Find(pdgc) ->Mass();
857  double M2_1 = pLib->Find(t1code)->Mass();
858  double M2_2 = pLib->Find(t2code)->Mass();
859  //double M2 = M2_1 + M2_2;
860  double M3 = pLib->Find(scode) ->Mass();
861  double M4 = pLib->Find(s2code)->Mass();
862 
863  // handle fermi momentum
864  double E2_1L, E2_2L;
865  TVector3 tP2_1L, tP2_2L;
866  //TLorentzVector dNucl_P4;
867  Target target(ev->TargetNucleus()->Pdg());
868  if(fDoFermi)
869  {
870  target.SetHitNucPdg(t1code);
871  fNuclmodel->GenerateNucleon(target);
872  tP2_1L=fFermiFac * fNuclmodel->Momentum3();
873  E2_1L = TMath::Sqrt(tP2_1L.Mag2() + M2_1*M2_1);
874 
875  target.SetHitNucPdg(t2code);
876  fNuclmodel->GenerateNucleon(target);
877  tP2_2L=fFermiFac * fNuclmodel->Momentum3();
878  E2_2L = TMath::Sqrt(tP2_2L.Mag2() + M2_2*M2_2);
879 
880  //dNucl_P4=TLorentzVector(tP2_1L+tP2_2L,E2_1L+E2_2L);
881  }
882  else
883  {
884  tP2_1L.SetXYZ(0.0, 0.0, 0.0);
885  E2_1L = M2_1;
886 
887  tP2_2L.SetXYZ(0.0, 0.0, 0.0);
888  E2_2L = M2_2;
889  }
890  TLorentzVector dNucl_P4=TLorentzVector(tP2_1L+tP2_2L,E2_1L+E2_2L);
891 
892  double E2L = E2_1L + E2_2L;
893 
894  // adjust p to reflect scattering
895  // get random scattering angle
896  double C3CM = fHadroData->IntBounce(p,t1code,scode,fate_hN);
897  if (C3CM<-1.)
898  {
899  LOG("HAIntranuke", pWARN) << "Inelastic() failed: IntBounce returned bad angle - try again";
900  exceptions::INukeException exception;
901  exception.SetReason("unphysical angle for hN scattering");
902  throw exception;
903  return;
904  }
905 
906  TLorentzVector t4P1L,t4P2L,t4P3L,t4P4L;
907  t4P1L=*p->P4();
908  t4P2L=TLorentzVector(TVector3(tP2_1L+tP2_2L),E2L);
909  double bindE=0.075; // set to fit McKeown data
910  //double bindE=0.0;
911  if (utils::intranuke::TwoBodyKinematics(M3,M4,t4P1L,t4P2L,t4P3L,t4P4L,C3CM,fRemnP4,bindE))
912  {
913  if (pdgc==kPdgPiP || pdgc==kPdgKP) fRemnZ++;
914  if (pdgc==kPdgPiM || pdgc==kPdgKM) fRemnZ--;
915  if (t1code==kPdgProton) fRemnZ--;
916  if (t2code==kPdgProton) fRemnZ--;
917  fRemnA-=2;
918 
919  fRemnP4-=dNucl_P4;
920 
921  // create t particles w/ appropriate momenta, code, and status
922  // Set target's mom to be the mom of the hadron that was cloned
923  GHepParticle* t1 = new GHepParticle(*p);
924  GHepParticle* t2 = new GHepParticle(*p);
925  t1->SetFirstMother(p->FirstMother());
926  t1->SetLastMother(p->LastMother());
927  t2->SetFirstMother(p->FirstMother());
928  t2->SetLastMother(p->LastMother());
929 
930  // adjust p to reflect scattering
931  t1->SetPdgCode(scode);
932  t1->SetMomentum(t4P3L);
933 
934  t2->SetPdgCode(s2code);
935  t2->SetMomentum(t4P4L);
936 
937  t1->SetStatus(kIStStableFinalState);
938  t2->SetStatus(kIStStableFinalState);
939 
940  ev->AddParticle(*t1);
941  ev->AddParticle(*t2);
942  delete t1;
943  delete t2;
944 
945  return;
946  }
947  else
948  {
949  LOG("HAIntranuke", pNOTICE) << "Inelastic in hA failed calling TwoBodyKineamtics";
950  exceptions::INukeException exception;
951  exception.SetReason("Pion absorption kinematics through TwoBodyKinematics failed");
952  throw exception;
953 
954  }
955 
956  } // end pi d -> N N
957  else // multi-nucleon
958  {
959 
960  // declare some parameters for double gaussian and determine values chosen
961  // parameters for proton and pi+, others come from isospin transformations
962 
963  double ns0=0; // mean - sum of nucleons
964  double nd0=0; // mean - difference of nucleons
965  double Sig_ns=0; // std dev - sum
966  double Sig_nd=0; // std dev - diff
967  double gam_ns=0; // exponential decay rate (for nucleons)
968 
969  if ( pdg::IsNeutronOrProton (pdgc) ) // nucleon probe
970  {
971  // antisymmetric about Z=N
972  if (fRemnA-fRemnZ > fRemnZ)
973  nd0 = 135.227 * TMath::Exp(-7.124*(fRemnA-fRemnZ)/double(fRemnA)) - 2.762;
974  else
975  nd0 = -135.227 * TMath::Exp(-7.124* fRemnZ /double(fRemnA)) + 4.914;
976 
977  Sig_nd = 2.034 + fRemnA * 0.007846;
978 
979  double c1 = 0.041 + ke * 0.0001525;
980  double c2 = -0.003444 - ke * 0.00002324;
981 //change last factor from 30 to 15 so that gam_ns always larger than 0
982 //add check to be certain
983  double c3 = 0.064 - ke * 0.000015;
984  gam_ns = c1 * TMath::Exp(c2*fRemnA) + c3;
985  if(gam_ns<0.002) gam_ns = 0.002;
986  //gam_ns = 10.;
987  LOG("HAIntranuke", pINFO) << "nucleon absorption";
988  LOG("HAIntranuke", pINFO) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
989  // LOG("HAIntranuke", pINFO) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
990  LOG("HAIntranuke", pINFO) << "--> gam_ns = " << gam_ns;
991  }
992  else if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) //pion probe
993  {
994  ns0 = .0001*(1.+ke/250.) * (fRemnA-10)*(fRemnA-10) + 3.5;
995  nd0 = (1.+ke/250.) - ((fRemnA/200.)*(1. + 2.*ke/250.));
996  Sig_ns = (10. + 4. * ke/250.)*TMath::Power(fRemnA/250.,0.9); //(1. - TMath::Exp(-0.02*fRemnA));
997  Sig_nd = 4*(1 - TMath::Exp(-0.03*ke));
998  LOG("HAIntranuke", pINFO) << "pion absorption";
999  LOG("HAIntranuke", pINFO) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
1000  LOG("HAIntranuke", pINFO) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
1001  }
1002  else if (pdgc==kPdgKP || pdgc==kPdgKM) // kaon probe
1003  {
1004  ns0 = (rnd->RndFsi().Rndm()>0.5?3:2);
1005  nd0 = 1.;
1006  Sig_ns = 0.1;
1007  Sig_nd = 0.1;
1008  LOG("HAIntranuke", pINFO) << "kaon absorption - set ns, nd later";
1009  // LOG("HAIntranuke", pINFO) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
1010  // LOG("HAIntranuke", pINFO) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
1011  }
1012  else
1013  {
1014  LOG("HAIntranuke", pWARN) << "Inelastic() cannot handle absorption reaction for " << p->Name();
1015  }
1016 
1017  // account for different isospin
1018  if (pdgc==kPdgPi0 || pdgc==kPdgNeutron) nd0-=2.;
1019  if (pdgc==kPdgPiM) nd0-=4.;
1020 
1021  int iter=0; // counter
1022  int np=0,nn=0; // # of p, # of n
1023  bool not_done=true;
1024  double u1 = 0, u2 = 0;
1025 
1026  while (not_done)
1027  {
1028  // infinite loop check
1029  if (iter>=100) {
1030  LOG("HAIntranuke", pNOTICE) << "Error: could not choose absorption final state";
1031  LOG("HAIntranuke", pNOTICE) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
1032  LOG("HAIntranuke", pNOTICE) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
1033  LOG("HAIntranuke", pNOTICE) << "--> gam_ns = " << gam_ns;
1034  LOG("HAIntranuke", pNOTICE) << "--> A = " << fRemnA << ", Z = " << fRemnZ << ", Energy = " << ke;
1035  exceptions::INukeException exception;
1036  exception.SetReason("Absorption choice of # of p,n failed");
1037  throw exception;
1038  }
1039  //here??
1040 
1041  // Box-Muller transform
1042  // Takes two uniform distribution random variables on (0,1]
1043  // Creates two normally distributed random variables on (0,inf)
1044 
1045  u1 = rnd->RndFsi().Rndm(); // uniform random variable 1
1046  u2 = rnd->RndFsi().Rndm(); // " " 2
1047  if (u1==0) u1 = rnd->RndFsi().Rndm();
1048  if (u2==0) u2 = rnd->RndFsi().Rndm(); // Just in case
1049 
1050  // normally distributed random variable
1051  double x2 = TMath::Sqrt(-2*TMath::Log(u1))*TMath::Sin(2*kPi*u2);
1052 
1053  double ns = 0;
1054 
1055  if ( pdg::IsNeutronOrProton (pdgc) ) //nucleon probe
1056  {
1057  ns = -TMath::Log(rnd->RndFsi().Rndm())/gam_ns; // exponential random variable
1058  }
1059  if ( pdg::IsKaon (pdgc) ) //charged kaon probe - either 2 or 3 nucleons to stay simple
1060  {
1061  ns = (rnd->RndFsi().Rndm()<0.5?2:3);
1062  }
1063  else if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) //pion probe
1064  {
1065  // Pion fit for sum takes for xs*exp((xs-x0)^2 / 2*sig_xs0)
1066  // Find random variable by first finding gaussian random variable
1067  // then throwing the value against a linear P.D.F.
1068  //
1069  // max is the maximum value allowed for the random variable (10 std + mean)
1070  // minimum allowed value is 0
1071 
1072  double max = ns0 + Sig_ns * 10;
1073  if(max>fRemnA) max=fRemnA;
1074  double x1 = 0;
1075  bool not_found = true;
1076  int iter2 = 0;
1077 
1078  while (not_found)
1079  {
1080  // infinite loop check
1081  if (iter2>=100)
1082  {
1083  LOG("HAIntranuke", pNOTICE) << "Error: stuck in random variable loop for ns";
1084  LOG("HAIntranuke", pNOTICE) << "--> mean of sum parent distr = " << ns0 << ", Stand dev = " << Sig_ns;
1085  LOG("HAIntranuke", pNOTICE) << "--> A = " << fRemnA << ", Z = " << fRemnZ << ", Energy = " << ke;
1086 
1087  exceptions::INukeException exception;
1088  exception.SetReason("Random number generator for choice of #p,n final state failed - unusual - redo kinematics");
1089  throw exception;
1090  }
1091 
1092  // calculate exponential random variable
1093  u1 = rnd->RndFsi().Rndm();
1094  u2 = rnd->RndFsi().Rndm();
1095  if (u1==0) u1 = rnd->RndFsi().Rndm();
1096  if (u2==0) u2 = rnd->RndFsi().Rndm();
1097  x1 = TMath::Sqrt(-2*TMath::Log(u1))*TMath::Cos(2*kPi*u2);
1098 
1099  ns = ns0 + Sig_ns * x1;
1100  if ( ns>max || ns<0 ) {iter2++; continue;}
1101  else if ( rnd->RndFsi().Rndm() > (ns/max) ) {iter2++; continue;}
1102  else {
1103  // accept this sum value
1104  not_found=false;
1105  }
1106  } //while(not_found)
1107  }//else pion
1108 
1109  double nd = nd0 + Sig_nd * x2; // difference (p-n) for both pion, nucleon probe
1110  if (pdgc==kPdgKP) // special for KP
1111  { if (ns==2) nd=0;
1112  if (ns>2) nd=1; }
1113 
1114  np = int((ns+nd)/2.+.5); // Addition of .5 for rounding correction
1115  nn = int((ns-nd)/2.+.5);
1116 
1117  LOG("HAIntranuke", pINFO) << "ns = "<<ns<<", nd = "<<nd<<", np = "<<np<<", nn = "<<nn;
1118  //LOG("HAIntranuke", pNOTICE) << "RemA = "<<fRemnA<<", RemZ = "<<fRemnZ<<", probe = "<<pdgc;
1119 
1120  /*if ((ns+nd)/2. < 0 || (ns-nd)/2. < 0) {iter++; continue;}
1121  else */
1122  //check for problems befor executing phase space 'decay'
1123  if (np < 0 || nn < 0 ) {iter++; continue;}
1124  else if (np + nn < 2. ) {iter++; continue;}
1125  else if ((np + nn == 2.) && pdg::IsNeutronOrProton (pdgc)) {iter++; continue;}
1126  else if (np > fRemnZ + ((pdg::IsProton(pdgc) || pdgc==kPdgPiP || pdgc==kPdgKP)?1:0)
1127  - ((pdgc==kPdgPiM || pdgc==kPdgKM)?1:0)) {iter++; continue;}
1128  else if (nn > fRemnA-fRemnZ + ((pdg::IsNeutron(pdgc)||pdgc==kPdgPiM||pdgc==kPdgKM)?1:0)
1129  - ((pdgc==kPdgPiP||pdgc==kPdgKP)?1:0)) {iter++; continue;}
1130  else {
1131  not_done=false; //success
1132  LOG("HAIntranuke",pINFO) << "success, iter = " << iter << " np, nn = " << np << " " << nn;
1133  if (np+nn>86) // too many particles, scale down
1134  {
1135  double frac = 85./double(np+nn);
1136  np = int(np*frac);
1137  nn = int(nn*frac);
1138  }
1139 
1140  if ( (np==fRemnZ +((pdg::IsProton (pdgc)||pdgc==kPdgPiP||pdgc==kPdgKP)?1:0)-(pdgc==kPdgPiM||pdgc==kPdgKM?1:0))
1141  &&(nn==fRemnA-fRemnZ+((pdg::IsNeutron(pdgc)||pdgc==kPdgPiM||pdgc==kPdgKM)?1:0)-(pdgc==kPdgPiP||pdgc==kPdgKP?1:0)) )
1142  { // leave at least one nucleon in the nucleus to prevent excess momentum
1143  if (rnd->RndFsi().Rndm()<np/(double)(np+nn)) np--;
1144  else nn--;
1145  }
1146 
1147  LOG("HAIntranuke", pNOTICE) << "Final state chosen; # protons : "
1148  << np << ", # neutrons : " << nn;
1149  }
1150  } //while(not_done)
1151 
1152  // change remnants to reflect probe
1153  if ( pdgc==kPdgProton || pdgc==kPdgPiP || pdgc==kPdgKP) fRemnZ++;
1154  if ( pdgc==kPdgPiM || pdgc==kPdgKM) fRemnZ--;
1155  if ( pdg::IsNeutronOrProton (pdgc) ) fRemnA++;
1156 
1157  // PhaseSpaceDecay forbids anything over 18 particles
1158  //
1159  // If over 18 particles, split into 5 groups and perform decay on each group
1160  // Anything over 85 particles reduced to 85 in previous step
1161  //
1162  // 4 of the nucleons are used as "probes" as well as original probe,
1163  // with each one sharing 1/5 of the total incident momentum
1164  //
1165  // Note: choosing 5 groups and distributing momentum evenly was arbitrary
1166  // Needs to be revised later
1167 
1168  if ((np+nn)>18)
1169  {
1170  // code lists
1171  PDGCodeList list0(allow_dup);
1172  PDGCodeList list1(allow_dup);
1173  PDGCodeList list2(allow_dup);
1174  PDGCodeList list3(allow_dup);
1175  PDGCodeList list4(allow_dup);
1176  PDGCodeList* listar[5] = {&list0, &list1, &list2, &list3, &list4};
1177 
1178  //set up HadronClusters
1179  // simple for now, each (of 5) hadron cluster has 1/5 of mom and KE
1180 
1181  double probM = pLib->Find(pdgc) ->Mass();
1182  TVector3 pP3 = p->P4()->Vect() * (1./5.);
1183  double probKE = p->P4()->E() -probM;
1184  double clusKE = probKE * (1./5.);
1185  TLorentzVector clusP4(pP3,clusKE); //no mass
1186 
1187  TLorentzVector X4(*p->X4());
1188  GHepStatus_t ist = kIStNucleonClusterTarget;
1189 
1190  int mom = p->FirstMother();
1191 
1192  GHepParticle * p0 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1193  GHepParticle * p1 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1194  GHepParticle * p2 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1195  GHepParticle * p3 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1196  GHepParticle * p4 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1197 
1198  // To conserve 4-momenta
1199  // fRemnP4 -= probP4 + protP4*np_p + neutP4*(4-np_p) - *p->P4();
1200  fRemnP4 -= 5.*clusP4 - *p->P4();
1201 
1202  for (int i=0;i<(np+nn);i++)
1203  {
1204  if (i<np)
1205  {
1206  listar[i%5]->push_back(kPdgProton);
1207  fRemnZ--;
1208  }
1209  else listar[i%5]->push_back(kPdgNeutron);
1210  fRemnA--;
1211  }
1212  for (int i=0;i<5;i++)
1213  {
1214  LOG("HAIntranuke", pINFO) << "List" << i << " size: " << listar[i]->size();
1215  if (listar[i]->size() <2)
1216  {
1217  exceptions::INukeException exception;
1218  exception.SetReason("too few particles for Phase Space decay - try again");
1219  throw exception;
1220  }
1221  }
1222 
1223  // commented out to better fit with absorption reactions
1224  // Add the fermi energy of the nucleons to the phase space
1225  /*if(fDoFermi)
1226  {
1227  GHepParticle* p_ar[5] = {cl, p1, p2, p3, p4};
1228  for (int i=0;i<5;i++)
1229  {
1230  Target target(ev->TargetNucleus()->Pdg());
1231  TVector3 pBuf = p_ar[i]->P4()->Vect();
1232  double mBuf = p_ar[i]->Mass();
1233  double eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1234  TLorentzVector tSum(pBuf,eBuf);
1235  double mSum = 0.0;
1236  vector<int>::const_iterator pdg_iter;
1237  for(pdg_iter=++(listar[i]->begin());pdg_iter!=listar[i]->end();++pdg_iter)
1238  {
1239  target.SetHitNucPdg(*pdg_iter);
1240  fNuclmodel->GenerateNucleon(target);
1241  mBuf = pLib->Find(*pdg_iter)->Mass();
1242  mSum += mBuf;
1243  pBuf = fFermiFac * fNuclmodel->Momentum3();
1244  eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1245  tSum += TLorentzVector(pBuf,eBuf);
1246  fRemnP4 -= TLorentzVector(pBuf,eBuf-mBuf);
1247  }
1248  TLorentzVector dP4 = tSum + TLorentzVector(TVector3(0,0,0),-mSum);
1249  p_ar[i]->SetMomentum(dP4);
1250  }
1251  }*/
1252 
1253  bool success1 = utils::intranuke::PhaseSpaceDecay(ev,p0,*listar[0],fRemnP4,fNucRmvE,kIMdHA);
1254  bool success2 = utils::intranuke::PhaseSpaceDecay(ev,p1,*listar[1],fRemnP4,fNucRmvE,kIMdHA);
1255  bool success3 = utils::intranuke::PhaseSpaceDecay(ev,p2,*listar[2],fRemnP4,fNucRmvE,kIMdHA);
1256  bool success4 = utils::intranuke::PhaseSpaceDecay(ev,p3,*listar[3],fRemnP4,fNucRmvE,kIMdHA);
1257  bool success5 = utils::intranuke::PhaseSpaceDecay(ev,p4,*listar[4],fRemnP4,fNucRmvE,kIMdHA);
1258  if(success1 && success2 && success3 && success4 && success5)
1259  {
1260  LOG("HAIntranuke", pINFO)<<"Successful many-body absorption - n>=18";
1261  }
1262  else
1263  {
1264  // try to recover
1265  LOG("HAIntranuke", pWARN) << "PhaseSpace decay fails for HadrCluster- recovery likely incorrect - rethrow event";
1266  p->SetStatus(kIStStableFinalState);
1267  ev->AddParticle(*p);
1268  fRemnA+=np+nn;
1269  fRemnZ+=np;
1270  if ( pdgc==kPdgProton || pdgc==kPdgPiP ) fRemnZ--;
1271  if ( pdgc==kPdgPiM ) fRemnZ++;
1272  if ( pdg::IsNeutronOrProton (pdgc) ) fRemnA--;
1273  /* exceptions::INukeException exception;
1274  exception.SetReason("Phase space generation of absorption final state failed");
1275  throw exception;
1276  */
1277  }
1278 
1279  // delete cl;
1280  delete p0;
1281  delete p1;
1282  delete p2;
1283  delete p3;
1284  delete p4;
1285 
1286  }
1287  else // less than 18 particles pion
1288  {
1289  if (pdgc==kPdgKP) list.push_back(kPdgKP); //normally conserve strangeness
1290  if (pdgc==kPdgKM) list.push_back(kPdgKM);
1291  for (int i=0;i<np;i++)
1292  {
1293  list.push_back(kPdgProton);
1294  fRemnA--;
1295  fRemnZ--;
1296  }
1297  for (int i=0;i<nn;i++)
1298  {
1299  list.push_back(kPdgNeutron);
1300  fRemnA--;
1301  }
1302 
1303  // Library instance for reference
1304  //PDGLibrary * pLib = PDGLibrary::Instance();
1305 
1306  // commented out to better fit with absorption reactions
1307  // Add the fermi energy of the nucleons to the phase space
1308  /*if(fDoFermi)
1309  {
1310  Target target(ev->TargetNucleus()->Pdg());
1311  TVector3 pBuf = p->P4()->Vect();
1312  double mBuf = p->Mass();
1313  double eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1314  TLorentzVector tSum(pBuf,eBuf);
1315  double mSum = 0.0;
1316  vector<int>::const_iterator pdg_iter;
1317  for(pdg_iter=++(list.begin());pdg_iter!=list.end();++pdg_iter)
1318  {
1319  target.SetHitNucPdg(*pdg_iter);
1320  fNuclmodel->GenerateNucleon(target);
1321  mBuf = pLib->Find(*pdg_iter)->Mass();
1322  mSum += mBuf;
1323  pBuf = fFermiFac * fNuclmodel->Momentum3();
1324  eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1325  tSum += TLorentzVector(pBuf,eBuf);
1326  fRemnP4 -= TLorentzVector(pBuf,eBuf-mBuf);
1327  }
1328  TLorentzVector dP4 = tSum + TLorentzVector(TVector3(0,0,0),-mSum);
1329  p->SetMomentum(dP4);
1330  }*/
1331 
1332  LOG("HAIntranuke", pDEBUG)
1333  << "Remnant nucleus (A,Z) = (" << fRemnA << ", " << fRemnZ << ")";
1334  LOG("HAIntranuke", pINFO) << " list size: " << np+nn;
1335  if (np+nn <2)
1336  {
1337  exceptions::INukeException exception;
1338  exception.SetReason("too few particles for Phase Space decay - try again");
1339  throw exception;
1340  }
1341  // GHepParticle * cl = new GHepParticle(*p);
1342  // cl->SetPdgCode(kPdgDecayNuclCluster);
1343  bool success = utils::intranuke::PhaseSpaceDecay(ev,p,list,fRemnP4,fNucRmvE,kIMdHA);
1344  if (success)
1345  {
1346  LOG ("HAIntranuke",pINFO) << "Successful many-body absorption, n<=18";
1347  }
1348  else {
1349  // recover
1350  p->SetStatus(kIStStableFinalState);
1351  ev->AddParticle(*p);
1352  fRemnA+=np+nn;
1353  fRemnZ+=np;
1354  if ( pdgc==kPdgProton || pdgc==kPdgPiP ) fRemnZ--;
1355  if ( pdgc==kPdgPiM ) fRemnZ++;
1356  if ( pdg::IsNeutronOrProton (pdgc) ) fRemnA--;
1357  exceptions::INukeException exception;
1358  exception.SetReason("Phase space generation of absorption final state failed");
1359  throw exception;
1360  }
1361  }
1362  } // end multi-nucleon FS
1363  }
1364  else // not absorption/pipro
1365  {
1366  LOG("HAIntranuke", pWARN)
1367  << "Inelastic() can not handle fate: " << INukeHadroFates::AsString(fate);
1368  return;
1369  }
1370 }
genie::Intranuke::fFermiMomentum
double fFermiMomentum
whether or not particle collision is pauli blocked
Definition: Intranuke.h:112
genie::exceptions::INukeException
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
Definition: INukeException.h:38
genie::GHepParticle::SetFirstMother
void SetFirstMother(int m)
Definition: GHepParticle.h:132
genie::pdg::IsPion
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:326
genie::RandomGen::RndFsi
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:59
genie::Intranuke::fHadroData
INukeHadroData * fHadroData
a collection of h+N,h+A data &amp; calculations
Definition: Intranuke.h:93
genie::Intranuke::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke.h:98
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::GHepParticle::P4
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
genie::GHepParticle::SetMomentum
void SetMomentum(const TLorentzVector &p4)
Definition: GHepParticle.cxx:236
genie::units::ns
static constexpr double ns
Definition: Units.h:98
genie::pdg::IsKaon
bool IsKaon(int pdgc)
Definition: PDGUtils.cxx:331
genie::units::MeV
static constexpr double MeV
Definition: Units.h:129
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
genie::NuclearModelI::Momentum3
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:75
genie::PDGCodeList
A list of PDG codes.
Definition: PDGCodeList.h:32
genie::GHepParticle::LastMother
int LastMother(void) const
Definition: GHepParticle.h:67
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:63
genie::GHepParticle::FirstMother
int FirstMother(void) const
Definition: GHepParticle.h:66
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:135
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
genie::kPdgCompNuclCluster
const int kPdgCompNuclCluster
Definition: PDGCodes.h:217
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::INukeFateHN_t
enum genie::EINukeFateHN_t INukeFateHN_t
genie::utils::intranuke::PhaseSpaceDecay
bool PhaseSpaceDecay(GHepRecord *ev, GHepParticle *p, const PDGCodeList &pdgv, TLorentzVector &RemnP4, double NucRmvE, EINukeMode mode=kIMdHA)
general phase space decay method
Definition: INukeUtils.cxx:1555
genie::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:173
genie::utils::intranuke::PionProduction
bool PionProduction(GHepRecord *ev, GHepParticle *p, GHepParticle *s1, GHepParticle *s2, GHepParticle *s3, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI *Nuclmodel)
Definition: INukeUtils.cxx:1178
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::utils::intranuke::TwoBodyKinematics
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition: INukeUtils.cxx:754
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::GHepRecord::TargetNucleus
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
pINFO
#define pINFO
Definition: Messenger.h:62
genie::Intranuke::fRemnA
int fRemnA
remnant nucleus A
Definition: Intranuke.h:96
pWARN
#define pWARN
Definition: Messenger.h:60
genie::Intranuke::fNucRmvE
double fNucRmvE
binding energy to subtract from cascade nucleons
Definition: Intranuke.h:104
genie::INukeHadroData::IntBounce
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)
Definition: INukeHadroData.cxx:1294
genie::GHepParticle::KinE
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
Definition: GHepParticle.cxx:159
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::GHepParticle::SetLastMother
void SetLastMother(int m)
Definition: GHepParticle.h:133
genie::PDGLibrary
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
genie::GHepParticle::X4
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
genie::GHepParticle::SetStatus
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:126
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::Target::SetHitNucPdg
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:171
genie::pdg::IsNeutronOrProton
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:351
genie::Intranuke::fRemnZ
int fRemnZ
remnant nucleus Z
Definition: Intranuke.h:97
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
genie::Intranuke::fDoFermi
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke.h:113
genie::GHepRecord::AddParticle
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
genie::Intranuke::fFermiFac
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke.h:111
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::NuclearModelI::GenerateNucleon
virtual bool GenerateNucleon(const Target &) const =0
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
genie::GHepParticle
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
genie::Intranuke::fNuclmodel
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke.h:95
genie::PDGCodeList::push_back
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
genie::GHepStatus_t
enum genie::EGHepStatus GHepStatus_t
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
genie::GHepParticle::SetPdgCode
void SetPdgCode(int c)
Definition: GHepParticle.cxx:230
void HAIntranuke::InelasticHA ( GHepRecord *  ev,
GHepParticle *  p,
INukeFateHA_t  fate 
) const
private

Definition at line 542 of file HAIntranuke.cxx.

References genie::GHepRecord::AddParticle(), genie::INukeHadroFates::AsString(), genie::GHepParticle::E(), genie::Intranuke::fDoFermi, genie::Intranuke::fFermiFac, genie::Intranuke::fHadroData, genie::Intranuke::fNuclmodel, genie::Intranuke::fRemnA, genie::Intranuke::fRemnP4, genie::Intranuke::fRemnZ, genie::NuclearModelI::GenerateNucleon(), genie::RandomGen::Instance(), genie::INukeHadroData::IntBounce(), genie::kIHAFtCEx, genie::kIHAFtInelas, genie::kIHNFtCEx, genie::kIHNFtElas, genie::kIMdHA, genie::GHepParticle::KinE(), genie::kIStStableFinalState, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, LOG, genie::GHepParticle::Mass(), genie::NuclearModelI::Momentum3(), genie::GHepParticle::Name(), genie::GHepParticle::P4(), pDEBUG, genie::GHepParticle::Pdg(), pINFO, pNOTICE, genie::GHepRecord::Probe(), pWARN, genie::GHepParticle::Px(), genie::GHepParticle::Py(), genie::GHepParticle::Pz(), genie::RandomGen::RndFsi(), genie::GHepParticle::SetMomentum(), genie::GHepParticle::SetPdgCode(), genie::exceptions::INukeException::SetReason(), genie::GHepParticle::SetStatus(), genie::GHepRecord::TargetNucleus(), and genie::utils::intranuke::TwoBodyCollision().

Referenced by SimulateHadronicFinalStateKinematics().

544 {
545  // charge exch and inelastic - scatters particle within nucleus, hA version
546  // each are treated as quasielastic, particle scatters off single nucleon
547 
548 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
549  LOG("HAIntranuke", pDEBUG)
550  << "InelasticHA() is invoked for a : " << p->Name()
551  << " whose fate is : " << INukeHadroFates::AsString(fate);
552 #endif
553  if(ev->Probe() ) {
554  LOG("HAIntranuke", pINFO) << " probe KE = " << ev->Probe()->KinE();
555  }
556  if(fate!=kIHAFtCEx && fate!=kIHAFtInelas)
557  {
558  LOG("HAIntranuke", pWARN)
559  << "InelasticHA() cannot handle fate: " << INukeHadroFates::AsString(fate);
560  return;
561  }
562 
563  // Random number generator
564  RandomGen * rnd = RandomGen::Instance();
565 
566  // vars for incoming particle, target, and scattered pdg codes
567  int pcode = p->Pdg();
568  int tcode, scode, s2code;
569  double ppcnt = (double) fRemnZ / (double) fRemnA; // % of protons
570 
571  // Select a hadron fate in HN mode
572  INukeFateHN_t h_fate;
573  if (fate == kIHAFtCEx) h_fate = kIHNFtCEx;
574  else h_fate = kIHNFtElas;
575 
576  // Select a target randomly, weighted to #
577  // -- Unless, of course, the fate is CEx,
578  // -- in which case the target may be deterministic
579  // Also assign scattered particle code
580  if(fate==kIHAFtCEx)
581  {
582  if(pcode==kPdgPiP) {tcode = kPdgNeutron; scode = kPdgPi0; s2code = kPdgProton;}
583  else if(pcode==kPdgPiM) {tcode = kPdgProton; scode = kPdgPi0; s2code = kPdgNeutron;}
584  else if(pcode==kPdgPi0)
585  {
586  // for pi0
587  tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton) :(kPdgNeutron);
588  scode = (tcode == kPdgProton) ?(kPdgPiP) :(kPdgPiM);
589  s2code = (tcode == kPdgProton) ?(kPdgNeutron):(kPdgProton);
590  }
591  else if(pcode==kPdgProton) {tcode = kPdgNeutron; scode = kPdgNeutron; s2code = kPdgProton;}
592  else if(pcode==kPdgNeutron){tcode = kPdgProton; scode = kPdgProton; s2code = kPdgNeutron;}
593  else
594  { LOG("HAIntranuke", pWARN) << "InelasticHA() cannot handle fate: "
595  << INukeHadroFates::AsString(fate)
596  << " for particle " << p->Name();
597  return;
598  }
599  }
600  else
601  {
602  tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton):(kPdgNeutron);
603  // if(pcode == kPdgKP || pcode == kPdgKM) tcode = kPdgProton;
604  scode = pcode;
605  s2code = tcode;
606  }
607 
608  // check remnants
609  if ( fRemnA < 1 ) //we've blown nucleus apart, no need to retry anything - exit
610  {
611  LOG("HAIntranuke",pNOTICE) << "InelasticHA() stops : not enough nucleons";
612  p->SetStatus(kIStStableFinalState);
613  ev->AddParticle(*p);
614  return;
615  }
616  else if ( fRemnZ + (((pcode==kPdgProton)||(pcode==kPdgPiP))?1:0) - (pcode==kPdgPiM?1:0)
617  < ((( scode==kPdgProton)||( scode==kPdgPiP)) ?1:0) - (scode ==kPdgPiM ?1:0)
618  + (((s2code==kPdgProton)||(s2code==kPdgPiP)) ?1:0) - (s2code==kPdgPiM ?1:0) )
619  {
620  LOG("HAIntranuke",pWARN) << "InelasticHA() failed : too few protons in nucleus";
621  p->SetStatus(kIStStableFinalState);
622  ev->AddParticle(*p);
623  return; // another extreme case, best strategy is to exit and go to next event
624  }
625 
626  GHepParticle t(*p);
627  t.SetPdgCode(tcode);
628 
629  // set up fermi target
630  Target target(ev->TargetNucleus()->Pdg());
631  double tM = t.Mass();
632 
633  // handle fermi momentum
634  if(fDoFermi)
635  {
636  target.SetHitNucPdg(tcode);
637  fNuclmodel->GenerateNucleon(target);
638  TVector3 tP3 = fFermiFac * fNuclmodel->Momentum3();
639  double tE = TMath::Sqrt(tP3.Mag2()+ tM*tM);
640  t.SetMomentum(TLorentzVector(tP3,tE));
641  }
642  else
643  {
644  t.SetMomentum(0,0,0,tM);
645  }
646 
647  GHepParticle * cl = new GHepParticle(*p); // clone particle, to run IntBounce at proper energy
648  // calculate energy and momentum using invariant mass
649  double pM = p->Mass();
650  double E_p = ((*p->P4() + *t.P4()).Mag2() - tM*tM - pM*pM)/(2.0*tM);
651  double P_p = TMath::Sqrt(E_p*E_p - pM*pM);
652  cl->SetMomentum(TLorentzVector(P_p,0,0,E_p));
653  // momentum doesn't have to be in right direction, only magnitude
654  double C3CM = fHadroData->IntBounce(cl,tcode,scode,h_fate);
655  delete cl;
656  if (C3CM<-1.) // hope this doesn't occur too often - unphysical but we just pass it on
657  {
658  LOG("HAIntranuke", pWARN) << "unphysical angle chosen in InelasicHA - put particle outside nucleus";
659  p->SetStatus(kIStStableFinalState);
660  ev->AddParticle(*p);
661  return;
662  }
663  double KE1L = p->KinE();
664  double KE2L = t.KinE();
665  LOG("HAIntranuke",pINFO)
666  << " KE1L = " << KE1L << " " << KE1L << " KE2L = " << KE2L;
667  GHepParticle cl1(*p);
668  GHepParticle cl2(t);
669  bool success = utils::intranuke::TwoBodyCollision(ev,pcode,tcode,scode,s2code,C3CM,
670  &cl1,&cl2,fRemnA,fRemnZ,fRemnP4,kIMdHA);
671  if(success)
672  {
673  double P3L = TMath::Sqrt(cl1.Px()*cl1.Px() + cl1.Py()*cl1.Py() + cl1.Pz()*cl1.Pz());
674  double P4L = TMath::Sqrt(cl2.Px()*cl2.Px() + cl2.Py()*cl2.Py() + cl2.Pz()*cl2.Pz());
675  double E3L = cl1.KinE();
676  double E4L = cl2.KinE();
677  LOG ("HAIntranuke",pINFO) << "Successful quasielastic scattering or charge exchange";
678  LOG("HAIntranuke",pINFO)
679  << "C3CM = " << C3CM << "\n P3L, E3L = "
680  << P3L << " " << E3L << " P4L, E4L = "<< P4L << " " << E4L ;
681  if(ev->Probe() ) { LOG("HAIntranuke",pINFO)
682  << "P4L = " << P4L << " ;E4L= " << E4L << "\n probe KE = " << ev->Probe()->KinE() << "\n";
683  if (ev->Probe() && (E3L>ev->Probe()->E()||E4L>ev->Probe()->E())) //is this redundant?
684  {
685  exceptions::INukeException exception;
686  exception.SetReason("TwoBodyCollison gives KE> probe KE in hA simulation");
687  throw exception;
688  }
689  }
690  // always get here, even nucleon decay
691  ev->AddParticle(cl1);
692  ev->AddParticle(cl2);
693  LOG("HAIntranuke", pDEBUG) << "Nucleus : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
694 
695  } else
696  {
697  exceptions::INukeException exception;
698  exception.SetReason("TwoBodyCollison failed in hA simulation");
699  throw exception;
700  }
701 
702 }
genie::exceptions::INukeException
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
Definition: INukeException.h:38
genie::RandomGen::RndFsi
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:59
genie::Intranuke::fHadroData
INukeHadroData * fHadroData
a collection of h+N,h+A data &amp; calculations
Definition: Intranuke.h:93
genie::GHepParticle::E
double E(void) const
Get energy.
Definition: GHepParticle.h:91
genie::Intranuke::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke.h:98
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::GHepParticle::P4
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
genie::GHepParticle::SetMomentum
void SetMomentum(const TLorentzVector &p4)
Definition: GHepParticle.cxx:236
genie::GHepParticle::Mass
double Mass(void) const
Mass that corresponds to the PDG code.
Definition: GHepParticle.cxx:143
genie::Target::Mass
double Mass(void) const
Definition: Target.cxx:224
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
genie::NuclearModelI::Momentum3
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:75
genie::GHepRecord::Probe
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:284
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:63
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:135
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::INukeFateHN_t
enum genie::EINukeFateHN_t INukeFateHN_t
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::GHepRecord::TargetNucleus
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
pINFO
#define pINFO
Definition: Messenger.h:62
genie::Intranuke::fRemnA
int fRemnA
remnant nucleus A
Definition: Intranuke.h:96
pWARN
#define pWARN
Definition: Messenger.h:60
genie::INukeHadroData::IntBounce
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)
Definition: INukeHadroData.cxx:1294
genie::GHepParticle::KinE
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
Definition: GHepParticle.cxx:159
genie::GHepParticle::SetStatus
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:126
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::Intranuke::fRemnZ
int fRemnZ
remnant nucleus Z
Definition: Intranuke.h:97
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
genie::Intranuke::fDoFermi
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke.h:113
genie::GHepRecord::AddParticle
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
genie::Intranuke::fFermiFac
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke.h:111
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::NuclearModelI::GenerateNucleon
virtual bool GenerateNucleon(const Target &) const =0
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
genie::GHepParticle
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
genie::utils::intranuke::TwoBodyCollision
bool TwoBodyCollision(GHepRecord *ev, int pcode, int tcode, int scode, int s2code, double C3CM, GHepParticle *p, GHepParticle *t, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, EINukeMode mode=kIMdHA)
Intranuke utility functions.
Definition: INukeUtils.cxx:625
genie::Intranuke::fNuclmodel
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke.h:95
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
void HAIntranuke::LoadConfig ( void  )
privatevirtual

Implements genie::Intranuke.

Definition at line 1385 of file HAIntranuke.cxx.

References genie::INukeMode::AsString(), genie::Intranuke::fChPionMFPScale, genie::Intranuke::fDelRNucleon, genie::Intranuke::fDelRPion, genie::Intranuke::fDoCompoundNucleus, genie::Intranuke::fDoFermi, genie::Intranuke::fEPreEq, genie::Intranuke::fFermiFac, genie::Intranuke::fFermiMomentum, genie::Intranuke::fHadroData, genie::Intranuke::fHadStep, genie::Intranuke::fNeutralPionMFPScale, genie::Intranuke::fNR, genie::Intranuke::fNucAbsFac, genie::Intranuke::fNucCEXFac, genie::Intranuke::fNucleonFracAbsScale, genie::Intranuke::fNucleonFracCExScale, genie::Intranuke::fNucleonFracElasScale, genie::Intranuke::fNucleonFracInelScale, genie::Intranuke::fNucleonFracPiProdScale, genie::Intranuke::fNucleonMFPScale, genie::Intranuke::fNuclmodel, genie::Intranuke::fNucRmvE, genie::Intranuke::fPionFracAbsScale, genie::Intranuke::fPionFracCExScale, genie::Intranuke::fPionFracElasScale, genie::Intranuke::fPionFracInelScale, genie::Intranuke::fPionFracPiProdScale, genie::Intranuke::fR0, genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), genie::INukeHadroData::Instance(), genie::kIMdHA, LOG, pINFO, and genie::Algorithm::SubAlg().

Referenced by Configure().

1386 {
1387  // load hadronic cross sections
1388  fHadroData = INukeHadroData::Instance();
1389 
1390  // fermi momentum setup
1391  fNuclmodel = dynamic_cast<const NuclearModelI *>( this -> SubAlg("NuclearModel") ) ;
1392 
1393  // other intranuke config params
1394  GetParam( "NUCL-R0", fR0 ); // fm
1395  GetParam( "NUCL-NR", fNR );
1396 
1397  GetParam( "INUKE-NucRemovalE", fNucRmvE ); // GeV
1398  GetParam( "INUKE-HadStep", fHadStep ) ;
1399  GetParam( "INUKE-NucAbsFac", fNucAbsFac ) ;
1400  GetParam( "INUKE-NucCEXFac", fNucCEXFac ) ;
1401  GetParam( "INUKE-Energy_Pre_Eq", fEPreEq ) ;
1402  GetParam( "INUKE-FermiFac", fFermiFac ) ;
1403  GetParam( "INUKE-FermiMomentum", fFermiMomentum ) ;
1404 
1405  GetParam( "INUKE-DoCompoundNucleus", fDoCompoundNucleus ) ;
1406  GetParam( "INUKE-DoFermi", fDoFermi ) ;
1407 
1408  GetParam( "HAINUKE-DelRPion", fDelRPion ) ;
1409  GetParam( "HAINUKE-DelRNucleon", fDelRNucleon ) ;
1410 
1411  GetParamDef( "FSI-ChargedPion-MFPScale", fChPionMFPScale, 1.0 ) ;
1412  GetParamDef( "FSI-NeutralPion-MFPScale", fNeutralPionMFPScale, 1.0 ) ;
1413  GetParamDef( "FSI-Pion-FracCExScale", fPionFracCExScale, 1.0 ) ;
1414  GetParamDef( "FSI-Pion-FracAbsScale", fPionFracAbsScale, 1.0 ) ;
1415  GetParamDef( "FSI-Pion-FracElasScale", fPionFracElasScale, 1.0 ) ;
1416  GetParamDef( "FSI-Pion-FracInelScale", fPionFracInelScale, 1.0 ) ;
1417  GetParamDef( "FSI-Pion-FracPiProdScale", fPionFracPiProdScale, 1.0 ) ;
1418  GetParamDef( "FSI-Nucleon-MFPScale", fNucleonMFPScale, 1.0 ) ;
1419  GetParamDef( "FSI-Nucleon-FracCExScale", fNucleonFracCExScale , 1.0 ) ;
1420  GetParamDef( "FSI-Nucleon-FracElasScale", fNucleonFracElasScale, 1.0 ) ;
1421  GetParamDef( "FSI-Nucleon-FracInelScale", fNucleonFracInelScale, 1.0 ) ;
1422  GetParamDef( "FSI-Nucleon-FracAbsScale", fNucleonFracAbsScale, 1.0 ) ;
1423  GetParamDef( "FSI-Nucleon-FracPiProdScale", fNucleonFracPiProdScale, 1.0 ) ;
1424 
1425  // report
1426  LOG("HAIntranuke", pINFO) << "Settings for INTRANUKE mode: " << INukeMode::AsString(kIMdHA);
1427  LOG("HAIntranuke", pINFO) << "R0 = " << fR0 << " fermi";
1428  LOG("HAIntranuke", pINFO) << "NR = " << fNR;
1429  LOG("HAIntranuke", pINFO) << "DelRPion = " << fDelRPion;
1430  LOG("HAIntranuke", pINFO) << "DelRNucleon = " << fDelRNucleon;
1431  LOG("HAIntranuke", pINFO) << "HadStep = " << fHadStep << " fermi";
1432  LOG("HAIntranuke", pINFO) << "EPreEq = " << fHadStep << " fermi";
1433  LOG("HAIntranuke", pINFO) << "NucAbsFac = " << fNucAbsFac;
1434  LOG("HAIntranuke", pINFO) << "NucCEXFac = " << fNucCEXFac;
1435  LOG("HAIntranuke", pINFO) << "FermiFac = " << fFermiFac;
1436  LOG("HAIntranuke", pINFO) << "FermiMomtm = " << fFermiMomentum;
1437  LOG("HAIntranuke", pINFO) << "DoFermi? = " << ((fDoFermi)?(true):(false));
1438  LOG("HAIntranuke", pINFO) << "DoCmpndNuc? = " << ((fDoCompoundNucleus)?(true):(false));
1439 }
genie::Intranuke::fFermiMomentum
double fFermiMomentum
whether or not particle collision is pauli blocked
Definition: Intranuke.h:112
genie::INukeMode::AsString
static string AsString(INukeMode_t mode)
Definition: INukeMode.h:41
genie::Intranuke::fNucleonFracInelScale
double fNucleonFracInelScale
Definition: Intranuke.h:127
genie::Intranuke::fHadroData
INukeHadroData * fHadroData
a collection of h+N,h+A data &amp; calculations
Definition: Intranuke.h:93
genie::Intranuke::fHadStep
double fHadStep
step size for intranuclear hadron transport
Definition: Intranuke.h:107
genie::Intranuke::fNucleonFracCExScale
double fNucleonFracCExScale
Definition: Intranuke.h:125
genie::NuclearModelI
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
genie::Intranuke::fNucleonFracAbsScale
double fNucleonFracAbsScale
Definition: Intranuke.h:128
genie::Intranuke::fDelRNucleon
double fDelRNucleon
factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke.h:106
genie::Intranuke::fNucleonMFPScale
double fNucleonMFPScale
Definition: Intranuke.h:124
genie::Intranuke::fDelRPion
double fDelRPion
factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke.h:105
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::Intranuke::fPionFracPiProdScale
double fPionFracPiProdScale
Definition: Intranuke.h:123
genie::Intranuke::fPionFracInelScale
double fPionFracInelScale
Definition: Intranuke.h:121
genie::INukeHadroData::Instance
static INukeHadroData * Instance(void)
Definition: INukeHadroData.cxx:172
genie::Intranuke::fDoCompoundNucleus
bool fDoCompoundNucleus
whether or not to do compound nucleus considerations
Definition: Intranuke.h:115
genie::Intranuke::fR0
double fR0
effective nuclear size param
Definition: Intranuke.h:102
pINFO
#define pINFO
Definition: Messenger.h:62
genie::Intranuke::fNucRmvE
double fNucRmvE
binding energy to subtract from cascade nucleons
Definition: Intranuke.h:104
genie::Intranuke::fNucAbsFac
double fNucAbsFac
absorption xsec correction factor (hN Mode)
Definition: Intranuke.h:108
genie::Intranuke::fEPreEq
double fEPreEq
threshold for pre-equilibrium reaction
Definition: Intranuke.h:110
genie::Intranuke::fPionFracAbsScale
double fPionFracAbsScale
Definition: Intranuke.h:122
genie::Intranuke::fPionFracCExScale
double fPionFracCExScale
Definition: Intranuke.h:119
genie::Intranuke::fPionFracElasScale
double fPionFracElasScale
Definition: Intranuke.h:120
genie::Intranuke::fNucCEXFac
double fNucCEXFac
charge exchange xsec correction factor (hN Mode)
Definition: Intranuke.h:109
genie::Intranuke::fChPionMFPScale
double fChPionMFPScale
tweaking factors for tuning
Definition: Intranuke.h:117
genie::Intranuke::fDoFermi
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke.h:113
genie::Intranuke::fFermiFac
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke.h:111
genie::Algorithm::GetParamDef
bool GetParamDef(const RgKey &name, T &p, const T &def) const
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::Intranuke::fNucleonFracPiProdScale
double fNucleonFracPiProdScale
Definition: Intranuke.h:129
genie::Intranuke::fNeutralPionMFPScale
double fNeutralPionMFPScale
Definition: Intranuke.h:118
genie::Intranuke::fNucleonFracElasScale
double fNucleonFracElasScale
Definition: Intranuke.h:126
genie::Intranuke::fNuclmodel
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke.h:95
genie::Intranuke::fNR
double fNR
param multiplying the nuclear radius, determining how far to track hadrons beyond the &quot;nuclear bounda...
Definition: Intranuke.h:103
genie::Algorithm::SubAlg
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
double HAIntranuke::PiBounce ( void  ) const
private

Definition at line 363 of file HAIntranuke.cxx.

References genie::RandomGen::Instance(), LOG, pNOTICE, and genie::RandomGen::RndFsi().

Referenced by ElasHA().

364 {
365 // [adapted from neugen3 intranuke_bounce.F]
366 // [is a fortran stub / difficult to understand - needs to be improved]
367 //
368 // Generates theta in radians for elastic pion-nucleus scattering/
369 // Lookup table is based on Fig 17 of Freedman, Miller and Henley, Nucl.Phys.
370 // A389, 457 (1982)
371 //
372  const int nprob = 25;
373  double dintor = 0.0174533;
374  double denom = 47979.453;
375  double rprob[nprob] = {
376  5000., 4200., 3000., 2600., 2100., 1800., 1200., 750., 500., 230., 120.,
377  35., 9., 3., 11., 18., 29., 27., 20., 14., 10., 6., 2., 0.14, 0.19 };
378 
379  double angles[nprob];
380  for(int i=0; i<nprob; i++) angles[i] = 2.5*i;
381 
382  RandomGen * rnd = RandomGen::Instance();
383  double r = rnd->RndFsi().Rndm();
384 
385  double xsum = 0.;
386  double theta = 0.;
387  double binl = 0.;
388  double binh = 0.;
389  int tj = 0;
390  for(int i=0; i<60; i++) {
391  theta = i+0.5;
392  for(int j=0; j < nprob-1; j++) {
393  binl = angles[j];
394  binh = angles[j+1];
395  tj=j;
396  if(binl<=theta && binh>=theta) break;
397  tj=0;
398  }//j
399  int itj = tj;
400  double tfract = (theta-binl)/2.5;
401  double delp = rprob[itj+1] - rprob[itj];
402  xsum += (rprob[itj] + tfract*delp)/denom;
403  if(xsum>r) break;
404  theta = 0.;
405  }//i
406 
407  theta *= dintor;
408 
409  LOG("HAIntranuke", pNOTICE)
410  << "Generated pi+A elastic scattering angle = " << theta << " radians";
411 
412  return theta;
413 }
genie::RandomGen::RndFsi
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:59
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
double HAIntranuke::PnBounce ( void  ) const
private

Definition at line 415 of file HAIntranuke.cxx.

References genie::RandomGen::Instance(), LOG, pNOTICE, and genie::RandomGen::RndFsi().

Referenced by ElasHA().

416 {
417 // [adapted from neugen3 intranuke_pnbounce.F]
418 // [is a fortran stub / difficult to understand - needs to be improved]
419 //
420 // Generates theta in radians for elastic nucleon-nucleus scattering.
421 // Use 800 MeV p+O16 as template in same (highly simplified) spirit as pi+A
422 // from table in Adams et al., PRL 1979. Guess value at 0-2 deg based on Ni
423 // data.
424 //
425  const int nprob = 20;
426  double dintor = 0.0174533;
427  double denom = 11967.0;
428  double rprob[nprob] = {
429  2400., 2350., 2200., 2000., 1728., 1261., 713., 312., 106., 35.,
430  6., 5., 10., 12., 11., 9., 6., 1., 1., 1. };
431 
432  double angles[nprob];
433  for(int i=0; i<nprob; i++) angles[i] = 1.0*i;
434 
435  RandomGen * rnd = RandomGen::Instance();
436  double r = rnd->RndFsi().Rndm();
437 
438  double xsum = 0.;
439  double theta = 0.;
440  double binl = 0.;
441  double binh = 0.;
442  int tj = 0;
443  for(int i=0; i< nprob; i++) {
444  theta = i+0.5;
445  for(int j=0; j < nprob-1; j++) {
446  binl = angles[j];
447  binh = angles[j+1];
448  tj=j;
449  if(binl<=theta && binh>=theta) break;
450  tj=0;
451  }//j
452  int itj = tj;
453  double tfract = (theta-binl)/2.5;
454  double delp = rprob[itj+1] - rprob[itj];
455  xsum += (rprob[itj] + tfract*delp)/denom;
456  if(xsum>r) break;
457  theta = 0.;
458  }//i
459 
460  theta *= dintor;
461 
462  LOG("HAIntranuke", pNOTICE)
463  << "Generated N+A elastic scattering angle = " << theta << " radians";
464 
465  return theta;
466 }
genie::RandomGen::RndFsi
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:59
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
void HAIntranuke::ProcessEventRecord ( GHepRecord *  event_rec) const
virtual

Reimplemented from genie::Intranuke.

Definition at line 108 of file HAIntranuke.cxx.

References LOG, pINFO, pNOTICE, and genie::Intranuke::ProcessEventRecord().

109 {
110  LOG("HAIntranuke", pNOTICE)
111  << "************ Running HA MODE INTRANUKE ************";
112 
113  Intranuke::ProcessEventRecord(evrec);
114 
115  LOG("HAIntranuke", pINFO) << "Done with this event";
116 }
genie::Intranuke::ProcessEventRecord
virtual void ProcessEventRecord(GHepRecord *event_rec) const
Definition: Intranuke.cxx:114
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
pINFO
#define pINFO
Definition: Messenger.h:62
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
void HAIntranuke::SimulateHadronicFinalState ( GHepRecord *  ev,
GHepParticle *  p 
) const
privatevirtual

Implements genie::Intranuke.

Definition at line 118 of file HAIntranuke.cxx.

References genie::GHepRecord::AddParticle(), genie::INukeHadroFates::AsString(), genie::GHepParticle::FirstMother(), fNumIterations, HadronFateHA(), genie::kIHAFtUndefined, genie::kIStStableFinalState, genie::kPdgGamma, genie::kPdgKM, genie::kPdgKP, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, LOG, genie::GHepParticle::Name(), genie::GHepRecord::Particle(), genie::GHepParticle::Pdg(), pERROR, pNOTICE, genie::GHepParticle::SetStatus(), and SimulateHadronicFinalStateKinematics().

Referenced by SimulateHadronicFinalStateKinematics().

120 {
121 // Simulate a hadron interaction for the input particle p in HA mode
122 //
123  // check inputs
124  if(!p || !ev) {
125  LOG("HAIntranuke", pERROR) << "** Null input!";
126  return;
127  }
128 
129  // get particle code and check whether this particle can be handled
130  int pdgc = p->Pdg();
131  bool is_gamma = (pdgc==kPdgGamma);
132  bool is_pion = (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0);
133  bool is_kaon = (pdgc==kPdgKP || pdgc==kPdgKM);
134  bool is_baryon = (pdgc==kPdgProton || pdgc==kPdgNeutron);
135  bool is_handled = (is_baryon || is_pion || is_kaon || is_gamma);
136  if(!is_handled) {
137  LOG("HAIntranuke", pERROR) << "** Can not handle particle: " << p->Name();
138  return;
139  }
140 
141  // select a fate for the input particle
142  INukeFateHA_t fate = this->HadronFateHA(p);
143 
144  // store the fate
145  ev->Particle(p->FirstMother())->SetRescatterCode((int)fate);
146 
147  if(fate == kIHAFtUndefined) {
148  LOG("HAIntranuke", pERROR) << "** Couldn't select a fate";
149  p->SetStatus(kIStStableFinalState);
150  ev->AddParticle(*p);
151  return;
152  }
153  LOG("HAIntranuke", pNOTICE)
154  << "Selected "<< p->Name() << " fate: "<< INukeHadroFates::AsString(fate);
155 
156  // try to generate kinematics - repeat till is done (should seldom need >2)
157  fNumIterations = 0;
158  this->SimulateHadronicFinalStateKinematics(ev,p);
159 }
genie::GHepRecord::Particle
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
pERROR
#define pERROR
Definition: Messenger.h:59
genie::HAIntranuke::fNumIterations
unsigned int fNumIterations
Definition: HAIntranuke.h:79
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:63
genie::GHepParticle::FirstMother
int FirstMother(void) const
Definition: GHepParticle.h:66
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:135
genie::HAIntranuke::SimulateHadronicFinalStateKinematics
void SimulateHadronicFinalStateKinematics(GHepRecord *ev, GHepParticle *p) const
Definition: HAIntranuke.cxx:161
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:173
genie::kPdgGamma
const int kPdgGamma
Definition: PDGCodes.h:189
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::HAIntranuke::HadronFateHA
INukeFateHA_t HadronFateHA(const GHepParticle *p) const
Definition: HAIntranuke.cxx:212
genie::GHepParticle::SetStatus
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:126
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
genie::GHepRecord::AddParticle
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
genie::INukeFateHA_t
enum genie::EINukeFateHA_t INukeFateHA_t
void HAIntranuke::SimulateHadronicFinalStateKinematics ( GHepRecord *  ev,
GHepParticle *  p 
) const
private

Definition at line 161 of file HAIntranuke.cxx.

References genie::INukeHadroFates::AsString(), ElasHA(), genie::GHepParticle::FirstMother(), fNumIterations, Inelastic(), InelasticHA(), genie::kIHAFtAbs, genie::kIHAFtCEx, genie::kIHAFtElas, genie::kIHAFtInelas, genie::kIHAFtPiProd, LOG, genie::GHepParticle::Name(), genie::GHepRecord::Particle(), pINFO, pNOTICE, and SimulateHadronicFinalState().

Referenced by SimulateHadronicFinalState().

163 {
164  // get stored fate
165  INukeFateHA_t fate = (INukeFateHA_t)
166  ev->Particle(p->FirstMother())->RescatterCode();
167 
168  LOG("HAIntranuke", pINFO)
169  << "Generating kinematics for " << p->Name()
170  << " fate: "<< INukeHadroFates::AsString(fate);
171 
172  // try to generate kinematics for the selected fate
173 
174  try
175  {
176  fNumIterations++;
177  if (fate == kIHAFtElas)
178  {
179  this->ElasHA(ev,p,fate);
180  }
181  else
182  if (fate == kIHAFtInelas || fate == kIHAFtCEx)
183  {
184  this->InelasticHA(ev,p,fate);
185  }
186  else if (fate == kIHAFtAbs || fate == kIHAFtPiProd)
187  {
188  this->Inelastic(ev,p,fate);
189  }
190  }
191  catch(exceptions::INukeException exception)
192  {
193  LOG("HAIntranuke", pNOTICE)
194  << exception;
195  if(fNumIterations <= 100) {
196  LOG("HAIntranuke", pNOTICE)
197  << "Failed attempt to generate kinematics for "
198  << p->Name() << " fate: " << INukeHadroFates::AsString(fate)
199  << " - After " << fNumIterations << " tries, still retrying...";
200  this->SimulateHadronicFinalStateKinematics(ev,p);
201  } else {
202  LOG("HAIntranuke", pNOTICE)
203  << "Failed attempt to generate kinematics for "
204  << p->Name() << " fate: " << INukeHadroFates::AsString(fate)
205  << " after " << fNumIterations-1
206  << " attempts. Trying a new fate...";
207  this->SimulateHadronicFinalState(ev,p);
208  }
209  }
210 }
genie::exceptions::INukeException
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
Definition: INukeException.h:38
genie::GHepRecord::Particle
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
genie::HAIntranuke::fNumIterations
unsigned int fNumIterations
Definition: HAIntranuke.h:79
genie::HAIntranuke::InelasticHA
void InelasticHA(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition: HAIntranuke.cxx:542
genie::GHepParticle::FirstMother
int FirstMother(void) const
Definition: GHepParticle.h:66
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:135
genie::HAIntranuke::SimulateHadronicFinalStateKinematics
void SimulateHadronicFinalStateKinematics(GHepRecord *ev, GHepParticle *p) const
Definition: HAIntranuke.cxx:161
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::HAIntranuke::Inelastic
void Inelastic(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition: HAIntranuke.cxx:704
pINFO
#define pINFO
Definition: Messenger.h:62
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::HAIntranuke::ElasHA
void ElasHA(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition: HAIntranuke.cxx:468
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::INukeFateHA_t
enum genie::EINukeFateHA_t INukeFateHA_t
genie::HAIntranuke::SimulateHadronicFinalState
void SimulateHadronicFinalState(GHepRecord *ev, GHepParticle *p) const
Definition: HAIntranuke.cxx:118

Friends And Related Function Documentation

friend class IntranukeTester
friend

Definition at line 53 of file HAIntranuke.h.

Member Data Documentation

unsigned int genie::HAIntranuke::fNumIterations
mutableprivate

Definition at line 79 of file HAIntranuke.h.

Referenced by SimulateHadronicFinalState(), and SimulateHadronicFinalStateKinematics().

The documentation for this class was generated from the following files:
Generated on Wed Jul 17 2024 14:46:58 for GENIEGenerator by   doxygen 1.8.5