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

#include <HAIntranuke2018.h>

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

Public Member Functions

 HAIntranuke2018 ()
 
 HAIntranuke2018 (string config)
 
 ~HAIntranuke2018 ()
 
void ProcessEventRecord (GHepRecord *event_rec) const
 
virtual string GetINukeMode () const
 
virtual string GetGenINukeMode () const
 
- Public Member Functions inherited from genie::Intranuke2018
 Intranuke2018 ()
 
 Intranuke2018 (string name)
 
 Intranuke2018 (string name, string config)
 
 ~Intranuke2018 ()
 
virtual void Configure (const Registry &config)
 
virtual void Configure (string param_set)
 
void SetRemnA (int A)
 
void SetRemnZ (int Z)
 
double GetRemnA () const
 
double GetRemnZ () const
 
double GetR0 () const
 
double GetNR () const
 
double GetDelRPion () const
 
double GetDelRNucleon () const
 
double GetNucRmvE () const
 
double GetHadStep () const
 
bool GetUseOset () const
 
bool GetAltOset () const
 
bool GetXsecNNCorr () const
 
- 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
 
int HandleCompoundNucleus (GHepRecord *ev, GHepParticle *p, int mom) const
 

Private Attributes

int nuclA
 value of A for the target nucleus in hA mode More...
 
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::Intranuke2018
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::Intranuke2018
double fTrackingRadius
 tracking radius for the nucleus in the current event More...
 
TGenPhaseSpace fGenPhaseSpace
 a phase space generator More...
 
INukeHadroData2018 * fHadroData2018
 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...
 
bool fUseOset
 Oset model for low energy pion in hN. More...
 
bool fAltOset
 NuWro's table-based implementation (not recommended) More...
 
bool fXsecNNCorr
 use nuclear medium correction for NN cross section More...
 
double fChPionMFPScale
 tweaking factors for tuning More...
 
double fNeutralPionMFPScale
 
double fPionFracCExScale
 
double fPionFracInelScale
 
double fChPionFracAbsScale
 
double fNeutralPionFracAbsScale
 
double fPionFracPiProdScale
 
double fNucleonMFPScale
 
double fNucleonFracCExScale
 
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 HAIntranuke2018.h.

Constructor & Destructor Documentation

HAIntranuke2018::HAIntranuke2018 ( )

Definition at line 96 of file HAIntranuke2018.cxx.

96  :
97 Intranuke2018("genie::HAIntranuke2018")
98 {
99 
100 }
HAIntranuke2018::HAIntranuke2018 ( string  config)

Definition at line 102 of file HAIntranuke2018.cxx.

102  :
103 Intranuke2018("genie::HAIntranuke2018",config)
104 {
105 
106 }
HAIntranuke2018::~HAIntranuke2018 ( )

Definition at line 108 of file HAIntranuke2018.cxx.

109 {
110 
111 }

Member Function Documentation

void HAIntranuke2018::ElasHA ( GHepRecord *  ev,
GHepParticle *  p,
INukeFateHA_t  fate 
) const
private

Definition at line 481 of file HAIntranuke2018.cxx.

References genie::GHepParticle::A(), genie::GHepRecord::AddParticle(), genie::INukeHadroFates::AsString(), genie::PDGLibrary::Find(), genie::Intranuke2018::fRemnA, genie::Intranuke2018::fRemnP4, genie::Intranuke2018::fRemnZ, genie::PDGLibrary::Instance(), 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::intranuke2018::TwoBodyKinematics().

483 {
484  // scatters particle within nucleus, copy of hN code meant to run only once
485  // in hA mode
486 
487  LOG("HAIntranuke2018", pDEBUG)
488  << "ElasHA() is invoked for a : " << p->Name()
489  << " whose fate is : " << INukeHadroFates::AsString(fate);
490 
491  /* if(fate!=kIHAFtElas)
492  {
493  LOG("HAIntranuke2018", pWARN)
494  << "ElasHA() cannot handle fate: " << INukeHadroFates::AsString(fate);
495  return;
496  } */
497 
498  // check remnants
499  if(fRemnA<0 || fRemnZ<0) // best to stop it here and not try again.
500  {
501  LOG("HAIntranuke2018", pWARN) << "Invalid Nucleus! : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
502  p->SetStatus(kIStStableFinalState);
503  ev->AddParticle(*p);
504  return;
505  }
506 
507  // vars for incoming particle, target, and scattered pdg codes
508  int pcode = p->Pdg();
509  double Mp = p->Mass();
510  double Mt = 0.;
511  if (ev->TargetNucleus()->A()==fRemnA)
512  { Mt = PDGLibrary::Instance()->Find(ev->TargetNucleus()->Pdg())->Mass(); }
513  else
514  {
515  Mt = fRemnP4.M();
516  }
517  TLorentzVector t4PpL = *p->P4();
518  TLorentzVector t4PtL = fRemnP4;
519  double C3CM = 0.0;
520 
521  // calculate scattering angle
522  if(pcode==kPdgNeutron||pcode==kPdgProton) C3CM = TMath::Cos(this->PnBounce());
523  else C3CM = TMath::Cos(this->PiBounce());
524 
525  // calculate final 4 momentum of probe
526  TLorentzVector t4P3L, t4P4L;
527 
528  if (!utils::intranuke2018::TwoBodyKinematics(Mp,Mt,t4PpL,t4PtL,t4P3L,t4P4L,C3CM,fRemnP4))
529  {
530  LOG("HAIntranuke2018", pNOTICE) << "ElasHA() failed";
531  exceptions::INukeException exception;
532  exception.SetReason("TwoBodyKinematics failed in ElasHA");
533  throw exception;
534  }
535 
536  // Update probe particle
537  p->SetMomentum(t4P3L);
538  p->SetStatus(kIStStableFinalState);
539 
540  // Update Remnant nucleus
541  fRemnP4 = t4P4L;
542  LOG("HAIntranuke2018",pINFO)
543  << "C3cm = " << C3CM;
544  LOG("HAIntranuke2018",pINFO)
545  << "|p3| = " << t4P3L.Vect().Mag() << ", E3 = " << t4P3L.E() << ",Mp = " << Mp;
546  LOG("HAIntranuke2018",pINFO)
547  << "|p4| = " << fRemnP4.Vect().Mag() << ", E4 = " << fRemnP4.E() << ",Mt = " << Mt;
548 
549  ev->AddParticle(*p);
550 
551 }
genie::exceptions::INukeException
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
Definition: INukeException.h:38
genie::GHepParticle::P4
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
genie::Intranuke2018::fRemnA
int fRemnA
remnant nucleus A
Definition: Intranuke2018.h:119
genie::HAIntranuke2018::PiBounce
double PiBounce(void) const
Definition: HAIntranuke2018.cxx:376
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::intranuke2018::TwoBodyKinematics
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition: INukeUtils2018.cxx:873
genie::GHepRecord::TargetNucleus
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
pINFO
#define pINFO
Definition: Messenger.h:62
pWARN
#define pWARN
Definition: Messenger.h:60
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::HAIntranuke2018::PnBounce
double PnBounce(void) const
Definition: HAIntranuke2018.cxx:428
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::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
genie::Intranuke2018::fRemnZ
int fRemnZ
remnant nucleus Z
Definition: Intranuke2018.h:120
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
genie::Intranuke2018::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke2018.h:121
virtual string genie::HAIntranuke2018::GetGenINukeMode ( ) const
inlinevirtual

Reimplemented from genie::Intranuke2018.

Definition at line 63 of file HAIntranuke2018.h.

63 {return "hA";};
virtual string genie::HAIntranuke2018::GetINukeMode ( ) const
inlinevirtual

Reimplemented from genie::Intranuke2018.

Definition at line 62 of file HAIntranuke2018.h.

62 {return "hA2018";};
INukeFateHA_t HAIntranuke2018::HadronFateHA ( const GHepParticle *  p) const
private

Definition at line 224 of file HAIntranuke2018.cxx.

References genie::INukeHadroFates::AsString(), genie::Intranuke2018::fChPionFracAbsScale, genie::Intranuke2018::fHadroData2018, genie::Intranuke2018::fNeutralPionFracAbsScale, genie::Intranuke2018::fNucleonFracAbsScale, genie::Intranuke2018::fNucleonFracCExScale, genie::Intranuke2018::fNucleonFracInelScale, genie::Intranuke2018::fNucleonFracPiProdScale, genie::Intranuke2018::fPionFracCExScale, genie::Intranuke2018::fPionFracInelScale, genie::Intranuke2018::fPionFracPiProdScale, genie::INukeHadroData2018::FracADep(), genie::INukeHadroData2018::FracAIndep(), genie::RandomGen::Instance(), genie::kIHAFtAbs, genie::kIHAFtCEx, genie::kIHAFtCmp, 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(), nuclA, pDEBUG, genie::GHepParticle::Pdg(), pINFO, pWARN, and genie::RandomGen::RndFsi().

Referenced by SimulateHadronicFinalState().

