GENIEGenerator
Loading...
Searching...
No Matches
genie::HNIntranuke2025 Class Reference

#include <HNIntranuke2025.h>

Inheritance diagram for genie::HNIntranuke2025:
[legend]
Collaboration diagram for genie::HNIntranuke2025:
[legend]

Public Member Functions

 HNIntranuke2025 ()
 HNIntranuke2025 (string config)
 ~HNIntranuke2025 ()
void ProcessEventRecord (GHepRecord *event_rec) const
virtual string GetINukeMode () const
virtual string GetGenINukeMode () const
Public Member Functions inherited from genie::Intranuke2025
 Intranuke2025 ()
 Intranuke2025 (string name)
 Intranuke2025 (string name, string config)
 ~Intranuke2025 ()
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 RegistryGetConfig (void) const
RegistryGetOwnedConfig (void)
virtual const AlgIdId (void) const
 Get algorithm ID.
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status.
virtual bool AllowReconfig (void) const
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm.
virtual void SetId (const AlgId &id)
 Set algorithm ID.
virtual void SetId (string name, string config)
const AlgorithmSubAlg (const RgKey &registry_key) const
void AdoptConfig (void)
void AdoptSubstructure (void)
virtual void Print (ostream &stream) const
 Print algorithm info.

Private Member Functions

void LoadConfig (void)
void SimulateHadronicFinalState (GHepRecord *ev, GHepParticle *p) const
INukeFateHN_t HadronFateHN (const GHepParticle *p) const
INukeFateHN_t HadronFateOset () const
double FateWeight (int pdgc, INukeFateHN_t fate) const
void ElasHN (GHepRecord *ev, GHepParticle *p, INukeFateHN_t fate) const
void AbsorbHN (GHepRecord *ev, GHepParticle *p, INukeFateHN_t fate) const
void InelasticHN (GHepRecord *ev, GHepParticle *p) const
void GammaInelasticHN (GHepRecord *ev, GHepParticle *p, INukeFateHN_t fate) const
bool HandleCompoundNucleusHN (GHepRecord *ev, GHepParticle *p) const
int HandleCompoundNucleus (GHepRecord *ev, GHepParticle *p, int mom) const

Private Attributes

int nuclA
 value of A for the target nucleus in hA mode
double fNucQEFac

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::Intranuke2025
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)
RegistryExtractLocalConfig (const Registry &in) const
RegistryExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key.
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.
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.
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
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership
int MergeTopRegistry (const Registry &r)
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships.
Protected Attributes inherited from genie::Intranuke2025
double fTrackingRadius
 tracking radius for the nucleus in the current event
TGenPhaseSpace fGenPhaseSpace
 a phase space generator
INukeHadroData2025fHadroData2025
 a collection of h+N,h+A data & calculations
AlgFactoryfAlgf
 algorithm factory instance
const NuclearModelIfNuclmodel
 nuclear model used to generate fermi momentum
int fRemnA
 remnant nucleus A
int fRemnZ
 remnant nucleus Z
TLorentzVector fRemnP4
 P4 of remnant system.
GEvGenMode_t fGMode
 event generation mode (lepton+A, hadron+A, ...)
double fR0
 effective nuclear size param
double fNR
 param multiplying the nuclear radius, determining how far to track hadrons beyond the "nuclear boundary"
double fNucRmvE
 binding energy to subtract from cascade nucleons
double fDelRPion
 factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement
double fDelRNucleon
 factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement
double fHadStep
 step size for intranuclear hadron transport
double fNucAbsFac
 absorption xsec correction factor (hN Mode)
double fNucCEXFac
 charge exchange xsec correction factor (hN Mode)
double fEPreEq
 threshold for pre-equilibrium reaction
double fFermiFac
 testing parameter to modify fermi momentum
double fFermiMomentum
 whether or not particle collision is pauli blocked
bool fDoFermi
 whether or not to do fermi mom.
bool fDoMassDiff
 whether or not to do mass diff. mode
bool fDoCompoundNucleus
 whether or not to do compound nucleus considerations
bool fUseOset
 Oset model for low energy pion in hN.
bool fAltOset
 NuWro's table-based implementation (not recommended)
bool fXsecNNCorr
 use nuclear medium correction for NN cross section
double fChPionMFPScale
 tweaking factors for tuning
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,...)
AlgId fID
 algorithm name and configuration set
vector< Registry * > fConfVect
vector< bool > fOwnerships
 ownership for every registry in fConfVect
AlgStatus_t fStatus
 algorithm execution status
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool)

Detailed Description

Definition at line 51 of file HNIntranuke2025.h.

Constructor & Destructor Documentation

◆ HNIntranuke2025() [1/2]

HNIntranuke2025::HNIntranuke2025 ( )

Definition at line 97 of file HNIntranuke2025.cxx.

97 :
98Intranuke2025("genie::HNIntranuke2025")
99{
100
101}

References genie::Intranuke2025::Intranuke2025().

◆ HNIntranuke2025() [2/2]

HNIntranuke2025::HNIntranuke2025 ( string config)

Definition at line 103 of file HNIntranuke2025.cxx.

103 :
104Intranuke2025("genie::HNIntranuke2025",config)
105{
106
107}

References genie::Intranuke2025::Intranuke2025().

◆ ~HNIntranuke2025()

HNIntranuke2025::~HNIntranuke2025 ( )

Definition at line 109 of file HNIntranuke2025.cxx.

110{
111
112}

Member Function Documentation

◆ AbsorbHN()

void HNIntranuke2025::AbsorbHN ( GHepRecord * ev,
GHepParticle * p,
INukeFateHN_t fate ) const
private

Definition at line 385 of file HNIntranuke2025.cxx.

