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

#include <HAIntranuke2015.h>

Inheritance diagram for genie::HAIntranuke2015:
genie::Intranuke2015 genie::EventRecordVisitorI genie::Algorithm

Public Member Functions

 HAIntranuke2015 ()
 
 HAIntranuke2015 (string config)
 
 ~HAIntranuke2015 ()
 
void ProcessEventRecord (GHepRecord *event_rec) const
 
virtual string GetINukeMode () const
 
- Public Member Functions inherited from genie::Intranuke2015
 Intranuke2015 ()
 
 Intranuke2015 (string name)
 
 Intranuke2015 (string name, string config)
 
 ~Intranuke2015 ()
 
void Configure (const Registry &config)
 Configure the algorithm. More...
 
void Configure (string param_set)
 Configure the algorithm. More...
 
- Public Member Functions inherited from genie::EventRecordVisitorI
virtual ~EventRecordVisitorI ()
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 Lookup configuration from the config pool. More...
 
virtual const Registry & GetConfig (void) const
 Get configuration registry. More...
 
Registry * GetOwnedConfig (void)
 Get a writeable version of an owned configuration Registry. More...
 
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

- Protected Member Functions inherited from genie::Intranuke2015
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)
 
- Protected Attributes inherited from genie::Intranuke2015
double fTrackingRadius
 tracking radius for the nucleus in the current event More...
 
TGenPhaseSpace fGenPhaseSpace
 a phase space generator More...
 
INukeHadroData2015 * fHadroData2015
 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 fFreeStep
 produced particle free stem, in fm 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...
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsConfig
 true if it owns its config. registry More...
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
Registry * fConfig
 config. (either owned or pointing to config pool) 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 HAIntranuke2015.h.

Constructor & Destructor Documentation

HAIntranuke2015::HAIntranuke2015 ( )

Definition at line 98 of file HAIntranuke2015.cxx.

98  :
99 Intranuke2015("genie::HAIntranuke2015")
100 {
101 
102 }
HAIntranuke2015::HAIntranuke2015 ( string  config)

Definition at line 104 of file HAIntranuke2015.cxx.

104  :
105 Intranuke2015("genie::HAIntranuke2015",config)
106 {
107 
108 }
generate_config.config
config
Definition: generate_config.py:3
HAIntranuke2015::~HAIntranuke2015 ( )

Definition at line 110 of file HAIntranuke2015.cxx.

111 {
112 
113 }

Member Function Documentation

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

Definition at line 455 of file HAIntranuke2015.cxx.

457 {
458  // scatters particle within nucleus, copy of hN code meant to run only once
459  // in hA mode
460 
461 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
462  LOG("HAIntranuke2015", pDEBUG)
463  << "ElasHA() is invoked for a : " << p->Name()
464  << " whose fate is : " << INukeHadroFates::AsString(fate);
465 #endif
466 
467  if(fate!=kIHAFtElas)
468  {
469  LOG("HAIntranuke2015", pWARN)
470  << "ElasHA() cannot handle fate: " << INukeHadroFates::AsString(fate);
471  return;
472  }
473 
474  // check remnants
475  if(fRemnA<0 || fRemnZ<0) // best to stop it here and not try again.
476  {
477  LOG("HAIntranuke2015", pWARN) << "Invalid Nucleus! : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
478  p->SetStatus(kIStStableFinalState);
479  ev->AddParticle(*p);
480  return;
481  }
482 
483  // vars for incoming particle, target, and scattered pdg codes
484  int pcode = p->Pdg();
485  double Mp = p->Mass();
486  double Mt = 0.;
487  if (ev->TargetNucleus()->A()==fRemnA)
488  { Mt = PDGLibrary::Instance()->Find(ev->TargetNucleus()->Pdg())->Mass(); }
489  else
490  {
491  Mt = fRemnP4.M();
492  }
493  TLorentzVector t4PpL = *p->P4();
494  TLorentzVector t4PtL = fRemnP4;
495  double C3CM = 0.0;
496 
497  // calculate scattering angle
498  if(pcode==kPdgNeutron||pcode==kPdgProton) C3CM = TMath::Cos(this->PnBounce());
499  else C3CM = TMath::Cos(this->PiBounce());
500 
501  // calculate final 4 momentum of probe
502  TLorentzVector t4P3L, t4P4L;
503 
504  if (!utils::intranuke2015::TwoBodyKinematics(Mp,Mt,t4PpL,t4PtL,t4P3L,t4P4L,C3CM,fRemnP4))
505  {
506  LOG("HAIntranuke2015", pNOTICE) << "ElasHA() failed";
507  exceptions::INukeException exception;
508  exception.SetReason("TwoBodyKinematics failed in ElasHA");
509  throw exception;
510  }
511 
512  // Update probe particle
513  p->SetMomentum(t4P3L);
514  p->SetStatus(kIStStableFinalState);
515 
516  // Update Remnant nucleus
517  fRemnP4 = t4P4L;
518  LOG("HAIntranuke2015",pINFO)
519  << "C3cm = " << C3CM;
520  LOG("HAIntranuke2015",pINFO)
521  << "|p3| = " << t4P3L.Vect().Mag() << ", E3 = " << t4P3L.E() << ",Mp = " << Mp;
522  LOG("HAIntranuke2015",pINFO)
523  << "|p4| = " << fRemnP4.Vect().Mag() << ", E4 = " << fRemnP4.E() << ",Mt = " << Mt;
524 
525  ev->AddParticle(*p);
526 
527 }
genie::exceptions::INukeException
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
Definition: INukeException.h:38
genie::utils::intranuke2015::TwoBodyKinematics
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
Definition: INukeUtils2015.cxx:886
genie::utils::res::Mass
double Mass(Resonance_t res)
resonance mass (GeV)
Definition: BaryonResUtils.cxx:435
genie::GHepParticle::SetMomentum
void SetMomentum(const TLorentzVector &p4)
Definition: GHepParticle.cxx:254
genie::GHepParticle::Mass
double Mass(void) const
Mass that corresponds to the PDG code.
Definition: GHepParticle.cxx:161
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:64
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:153
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:87
genie::HAIntranuke2015::PiBounce
double PiBounce(void) const
Definition: HAIntranuke2015.cxx:350
genie::GHepRecord::TargetNucleus
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:323
pINFO
#define pINFO
Definition: Messenger.h:53
pWARN
#define pWARN
Definition: Messenger.h:51
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:45
genie::GHepParticle::SetStatus
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:125
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:63
genie::Intranuke2015::fRemnZ
int fRemnZ
remnant nucleus Z
Definition: Intranuke2015.h:99
genie::Intranuke2015::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke2015.h:100
genie::GHepRecord::AddParticle
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:526
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:62
genie::GHepParticle::A
int A(void) const
Definition: GHepParticle.cxx:202
pNOTICE
#define pNOTICE
Definition: Messenger.h:52
genie::GHepParticle::P4
TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
genie::HAIntranuke2015::PnBounce
double PnBounce(void) const
Definition: HAIntranuke2015.cxx:402
genie::Intranuke2015::fRemnA
int fRemnA
remnant nucleus A
Definition: Intranuke2015.h:98
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:64
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
pDEBUG
#define pDEBUG
Definition: Messenger.h:54
virtual string genie::HAIntranuke2015::GetINukeMode ( ) const
inlinevirtual

