DISHadronicSystemGenerator.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  @ Dec 03, 2007 - CA
14  If the formation zone is too large then the particle is placed back,
15  just "outside the nucleus". Update the definition of "outside the
16  nucleus" to be the maximum distance that a particle can be tracked
17  by the intranuclear cascade + a couple of fermis. Brings the code
18  in sync with changes in the intranuclear cascade tracking algorithm.
19  @ Feb 07, 2009 - CA
20  Removed call to AddTargetNucleusRemnant(). This simulation step is now
21  performed further upstream in the processing chain.
22  @ Mar 03, 2009 - CA
23  Moved into the new DIS package from its previous location (EVGModules).
24  @ Sep 15, 2009 - CA
25  IsNucleus() is no longer available in GHepParticle. Use pdg::IsIon().
26  @ Feb 08, 2013 - CA
27  Use the formation zone code from PhysUtils (also used by reweighting)
28  rather than having own implementation here
29 */
30 //____________________________________________________________________________
31 
32 #include <RVersion.h>
33 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6)
34 #include <TMCParticle.h>
35 #else
36 #include <TMCParticle6.h>
37 #endif
38 
57 
58 using namespace genie;
59 using namespace genie::controls;
60 using namespace genie::constants;
61 using namespace genie::utils;
62 
63 //___________________________________________________________________________
65 HadronicSystemGenerator("genie::DISHadronicSystemGenerator")
66 {
67 
68 }
69 //___________________________________________________________________________
71 HadronicSystemGenerator("genie::DISHadronicSystemGenerator", config)
72 {
73 
74 }
75 //___________________________________________________________________________
77 {
78 
79 }
80 //___________________________________________________________________________
82 {
83 // This method generates the final state hadronic system
84 
85  //-- Add an entry for the DIS Pre-Fragm. Hadronic State
86  this->AddFinalHadronicSyst(evrec);
87 
88  //-- Add the fragmentation products
89  this->AddFragmentationProducts(evrec);
90 
91  //-- Simulate the formation zone if not taken directly from the
92  // hadronization model
93  this->SimulateFormationZone(evrec);
94 }
95 //___________________________________________________________________________
97  GHepRecord * evrec) const
98 {
99 // Calls a hadronizer and adds the fragmentation products at the GHEP
100 
101  GHepParticle * neutrino = evrec->Probe();
102  const TLorentzVector & vtx = *(neutrino->X4());
103 
104  //-- Compute the hadronic system invariant mass
105  TLorentzVector p4Had = this->Hadronic4pLAB(evrec);
106  double W = p4Had.M();
107 
108  Interaction * interaction = evrec->Summary();
109  interaction->KinePtr()->SetW(W);
110 
111  //-- Run the hadronization model and get the fragmentation products:
112  // A collection of ROOT TMCParticles (equiv. to a LUJETS record)
113 
114  TClonesArray * plist = fHadronizationModel->Hadronize(interaction);
115  if(!plist) {
116  LOG("DISHadronicVtx", pWARN)
117  << "Got an empty particle list. Hadronizer failed!";
118  LOG("DISHadronicVtx", pWARN)
119  << "Quitting the current event generation thread";
120 
121  evrec->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
122 
124  exception.SetReason("Could not simulate the hadronic system");
125  exception.SwitchOnFastForward();
126  throw exception;
127 
128  return;
129  }
130 
131  //-- Take the hadronic system weight to handle cases that the hadronizer
132  // was asked to produce weighted events
133  double wght = fHadronizationModel->Weight();
134 
135  //-- Translate the fragmentation products from TMCParticles to
136  // GHepParticles and copy them to the event record.
137 
138  int mom = evrec->FinalStateHadronicSystemPosition();
139  assert(mom!=-1);
140 
141  TMCParticle * p = 0;
142  TIter particle_iter(plist);
143 
144  bool is_nucleus = interaction->InitState().Tgt().IsNucleus();
145  GHepStatus_t istfin = (is_nucleus) ?
147 
148  // Vector defining rotation from LAB to LAB' (z:= \vec{phad})
149  TVector3 unitvq = p4Had.Vect().Unit();
150 
151  // Boost velocity LAB' -> HCM
152  TVector3 beta(0,0,p4Had.P()/p4Had.Energy());
153 
154  while( (p = (TMCParticle *) particle_iter.Next()) ) {
155 
156  int pdgc = p->GetKF();
157  int ks = p->GetKS();
158 
159  if(fFilterPreFragmEntries && ks!=1) continue;
160 
161  // The fragmentation products are generated in the hadronic CM frame
162  // where the z>0 axis is the \vec{phad} direction. For each particle
163  // returned by the hadronizer:
164  // - boost it back to LAB' frame {z:=\vec{phad}} / doesn't affect pT
165  // - rotate its 3-momentum from LAB' to LAB
166 
167  TLorentzVector p4o(p->GetPx(), p->GetPy(), p->GetPz(), p->GetEnergy());
168  p4o.Boost(beta);
169  TVector3 p3 = p4o.Vect();
170  p3.RotateUz(unitvq);
171  TLorentzVector p4(p3,p4o.Energy());
172 
173  // copy final state particles to the event record
174  GHepStatus_t ist = (ks==1) ? istfin : kIStDISPreFragmHadronicState;
175 
176  // handle gammas, and leptons that might come from internal pythia decays
177  // mark them as final state particles
178  bool not_hadron = (pdgc==kPdgGamma ||
180  if(not_hadron) { ist = kIStStableFinalState; }
181 
182  int im = mom + 1 + p->GetParent();
183  int ifc = (p->GetFirstChild() == -1) ? -1 : mom + 1 + p->GetFirstChild();
184  int ilc = (p->GetLastChild() == -1) ? -1 : mom + 1 + p->GetLastChild();
185 
186  evrec->AddParticle(pdgc, ist, im,-1, ifc, ilc, p4,vtx);
187 
188  } // fragmentation-products-iterator
189 
190  //-- Handle the case that the hadronizer produced weighted events and
191  // take into account that the current event might be already weighted
192  evrec->SetWeight (wght * evrec->Weight());
193 
194  plist->Delete();
195  delete plist;
196 }
197 //___________________________________________________________________________
199  GHepRecord * evrec) const
200 {
201  LOG("DISHadronicVtx", pDEBUG)
202  << "Simulating formation zone for the DIS hadronic system";
203 
204  GHepParticle * nucltgt = evrec->TargetNucleus();
205  if (!nucltgt) {
206  LOG("DISHadronicVtx", pDEBUG)
207  << "No nuclear target was found - No need to simulate formation zones";
208  return;
209  }
210  // Compute the nuclear radius & how far away a particle is being tracked by
211  // the intranuclear hadron transport
212  assert(nucltgt && pdg::IsIon(nucltgt->Pdg()));
213  double A = nucltgt->A();
214  double R = fR0 * TMath::Power(A, 1./3.);
215  R *= TMath::Max(fNR,1.); // particle is tracked much further outside the nuclear boundary as the density is non-zero
216 
217  // Decay very short living particles so that we can give the formation
218  // zone to the daughters
219  this->PreHadronTransportDecays(evrec);
220 
221  // Get hadronic system's 3-momentum
222  GHepParticle * hadronic_system = evrec->FinalStateHadronicSystem();
223  TVector3 p3hadr = hadronic_system->P4()->Vect(); // (px,py,pz)
224 
225  // Loop over GHEP and set the formation zone to the right particles
226  // Limit the maximum formation zone so that particles escaping the
227  // nucleus are placed right outside...
228 
229  TObjArrayIter piter(evrec);
230  GHepParticle * p = 0;
231  int icurr = -1;
232 
233  while( (p = (GHepParticle *) piter.Next()) )
234  {
235  icurr++;
236 
237  GHepStatus_t ist = p->Status();
238  bool apply_formation_zone = (ist==kIStHadronInTheNucleus);
239 
240  if(!apply_formation_zone) continue;
241 
242  LOG("DISHadronicVtx", pINFO)
243  << "Applying formation-zone to " << p->Name();
244 
245  double m = p->Mass();
246  int pdgc = p->Pdg();
247  const TLorentzVector & p4 = *(p->P4());
248  double ct0=0.;
249  pdg::IsNucleon(pdgc) ? ct0=fct0nucleon : ct0=fct0pion;
250  double fz = phys::FormationZone(m,p4,p3hadr,ct0,fK);
251 
252  //-- Apply the formation zone step
253 
254  double step = fz;
255  TVector3 dr = p->P4()->Vect().Unit(); // unit vector along its direction
256  // double c = kLightSpeed / (units::fm/units::ns); // c in fm/nsec
257  dr.SetMag(step); // spatial step size
258  // double dt = step/c; // temporal step:
259  double dt = 0;
260  TLorentzVector dx4(dr,dt); // 4-vector step
261  TLorentzVector x4new = *(p->X4()) + dx4; // new position
262 
263  //-- If the formation zone was large enough that the particle is now outside
264  // the nucleus make sure that it is not placed further away from the
265  // (max distance particles tracked by intranuclear cascade) + ~2 fm
266  double epsilon = 2; // fm
267  double r = x4new.Vect().Mag(); // fm
268  double rmax = R+epsilon;
269  if(r > rmax) {
270  LOG("DISHadronicVtx", pINFO)
271  << "Particle was stepped too far away (r = " << r << " fm)";
272  LOG("DISHadronicVtx", pINFO)
273  << "Placing it ~2 fm away from the furthermost position tracked "
274  << "by intranuclear cascades (r' = " << rmax << " fm)";
275  double scale = rmax/r;
276  x4new *= scale;
277  }
278 
279  p->SetPosition(x4new);
280  }
281 }
282 //___________________________________________________________________________
284 {
285  Algorithm::Configure(config);
286  this->LoadConfig();
287 }
288 //____________________________________________________________________________
290 {
291  Algorithm::Configure(config);
292  this->LoadConfig();
293 }
294 //____________________________________________________________________________
296 {
298  fPreINukeDecayer = 0;
299 
300  //-- Get the requested hadronization model
302  dynamic_cast<const HadronizationModelI *> (this->SubAlg("Hadronizer"));
304 
305  //-- Handle pre-intranuke decays
307  dynamic_cast<const EventRecordVisitorI *> (this->SubAlg("PreTransportDecayer"));
309 
310  //-- Get flag to determine whether we copy all fragmentation record entries
311  // into the GHEP record or just the ones marked with kf=1
312  GetParamDef( "FilterPreFragm", fFilterPreFragmEntries, false ) ;
313 
314  //-- Get parameters controlling the nuclear sizes
315  GetParam( "NUCL-R0", fR0 ) ;
316  GetParam( "NUCL-NR", fNR ) ;
317 
318  //-- Get parameters controlling the formation zone simulation
319  GetParam( "FZONE-ct0pion", fct0pion ) ;
320  GetParam( "FZONE-ct0nucleon",fct0nucleon ) ;
321  GetParam( "FZONE-KPt2", fK ) ;
322 
323  LOG("DISHadronicVtx", pDEBUG) << "ct0pion = " << fct0pion << " fermi";
324  LOG("DISHadronicVtx", pDEBUG) << "ct0nucleon = " << fct0nucleon << " fermi";
325  LOG("DISHadronicVtx", pDEBUG) << "K(pt^2) = " << fK;
326 }
327 //____________________________________________________________________________
328 
Basic constants.
virtual void SetWeight(double wght)
Definition: GHepRecord.h:131
double fct0pion
formation zone (c * formation time) - for pions
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
TLorentzVector Hadronic4pLAB(GHepRecord *event_rec) const
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:135
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
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:309
enum genie::EGHepStatus GHepStatus_t
const char * p
Definition: xmltok.h:285
bool IsNucleus(void) const
Definition: Target.cxx:289
const EventRecordVisitorI * fPreINukeDecayer
double fK
param multiplying pT^2 in formation zone calculation
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:99
virtual double Weight(void) const
Definition: GHepRecord.h:125
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
Definition: config.py:1
double Mass(void) const
Mass that corresponds to the PDG code.
Double_t beta
GHepStatus_t Status(void) const
Definition: GHepParticle.h:65
double FormationZone(double m, const TLorentzVector &p, const TVector3 &p3hadr, double ct0, double K)
Definition: PhysUtils.cxx:29
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
int Pdg(void) const
Definition: GHepParticle.h:64
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
Definition: Cand.cxx:23
Double_t scale
Definition: plot.C:25
void SetPosition(const TLorentzVector &v4)
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const HadronizationModelI * fHadronizationModel
const int kPdgGamma
Definition: PDGCodes.h:166
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
double fct0nucleon
formation zone (c * formation time) - for nucleons
#define R(x)
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
virtual TClonesArray * Hadronize(const Interaction *) const =0
#define pINFO
Definition: Messenger.h:63
virtual GHepParticle * FinalStateHadronicSystem(void) const
Definition: GHepRecord.cxx:379
Misc GENIE control constants.
#define pWARN
Definition: Messenger.h:61
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:93
static const double A
Definition: Units.h:82
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:118
bool IsIon(int pdgc)
Definition: PDGUtils.cxx:40
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
void AddFragmentationProducts(GHepRecord *event_rec) const
double epsilon
Pure abstract base class. Defines the HadronizationModelI interface to be implemented by any algorith...
void SimulateFormationZone(GHepRecord *event_rec) const
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
const InitialState & InitState(void) const
Definition: Interaction.h:69
Abstract class. Is used to pass some commonly recurring methods to all concrete implementations of th...
std::initializer_list< hep_hpc::hdf5::PropertyList > plist
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
const Target & Tgt(void) const
Definition: InitialState.h:67
void AddFinalHadronicSyst(GHepRecord *event_rec) const
double fR0
param controling nuclear size
#define W(x)
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
virtual int FinalStateHadronicSystemPosition(void) const
Definition: GHepRecord.cxx:507
Root of GENIE utility namespaces.
virtual double Weight(void) const =0
#define pDEBUG
Definition: Messenger.h:64
void PreHadronTransportDecays(GHepRecord *event_rec) const
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353