225 {
226 // Select a hadron fate in HA mode
227 //
228  RandomGen * rnd = RandomGen::Instance();
229 
230  // get pdgc code & kinetic energy in MeV
231  int pdgc = p->Pdg();
232  double ke = p->KinE() / units::MeV;
233 
234  //bool isPion = (pdgc == kPdgPiP or pdgc == kPdgPi0 or pdgc == kPdgPiM);
235  //if (isPion and fUseOset and ke < 350.0) return this->HadronFateOset();
236 
237  LOG("HAIntranuke2018", pINFO)
238  << "Selecting hA fate for " << p->Name() << " with KE = " << ke << " MeV";
239 
240  // try to generate a hadron fate
241  unsigned int iter = 0;
242  while(iter++ < kRjMaxIterations) {
243 
244  // handle pions
245  //
246  if (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0) {
247 
248  double frac_cex = fHadroData2018->FracADep(pdgc, kIHAFtCEx, ke, nuclA);
249  // double frac_elas = fHadroData2018->FracADep(pdgc, kIHAFtElas, ke, nuclA);
250  double frac_inel = fHadroData2018->FracADep(pdgc, kIHAFtInelas, ke, nuclA);
251  double frac_abs = fHadroData2018->FracADep(pdgc, kIHAFtAbs, ke, nuclA);
252  double frac_piprod = fHadroData2018->FracADep(pdgc, kIHAFtPiProd, ke, nuclA);
253  LOG("HAIntranuke2018", pDEBUG)
254  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCEx) << "} = " << frac_cex
255  // << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
256  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
257  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs
258  << "\n frac{" << INukeHadroFates::AsString(kIHAFtPiProd) << "} = " << frac_piprod;
259 
260  // apply external tweaks to fractions
261  frac_cex *= fPionFracCExScale;
262  frac_inel *= fPionFracInelScale;
263  if (pdgc==kPdgPiP || pdgc==kPdgPiM) frac_abs *= fChPionFracAbsScale;
264  if (pdgc==kPdgPi0) frac_abs *= fNeutralPionFracAbsScale;
265  frac_piprod *= fPionFracPiProdScale;
266 
267  double frac_rescale = 1./(frac_cex + frac_inel + frac_abs + frac_piprod);
268 
269  frac_cex *= frac_rescale;
270  frac_inel *= frac_rescale;
271  frac_abs *= frac_rescale;
272  frac_piprod *= frac_rescale;
273 
274 
275  // compute total fraction (can be <1 if fates have been switched off)
276  double tf = frac_cex +
277  // frac_elas +
278  frac_inel +
279  frac_abs +
280  frac_piprod;
281 
282  double r = tf * rnd->RndFsi().Rndm();
283 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
284  LOG("HAIntranuke2018", pDEBUG) << "r = " << r << " (max = " << tf << ")";
285 #endif
286  double cf=0; // current fraction
287  if(r < (cf += frac_cex )) return kIHAFtCEx; // cex
288  // if(r < (cf += frac_elas )) return kIHAFtElas; // elas
289  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
290  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
291  if(r < (cf += frac_piprod )) return kIHAFtPiProd; // pi prod
292 
293  LOG("HAIntranuke2018", pWARN)
294  << "No selection after going through all fates! "
295  << "Total fraction = " << tf << " (r = " << r << ")";
296  }
297 
298  // handle nucleons
299  else if (pdgc==kPdgProton || pdgc==kPdgNeutron) {
300  double frac_cex = fHadroData2018->FracAIndep(pdgc, kIHAFtCEx, ke);
301  //double frac_elas = fHadroData2018->FracAIndep(pdgc, kIHAFtElas, ke);
302  double frac_inel = fHadroData2018->FracAIndep(pdgc, kIHAFtInelas, ke);
303  double frac_abs = fHadroData2018->FracAIndep(pdgc, kIHAFtAbs, ke);
304  double frac_pipro = fHadroData2018->FracAIndep(pdgc, kIHAFtPiProd, ke);
305  double frac_cmp = fHadroData2018->FracAIndep(pdgc, kIHAFtCmp , ke);
306 
307  LOG("HAIntranuke2018", pINFO)
308  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCEx) << "} = " << frac_cex
309  // << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
310  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
311  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs
312  << "\n frac{" << INukeHadroFates::AsString(kIHAFtPiProd) << "} = " << frac_pipro
313  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCmp) << "} = " << frac_cmp; //suarez edit, cmp
314 
315  // apply external tweaks to fractions
316  frac_cex *= fNucleonFracCExScale;
317  frac_inel *= fNucleonFracInelScale;
318  frac_abs *= fNucleonFracAbsScale;
319  frac_pipro *= fNucleonFracPiProdScale;
320 
321  double frac_rescale = 1./(frac_cex + frac_inel + frac_abs + frac_pipro);
322 
323  frac_cex *= frac_rescale;
324  frac_inel *= frac_rescale;
325  frac_abs *= frac_rescale;
326  frac_pipro *= frac_rescale;
327 
328  // compute total fraction (can be <1 if fates have been switched off)
329  double tf = frac_cex +
330  //frac_elas +
331  frac_inel +
332  frac_abs +
333  frac_pipro +
334  frac_cmp; //suarez edit, cmp
335 
336  double r = tf * rnd->RndFsi().Rndm();
337 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
338  LOG("HAIntranuke2018", pDEBUG) << "r = " << r << " (max = " << tf << ")";
339 #endif
340  double cf=0; // current fraction
341  if(r < (cf += frac_cex )) return kIHAFtCEx; // cex
342  //if(r < (cf += frac_elas )) return kIHAFtElas; // elas
343  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
344  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
345  if(r < (cf += frac_pipro )) return kIHAFtPiProd; // pi prod
346  if(r < (cf += frac_cmp )) return kIHAFtCmp; //suarez edit, cmp
347 
348  LOG("HAIntranuke2018", pWARN)
349  << "No selection after going through all fates! "
350  << "Total fraction = " << tf << " (r = " << r << ")";
351  }
352  // handle kaons
353  else if (pdgc==kPdgKP || pdgc==kPdgKM) {
354  double frac_inel = fHadroData2018->FracAIndep(pdgc, kIHAFtInelas, ke);
355  double frac_abs = fHadroData2018->FracAIndep(pdgc, kIHAFtAbs, ke);
356 
357  LOG("HAIntranuke2018", pDEBUG)
358  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
359  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs;
360  // compute total fraction (can be <1 if fates have been switched off)
361  double tf = frac_inel +
362  frac_abs;
363  double r = tf * rnd->RndFsi().Rndm();
364 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
365  LOG("HAIntranuke2018", pDEBUG) << "r = " << r << " (max = " << tf << ")";
366 #endif
367  double cf=0; // current fraction
368  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
369  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
370  }
371  }//iterations
372 
373  return kIHAFtUndefined;
374 }
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::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::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:63
genie::INukeHadroData2018::FracADep
double FracADep(int hpdgc, INukeFateHA_t fate, double ke, int targA) const
Definition: INukeHadroData2018.cxx:1377
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::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:173
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::INukeHadroData2018::FracAIndep
double FracAIndep(int hpdgc, INukeFateHA_t fate, double ke) const
Definition: INukeHadroData2018.cxx:1420
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::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::Intranuke2018::fHadroData2018
INukeHadroData2018 * fHadroData2018
a collection of h+N,h+A data &amp; calculations
Definition: Intranuke2018.h:116
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
genie::HAIntranuke2018::nuclA
int nuclA
value of A for the target nucleus in hA mode
Definition: HAIntranuke2018.h:81
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::controls::kRjMaxIterations
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
int HAIntranuke2018::HandleCompoundNucleus ( GHepRecord *  ev,
GHepParticle *  p,
int  mom 
) const
privatevirtual

Implements genie::Intranuke2018.

Definition at line 1513 of file HAIntranuke2018.cxx.

1514 {
1515 
1516  // only relevant for hN mode - not anymore.
1517  return false;
1518 
1519 }
void HAIntranuke2018::Inelastic ( GHepRecord *  ev,
GHepParticle *  p,
INukeFateHA_t  fate 
) const
private

Definition at line 739 of file HAIntranuke2018.cxx.

