GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DISHadronicSystemGenerator.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <RVersion.h>
12 
29 
30 using namespace genie;
31 using namespace genie::controls;
32 using namespace genie::constants;
33 using namespace genie::utils;
34 
35 //___________________________________________________________________________
37 HadronicSystemGenerator("genie::DISHadronicSystemGenerator")
38 {
39 
40 }
41 //___________________________________________________________________________
43 HadronicSystemGenerator("genie::DISHadronicSystemGenerator", config)
44 {
45 
46 }
47 //___________________________________________________________________________
49 {
50 
51 }
52 //___________________________________________________________________________
54 {
55 // This method generates the final state hadronic system
56 
57  //-- Add an entry for the DIS Pre-Fragm. Hadronic State
58  this->AddFinalHadronicSyst(evrec);
59 
60  // set W in the Event Summary
61  //-- Compute the hadronic system invariant mass
62  TLorentzVector p4Had = this->Hadronic4pLAB(evrec);
63 
64  Interaction * interaction = evrec->Summary();
65  interaction->KinePtr()->SetW( p4Had.M() );
66 
67  //-- Add the fragmentation products
69 
70 
71  //-- Simulate the formation zone if not taken directly from the
72  // hadronization model
73  this->SimulateFormationZone(evrec);
74 }
75 //___________________________________________________________________________
77  GHepRecord * evrec) const
78 {
79  LOG("DISHadronicVtx", pDEBUG)
80  << "Simulating formation zone for the DIS hadronic system";
81 
82  GHepParticle * nucltgt = evrec->TargetNucleus();
83  if (!nucltgt) {
84  LOG("DISHadronicVtx", pDEBUG)
85  << "No nuclear target was found - No need to simulate formation zones";
86  return;
87  }
88  // Compute the nuclear radius & how far away a particle is being tracked by
89  // the intranuclear hadron transport
90  assert(nucltgt && pdg::IsIon(nucltgt->Pdg()));
91  double A = nucltgt->A();
92  double R = fR0 * TMath::Power(A, 1./3.);
93  R *= TMath::Max(fNR,1.); // particle is tracked much further outside the nuclear boundary as the density is non-zero
94 
95  // Decay very short living particles so that we can give the formation
96  // zone to the daughters
97  this->PreHadronTransportDecays(evrec);
98 
99  // Get hadronic system's 3-momentum
100  GHepParticle * hadronic_system = evrec->FinalStateHadronicSystem();
101  TVector3 p3hadr = hadronic_system->P4()->Vect(); // (px,py,pz)
102 
103  // Loop over GHEP and set the formation zone to the right particles
104  // Limit the maximum formation zone so that particles escaping the
105  // nucleus are placed right outside...
106 
107  TObjArrayIter piter(evrec);
108  GHepParticle * p = 0;
109  int icurr = -1;
110 
111  while( (p = (GHepParticle *) piter.Next()) )
112  {
113  icurr++;
114 
115  GHepStatus_t ist = p->Status();
116  bool apply_formation_zone = (ist==kIStHadronInTheNucleus);
117 
118  if(!apply_formation_zone) continue;
119 
120  LOG("DISHadronicVtx", pINFO)
121  << "Applying formation-zone to " << p->Name();
122 
123  double m = p->Mass();
124  int pdgc = p->Pdg();
125  const TLorentzVector & p4 = *(p->P4());
126  double ct0=0.;
127  pdg::IsNucleon(pdgc) ? ct0=fct0nucleon : ct0=fct0pion;
128  double fz = phys::FormationZone(m,p4,p3hadr,ct0,fK);
129 
130  //-- Apply the formation zone step
131 
132  double step = fz;
133  TVector3 dr = p->P4()->Vect().Unit(); // unit vector along its direction
134  // double c = kLightSpeed / (units::fm/units::ns); // c in fm/nsec
135  dr.SetMag(step); // spatial step size
136  // double dt = step/c; // temporal step:
137  double dt = 0;
138  TLorentzVector dx4(dr,dt); // 4-vector step
139  TLorentzVector x4new = *(p->X4()) + dx4; // new position
140 
141  //-- If the formation zone was large enough that the particle is now outside
142  // the nucleus make sure that it is not placed further away from the
143  // (max distance particles tracked by intranuclear cascade) + ~2 fm
144  double epsilon = 2; // fm
145  double r = x4new.Vect().Mag(); // fm
146  double rmax = R+epsilon;
147  if(r > rmax) {
148  LOG("DISHadronicVtx", pINFO)
149  << "Particle was stepped too far away (r = " << r << " fm)";
150  LOG("DISHadronicVtx", pINFO)
151  << "Placing it ~2 fm away from the furthermost position tracked "
152  << "by intranuclear cascades (r' = " << rmax << " fm)";
153  double scale = rmax/r;
154  x4new *= scale;
155  }
156 
157  p->SetPosition(x4new);
158  }
159 }
160 //___________________________________________________________________________
162 {
163  Algorithm::Configure(config);
164  this->LoadConfig();
165 }
166 //____________________________________________________________________________
168 {
169  Algorithm::Configure(config);
170  this->LoadConfig();
171 }
172 //____________________________________________________________________________
174 {
176  fPreINukeDecayer = 0;
177 
178  //-- Get the requested hadronization model
180  dynamic_cast<const EventRecordVisitorI *> (this->SubAlg("Hadronizer"));
181  assert(fHadronizationModel);
182 
183  //-- Handle pre-intranuke decays
185  dynamic_cast<const EventRecordVisitorI *> (this->SubAlg("PreTransportDecayer"));
186  assert(fPreINukeDecayer);
187 
188  //-- Get flag to determine whether we copy all fragmentation record entries
189  // into the GHEP record or just the ones marked with kf=1
190  GetParamDef( "FilterPreFragm", fFilterPreFragmEntries, false ) ;
191 
192  //-- Get parameters controlling the nuclear sizes
193  GetParam( "NUCL-R0", fR0 ) ;
194  GetParam( "NUCL-NR", fNR ) ;
195 
196  //-- Get parameters controlling the formation zone simulation
197  GetParam( "FZONE-ct0pion", fct0pion ) ;
198  GetParam( "FZONE-ct0nucleon",fct0nucleon ) ;
199  GetParam( "FZONE-KPt2", fK ) ;
200 
201  LOG("DISHadronicVtx", pDEBUG) << "ct0pion = " << fct0pion << " fermi";
202  LOG("DISHadronicVtx", pDEBUG) << "ct0nucleon = " << fct0nucleon << " fermi";
203  LOG("DISHadronicVtx", pDEBUG) << "K(pt^2) = " << fK;
204 }
205 //____________________________________________________________________________
double fct0pion
formation zone (c * formation time) - for pions
TLorentzVector Hadronic4pLAB(GHepRecord *event_rec) const
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
void ProcessEventRecord(GHepRecord *event_rec) const
bool IsNucleon(int pdgc)
Definition: PDGUtils.cxx:346
const EventRecordVisitorI * fPreINukeDecayer
double fK
param multiplying pT^2 in formation zone calculation
double Mass(void) const
Mass that corresponds to the PDG code.
const EventRecordVisitorI * fHadronizationModel
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
double FormationZone(double m, const TLorentzVector &p, const TVector3 &p3hadr, double ct0, double K)
Definition: PhysUtils.cxx:18
const double epsilon
int Pdg(void) const
Definition: GHepParticle.h:63
double fNR
how far beyond the nuclear boundary does the particle tracker goes?
string Name(void) const
Name that corresponds to the PDG code.
Summary information for an interaction.
Definition: Interaction.h:56
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:96
static constexpr double A
Definition: Units.h:74
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
double fct0nucleon
formation zone (c * formation time) - for nucleons
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
#define pINFO
Definition: Messenger.h:62
virtual GHepParticle * FinalStateHadronicSystem(void) const
Definition: GHepRecord.cxx:342
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
bool IsIon(int pdgc)
Definition: PDGUtils.cxx:42
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void SimulateFormationZone(GHepRecord *event_rec) const
Abstract class. Is used to pass some commonly recurring methods to all concrete implementations of th...
int A(void) const
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
void AddFinalHadronicSyst(GHepRecord *event_rec) const
double fR0
param controling nuclear size
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
static constexpr double m
Definition: Units.h:71
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:39
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63
void PreHadronTransportDecays(GHepRecord *event_rec) const
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345