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

#include <HNIntranuke2018.h>

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

Public Member Functions

 HNIntranuke2018 ()
 HNIntranuke2018 (string config)
 ~HNIntranuke2018 ()
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 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::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)
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::Intranuke2018
double fTrackingRadius
 tracking radius for the nucleus in the current event
TGenPhaseSpace fGenPhaseSpace
 a phase space generator
INukeHadroData2018fHadroData2018
 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 HNIntranuke2018.h.

Constructor & Destructor Documentation

◆ HNIntranuke2018() [1/2]

HNIntranuke2018::HNIntranuke2018 ( )

Definition at line 96 of file HNIntranuke2018.cxx.

96 :
97Intranuke2018("genie::HNIntranuke2018")
98{
99
100}

References genie::Intranuke2018::Intranuke2018().

◆ HNIntranuke2018() [2/2]

HNIntranuke2018::HNIntranuke2018 ( string config)

Definition at line 102 of file HNIntranuke2018.cxx.

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

References genie::Intranuke2018::Intranuke2018().

◆ ~HNIntranuke2018()

HNIntranuke2018::~HNIntranuke2018 ( )

Definition at line 108 of file HNIntranuke2018.cxx.

109{
110
111}

Member Function Documentation

◆ AbsorbHN()

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

Definition at line 388 of file HNIntranuke2018.cxx.

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

Definition at line 659 of file HNIntranuke2018.cxx.

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

Referenced by SimulateHadronicFinalState().

◆ FateWeight()

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

Definition at line 363 of file HNIntranuke2018.cxx.

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

References genie::Intranuke2018::fDoCompoundNucleus, genie::Intranuke2018::fRemnA, genie::Intranuke2018::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 HNIntranuke2018::GammaInelasticHN ( GHepRecord * ev,
GHepParticle * p,
INukeFateHN_t fate ) const
private

Definition at line 804 of file HNIntranuke2018.cxx.

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

References genie::GHepRecord::AddParticle(), genie::INukeHadroFates::AsString(), genie::Intranuke2018::fDoFermi, genie::Intranuke2018::fFermiFac, genie::Intranuke2018::fHadroData2018, genie::Intranuke2018::fNuclmodel, genie::Intranuke2018::fRemnA, genie::Intranuke2018::fRemnP4, genie::Intranuke2018::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::intranuke2018::TwoBodyCollision().

Referenced by SimulateHadronicFinalState().

◆ GetGenINukeMode()

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

Reimplemented from genie::Intranuke2018.

Definition at line 63 of file HNIntranuke2018.h.

63{return "hN";};

◆ GetINukeMode()

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

Reimplemented from genie::Intranuke2018.

Definition at line 62 of file HNIntranuke2018.h.

62{return "hN2018";};

◆ HadronFateHN()

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

Definition at line 208 of file HNIntranuke2018.cxx.

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

Definition at line 1012 of file HNIntranuke2018.cxx.

1013{
1014 //LOG("HNIntranuke2018", pWARN) << "IN HadronFateOset";
1015
1016 //LOG("HNIntranuke2018", pWARN) << "{ frac abs = " << osetUtils::currentInstance->getAbsorptionFraction();
1017 //LOG("HNIntranuke2018", pWARN) << " frac cex = " << osetUtils::currentInstance->getCexFraction() << " }";
1018
1019 double fractionAbsorption = osetUtils::currentInstance->getAbsorptionFraction();
1020 double fractionCex = osetUtils::currentInstance->getCexFraction ();
1021 double fractionElas = 1 - (fractionAbsorption + fractionCex);
1022
1023 fractionCex *= fNucCEXFac; // scaling factors
1024 fractionAbsorption *= fNucAbsFac;
1025 fractionElas *= fNucQEFac;
1026
1027 double totalFrac = fractionCex + fractionAbsorption + fractionElas;
1028
1029 RandomGen *randomGenerator = RandomGen::Instance();
1030 const double randomNumber = randomGenerator->RndFsi().Rndm() * totalFrac;
1031
1032 LOG("HNIntranuke2018", pNOTICE) << "{ frac abs = " << fractionAbsorption;
1033 LOG("HNIntranuke2018", pNOTICE) << " frac cex = " << fractionCex;
1034 LOG("HNIntranuke2018", pNOTICE) << " frac elas = " << fractionElas << " }";
1035
1036 if (randomNumber < fractionAbsorption && fRemnA > 1) return kIHNFtAbs;
1037 else if (randomNumber < fractionAbsorption + fractionCex) return kIHNFtCEx;
1038 else return kIHNFtElas;
1039}
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::Intranuke2018::fNucAbsFac, genie::Intranuke2018::fNucCEXFac, fNucQEFac, genie::RandomGen::Instance(), genie::kIHNFtAbs, genie::kIHNFtCEx, genie::kIHNFtElas, LOG, pNOTICE, and genie::RandomGen::RndFsi().

Referenced by HadronFateHN().

◆ HandleCompoundNucleus()

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

Implements genie::Intranuke2018.

Definition at line 908 of file HNIntranuke2018.cxx.