Reimplemented from genie::Intranuke2015.

Definition at line 62 of file HAIntranuke2015.h.

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

Definition at line 226 of file HAIntranuke2015.cxx.

227 {
228 // Select a hadron fate in HA mode
229 //
230  RandomGen * rnd = RandomGen::Instance();
231 
232  // get pdgc code & kinetic energy in MeV
233  int pdgc = p->Pdg();
234  double ke = p->KinE() / units::MeV;
235 
236  //bool isPion = (pdgc == kPdgPiP or pdgc == kPdgPi0 or pdgc == kPdgPiM);
237  //if (isPion and fUseOset and ke < 350.0) return this->HadronFateOset();
238 
239  LOG("HAIntranuke2015", pINFO)
240  << "Selecting hA fate for " << p->Name() << " with KE = " << ke << " MeV";
241 
242  // try to generate a hadron fate
243  unsigned int iter = 0;
244  while(iter++ < kRjMaxIterations) {
245 
246  // handle pions
247  //
248  if (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0) {
249 
250  double frac_cex = fHadroData2015->FracADep(pdgc, kIHAFtCEx, ke, nuclA);
251  // double frac_elas = fHadroData2015->FracADep(pdgc, kIHAFtElas, ke, nuclA);
252  double frac_inel = fHadroData2015->FracADep(pdgc, kIHAFtInelas, ke, nuclA);
253  double frac_abs = fHadroData2015->FracADep(pdgc, kIHAFtAbs, ke, nuclA);
254  double frac_piprod = fHadroData2015->FracADep(pdgc, kIHAFtPiProd, ke, nuclA);
255  LOG("HAIntranuke2015", pDEBUG)
256  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCEx) << "} = " << frac_cex
257  // << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
258  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
259  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs
260  << "\n frac{" << INukeHadroFates::AsString(kIHAFtPiProd) << "} = " << frac_piprod;
261 
262  // compute total fraction (can be <1 if fates have been switched off)
263  double tf = frac_cex +
264  // frac_elas +
265  frac_inel +
266  frac_abs +
267  frac_piprod;
268 
269  double r = tf * rnd->RndFsi().Rndm();
270 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
271  LOG("HAIntranuke2015", pDEBUG) << "r = " << r << " (max = " << tf << ")";
272 #endif
273  double cf=0; // current fraction
274  if(r < (cf += frac_cex )) return kIHAFtCEx; // cex
275  // if(r < (cf += frac_elas )) return kIHAFtElas; // elas
276  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
277  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
278  if(r < (cf += frac_piprod )) return kIHAFtPiProd; // pi prod
279 
280  LOG("HAIntranuke2015", pWARN)
281  << "No selection after going through all fates! "
282  << "Total fraction = " << tf << " (r = " << r << ")";
283  }
284 
285  // handle nucleons
286  else if (pdgc==kPdgProton || pdgc==kPdgNeutron) {
287  double frac_cex = fHadroData2015->FracAIndep(pdgc, kIHAFtCEx, ke);
288  //double frac_elas = fHadroData2015->FracAIndep(pdgc, kIHAFtElas, ke);
289  double frac_inel = fHadroData2015->FracAIndep(pdgc, kIHAFtInelas, ke);
290  double frac_abs = fHadroData2015->FracAIndep(pdgc, kIHAFtAbs, ke);
291  double frac_pipro = fHadroData2015->FracAIndep(pdgc, kIHAFtPiProd, ke);
292  double frac_cmp = fHadroData2015->FracAIndep(pdgc, kIHAFtCmp , ke);
293 
294  LOG("HAIntranuke2015", pINFO)
295  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCEx) << "} = " << frac_cex
296  // << "\n frac{" << INukeHadroFates::AsString(kIHAFtElas) << "} = " << frac_elas
297  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
298  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs
299  << "\n frac{" << INukeHadroFates::AsString(kIHAFtPiProd) << "} = " << frac_pipro
300  << "\n frac{" << INukeHadroFates::AsString(kIHAFtCmp) << "} = " << frac_cmp; //suarez edit, cmp
301 
302  // compute total fraction (can be <1 if fates have been switched off)
303  double tf = frac_cex +
304  //frac_elas +
305  frac_inel +
306  frac_abs +
307  frac_pipro +
308  frac_cmp; //suarez edit, cmp
309 
310  double r = tf * rnd->RndFsi().Rndm();
311 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
312  LOG("HAIntranuke2015", pDEBUG) << "r = " << r << " (max = " << tf << ")";
313 #endif
314  double cf=0; // current fraction
315  if(r < (cf += frac_cex )) return kIHAFtCEx; // cex
316  //if(r < (cf += frac_elas )) return kIHAFtElas; // elas
317  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
318  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
319  if(r < (cf += frac_pipro )) return kIHAFtPiProd; // pi prod
320  if(r < (cf += frac_cmp )) return kIHAFtCmp; //suarez edit, cmp
321 
322  LOG("HAIntranuke2015", pWARN)
323  << "No selection after going through all fates! "
324  << "Total fraction = " << tf << " (r = " << r << ")";
325  }
326  // handle kaons
327  else if (pdgc==kPdgKP || pdgc==kPdgKM) {
328  double frac_inel = fHadroData2015->FracAIndep(pdgc, kIHAFtInelas, ke);
329  double frac_abs = fHadroData2015->FracAIndep(pdgc, kIHAFtAbs, ke);
330 
331  LOG("HAIntranuke2015", pDEBUG)
332  << "\n frac{" << INukeHadroFates::AsString(kIHAFtInelas) << "} = " << frac_inel
333  << "\n frac{" << INukeHadroFates::AsString(kIHAFtAbs) << "} = " << frac_abs;
334  // compute total fraction (can be <1 if fates have been switched off)
335  double tf = frac_inel +
336  frac_abs;
337  double r = tf * rnd->RndFsi().Rndm();
338 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
339  LOG("HAIntranuke2015", pDEBUG) << "r = " << r << " (max = " << tf << ")";
340 #endif
341  double cf=0; // current fraction
342  if(r < (cf += frac_inel )) return kIHAFtInelas; // inelas
343  if(r < (cf += frac_abs )) return kIHAFtAbs; // abs
344  }
345  }//iterations
346 
347  return kIHAFtUndefined;
348 }
genie::RandomGen::RndFsi
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:60
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
genie::units::MeV
static const double MeV
Definition: Units.h:130
genie::INukeHadroData2015::FracAIndep
double FracAIndep(int hpdgc, INukeFateHA_t fate, double ke) const
Definition: INukeHadroData2015.cxx:1383
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:30
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:64
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:153
genie::Intranuke2015::fHadroData2015
INukeHadroData2015 * fHadroData2015
a collection of h+N,h+A data & calculations
Definition: Intranuke2015.h:95
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:87
genie::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:147
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:146
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:132
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:134
pINFO
#define pINFO
Definition: Messenger.h:53
pWARN
#define pWARN
Definition: Messenger.h:51
genie::GHepParticle::KinE
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
Definition: GHepParticle.cxx:177
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::HAIntranuke2015::nuclA
int nuclA
value of A for the target nucleus in hA mode
Definition: HAIntranuke2015.h:80
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:133
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:62
genie::controls::kRjMaxIterations
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:64
genie::INukeHadroData2015::FracADep
double FracADep(int hpdgc, INukeFateHA_t fate, double ke, int targA) const
Definition: INukeHadroData2015.cxx:1326
plot_model.r
r
Definition: plot_model.py:48
pDEBUG
#define pDEBUG
Definition: Messenger.h:54
int HAIntranuke2015::HandleCompoundNucleus ( GHepRecord *  ev,
GHepParticle *  p,
int  mom 
) const
privatevirtual