387{
388 // handles pi+d->2p, pi-d->nn, pi0 d->pn absorbtion, all using pi+d values
389
390 int pdgc = p->Pdg();
391
392#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
393 LOG("HNIntranuke2025", pDEBUG)
394 << "AbsorbHN() is invoked for a : " << p->Name()
395 << " whose fate is : " << INukeHadroFates2025::AsString(fate);
396#endif
397
398 // check fate
399 if(fate!=kIHNFtAbs)
400 {
401 LOG("HNIntranuke2025", pWARN)
402 << "AbsorbHN() cannot handle fate: " << INukeHadroFates2025::AsString(fate);
403 return;
404 }
405
406 // random number generator
407 RandomGen * rnd = RandomGen::Instance();
408
409 // Notes on the kinematics
410 // -- Simple variables are used for efficiency
411 // -- Variables are numbered according to particle
412 // -- -- #1 -> incoming particle
413 // -- -- #2 -> target (here, 2_1 and 2_2 for individual particles)
414 // -- -- #3 -> scattered incoming (Particle tracked in hA mode)
415 // -- -- #4 -> other scattered particle
416 // -- Suffix "L" is for lab frame, suffix "CM" is for center of mass frame
417 // -- Subscript "z" is for parallel component, "t" is for transverse
418
419 int pcode, t1code, t2code, scode, s2code; // particles
420 double M1, M2_1, M2_2, M3, M4; // rest energies, in GeV
421 double E1L, P1L, E2L, P2L, E3L, P3L, E4L, P4L;
422 double P1zL, P2zL;
423 double beta, gm; // speed and gamma for CM frame in lab
424 double Et, E2CM;
425 double C3CM, S3CM; // cos and sin of scattering angle
426 double Theta1, Theta2, theta5;
427 double PHI3; // transverse scattering angle
428 double E1CM, E3CM, E4CM, P3CM;
429 double P3zL, P3tL, P4zL, P4tL;
430 double E2_1L, E2_2L;
431 TVector3 tP2_1L, tP2_2L, tP1L, tP2L, tPtot, P1zCM, P2zCM;
432 TVector3 tP3L, tP4L;
433 TVector3 bDir, tTrans, tbeta, tVect;
434
435 // Library instance for reference
436 PDGLibrary * pLib = PDGLibrary::Instance();
437
438 // Handle fermi target
439 Target target(ev->TargetNucleus()->Pdg());
440
441 // Target should be a deuteron, but for now
442 // handling it as seperate nucleons
443 if(pdgc==211) // pi-plus
444 {
445 pcode = 211;
446 t1code = 2212; // proton
447 t2code = 2112; // neutron
448 scode = 2212;
449 s2code = 2212;
450 }
451 else if(pdgc==-211) // pi-minus
452 {
453 pcode = -211;
454 t1code = 2212;
455 t2code = 2112;
456 scode = 2112;
457 s2code = 2112;
458 }
459 else if(pdgc==111) // pi-zero
460 {
461 pcode = 111;
462 t1code = 2212;
463 t2code = 2112;
464 scode = 2212;
465 s2code = 2112;
466 }
467 else
468 {
469 LOG("HNIntranuke2025", pWARN)
470 << "AbsorbHN() cannot handle probe: " << pdgc;
471 return;
472 }
473
474 // assign proper masses
475 M1 = pLib->Find(pcode) ->Mass();
476 M2_1 = pLib->Find(t1code)->Mass();
477 M2_2 = pLib->Find(t2code)->Mass();
478 M3 = pLib->Find(scode) ->Mass();
479 M4 = pLib->Find(s2code)->Mass();
480
481 // handle fermi momentum
482 if(fDoFermi)
483 {
484 target.SetHitNucPdg(t1code);
485 fNuclmodel->GenerateNucleon(target);
486 tP2_1L=fFermiFac * fNuclmodel->Momentum3();
487 E2_1L = TMath::Sqrt(tP2_1L.Mag2() + M2_1*M2_1);
488
489 target.SetHitNucPdg(t2code);
490 fNuclmodel->GenerateNucleon(target);
491 tP2_2L=fFermiFac * fNuclmodel->Momentum3();
492 E2_2L = TMath::Sqrt(tP2_2L.Mag2() + M2_2*M2_2);
493 }
494 else
495 {
496 tP2_1L.SetXYZ(0.0, 0.0, 0.0);
497 E2_1L = M2_1;
498
499 tP2_2L.SetXYZ(0.0, 0.0, 0.0);
500 E2_2L = M2_2;
501 }
502
503 E2L = E2_1L + E2_2L;
504
505 // adjust p to reflect scattering
506 // get random scattering angle
507 C3CM = fHadroData2025->IntBounce(p,t1code,scode,fate);
508 if (C3CM<-1.)
509 {
511 ev->AddParticle(*p);
512 return;
513 }
514 S3CM = TMath::Sqrt(1.0 - C3CM*C3CM);
515
516 // Get lab energy and momenta
517 E1L = p->E();
518 if(E1L<0.001) E1L=0.001;
519 P1L = TMath::Sqrt(E1L*E1L - M1*M1);
520 tP1L = p->P4()->Vect();
521 tP2L = tP2_1L + tP2_2L;
522 P2L = tP2L.Mag();
523 tPtot = tP1L + tP2L;
524
525 // get unit vectors and angles needed for later
526 bDir = tPtot.Unit();
527 Theta1 = tP1L.Angle(bDir);
528 Theta2 = tP2L.Angle(bDir);
529
530 // get parallel and transverse components
531 P1zL = P1L*TMath::Cos(Theta1);
532 P2zL = P2L*TMath::Cos(Theta2);
533 tVect.SetXYZ(1,0,0);
534 if(TMath::Abs((tVect - bDir).Mag())<.01) tVect.SetXYZ(0,1,0);
535 theta5 = tVect.Angle(bDir);
536 tTrans = (tVect - TMath::Cos(theta5)*bDir).Unit();
537
538 // calculate beta and gamma
539 tbeta = tPtot * (1.0 / (E1L + E2L));
540 beta = tbeta.Mag();
541 gm = 1.0 / TMath::Sqrt(1.0 - beta*beta);
542
543 // boost to CM frame to get scattered particle momenta
544 E1CM = gm*E1L - gm*beta*P1zL;
545 P1zCM = gm*P1zL*bDir - gm*tbeta*E1L;
546 E2CM = gm*E2L - gm*beta*P2zL;
547 P2zCM = gm*P2zL*bDir - gm*tbeta*E2L;
548 Et = E1CM + E2CM;
549 E3CM = (Et*Et + (M3*M3) - (M4*M4)) / (2.0*Et);
550 E4CM = Et - E3CM;
551 P3CM = TMath::Sqrt(E3CM*E3CM - M3*M3);
552
553 // boost back to lab
554 P3zL = gm*beta*E3CM + gm*P3CM*C3CM;
555 P3tL = P3CM*S3CM;
556 P4zL = gm*beta*E4CM + gm*P3CM*(-C3CM);
557 P4tL = P3CM*(-S3CM);
558
559 P3L = TMath::Sqrt(P3zL*P3zL + P3tL*P3tL);
560 P4L = TMath::Sqrt(P4zL*P4zL + P4tL*P4tL);
561
562 // check for too small values
563 // may introduce error, so warn if it occurs
564 if(!(TMath::Finite(P3L))||P3L<.001)
565 {
566 LOG("HNIntranuke2025",pINFO)
567 << "Particle 3 " << M3 << " momentum small or non-finite: " << P3L
568 << "\n" << "--> Assigning .001 as new momentum";
569 P3tL = 0;
570 P3zL = .001;
571 P3L = .001;
572 E3L = TMath::Sqrt(P3L*P3L + M3*M3);
573 }
574
575 if(!(TMath::Finite(P4L))||P4L<.001)
576 {
577 LOG("HNIntranuke2025",pINFO)
578 << "Particle 4 " << M4 << " momentum small or non-finite: " << P4L
579 << "\n" << "--> Assigning .001 as new momentum";
580 P4tL = 0;
581 P4zL = .001;
582 P4L = .001;
583 E4L = TMath::Sqrt(P4L*P4L + M4*M4);
584 }
585
586 // pauli blocking (do not apply PB for Oset)
587 //if(!fUseOset && (P3L < fFermiMomentum || P4L < fFermiMomentum))
588 double ke = p->KinE() / units::MeV;
589 if((!fUseOset || ke > 350.0 ) && (P3L < fFermiMomentum || P4L < fFermiMomentum))
590 {
591#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
592 LOG("HNIntranuke2025",pINFO) << "AbsorbHN failed: Pauli blocking";
593#endif
594 /*
595 p->SetStatus(kIStHadronInTheNucleus);
596 //disable until needed
597 // utils::intranuke2025::StepParticle(p,fFreeStep,fTrackingRadius);
598 ev->AddParticle(*p);
599 return;
600 */
601 // new attempt at error handling:
602 LOG("HNIntranuke2025", pINFO) << "AbsorbHN failed: Pauli blocking";
603 exceptions::INukeException exception;
604 exception.SetReason("hN absorption failed");
605 throw exception;
606 }
607
608 // handle remnant nucleus updates
609 fRemnZ--;
610 fRemnA -=2;
611 fRemnP4 -= TLorentzVector(tP2_1L,E2_1L);
612 fRemnP4 -= TLorentzVector(tP2_2L,E2_2L);
613
614 // get random phi angle, distributed uniformally in 360 deg
615 PHI3 = 2 * kPi * rnd->RndFsi().Rndm();
616
617 tP3L = P3zL*bDir + P3tL*tTrans;
618 tP4L = P4zL*bDir + P4tL*tTrans;
619
620 tP3L.Rotate(PHI3,bDir); // randomize transverse components
621 tP4L.Rotate(PHI3,bDir);
622
623 E3L = TMath::Sqrt(P3L*P3L + M3*M3);
624 E4L = TMath::Sqrt(P4L*P4L + M4*M4);
625
626 // create t particle w/ appropriate momenta, code, and status
627 // set target's mom to be the mom of the hadron that was cloned
628 GHepParticle * t = new GHepParticle(*p);
630 t->SetLastMother(p->LastMother());
631
632 TLorentzVector t4P4L(tP4L,E4L);
633 t->SetPdgCode(s2code);
634 t->SetMomentum(t4P4L);
636
637 // adjust p to reflect scattering
638 TLorentzVector t4P3L(tP3L,E3L);
639 p->SetPdgCode(scode);
640 p->SetMomentum(t4P3L);
642
643#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
644 LOG("HNIntranuke2025", pDEBUG)
645 << "|p3| = " << (P3L) << ", E3 = " << (E3L);
646 LOG("HNIntranuke2025", pDEBUG)
647 << "|p4| = " << (P4L) << ", E4 = " << (E4L);
648#endif
649
650 ev->AddParticle(*p);
651 ev->AddParticle(*t);
652
653 delete t; // delete particle clone
654}
#define pINFO
Definition Messenger.h:62
#define pDEBUG
Definition Messenger.h:63
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
#define pWARN
Definition Messenger.h:60
string Name(void) const
Name that corresponds to the PDG code.
int FirstMother(void) const
void SetMomentum(const TLorentzVector &p4)
int Pdg(void) const
void SetFirstMother(int m)
void SetLastMother(int m)
const TLorentzVector * P4(void) const
int LastMother(void) const
void SetStatus(GHepStatus_t s)
double E(void) const
Get energy.
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
virtual GHepParticle * TargetNucleus(void) const
virtual void AddParticle(const GHepParticle &p)
static string AsString(INukeFateHN_t fate)
int fRemnA
remnant nucleus A
int fRemnZ
remnant nucleus Z
double fFermiMomentum
whether or not particle collision is pauli blocked
bool fUseOset
Oset model for low energy pion in hN.
TLorentzVector fRemnP4
P4 of remnant system.
double fFermiFac
testing parameter to modify fermi momentum
INukeHadroData2025 * fHadroData2025
a collection of h+N,h+A data & calculations
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
bool fDoFermi
whether or not to do fermi mom.
static PDGLibrary * Instance(void)
TParticlePDG * Find(int pdgc, bool must_exist=true)
static RandomGen * Instance()
Access instance.
Definition RandomGen.cxx:74
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition RandomGen.h:59
static constexpr double MeV
Definition Units.h:129
@ kIStHadronInTheNucleus
Definition GHepStatus.h:37
@ kIStStableFinalState
Definition GHepStatus.h:30

