Intranuke.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 
52 */
53 //____________________________________________________________________________
54 
55 #include <cstdlib>
56 #include <sstream>
57 
58 #include <TMath.h>
59 
82 
83 using std::ostringstream;
84 
85 using namespace genie;
86 using namespace genie::utils;
87 using namespace genie::constants;
88 using namespace genie::controls;
89 
90 //___________________________________________________________________________
93 {
94 
95 }
96 //___________________________________________________________________________
99 {
100 
101 }
102 //___________________________________________________________________________
104 EventRecordVisitorI(name, config)
105 {
106 
107 }
108 //___________________________________________________________________________
110 {
111 
112 }
113 //___________________________________________________________________________
115 {
116  // Do not continue if there is no nuclear target
117  GHepParticle * nucltgt = evrec->TargetNucleus();
118  if (!nucltgt) {
119  LOG("HNIntranuke", pINFO) << "No nuclear target found - INTRANUKE exits";
120  return;
121  }
122 
123  // Decide tracking radius for the current nucleus (few * R0 * A^1/3)
124  this->SetTrackingRadius(nucltgt);
125 
126  // Understand what the event generation mode is (hadron/photon-nucleus,
127  // lepton-nucleus, nucleon decay) from the input event.
128  // The determined mode has an effect on INTRANUKE behaviour (how to lookup
129  // the residual nucleus, whether to set an intranuclear vtx etc) but it
130  // does not affect the INTRANUKE physics.
131  fGMode = evrec->EventGenerationMode();
132 
133  // For lepton-nucleus scattering and for nucleon decay intranuclear vtx
134  // position (in the target nucleus coord system) is set elsewhere.
135  // This method only takes effect in hadron/photon-nucleus interactions.
136  // In this special mode, an interaction vertex is set at the periphery
137  // of the target nucleus.
138  if(fGMode == kGMdHadronNucleus ||
140  {
141  this->GenerateVertex(evrec);
142  }
143 
144  // Now transport all hadrons outside the tracking radius.
145  // Stepping part is common for both HA and HN.
146  // Once it has been estabished that an interaction takes place then
147  // HA and HN specific code takes over in order to simulate the final state.
148  this->TransportHadrons(evrec);
149 }
150 //___________________________________________________________________________
152 {
153 // Sets a vertex in the nucleus periphery
154 // Called onlt in hadron/photon-nucleus interactions.
155 
156  GHepParticle * nucltgt = evrec->TargetNucleus();
157  assert(nucltgt);
158 
160  TVector3 vtx(999999.,999999.,999999.);
161 
162  // *** For h+A events (test mode):
163  // Assume a hadron beam with uniform intensity across an area,
164  // so we need to choose events uniformly within that area.
165  double x=999999., y=999999., epsilon = 0.001;
166  double R2 = TMath::Power(fTrackingRadius,2.);
167  double rp2 = TMath::Power(x,2.) + TMath::Power(y,2.);
168  while(rp2 > R2-epsilon) {
169  x = (fTrackingRadius-epsilon) * rnd->RndFsi().Rndm();
170  y = -fTrackingRadius + 2*fTrackingRadius * rnd->RndFsi().Rndm();
171  y -= ((y>0) ? epsilon : -epsilon);
172  rp2 = TMath::Power(x,2.) + TMath::Power(y,2.);
173  }
174  vtx.SetXYZ(x,y, -1.*TMath::Sqrt(TMath::Max(0.,R2-rp2)) + epsilon);
175 
176  // get the actual unit vector along the incoming hadron direction
177  TVector3 direction = evrec->Probe()->P4()->Vect().Unit();
178 
179  // rotate the vtx position
180  vtx.RotateUz(direction);
181 
182  LOG("Intranuke", pNOTICE)
183  << "Generated vtx @ R = " << vtx.Mag() << " fm / "
184  << print::Vec3AsString(&vtx);
185 
186  TObjArrayIter piter(evrec);
187  GHepParticle * p = 0;
188  while( (p = (GHepParticle *) piter.Next()) )
189  {
190  if(pdg::IsPseudoParticle(p->Pdg())) continue;
191  if(pdg::IsIon (p->Pdg())) continue;
192 
193  p->SetPosition(vtx.x(), vtx.y(), vtx.z(), 0.);
194  }
195 }
196 //___________________________________________________________________________
198 {
199  assert(p && pdg::IsIon(p->Pdg()));
200  double A = p->A();
201  fTrackingRadius = fR0 * TMath::Power(A, 1./3.);
202 
203  // multiply that by some input factor so that hadrons are tracked
204  // beyond the nuclear 'boundary' since the nuclear density distribution
205  // is not zero there
206  fTrackingRadius *= fNR;
207 
208  LOG("Intranuke", pNOTICE)
209  << "Setting tracking radius to R = " << fTrackingRadius;
210 }
211 //___________________________________________________________________________
213 {
214 // checks whether the particle should be rescattered
215 
216  assert(p);
217 
218  if(fGMode == kGMdHadronNucleus ||
220  // hadron/photon-nucleus scattering propagate the incoming particle
221  return (
223  && !pdg::IsIon(p->Pdg()));
224  }
225  else {
226  // attempt to rescatter anything marked as 'hadron in the nucleus'
227  return (p->Status() == kIStHadronInTheNucleus);
228  }
229 }
230 //___________________________________________________________________________
232 {
233 // checks whether a particle that needs to be rescattered, can in fact be
234 // rescattered by this cascade MC
235 
236  assert(p);
237  return ( p->Pdg() == kPdgPiP ||
238  p->Pdg() == kPdgPiM ||
239  p->Pdg() == kPdgPi0 ||
240  p->Pdg() == kPdgProton ||
241  p->Pdg() == kPdgNeutron ||
242  // p->Pdg() == kPdgGamma ||
243  p->Pdg() == kPdgKP //||
244  // p->Pdg() == kPdgKM
245  );
246 }
247 //___________________________________________________________________________
249 {
250 // check whether the input particle is still within the nucleus
251 //
252  return (p->X4()->Vect().Mag() < fTrackingRadius + fHadStep);
253 }
254 //___________________________________________________________________________
256 {
257 // transport all hadrons outside the nucleus
258 
259  int inucl = -1;
260  fRemnA = -1;
261  fRemnZ = -1;
262 
263  // Get 'nuclear environment' at the beginning of hadron transport
264  // and keep track of the remnant nucleus A,Z
265 
266  if(fGMode == kGMdHadronNucleus ||
268  {
269  inucl = evrec->TargetNucleusPosition();
270  }
271  else
272  if(fGMode == kGMdLeptonNucleus ||
276  {
277  inucl = evrec->RemnantNucleusPosition();
278  }
279 
280  LOG("Intranuke", pNOTICE)
281  << "Propagating hadrons within nucleus found in position = " << inucl;
282  GHepParticle * nucl = evrec->Particle(inucl);
283  if(!nucl) {
284  LOG("Intranuke", pERROR)
285  << "No nucleus found in position = " << inucl;
286  LOG("Intranuke", pERROR)
287  << *evrec;
288  return;
289  }
290 
291  fRemnA = nucl->A();
292  fRemnZ = nucl->Z();
293 
294  LOG("Intranuke", pNOTICE)
295  << "Nucleus (A,Z) = (" << fRemnA << ", " << fRemnZ << ")";
296 
297  const TLorentzVector & p4nucl = *(nucl->P4());
298  fRemnP4 = p4nucl;
299 
300  // Loop over GHEP and run intranuclear rescattering on handled particles
301  TObjArrayIter piter(evrec);
302  GHepParticle * p = 0;
303  int icurr = -1;
304 
305  while( (p = (GHepParticle *) piter.Next()) )
306  {
307  icurr++;
308 
309  // Check whether the particle needs rescattering, otherwise skip it
310  if( ! this->NeedsRescattering(p) ) continue;
311 
312  if(this->HandleCompoundNucleus(evrec,p,icurr)) continue;
313 
314  LOG("Intranuke", 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("Intranuke", 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  if(has_interacted && fRemnA>0) {
351  // the particle interacts - simulate the hadronic interaction
352  LOG("Intranuke", pNOTICE)
353  << "Particle has interacted at location: "
354  << sp->X4()->Vect().Mag() << " / nucl rad= " << fTrackingRadius;
355  this->SimulateHadronicFinalState(evrec,sp);
356  } else if(has_interacted && fRemnA<=0) {
357  // nothing left to interact with!
358  LOG("Intranuke", pNOTICE)
359  << "*** Nothing left to interact with, escaping.";
361  evrec->AddParticle(*sp);
362  evrec->Particle(sp->FirstMother())->SetRescatterCode(1);
363  } else {
364  // the exits the nucleus without interacting - Done with it!
365  LOG("Intranuke", pNOTICE)
366  << "*** Hadron escaped the nucleus! Done with it.";
368  evrec->AddParticle(*sp);
369  evrec->Particle(sp->FirstMother())->SetRescatterCode(1);
370  }
371  delete sp;
372 
373  // Current snapshot
374  //LOG("Intranuke", pINFO) << "Current event record snapshot: " << *evrec;
375 
376  }// GHEP entries
377 
378  // Add remnant nucleus - that 'hadronic blob' has all the remaining hadronic
379  // 4p not put explicitly into the simulated particles
380  TLorentzVector v4(0.,0.,0.,0.);
381  GHepParticle remnant_nucleus(
383  evrec->AddParticle(remnant_nucleus);
384  // Mark the initial remnant nucleus as an intermediate state
385  // Don't do that in the hadron/photon-nucleus scatterig mode since the initial
386  // remnant nucleus and the target nucleus coincide.
387  if(fGMode != kGMdHadronNucleus &&
389  evrec->Particle(inucl)->SetStatus(kIStIntermediateState);
390  }
391 }
392 //___________________________________________________________________________
394 {
395 // Generate a step (in fermis) for particle p in the input event.
396 // Computes the mean free path L and generate an 'interaction' distance d
397 // from an exp(-d/L) distribution
398 
400 
401  int pdgc = p->Pdg();
402 
403  double scale = 1.;
404  if (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0) {
405  scale = fPionMFPScale;
406  }
407  else if (pdgc==kPdgProton || pdgc==kPdgNeutron) {
408  scale = fNucleonMFPScale;
409  }
410 
411  double L = utils::intranuke::MeanFreePath(pdgc, *p->X4(), *p->P4(), fRemnA,
413  L *= scale;
414  double d = -1.*L * TMath::Log(rnd->RndFsi().Rndm());
415 
416  LOG("Intranuke", pDEBUG)
417  << "Mean free path = " << L << " fm / "
418  << "Generated path length = " << d << " fm";
419  return d;
420 }
421 //___________________________________________________________________________
423 {
424  Algorithm::Configure(config);
425  this->LoadConfig();
426 }
427 //___________________________________________________________________________
428 void Intranuke::Configure(string param_set)
429 {
430  Algorithm::Configure(param_set);
431  this->LoadConfig();
432 }
433 //___________________________________________________________________________
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
TLorentzVector fRemnP4
P4 of remnant system.
Definition: Intranuke.h:98
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...
const int kPdgHadronicBlob
Definition: PDGCodes.h:188
virtual int RemnantNucleusPosition(void) const
Definition: GHepRecord.cxx:433
const char * p
Definition: xmltok.h:285
double fHadStep
step size for intranuclear hadron transport
Definition: Intranuke.h:107
double fTrackingRadius
tracking radius for the nucleus in the current event
Definition: Intranuke.h:91
virtual void ProcessEventRecord(GHepRecord *event_rec) const
Definition: Intranuke.cxx:114
Definition: config.py:1
void SetTrackingRadius(const GHepParticle *p) const
Definition: Intranuke.cxx:197
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
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
int Pdg(void) const
Definition: GHepParticle.h:64
double fDelRNucleon
factor by which Nucleon Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke.h:106
int FirstMother(void) const
Definition: GHepParticle.h:67
string Name(void) const
Name that corresponds to the PDG code.
double fNucleonMFPScale
Definition: Intranuke.h:123
GEvGenMode_t fGMode
event generation mode (lepton+A, hadron+A, ...)
Definition: Intranuke.h:99
Double_t scale
Definition: plot.C:25
double fDelRPion
factor by which Pion Compton wavelength gets multiplied to become nuclear size enhancement ...
Definition: Intranuke.h:105
void SetPosition(const TLorentzVector &v4)
virtual void SimulateHadronicFinalState(GHepRecord *ev, GHepParticle *p) const =0
#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
bool NeedsRescattering(const GHepParticle *p) const
Definition: Intranuke.cxx:212
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
Float_t d
Definition: plot.C:236
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
virtual void LoadConfig(void)=0
double fR0
effective nuclear size param
Definition: Intranuke.h:102
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
#define pINFO
Definition: Messenger.h:63
int fRemnA
remnant nucleus A
Definition: Intranuke.h:96
bool CanRescatter(const GHepParticle *p) const
Definition: Intranuke.cxx:231
virtual bool HandleCompoundNucleus(GHepRecord *ev, GHepParticle *p, int mom) const =0
Misc GENIE control constants.
void GenerateVertex(GHepRecord *ev) const
Definition: Intranuke.cxx:151
double MeanFreePath(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double nRpi=0.5, double nRnuc=1.0)
Mean free path (pions, nucleons)
Definition: INukeUtils.cxx:70
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
bool IsIon(int pdgc)
Definition: PDGUtils.cxx:40
void TransportHadrons(GHepRecord *ev) const
Definition: Intranuke.cxx:255
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
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
bool IsPseudoParticle(int pdgc)
Definition: PDGUtils.cxx:25
double epsilon
int fRemnZ
remnant nucleus Z
Definition: Intranuke.h:97
assert(nhit_max >=nhit_nbins)
const int kPdgPiM
Definition: PDGCodes.h:136
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
const int kPdgProton
Definition: PDGCodes.h:65
double GenerateStep(GHepRecord *ev, GHepParticle *p) const
Definition: Intranuke.cxx:393
int A(void) const
#define pNOTICE
Definition: Messenger.h:62
double fPionMFPScale
tweaking factors for tuning
Definition: Intranuke.h:117
const int kPdgNeutron
Definition: PDGCodes.h:67
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
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
bool IsInNucleus(const GHepParticle *p) const
Definition: Intranuke.cxx:248
virtual int TargetNucleusPosition(void) const
Definition: GHepRecord.cxx:406
Root of GENIE utility namespaces.
void Configure(const Registry &config)
Definition: Intranuke.cxx:422
double fNR
param multiplying the nuclear radius, determining how far to track hadrons beyond the "nuclear bounda...
Definition: Intranuke.h:103
#define pDEBUG
Definition: Messenger.h:64