Implements genie::Intranuke2015.

Definition at line 1428 of file HAIntranuke2015.cxx.

1429 {
1430 
1431  // only relevant for hN mode - not anymore.
1432  return false;
1433 
1434 }
void HAIntranuke2015::Inelastic ( GHepRecord *  ev,
GHepParticle *  p,
INukeFateHA_t  fate 
) const
private

Definition at line 715 of file HAIntranuke2015.cxx.

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

Definition at line 529 of file HAIntranuke2015.cxx.

531 {
532  // charge exch and inelastic - scatters particle within nucleus, hA version
533  // each are treated as quasielastic, particle scatters off single nucleon
534 
535 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
536  LOG("HAIntranuke2015", pDEBUG)
537  << "InelasticHA() is invoked for a : " << p->Name()
538  << " whose fate is : " << INukeHadroFates::AsString(fate);
539 #endif
540  if(ev->Probe() ) {
541  LOG("HAIntranuke2015", pINFO) << " probe KE = " << ev->Probe()->KinE();
542  }
543  if(fate!=kIHAFtCEx && fate!=kIHAFtInelas)
544  {
545  LOG("HAIntranuke2015", pWARN)
546  << "InelasticHA() cannot handle fate: " << INukeHadroFates::AsString(fate);
547  return;
548  }
549 
550  // Random number generator
551  RandomGen * rnd = RandomGen::Instance();
552 
553  // vars for incoming particle, target, and scattered pdg codes
554  int pcode = p->Pdg();
555  int tcode, scode, s2code;
556  double ppcnt = (double) fRemnZ / (double) fRemnA; // % of protons
557 
558  // Select a hadron fate in HN mode
559  INukeFateHN_t h_fate;
560  if (fate == kIHAFtCEx) h_fate = kIHNFtCEx;
561  else h_fate = kIHNFtElas;
562 
563  // Select a target randomly, weighted to #
564  // -- Unless, of course, the fate is CEx,
565  // -- in which case the target may be deterministic
566  // Also assign scattered particle code
567  if(fate==kIHAFtCEx)
568  {
569  if(pcode==kPdgPiP) {tcode = kPdgNeutron; scode = kPdgPi0; s2code = kPdgProton;}
570  else if(pcode==kPdgPiM) {tcode = kPdgProton; scode = kPdgPi0; s2code = kPdgNeutron;}
571  else if(pcode==kPdgPi0)
572  {
573  // for pi0
574  tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton) :(kPdgNeutron);
575  scode = (tcode == kPdgProton) ?(kPdgPiP) :(kPdgPiM);
576  s2code = (tcode == kPdgProton) ?(kPdgNeutron):(kPdgProton);
577  }
578  else if(pcode==kPdgProton) {tcode = kPdgNeutron; scode = kPdgNeutron; s2code = kPdgProton;}
579  else if(pcode==kPdgNeutron){tcode = kPdgProton; scode = kPdgProton; s2code = kPdgNeutron;}
580  else
581  { LOG("HAIntranuke2015", pWARN) << "InelasticHA() cannot handle fate: "
582  << INukeHadroFates::AsString(fate)
583  << " for particle " << p->Name();
584  return;
585  }
586  }
587  else
588  {
589  tcode = (rnd->RndFsi().Rndm()<=ppcnt)?(kPdgProton):(kPdgNeutron);
590  // if(pcode == kPdgKP || pcode == kPdgKM) tcode = kPdgProton;
591  scode = pcode;
592  s2code = tcode;
593  }
594 
595  // check remnants
596  if ( fRemnA < 1 ) //we've blown nucleus apart, no need to retry anything - exit
597  {
598  LOG("HAIntranuke2015",pNOTICE) << "InelasticHA() stops : not enough nucleons";
599  p->SetStatus(kIStStableFinalState);
600  ev->AddParticle(*p);
601  return;
602  }
603  else if ( fRemnZ + (((pcode==kPdgProton)||(pcode==kPdgPiP))?1:0) - (pcode==kPdgPiM?1:0)
604  < ((( scode==kPdgProton)||( scode==kPdgPiP)) ?1:0) - (scode ==kPdgPiM ?1:0)
605  + (((s2code==kPdgProton)||(s2code==kPdgPiP)) ?1:0) - (s2code==kPdgPiM ?1:0) )
606  {
607  LOG("HAIntranuke2015",pWARN) << "InelasticHA() failed : too few protons in nucleus";
608  p->SetStatus(kIStStableFinalState);
609  ev->AddParticle(*p);
610  return; // another extreme case, best strategy is to exit and go to next event
611  }
612 
613  GHepParticle t(*p);
614  t.SetPdgCode(tcode);
615 
616  // set up fermi target
617  Target target(ev->TargetNucleus()->Pdg());
618  double tM = t.Mass();
619 
620  // handle fermi momentum
621  if(fDoFermi)
622  {
623  target.SetHitNucPdg(tcode);
624  fNuclmodel->GenerateNucleon(target);
625  TVector3 tP3 = fFermiFac * fNuclmodel->Momentum3();
626  double tE = TMath::Sqrt(tP3.Mag2()+ tM*tM);
627  t.SetMomentum(TLorentzVector(tP3,tE));
628  }
629  else
630  {
631  t.SetMomentum(0,0,0,tM);
632  }
633 
634  GHepParticle * cl = new GHepParticle(*p); // clone particle, to run IntBounce at proper energy
635  // calculate energy and momentum using invariant mass
636  double pM = p->Mass();
637  double E_p = ((*p->P4() + *t.P4()).Mag2() - tM*tM - pM*pM)/(2.0*tM);
638  double P_p = TMath::Sqrt(E_p*E_p - pM*pM);
639  cl->SetMomentum(TLorentzVector(P_p,0,0,E_p));
640  // momentum doesn't have to be in right direction, only magnitude
641  double C3CM = fHadroData2015->IntBounce(cl,tcode,scode,h_fate);
642  delete cl;
643  if (C3CM<-1.) // hope this doesn't occur too often - unphysical but we just pass it on
644  {
645  LOG("HAIntranuke2015", pWARN) << "unphysical angle chosen in InelasicHA - put particle outside nucleus";
646  p->SetStatus(kIStStableFinalState);
647  ev->AddParticle(*p);
648  return;
649  }
650  double KE1L = p->KinE();
651  double KE2L = t.KinE();
652  LOG("HAIntranuke2015",pINFO)
653  << " KE1L = " << KE1L << " " << KE1L << " KE2L = " << KE2L;
654  GHepParticle cl1(*p);
655  GHepParticle cl2(t);
656  bool success = utils::intranuke2015::TwoBodyCollision(ev,pcode,tcode,scode,s2code,C3CM,
657  &cl1,&cl2,fRemnA,fRemnZ,fRemnP4,kIMdHA);
658  if(success)
659  {
660  double P3L = TMath::Sqrt(cl1.Px()*cl1.Px() + cl1.Py()*cl1.Py() + cl1.Pz()*cl1.Pz());
661  double P4L = TMath::Sqrt(cl2.Px()*cl2.Px() + cl2.Py()*cl2.Py() + cl2.Pz()*cl2.Pz());
662  double E3L = cl1.KinE();
663  double E4L = cl2.KinE();
664  LOG ("HAIntranuke2015",pINFO) << "Successful quasielastic scattering or charge exchange";
665  LOG("HAIntranuke",pINFO)
666  << "C3CM = " << C3CM << "\n P3L, E3L = "
667  << P3L << " " << E3L << " P4L, E4L = "<< P4L << " " << E4L ;
668  if(ev->Probe() ) { LOG("HAIntranuke",pINFO)
669  << "P4L = " << P4L << " ;E4L= " << E4L << "\n probe KE = " << ev->Probe()->KinE() << "\n";
670  LOG("HAIntranuke2015", pINFO) << "Nucleus : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
671  TParticlePDG * remn = 0;
672  double MassRem = 0.;
673  int ipdgc = pdg::IonPdgCode(fRemnA, fRemnZ);
674  remn = PDGLibrary::Instance()->Find(ipdgc);
675  if(!remn)
676  {
677  LOG("HAIntranuke2015", pINFO)
678  << "NO Particle with [A = " << fRemnA << ", Z = " << fRemnZ
679  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
680  }
681  else
682  {
683  MassRem = remn->Mass();
684  LOG("HAIntranuke2015", pINFO)
685  << "Particle with [A = " << fRemnA << ", Z = " << fRemnZ
686  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
687  }
688  double ERemn = fRemnP4.E();
689  double PRemn = TMath::Sqrt(fRemnP4.Px()*fRemnP4.Px() + fRemnP4.Py()*fRemnP4.Py() + fRemnP4.Pz()*fRemnP4.Pz());
690  double MRemn = TMath::Sqrt(ERemn*ERemn - PRemn*PRemn);
691  LOG("HAIntranuke2015",pINFO) << "PRemn = " << PRemn << " ;ERemn= " << ERemn;
692  LOG("HAIntranuke2015",pINFO) << "MRemn= " << MRemn << " ;true Mass= " << MassRem << " ; excitation energy= " << (MRemn-MassRem)*1000. << " MeV";
693  }
694  if (ev->Probe() && (E3L>ev->Probe()->KinE())) //assuming E3 is most important, definitely for pion. what about pp?
695  {
696  // LOG("HAIntranuke",pINFO)
697  // << "E3Lagain = " << E3L << " ;E4L= " << E4L << "\n probe KE = " << ev->Probe()->KinE() << "\n";
698  exceptions::INukeException exception;
699  exception.SetReason("TwoBodyCollison gives KE> probe KE in hA simulation");
700  throw exception;
701  }
702  ev->AddParticle(cl1);
703  ev->AddParticle(cl2);
704 
705  LOG("HAIntranuke2015", pDEBUG) << "Nucleus : (A,Z) = ("<<fRemnA<<','<<fRemnZ<<')';
706  } else
707  {
708  exceptions::INukeException exception;
709  exception.SetReason("TwoBodyCollison failed in hA simulation");
710  throw exception;
711  }
712 
713 }
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:60
genie::INukeHadroData2015::IntBounce
double IntBounce(const GHepParticle *p, int target, int s1, INukeFateHN_t fate)
Definition: INukeHadroData2015.cxx:1605
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
generate_CCQE_events.t
t
Definition: generate_CCQE_events.py:68
genie::GHepParticle::SetMomentum
void SetMomentum(const TLorentzVector &p4)
Definition: GHepParticle.cxx:254
genie::GHepParticle::Mass
double Mass(void) const
Mass that corresponds to the PDG code.
Definition: GHepParticle.cxx:161
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
generate_CCQE_events.target
target
Definition: generate_CCQE_events.py:44
genie::GHepRecord::Probe
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:314
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:64
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:153
genie::Intranuke2015::fHadroData2015
INukeHadroData2015 * fHadroData2015
a collection of h+N,h+A data & calculations
Definition: Intranuke2015.h:95
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:87
genie::Intranuke2015::fNuclmodel
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke2015.h:97
genie::Intranuke2015::fDoFermi
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke2015.h:116
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:132
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:134
genie::Intranuke2015::fFermiFac
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke2015.h:113
genie::GHepRecord::TargetNucleus
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:323
pINFO
#define pINFO
Definition: Messenger.h:53
genie::INukeFateHN_t
enum genie::EINukeFateHN_t INukeFateHN_t
pWARN
#define pWARN
Definition: Messenger.h:51
genie::GHepParticle::KinE
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
Definition: GHepParticle.cxx:177
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:45
genie::GHepParticle::SetStatus
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:125
genie::NuclearModelI::Momentum3
virtual TVector3 Momentum3(void) const
Definition: NuclearModelI.h:59
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:63
genie::pdg::IonPdgCode
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:82
genie::Intranuke2015::fRemnZ
int fRemnZ
remnant nucleus Z
Definition: Intranuke2015.h:99
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:133
genie::Intranuke2015::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke2015.h:100
genie::NuclearModelI::GenerateNucleon
virtual bool GenerateNucleon(const Target &) const =0
genie::GHepRecord::AddParticle
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:526
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:62
genie::utils::intranuke2015::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: INukeUtils2015.cxx:742
pNOTICE
#define pNOTICE
Definition: Messenger.h:52
genie::GHepParticle::P4
TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
genie::Intranuke2015::fRemnA
int fRemnA
remnant nucleus A
Definition: Intranuke2015.h:98
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:64
genie::GHepParticle
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
pDEBUG
#define pDEBUG
Definition: Messenger.h:54
void HAIntranuke2015::LoadConfig ( void  )
privatevirtual