References genie::GHepRecord::AddParticle(), genie::INukeHadroFates2025::AsString(), genie::GHepParticle::E(), genie::Intranuke2025::fDoFermi, genie::Intranuke2025::fFermiFac, genie::Intranuke2025::fFermiMomentum, genie::Intranuke2025::fHadroData2025, genie::PDGLibrary::Find(), genie::GHepParticle::FirstMother(), genie::Intranuke2025::fNuclmodel, genie::Intranuke2025::fRemnA, genie::Intranuke2025::fRemnP4, genie::Intranuke2025::fRemnZ, genie::Intranuke2025::fUseOset, genie::PDGLibrary::Instance(), genie::RandomGen::Instance(), genie::kIHNFtAbs, genie::GHepParticle::KinE(), genie::kIStHadronInTheNucleus, genie::kIStStableFinalState, genie::constants::kPi, genie::GHepParticle::LastMother(), LOG, genie::units::MeV, genie::GHepParticle::Name(), genie::GHepParticle::P4(), pDEBUG, genie::GHepParticle::Pdg(), pINFO, 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(), and genie::GHepRecord::TargetNucleus().

Referenced by SimulateHadronicFinalState().

◆ ElasHN()

void HNIntranuke2025::ElasHN ( GHepRecord * ev,
GHepParticle * p,
INukeFateHN_t fate ) const
private

Definition at line 656 of file HNIntranuke2025.cxx.

658{
659 // scatters particle within nucleus
660
661#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
662 LOG("HNIntranuke2025", pDEBUG)
663 << "ElasHN() is invoked for a : " << p->Name()
664 << " whose fate is : " << INukeHadroFates2025::AsString(fate);
665#endif
666
667 if(fate!=kIHNFtCEx && fate!=kIHNFtElas)
668 {
669 LOG("HNIntranuke2025", pWARN)
670 << "ElasHN() cannot handle fate: " << INukeHadroFates2025::AsString(fate);
671 return;
672 }
673
674 // Random number generator
675 RandomGen * rnd = RandomGen::Instance();
676
677 // vars for incoming particle, target, and scattered pdg codes
678 int pcode = p->Pdg();
679 int tcode, scode, s2code;
680 double ppcnt = (double) fRemnZ / (double) fRemnA; // % of protons
681
682 // Select a target randomly, weighted to #
683 // -- Unless, of course, the fate is CEx,
684 // -- in which case the target may be deterministic
685 // Also assign scattered particle code
686 if(fate==kIHNFtCEx)
687 {
688 if(pcode==kPdgPiP) {tcode = kPdgNeutron; scode = kPdgPi0; s2code = kPdgProton;}
689 else if(pcode==kPdgPiM) {tcode = kPdgProton; scode = kPdgPi0; s2code = kPdgNeutron;}
690 else if(pcode==kPdgKP) {tcode = kPdgNeutron; scode = kPdgK0; s2code = kPdgProton;}
691 else
692 {
693 // for pi0
694 tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton) :(kPdgNeutron);
695 scode = (tcode == kPdgProton) ?(kPdgPiP) :(kPdgPiM);
696 s2code = (tcode == kPdgProton) ?(kPdgNeutron):(kPdgProton);
697 }
698 }
699 else
700 {
701 tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton):(kPdgNeutron);
702 scode = pcode;
703 s2code = tcode;
704 }
705
706 // get random scattering angle
707 double C3CM = fHadroData2025->IntBounce(p,tcode,scode,fate);
708 if (C3CM<-1.)
709 {
711 ev->AddParticle(*p);
712 return;
713 }
714
715 // create scattered particle
716 GHepParticle * t = new GHepParticle(*p);
717 t->SetPdgCode(tcode);
718 double Mt = t->Mass();
719 //t->SetMomentum(TLorentzVector(0,0,0,Mt));
720 t->SetRemovalEnergy(0);
721 // handle fermi momentum
722 if(fDoFermi)
723 {
724 // Handle fermi target
725 Target target(ev->TargetNucleus()->Pdg());
726 //LOG("HAIntranuke2025", pNOTICE) << "Nuclmodel= " << fNuclmodel->ModelType(target) ;
727 target.SetHitNucPdg(tcode);
728 fNuclmodel->GenerateNucleon(target);
729 TVector3 tP3L = fFermiFac * fNuclmodel->Momentum3();
730 double tE = TMath::Sqrt(tP3L.Mag2() + Mt*Mt);
731 t->SetMomentum(TLorentzVector(tP3L,tE));
732 }
733 else
734 {
735 t->SetMomentum(TLorentzVector(0,0,0,Mt));
736 }
737
738 bool pass = utils::intranuke2025::TwoBodyCollision(ev,pcode,tcode,scode,s2code,C3CM,
740
741#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
742 LOG("HNIntranuke2025",pDEBUG)
743 << "|p3| = " << P3L << ", E3 = " << E3L;
744 LOG("HNIntranuke2025",pDEBUG)
745 << "|p4| = " << P4L << ", E4 = " << E4L;
746#endif
747
748 if (pass==true)
749 {
750 ev->AddParticle(*p);
751 ev->AddParticle(*t);
752 } else
753 {
754 delete t; //fixes memory leak
755 LOG("HNIntranuke2025", pINFO) << "Elastic in hN failed calling TwoBodyCollision";
756 exceptions::INukeException exception;
757 exception.SetReason("hN scattering kinematics through TwoBodyCollision failed");
758 throw exception;
759 }
760
761 delete t;
762
763}
void SetRemovalEnergy(double Erm)
double Mass(void) const
Mass that corresponds to the PDG code.
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.
const int kPdgPiM
Definition PDGCodes.h:159
const int kPdgProton
Definition PDGCodes.h:81
const int kPdgPi0
Definition PDGCodes.h:160
@ kIMdHN
Definition INukeMode.h:32
const int kPdgKP
Definition PDGCodes.h:172
const int kPdgNeutron
Definition PDGCodes.h:83
const int kPdgPiP
Definition PDGCodes.h:158
const int kPdgK0
Definition PDGCodes.h:174

