INukeDeltaPropg.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Oct 01, 2009 - CA
14  Was first added in v2.5.1
15 
16 */
17 //____________________________________________________________________________
18 
19 #include <TLorentzVector.h>
20 #include <TVector3.h>
21 
32 
33 using namespace genie;
34 
35 //___________________________________________________________________________
37 EventRecordVisitorI("genie::INukeDeltaPropg")
38 {
39 
40 }
41 //___________________________________________________________________________
43 EventRecordVisitorI("genie::INukeDeltaPropg", config)
44 {
45 
46 }
47 //___________________________________________________________________________
49 {
50 
51 }
52 //___________________________________________________________________________
54 {
55  // Check that we have an interaction with a nuclear target. If not skip...
56  GHepParticle * nucltgt = event->TargetNucleus();
57  if (!nucltgt) {
58  LOG("INukeDelta", pINFO)
59  << "No nuclear target. Skipping....";
60  return;
61  }
62 
63  // mass number, nuclear radius, step size
64  int A = nucltgt->A();
65  double step_sz = fHadStep;
66  double nucl_radius = utils::nuclear::Radius(A,fR0);
67  nucl_radius *= fNR; // track the particle further out
68 
70 
71  // Loop over GHEP rescatter handled particles
72  TObjArrayIter piter(event);
73  GHepParticle * p = 0;
74  int icurr = -1;
75 
76  while( (p = (GHepParticle *) piter.Next()) )
77  {
78  icurr++;
79 
80  // handle?
81  int pdgc = p->Pdg();
82  bool delta = (pdgc == kPdgP33m1232_DeltaPP);
83  if (!delta) return;
84  GHepStatus_t ist = p->Status();
85  bool in_nucleus = (ist == kIStHadronInTheNucleus);
86  if (!in_nucleus) return;
87 
88  LOG("INukeDelta", pNOTICE)
89  << " >> Stepping a " << p->Name()
90  << " with kinetic E = " << p->KinE() << " GeV";
91 
92  // Rescatter a clone, not the original particle
93  GHepParticle * sp = new GHepParticle(*p);
94 
95  // Set clone's mom to be the hadron that was cloned
96  sp->SetFirstMother(icurr);
97 
98  // Start stepping particle out of the nucleus
99  bool has_interacted = false;
100  bool has_decayed = false;
101  while (1) {
102 
103  const TLorentzVector & p4 = *(sp->P4());
104  const TLorentzVector & x4 = *(sp->X4());
105 
106  bool is_in = (x4.Vect().Mag() < nucl_radius + step_sz);
107  if (!is_in) break;
108 
109  // step
110  utils::intranuke::StepParticle(sp, step_sz, nucl_radius);
111 
112  // check whether it decayed at this step
113  double Ldec = 0.; // calculate
114  double ddec = -1. * Ldec * TMath::Log(rnd->RndFsi().Rndm());
115  has_decayed = (ddec < step_sz);
116  if(has_decayed) break;
117 
118  // check whether it interacted at this step
119  double Lint = utils::intranuke::MeanFreePath_Delta(pdgc,x4,p4,A);
120  double dint = -1. * Lint * TMath::Log(rnd->RndFsi().Rndm());
121  has_interacted = (dint < step_sz);
122  if(has_interacted) break;
123 
124  }//stepping
125 
126  if(has_decayed) {
127  // the particle decays
128 
129  }
130  else
131  if(has_interacted) {
132  // the particle interacts - simulate the hadronic interaction
133  LOG("INukeDelta", pINFO)
134  << "Particle has interacted at location: "
135  << sp->X4()->Vect().Mag() << " / nucl radius = " << nucl_radius;
136 
137  //
138  // *** Temporary / Used in test mode only ***
139  // For now, just change the Delta++ to a proton to prevent the
140  // Delta++ decay later on.
141  // Causes energy and charge non-conservation
142  //
143  sp->SetPdgCode(kPdgProton);
145  event->AddParticle(*sp);
146  }
147  else {
148  // the particle escapes the nucleus
149  LOG("Intranuke", pNOTICE)
150  << "*** Hadron escaped the nucleus! Done with it.";
152  event->AddParticle(*sp);
153  }//decay? interacts? escapes?
154 
155  }//particle loop
156 }
157 //___________________________________________________________________________
159 {
160  Algorithm::Configure(config);
161  this->LoadConfig();
162 }
163 //___________________________________________________________________________
165 {
166  Algorithm::Configure(config);
167  this->LoadConfig();
168 }
169 //___________________________________________________________________________
171 {
172 
173  GetParam( "NUCL-R0", fR0 ) ; //fm
174  GetParam( "NUCL-NR", fNR ) ;
175  GetParam( "INUKE-HadStep", fHadStep ) ; //fm
176 
177 }
178 //___________________________________________________________________________
179 
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:91
void SetFirstMother(int m)
Definition: GHepParticle.h:133
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:60
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
double delta
Definition: runWimpSim.h:98
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
enum genie::EGHepStatus GHepStatus_t
const char * p
Definition: xmltok.h:285
double MeanFreePath_Delta(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A)
Mean free path (Delta++ test)
Definition: INukeUtils.cxx:170
Definition: config.py:1
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
GHepStatus_t Status(void) const
Definition: GHepParticle.h:65
int Pdg(void) const
Definition: GHepParticle.h:64
double Radius(int A, double Ro=constants::kNucRo)
string Name(void) const
Name that corresponds to the PDG code.
double fR0
effective nuclear size param
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void Configure(const Registry &config)
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
#define pINFO
Definition: Messenger.h:63
double fNR
param multiplying the nuclear radius, determining how far to track hadrons beyond the "nuclear bounda...
static const double A
Definition: Units.h:82
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
void StepParticle(GHepParticle *p, double step, double nuclear_radius=-1.)
Step particle.
Definition: INukeUtils.cxx:341
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:127
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const int kPdgProton
Definition: PDGCodes.h:65
void ProcessEventRecord(GHepRecord *event_rec) const
int A(void) const
#define pNOTICE
Definition: Messenger.h:62
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
double fHadStep
step size for intranuclear hadron transport
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
void SetPdgCode(int c)