Implements genie::Intranuke2015.

Definition at line 1436 of file HAIntranuke2015.cxx.

1437 {
1438 
1439  AlgConfigPool * confp = AlgConfigPool::Instance();
1440  const Registry * gc = confp->GlobalParameterList();
1441 
1442  // load hadronic cross sections
1443  fHadroData2015 = INukeHadroData2015::Instance();
1444 
1445  // fermi momentum setup
1446  fAlgf = AlgFactory::Instance();
1447  fNuclmodel = dynamic_cast<const NuclearModelI *>
1448  (fAlgf->GetAlgorithm("genie::FGMBodekRitchie","Default"));
1449 
1450  // other intranuke config params
1451  fR0 = fConfig->GetDoubleDef ("R0", gc->GetDouble("NUCL-R0")); // fm
1452  fNR = fConfig->GetDoubleDef ("NR", gc->GetDouble("NUCL-NR"));
1453  fNucRmvE = fConfig->GetDoubleDef ("NucRmvE", gc->GetDouble("INUKE-NucRemovalE")); // GeV
1454  fDelRPion = fConfig->GetDoubleDef ("DelRPion", gc->GetDouble("HAINUKE-DelRPion"));
1455  fDelRNucleon = fConfig->GetDoubleDef ("DelRNucleon", gc->GetDouble("HAINUKE-DelRNucleon"));
1456  fHadStep = fConfig->GetDoubleDef ("HadStep", gc->GetDouble("INUKE-HadStep")); // fm
1457  fEPreEq = fConfig->GetDoubleDef ("EPreEq", gc->GetDouble("INUKE-Energy_Pre_Eq"));
1458  fNucAbsFac = fConfig->GetDoubleDef ("NucAbsFac", gc->GetDouble("INUKE-NucAbsFac"));
1459  fNucCEXFac = fConfig->GetDoubleDef ("NucCEXFac", gc->GetDouble("INUKE-NucCEXFac"));
1460  fFermiFac = fConfig->GetDoubleDef ("FermiFac", gc->GetDouble("INUKE-FermiFac"));
1461  fFermiMomentum = fConfig->GetDoubleDef ("FermiMomentum",gc->GetDouble("INUKE-FermiMomentum"));
1462  fDoFermi = fConfig->GetBoolDef ("DoFermi", gc->GetBool("INUKE-DoFermi"));
1463  fFreeStep = fConfig->GetDoubleDef ("FreeStep", gc->GetDouble("INUKE-FreeStep"));
1464  fDoCompoundNucleus = fConfig->GetBoolDef ("DoCompoundNucleus", gc->GetBool("INUKE-DoCompoundNucleus"));
1465  fUseOset = fConfig->GetBoolDef ("UseOset", false);
1466  fAltOset = fConfig->GetBoolDef ("AltOset", false);
1467  fXsecNNCorr = fConfig->GetBoolDef ("XsecNNCorr", gc->GetBool("INUKE-XsecNNCorr"));
1468 
1469  // report
1470  LOG("HAIntranuke2015", pINFO) << "Settings for INTRANUKE mode: " << INukeMode::AsString(kIMdHA);
1471  LOG("HAIntranuke2015", pINFO) << "R0 = " << fR0 << " fermi";
1472  LOG("HAIntranuke2015", pINFO) << "NR = " << fNR;
1473  LOG("HAIntranuke2015", pINFO) << "DelRPion = " << fDelRPion;
1474  LOG("HAIntranuke2015", pINFO) << "DelRNucleon = " << fDelRNucleon;
1475  LOG("HAIntranuke2015", pINFO) << "HadStep = " << fHadStep << " fermi";
1476  LOG("HAIntranuke2015", pINFO) << "EPreEq = " << fHadStep << " fermi";
1477  LOG("HAIntranuke2015", pINFO) << "NucAbsFac = " << fNucAbsFac;
1478  LOG("HAIntranuke2015", pINFO) << "NucCEXFac = " << fNucCEXFac;
1479  LOG("HAIntranuke2015", pINFO) << "FermiFac = " << fFermiFac;
1480  LOG("HAIntranuke2015", pINFO) << "FreeStep = " << fFreeStep; // free step in fm
1481  LOG("HAIntranuke2015", pINFO) << "FermiMomtm = " << fFermiMomentum;
1482  LOG("HAIntranuke2015", pINFO) << "DoFermi? = " << ((fDoFermi)?(true):(false));
1483  LOG("HAIntranuke2015", pINFO) << "DoCmpndNuc? = " << ((fDoCompoundNucleus)?(true):(false));
1484  LOG("HAIntranuke2015", pINFO) << "XsecNNCorr? = " << ((fXsecNNCorr)?(true):(false));
1485 }
genie::INukeMode::AsString
static string AsString(INukeMode_t mode)
Definition: INukeMode.h:41
genie::Registry::GetDoubleDef
RgDbl GetDoubleDef(RgKey key, RgDbl def_opt, bool set_def=true)
Definition: Registry.cxx:549
genie::Intranuke2015::fUseOset
bool fUseOset
Oset model for low energy pion in hN.
Definition: Intranuke2015.h:119
genie::INukeHadroData2015::Instance
static INukeHadroData2015 * Instance(void)
Definition: INukeHadroData2015.cxx:157
genie::Intranuke2015::fR0
double fR0
effective nuclear size param
Definition: Intranuke2015.h:104
genie::Intranuke2015::fNucRmvE
double fNucRmvE
binding energy to subtract from cascade nucleons
Definition: Intranuke2015.h:106
genie::AlgConfigPool
A singleton class holding all configuration registries built while parsing all loaded XML configurati...
Definition: AlgConfigPool.h:41
genie::Intranuke2015::fXsecNNCorr
bool fXsecNNCorr
use nuclear medium correction for NN cross section
Definition: Intranuke2015.h:121
genie::Intranuke2015::fFreeStep
double fFreeStep
produced particle free stem, in fm
Definition: Intranuke2015.h:114
genie::NuclearModelI
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
genie::Intranuke2015::fEPreEq
double fEPreEq
threshold for pre-equilibrium reaction
Definition: Intranuke2015.h:112
genie::Registry::GetDouble
RgDbl GetDouble(RgKey key) const
Definition: Registry.cxx:488
genie::Intranuke2015::fFermiMomentum
double fFermiMomentum
whether or not particle collision is pauli blocked
Definition: Intranuke2015.h:115
genie::Intranuke2015::fNR
double fNR
param multiplying the nuclear radius, determining how far to track hadrons beyond the "nuclear bounda...
Definition: Intranuke2015.h:105
genie::Intranuke2015::fHadroData2015
INukeHadroData2015 * fHadroData2015
a collection of h+N,h+A data & calculations
Definition: Intranuke2015.h:95
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:87
genie::Intranuke2015::fDoCompoundNucleus
bool fDoCompoundNucleus
whether or not to do compound nucleus considerations
Definition: Intranuke2015.h:118
genie::Intranuke2015::fNuclmodel
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke2015.h:97
genie::Intranuke2015::fDoFermi
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke2015.h:116
genie::Intranuke2015::fAlgf
AlgFactory * fAlgf
algorithm factory instance
Definition: Intranuke2015.h:96
genie::AlgFactory::GetAlgorithm
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:86
genie::Intranuke2015::fFermiFac
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke2015.h:113
pINFO
#define pINFO
Definition: Messenger.h:53
genie::AlgConfigPool::GlobalParameterList
Registry * GlobalParameterList(void) const
Definition: AlgConfigPool.cxx:495
genie::AlgFactory::Instance
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
genie::Registry
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
genie::Registry::GetBool
RgBool GetBool(RgKey key) const
Definition: Registry.cxx:474
genie::Intranuke2015::fAltOset
bool fAltOset
NuWro&#39;s table-based implementation (not recommended)
Definition: Intranuke2015.h:120
genie::Algorithm::fConfig
Registry * fConfig
config. (either owned or pointing to config pool)
Definition: Algorithm.h:122
genie::Intranuke2015::fDelRNucleon
double fDelRNucleon
factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke2015.h:108
genie::Registry::GetBoolDef
RgBool GetBoolDef(RgKey key, RgBool def_opt, bool set_def=true)
Definition: Registry.cxx:539
genie::Intranuke2015::fHadStep
double fHadStep
step size for intranuclear hadron transport
Definition: Intranuke2015.h:109
genie::Intranuke2015::fNucAbsFac
double fNucAbsFac
absorption xsec correction factor (hN Mode)
Definition: Intranuke2015.h:110
genie::AlgConfigPool::Instance
static AlgConfigPool * Instance()
Definition: AlgConfigPool.cxx:98
genie::Intranuke2015::fDelRPion
double fDelRPion
factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke2015.h:107
genie::Intranuke2015::fNucCEXFac
double fNucCEXFac
charge exchange xsec correction factor (hN Mode)
Definition: Intranuke2015.h:111
double HAIntranuke2015::PiBounce ( void  ) const
private