References genie::GHepRecord::AddParticle(), genie::INukeHadroFates2025::AsString(), genie::Intranuke2025::fDoFermi, genie::Intranuke2025::fFermiFac, genie::Intranuke2025::fHadroData2025, genie::Intranuke2025::fNuclmodel, genie::Intranuke2025::fRemnA, genie::Intranuke2025::fRemnP4, genie::Intranuke2025::fRemnZ, genie::RandomGen::Instance(), genie::kIHNFtCEx, genie::kIHNFtElas, genie::kIMdHN, genie::kIStStableFinalState, genie::kPdgK0, genie::kPdgKP, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, LOG, genie::GHepParticle::Mass(), genie::GHepParticle::Name(), pDEBUG, genie::GHepParticle::Pdg(), pINFO, pWARN, genie::RandomGen::RndFsi(), genie::Target::SetHitNucPdg(), genie::GHepParticle::SetMomentum(), genie::GHepParticle::SetPdgCode(), genie::exceptions::INukeException::SetReason(), genie::GHepParticle::SetRemovalEnergy(), genie::GHepParticle::SetStatus(), genie::GHepRecord::TargetNucleus(), and genie::utils::intranuke2025::TwoBodyCollision().

Referenced by SimulateHadronicFinalState().

◆ FateWeight()

double HNIntranuke2025::FateWeight ( int pdgc,
INukeFateHN_t fate ) const
private

Definition at line 360 of file HNIntranuke2025.cxx.

361{
362 // turn fates off if the remnant nucleus does not have the number of p,n
363 // required
364
365 int np = fRemnZ;
366 int nn = fRemnA - fRemnZ;
367
368 if (np < 1 && nn < 1)
369 {
370 LOG("HNIntranuke2025", pERROR) << "** Nothing left in nucleus!!! **";
371 return 0;
372 }
373 else
374 {
375 if (fate == kIHNFtCEx && pdgc==kPdgPiP ) { return (nn>=1) ? 1. : 0.; }
376 if (fate == kIHNFtCEx && pdgc==kPdgPiM ) { return (np>=1) ? 1. : 0.; }
377 if (fate == kIHNFtCEx && pdgc==kPdgKP ) { return (nn>=1) ? 1. : 0.; } //Added, changed np to nn
378 if (fate == kIHNFtAbs) { return ((nn>=1) && (np>=1)) ? 1. : 0.; }
379 if (fate == kIHNFtCmp ) { return ((pdgc==kPdgProton||pdgc==kPdgNeutron)&&fDoCompoundNucleus&&fRemnA>5) ? 1. : 0.; }
380
381 }
382 return 1.;
383}
#define pERROR
Definition Messenger.h:59
bool fDoCompoundNucleus
whether or not to do compound nucleus considerations

References genie::Intranuke2025::fDoCompoundNucleus, genie::Intranuke2025::fRemnA, genie::Intranuke2025::fRemnZ, genie::kIHNFtAbs, genie::kIHNFtCEx, genie::kIHNFtCmp, genie::kPdgKP, genie::kPdgNeutron, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, LOG, and pERROR.

Referenced by HadronFateHN().

◆ GammaInelasticHN()

void HNIntranuke2025::GammaInelasticHN ( GHepRecord * ev,
GHepParticle * p,
INukeFateHN_t fate ) const
private

Definition at line 801 of file HNIntranuke2025.cxx.

802{
803 // This function handles pion photoproduction reactions
804
805#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
806 LOG("HNIntranuke2025", pDEBUG)
807 << "GammaInelasticHN() is invoked for a : " << p->Name()
808 << " whose fate is : " << INukeHadroFates2025::AsString(fate);
809#endif
810
811 if(fate!=kIHNFtInelas && p->Pdg()!=kPdgGamma)
812 {
813 LOG("HNIntranuke2025", pWARN)
814 << "GammaInelasticHN() cannot handle fate: " << INukeHadroFates2025::AsString(fate);
815 return;
816 }
817
818 // random number generator
819 RandomGen * rnd = RandomGen::Instance();
820
821 // vars for incoming particle, target, and scattered reaction products
822 double ppcnt = (double) fRemnZ / (double) fRemnA; // % of protons
823 int pcode = p->Pdg();
824 int tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton):(kPdgNeutron);
825 int scode, s2code;
826 double ke = p->KinE() / units::MeV;
827
828 LOG("HNIntranuke2025", pNOTICE)
829 << "Particle code: " << pcode << ", target: " << tcode;
830
831
832 if (rnd->RndFsi().Rndm() * (fHadroData2025 -> XSecGamp_fs() -> Evaluate(ke) +
833 fHadroData2025 -> XSecGamn_fs() -> Evaluate(ke) )
834 <= fHadroData2025 -> XSecGamp_fs() -> Evaluate(ke) ) { scode = kPdgProton; }
835 else { scode = kPdgNeutron; }
836
837 //scode=fHadroData2025->AngleAndProduct(p,tcode,C3CM,fate);
838 //double C3CM = 0.0; // cos of scattering angle
839 double C3CM = fHadroData2025->IntBounce(p,tcode,scode,fate);
840
841 if ((tcode == kPdgProton ) && (scode==kPdgProton )) {s2code=kPdgPi0;}
842 else if ((tcode == kPdgProton ) && (scode==kPdgNeutron)) {s2code=kPdgPiP;}
843 else if ((tcode == kPdgNeutron) && (scode==kPdgProton )) {s2code=kPdgPiM;}
844 else if ((tcode == kPdgNeutron) && (scode==kPdgNeutron)) {s2code=kPdgPi0;}
845 else {
846 LOG("HNIntranuke2025", pERROR)
847 << "Error: could not determine particle final states";
848 ev->AddParticle(*p);
849 return;
850 }
851
852 LOG("HNIntranuke2025", pNOTICE)
853 << "GammaInelastic fate: " << INukeHadroFates2025::AsString(fate);
854 LOG("HNIntranuke2025", pNOTICE)
855 << " final state: " << scode << " and " << s2code;
856 LOG("HNIntranuke2025", pNOTICE)
857 << " scattering angle: " << C3CM;
858
859 GHepParticle * t = new GHepParticle(*p);
860 t->SetPdgCode(tcode);
861 double Mt = t->Mass();
862
863 // handle fermi momentum
864 if(fDoFermi)
865 {
866 // Handle fermi target
867 Target target(ev->TargetNucleus()->Pdg());
868
869 target.SetHitNucPdg(tcode);
870 fNuclmodel->GenerateNucleon(target);
871 TVector3 tP3L = fFermiFac * fNuclmodel->Momentum3();
872 double tE = TMath::Sqrt(tP3L.Mag2() + Mt*Mt);
873 t->SetMomentum(TLorentzVector(tP3L,tE));
874 }
875 else
876 {
877 t->SetMomentum(TLorentzVector(0,0,0,Mt));
878 }
879
880 bool pass = utils::intranuke2025::TwoBodyCollision(ev,pcode,tcode,scode,s2code,C3CM,
882
883#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
884 LOG("HNIntranuke2025",pDEBUG)
885 << "|p3| = " << P3L << ", E3 = " << E3L;
886 LOG("HNIntranuke2025",pDEBUG)
887 << "|p4| = " << P4L << ", E4 = " << E4L;
888#endif
889
890 if (pass==true)
891 {
892 //p->SetStatus(kIStStableFinalState);
893 //t->SetStatus(kIStStableFinalState);
894 ev->AddParticle(*p);
895 ev->AddParticle(*t);
896 } else
897 {
898 ev->AddParticle(*p);
899 }
900
901 delete t;
902
903}
#define pNOTICE
Definition Messenger.h:61
const int kPdgGamma
Definition PDGCodes.h:189