References genie::GHepRecord::AddParticle(), genie::INukeHadroFates::AsString(), genie::Intranuke2018::fDoFermi, genie::Intranuke2018::fFermiFac, genie::Intranuke2018::fFermiMomentum, genie::Intranuke2018::fHadroData2018, genie::PDGLibrary::Find(), genie::GHepParticle::FirstMother(), genie::Intranuke2018::fNuclmodel, genie::Intranuke2018::fNucRmvE, genie::Intranuke2018::fRemnA, genie::Intranuke2018::fRemnP4, genie::Intranuke2018::fRemnZ, genie::NuclearModelI::GenerateNucleon(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::INukeHadroData2018::IntBounce(), genie::pdg::IonPdgCode(), 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::intranuke2018::PhaseSpaceDecay(), pINFO, genie::utils::intranuke2018::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::intranuke2018::TwoBodyKinematics(), and genie::GHepParticle::X4().

Referenced by SimulateHadronicFinalStateKinematics().

741 {
742 
743  // Aaron Meyer (05/25/10)
744  //
745  // Called to handle all absorption and pi production reactions
746  //
747  // Nucleons -> Reaction approximated by exponential decay in p+n (sum) space,
748  // gaussian in p-n (difference) space
749  // -fit to hN simulations p C, Fe, Pb at 200, 800 MeV
750  // -get n from isospin, np-nn smaller by 2
751  // Pions -> Reaction approximated with a modified gaussian in p+n space,
752  // normal gaussian in p-n space
753  // -based on fits to multiplicity distributions of hN model
754  // for pi+ C, Fe, Pb at 250, 500 MeV
755  // -fit sum and diff of nn, np to Gaussian
756  // -get pi0 from isospin, np-nn smaller by 2
757  // -get pi- from isospin, np-nn smaller by 4
758  // -add 2-body absorption to better match McKeown data
759  // Kaons -> no guidance, use same code as pions.
760  //
761  // Normally distributed random number generated using Box-Muller transformation
762  //
763  // Pion production reactions rescatter pions in nucleus, otherwise unchanged from
764  // older versions of GENIE
765  //
766 
767 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
768  LOG("HAIntranuke2018", pDEBUG)
769  << "Inelastic() is invoked for a : " << p->Name()
770  << " whose fate is : " << INukeHadroFates::AsString(fate);
771 #endif
772 
773  bool allow_dup = true;
774  PDGCodeList list(allow_dup); // list of final state particles
775 
776  // only absorption/pipro fates allowed
777  if (fate == kIHAFtPiProd) {
778 
779  GHepParticle s1(*p);
780  GHepParticle s2(*p);
781  GHepParticle s3(*p);
782 
783  bool success = utils::intranuke2018::PionProduction(
784  ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4, fDoFermi,fFermiFac,fFermiMomentum,fNuclmodel);
785 
786  if (success){
787  LOG ("HAIntranuke2018",pINFO) << " successful pion production fate";
788  // set status of particles and conserve charge/baryon number
789  s1.SetStatus(kIStStableFinalState); //should be redundant
790  // if (pdg::IsPion(s2->Pdg())) s2->SetStatus(kIStHadronInTheNucleus);
791  s2.SetStatus(kIStStableFinalState);
792  // if (pdg::IsPion(s3->Pdg())) s3->SetStatus(kIStHadronInTheNucleus);
793  s3.SetStatus(kIStStableFinalState);
794 
795  ev->AddParticle(s1);
796  ev->AddParticle(s2);
797  ev->AddParticle(s3);
798 
799  return;
800  }
801  else {
802  LOG("HAIntranuke2018", pNOTICE) << "Error: could not create pion production final state";
803  exceptions::INukeException exception;
804  exception.SetReason("PionProduction kinematics failed - retry kinematics");
805  throw exception;
806  }
807  }
808 
809  else if (fate==kIHAFtAbs)
810 // tuned for pions - mixture of 2-body and many-body
811 // use same for kaons as there is no guidance
812  {
813  // Instances for reference
814  PDGLibrary * pLib = PDGLibrary::Instance();
815  RandomGen * rnd = RandomGen::Instance();
816 
817  double ke = p->KinE() / units::MeV;
818  int pdgc = p->Pdg();
819 
820  if (fRemnA<2)
821  {
822  LOG("HAIntranuke2018", pNOTICE) << "stop propagation - could not create absorption final state: too few particles";
823  p->SetStatus(kIStStableFinalState);
824  ev->AddParticle(*p);
825  return;
826  }
827  if (fRemnZ<1 && (pdgc==kPdgPiM || pdgc==kPdgKM))
828  {
829  LOG("HAIntranuke2018", pNOTICE) << "stop propagation - could not create absorption final state: Pi- or K- cannot be absorbed by only neutrons";
830  p->SetStatus(kIStStableFinalState);
831  ev->AddParticle(*p);
832  return;
833  }
834  if (fRemnA-fRemnZ<1 && (pdgc==kPdgPiP || pdgc==kPdgKP))
835  {
836  LOG("HAIntranuke2018", pINFO) << "stop propagation - could not create absorption final state: Pi+ or K+ cannot be absorbed by only protons";
837  p->SetStatus(kIStStableFinalState);
838  ev->AddParticle(*p);
839  return;
840  }
841 
842  // for now, empirical split between multi-nucleon absorption and pi d -> N N
843  //
844  // added 03/21/11 - Aaron Meyer
845  //
846  if (pdg::IsPion(pdgc) && rnd->RndFsi().Rndm()<1.14*(.903-0.00189*fRemnA)*(1.35-0.00467*ke))
847  { // pi d -> N N, probability determined empirically with McKeown data
848 
849  INukeFateHN_t fate_hN=kIHNFtAbs;
850  int t1code,t2code,scode,s2code;
851  double ppcnt = (double) fRemnZ / (double) fRemnA; // % of protons
852 
853  // choose target nucleon
854  // -- fates weighted by values from Engel, Mosel...
855  if (pdgc==kPdgPiP) {
856  double Prob_pipd_pp=2.*ppcnt*(1.-ppcnt);
857  double Prob_pipnn_pn=.083*(1.-ppcnt)*(1.-ppcnt);
858  if (rnd->RndFsi().Rndm()*(Prob_pipd_pp+Prob_pipnn_pn)<Prob_pipd_pp){
859  t1code=kPdgNeutron; t2code=kPdgProton;
860  scode=kPdgProton; s2code=kPdgProton;}
861  else{
862  t1code=kPdgNeutron; t2code=kPdgNeutron;
863  scode=kPdgProton; s2code=kPdgNeutron;}
864  }
865  else if (pdgc==kPdgPiM) {
866  double Prob_pimd_nn=2.*ppcnt*(1.-ppcnt);
867  double Prob_pimpp_pn=.083*ppcnt*ppcnt;
868  if (rnd->RndFsi().Rndm()*(Prob_pimd_nn+Prob_pimpp_pn)<Prob_pimd_nn){
869  t1code=kPdgProton; t2code=kPdgNeutron;
870  scode=kPdgNeutron; s2code=kPdgNeutron; }
871  else{
872  t1code=kPdgProton; t2code=kPdgProton;
873  scode=kPdgProton; s2code=kPdgNeutron;}
874  }
875  else { // pi0
876  double Prob_pi0d_pn=0.88*ppcnt*(1.-ppcnt); // 2 * .44
877  double Prob_pi0pp_pp=.14*ppcnt*ppcnt;
878  double Prob_pi0nn_nn=.14*(1.-ppcnt)*(1.-ppcnt);
879  double random_number = rnd->RndFsi().Rndm();
880  if (random_number*(Prob_pi0d_pn+Prob_pi0pp_pp+Prob_pi0nn_nn)<Prob_pi0d_pn){
881  t1code=kPdgNeutron; t2code=kPdgProton;
882  scode=kPdgNeutron; s2code=kPdgProton; }
883  else if (random_number*(Prob_pi0d_pn+Prob_pi0pp_pp+Prob_pi0nn_nn)<(Prob_pi0d_pn+Prob_pi0pp_pp)){
884  t1code=kPdgProton; t2code=kPdgProton;
885  scode=kPdgProton; s2code=kPdgProton; }
886  else {
887  t1code=kPdgNeutron; t2code=kPdgNeutron;
888  scode=kPdgNeutron; s2code=kPdgNeutron; }
889  }
890  LOG("HAIntranuke2018",pINFO) << "choose 2 body absorption, probe, fs = " << pdgc <<" "<< scode <<" "<<s2code;
891  // assign proper masses
892  //double M1 = pLib->Find(pdgc) ->Mass();
893  double M2_1 = pLib->Find(t1code)->Mass();
894  double M2_2 = pLib->Find(t2code)->Mass();
895  //double M2 = M2_1 + M2_2;
896  double M3 = pLib->Find(scode) ->Mass();
897  double M4 = pLib->Find(s2code)->Mass();
898 
899  // handle fermi momentum
900  double E2_1L, E2_2L;
901  TVector3 tP2_1L, tP2_2L;
902  //TLorentzVector dNucl_P4;
903  Target target(ev->TargetNucleus()->Pdg());
904  if(fDoFermi)
905  {
906  target.SetHitNucPdg(t1code);
907  fNuclmodel->GenerateNucleon(target);
908  //LOG("HAIntranuke2018", pNOTICE) << "Nuclmodel= " << fNuclmodel->ModelType(target) ;
909  tP2_1L=fFermiFac * fNuclmodel->Momentum3();
910  E2_1L = TMath::Sqrt(tP2_1L.Mag2() + M2_1*M2_1);
911 
912  target.SetHitNucPdg(t2code);
913  fNuclmodel->GenerateNucleon(target);
914  tP2_2L=fFermiFac * fNuclmodel->Momentum3();
915  E2_2L = TMath::Sqrt(tP2_2L.Mag2() + M2_2*M2_2);
916 
917  //dNucl_P4=TLorentzVector(tP2_1L+tP2_2L,E2_1L+E2_2L);
918  }
919  else
920  {
921  tP2_1L.SetXYZ(0.0, 0.0, 0.0);
922  E2_1L = M2_1;
923 
924  tP2_2L.SetXYZ(0.0, 0.0, 0.0);
925  E2_2L = M2_2;
926  }
927  TLorentzVector dNucl_P4=TLorentzVector(tP2_1L+tP2_2L,E2_1L+E2_2L);
928 
929  double E2L = E2_1L + E2_2L;
930 
931  // adjust p to reflect scattering
932  // get random scattering angle
933  double C3CM = fHadroData2018->IntBounce(p,t1code,scode,fate_hN);
934  if (C3CM<-1.)
935  {
936  LOG("HAIntranuke2018", pWARN) << "Inelastic() failed: IntBounce returned bad angle - try again";
937  exceptions::INukeException exception;
938  exception.SetReason("unphysical angle for hN scattering");
939  throw exception;
940  return;
941  }
942 
943  TLorentzVector t4P1L,t4P2L,t4P3L,t4P4L;
944  t4P1L=*p->P4();
945  t4P2L=TLorentzVector(TVector3(tP2_1L+tP2_2L),E2L);
946  double bindE=0.050; // set to fit McKeown data, updated aug 18
947  //double bindE=0.0;
948  if (utils::intranuke2018::TwoBodyKinematics(M3,M4,t4P1L,t4P2L,t4P3L,t4P4L,C3CM,fRemnP4,bindE))
949  {
950  //construct remnant nucleus and its mass
951 
952  if (pdgc==kPdgPiP || pdgc==kPdgKP) fRemnZ++;
953  if (pdgc==kPdgPiM || pdgc==kPdgKM) fRemnZ--;
954  if (t1code==kPdgProton) fRemnZ--;
955  if (t2code==kPdgProton) fRemnZ--;
956  fRemnA-=2;
957 
958  fRemnP4-=dNucl_P4;
959 
960  TParticlePDG * remn = 0;
961  double MassRem = 0.;
962  int ipdgc = pdg::IonPdgCode(fRemnA, fRemnZ);
963  remn = PDGLibrary::Instance()->Find(ipdgc);
964  if(!remn)
965  {
966  LOG("HAIntranuke2018", pINFO)
967  << "NO Particle with [A = " << fRemnA << ", Z = " << fRemnZ
968  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
969  }
970  else
971  {
972  MassRem = remn->Mass();
973  LOG("HAIntranuke2018", pINFO)
974  << "Particle with [A = " << fRemnA << ", Z = " << fRemnZ
975  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
976  }
977  double ERemn = fRemnP4.E();
978  double PRemn = TMath::Sqrt(fRemnP4.Px()*fRemnP4.Px() + fRemnP4.Py()*fRemnP4.Py() + fRemnP4.Pz()*fRemnP4.Pz());
979  double MRemn = TMath::Sqrt(ERemn*ERemn - PRemn*PRemn);
980  LOG("HAIntranuke2018",pINFO) << "PRemn = " << PRemn << " ;ERemn= " << ERemn;
981  LOG("HAIntranuke2018",pINFO) << "expt MRemn= " << MRemn << " ;true Mass= " << MassRem << " ; excitation energy (>0 good)= " << (MRemn-MassRem)*1000. << " MeV";
982 
983  // create t particles w/ appropriate momenta, code, and status
984  // Set target's mom to be the mom of the hadron that was cloned
985  GHepParticle* t1 = new GHepParticle(*p);
986  GHepParticle* t2 = new GHepParticle(*p);
987  t1->SetFirstMother(p->FirstMother());
988  t1->SetLastMother(p->LastMother());
989  t2->SetFirstMother(p->FirstMother());
990  t2->SetLastMother(p->LastMother());
991 
992  // adjust p to reflect scattering
993  t1->SetPdgCode(scode);
994  t1->SetMomentum(t4P3L);
995 
996  t2->SetPdgCode(s2code);
997  t2->SetMomentum(t4P4L);
998 
999  t1->SetStatus(kIStStableFinalState);
1000  t2->SetStatus(kIStStableFinalState);
1001 
1002  ev->AddParticle(*t1);
1003  ev->AddParticle(*t2);
1004  delete t1;
1005  delete t2;
1006 
1007  return;
1008  }
1009  else
1010  {
1011  LOG("HAIntranuke2018", pNOTICE) << "Inelastic in hA failed calling TwoBodyKineamtics";
1012  exceptions::INukeException exception;
1013  exception.SetReason("Pion absorption kinematics through TwoBodyKinematics failed");
1014  throw exception;
1015 
1016  }
1017 
1018  } // end pi d -> N N
1019  else // multi-nucleon
1020  {
1021 
1022  // declare some parameters for double gaussian and determine values chosen
1023  // parameters for proton and pi+, others come from isospin transformations
1024 
1025  double ns0=0; // mean - sum of nucleons
1026  double nd0=0; // mean - difference of nucleons
1027  double Sig_ns=0; // std dev - sum
1028  double Sig_nd=0; // std dev - diff
1029  double gam_ns=0; // exponential decay rate (for nucleons)
1030 
1031  if ( pdg::IsNeutronOrProton (pdgc) ) // nucleon probe
1032  {
1033  // antisymmetric about Z=N
1034  if (fRemnA-fRemnZ > fRemnZ)
1035  nd0 = 135.227 * TMath::Exp(-7.124*(fRemnA-fRemnZ)/double(fRemnA)) - 2.762;
1036  else
1037  nd0 = -135.227 * TMath::Exp(-7.124* fRemnZ /double(fRemnA)) + 4.914;
1038 
1039  Sig_nd = 2.034 + fRemnA * 0.007846;
1040 
1041  double c1 = 0.041 + ke * 0.0001525;
1042  double c2 = -0.003444 - ke * 0.00002324;
1043 //change last factor from 30 to 15 so that gam_ns always larger than 0
1044 //add check to be certain
1045  double c3 = 0.064 - ke * 0.000015;
1046  gam_ns = c1 * TMath::Exp(c2*fRemnA) + c3;
1047  if(gam_ns<0.002) gam_ns = 0.002;
1048  //gam_ns = 10.;
1049  LOG("HAIntranuke2018", pINFO) << "nucleon absorption";
1050  LOG("HAIntranuke2018", pINFO) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
1051  // LOG("HAIntranuke2018", pINFO) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
1052  LOG("HAIntranuke2018", pINFO) << "--> gam_ns = " << gam_ns;
1053  }
1054  else if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) //pion probe
1055  {
1056  ns0 = .0001*(1.+ke/250.) * (fRemnA-10)*(fRemnA-10) + 3.5;
1057  nd0 = (1.+ke/250.) - ((fRemnA/200.)*(1. + 2.*ke/250.));
1058  Sig_ns = (10. + 4. * ke/250.)*TMath::Power(fRemnA/250.,0.9); //(1. - TMath::Exp(-0.02*fRemnA));
1059  Sig_nd = 4*(1 - TMath::Exp(-0.03*ke));
1060  LOG("HAIntranuke2018", pINFO) << "pion absorption";
1061  LOG("HAIntranuke2018", pINFO) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
1062  LOG("HAIntranuke2018", pINFO) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
1063  }
1064  else if (pdgc==kPdgKP || pdgc==kPdgKM) // kaon probe
1065  {
1066  ns0 = (rnd->RndFsi().Rndm()>0.5?3:2);
1067  nd0 = 1.;
1068  Sig_ns = 0.1;
1069  Sig_nd = 0.1;
1070  LOG("HAIntranuke2018", pINFO) << "kaon absorption - set ns, nd later";
1071  // LOG("HAIntranuke2018", pINFO) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
1072  // LOG("HAIntranuke2018", pINFO) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
1073  }
1074  else
1075  {
1076  LOG("HAIntranuke2018", pWARN) << "Inelastic() cannot handle absorption reaction for " << p->Name();
1077  }
1078 
1079  // account for different isospin
1080  if (pdgc==kPdgPi0 || pdgc==kPdgNeutron) nd0-=2.;
1081  if (pdgc==kPdgPiM) nd0-=4.;
1082 
1083  int iter=0; // counter
1084  int np=0,nn=0; // # of p, # of n
1085  bool not_done=true;
1086  double u1 = 0, u2 = 0;
1087 
1088  while (not_done)
1089  {
1090  // infinite loop check
1091  if (iter>=100) {
1092  LOG("HAIntranuke2018", pNOTICE) << "Error: could not choose absorption final state";
1093  LOG("HAIntranuke2018", pNOTICE) << "--> mean diff distr = " << nd0 << ", stand dev = " << Sig_nd;
1094  LOG("HAIntranuke2018", pNOTICE) << "--> mean sum distr = " << ns0 << ", Stand dev = " << Sig_ns;
1095  LOG("HAIntranuke2018", pNOTICE) << "--> gam_ns = " << gam_ns;
1096  LOG("HAIntranuke2018", pNOTICE) << "--> A = " << fRemnA << ", Z = " << fRemnZ << ", Energy = " << ke;
1097  exceptions::INukeException exception;
1098  exception.SetReason("Absorption choice of # of p,n failed");
1099  throw exception;
1100  }
1101  //here??
1102 
1103  // Box-Muller transform
1104  // Takes two uniform distribution random variables on (0,1]
1105  // Creates two normally distributed random variables on (0,inf)
1106 
1107  u1 = rnd->RndFsi().Rndm(); // uniform random variable 1
1108  u2 = rnd->RndFsi().Rndm(); // " " 2
1109  if (u1==0) u1 = rnd->RndFsi().Rndm();
1110  if (u2==0) u2 = rnd->RndFsi().Rndm(); // Just in case
1111 
1112  // normally distributed random variable
1113  double x2 = TMath::Sqrt(-2*TMath::Log(u1))*TMath::Sin(2*kPi*u2);
1114 
1115  double ns = 0;
1116 
1117  if ( pdg::IsNeutronOrProton (pdgc) ) //nucleon probe
1118  {
1119  ns = -TMath::Log(rnd->RndFsi().Rndm())/gam_ns; // exponential random variable
1120  }
1121  if ( pdg::IsKaon (pdgc) ) //charged kaon probe - either 2 or 3 nucleons to stay simple
1122  {
1123  ns = (rnd->RndFsi().Rndm()<0.5?2:3);
1124  }
1125  else if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) //pion probe
1126  {
1127  // Pion fit for sum takes for xs*exp((xs-x0)^2 / 2*sig_xs0)
1128  // Find random variable by first finding gaussian random variable
1129  // then throwing the value against a linear P.D.F.
1130  //
1131  // max is the maximum value allowed for the random variable (10 std + mean)
1132  // minimum allowed value is 0
1133 
1134  double max = ns0 + Sig_ns * 10;
1135  if(max>fRemnA) max=fRemnA;
1136  double x1 = 0;
1137  bool not_found = true;
1138  int iter2 = 0;
1139 
1140  while (not_found)
1141  {
1142  // infinite loop check
1143  if (iter2>=100)
1144  {
1145  LOG("HAIntranuke2018", pNOTICE) << "Error: stuck in random variable loop for ns";
1146  LOG("HAIntranuke2018", pNOTICE) << "--> mean of sum parent distr = " << ns0 << ", Stand dev = " << Sig_ns;
1147  LOG("HAIntranuke2018", pNOTICE) << "--> A = " << fRemnA << ", Z = " << fRemnZ << ", Energy = " << ke;
1148 
1149  exceptions::INukeException exception;
1150  exception.SetReason("Random number generator for choice of #p,n final state failed - unusual - redo kinematics");
1151  throw exception;
1152  }
1153 
1154  // calculate exponential random variable
1155  u1 = rnd->RndFsi().Rndm();
1156  u2 = rnd->RndFsi().Rndm();
1157  if (u1==0) u1 = rnd->RndFsi().Rndm();
1158  if (u2==0) u2 = rnd->RndFsi().Rndm();
1159  x1 = TMath::Sqrt(-2*TMath::Log(u1))*TMath::Cos(2*kPi*u2);
1160 
1161  ns = ns0 + Sig_ns * x1;
1162  if ( ns>max || ns<0 ) {iter2++; continue;}
1163  else if ( rnd->RndFsi().Rndm() > (ns/max) ) {iter2++; continue;}
1164  else {
1165  // accept this sum value
1166  not_found=false;
1167  }
1168  } //while(not_found)
1169  }//else pion
1170 
1171  double nd = nd0 + Sig_nd * x2; // difference (p-n) for both pion, nucleon probe
1172  if (pdgc==kPdgKP) // special for KP
1173  { if (ns==2) nd=0;
1174  if (ns>2) nd=1; }
1175 
1176  np = int((ns+nd)/2.+.5); // Addition of .5 for rounding correction
1177  nn = int((ns-nd)/2.+.5);
1178 
1179  LOG("HAIntranuke2018", pINFO) << "ns = "<<ns<<", nd = "<<nd<<", np = "<<np<<", nn = "<<nn;
1180  //LOG("HAIntranuke2018", pNOTICE) << "RemA = "<<fRemnA<<", RemZ = "<<fRemnZ<<", probe = "<<pdgc;
1181 
1182  /*if ((ns+nd)/2. < 0 || (ns-nd)/2. < 0) {iter++; continue;}
1183  else */
1184  //check for problems befor executing phase space 'decay'
1185  if (np < 0 || nn < 0 ) {iter++; continue;}
1186  else if (np + nn < 2. ) {iter++; continue;}
1187  else if ((np + nn == 2.) && pdg::IsNeutronOrProton (pdgc)) {iter++; continue;}
1188  else if (np > fRemnZ + ((pdg::IsProton(pdgc) || pdgc==kPdgPiP || pdgc==kPdgKP)?1:0)
1189  - ((pdgc==kPdgPiM || pdgc==kPdgKM)?1:0)) {iter++; continue;}
1190  else if (nn > fRemnA-fRemnZ + ((pdg::IsNeutron(pdgc)||pdgc==kPdgPiM||pdgc==kPdgKM)?1:0)
1191  - ((pdgc==kPdgPiP||pdgc==kPdgKP)?1:0)) {iter++; continue;}
1192  else {
1193  not_done=false; //success
1194  LOG("HAIntranuke2018",pINFO) << "success, iter = " << iter << " np, nn = " << np << " " << nn;
1195  if (np+nn>86) // too many particles, scale down
1196  {
1197  double frac = 85./double(np+nn);
1198  np = int(np*frac);
1199  nn = int(nn*frac);
1200  }
1201 
1202  if ( (np==fRemnZ +((pdg::IsProton (pdgc)||pdgc==kPdgPiP||pdgc==kPdgKP)?1:0)-(pdgc==kPdgPiM||pdgc==kPdgKM?1:0))
1203  &&(nn==fRemnA-fRemnZ+((pdg::IsNeutron(pdgc)||pdgc==kPdgPiM||pdgc==kPdgKM)?1:0)-(pdgc==kPdgPiP||pdgc==kPdgKP?1:0)) )
1204  { // leave at least one nucleon in the nucleus to prevent excess momentum
1205  if (rnd->RndFsi().Rndm()<np/(double)(np+nn)) np--;
1206  else nn--;
1207  }
1208 
1209  LOG("HAIntranuke2018", pNOTICE) << "Final state chosen; # protons : "
1210  << np << ", # neutrons : " << nn;
1211  }
1212  } //while(not_done)
1213 
1214  // change remnants to reflect probe
1215  if ( pdgc==kPdgProton || pdgc==kPdgPiP || pdgc==kPdgKP) fRemnZ++;
1216  if ( pdgc==kPdgPiM || pdgc==kPdgKM) fRemnZ--;
1217  if ( pdg::IsNeutronOrProton (pdgc) ) fRemnA++;
1218 
1219  // PhaseSpaceDecay forbids anything over 18 particles
1220  //
1221  // If over 18 particles, split into 5 groups and perform decay on each group
1222  // Anything over 85 particles reduced to 85 in previous step
1223  //
1224  // 4 of the nucleons are used as "probes" as well as original probe,
1225  // with each one sharing 1/5 of the total incident momentum
1226  //
1227  // Note: choosing 5 groups and distributing momentum evenly was arbitrary
1228  // Needs to be revised later
1229 
1230  if ((np+nn)>18)
1231  {
1232  // code lists
1233  PDGCodeList list0(allow_dup);
1234  PDGCodeList list1(allow_dup);
1235  PDGCodeList list2(allow_dup);
1236  PDGCodeList list3(allow_dup);
1237  PDGCodeList list4(allow_dup);
1238  PDGCodeList* listar[5] = {&list0, &list1, &list2, &list3, &list4};
1239 
1240  //set up HadronClusters
1241  // simple for now, each (of 5) in hadron cluster has 1/5 of mom and KE
1242 
1243  double probM = pLib->Find(pdgc) ->Mass();
1244  probM -= .025; // BE correction
1245  TVector3 pP3 = p->P4()->Vect() * (1./5.);
1246  double probKE = p->P4()->E() -probM;
1247  double clusKE = probKE * (1./5.);
1248  TLorentzVector clusP4(pP3,clusKE); //no mass
1249  LOG("HAIntranuke2018",pINFO) << "probM = " << probM << " ;clusKE= " << clusKE;
1250  TLorentzVector X4(*p->X4());
1251  GHepStatus_t ist = kIStNucleonClusterTarget;
1252 
1253  int mom = p->FirstMother();
1254 
1255  GHepParticle * p0 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1256  GHepParticle * p1 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1257  GHepParticle * p2 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1258  GHepParticle * p3 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1259  GHepParticle * p4 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1260 
1261  // To conserve 4-momenta
1262  // fRemnP4 -= probP4 + protP4*np_p + neutP4*(4-np_p) - *p->P4();
1263  fRemnP4 -= 5.*clusP4 - *p->P4();
1264 
1265  for (int i=0;i<(np+nn);i++)
1266  {
1267  if (i<np)
1268  {
1269  listar[i%5]->push_back(kPdgProton);
1270  fRemnZ--;
1271  }
1272  else listar[i%5]->push_back(kPdgNeutron);
1273  fRemnA--;
1274  }
1275  for (int i=0;i<5;i++)
1276  {
1277  LOG("HAIntranuke2018", pINFO) << "List" << i << " size: " << listar[i]->size();
1278  if (listar[i]->size() <2)
1279  {
1280  exceptions::INukeException exception;
1281  exception.SetReason("too few particles for Phase Space decay - try again");
1282  throw exception;
1283  }
1284  }
1285 
1286  // commented out to better fit with absorption reactions
1287  // Add the fermi energy of the nucleons to the phase space
1288  /*if(fDoFermi)
1289  {
1290  GHepParticle* p_ar[5] = {cl, p1, p2, p3, p4};
1291  for (int i=0;i<5;i++)
1292  {
1293  Target target(ev->TargetNucleus()->Pdg());
1294  TVector3 pBuf = p_ar[i]->P4()->Vect();
1295  double mBuf = p_ar[i]->Mass();
1296  double eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1297  TLorentzVector tSum(pBuf,eBuf);
1298  double mSum = 0.0;
1299  vector<int>::const_iterator pdg_iter;
1300  for(pdg_iter=++(listar[i]->begin());pdg_iter!=listar[i]->end();++pdg_iter)
1301  {
1302  target.SetHitNucPdg(*pdg_iter);
1303  fNuclmodel->GenerateNucleon(target);
1304  mBuf = pLib->Find(*pdg_iter)->Mass();
1305  mSum += mBuf;
1306  pBuf = fFermiFac * fNuclmodel->Momentum3();
1307  eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1308  tSum += TLorentzVector(pBuf,eBuf);
1309  fRemnP4 -= TLorentzVector(pBuf,eBuf-mBuf);
1310  }
1311  TLorentzVector dP4 = tSum + TLorentzVector(TVector3(0,0,0),-mSum);
1312  p_ar[i]->SetMomentum(dP4);
1313  }
1314  }*/
1315 
1316  bool success1 = utils::intranuke2018::PhaseSpaceDecay(ev,p0,*listar[0],fRemnP4,fNucRmvE,kIMdHA);
1317  bool success2 = utils::intranuke2018::PhaseSpaceDecay(ev,p1,*listar[1],fRemnP4,fNucRmvE,kIMdHA);
1318  bool success3 = utils::intranuke2018::PhaseSpaceDecay(ev,p2,*listar[2],fRemnP4,fNucRmvE,kIMdHA);
1319  bool success4 = utils::intranuke2018::PhaseSpaceDecay(ev,p3,*listar[3],fRemnP4,fNucRmvE,kIMdHA);
1320  bool success5 = utils::intranuke2018::PhaseSpaceDecay(ev,p4,*listar[4],fRemnP4,fNucRmvE,kIMdHA);
1321  if(success1 && success2 && success3 && success4 && success5)
1322  {
1323  LOG("HAIntranuke2018", pINFO)<<"Successful many-body absorption - n>=18";
1324  LOG("HAIntranuke2018", pDEBUG) << "Nucleus : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
1325  TParticlePDG * remn = 0;
1326  double MassRem = 0.;
1327  int ipdgc = pdg::IonPdgCode(fRemnA, fRemnZ);
1328  remn = PDGLibrary::Instance()->Find(ipdgc);
1329  if(!remn)
1330  {
1331  LOG("HAIntranuke2018", pINFO)
1332  << "NO Particle with [A = " << fRemnA << ", Z = " << fRemnZ
1333  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
1334  }
1335  else
1336  {
1337  MassRem = remn->Mass();
1338  LOG("HAIntranuke2018", pINFO)
1339  << "Particle with [A = " << fRemnA << ", Z = " << fRemnZ
1340  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
1341  }
1342  double ERemn = fRemnP4.E();
1343  double PRemn = TMath::Sqrt(fRemnP4.Px()*fRemnP4.Px() + fRemnP4.Py()*fRemnP4.Py() + fRemnP4.Pz()*fRemnP4.Pz());
1344  double MRemn = TMath::Sqrt(ERemn*ERemn - PRemn*PRemn);
1345  LOG("HAIntranuke2018",pINFO) << "PRemn = " << PRemn << " ;ERemn= " << ERemn;
1346  LOG("HAIntranuke2018",pINFO) << "MRemn= " << MRemn << " ;true Mass= " << MassRem << " ; excitation energy (>0 good)= " << (MRemn-MassRem)*1000. << " MeV";
1347 
1348  }
1349  else
1350  {
1351  // try to recover
1352  LOG("HAIntranuke2018", pWARN) << "PhaseSpace decay fails for HadrCluster- recovery likely incorrect - rethrow event";
1353  p->SetStatus(kIStStableFinalState);
1354  ev->AddParticle(*p);
1355  fRemnA+=np+nn;
1356  fRemnZ+=np;
1357  if ( pdgc==kPdgProton || pdgc==kPdgPiP ) fRemnZ--;
1358  if ( pdgc==kPdgPiM ) fRemnZ++;
1359  if ( pdg::IsNeutronOrProton (pdgc) ) fRemnA--;
1360  /* exceptions::INukeException exception;
1361  exception.SetReason("Phase space generation of absorption final state failed");
1362  throw exception;
1363  */
1364  }
1365 
1366  // delete cl;
1367  delete p0;
1368  delete p1;
1369  delete p2;
1370  delete p3;
1371  delete p4;
1372 
1373  }
1374  else // less than 18 particles pion
1375  {
1376  if (pdgc==kPdgKP) list.push_back(kPdgKP); //normally conserve strangeness
1377  if (pdgc==kPdgKM) list.push_back(kPdgKM);
1378  /*
1379  TParticlePDG * remn0 = 0;
1380  int ipdgc0 = pdg::IonPdgCode(fRemnA, fRemnZ);
1381  remn0 = PDGLibrary::Instance()->Find(ipdgc0);
1382  double Mass0 = remn0->Mass();
1383  TParticlePDG * remnt = 0;
1384  int ipdgct = pdg::IonPdgCode(fRemnA-(nn+np), fRemnZ-np);
1385  remnt = PDGLibrary::Instance()->Find(ipdgct);
1386  double MassRemt = remnt->Mass();
1387  LOG("HAIntranuke2018",pINFO) << "Mass0 = " << Mass0 << " ;Masst= " << MassRemt << " ; diff/nucleon= "<< (Mass0-MassRemt)/(np+nn);
1388  */
1389  //set up HadronCluster
1390 
1391  double probM = pLib->Find(pdgc) ->Mass();
1392  double probBE = (np+nn)*.005; // BE correction
1393  TVector3 pP3 = p->P4()->Vect();
1394  double probKE = p->P4()->E() - (probM - probBE);
1395  double clusKE = probKE; // + np*0.9383 + nn*.9396;
1396  TLorentzVector clusP4(pP3,clusKE); //no mass is correct
1397  LOG("HAIntranuke2018",pINFO) << "probM = " << probM << " ;clusKE= " << clusKE;
1398  TLorentzVector X4(*p->X4());
1399  GHepStatus_t ist = kIStNucleonClusterTarget;
1400  int mom = p->FirstMother();
1401 
1402  GHepParticle * p0 = new GHepParticle(kPdgCompNuclCluster,ist, mom,-1,-1,-1,clusP4,X4);
1403 
1404  //set up remnant nucleus
1405  fRemnP4 -= clusP4 - *p->P4();
1406 
1407  for (int i=0;i<np;i++)
1408  {
1409  list.push_back(kPdgProton);
1410  fRemnA--;
1411  fRemnZ--;
1412  }
1413  for (int i=0;i<nn;i++)
1414  {
1415  list.push_back(kPdgNeutron);
1416  fRemnA--;
1417  }
1418 
1419  // Library instance for reference
1420  //PDGLibrary * pLib = PDGLibrary::Instance();
1421 
1422  // commented out to better fit with absorption reactions
1423  // Add the fermi energy of the nucleons to the phase space
1424  /*if(fDoFermi)
1425  {
1426  Target target(ev->TargetNucleus()->Pdg());
1427  TVector3 pBuf = p->P4()->Vect();
1428  double mBuf = p->Mass();
1429  double eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1430  TLorentzVector tSum(pBuf,eBuf);
1431  double mSum = 0.0;
1432  vector<int>::const_iterator pdg_iter;
1433  for(pdg_iter=++(list.begin());pdg_iter!=list.end();++pdg_iter)
1434  {
1435  target.SetHitNucPdg(*pdg_iter);
1436  fNuclmodel->GenerateNucleon(target);
1437  mBuf = pLib->Find(*pdg_iter)->Mass();
1438  mSum += mBuf;
1439  pBuf = fFermiFac * fNuclmodel->Momentum3();
1440  eBuf = TMath::Sqrt(pBuf.Mag2() + mBuf*mBuf);
1441  tSum += TLorentzVector(pBuf,eBuf);
1442  fRemnP4 -= TLorentzVector(pBuf,eBuf-mBuf);
1443  }
1444  TLorentzVector dP4 = tSum + TLorentzVector(TVector3(0,0,0),-mSum);
1445  p->SetMomentum(dP4);
1446  }*/
1447 
1448  LOG("HAIntranuke2018", pDEBUG)
1449  << "Remnant nucleus (A,Z) = (" << fRemnA << ", " << fRemnZ << ")";
1450  LOG("HAIntranuke2018", pINFO) << " list size: " << np+nn;
1451  if (np+nn <2)
1452  {
1453  exceptions::INukeException exception;
1454  exception.SetReason("too few particles for Phase Space decay - try again");
1455  throw exception;
1456  }
1457  // GHepParticle * cl = new GHepParticle(*p);
1458  // cl->SetPdgCode(kPdgDecayNuclCluster);
1459  //bool success1 = utils::intranuke2018::PhaseSpaceDecay(ev,p0,*listar[0],fRemnP4,fNucRmvE,kIMdHA);
1460  bool success = utils::intranuke2018::PhaseSpaceDecay(ev,p0,list,fRemnP4,fNucRmvE,kIMdHA);
1461  if (success)
1462  {
1463  LOG ("HAIntranuke2018",pINFO) << "Successful many-body absorption, n<=18";
1464  LOG("HAIntranuke2018", pDEBUG) << "Nucleus : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
1465  TParticlePDG * remn = 0;
1466  double MassRem = 0.;
1467  int ipdgc = pdg::IonPdgCode(fRemnA, fRemnZ);
1468  remn = PDGLibrary::Instance()->Find(ipdgc);
1469  if(!remn)
1470  {
1471  LOG("HAIntranuke2018", pINFO)
1472  << "NO Particle with [A = " << fRemnA << ", Z = " << fRemnZ
1473  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
1474  }
1475  else
1476  {
1477  MassRem = remn->Mass();
1478  LOG("HAIntranuke2018", pINFO)
1479  << "Particle with [A = " << fRemnA << ", Z = " << fRemnZ
1480  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
1481  }
1482  double ERemn = fRemnP4.E();
1483  double PRemn = TMath::Sqrt(fRemnP4.Px()*fRemnP4.Px() + fRemnP4.Py()*fRemnP4.Py() + fRemnP4.Pz()*fRemnP4.Pz());
1484  double MRemn = TMath::Sqrt(ERemn*ERemn - PRemn*PRemn);
1485  LOG("HAIntranuke2018",pINFO) << "PRemn = " << PRemn << " ;ERemn= " << ERemn;
1486  LOG("HAIntranuke2018",pINFO) << "expt MRemn= " << MRemn << " ;true Mass= " << MassRem << " ; excitation energy (>0 good)= " << (MRemn-MassRem)*1000. << " MeV";
1487  }
1488  else {
1489  // recover
1490  p->SetStatus(kIStStableFinalState);
1491  ev->AddParticle(*p);
1492  fRemnA+=np+nn;
1493  fRemnZ+=np;
1494  if ( pdgc==kPdgProton || pdgc==kPdgPiP ) fRemnZ--;
1495  if ( pdgc==kPdgPiM ) fRemnZ++;
1496  if ( pdg::IsNeutronOrProton (pdgc) ) fRemnA--;
1497  exceptions::INukeException exception;
1498  exception.SetReason("Phase space generation of absorption final state failed");
1499  throw exception;
1500  }
1501  delete p0;
1502  }
1503  } // end multi-nucleon FS
1504  }
1505  else // not absorption/pipro
1506  {
1507  LOG("HAIntranuke2018", pWARN)
1508  << "Inelastic() can not handle fate: " << INukeHadroFates::AsString(fate);
1509  return;
1510  }
1511 }
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::Intranuke2018::fFermiMomentum
double fFermiMomentum
whether or not particle collision is pauli blocked
Definition: Intranuke2018.h:135
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::utils::intranuke2018::PhaseSpaceDecay
bool PhaseSpaceDecay(GHepRecord *ev, GHepParticle *p, const PDGCodeList &pdgv, TLorentzVector &RemnP4, double NucRmvE, EINukeMode mode=kIMdHA)
general phase space decay method
Definition: INukeUtils2018.cxx:1678
genie::Intranuke2018::fFermiFac
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke2018.h:134
genie::Intranuke2018::fRemnA
int fRemnA
remnant nucleus A
Definition: Intranuke2018.h:119
genie::utils::intranuke2018::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: INukeUtils2018.cxx:1301
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::Intranuke2018::fNuclmodel
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke2018.h:118
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::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:173
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::intranuke2018::TwoBodyKinematics
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition: INukeUtils2018.cxx:873
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
genie::Intranuke2018::fDoFermi
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke2018.h:136
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::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::Intranuke2018::fHadroData2018
INukeHadroData2018 * fHadroData2018
a collection of h+N,h+A data &amp; calculations
Definition: Intranuke2018.h:116
genie::pdg::IsNeutronOrProton
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:351
genie::pdg::IonPdgCode
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
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
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::Intranuke2018::fNucRmvE
double fNucRmvE
binding energy to subtract from cascade nucleons
Definition: Intranuke2018.h:127
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::INukeHadroData2018::IntBounce
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)
Definition: INukeHadroData2018.cxx:1667
genie::PDGCodeList::push_back
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
genie::GHepStatus_t
enum genie::EGHepStatus GHepStatus_t
genie::Intranuke2018::fRemnZ
int fRemnZ
remnant nucleus Z
Definition: Intranuke2018.h:120
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
genie::GHepParticle::SetPdgCode
void SetPdgCode(int c)
Definition: GHepParticle.cxx:230
genie::Intranuke2018::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke2018.h:121
void HAIntranuke2018::InelasticHA ( GHepRecord *  ev,
GHepParticle *  p,
INukeFateHA_t  fate 
) const
private