Definition at line 350 of file HAIntranuke2015.cxx.

351 {
352 // [adapted from neugen3 intranuke_bounce.F]
353 // [is a fortran stub / difficult to understand - needs to be improved]
354 //
355 // Generates theta in radians for elastic pion-nucleus scattering/
356 // Lookup table is based on Fig 17 of Freedman, Miller and Henley, Nucl.Phys.
357 // A389, 457 (1982)
358 //
359  const int nprob = 25;
360  double dintor = 0.0174533;
361  double denom = 47979.453;
362  double rprob[nprob] = {
363  5000., 4200., 3000., 2600., 2100., 1800., 1200., 750., 500., 230., 120.,
364  35., 9., 3., 11., 18., 29., 27., 20., 14., 10., 6., 2., 0.14, 0.19 };
365 
366  double angles[nprob];
367  for(int i=0; i<nprob; i++) angles[i] = 2.5*i;
368 
369  RandomGen * rnd = RandomGen::Instance();
370  double r = rnd->RndFsi().Rndm();
371 
372  double xsum = 0.;
373  double theta = 0.;
374  double binl = 0.;
375  double binh = 0.;
376  int tj = 0;
377  for(int i=0; i<60; i++) {
378  theta = i+0.5;
379  for(int j=0; j < nprob-1; j++) {
380  binl = angles[j];
381  binh = angles[j+1];
382  tj=j;
383  if(binl<=theta && binh>=theta) break;
384  tj=0;
385  }//j
386  int itj = tj;
387  double tfract = (theta-binl)/2.5;
388  double delp = rprob[itj+1] - rprob[itj];
389  xsum += (rprob[itj] + tfract*delp)/denom;
390  if(xsum>r) break;
391  theta = 0.;
392  }//i
393 
394  theta *= dintor;
395 
396  LOG("HAIntranuke2015", pNOTICE)
397  << "Generated pi+A elastic scattering angle = " << theta << " radians";
398 
399  return theta;
400 }
genie::RandomGen::RndFsi
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:60
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
sensitivity_comps.i
int i
Definition: sensitivity_comps.py:15
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:87
pNOTICE
#define pNOTICE
Definition: Messenger.h:52
plot_model.r
r
Definition: plot_model.py:48
double HAIntranuke2015::PnBounce ( void  ) const
private

