Public Member Functions | Private Member Functions | Private Attributes | Friends | List of all members
genie::HAIntranuke2018 Class Reference

#include <HAIntranuke2018.h>

Inheritance diagram for genie::HAIntranuke2018:
genie::Intranuke2018 genie::EventRecordVisitorI genie::Algorithm

Public Member Functions

 HAIntranuke2018 ()
 
 HAIntranuke2018 (string config)
 
 ~HAIntranuke2018 ()
 
void ProcessEventRecord (GHepRecord *event_rec) const
 
virtual string GetINukeMode () const
 
virtual string GetGenINukeMode () const
 
- Public Member Functions inherited from genie::Intranuke2018
 Intranuke2018 ()
 
 Intranuke2018 (string name)
 
 Intranuke2018 (string name, string config)
 
 ~Intranuke2018 ()
 
virtual void Configure (const Registry &config)
 
virtual void Configure (string param_set)
 
void SetRemnA (int A)
 
void SetRemnZ (int Z)
 
double GetRemnA () const
 
double GetRemnZ () const
 
double GetR0 () const
 
double GetNR () const
 
double GetDelRPion () const
 
double GetDelRNucleon () const
 
double GetNucRmvE () const
 
double GetHadStep () const
 
bool GetUseOset () const
 
bool GetAltOset () const
 
bool GetXsecNNCorr () const
 
- Public Member Functions inherited from genie::EventRecordVisitorI
virtual ~EventRecordVisitorI ()
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 
virtual const Registry & GetConfig (void) const
 
Registry * GetOwnedConfig (void)
 
virtual const AlgId & Id (void) const
 Get algorithm ID. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
virtual void SetId (string name, string config)
 
const Algorithm * SubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Private Member Functions

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

Private Attributes

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

Friends

class IntranukeTester
 

Additional Inherited Members

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

Detailed Description

Definition at line 51 of file HAIntranuke2018.h.

Constructor & Destructor Documentation

HAIntranuke2018::HAIntranuke2018 ( )

Definition at line 96 of file HAIntranuke2018.cxx.

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

Definition at line 102 of file HAIntranuke2018.cxx.

102  :
103 Intranuke2018("genie::HAIntranuke2018",config)
104 {
105 
106 }
config
static Config * config
Definition: config.cpp:1054
HAIntranuke2018::~HAIntranuke2018 ( )

Definition at line 108 of file HAIntranuke2018.cxx.

109 {
110 
111 }

Member Function Documentation

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

Definition at line 481 of file HAIntranuke2018.cxx.

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

Reimplemented from genie::Intranuke2018.

Definition at line 63 of file HAIntranuke2018.h.

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

Reimplemented from genie::Intranuke2018.

Definition at line 62 of file HAIntranuke2018.h.

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

Definition at line 224 of file HAIntranuke2018.cxx.

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

Implements genie::Intranuke2018.

Definition at line 1512 of file HAIntranuke2018.cxx.

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

Definition at line 739 of file HAIntranuke2018.cxx.

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

Definition at line 553 of file HAIntranuke2018.cxx.

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

Implements genie::Intranuke2018.

Definition at line 1520 of file HAIntranuke2018.cxx.

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

Definition at line 376 of file HAIntranuke2018.cxx.

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

Definition at line 428 of file HAIntranuke2018.cxx.

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

Reimplemented from genie::Intranuke2018.

Definition at line 113 of file HAIntranuke2018.cxx.

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

Implements genie::Intranuke2018.

Definition at line 125 of file HAIntranuke2018.cxx.

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

Definition at line 168 of file HAIntranuke2018.cxx.

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

Friends And Related Function Documentation

friend class IntranukeTester
friend

Definition at line 53 of file HAIntranuke2018.h.

Member Data Documentation

unsigned int genie::HAIntranuke2018::fNumIterations
mutableprivate

Definition at line 82 of file HAIntranuke2018.h.

int genie::HAIntranuke2018::nuclA
mutableprivate

value of A for the target nucleus in hA mode

Definition at line 81 of file HAIntranuke2018.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11