Definition at line 553 of file HAIntranuke2018.cxx.

References genie::GHepRecord::AddParticle(), genie::INukeHadroFates::AsString(), genie::Intranuke2018::fDoFermi, genie::Intranuke2018::fFermiFac, genie::Intranuke2018::fHadroData2018, genie::PDGLibrary::Find(), genie::Intranuke2018::fNuclmodel, genie::Intranuke2018::fRemnA, genie::Intranuke2018::fRemnP4, genie::Intranuke2018::fRemnZ, genie::NuclearModelI::GenerateNucleon(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::INukeHadroData2018::IntBounce(), genie::pdg::IonPdgCode(), 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::intranuke2018::TwoBodyCollision().

Referenced by SimulateHadronicFinalStateKinematics().

555 {
556  // charge exch and inelastic - scatters particle within nucleus, hA version
557  // each are treated as quasielastic, particle scatters off single nucleon
558 
559 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
560  LOG("HAIntranuke2018", pDEBUG)
561  << "InelasticHA() is invoked for a : " << p->Name()
562  << " whose fate is : " << INukeHadroFates::AsString(fate);
563 #endif
564  if(ev->Probe() ) {
565  LOG("HAIntranuke2018", pINFO) << " probe KE = " << ev->Probe()->KinE();
566  }
567  if(fate!=kIHAFtCEx && fate!=kIHAFtInelas)
568  {
569  LOG("HAIntranuke2018", pWARN)
570  << "InelasticHA() cannot handle fate: " << INukeHadroFates::AsString(fate);
571  return;
572  }
573 
574  // Random number generator
575  RandomGen * rnd = RandomGen::Instance();
576 
577  // vars for incoming particle, target, and scattered pdg codes
578  int pcode = p->Pdg();
579  int tcode, scode, s2code;
580  double ppcnt = (double) fRemnZ / (double) fRemnA; // % of protons
581 
582  // Select a hadron fate in HN mode
583  INukeFateHN_t h_fate;
584  if (fate == kIHAFtCEx) h_fate = kIHNFtCEx;
585  else h_fate = kIHNFtElas;
586 
587  // Select a target randomly, weighted to #
588  // -- Unless, of course, the fate is CEx,
589  // -- in which case the target may be deterministic
590  // Also assign scattered particle code
591  if(fate==kIHAFtCEx)
592  {
593  if(pcode==kPdgPiP) {tcode = kPdgNeutron; scode = kPdgPi0; s2code = kPdgProton;}
594  else if(pcode==kPdgPiM) {tcode = kPdgProton; scode = kPdgPi0; s2code = kPdgNeutron;}
595  else if(pcode==kPdgPi0)
596  {
597  // for pi0
598  tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton) :(kPdgNeutron);
599  scode = (tcode == kPdgProton) ?(kPdgPiP) :(kPdgPiM);
600  s2code = (tcode == kPdgProton) ?(kPdgNeutron):(kPdgProton);
601  }
602  else if(pcode==kPdgProton) {tcode = kPdgNeutron; scode = kPdgNeutron; s2code = kPdgProton;}
603  else if(pcode==kPdgNeutron){tcode = kPdgProton; scode = kPdgProton; s2code = kPdgNeutron;}
604  else
605  { LOG("HAIntranuke2018", pWARN) << "InelasticHA() cannot handle fate: "
606  << INukeHadroFates::AsString(fate)
607  << " for particle " << p->Name();
608  return;
609  }
610  }
611  else
612  {
613  tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton):(kPdgNeutron);
614  // if(pcode == kPdgKP || pcode == kPdgKM) tcode = kPdgProton;
615  scode = pcode;
616  s2code = tcode;
617  }
618 
619  // check remnants
620  if ( fRemnA < 1 ) //we've blown nucleus apart, no need to retry anything - exit
621  {
622  LOG("HAIntranuke2018",pNOTICE) << "InelasticHA() stops : not enough nucleons";
623  p->SetStatus(kIStStableFinalState);
624  ev->AddParticle(*p);
625  return;
626  }
627  else if ( fRemnZ + (((pcode==kPdgProton)||(pcode==kPdgPiP))?1:0) - (pcode==kPdgPiM?1:0)
628  < ((( scode==kPdgProton)||( scode==kPdgPiP)) ?1:0) - (scode ==kPdgPiM ?1:0)
629  + (((s2code==kPdgProton)||(s2code==kPdgPiP)) ?1:0) - (s2code==kPdgPiM ?1:0) )
630  {
631  LOG("HAIntranuke2018",pWARN) << "InelasticHA() failed : too few protons in nucleus";
632  p->SetStatus(kIStStableFinalState);
633  ev->AddParticle(*p);
634  return; // another extreme case, best strategy is to exit and go to next event
635  }
636 
637  GHepParticle t(*p);
638  t.SetPdgCode(tcode);
639 
640  // set up fermi target
641  Target target(ev->TargetNucleus()->Pdg());
642  double tM = t.Mass();
643 
644  // handle fermi momentum
645  if(fDoFermi)
646  {
647  target.SetHitNucPdg(tcode);
648  fNuclmodel->GenerateNucleon(target);
649  TVector3 tP3 = fFermiFac * fNuclmodel->Momentum3();
650  double tE = TMath::Sqrt(tP3.Mag2()+ tM*tM);
651  t.SetMomentum(TLorentzVector(tP3,tE));
652  }
653  else
654  {
655  t.SetMomentum(0,0,0,tM);
656  }
657 
658  GHepParticle * cl = new GHepParticle(*p); // clone particle, to run IntBounce at proper energy
659  // calculate energy and momentum using invariant mass
660  double pM = p->Mass();
661  double E_p = ((*p->P4() + *t.P4()).Mag2() - tM*tM - pM*pM)/(2.0*tM);
662  double P_p = TMath::Sqrt(E_p*E_p - pM*pM);
663  cl->SetMomentum(TLorentzVector(P_p,0,0,E_p));
664  // momentum doesn't have to be in right direction, only magnitude
665  double C3CM = fHadroData2018->IntBounce(cl,tcode,scode,h_fate);
666  delete cl;
667  if (C3CM<-1.) // hope this doesn't occur too often - unphysical but we just pass it on
668  {
669  LOG("HAIntranuke2018", pWARN) << "unphysical angle chosen in InelasicHA - put particle outside nucleus";
670  p->SetStatus(kIStStableFinalState);
671  ev->AddParticle(*p);
672  return;
673  }
674  double KE1L = p->KinE();
675  double KE2L = t.KinE();
676  LOG("HAIntranuke2018",pINFO)
677  << " KE1L = " << KE1L << " " << KE1L << " KE2L = " << KE2L;
678  GHepParticle cl1(*p);
679  GHepParticle cl2(t);
680  bool success = utils::intranuke2018::TwoBodyCollision(ev,pcode,tcode,scode,s2code,C3CM,
681  &cl1,&cl2,fRemnA,fRemnZ,fRemnP4,kIMdHA);
682  if(success)
683  {
684  double P3L = TMath::Sqrt(cl1.Px()*cl1.Px() + cl1.Py()*cl1.Py() + cl1.Pz()*cl1.Pz());
685  double P4L = TMath::Sqrt(cl2.Px()*cl2.Px() + cl2.Py()*cl2.Py() + cl2.Pz()*cl2.Pz());
686  double E3L = cl1.KinE();
687  double E4L = cl2.KinE();
688  LOG ("HAIntranuke2018",pINFO) << "Successful quasielastic scattering or charge exchange";
689  LOG("HAIntranuke",pINFO)
690  << "C3CM = " << C3CM << "\n P3L, E3L = "
691  << P3L << " " << E3L << " P4L, E4L = "<< P4L << " " << E4L ;
692  if(ev->Probe() ) { LOG("HAIntranuke",pINFO)
693  << "P4L = " << P4L << " ;E4L= " << E4L << "\n probe KE = " << ev->Probe()->KinE() << "\n";
694  LOG("HAIntranuke2018", pINFO) << "Nucleus : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
695  TParticlePDG * remn = 0;
696  double MassRem = 0.;
697  int ipdgc = pdg::IonPdgCode(fRemnA, fRemnZ);
698  remn = PDGLibrary::Instance()->Find(ipdgc);
699  if(!remn)
700  {
701  LOG("HAIntranuke2018", pINFO)
702  << "NO Particle with [A = " << fRemnA << ", Z = " << fRemnZ
703  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
704  }
705  else
706  {
707  MassRem = remn->Mass();
708  LOG("HAIntranuke2018", pINFO)
709  << "Particle with [A = " << fRemnA << ", Z = " << fRemnZ
710  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
711  }
712  double ERemn = fRemnP4.E();
713  double PRemn = TMath::Sqrt(fRemnP4.Px()*fRemnP4.Px() + fRemnP4.Py()*fRemnP4.Py() + fRemnP4.Pz()*fRemnP4.Pz());
714  double MRemn = TMath::Sqrt(ERemn*ERemn - PRemn*PRemn);
715  LOG("HAIntranuke2018",pINFO) << "PRemn = " << PRemn << " ;ERemn= " << ERemn;
716  LOG("HAIntranuke2018",pINFO) << "MRemn= " << MRemn << " ;true Mass= " << MassRem << " ; excitation energy= " << (MRemn-MassRem)*1000. << " MeV";
717  }
718  if (ev->Probe() && (E3L>ev->Probe()->KinE())) //assuming E3 is most important, definitely for pion. what about pp?
719  {
720  // LOG("HAIntranuke",pINFO)
721  // << "E3Lagain = " << E3L << " ;E4L= " << E4L << "\n probe KE = " << ev->Probe()->KinE() << "\n";
722  exceptions::INukeException exception;
723  exception.SetReason("TwoBodyCollison gives KE> probe KE in hA simulation");
724  throw exception;
725  }
726  ev->AddParticle(cl1);
727  ev->AddParticle(cl2);
728 
729  LOG("HAIntranuke2018", pDEBUG) << "Nucleus : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
730  } else
731  {
732  exceptions::INukeException exception;
733  exception.SetReason("TwoBodyCollison failed in hA simulation");
734  throw exception;
735  }
736 
737 }
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::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
genie::GHepParticle::P4
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
genie::Intranuke2018::fFermiFac
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke2018.h:134
genie::Intranuke2018::fRemnA
int fRemnA
remnant nucleus A
Definition: Intranuke2018.h:119
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::Intranuke2018::fNuclmodel
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke2018.h:118
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::utils::intranuke2018::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: INukeUtils2018.cxx:729
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
genie::Intranuke2018::fDoFermi
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke2018.h:136
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::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::Intranuke2018::fHadroData2018
INukeHadroData2018 * fHadroData2018
a collection of h+N,h+A data &amp; calculations
Definition: Intranuke2018.h:116
genie::pdg::IonPdgCode
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
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
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::INukeHadroData2018::IntBounce
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)
Definition: INukeHadroData2018.cxx:1667
genie::Intranuke2018::fRemnZ
int fRemnZ
remnant nucleus Z
Definition: Intranuke2018.h:120
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
genie::Intranuke2018::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke2018.h:121
void HAIntranuke2018::LoadConfig ( void  )
privatevirtual

