Intranuke2018.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: Steve Dytman <dytman+@pitt.edu>, Pittsburgh Univ.
8  Aaron Meyer <asm58@pitt.edu>, Pittsburgh Univ.
9  Alex Bell, Pittsburgh Univ.
10  Hugh Gallagher <gallag@minos.phy.tufts.edu>, Tufts Univ.
11  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>, Rutherford Lab.
12  September 20, 2005
13 
14  For the class documentation see the corresponding header file.
15 
16  Important revisions after version 2.0.0 :
17  @ Nov 30, 2007 - SD
18  Changed the hadron tracking algorithm to take into account the radial
19  nuclear density dependence. Using the somewhat empirical approach of
20  increasing the nuclear radius by a const (tunable) number times the tracked
21  particle's de Broglie wavelength as this helps getting the hadron+nucleus
22  cross sections right.
23  @ Mar 08, 2008 - CA
24  Fixed code retrieving the remnant nucleus which stopped working as soon as
25  simulation of nuclear de-excitation started pushing photons in the target
26  nucleus daughter list.
27  @ Jun 20, 2008 - CA
28  Fix a mem leak: The (clone of the) GHepParticle being re-scattered was not
29  deleted after it was added at the GHEP event record.
30  @ Jan 28, 2009 - CA
31  The nuclear remnant is now marked as a kIStFinalStateNuclearRemnant, not
32  as a kIStStableFinalState.
33  @ Sep 15, 2009 - CA
34  IsFake() and IsNucleus() are no longer available in GHepParticle.
35  Use pdg::IsPseudoParticle() and pdg::IsIon().
36  @ Sep 15, 2009 - CA
37  Store the rescattering code (hadron fate) in the GHEP record so as to
38  facilitate event reweighting.
39  @ Jul 15, 2010 - AM
40  Split Intranuke class into two separate classes, one for each interaction mode.
41  Intranuke.cxx now only contains methods common to both classes and associated
42  with the stepping of the hadrons through the nucleus and with configuration.
43  @ Nov 20, 2011 - CA
44  Tweaked the way TransportHadrons() looks-up the nuclear environment so that
45  it works for the nucleon decay mode as well.
46  @ Dec 08, 2011 - CA
47  Some minor structural changes. The new GEvGenMode_t is determined at the
48  start of the event processing and is used throughout. fInTestMode flag and
49  special INTRANUKE configs not needed. ProcessEventRecord() was added by
50  factoring out code from HNIntranuke and HAIntranuke. Some comments added.
51  @ Dec 23, 2014 - TG, SD
52  New 2014 class for latest Intranuke model
53  @ Apr 26, 2018 - SD
54  Change year 2015 to 2018
55 
56 */
57 //____________________________________________________________________________
58 
59 #include <cstdlib>
60 #include <sstream>
61 
62 #include <TMath.h>
63 
87 
88 using std::ostringstream;
89 
90 using namespace genie;
91 using namespace genie::utils;
92 using namespace genie::constants;
93 using namespace genie::controls;
94 
95 //___________________________________________________________________________
98 {
99 
100 }
101 //___________________________________________________________________________
104 {
105 
106 }
107 //___________________________________________________________________________
109 EventRecordVisitorI(name, config)
110 {
111 
112 }
113 //___________________________________________________________________________
115 {
116 
117 }
118 //___________________________________________________________________________
120 {
121  // Do not continue if there is no nuclear target
122  GHepParticle * nucltgt = evrec->TargetNucleus();
123  if (!nucltgt) {
124  LOG("HNIntranuke2018", pINFO) << "No nuclear target found - INTRANUKE exits";
125  return;
126  }
127 
128  // Decide tracking radius for the current nucleus (few * R0 * A^1/3)
129  this->SetTrackingRadius(nucltgt);
130 
131  // Understand what the event generation mode is (hadron/photon-nucleus,
132  // lepton-nucleus, nucleon decay) from the input event.
133  // The determined mode has an effect on INTRANUKE behaviour (how to lookup
134  // the residual nucleus, whether to set an intranuclear vtx etc) but it
135  // does not affect the INTRANUKE physics.
136  fGMode = evrec->EventGenerationMode();
137 
138  // For lepton-nucleus scattering and for nucleon decay intranuclear vtx
139  // position (in the target nucleus coord system) is set elsewhere.
140  // This method only takes effect in hadron/photon-nucleus interactions.
141  // In this special mode, an interaction vertex is set at the periphery
142  // of the target nucleus.
143  if(fGMode == kGMdHadronNucleus ||
145  {
146  this->GenerateVertex(evrec);
147  }
148 
149  // Now transport all hadrons outside the tracking radius.
150  // Stepping part is common for both HA and HN.
151  // Once it has been estabished that an interaction takes place then
152  // HA and HN specific code takes over in order to simulate the final state.
153  this->TransportHadrons(evrec);
154 }
155 //___________________________________________________________________________
157 {
158 // Sets a vertex in the nucleus periphery
159 // Called onlt in hadron/photon-nucleus interactions.
160 
161  GHepParticle * nucltgt = evrec->TargetNucleus();
162  assert(nucltgt);
163 
165  TVector3 vtx(999999.,999999.,999999.);
166 
167  // *** For h+A events (test mode):
168  // Assume a hadron beam with uniform intensity across an area,
169  // so we need to choose events uniformly within that area.
170  double x=999999., y=999999., epsilon = 0.001;
171  double R2 = TMath::Power(fTrackingRadius,2.);
172  double rp2 = TMath::Power(x,2.) + TMath::Power(y,2.);
173  while(rp2 > R2-epsilon) {
174  x = (fTrackingRadius-epsilon) * rnd->RndFsi().Rndm();
175  y = -fTrackingRadius + 2*fTrackingRadius * rnd->RndFsi().Rndm();
176  y -= ((y>0) ? epsilon : -epsilon);
177  rp2 = TMath::Power(x,2.) + TMath::Power(y,2.);
178  }
179  vtx.SetXYZ(x,y, -1.*TMath::Sqrt(TMath::Max(0.,R2-rp2)) + epsilon);
180 
181  // get the actual unit vector along the incoming hadron direction
182  TVector3 direction = evrec->Probe()->P4()->Vect().Unit();
183 
184  // rotate the vtx position
185  vtx.RotateUz(direction);
186 
187  LOG("Intranuke2018", pNOTICE)
188  << "Generated vtx @ R = " << vtx.Mag() << " fm / "
189  << print::Vec3AsString(&vtx);
190 
191  TObjArrayIter piter(evrec);
192  GHepParticle * p = 0;
193  while( (p = (GHepParticle *) piter.Next()) )
194  {
195  if(pdg::IsPseudoParticle(p->Pdg())) continue;
196  if(pdg::IsIon (p->Pdg())) continue;
197 
198  p->SetPosition(vtx.x(), vtx.y(), vtx.z(), 0.);
199  }
200 }
201 //___________________________________________________________________________
203 {
204  assert(p && pdg::IsIon(p->Pdg()));
205  double A = p->A();
206  fTrackingRadius = fR0 * TMath::Power(A, 1./3.);
207 
208  // multiply that by some input factor so that hadrons are tracked
209  // beyond the nuclear 'boundary' since the nuclear density distribution
210  // is not zero there
211  fTrackingRadius *= fNR;
212 
213  LOG("Intranuke2018", pNOTICE)
214  << "Setting tracking radius to R = " << fTrackingRadius;
215 }
216 //___________________________________________________________________________
218 {
219 // checks whether the particle should be rescattered
220 
221  assert(p);
222 
223  if(fGMode == kGMdHadronNucleus ||
225  // hadron/photon-nucleus scattering propagate the incoming particle
226  return (
228  && !pdg::IsIon(p->Pdg()));
229  }
230  else {
231  // attempt to rescatter anything marked as 'hadron in the nucleus'
232  return (p->Status() == kIStHadronInTheNucleus);
233  }
234 }
235 //___________________________________________________________________________
237 {
238 // checks whether a particle that needs to be rescattered, can in fact be
239 // rescattered by this cascade MC
240 
241  assert(p);
242  return ( p->Pdg() == kPdgPiP ||
243  p->Pdg() == kPdgPiM ||
244  p->Pdg() == kPdgPi0 ||
245  p->Pdg() == kPdgProton ||
246  p->Pdg() == kPdgNeutron ||
247  // p->Pdg() == kPdgGamma ||
248  p->Pdg() == kPdgKP //||
249  // p->Pdg() == kPdgKM
250  );
251 }
252 //___________________________________________________________________________
254 {
255 // check whether the input particle is still within the nucleus
256 //
257  return (p->X4()->Vect().Mag() < fTrackingRadius + fHadStep);
258 }
259 //___________________________________________________________________________
261 {
262 // transport all hadrons outside the nucleus
263 
264  int inucl = -1;
265  fRemnA = -1;
266  fRemnZ = -1;
267 
268  // Get 'nuclear environment' at the beginning of hadron transport
269  // and keep track of the remnant nucleus A,Z
270 
271  if(fGMode == kGMdHadronNucleus ||
273  {
274  inucl = evrec->TargetNucleusPosition();
275  }
276  else if(fGMode == kGMdLeptonNucleus ||
279  inucl = evrec->RemnantNucleusPosition();
280  }
281 
282  LOG("Intranuke2018", pNOTICE)
283  << "Propagating hadrons within nucleus found in position = " << inucl;
284  GHepParticle * nucl = evrec->Particle(inucl);
285  if(!nucl) {
286  LOG("Intranuke2018", pERROR)
287  << "No nucleus found in position = " << inucl;
288  LOG("Intranuke2018", pERROR)
289  << *evrec;
290  return;
291  }
292 
293  fRemnA = nucl->A();
294  fRemnZ = nucl->Z();
295 
296  LOG("Intranuke2018", pNOTICE)
297  << "Nucleus (A,Z) = (" << fRemnA << ", " << fRemnZ << ")";
298 
299  const TLorentzVector & p4nucl = *(nucl->P4());
300  fRemnP4 = p4nucl;
301 
302  // Loop over GHEP and run intranuclear rescattering on handled particles
303  TObjArrayIter piter(evrec);
304  GHepParticle * p = 0;
305  int icurr = -1;
306 
307  while( (p = (GHepParticle *) piter.Next()) )
308  {
309  icurr++;
310 
311  // Check whether the particle needs rescattering, otherwise skip it
312  if( ! this->NeedsRescattering(p) ) continue;
313 
314  LOG("Intranuke2018", pNOTICE)
315  << " >> Stepping a " << p->Name()
316  << " with kinetic E = " << p->KinE() << " GeV";
317 
318  // Rescatter a clone, not the original particle
319  GHepParticle * sp = new GHepParticle(*p);
320 
321  // Set clone's mom to be the hadron that was cloned
322  sp->SetFirstMother(icurr);
323 
324  // Check whether the particle can be rescattered
325  if(!this->CanRescatter(sp)) {
326 
327  // if I can't rescatter it, I will just take it out of the nucleus
328  LOG("Intranuke2018", pNOTICE)
329  << "... Current version can't rescatter a " << sp->Name();
330  sp->SetFirstMother(icurr);
332  evrec->AddParticle(*sp);
333  delete sp;
334  continue; // <-- skip to next GHEP entry
335  }
336 
337  // Start stepping particle out of the nucleus
338  bool has_interacted = false;
339  while ( this-> IsInNucleus(sp) )
340  {
341  // advance the hadron by a step
343 
344  // check whether it interacts
345  double d = this->GenerateStep(evrec,sp);
346  has_interacted = (d<fHadStep);
347  if(has_interacted) break;
348  }//stepping
349 
350  //updating the position of the original particle with the position of the clone
351  evrec->Particle(sp->FirstMother())->SetPosition(*(sp->X4()));
352 
353  if(has_interacted && fRemnA>0) {
354  // the particle interacts - simulate the hadronic interaction
355  LOG("Intranuke2018", pNOTICE)
356  << "Particle has interacted at location: "
357  << sp->X4()->Vect().Mag() << " / nucl rad= " << fTrackingRadius;
358  this->SimulateHadronicFinalState(evrec,sp);
359  } else if(has_interacted && fRemnA<=0) {
360  // nothing left to interact with!
361  LOG("Intranuke2018", pNOTICE)
362  << "*** Nothing left to interact with, escaping.";
364  evrec->AddParticle(*sp);
365  evrec->Particle(sp->FirstMother())->SetRescatterCode(1);
366  } else {
367  // the exits the nucleus without interacting - Done with it!
368  LOG("Intranuke2018", pNOTICE)
369  << "*** Hadron escaped the nucleus! Done with it.";
371  evrec->AddParticle(*sp);
372  evrec->Particle(sp->FirstMother())->SetRescatterCode(1);
373  }
374  delete sp;
375 
376  // Current snapshot
377  //LOG("Intranuke2018", pINFO) << "Current event record snapshot: " << *evrec;
378 
379  }// GHEP entries
380 
381  // Add remnant nucleus - that 'hadronic blob' has all the remaining hadronic
382  // 4p not put explicitly into the simulated particles
383  TLorentzVector v4(0.,0.,0.,0.);
384  GHepParticle remnant_nucleus(
386  evrec->AddParticle(remnant_nucleus);
387  // Mark the initial remnant nucleus as an intermediate state
388  // Don't do that in the hadron/photon-nucleus scatterig mode since the initial
389  // remnant nucleus and the target nucleus coincide.
390  if(fGMode != kGMdHadronNucleus &&
392  evrec->Particle(inucl)->SetStatus(kIStIntermediateState);
393  }
394 }
395 //___________________________________________________________________________
396 double Intranuke2018::GenerateStep(GHepRecord* /*evrec*/, GHepParticle* p) const //Added ev to get tgt argument//
397 {
398 // Generate a step (in fermis) for particle p in the input event.
399 // Computes the mean free path L and generate an 'interaction' distance d
400 // from an exp(-d/L) distribution
401 
402  int pdgc = p->Pdg();
403 
404  double scale = 1.;
405  if (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0) {
406  scale = fPionMFPScale;
407  }
408  else if (pdgc==kPdgProton || pdgc==kPdgNeutron) {
409  scale = fNucleonMFPScale;
410  }
411 
413 
414  string fINukeMode = this->GetINukeMode();
415  string fINukeModeGen = this->GetGenINukeMode();
416 
417  double L = utils::intranuke2018::MeanFreePath(p->Pdg(), *p->X4(), *p->P4(), fRemnA,
419 
420  LOG("Intranuke2018", pDEBUG) << "mode= " << fINukeModeGen;
421  if(fINukeModeGen == "hA") L *= scale;
422 
423  double d = -1.*L * TMath::Log(rnd->RndFsi().Rndm());
424 
425  /* LOG("Intranuke2018", pDEBUG)
426  << "mode= " << fINukeMode << "; Mean free path = " << L << " fm / "
427  << "Generated path length = " << d << " fm";
428  */
429  return d;
430 }
431 //___________________________________________________________________________
433 {
434  Algorithm::Configure(config);
435  this->LoadConfig();
436 }
437 //___________________________________________________________________________
438 void Intranuke2018::Configure(string param_set)
439 {
440  Algorithm::Configure(param_set);
441  this->LoadConfig();
442 }
443 //___________________________________________________________________________
const XML_Char * name
Definition: expat.h:151
int Z(void) const
void SetFirstMother(int m)
Definition: GHepParticle.h:133
Basic constants.
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:148
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
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
double fDelRPion
factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement ...
virtual string GetINukeMode() const
Definition: Intranuke2018.h:72
const int kPdgHadronicBlob
Definition: PDGCodes.h:188
virtual int RemnantNucleusPosition(void) const
Definition: GHepRecord.cxx:433
double fNR
param multiplying the nuclear radius, determining how far to track hadrons beyond the "nuclear bounda...
void StepParticle(GHepParticle *p, double step, double nuclear_radius=-1.)
Step particle.
const char * p
Definition: xmltok.h:285
virtual void Configure(const Registry &config)
int fRemnA
remnant nucleus A
Definition: Intranuke2018.h:99
bool fXsecNNCorr
use nuclear medium correction for NN cross section
Definition: config.py:1
bool IsInNucleus(const GHepParticle *p) const
bool CanRescatter(const GHepParticle *p) const
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
GEvGenMode_t fGMode
event generation mode (lepton+A, hadron+A, ...)
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
int Pdg(void) const
Definition: GHepParticle.h:64
int FirstMother(void) const
Definition: GHepParticle.h:67
string Name(void) const
Name that corresponds to the PDG code.
Double_t scale
Definition: plot.C:25
void SetPosition(const TLorentzVector &v4)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
static constexpr double L
double fPionMFPScale
tweaking factors for tuning
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
GEvGenMode_t EventGenerationMode(void) const
Definition: GHepRecord.cxx:253
const int kPdgKP
Definition: PDGCodes.h:149
bool NeedsRescattering(const GHepParticle *p) const
Float_t d
Definition: plot.C:236
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
void SetTrackingRadius(const GHepParticle *p) const
#define pINFO
Definition: Messenger.h:63
virtual void LoadConfig(void)=0
double MeanFreePath(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double nRpi=0.5, double nRnuc=1.0, const bool useOset=false, const bool altOset=false, const bool xsecNNCorr=false, string INukeMode="XX2018")
Mean free path (pions, nucleons)
Misc GENIE control constants.
void GenerateVertex(GHepRecord *ev) const
double fR0
effective nuclear size param
double GenerateStep(GHepRecord *ev, GHepParticle *p) const
virtual void ProcessEventRecord(GHepRecord *event_rec) const
virtual string GetGenINukeMode() const
Definition: Intranuke2018.h:73
static const double A
Definition: Units.h:82
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
bool IsIon(int pdgc)
Definition: PDGUtils.cxx:40
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:127
bool fAltOset
NuWro&#39;s table-based implementation (not recommended)
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
bool IsPseudoParticle(int pdgc)
Definition: PDGUtils.cxx:25
double epsilon
double fDelRNucleon
factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement ...
double fHadStep
step size for intranuclear hadron transport
assert(nhit_max >=nhit_nbins)
const int kPdgPiM
Definition: PDGCodes.h:136
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
void TransportHadrons(GHepRecord *ev) const
const int kPdgProton
Definition: PDGCodes.h:65
int A(void) const
#define pNOTICE
Definition: Messenger.h:62
double fTrackingRadius
tracking radius for the nucleus in the current event
Definition: Intranuke2018.h:94
const int kPdgNeutron
Definition: PDGCodes.h:67
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
virtual void SimulateHadronicFinalState(GHepRecord *ev, GHepParticle *p) const =0
bool fUseOset
Oset model for low energy pion in hN.
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
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:87
virtual int TargetNucleusPosition(void) const
Definition: GHepRecord.cxx:406
Root of GENIE utility namespaces.
int fRemnZ
remnant nucleus Z
#define pDEBUG
Definition: Messenger.h:64
TLorentzVector fRemnP4
P4 of remnant system.