909{
910
911 // handle compound nucleus option
912 // -- Call the PreEquilibrium function
914 { // random number generator
915 //unused var - quiet compiler warning//RandomGen * rnd = RandomGen::Instance();
916
917 if((p->KinE() < fEPreEq) )
918 {
919 if(fRemnA>4) //this needs to be matched to what is in PreEq and Eq
920 {
921 GHepParticle * sp = new GHepParticle(*p);
922 sp->SetFirstMother(mom);
923 // this was PreEquilibrium - now just used for hN
924 //same arguement lists for PreEq and Eq
927
928 delete sp;
929 return 2;
930 }
931 else
932 {
933 // nothing left to interact with!
934 LOG("HNIntranuke2018", pNOTICE)
935 << "*** Nothing left to interact with, escaping.";
936 GHepParticle * sp = new GHepParticle(*p);
937 sp->SetFirstMother(mom);
939 ev->AddParticle(*sp);
940 delete sp;
941 return 1;
942 }
943 }
944 }
945 return 0;
946}
double fNucRmvE
binding energy to subtract from cascade nucleons
double fEPreEq
threshold for pre-equilibrium reaction
bool IsInNucleus(const GHepParticle *p) const
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::intranuke2018::Equilibrium(), genie::Intranuke2018::fDoCompoundNucleus, genie::Intranuke2018::fDoFermi, genie::Intranuke2018::fEPreEq, genie::Intranuke2018::fFermiFac, genie::Intranuke2018::fNuclmodel, genie::Intranuke2018::fNucRmvE, genie::Intranuke2018::fRemnA, genie::Intranuke2018::fRemnP4, genie::Intranuke2018::fRemnZ, genie::Intranuke2018::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 HNIntranuke2018::HandleCompoundNucleusHN ( GHepRecord * ev,
GHepParticle * p ) const
private

Definition at line 948 of file HNIntranuke2018.cxx.

949{
950 return (this->HandleCompoundNucleus(ev,p,p->FirstMother())==2);
951}
int HandleCompoundNucleus(GHepRecord *ev, GHepParticle *p, int mom) const

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

◆ InelasticHN()

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

Definition at line 768 of file HNIntranuke2018.cxx.

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

Referenced by SimulateHadronicFinalState().

◆ LoadConfig()

void HNIntranuke2018::LoadConfig ( void )
privatevirtual

Implements genie::Intranuke2018.

Definition at line 954 of file HNIntranuke2018.cxx.

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

References genie::INukeMode::AsString(), genie::Intranuke2018::fAltOset, 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::fNeutralPionMFPScale, genie::Intranuke2018::fNR, genie::Intranuke2018::fNucAbsFac, genie::Intranuke2018::fNucCEXFac, genie::Intranuke2018::fNucleonMFPScale, genie::Intranuke2018::fNuclmodel, fNucQEFac, genie::Intranuke2018::fNucRmvE, genie::Intranuke2018::fR0, genie::Intranuke2018::fUseOset, genie::Intranuke2018::fXsecNNCorr, genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), genie::INukeHadroData2018::Instance(), genie::kIMdHN, LOG, pINFO, pWARN, and genie::Algorithm::SubAlg().

◆ ProcessEventRecord()

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

Reimplemented from genie::Intranuke2018.

Definition at line 113 of file HNIntranuke2018.cxx.

114{
115 LOG("HNIntranuke2018", pNOTICE)
116 << "************ Running hN2018 MODE INTRANUKE ************";
117
118 /* LOG("HNIntranuke2018", pWARN)
119 << print::PrintFramedMesg(
120 "Experimental code (INTRANUKE/hN model) - Run at your own risk");
121 */
122
124
125 LOG("HNIntranuke2018", pINFO) << "Done with this event";
126}
virtual void ProcessEventRecord(GHepRecord *event_rec) const

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

◆ SimulateHadronicFinalState()

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

Implements genie::Intranuke2018.

Definition at line 128 of file HNIntranuke2018.cxx.

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

References AbsorbHN(), genie::GHepRecord::AddParticle(), genie::INukeHadroFates::AsString(), ElasHN(), genie::Intranuke2018::fDoFermi, genie::Intranuke2018::fFermiFac, genie::GHepParticle::FirstMother(), genie::Intranuke2018::fNuclmodel, genie::Intranuke2018::fNucRmvE, genie::Intranuke2018::fRemnA, genie::Intranuke2018::fRemnP4, genie::Intranuke2018::fRemnZ, GammaInelasticHN(), HadronFateHN(), InelasticHN(), genie::kIHNFtAbs, genie::kIHNFtCEx, genie::kIHNFtCmp, genie::kIHNFtElas, genie::kIHNFtInelas, genie::kIHNFtNoInteraction, genie::kIHNFtUndefined, genie::kIMdHN, 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, pNOTICE, genie::utils::intranuke2018::PreEquilibrium(), genie::GHepParticle::SetRescatterCode(), genie::GHepParticle::SetStatus(), and SimulateHadronicFinalState().

Referenced by SimulateHadronicFinalState().

◆ IntranukeTester

friend class IntranukeTester
friend

Definition at line 53 of file HNIntranuke2018.h.

References IntranukeTester.

Referenced by IntranukeTester.

Member Data Documentation

◆ fNucQEFac

double genie::HNIntranuke2018::fNucQEFac
private

Definition at line 84 of file HNIntranuke2018.h.

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

◆ nuclA

int genie::HNIntranuke2018::nuclA
mutableprivate

value of A for the target nucleus in hA mode

Definition at line 81 of file HNIntranuke2018.h.


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