Implements genie::Intranuke2018.

Definition at line 1521 of file HAIntranuke2018.cxx.

References genie::INukeMode::AsString(), genie::Intranuke2018::fAltOset, genie::Intranuke2018::fChPionFracAbsScale, genie::Intranuke2018::fChPionMFPScale, genie::Intranuke2018::fDelRNucleon, genie::Intranuke2018::fDelRPion, genie::Intranuke2018::fDoCompoundNucleus, genie::Intranuke2018::fDoFermi, genie::Intranuke2018::fEPreEq, genie::Intranuke2018::fFermiFac, genie::Intranuke2018::fFermiMomentum, genie::Intranuke2018::fHadroData2018, genie::Intranuke2018::fHadStep, genie::Intranuke2018::fNeutralPionFracAbsScale, genie::Intranuke2018::fNeutralPionMFPScale, genie::Intranuke2018::fNR, genie::Intranuke2018::fNucAbsFac, genie::Intranuke2018::fNucCEXFac, genie::Intranuke2018::fNucleonFracAbsScale, genie::Intranuke2018::fNucleonFracCExScale, genie::Intranuke2018::fNucleonFracInelScale, genie::Intranuke2018::fNucleonFracPiProdScale, genie::Intranuke2018::fNucleonMFPScale, genie::Intranuke2018::fNuclmodel, genie::Intranuke2018::fNucRmvE, genie::Intranuke2018::fPionFracCExScale, genie::Intranuke2018::fPionFracInelScale, genie::Intranuke2018::fPionFracPiProdScale, genie::Intranuke2018::fR0, genie::Intranuke2018::fUseOset, genie::Intranuke2018::fXsecNNCorr, genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), genie::INukeHadroData2018::Instance(), genie::kIMdHA, LOG, pINFO, and genie::Algorithm::SubAlg().