References genie::GHepRecord::AddParticle(), genie::INukeHadroFates2025::AsString(), genie::Intranuke2025::fDoFermi, genie::Intranuke2025::fFermiFac, genie::Intranuke2025::fHadroData2025, genie::Intranuke2025::fNuclmodel, genie::Intranuke2025::fRemnA, genie::Intranuke2025::fRemnP4, genie::Intranuke2025::fRemnZ, genie::RandomGen::Instance(), genie::kIHNFtInelas, genie::kIMdHN, genie::GHepParticle::KinE(), genie::kPdgGamma, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, LOG, genie::GHepParticle::Mass(), genie::units::MeV, genie::GHepParticle::Name(), pDEBUG, genie::GHepParticle::Pdg(), pERROR, pNOTICE, pWARN, genie::RandomGen::RndFsi(), genie::Target::SetHitNucPdg(), genie::GHepParticle::SetMomentum(), genie::GHepParticle::SetPdgCode(), genie::GHepRecord::TargetNucleus(), and genie::utils::intranuke2025::TwoBodyCollision().

Referenced by SimulateHadronicFinalState().

◆ GetGenINukeMode()

virtual string genie::HNIntranuke2025::GetGenINukeMode ( ) const
inlinevirtual

Reimplemented from genie::Intranuke2025.

Definition at line 63 of file HNIntranuke2025.h.

63{return "hN";};

◆ GetINukeMode()

virtual string genie::HNIntranuke2025::GetINukeMode ( ) const
inlinevirtual

Reimplemented from genie::Intranuke2025.

Definition at line 62 of file HNIntranuke2025.h.

62{return "hN2025";};

◆ HadronFateHN()

INukeFateHN_t HNIntranuke2025::HadronFateHN ( const GHepParticle * p) const
private

Definition at line 205 of file HNIntranuke2025.cxx.

206{
207// Select a hadron fate in HN mode
208//
209 RandomGen * rnd = RandomGen::Instance();
210
211 // get pdgc code & kinetic energy in MeV
212 int pdgc = p->Pdg();
213 double ke = p->KinE() / units::MeV;
214
215 bool isPion = (pdgc == kPdgPiP or pdgc == kPdgPi0 or pdgc == kPdgPiM);
216
217 if (isPion and fUseOset and ke < 350.0) return HadronFateOset ();
218
219 LOG("HNIntranuke2025", pNOTICE)
220 << "Selecting hN fate for " << p->Name() << " with KE = " << ke << " MeV";
221
222 // try to generate a hadron fate
223 unsigned int iter = 0;
224 while(iter++ < kRjMaxIterations) {
225
226 // handle pions
227 //
228 if (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0) {
229
230 double frac_cex = this->FateWeight(pdgc, kIHNFtCEx)
231 * fHadroData2025->Frac(pdgc, kIHNFtCEx, ke, fRemnA, fRemnZ);
232 double frac_elas = this->FateWeight(pdgc, kIHNFtElas)
233 * fHadroData2025->Frac(pdgc, kIHNFtElas, ke, fRemnA, fRemnZ);
234 double frac_inel = this->FateWeight(pdgc, kIHNFtInelas)
235 * fHadroData2025->Frac(pdgc, kIHNFtInelas, ke, fRemnA, fRemnZ);
236 double frac_abs = this->FateWeight(pdgc, kIHNFtAbs)
237 * fHadroData2025->Frac(pdgc, kIHNFtAbs, ke, fRemnA, fRemnZ);
238
239 frac_cex *= fNucCEXFac; // scaling factors
240 frac_abs *= fNucAbsFac;
241 frac_elas *= fNucQEFac;
242 if(pdgc==kPdgPi0) frac_abs*= 0.665; //isospin factor
243
244 LOG("HNIntranuke2025", pNOTICE)
245 << "\n frac{" << INukeHadroFates2025::AsString(kIHNFtCEx) << "} = " << frac_cex
246 << "\n frac{" << INukeHadroFates2025::AsString(kIHNFtElas) << "} = " << frac_elas
247 << "\n frac{" << INukeHadroFates2025::AsString(kIHNFtInelas) << "} = " << frac_inel
248 << "\n frac{" << INukeHadroFates2025::AsString(kIHNFtAbs) << "} = " << frac_abs;
249
250 // compute total fraction (can be <1 if fates have been switched off)
251 double tf = frac_cex +
252 frac_elas +
253 frac_inel +
254 frac_abs;
255
256 double r = tf * rnd->RndFsi().Rndm();
257#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
258 LOG("HNIntranuke2025", pDEBUG) << "r = " << r << " (max = " << tf << ")";
259#endif
260
261 double cf=0; // current fraction
262
263 if(r < (cf += frac_cex )) return kIHNFtCEx; //cex
264 if(r < (cf += frac_elas )) return kIHNFtElas; //elas
265 if(r < (cf += frac_inel )) return kIHNFtInelas; //inelas
266 if(r < (cf += frac_abs )) return kIHNFtAbs; //abs
267
268 LOG("HNIntranuke2025", pWARN)
269 << "No selection after going through all fates! "
270 << "Total fraction = " << tf << " (r = " << r << ")";
271 ////////////////////////////
272 return kIHNFtUndefined;
273 }
274
275 // handle nucleons
276 else if (pdgc==kPdgProton || pdgc==kPdgNeutron) {
277
278 double frac_elas = this->FateWeight(pdgc, kIHNFtElas)
279 * fHadroData2025->Frac(pdgc, kIHNFtElas, ke, fRemnA, fRemnZ);
280 double frac_inel = this->FateWeight(pdgc, kIHNFtInelas)
281 * fHadroData2025->Frac(pdgc, kIHNFtInelas, ke, fRemnA, fRemnZ);
282 double frac_cmp = this->FateWeight(pdgc, kIHNFtCmp)
283 * fHadroData2025->Frac(pdgc, kIHNFtCmp, ke, fRemnA , fRemnZ);
284
285 LOG("HNIntranuke2025", pINFO)
286 << "\n frac{" << INukeHadroFates2025::AsString(kIHNFtElas) << "} = " << frac_elas
287 << "\n frac{" << INukeHadroFates2025::AsString(kIHNFtInelas) << "} = " << frac_inel;
288
289 // compute total fraction (can be <1 if fates have been switched off)
290 double tf = frac_elas +
291 frac_inel +
292 frac_cmp;
293
294 double r = tf * rnd->RndFsi().Rndm();
295
296#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
297 LOG("HNIntranuke2025", pDEBUG) << "r = " << r << " (max = " << tf << ")";
298#endif
299
300 double cf=0; // current fraction
301 if(r < (cf += frac_elas )) return kIHNFtElas; // elas
302 if(r < (cf += frac_inel )) return kIHNFtInelas; // inelas
303 if(r < (cf += frac_cmp )) return kIHNFtCmp; // cmp
304
305 LOG("HNIntranuke2025", pWARN)
306 << "No selection after going through all fates! "
307 << "Total fraction = " << tf << " (r = " << r << ")";
308 //////////////////////////
309 return kIHNFtUndefined;
310 }
311
312 // handle gamma -- does not currently consider the elastic case
313 else if (pdgc==kPdgGamma) return kIHNFtInelas;
314 // Handle kaon -- elastic + charge exchange
315 else if (pdgc==kPdgKP){
316 double frac_cex = this->FateWeight(pdgc, kIHNFtCEx)
317 * fHadroData2025->Frac(pdgc, kIHNFtCEx, ke, fRemnA, fRemnZ);
318 double frac_elas = this->FateWeight(pdgc, kIHNFtElas)
319 * fHadroData2025->Frac(pdgc, kIHNFtElas, ke, fRemnA, fRemnZ);
320
321 // frac_cex *= fNucCEXFac; // scaling factors
322 // frac_elas *= fNucQEFac; // Flor - Correct scaling factors?
323
324 LOG("HNIntranuke", pINFO)
325 << "\n frac{" << INukeHadroFates2025::AsString(kIHNFtCEx) << "} = " << frac_cex
326 << "\n frac{" << INukeHadroFates2025::AsString(kIHNFtElas) << "} = " << frac_elas;
327
328 // compute total fraction (can be <1 if fates have been switched off)
329 double tf = frac_cex +
330 frac_elas;
331
332 double r = tf * rnd->RndFsi().Rndm();
333#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
334 LOG("HNIntranuke", pDEBUG) << "r = " << r << " (max = " << tf << ")";
335#endif
336
337 double cf=0; // current fraction
338
339 if(r < (cf += frac_cex )) return kIHNFtCEx; //cex
340 if(r < (cf += frac_elas )) return kIHNFtElas; //elas
341
342 LOG("HNIntranuke", pWARN)
343 << "No selection after going through all fates! "
344 << "Total fraction = " << tf << " (r = " << r << ")";
345 ////////////////////////////
346 return kIHNFtUndefined;
347 }//End K+
348
349 /*else if (pdgc==kPdgKM){
350
351 return kIHNFtElas;
352 }//End K-*/
353
354 }//iterations
355
356 return kIHNFtUndefined;
357}
INukeFateHN_t HadronFateOset() const
double FateWeight(int pdgc, INukeFateHN_t fate) const
double fNucAbsFac
absorption xsec correction factor (hN Mode)
double fNucCEXFac
charge exchange xsec correction factor (hN Mode)
static const unsigned int kRjMaxIterations
Definition Controls.h:26
@ kIHNFtUndefined

