GENIEGenerator
Loading...
Searching...
No Matches
Public Member Functions | Private Member Functions | Private Attributes | Friends | List of all members
genie::HAIntranuke2025 Class Reference

#include <HAIntranuke2025.h>

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

Public Member Functions

 HAIntranuke2025 ()
 HAIntranuke2025 (string config)
 ~HAIntranuke2025 ()
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 Registry & GetConfig (void) const
Registry * GetOwnedConfig (void)
virtual const AlgId & Id (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 Algorithm * SubAlg (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
void SimulateHadronicFinalStateKinematics (GHepRecord *ev, GHepParticle *p) const
INukeFateHA_t HadronFateHA (const GHepParticle *p) const
void Inelastic (GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
void ElasHA (GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
void InelasticHA (GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
double PiBounce (void) const
double PnBounce (void) const
int HandleCompoundNucleus (GHepRecord *ev, GHepParticle *p, int mom) const

Private Attributes

int nuclA
 value of A for the target nucleus in hA mode
unsigned int fNumIterations

Friends

class IntranukeTester

Additional Inherited Members

Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
static string BuildParamVectSizeKey (const std::string &comm_name)
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
static string BuildParamMatRowSizeKey (const std::string &comm_name)
static string BuildParamMatColSizeKey (const std::string &comm_name)
Protected Member Functions inherited from genie::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)
Registry * ExtractLocalConfig (const Registry &in) const
Registry * ExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key.
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
INukeHadroData2025 * fHadroData2025
 a collection of h+N,h+A data & calculations
AlgFactory * fAlgf
 algorithm factory instance
const NuclearModelI * fNuclmodel
 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
AlgMap * fOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool)

Detailed Description

Definition at line 55 of file HAIntranuke2025.h.

Constructor & Destructor Documentation

◆ HAIntranuke2025() [1/2]

HAIntranuke2025::HAIntranuke2025 ( )

Definition at line 74 of file HAIntranuke2025.cxx.

74 :
75Intranuke2025("genie::HAIntranuke2025")
76{
77
78}

References genie::Intranuke2025::Intranuke2025().

◆ HAIntranuke2025() [2/2]

HAIntranuke2025::HAIntranuke2025 ( string config)

Definition at line 80 of file HAIntranuke2025.cxx.

80 :
81Intranuke2025("genie::HAIntranuke2025",config)
82{
83
84}

References genie::Intranuke2025::Intranuke2025().

◆ ~HAIntranuke2025()

HAIntranuke2025::~HAIntranuke2025 ( )

Definition at line 86 of file HAIntranuke2025.cxx.

87{
88
89}

Member Function Documentation

◆ ElasHA()

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

Definition at line 493 of file HAIntranuke2025.cxx.

495{
496 // scatters particle within nucleus, copy of hN code meant to run only once
497 // in hA mode
498
499 LOG("HAIntranuke2025", pDEBUG)
500 << "ElasHA() is invoked for a : " << p->Name()
501 << " whose fate is : " << INukeHadroFates2025::AsString(fate);
502
503 /* if(fate!=kIHAFtElas)
504 {
505 LOG("HAIntranuke2025", pWARN)
506 << "ElasHA() cannot handle fate: " << INukeHadroFates::AsString(fate);
507 return;
508 } */
509
510 // check remnants
511 if(fRemnA<0 || fRemnZ<0) // best to stop it here and not try again.
512 {
513 LOG("HAIntranuke2025", pWARN) << "Invalid Nucleus! : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
514 p->SetStatus(kIStStableFinalState);
515 ev->AddParticle(*p);
516 return;
517 }
518
519 // vars for incoming particle, target, and scattered pdg codes
520 int pcode = p->Pdg();
521 double Mp = p->Mass();
522 double Mt = 0.;
523 if (ev->TargetNucleus()->A()==fRemnA)
524 { Mt = PDGLibrary::Instance()->Find(ev->TargetNucleus()->Pdg())->Mass(); }
525 else
526 {
527 Mt = fRemnP4.M();
528 }
529 TLorentzVector t4PpL = *p->P4();
530 TLorentzVector t4PtL = fRemnP4;
531 double C3CM = 0.0;
532
533 // calculate scattering angle
534 if(pcode==kPdgNeutron||pcode==kPdgProton) C3CM = TMath::Cos(this->PnBounce());
535 else C3CM = TMath::Cos(this->PiBounce());
536
537 // calculate final 4 momentum of probe
538 TLorentzVector t4P3L, t4P4L;
539
540 if (!utils::intranuke2025::TwoBodyKinematics(Mp,Mt,t4PpL,t4PtL,t4P3L,t4P4L,C3CM,fRemnP4))
541 {
542 LOG("HAIntranuke2025", pNOTICE) << "ElasHA() failed";
543 exceptions::INukeException exception;
544 exception.SetReason("TwoBodyKinematics failed in ElasHA");
545 throw exception;
546 }
547
548 // Update probe particle
549 p->SetMomentum(t4P3L);
550 p->SetStatus(kIStStableFinalState);
551
552 // Update Remnant nucleus
553 fRemnP4 = t4P4L;
554 LOG("HAIntranuke2025",pINFO)
555 << "C3cm = " << C3CM;
556 LOG("HAIntranuke2025",pINFO)
557 << "|p3| = " << t4P3L.Vect().Mag() << ", E3 = " << t4P3L.E() << ",Mp = " << Mp;
558 LOG("HAIntranuke2025",pINFO)
559 << "|p4| = " << fRemnP4.Vect().Mag() << ", E4 = " << fRemnP4.E() << ",Mt = " << Mt;
560
561 ev->AddParticle(*p);
562
563}
pNOTICE
#define pNOTICE
Definition Messenger.h:61
pINFO
#define pINFO
Definition Messenger.h:62
pDEBUG
#define pDEBUG
Definition Messenger.h:63
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
pWARN
#define pWARN
Definition Messenger.h:60
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition GHepParticle.cxx:135
genie::GHepParticle::SetMomentum
void SetMomentum(const TLorentzVector &p4)
Definition GHepParticle.cxx:236
genie::GHepParticle::Pdg
int Pdg(void) const
Definition GHepParticle.h:63
genie::GHepParticle::P4
const TLorentzVector * P4(void) const
Definition GHepParticle.h:78
genie::GHepParticle::Mass
double Mass(void) const
Mass that corresponds to the PDG code.
Definition GHepParticle.cxx:143
genie::GHepParticle::SetStatus
void SetStatus(GHepStatus_t s)
Definition GHepParticle.h:126
genie::GHepParticle::A
int A(void) const
Definition GHepParticle.cxx:184
genie::GHepRecord::TargetNucleus
virtual GHepParticle * TargetNucleus(void) const
Definition GHepRecord.cxx:293
genie::GHepRecord::AddParticle
virtual void AddParticle(const GHepParticle &p)
Definition GHepRecord.cxx:499
genie::HAIntranuke2025::PiBounce
double PiBounce(void) const
Definition HAIntranuke2025.cxx:388
genie::HAIntranuke2025::PnBounce
double PnBounce(void) const
Definition HAIntranuke2025.cxx:440
genie::INukeHadroFates2025::AsString
static string AsString(INukeFateHN_t fate)
Definition INukeHadroFates2025.h:71
genie::Intranuke2025::fRemnA
int fRemnA
remnant nucleus A
Definition Intranuke2025.h:121
genie::Intranuke2025::fRemnZ
int fRemnZ
remnant nucleus Z
Definition Intranuke2025.h:122
genie::Intranuke2025::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition Intranuke2025.h:123
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition PDGLibrary.cxx:68
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition PDGLibrary.cxx:86
genie::utils::intranuke2025::TwoBodyKinematics
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition INukeUtils2025.cxx:858
genie::kIStStableFinalState
@ kIStStableFinalState
Definition GHepStatus.h:30
genie::kPdgProton
const int kPdgProton
Definition PDGCodes.h:81
genie::kPdgNeutron
const int kPdgNeutron
Definition PDGCodes.h:83

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

◆ GetGenINukeMode()

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

Reimplemented from genie::Intranuke2025.

Definition at line 67 of file HAIntranuke2025.h.

67{return "hA";};

◆ GetINukeMode()

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

Reimplemented from genie::Intranuke2025.

Definition at line 66 of file HAIntranuke2025.h.

66{return "hA2025";};

◆ HadronFateHA()

INukeFateHA_t HAIntranuke2025::HadronFateHA ( const GHepParticle * p) const
private

Definition at line 203 of file HAIntranuke2025.cxx.

204{
205// Select a hadron fate in HA mode
206//
207 RandomGen * rnd = RandomGen::Instance();
208
209 // get pdgc code & kinetic energy in MeV
210 int pdgc = p->Pdg();
211 double ke = p->KinE() / units::MeV;
212
213 //bool isPion = (pdgc == kPdgPiP or pdgc == kPdgPi0 or pdgc == kPdgPiM);
214 //if (isPion and fUseOset and ke < 350.0) return this->HadronFateOset();
215
216 LOG("HAIntranuke2025", pINFO)
217 << "Selecting hA fate for " << p->Name() << " with KE = " << ke << " MeV";
218
219 // try to generate a hadron fate
220 unsigned int iter = 0;
221 while(iter++ < kRjMaxIterations) {
222
223 // handle pions
224 //
225 if (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0) {
226
227 double frac_cex = fHadroData2025->FracADep(pdgc, kIHAFtCEx, ke, nuclA);
228 // double frac_elas = fHadroData2025->FracADep(pdgc, kIHAFtElas, ke, nuclA);
229 double frac_inel = fHadroData2025->FracADep(pdgc, kIHAFtInelas, ke, nuclA);
230 double frac_abs = fHadroData2025->FracADep(pdgc, kIHAFtAbs, ke, nuclA);
231 double frac_piprod = fHadroData2025->FracADep(pdgc, kIHAFtPiProd, ke, nuclA);
232 LOG("HAIntranuke2025", pDEBUG)
233 << "\n frac{" << INukeHadroFates2025::AsString(kIHAFtCEx) << "} = " << frac_cex
234 // << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
235 << "\n frac{" << INukeHadroFates2025::AsString(kIHAFtInelas) << "} = " << frac_inel
236 << "\n frac{" << INukeHadroFates2025::AsString(kIHAFtAbs) << "} = " << frac_abs
237 << "\n frac{" << INukeHadroFates2025::AsString(kIHAFtPiProd) << "} = " << frac_piprod;
238
239 // apply external tweaks to fractions
240 frac_cex *= fPionFracCExScale;
241 frac_inel *= fPionFracInelScale;
242 if (pdgc==kPdgPiP || pdgc==kPdgPiM) frac_abs *= fChPionFracAbsScale;
243 if (pdgc==kPdgPi0) frac_abs *= fNeutralPionFracAbsScale;
244 frac_piprod *= fPionFracPiProdScale;
245
246
247 // Flag to enable or disable π0/π+ ratio corrections
248 bool apply_pi0_ratio_correction = true;
249
250 if (apply_pi0_ratio_correction && pdgc == kPdgPi0) {
251
252 // Cap the kinetic energy at 1000 MeV for the correction calculation
253 // decision made by Steve Dytman and Hugh Gallagher (2006) as best for MINOS. Looks OK for DUNE in 1st simulations.
254 // This is reasonable because cross ections for KE above 1 GeV tend to be constand.
255 double ke_ratio = (ke > 1000.0) ? 1000.0 : ke;
256
257 // Define correction scale factors as a function of kinetic energy (in MeV) to match reslts of hN sim. MI Nov, 2025.
258 double ratio_cex = 0.0008702 * ke_ratio + 1.9047;
259 double ratio_abs = 0.0003291 * ke_ratio + 0.82617;
260 double ratio_inel = -0.0003209 * ke_ratio + 0.837764;
261 double ratio_piprod = 0.0004402 * ke_ratio + 0.47418;
262
263 // Apply the corrections only if kinetic energy is below 1000 MeV
264
265 frac_cex *= ratio_cex;
266 frac_abs *= ratio_abs;
267 frac_inel *= ratio_inel;
268
269 // Apply piprod correction only if KE is above 400 MeV
270 if (ke > 400.0) {
271 frac_piprod *= ratio_piprod;
272 }
273
274 }
275
276
277
278
279 double frac_rescale = 1./(frac_cex + frac_inel + frac_abs + frac_piprod);
280
281 frac_cex *= frac_rescale;
282 frac_inel *= frac_rescale;
283 frac_abs *= frac_rescale;
284 frac_piprod *= frac_rescale;
285
286
287 // compute total fraction (can be <1 if fates have been switched off)
288 double tf = frac_cex +
289 // frac_elas +
290 frac_inel +
291 frac_abs +
292 frac_piprod;
293
294 double r = tf * rnd->RndFsi().Rndm();
295#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
296 LOG("HAIntranuke2025", pDEBUG) << "r = " << r << " (max = " << tf << ")";
297#endif
298 double cf=0; // current fraction
299 if(r < (cf += frac_cex )) return kIHAFtCEx; // cex
300 // if(r < (cf += frac_elas )) return kIHAFtElas; // elas
301 if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
302 if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
303 if(r < (cf += frac_piprod )) return kIHAFtPiProd; // pi prod
304
305 LOG("HAIntranuke2025", pWARN)
306 << "No selection after going through all fates! "
307 << "Total fraction = " << tf << " (r = " << r << ")";
308 }
309
310 // handle nucleons
311 else if (pdgc==kPdgProton || pdgc==kPdgNeutron) {
312 double frac_cex = fHadroData2025->FracAIndep(pdgc, kIHAFtCEx, ke);
313 //double frac_elas = fHadroData2025->FracAIndep(pdgc, kIHAFtElas, ke);
314 double frac_inel = fHadroData2025->FracAIndep(pdgc, kIHAFtInelas, ke);
315 double frac_abs = fHadroData2025->FracAIndep(pdgc, kIHAFtAbs, ke);
316 double frac_pipro = fHadroData2025->FracAIndep(pdgc, kIHAFtPiProd, ke);
317 double frac_cmp = fHadroData2025->FracAIndep(pdgc, kIHAFtCmp , ke);
318
319 LOG("HAIntranuke2025", pINFO)
320 << "\n frac{" << INukeHadroFates2025::AsString(kIHAFtCEx) << "} = " << frac_cex
321 // << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
322 << "\n frac{" << INukeHadroFates2025::AsString(kIHAFtInelas) << "} = " << frac_inel
323 << "\n frac{" << INukeHadroFates2025::AsString(kIHAFtAbs) << "} = " << frac_abs
324 << "\n frac{" << INukeHadroFates2025::AsString(kIHAFtPiProd) << "} = " << frac_pipro
325 << "\n frac{" << INukeHadroFates2025::AsString(kIHAFtCmp) << "} = " << frac_cmp; //suarez edit, cmp
326
327 // apply external tweaks to fractions
328 frac_cex *= fNucleonFracCExScale;
329 frac_inel *= fNucleonFracInelScale;
330 frac_abs *= fNucleonFracAbsScale;
331 frac_pipro *= fNucleonFracPiProdScale;
332
333 double frac_rescale = 1./(frac_cex + frac_inel + frac_abs + frac_pipro);
334
335 frac_cex *= frac_rescale;
336 frac_inel *= frac_rescale;
337 frac_abs *= frac_rescale;
338 frac_pipro *= frac_rescale;
339
340 // compute total fraction (can be <1 if fates have been switched off)
341 double tf = frac_cex +
342 //frac_elas +
343 frac_inel +
344 frac_abs +
345 frac_pipro +
346 frac_cmp; //suarez edit, cmp
347
348 double r = tf * rnd->RndFsi().Rndm();
349#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
350 LOG("HAIntranuke2025", pDEBUG) << "r = " << r << " (max = " << tf << ")";
351#endif
352 double cf=0; // current fraction
353 if(r < (cf += frac_cex )) return kIHAFtCEx; // cex
354 //if(r < (cf += frac_elas )) return kIHAFtElas; // elas
355 if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
356 if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
357 if(r < (cf += frac_pipro )) return kIHAFtPiProd; // pi prod
358 if(r < (cf += frac_cmp )) return kIHAFtCmp; //suarez edit, cmp
359
360 LOG("HAIntranuke2025", pWARN)
361 << "No selection after going through all fates! "
362 << "Total fraction = " << tf << " (r = " << r << ")";
363 }
364 // handle kaons SD: only K+, remove K- in following statement
365 else if (pdgc==kPdgKP) {
366 double frac_inel = fHadroData2025->FracAIndep(pdgc, kIHAFtInelas, ke);
367 double frac_abs = fHadroData2025->FracAIndep(pdgc, kIHAFtAbs, ke);
368
369 LOG("HAIntranuke2025", pDEBUG)
370 << "\n frac{" << INukeHadroFates2025::AsString(kIHAFtInelas) << "} = " << frac_inel
371 << "\n frac{" << INukeHadroFates2025::AsString(kIHAFtAbs) << "} = " << frac_abs;
372 // compute total fraction (can be <1 if fates have been switched off)
373 double tf = frac_inel +
374 frac_abs;
375 double r = tf * rnd->RndFsi().Rndm();
376#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
377 LOG("HAIntranuke2025", pDEBUG) << "r = " << r << " (max = " << tf << ")";
378#endif
379 double cf=0; // current fraction
380 if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
381 if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
382 }
383 }//iterations
384
385 return kIHAFtUndefined;
386}
genie::GHepParticle::KinE
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
Definition GHepParticle.cxx:159
genie::HAIntranuke2025::nuclA
int nuclA
value of A for the target nucleus in hA mode
Definition HAIntranuke2025.h:85
genie::Intranuke2025::fHadroData2025
INukeHadroData2025 * fHadroData2025
a collection of h+N,h+A data & calculations
Definition Intranuke2025.h:118
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition RandomGen.cxx:74
genie::RandomGen::RndFsi
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition RandomGen.h:59
genie::controls::kRjMaxIterations
static const unsigned int kRjMaxIterations
Definition Controls.h:26
genie::units::MeV
static constexpr double MeV
Definition Units.h:129
genie::kPdgPiM
const int kPdgPiM
Definition PDGCodes.h:159
genie::kPdgPi0
const int kPdgPi0
Definition PDGCodes.h:160
genie::kPdgKP
const int kPdgKP
Definition PDGCodes.h:172
genie::kIHAFtUndefined
@ kIHAFtUndefined
Definition INukeHadroFates.h:48
genie::kPdgPiP
const int kPdgPiP
Definition PDGCodes.h:158

References genie::INukeHadroFates2025::AsString(), genie::Intranuke2025::fChPionFracAbsScale, genie::Intranuke2025::fHadroData2025, genie::Intranuke2025::fNeutralPionFracAbsScale, genie::Intranuke2025::fNucleonFracAbsScale, genie::Intranuke2025::fNucleonFracCExScale, genie::Intranuke2025::fNucleonFracInelScale, genie::Intranuke2025::fNucleonFracPiProdScale, genie::Intranuke2025::fPionFracCExScale, genie::Intranuke2025::fPionFracInelScale, genie::Intranuke2025::fPionFracPiProdScale, genie::RandomGen::Instance(), genie::kIHAFtAbs, genie::kIHAFtCEx, genie::kIHAFtCmp, genie::kIHAFtInelas, genie::kIHAFtPiProd, genie::kIHAFtUndefined, genie::GHepParticle::KinE(), genie::kPdgKP, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, genie::controls::kRjMaxIterations, LOG, genie::units::MeV, genie::GHepParticle::Name(), nuclA, pDEBUG, genie::GHepParticle::Pdg(), pINFO, pWARN, and genie::RandomGen::RndFsi().

Referenced by SimulateHadronicFinalState().

◆ HandleCompoundNucleus()

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

Implements genie::Intranuke2025.

Definition at line 1530 of file HAIntranuke2025.cxx.

1531{
1532
1533 // only relevant for hN mode - not anymore.
1534 return false;
1535
1536}

◆ Inelastic()

void HAIntranuke2025::Inelastic ( GHepRecord * ev,
GHepParticle * p,
INukeFateHA_t fate ) const
private

Definition at line 751 of file HAIntranuke2025.cxx.

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

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

Referenced by SimulateHadronicFinalStateKinematics().

◆ InelasticHA()

void HAIntranuke2025::InelasticHA ( GHepRecord * ev,
GHepParticle * p,
INukeFateHA_t fate ) const
private

Definition at line 565 of file HAIntranuke2025.cxx.

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

References genie::GHepRecord::AddParticle(), genie::INukeHadroFates2025::AsString(), genie::Intranuke2025::fDoFermi, genie::Intranuke2025::fFermiFac, genie::Intranuke2025::fHadroData2025, genie::PDGLibrary::Find(), genie::Intranuke2025::fNuclmodel, genie::Intranuke2025::fRemnA, genie::Intranuke2025::fRemnP4, genie::Intranuke2025::fRemnZ, genie::PDGLibrary::Instance(), genie::RandomGen::Instance(), genie::pdg::IonPdgCode(), genie::kIHAFtCEx, genie::kIHAFtInelas, genie::kIHNFtCEx, genie::kIHNFtElas, genie::kIMdHA, genie::GHepParticle::KinE(), genie::kIStStableFinalState, genie::kPdgNeutron, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, LOG, genie::GHepParticle::Mass(), genie::GHepParticle::Name(), genie::GHepParticle::P4(), pDEBUG, genie::GHepParticle::Pdg(), pINFO, pNOTICE, genie::GHepRecord::Probe(), pWARN, genie::GHepParticle::Px(), genie::GHepParticle::Py(), genie::GHepParticle::Pz(), genie::RandomGen::RndFsi(), genie::Target::SetHitNucPdg(), genie::GHepParticle::SetMomentum(), genie::GHepParticle::SetPdgCode(), genie::exceptions::INukeException::SetReason(), genie::GHepParticle::SetStatus(), genie::GHepRecord::TargetNucleus(), and genie::utils::intranuke2025::TwoBodyCollision().

Referenced by SimulateHadronicFinalStateKinematics().

◆ LoadConfig()

void HAIntranuke2025::LoadConfig ( void )
privatevirtual

Implements genie::Intranuke2025.

Definition at line 1538 of file HAIntranuke2025.cxx.

1539{
1540
1541 // load hadronic cross sections
1542 fHadroData2025 = INukeHadroData2025::Instance();
1543
1544 // fermi momentum setup
1545 // this is specifically set in Intranuke2025::Configure(string)
1546 fNuclmodel = dynamic_cast<const NuclearModelI *>( this -> SubAlg("NuclearModel") ) ;
1547
1548 // other intranuke config params
1549 GetParam( "NUCL-R0", fR0 ); // fm
1550 GetParam( "NUCL-NR", fNR );
1551
1552 GetParam( "INUKE-NucRemovalE", fNucRmvE ); // GeV
1553 GetParam( "INUKE-HadStep", fHadStep ) ;
1554 GetParam( "INUKE-NucAbsFac", fNucAbsFac ) ;
1555 GetParam( "INUKE-NucCEXFac", fNucCEXFac ) ;
1556 GetParam( "INUKE-Energy_Pre_Eq", fEPreEq ) ;
1557 GetParam( "INUKE-FermiFac", fFermiFac ) ;
1558 GetParam( "INUKE-FermiMomentum", fFermiMomentum ) ;
1559
1560 GetParam( "INUKE-DoCompoundNucleus", fDoCompoundNucleus ) ;
1561 GetParam( "INUKE-DoFermi", fDoFermi ) ;
1562 GetParam( "INUKE-XsecNNCorr", fXsecNNCorr ) ;
1563 GetParamDef( "UseOset", fUseOset, false ) ;
1564 GetParamDef( "AltOset", fAltOset, false ) ;
1565
1566 GetParam( "HAINUKE-DelRPion", fDelRPion ) ;
1567 GetParam( "HAINUKE-DelRNucleon", fDelRNucleon ) ;
1568
1569 GetParamDef( "FSI-ChargedPion-MFPScale", fChPionMFPScale, 1.0 ) ;
1570 GetParamDef( "FSI-NeutralPion-MFPScale", fNeutralPionMFPScale, 1.0 ) ;
1571 GetParamDef( "FSI-Pion-FracCExScale", fPionFracCExScale, 1.0 ) ;
1572 GetParamDef( "FSI-ChargedPion-FracAbsScale", fChPionFracAbsScale, 1.0 ) ;
1573 GetParamDef( "FSI-NeutralPion-FracAbsScale", fNeutralPionFracAbsScale,1.0 ) ;
1574 GetParamDef( "FSI-Pion-FracInelScale", fPionFracInelScale, 1.0 ) ;
1575 GetParamDef( "FSI-Pion-FracPiProdScale", fPionFracPiProdScale, 1.0 ) ;
1576 GetParamDef( "FSI-Nucleon-MFPScale", fNucleonMFPScale, 1.0 ) ;
1577 GetParamDef( "FSI-Nucleon-FracCExScale", fNucleonFracCExScale , 1.0 ) ;
1578 GetParamDef( "FSI-Nucleon-FracInelScale", fNucleonFracInelScale, 1.0 ) ;
1579 GetParamDef( "FSI-Nucleon-FracAbsScale", fNucleonFracAbsScale, 1.0 ) ;
1580 GetParamDef( "FSI-Nucleon-FracPiProdScale", fNucleonFracPiProdScale, 1.0 ) ;
1581
1582 // report
1583 LOG("HAIntranuke2025", pINFO) << "Settings for INTRANUKE mode: " << INukeMode::AsString(kIMdHA);
1584 LOG("HAIntranuke2025", pINFO) << "R0 = " << fR0 << " fermi";
1585 LOG("HAIntranuke2025", pINFO) << "NR = " << fNR;
1586 LOG("HAIntranuke2025", pINFO) << "DelRPion = " << fDelRPion;
1587 LOG("HAIntranuke2025", pINFO) << "DelRNucleon = " << fDelRNucleon;
1588 LOG("HAIntranuke2025", pINFO) << "HadStep = " << fHadStep << " fermi";
1589 LOG("HAIntranuke2025", pINFO) << "EPreEq = " << fHadStep << " fermi";
1590 LOG("HAIntranuke2025", pINFO) << "NucAbsFac = " << fNucAbsFac;
1591 LOG("HAIntranuke2025", pINFO) << "NucCEXFac = " << fNucCEXFac;
1592 LOG("HAIntranuke2025", pINFO) << "FermiFac = " << fFermiFac;
1593 LOG("HAIntranuke2025", pINFO) << "FermiMomtm = " << fFermiMomentum;
1594 LOG("HAIntranuke2025", pINFO) << "DoFermi? = " << ((fDoFermi)?(true):(false));
1595 LOG("HAIntranuke2025", pINFO) << "DoCmpndNuc? = " << ((fDoCompoundNucleus)?(true):(false));
1596 LOG("HAIntranuke2025", pINFO) << "XsecNNCorr? = " << ((fXsecNNCorr)?(true):(false));
1597}
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::Algorithm::GetParamDef
bool GetParamDef(const RgKey &name, T &p, const T &def) const
genie::Algorithm::SubAlg
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition Algorithm.cxx:345
genie::INukeHadroData2025::Instance
static INukeHadroData2025 * Instance(void)
Definition INukeHadroData2025.cxx:184
genie::INukeMode::AsString
static string AsString(INukeMode_t mode)
Definition INukeMode.h:41
genie::Intranuke2025::fChPionMFPScale
double fChPionMFPScale
tweaking factors for tuning
Definition Intranuke2025.h:145
genie::Intranuke2025::fDelRPion
double fDelRPion
factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement
Definition Intranuke2025.h:130
genie::Intranuke2025::fNucAbsFac
double fNucAbsFac
absorption xsec correction factor (hN Mode)
Definition Intranuke2025.h:133
genie::Intranuke2025::fHadStep
double fHadStep
step size for intranuclear hadron transport
Definition Intranuke2025.h:132
genie::Intranuke2025::fNR
double fNR
param multiplying the nuclear radius, determining how far to track hadrons beyond the "nuclear bounda...
Definition Intranuke2025.h:128
genie::Intranuke2025::fUseOset
bool fUseOset
Oset model for low energy pion in hN.
Definition Intranuke2025.h:141
genie::Intranuke2025::fAltOset
bool fAltOset
NuWro's table-based implementation (not recommended)
Definition Intranuke2025.h:142
genie::Intranuke2025::fXsecNNCorr
bool fXsecNNCorr
use nuclear medium correction for NN cross section
Definition Intranuke2025.h:143
genie::Intranuke2025::fDoCompoundNucleus
bool fDoCompoundNucleus
whether or not to do compound nucleus considerations
Definition Intranuke2025.h:140
genie::Intranuke2025::fDelRNucleon
double fDelRNucleon
factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement
Definition Intranuke2025.h:131
genie::Intranuke2025::fNucCEXFac
double fNucCEXFac
charge exchange xsec correction factor (hN Mode)
Definition Intranuke2025.h:134
genie::Intranuke2025::fR0
double fR0
effective nuclear size param
Definition Intranuke2025.h:127
genie::Intranuke2025::fEPreEq
double fEPreEq
threshold for pre-equilibrium reaction
Definition Intranuke2025.h:135

References genie::INukeMode::AsString(), genie::Intranuke2025::fAltOset, genie::Intranuke2025::fChPionFracAbsScale, 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::fNeutralPionFracAbsScale, genie::Intranuke2025::fNeutralPionMFPScale, genie::Intranuke2025::fNR, genie::Intranuke2025::fNucAbsFac, genie::Intranuke2025::fNucCEXFac, genie::Intranuke2025::fNucleonFracAbsScale, genie::Intranuke2025::fNucleonFracCExScale, genie::Intranuke2025::fNucleonFracInelScale, genie::Intranuke2025::fNucleonFracPiProdScale, genie::Intranuke2025::fNucleonMFPScale, genie::Intranuke2025::fNuclmodel, genie::Intranuke2025::fNucRmvE, genie::Intranuke2025::fPionFracCExScale, genie::Intranuke2025::fPionFracInelScale, genie::Intranuke2025::fPionFracPiProdScale, genie::Intranuke2025::fR0, genie::Intranuke2025::fUseOset, genie::Intranuke2025::fXsecNNCorr, genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), genie::INukeHadroData2025::Instance(), genie::kIMdHA, LOG, pINFO, and genie::Algorithm::SubAlg().

◆ PiBounce()

double HAIntranuke2025::PiBounce ( void ) const
private

Definition at line 388 of file HAIntranuke2025.cxx.

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

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

Referenced by ElasHA().

◆ PnBounce()

double HAIntranuke2025::PnBounce ( void ) const
private

Definition at line 440 of file HAIntranuke2025.cxx.

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

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

Referenced by ElasHA().

◆ ProcessEventRecord()

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

Reimplemented from genie::Intranuke2025.

Definition at line 91 of file HAIntranuke2025.cxx.

92{
93 LOG("HAIntranuke2025", pNOTICE)
94 << "************ Running hA2025 MODE INTRANUKE ************";
95 GHepParticle * nuclearTarget = evrec -> TargetNucleus();
96 nuclA = nuclearTarget -> A();
97
98 Intranuke2025::ProcessEventRecord(evrec);
99
100 LOG("HAIntranuke2025", pINFO) << "Done with this event";
101}
genie::Intranuke2025::ProcessEventRecord
virtual void ProcessEventRecord(GHepRecord *event_rec) const
Definition Intranuke2025.cxx:82

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

◆ SimulateHadronicFinalState()

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

Implements genie::Intranuke2025.

Definition at line 103 of file HAIntranuke2025.cxx.

105{
106// Simulate a hadron interaction for the input particle p in HA mode
107//
108 // check inputs
109 if(!p || !ev) {
110 LOG("HAIntranuke2025", pERROR) << "** Null input!";
111 return;
112 }
113
114 // get particle code and check whether this particle can be handled
115 int pdgc = p->Pdg();
116 bool is_gamma = (pdgc==kPdgGamma);
117 bool is_pion = (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0);
118 bool is_kaon = (pdgc==kPdgKP); // sd: only K+. need many changes for K-
119 bool is_baryon = (pdgc==kPdgProton || pdgc==kPdgNeutron);
120 bool is_handled = (is_baryon || is_pion || is_kaon || is_gamma);
121 if(!is_handled) {
122 LOG("HAIntranuke2025", pERROR) << "** Can not handle particle: " << p->Name();
123 return;
124 }
125
126 // select a fate for the input particle
127 INukeFateHA_t fate = this->HadronFateHA(p);
128
129 // store the fate
130 ev->Particle(p->FirstMother())->SetRescatterCode((int)fate);
131
132 if(fate == kIHAFtUndefined) {
133 LOG("HAIntranuke2025", pERROR) << "** Couldn't select a fate";
134 p->SetStatus(kIStStableFinalState);
135 ev->AddParticle(*p);
136 return;
137 }
138 LOG("HAIntranuke2025", pNOTICE)
139 << "Selected "<< p->Name() << " fate: "<< INukeHadroFates2025::AsString(fate);
140
141 // try to generate kinematics - repeat till is done (should seldom need >2)
142 fNumIterations = 0;
143 this->SimulateHadronicFinalStateKinematics(ev,p);
144}
pERROR
#define pERROR
Definition Messenger.h:59
genie::GHepParticle::SetRescatterCode
void SetRescatterCode(int code)
Definition GHepParticle.h:129
genie::GHepRecord::Particle
virtual GHepParticle * Particle(int position) const
Definition GHepRecord.cxx:104
genie::HAIntranuke2025::HadronFateHA
INukeFateHA_t HadronFateHA(const GHepParticle *p) const
Definition HAIntranuke2025.cxx:203
genie::HAIntranuke2025::SimulateHadronicFinalStateKinematics
void SimulateHadronicFinalStateKinematics(GHepRecord *ev, GHepParticle *p) const
Definition HAIntranuke2025.cxx:146
genie::INukeFateHA_t
enum genie::EINukeFateHA_t INukeFateHA_t
genie::kPdgGamma
const int kPdgGamma
Definition PDGCodes.h:189

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

Referenced by SimulateHadronicFinalStateKinematics().

◆ SimulateHadronicFinalStateKinematics()

void HAIntranuke2025::SimulateHadronicFinalStateKinematics ( GHepRecord * ev,
GHepParticle * p ) const
private

Definition at line 146 of file HAIntranuke2025.cxx.

148{
149 // get stored fate
150 INukeFateHA_t fate = (INukeFateHA_t)
151 ev->Particle(p->FirstMother())->RescatterCode();
152
153 LOG("HAIntranuke2025", pINFO)
154 << "Generating kinematics for " << p->Name()
155 << " fate: "<< INukeHadroFates2025::AsString(fate);
156
157 // try to generate kinematics for the selected fate
158
159 try
160 {
161 fNumIterations++;
162 /* if (fate == kIHAFtElas)
163 {
164 this->ElasHA(ev,p,fate);
165 }
166 else */
167 if (fate == kIHAFtInelas || fate == kIHAFtCEx)
168 {
169 this->InelasticHA(ev,p,fate);
170 }
171 else if (fate == kIHAFtAbs || fate == kIHAFtPiProd)
172 {
173 this->Inelastic(ev,p,fate);
174 }
175 else if (fate == kIHAFtCmp) //(suarez edit, 17 July, 2017: cmp)
176 {
177 LOG("HAIntranuke2025", pWARN) << "Running PreEquilibrium for kIHAFtCmp";
178 LOG("HAIntranuke2025", pFATAL) << "The PreEquilibrium and Compound Nucleus code are experimental, should not be used";
179 utils::intranuke2025::PreEquilibrium(ev,p,fRemnA,fRemnZ,fRemnP4,fDoFermi,fFermiFac,fNuclmodel,fNucRmvE,kIMdHA); //should be kiMdHA or HN?
180 }
181 }
182 catch(exceptions::INukeException exception)
183 {
184 LOG("HAIntranuke2025", pNOTICE)
185 << exception;
186 if(fNumIterations <= 100) {
187 LOG("HAIntranuke2025", pNOTICE)
188 << "Failed attempt to generate kinematics for "
189 << p->Name() << " fate: " << INukeHadroFates2025::AsString(fate)
190 << " - After " << fNumIterations << " tries, still retrying...";
191 this->SimulateHadronicFinalStateKinematics(ev,p);
192 } else {
193 LOG("HAIntranuke2025", pNOTICE)
194 << "Failed attempt to generate kinematics for "
195 << p->Name() << " fate: " << INukeHadroFates2025::AsString(fate)
196 << " after " << fNumIterations-1
197 << " attempts. Trying a new fate...";
198 this->SimulateHadronicFinalState(ev,p);
199 }
200 }
201}
pFATAL
#define pFATAL
Definition Messenger.h:56
genie::GHepParticle::RescatterCode
int RescatterCode(void) const
Definition GHepParticle.h:65
genie::HAIntranuke2025::Inelastic
void Inelastic(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition HAIntranuke2025.cxx:751
genie::HAIntranuke2025::InelasticHA
void InelasticHA(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition HAIntranuke2025.cxx:565
genie::HAIntranuke2025::SimulateHadronicFinalState
void SimulateHadronicFinalState(GHepRecord *ev, GHepParticle *p) const
Definition HAIntranuke2025.cxx:103
genie::utils::intranuke2025::PreEquilibrium
void PreEquilibrium(GHepRecord *ev, GHepParticle *p, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, const NuclearModelI *Nuclmodel, double NucRmvE, EINukeMode mode=kIMdHN)
Definition INukeUtils2025.cxx:476

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

Referenced by SimulateHadronicFinalState(), and SimulateHadronicFinalStateKinematics().

◆ IntranukeTester

friend class IntranukeTester
friend

Definition at line 57 of file HAIntranuke2025.h.

References IntranukeTester.

Referenced by IntranukeTester.

Member Data Documentation

◆ fNumIterations

unsigned int genie::HAIntranuke2025::fNumIterations
mutableprivate

Definition at line 86 of file HAIntranuke2025.h.

Referenced by SimulateHadronicFinalState(), and SimulateHadronicFinalStateKinematics().

◆ nuclA

int genie::HAIntranuke2025::nuclA
mutableprivate

value of A for the target nucleus in hA mode

Definition at line 85 of file HAIntranuke2025.h.

Referenced by HadronFateHA(), and ProcessEventRecord().

The documentation for this class was generated from the following files:
Generated on for GENIEGenerator by doxygen 1.14.0