Definition at line 402 of file HAIntranuke2015.cxx.

403 {
404 // [adapted from neugen3 intranuke_pnbounce.F]
405 // [is a fortran stub / difficult to understand - needs to be improved]
406 //
407 // Generates theta in radians for elastic nucleon-nucleus scattering.
408 // Use 800 MeV p+O16 as template in same (highly simplified) spirit as pi+A
409 // from table in Adams et al., PRL 1979. Guess value at 0-2 deg based on Ni
410 // data.
411 //
412  const int nprob = 20;
413  double dintor = 0.0174533;
414  double denom = 11967.0;
415  double rprob[nprob] = {
416  2400., 2350., 2200., 2000., 1728., 1261., 713., 312., 106., 35.,
417  6., 5., 10., 12., 11., 9., 6., 1., 1., 1. };
418 
419  double angles[nprob];
420  for(int i=0; i<nprob; i++) angles[i] = 1.0*i;
421 
422  RandomGen * rnd = RandomGen::Instance();
423  double r = rnd->RndFsi().Rndm();
424 
425  double xsum = 0.;
426  double theta = 0.;
427  double binl = 0.;
428  double binh = 0.;
429  int tj = 0;
430  for(int i=0; i< nprob; i++) {
431  theta = i+0.5;
432  for(int j=0; j < nprob-1; j++) {
433  binl = angles[j];
434  binh = angles[j+1];
435  tj=j;
436  if(binl<=theta && binh>=theta) break;
437  tj=0;
438  }//j
439  int itj = tj;
440  double tfract = (theta-binl)/2.5;
441  double delp = rprob[itj+1] - rprob[itj];
442  xsum += (rprob[itj] + tfract*delp)/denom;
443  if(xsum>r) break;
444  theta = 0.;
445  }//i
446 
447  theta *= dintor;
448 
449  LOG("HAIntranuke2015", pNOTICE)
450  << "Generated N+A elastic scattering angle = " << theta << " radians";
451 
452  return theta;
453 }
genie::RandomGen::RndFsi
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:60
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
genie::RandomGen
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
sensitivity_comps.i
int i
Definition: sensitivity_comps.py:15
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:87
pNOTICE
#define pNOTICE
Definition: Messenger.h:52
plot_model.r
r
Definition: plot_model.py:48
void HAIntranuke2015::ProcessEventRecord ( GHepRecord *  event_rec) const
virtual