References genie::INukeHadroFates2025::AsString(), FateWeight(), genie::Intranuke2025::fHadroData2025, genie::Intranuke2025::fNucAbsFac, genie::Intranuke2025::fNucCEXFac, fNucQEFac, genie::Intranuke2025::fRemnA, genie::Intranuke2025::fRemnZ, genie::Intranuke2025::fUseOset, HadronFateOset(), genie::RandomGen::Instance(), genie::kIHNFtAbs, genie::kIHNFtCEx, genie::kIHNFtCmp, genie::kIHNFtElas, genie::kIHNFtInelas, genie::kIHNFtUndefined, genie::GHepParticle::KinE(), genie::kPdgGamma, genie::kPdgKP, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, genie::controls::kRjMaxIterations, LOG, genie::units::MeV, genie::GHepParticle::Name(), pDEBUG, genie::GHepParticle::Pdg(), pINFO, pNOTICE, pWARN, and genie::RandomGen::RndFsi().

Referenced by SimulateHadronicFinalState().

◆ HadronFateOset()

INukeFateHN_t HNIntranuke2025::HadronFateOset ( ) const
private

Definition at line 1009 of file HNIntranuke2025.cxx.

1010{
1011 //LOG("HNIntranuke2025", pWARN) << "IN HadronFateOset";
1012
1013 //LOG("HNIntranuke2025", pWARN) << "{ frac abs = " << osetUtils::currentInstance->getAbsorptionFraction();
1014 //LOG("HNIntranuke2025", pWARN) << " frac cex = " << osetUtils::currentInstance->getCexFraction() << " }";
1015
1016 double fractionAbsorption = osetUtils::currentInstance->getAbsorptionFraction();
1017 double fractionCex = osetUtils::currentInstance->getCexFraction ();
1018 double fractionElas = 1 - (fractionAbsorption + fractionCex);
1019
1020 fractionCex *= fNucCEXFac; // scaling factors
1021 fractionAbsorption *= fNucAbsFac;
1022 fractionElas *= fNucQEFac;
1023
1024 double totalFrac = fractionCex + fractionAbsorption + fractionElas;
1025
1026 RandomGen *randomGenerator = RandomGen::Instance();
1027 const double randomNumber = randomGenerator->RndFsi().Rndm() * totalFrac;
1028
1029 LOG("HNIntranuke2025", pNOTICE) << "{ frac abs = " << fractionAbsorption;
1030 LOG("HNIntranuke2025", pNOTICE) << " frac cex = " << fractionCex;
1031 LOG("HNIntranuke2025", pNOTICE) << " frac elas = " << fractionElas << " }";
1032
1033 if (randomNumber < fractionAbsorption && fRemnA > 1) return kIHNFtAbs;
1034 else if (randomNumber < fractionAbsorption + fractionCex) return kIHNFtCEx;
1035 else return kIHNFtElas;
1036}
double getCexFraction() const
return fraction of cex events
Definition INukeOset.h:53
double getAbsorptionFraction() const
return fraction of absorption events
Definition INukeOset.h:59
INukeOset * currentInstance
Definition INukeOset.cxx:5

References osetUtils::currentInstance, genie::Intranuke2025::fNucAbsFac, genie::Intranuke2025::fNucCEXFac, fNucQEFac, genie::RandomGen::Instance(), genie::kIHNFtAbs, genie::kIHNFtCEx, genie::kIHNFtElas, LOG, pNOTICE, and genie::RandomGen::RndFsi().

Referenced by HadronFateHN().

◆ HandleCompoundNucleus()

int HNIntranuke2025::HandleCompoundNucleus ( GHepRecord * ev,
GHepParticle * p,
int mom ) const
privatevirtual

Implements genie::Intranuke2025.

Definition at line 905 of file HNIntranuke2025.cxx.

906{
907
908 // handle compound nucleus option
909 // -- Call the PreEquilibrium function
911 { // random number generator
912 //unused var - quiet compiler warning//RandomGen * rnd = RandomGen::Instance();
913
914 if((p->KinE() < fEPreEq) )
915 {
916 if(fRemnA>4) //this needs to be matched to what is in PreEq and Eq
917 {
918 GHepParticle * sp = new GHepParticle(*p);
919 sp->SetFirstMother(mom);
920 // this was PreEquilibrium - now just used for hN
921 //same arguement lists for PreEq and Eq
924
925 delete sp;
926 return 2;
927 }
928 else
929 {
930 // nothing left to interact with!
931 LOG("HNIntranuke2025", pNOTICE)
932 << "*** Nothing left to interact with, escaping.";
933 GHepParticle * sp = new GHepParticle(*p);
934 sp->SetFirstMother(mom);
936 ev->AddParticle(*sp);
937 delete sp;
938 return 1;
939 }
940 }
941 }
942 return 0;
943}
double fNucRmvE
binding energy to subtract from cascade nucleons
bool IsInNucleus(const GHepParticle *p) const
double fEPreEq
threshold for pre-equilibrium reaction
bool IsNeutronOrProton(int pdgc)
Definition PDGUtils.cxx:351
void Equilibrium(GHepRecord *ev, GHepParticle *p, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, const NuclearModelI *Nuclmodel, double NucRmvE, EINukeMode mode=kIMdHN)

References genie::GHepRecord::AddParticle(), genie::utils::intranuke2025::Equilibrium(), genie::Intranuke2025::fDoCompoundNucleus, genie::Intranuke2025::fDoFermi, genie::Intranuke2025::fEPreEq, genie::Intranuke2025::fFermiFac, genie::Intranuke2025::fNuclmodel, genie::Intranuke2025::fNucRmvE, genie::Intranuke2025::fRemnA, genie::Intranuke2025::fRemnP4, genie::Intranuke2025::fRemnZ, genie::Intranuke2025::IsInNucleus(), genie::pdg::IsNeutronOrProton(), genie::kIMdHN, genie::GHepParticle::KinE(), genie::kIStStableFinalState, LOG, genie::GHepParticle::Pdg(), pNOTICE, genie::GHepParticle::SetFirstMother(), and genie::GHepParticle::SetStatus().