1522 {
1523 
1524  // load hadronic cross sections
1525  fHadroData2018 = INukeHadroData2018::Instance();
1526 
1527  // fermi momentum setup
1528  // this is specifically set in Intranuke2018::Configure(string)
1529  fNuclmodel = dynamic_cast<const NuclearModelI *>( this -> SubAlg("NuclearModel") ) ;
1530 
1531  // other intranuke config params
1532  GetParam( "NUCL-R0", fR0 ); // fm
1533  GetParam( "NUCL-NR", fNR );
1534 
1535  GetParam( "INUKE-NucRemovalE", fNucRmvE ); // GeV
1536  GetParam( "INUKE-HadStep", fHadStep ) ;
1537  GetParam( "INUKE-NucAbsFac", fNucAbsFac ) ;
1538  GetParam( "INUKE-NucCEXFac", fNucCEXFac ) ;
1539  GetParam( "INUKE-Energy_Pre_Eq", fEPreEq ) ;
1540  GetParam( "INUKE-FermiFac", fFermiFac ) ;
1541  GetParam( "INUKE-FermiMomentum", fFermiMomentum ) ;
1542 
1543  GetParam( "INUKE-DoCompoundNucleus", fDoCompoundNucleus ) ;
1544  GetParam( "INUKE-DoFermi", fDoFermi ) ;
1545  GetParam( "INUKE-XsecNNCorr", fXsecNNCorr ) ;
1546  GetParamDef( "UseOset", fUseOset, false ) ;
1547  GetParamDef( "AltOset", fAltOset, false ) ;
1548 
1549  GetParam( "HAINUKE-DelRPion", fDelRPion ) ;
1550  GetParam( "HAINUKE-DelRNucleon", fDelRNucleon ) ;
1551 
1552  GetParamDef( "FSI-ChargedPion-MFPScale", fChPionMFPScale, 1.0 ) ;
1553  GetParamDef( "FSI-NeutralPion-MFPScale", fNeutralPionMFPScale, 1.0 ) ;
1554  GetParamDef( "FSI-Pion-FracCExScale", fPionFracCExScale, 1.0 ) ;
1555  GetParamDef( "FSI-ChargedPion-FracAbsScale", fChPionFracAbsScale, 1.0 ) ;
1556  GetParamDef( "FSI-NeutralPion-FracAbsScale", fNeutralPionFracAbsScale,1.0 ) ;
1557  GetParamDef( "FSI-Pion-FracInelScale", fPionFracInelScale, 1.0 ) ;
1558  GetParamDef( "FSI-Pion-FracPiProdScale", fPionFracPiProdScale, 1.0 ) ;
1559  GetParamDef( "FSI-Nucleon-MFPScale", fNucleonMFPScale, 1.0 ) ;
1560  GetParamDef( "FSI-Nucleon-FracCExScale", fNucleonFracCExScale , 1.0 ) ;
1561  GetParamDef( "FSI-Nucleon-FracInelScale", fNucleonFracInelScale, 1.0 ) ;
1562  GetParamDef( "FSI-Nucleon-FracAbsScale", fNucleonFracAbsScale, 1.0 ) ;
1563  GetParamDef( "FSI-Nucleon-FracPiProdScale", fNucleonFracPiProdScale, 1.0 ) ;
1564 
1565  // report
1566  LOG("HAIntranuke2018", pINFO) << "Settings for INTRANUKE mode: " << INukeMode::AsString(kIMdHA);
1567  LOG("HAIntranuke2018", pINFO) << "R0 = " << fR0 << " fermi";
1568  LOG("HAIntranuke2018", pINFO) << "NR = " << fNR;
1569  LOG("HAIntranuke2018", pINFO) << "DelRPion = " << fDelRPion;
1570  LOG("HAIntranuke2018", pINFO) << "DelRNucleon = " << fDelRNucleon;
1571  LOG("HAIntranuke2018", pINFO) << "HadStep = " << fHadStep << " fermi";
1572  LOG("HAIntranuke2018", pINFO) << "EPreEq = " << fHadStep << " fermi";
1573  LOG("HAIntranuke2018", pINFO) << "NucAbsFac = " << fNucAbsFac;
1574  LOG("HAIntranuke2018", pINFO) << "NucCEXFac = " << fNucCEXFac;
1575  LOG("HAIntranuke2018", pINFO) << "FermiFac = " << fFermiFac;
1576  LOG("HAIntranuke2018", pINFO) << "FermiMomtm = " << fFermiMomentum;
1577  LOG("HAIntranuke2018", pINFO) << "DoFermi? = " << ((fDoFermi)?(true):(false));
1578  LOG("HAIntranuke2018", pINFO) << "DoCmpndNuc? = " << ((fDoCompoundNucleus)?(true):(false));
1579  LOG("HAIntranuke2018", pINFO) << "XsecNNCorr? = " << ((fXsecNNCorr)?(true):(false));
1580 }
genie::INukeMode::AsString
static string AsString(INukeMode_t mode)
Definition: INukeMode.h:41
genie::Intranuke2018::fFermiMomentum
double fFermiMomentum
whether or not particle collision is pauli blocked
Definition: Intranuke2018.h:135
genie::Intranuke2018::fDelRPion
double fDelRPion
factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke2018.h:128
genie::Intranuke2018::fEPreEq
double fEPreEq
threshold for pre-equilibrium reaction
Definition: Intranuke2018.h:133
genie::Intranuke2018::fNR
double fNR
param multiplying the nuclear radius, determining how far to track hadrons beyond the &quot;nuclear bounda...
Definition: Intranuke2018.h:126
genie::Intranuke2018::fFermiFac
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke2018.h:134
genie::Intranuke2018::fChPionMFPScale
double fChPionMFPScale
tweaking factors for tuning
Definition: Intranuke2018.h:143
genie::NuclearModelI
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
genie::Intranuke2018::fXsecNNCorr
bool fXsecNNCorr
use nuclear medium correction for NN cross section
Definition: Intranuke2018.h:141
genie::Intranuke2018::fNucAbsFac
double fNucAbsFac
absorption xsec correction factor (hN Mode)
Definition: Intranuke2018.h:131
genie::Intranuke2018::fNuclmodel
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke2018.h:118
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::Intranuke2018::fDoFermi
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke2018.h:136
pINFO
#define pINFO
Definition: Messenger.h:62
genie::Intranuke2018::fNucCEXFac
double fNucCEXFac
charge exchange xsec correction factor (hN Mode)
Definition: Intranuke2018.h:132
genie::Intranuke2018::fR0
double fR0
effective nuclear size param
Definition: Intranuke2018.h:125
genie::Intranuke2018::fDoCompoundNucleus
bool fDoCompoundNucleus
whether or not to do compound nucleus considerations
Definition: Intranuke2018.h:138
genie::INukeHadroData2018::Instance
static INukeHadroData2018 * Instance(void)
Definition: INukeHadroData2018.cxx:175
genie::Intranuke2018::fAltOset
bool fAltOset
NuWro&#39;s table-based implementation (not recommended)
Definition: Intranuke2018.h:140
genie::Intranuke2018::fHadroData2018
INukeHadroData2018 * fHadroData2018
a collection of h+N,h+A data &amp; calculations
Definition: Intranuke2018.h:116
genie::Intranuke2018::fDelRNucleon
double fDelRNucleon
factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke2018.h:129
genie::Intranuke2018::fHadStep
double fHadStep
step size for intranuclear hadron transport
Definition: Intranuke2018.h:130
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::Intranuke2018::fNucRmvE
double fNucRmvE
binding energy to subtract from cascade nucleons
Definition: Intranuke2018.h:127
genie::Intranuke2018::fUseOset
bool fUseOset
Oset model for low energy pion in hN.
Definition: Intranuke2018.h:139
genie::Algorithm::SubAlg
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
double HAIntranuke2018::PiBounce ( void  ) const
private