Reimplemented from genie::Intranuke2015.

Definition at line 115 of file HAIntranuke2015.cxx.

116 {
117  LOG("HAIntranuke2015", pNOTICE)
118  << "************ Running hA2015 MODE INTRANUKE ************";
119  GHepParticle * nuclearTarget = evrec -> TargetNucleus();
120  nuclA = nuclearTarget -> A();
121 
122  Intranuke2015::ProcessEventRecord(evrec);
123 
124  LOG("HAIntranuke2015", pINFO) << "Done with this event";
125 }
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:87
pINFO
#define pINFO
Definition: Messenger.h:53
genie::units::A
static const double A
Definition: Units.h:82
genie::Intranuke2015::ProcessEventRecord
virtual void ProcessEventRecord(GHepRecord *event_rec) const
Definition: Intranuke2015.cxx:116
genie::HAIntranuke2015::nuclA
int nuclA
value of A for the target nucleus in hA mode
Definition: HAIntranuke2015.h:80
pNOTICE
#define pNOTICE
Definition: Messenger.h:52
genie::GHepParticle
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
void HAIntranuke2015::SimulateHadronicFinalState ( GHepRecord *  ev,
GHepParticle *  p 
) const
privatevirtual

Implements genie::Intranuke2015.

Definition at line 127 of file HAIntranuke2015.cxx.

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