Referenced by HandleCompoundNucleusHN().

◆ HandleCompoundNucleusHN()

bool HNIntranuke2025::HandleCompoundNucleusHN ( GHepRecord * ev,
GHepParticle * p ) const
private

Definition at line 945 of file HNIntranuke2025.cxx.

946{
947 return (this->HandleCompoundNucleus(ev,p,p->FirstMother())==2);
948}
int HandleCompoundNucleus(GHepRecord *ev, GHepParticle *p, int mom) const

References genie::GHepParticle::FirstMother(), and HandleCompoundNucleus().

◆ InelasticHN()

void HNIntranuke2025::InelasticHN ( GHepRecord * ev,
GHepParticle * p ) const
private

Definition at line 765 of file HNIntranuke2025.cxx.

766{
767 // Aaron Meyer (Jan 2010)
768 // Updated version of InelasticHN
769
770 GHepParticle s1(*p);
771 GHepParticle s2(*p);
772 GHepParticle s3(*p);
773 s2.SetRemovalEnergy(0);
774 s3.SetRemovalEnergy(0);
775
776
777
779 {
780 // set status of particles and return
781
782 s1.SetStatus(kIStHadronInTheNucleus);
783 s2.SetStatus(kIStHadronInTheNucleus);
784 s3.SetStatus(kIStHadronInTheNucleus);
785
786 ev->AddParticle(s1);
787 ev->AddParticle(s2);
788 ev->AddParticle(s3);
789 }
790 else
791 {
792 LOG("HNIntranuke2025", pNOTICE) << "Error: could not create pion production final state";
793 exceptions::INukeException exception;
794 exception.SetReason("PionProduction in hN failed");
795 throw exception;
796 }
797 return;
798
799}
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)

References genie::GHepRecord::AddParticle(), genie::Intranuke2025::fDoFermi, genie::Intranuke2025::fFermiFac, genie::Intranuke2025::fFermiMomentum, genie::Intranuke2025::fNuclmodel, genie::Intranuke2025::fRemnA, genie::Intranuke2025::fRemnP4, genie::Intranuke2025::fRemnZ, genie::kIStHadronInTheNucleus, LOG, genie::utils::intranuke2025::PionProduction(), pNOTICE, genie::exceptions::INukeException::SetReason(), genie::GHepParticle::SetRemovalEnergy(), and genie::GHepParticle::SetStatus().

Referenced by SimulateHadronicFinalState().

◆ LoadConfig()

void HNIntranuke2025::LoadConfig ( void )
privatevirtual

Implements genie::Intranuke2025.

Definition at line 951 of file HNIntranuke2025.cxx.

952{
953 // load hadronic cross sections
955
956 // fermi momentum setup
957 // this is specifically set in Intranuke2025::Configure(string)
958 fNuclmodel = dynamic_cast<const NuclearModelI *>( this -> SubAlg("NuclearModel") ) ;
959
960 // other intranuke config params
961 GetParam( "NUCL-R0", fR0 ); // fm
962 GetParam( "NUCL-NR", fNR );
963
964 GetParam( "INUKE-NucRemovalE", fNucRmvE ); // GeV
965 GetParam( "INUKE-HadStep", fHadStep ) ;
966 GetParam( "INUKE-NucAbsFac", fNucAbsFac ) ;
967 GetParam( "INUKE-NucQEFac", fNucQEFac ) ;
968 GetParam( "INUKE-NucCEXFac", fNucCEXFac ) ;
969 GetParam( "INUKE-Energy_Pre_Eq", fEPreEq ) ;
970 GetParam( "INUKE-FermiFac", fFermiFac ) ;
971 GetParam( "INUKE-FermiMomentum", fFermiMomentum ) ;
972
973 GetParam( "INUKE-DoCompoundNucleus", fDoCompoundNucleus ) ;
974 GetParam( "INUKE-DoFermi", fDoFermi ) ;
975 GetParam( "INUKE-XsecNNCorr", fXsecNNCorr ) ;
976 GetParamDef( "AltOset", fAltOset, false ) ;
977
978 GetParam( "HNINUKE-UseOset", fUseOset ) ;
979 GetParam( "HNINUKE-DelRPion", fDelRPion ) ;
980 GetParam( "HNINUKE-DelRNucleon", fDelRNucleon ) ;
981
982 GetParamDef( "FSI-ChargedPion-MFPScale", fChPionMFPScale, 1.0 ) ;
983 GetParamDef( "FSI-NeutralPion-MFPScale", fNeutralPionMFPScale, 1.0 ) ;
984 GetParamDef( "FSI-Nucleon-MFPScale", fNucleonMFPScale, 1.0 ) ;
985
986 // report
987 LOG("HNIntranuke2025", pINFO) << "Settings for Intranuke2025 mode: " << INukeMode::AsString(kIMdHN);
988 LOG("HNIntranuke2025", pWARN) << "R0 = " << fR0 << " fermi";
989 LOG("HNIntranuke2025", pWARN) << "NR = " << fNR;
990 LOG("HNIntranuke2025", pWARN) << "DelRPion = " << fDelRPion;
991 LOG("HNIntranuke2025", pWARN) << "DelRNucleon = " << fDelRNucleon;
992 LOG("HNIntranuke2025", pWARN) << "HadStep = " << fHadStep << " fermi";
993 LOG("HNIntranuke2025", pWARN) << "NucAbsFac = " << fNucAbsFac;
994 LOG("HNIntranuke2025", pWARN) << "NucQEFac = " << fNucQEFac;
995 LOG("HNIntranuke2025", pWARN) << "NucCEXFac = " << fNucCEXFac;
996 LOG("HNIntranuke2025", pWARN) << "EPreEq = " << fEPreEq;
997 LOG("HNIntranuke2025", pWARN) << "FermiFac = " << fFermiFac;
998 LOG("HNIntranuke2025", pWARN) << "FermiMomtm = " << fFermiMomentum;
999 LOG("HNIntranuke2025", pWARN) << "DoFermi? = " << ((fDoFermi)?(true):(false));
1000 LOG("HNIntranuke2025", pWARN) << "DoCmpndNuc? = " << ((fDoCompoundNucleus)?(true):(false));
1001 LOG("HNIntranuke2025", pWARN) << "useOset = " << fUseOset;
1002 LOG("HNIntranuke2025", pWARN) << "altOset = " << fAltOset;
1003 LOG("HNIntranuke2025", pWARN) << "XsecNNCorr? = " << ((fXsecNNCorr)?(true):(false));
1004 LOG("HNIntranuke2025", pWARN) << "FSI-ChargedPion-MFPScale = " << fChPionMFPScale;
1005 LOG("HNIntranuke2025", pWARN) << "FSI-NeutralPion-MFPScale = " << fNeutralPionMFPScale;
1006}
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Algorithm * SubAlg(const RgKey &registry_key) const
static INukeHadroData2025 * Instance(void)
static string AsString(INukeMode_t mode)
Definition INukeMode.h:41
double fChPionMFPScale
tweaking factors for tuning
double fDelRPion
factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement
double fHadStep
step size for intranuclear hadron transport
double fNR
param multiplying the nuclear radius, determining how far to track hadrons beyond the "nuclear bounda...
bool fAltOset
NuWro's table-based implementation (not recommended)
bool fXsecNNCorr
use nuclear medium correction for NN cross section
double fDelRNucleon
factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement
double fR0
effective nuclear size param