Definition at line 376 of file HAIntranuke2018.cxx.

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

Referenced by ElasHA().

377 {
378 // [adapted from neugen3 intranuke_bounce.F]
379 // [is a fortran stub / difficult to understand - needs to be improved]
380 //
381 // Generates theta in radians for elastic pion-nucleus scattering/
382 // Lookup table is based on Fig 17 of Freedman, Miller and Henley, Nucl.Phys.
383 // A389, 457 (1982)
384 //
385  const int nprob = 25;
386  double dintor = 0.0174533;
387  double denom = 47979.453;
388  double rprob[nprob] = {
389  5000., 4200., 3000., 2600., 2100., 1800., 1200., 750., 500., 230., 120.,
390  35., 9., 3., 11., 18., 29., 27., 20., 14., 10., 6., 2., 0.14, 0.19 };
391 
392  double angles[nprob];
393  for(int i=0; i<nprob; i++) angles[i] = 2.5*i;
394 
395  RandomGen * rnd = RandomGen::Instance();
396  double r = rnd->RndFsi().Rndm();
397 
398  double xsum = 0.;
399  double theta = 0.;
400  double binl = 0.;
401  double binh = 0.;
402  int tj = 0;
403  for(int i=0; i<60; i++) {
404  theta = i+0.5;
405  for(int j=0; j < nprob-1; j++) {
406  binl = angles[j];
407  binh = angles[j+1];
408  tj=j;
409  if(binl<=theta && binh>=theta) break;
410  tj=0;
411  }//j
412  int itj = tj;
413  double tfract = (theta-binl)/2.5;
414  double delp = rprob[itj+1] - rprob[itj];
415  xsum += (rprob[itj] + tfract*delp)/denom;
416  if(xsum>r) break;
417  theta = 0.;
418  }//i
419 
420  theta *= dintor;
421 
422  LOG("HAIntranuke2018", pNOTICE)
423  << "Generated pi+A elastic scattering angle = " << theta << " radians";
424 
425  return theta;
426 }
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 HAIntranuke2018::PnBounce ( void  ) const
private