Definition at line 170 of file HAIntranuke2015.cxx.

172 {
173  // get stored fate
174  INukeFateHA_t fate = (INukeFateHA_t)
175  ev->Particle(p->FirstMother())->RescatterCode();
176 
177  LOG("HAIntranuke2015", pINFO)
178  << "Generating kinematics for " << p->Name()
179  << " fate: "<< INukeHadroFates::AsString(fate);
180 
181  // try to generate kinematics for the selected fate
182 
183  try
184  {
185  fNumIterations++;
186  if (fate == kIHAFtElas)
187  {
188  this->ElasHA(ev,p,fate);
189  }
190  else
191  if (fate == kIHAFtInelas || fate == kIHAFtCEx)
192  {
193  this->InelasticHA(ev,p,fate);
194  }
195  else if (fate == kIHAFtAbs || fate == kIHAFtPiProd)
196  {
197  this->Inelastic(ev,p,fate);
198  }
199  else if (fate == kIHAFtCmp) //(suarez edit, 17 July, 2017: cmp)
200  {
201  LOG("HAIntranuke2015", pWARN) << "Running PreEquilibrium for kIHAFtCmp";
202  utils::intranuke2015::PreEquilibrium(ev,p,fRemnA,fRemnZ,fRemnP4,fDoFermi,fFermiFac,fNuclmodel,fNucRmvE,kIMdHA); //should be kiMdHA or HN?
203  }
204  }
205  catch(exceptions::INukeException exception)
206  {
207  LOG("HAIntranuke2015", pNOTICE)
208  << exception;
209  if(fNumIterations <= 100) {
210  LOG("HAIntranuke2015", pNOTICE)
211  << "Failed attempt to generate kinematics for "
212  << p->Name() << " fate: " << INukeHadroFates::AsString(fate)
213  << " - After " << fNumIterations << " tries, still retrying...";
214  this->SimulateHadronicFinalStateKinematics(ev,p);
215  } else {
216  LOG("HAIntranuke2015", pNOTICE)
217  << "Failed attempt to generate kinematics for "
218  << p->Name() << " fate: " << INukeHadroFates::AsString(fate)
219  << " after " << fNumIterations-1
220  << " attempts. Trying a new fate...";
221  this->SimulateHadronicFinalState(ev,p);
222  }
223  }
224 }
genie::INukeFateHA_t
enum genie::EINukeFateHA_t INukeFateHA_t
genie::HAIntranuke2015::InelasticHA
void InelasticHA(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition: HAIntranuke2015.cxx:529
genie::exceptions::INukeException
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
Definition: INukeException.h:38
genie::GHepRecord::Particle
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:148
genie::HAIntranuke2015::SimulateHadronicFinalStateKinematics
void SimulateHadronicFinalStateKinematics(GHepRecord *ev, GHepParticle *p) const
Definition: HAIntranuke2015.cxx:170
genie::HAIntranuke2015::ElasHA
void ElasHA(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition: HAIntranuke2015.cxx:455
genie::Intranuke2015::fNucRmvE
double fNucRmvE
binding energy to subtract from cascade nucleons
Definition: Intranuke2015.h:106
genie::HAIntranuke2015::SimulateHadronicFinalState
void SimulateHadronicFinalState(GHepRecord *ev, GHepParticle *p) const
Definition: HAIntranuke2015.cxx:127
genie::HAIntranuke2015::Inelastic
void Inelastic(GHepRecord *ev, GHepParticle *p, INukeFateHA_t fate) const
Definition: HAIntranuke2015.cxx:715
genie::GHepParticle::FirstMother
int FirstMother(void) const
Definition: GHepParticle.h:67
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:153
genie::utils::intranuke2015::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: INukeUtils2015.cxx:457
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:87
genie::Intranuke2015::fNuclmodel
const NuclearModelI * fNuclmodel
nuclear model used to generate fermi momentum
Definition: Intranuke2015.h:97
genie::Intranuke2015::fDoFermi
bool fDoFermi
whether or not to do fermi mom.
Definition: Intranuke2015.h:116
genie::Intranuke2015::fFermiFac
double fFermiFac
testing parameter to modify fermi momentum
Definition: Intranuke2015.h:113
genie::HAIntranuke2015::fNumIterations
unsigned int fNumIterations
Definition: HAIntranuke2015.h:81
pINFO
#define pINFO
Definition: Messenger.h:53
pWARN
#define pWARN
Definition: Messenger.h:51
genie::INukeHadroFates::AsString
static string AsString(INukeFateHN_t fate)
Definition: INukeHadroFates.h:68
genie::Intranuke2015::fRemnZ
int fRemnZ
remnant nucleus Z
Definition: Intranuke2015.h:99
genie::Intranuke2015::fRemnP4
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke2015.h:100
pNOTICE
#define pNOTICE
Definition: Messenger.h:52
genie::Intranuke2015::fRemnA
int fRemnA
remnant nucleus A
Definition: Intranuke2015.h:98
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Friends And Related Function Documentation

friend class IntranukeTester
friend

Definition at line 53 of file HAIntranuke2015.h.

Member Data Documentation

unsigned int genie::HAIntranuke2015::fNumIterations
mutableprivate

Definition at line 81 of file HAIntranuke2015.h.

int genie::HAIntranuke2015::nuclA
mutableprivate

value of A for the target nucleus in hA mode

Definition at line 80 of file HAIntranuke2015.h.

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