References genie::INukeMode::AsString(), genie::Intranuke2025::fAltOset, genie::Intranuke2025::fChPionMFPScale, genie::Intranuke2025::fDelRNucleon, genie::Intranuke2025::fDelRPion, genie::Intranuke2025::fDoCompoundNucleus, genie::Intranuke2025::fDoFermi, genie::Intranuke2025::fEPreEq, genie::Intranuke2025::fFermiFac, genie::Intranuke2025::fFermiMomentum, genie::Intranuke2025::fHadroData2025, genie::Intranuke2025::fHadStep, genie::Intranuke2025::fNeutralPionMFPScale, genie::Intranuke2025::fNR, genie::Intranuke2025::fNucAbsFac, genie::Intranuke2025::fNucCEXFac, genie::Intranuke2025::fNucleonMFPScale, genie::Intranuke2025::fNuclmodel, fNucQEFac, genie::Intranuke2025::fNucRmvE, genie::Intranuke2025::fR0, genie::Intranuke2025::fUseOset, genie::Intranuke2025::fXsecNNCorr, genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), genie::INukeHadroData2025::Instance(), genie::kIMdHN, LOG, pINFO, pWARN, and genie::Algorithm::SubAlg().

◆ ProcessEventRecord()

void HNIntranuke2025::ProcessEventRecord ( GHepRecord * event_rec) const
virtual

Reimplemented from genie::Intranuke2025.

Definition at line 114 of file HNIntranuke2025.cxx.

115{
116 LOG("HNIntranuke2025", pNOTICE)
117 << "************ Running hN2025 MODE INTRANUKE ************";
118
120
121 LOG("HNIntranuke2025", pINFO) << "Done with this event";
122}
virtual void ProcessEventRecord(GHepRecord *event_rec) const

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

◆ SimulateHadronicFinalState()

void HNIntranuke2025::SimulateHadronicFinalState ( GHepRecord * ev,
GHepParticle * p ) const
privatevirtual

Implements genie::Intranuke2025.

Definition at line 124 of file HNIntranuke2025.cxx.

125{
126// Simulate a hadron interaction for the input particle p in HN mode
127//
128 if(!p || !ev)
129 {
130 LOG("HNIntranuke2025", pERROR) << "** Null input!";
131 return;
132 }
133
134 // check particle id
135 int pdgc = p->Pdg();
136 bool is_pion = (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0);
137 bool is_kaon = (pdgc==kPdgKP);
138 bool is_baryon = (pdgc==kPdgProton || pdgc==kPdgNeutron);
139 bool is_gamma = (pdgc==kPdgGamma);
140 if(!(is_pion || is_baryon || is_gamma || is_kaon))
141 {
142 LOG("HNIntranuke2025", pERROR) << "** Cannot handle particle: " << p->Name();
143 return;
144 }
145 try
146 {
147 // select a fate for the input particle
148 INukeFateHN_t fate = this->HadronFateHN(p);
149
150 // store the fate
151 ev->Particle(p->FirstMother())->SetRescatterCode((int)fate);
152
153 if(fate == kIHNFtUndefined)
154 {
155 LOG("HNIntranuke2025", pERROR) << "** Couldn't select a fate";
156 LOG("HNIntranuke2025", pERROR) << "** Num Protons: " << fRemnZ
157 << ", Num Neutrons: "<<(fRemnA-fRemnZ);
158 LOG("HNIntranuke2025", pERROR) << "** Particle: " << "\n" << (*p);
159 //LOG("HNIntranuke2025", pERROR) << "** Event Record: " << "\n" << (*ev);
160 //p->SetStatus(kIStUndefined);
161 p->SetStatus(kIStStableFinalState);
162 ev->AddParticle(*p);
163 return;
164 }
165
166 LOG("HNIntranuke2025", pNOTICE)
167 << "Selected " << p->Name() << " fate: " << INukeHadroFates2025::AsString(fate);
168
169 // handle the reaction
170 if(fate == kIHNFtCEx || fate == kIHNFtElas)
171 {
172 this->ElasHN(ev,p,fate);
173 }
174 else if(fate == kIHNFtAbs) {this-> AbsorbHN(ev,p,fate);}
175 else if(fate == kIHNFtInelas && pdgc != kPdgGamma)
176 {
177#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
178 LOG("HNIntranuke2025", pDEBUG)
179 << "Invoking InelasticHN() for a : " << p->Name()
180 << " whose fate is : " << INukeHadroFates2025::AsString(fate);
181#endif
182
183 this-> InelasticHN(ev,p);
184 }
185 else if(fate == kIHNFtInelas && pdgc == kPdgGamma) {this-> GammaInelasticHN(ev,p,fate);}
186 else if(fate == kIHNFtCmp) LOG("HNIntranuke2025", pFATAL) << "The PreEquilibrium and Compound Nucleus code are experimental, should not be used";
187//{utils::intranuke2025::PreEquilibrium(ev,p,fRemnA,fRemnZ,fRemnP4,fDoFermi,fFermiFac,fNuclmodel,fNucRmvE,kIMdHN);}
188 else if(fate == kIHNFtNoInteraction)
189 {
190 p->SetStatus(kIStStableFinalState);
191 ev->AddParticle(*p);
192 return;
193 }
194 }
195 catch(exceptions::INukeException exception)
196 {
197 this->SimulateHadronicFinalState(ev,p);
198 LOG("HNIntranuke2025", pNOTICE)
199 << "retry call to SimulateHadronicFinalState ";
200 LOG("HNIntranuke2025", pNOTICE) << exception;
201
202 }
203}
#define pFATAL
Definition Messenger.h:56
void SetRescatterCode(int code)
virtual GHepParticle * Particle(int position) const
INukeFateHN_t HadronFateHN(const GHepParticle *p) const
void SimulateHadronicFinalState(GHepRecord *ev, GHepParticle *p) const
void ElasHN(GHepRecord *ev, GHepParticle *p, INukeFateHN_t fate) const
void AbsorbHN(GHepRecord *ev, GHepParticle *p, INukeFateHN_t fate) const
void InelasticHN(GHepRecord *ev, GHepParticle *p) const
void GammaInelasticHN(GHepRecord *ev, GHepParticle *p, INukeFateHN_t fate) const
@ kIHNFtNoInteraction
enum genie::EINukeFateHN_t INukeFateHN_t

References AbsorbHN(), genie::GHepRecord::AddParticle(), genie::INukeHadroFates2025::AsString(), ElasHN(), genie::GHepParticle::FirstMother(), genie::Intranuke2025::fRemnA, genie::Intranuke2025::fRemnZ, GammaInelasticHN(), HadronFateHN(), InelasticHN(), genie::kIHNFtAbs, genie::kIHNFtCEx, genie::kIHNFtCmp, genie::kIHNFtElas, genie::kIHNFtInelas, genie::kIHNFtNoInteraction, genie::kIHNFtUndefined, genie::kIStStableFinalState, genie::kPdgGamma, genie::kPdgKP, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, LOG, genie::GHepParticle::Name(), genie::GHepRecord::Particle(), pDEBUG, genie::GHepParticle::Pdg(), pERROR, pFATAL, pNOTICE, genie::GHepParticle::SetRescatterCode(), genie::GHepParticle::SetStatus(), and SimulateHadronicFinalState().

Referenced by SimulateHadronicFinalState().

◆ IntranukeTester

friend class IntranukeTester
friend

Definition at line 53 of file HNIntranuke2025.h.

References IntranukeTester.

Referenced by IntranukeTester.

Member Data Documentation

◆ fNucQEFac

double genie::HNIntranuke2025::fNucQEFac
private

Definition at line 84 of file HNIntranuke2025.h.

Referenced by HadronFateHN(), HadronFateOset(), and LoadConfig().

◆ nuclA

int genie::HNIntranuke2025::nuclA
mutableprivate

value of A for the target nucleus in hA mode

Definition at line 81 of file HNIntranuke2025.h.


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