Definition at line 428 of file HAIntranuke2018.cxx.

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

Referenced by ElasHA().

429 {
430 // [adapted from neugen3 intranuke_pnbounce.F]
431 // [is a fortran stub / difficult to understand - needs to be improved]
432 //
433 // Generates theta in radians for elastic nucleon-nucleus scattering.
434 // Use 800 MeV p+O16 as template in same (highly simplified) spirit as pi+A
435 // from table in Adams et al., PRL 1979. Guess value at 0-2 deg based on Ni
436 // data.
437 //
438  const int nprob = 20;
439  double dintor = 0.0174533;
440  double denom = 11967.0;
441  double rprob[nprob] = {
442  2400., 2350., 2200., 2000., 1728., 1261., 713., 312., 106., 35.,
443  6., 5., 10., 12., 11., 9., 6., 1., 1., 1. };
444 
445  double angles[nprob];
446  for(int i=0; i<nprob; i++) angles[i] = 1.0*i;
447 
448  RandomGen * rnd = RandomGen::Instance();
449  double r = rnd->RndFsi().Rndm();
450 
451  double xsum = 0.;
452  double theta = 0.;
453  double binl = 0.;
454  double binh = 0.;
455  int tj = 0;
456  for(int i=0; i< nprob; i++) {
457  theta = i+0.5;
458  for(int j=0; j < nprob-1; j++) {
459  binl = angles[j];
460  binh = angles[j+1];
461  tj=j;
462  if(binl<=theta && binh>=theta) break;
463  tj=0;
464  }//j
465  int itj = tj;
466  double tfract = (theta-binl)/2.5;
467  double delp = rprob[itj+1] - rprob[itj];
468  xsum += (rprob[itj] + tfract*delp)/denom;
469  if(xsum>r) break;
470  theta = 0.;
471  }//i
472 
473  theta *= dintor;
474 
475  LOG("HAIntranuke2018", pNOTICE)
476  << "Generated N+A elastic scattering angle = " << theta << " radians";
477 
478  return theta;
479 }
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 HAIntranuke2018::ProcessEventRecord ( GHepRecord *  event_rec) const
virtual

Reimplemented from genie::Intranuke2018.

Definition at line 113 of file HAIntranuke2018.cxx.

References genie::units::A, LOG, nuclA, pINFO, pNOTICE, and genie::Intranuke2018::ProcessEventRecord().

114 {
115  LOG("HAIntranuke2018", pNOTICE)
116  << "************ Running hA2018 MODE INTRANUKE ************";
117  GHepParticle * nuclearTarget = evrec -> TargetNucleus();
118  nuclA = nuclearTarget -> A();
119 
120  Intranuke2018::ProcessEventRecord(evrec);
121 
122  LOG("HAIntranuke2018", pINFO) << "Done with this event";
123 }
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::units::A
static constexpr double A
Definition: Units.h:74
pINFO
#define pINFO
Definition: Messenger.h:62
genie::Intranuke2018::ProcessEventRecord
virtual void ProcessEventRecord(GHepRecord *event_rec) const
Definition: Intranuke2018.cxx:119
genie::HAIntranuke2018::nuclA
int nuclA
value of A for the target nucleus in hA mode
Definition: HAIntranuke2018.h:81
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::GHepParticle
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
void HAIntranuke2018::SimulateHadronicFinalState ( GHepRecord *  ev,
GHepParticle *  p 
) const
privatevirtual

Implements genie::Intranuke2018.

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

127 {
128 // Simulate a hadron interaction for the input particle p in HA mode
129 //
130  // check inputs
131  if(!p || !ev) {
132  LOG("HAIntranuke2018", pERROR) << "** Null input!";
133  return;
134  }
135 
136  // get particle code and check whether this particle can be handled
137  int pdgc = p->Pdg();
138  bool is_gamma = (pdgc==kPdgGamma);
139  bool is_pion = (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0);
140  bool is_kaon = (pdgc==kPdgKP || pdgc==kPdgKM);
141  bool is_baryon = (pdgc==kPdgProton || pdgc==kPdgNeutron);
142  bool is_handled = (is_baryon || is_pion || is_kaon || is_gamma);
143  if(!is_handled) {
144  LOG("HAIntranuke2018", pERROR) << "** Can not handle particle: " << p->Name();
145  return;
146  }
147 
148  // select a fate for the input particle
149  INukeFateHA_t fate = this->HadronFateHA(p);
150 
151  // store the fate
152  ev->Particle(p->FirstMother())->SetRescatterCode((int)fate);
153 
154  if(fate == kIHAFtUndefined) {
155  LOG("HAIntranuke2018", pERROR) << "** Couldn't select a fate";
156  p->SetStatus(kIStStableFinalState);
157  ev->AddParticle(*p);
158  return;
159  }
160  LOG("HAIntranuke2018", pNOTICE)
161  << "Selected "<< p->Name() << " fate: "<< INukeHadroFates::AsString(fate);
162 
163  // try to generate kinematics - repeat till is done (should seldom need >2)
164  fNumIterations = 0;
165  this->SimulateHadronicFinalStateKinematics(ev,p);
166 }
genie::GHepRecord::Particle
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
pERROR
#define pERROR
Definition: Messenger.h:59
genie::HAIntranuke2018::fNumIterations
unsigned int fNumIterations
Definition: HAIntranuke2018.h:82
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
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::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::HAIntranuke2018::SimulateHadronicFinalStateKinematics
void SimulateHadronicFinalStateKinematics(GHepRecord *ev, GHepParticle *p) const
Definition: HAIntranuke2018.cxx:168
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::HAIntranuke2018::HadronFateHA
INukeFateHA_t HadronFateHA(const GHepParticle *p) const
Definition: HAIntranuke2018.cxx:224
genie::INukeFateHA_t
enum genie::EINukeFateHA_t INukeFateHA_t
void HAIntranuke2018::SimulateHadronicFinalStateKinematics ( GHepRecord *  ev,
GHepParticle *  p 
) const
private

Definition at line 168 of file HAIntranuke2018.cxx.

References genie::INukeHadroFates::AsString(), genie::Intranuke2018::fDoFermi, genie::Intranuke2018::fFermiFac, genie::GHepParticle::FirstMother(), genie::Intranuke2018::fNuclmodel, genie::Intranuke2018::fNucRmvE, fNumIterations, genie::Intranuke2018::fRemnA, genie::Intranuke2018::fRemnP4, genie::Intranuke2018::fRemnZ, Inelastic(), InelasticHA(), genie::kIHAFtAbs, genie::kIHAFtCEx, genie::kIHAFtCmp, genie::kIHAFtInelas, genie::kIHAFtPiProd, genie::kIMdHA, LOG, genie::GHepParticle::Name(), genie::GHepRecord::Particle(), pINFO, pNOTICE, genie::utils::intranuke2018::PreEquilibrium(), pWARN, and SimulateHadronicFinalState().

Referenced by SimulateHadronicFinalState().

170 {
171  // get stored fate
172  INukeFateHA_t fate = (INukeFateHA_t)
173  ev->Particle(p->FirstMother())->RescatterCode();
174 
175  LOG("HAIntranuke2018", pINFO)
176  << "Generating kinematics for " << p->Name()
177  << " fate: "<< INukeHadroFates::AsString(fate);
178 
179  // try to generate kinematics for the selected fate
180 
181  try
182  {
183  fNumIterations++;
184  /* if (fate == kIHAFtElas)
185  {
186  this->ElasHA(ev,p,fate);
187  }
188  else */
189  if (fate == kIHAFtInelas || fate == kIHAFtCEx)
190  {
191  this->InelasticHA(ev,p,fate);
192  }
193  else if (fate == kIHAFtAbs || fate == kIHAFtPiProd)
194  {
195  this->Inelastic(ev,p,fate);
196  }
197  else if (fate == kIHAFtCmp) //(suarez edit, 17 July, 2017: cmp)
198  {
199  LOG("HAIntranuke2018", pWARN) << "Running PreEquilibrium for kIHAFtCmp";
200  utils::intranuke2018::PreEquilibrium(ev,p,fRemnA,fRemnZ,fRemnP4,fDoFermi,fFermiFac,fNuclmodel,fNucRmvE,kIMdHA); //should be kiMdHA or HN?
201  }
202  }
203  catch(exceptions::INukeException exception)
204  {
205  LOG("HAIntranuke2018", pNOTICE)
206  << exception;
207  if(fNumIterations <= 100) {
208  LOG("HAIntranuke2018", pNOTICE)
209  << "Failed attempt to generate kinematics for "
210  << p->Name() << " fate: " << INukeHadroFates::AsString(fate)
211  << " - After " << fNumIterations << " tries, still retrying...";
212  this->SimulateHadronicFinalStateKinematics(ev,p);
213  } else {
214  LOG("HAIntranuke2018", pNOTICE)
215  << "Failed attempt to generate kinematics for "
216  << p->Name() << " fate: " << INukeHadroFates::AsString(fate)
217  << " after " << fNumIterations-1
218  << " attempts. Trying a new fate...";
219  this->SimulateHadronicFinalState(ev,p);
220  }
221  }
222 }
genie::exceptions::INukeException
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
Definition: INukeException.h:38
genie::utils::intranuke2018::PreEquilibrium
void PreEquilibrium(GHepRecord *ev, GHepParticle *p, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, const NuclearModelI *Nuclmodel, double NucRmvE, EINukeMode mode=kIMdHN)
Definition: INukeUtils2018.cxx:491
genie::GHepRecord::Particle
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
genie::HAIntranuke2018::SimulateHadronicFinalState
void SimulateHadronicFinalState(GHepRecord *ev, GHepParticle *p) const
Definition: HAIntranuke2018.cxx:125
genie::Intranuke2018::fFermiFac
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke2018.h:134
genie::HAIntranuke2018::fNumIterations
unsigned int fNumIterations
Definition: HAIntranuke2018.h:82
genie::Intranuke2018::fRemnA
int fRemnA
remnant nucleus A
Definition: Intranuke2018.h:119
genie::HAIntranuke2018::Inelastic
void Inelastic(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition: HAIntranuke2018.cxx:739
genie::HAIntranuke2018::InelasticHA
void InelasticHA(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition: HAIntranuke2018.cxx:553
genie::Intranuke2018::fNuclmodel
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke2018.h:118
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
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::Intranuke2018::fDoFermi
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke2018.h:136
pINFO
#define pINFO
Definition: Messenger.h:62
pWARN
#define pWARN
Definition: Messenger.h:60
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::HAIntranuke2018::SimulateHadronicFinalStateKinematics
void SimulateHadronicFinalStateKinematics(GHepRecord *ev, GHepParticle *p) const
Definition: HAIntranuke2018.cxx:168
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::Intranuke2018::fNucRmvE
double fNucRmvE
binding energy to subtract from cascade nucleons
Definition: Intranuke2018.h:127
genie::INukeFateHA_t
enum genie::EINukeFateHA_t INukeFateHA_t
genie::Intranuke2018::fRemnZ
int fRemnZ
remnant nucleus Z
Definition: Intranuke2018.h:120
genie::Intranuke2018::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke2018.h:121

Friends And Related Function Documentation

friend class IntranukeTester
friend

Definition at line 53 of file HAIntranuke2018.h.

Member Data Documentation

unsigned int genie::HAIntranuke2018::fNumIterations
mutableprivate

Definition at line 82 of file HAIntranuke2018.h.

Referenced by SimulateHadronicFinalState(), and SimulateHadronicFinalStateKinematics().

int genie::HAIntranuke2018::nuclA
mutableprivate

value of A for the target nucleus in hA mode

Definition at line 81 of file HAIntranuke2018.h.

Referenced by HadronFateHA(), and ProcessEventRecord().

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