Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TruthTrackingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TruthTrackingAction.cc
2 
3 #include "PHG4Particle.h" // for PHG4Particle
4 #include "PHG4Particlev2.h"
5 #include "PHG4Particlev3.h"
6 #include "PHG4Shower.h" // for PHG4Shower
7 #include "PHG4Showerv1.h"
8 #include "PHG4TrackUserInfoV1.h"
9 #include "PHG4TruthEventAction.h"
10 #include "PHG4TruthInfoContainer.h"
12 #include "PHG4VtxPointv1.h"
13 
14 #include <phool/getClass.h>
15 
16 #include <Geant4/G4DynamicParticle.hh> // for G4DynamicParticle
17 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
18 #include <Geant4/G4PrimaryParticle.hh>
19 #include <Geant4/G4SystemOfUnits.hh>
20 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
21 #include <Geant4/G4Track.hh>
22 #include <Geant4/G4TrackVector.hh> // for G4TrackVector
23 #include <Geant4/G4TrackingManager.hh>
24 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
25 
26 #include <cmath> // for sqrt
27 #include <cstddef> // for size_t
28 #include <iostream> // for operator<<, endl
29 #include <utility> // for pair
30 
31 using namespace std;
32 
33 //________________________________________________________
35  : m_EventAction(eventAction)
36  , m_TruthInfoList(nullptr)
37  , m_G4ParticleStack()
38  , m_CurrG4Particle()
39 {
40 }
41 
43 {
44  // insert particle into the output
45  PHG4Particle* ti = AddParticle(*m_TruthInfoList, *const_cast<G4Track*>(track));
46 
47  // Negative G4 track id values indicate unwanted tracks to be deleted
48  // Initially all tracks except primary ones flagged as unwanted
49  int track_id_g4 = track->GetTrackID() * (track->GetParentID() ? -1 : +1);
50  int trackid = ti->get_track_id();
51 
52  // add the user id to the geant4 user info
53  PHG4TrackUserInfo::SetUserTrackId(const_cast<G4Track*>(track), trackid);
54 
55  if (PHG4TrackUserInfoV1* p = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation()))
56  {
57  ti->set_parent_id(p->GetUserParentId());
58 
59  if (track->GetParentID())
60  {
61  ti->set_primary_id(p->GetUserPrimaryId());
62  }
63  else
64  {
65  PHG4TrackUserInfo::SetUserPrimaryId(const_cast<G4Track*>(track), trackid);
66  ti->set_primary_id(trackid);
67  }
68  }
69 
70  if (!track->GetParentID())
71  {
72  // primary track - propagate the barcode information
73  PHG4UserPrimaryParticleInformation* userdata = static_cast<PHG4UserPrimaryParticleInformation*>(track->GetDynamicParticle()->GetPrimaryParticle()->GetUserInformation());
74  if (userdata) ti->set_barcode(userdata->get_user_barcode());
75  }
76 
77  int vtxindex = ti->get_vtx_id();
78 
79  m_CurrG4Particle = {track_id_g4, trackid, vtxindex};
80 
81  // create or add to a new shower object --------------------------------------
82  if (!track->GetParentID())
83  {
85  PHG4TrackUserInfo::SetShower(const_cast<G4Track*>(track), shower);
86  m_TruthInfoList->AddShower(trackid, shower);
87  shower->set_id(trackid); // fyi, secondary showers may not share these ids
88  shower->set_parent_particle_id(trackid);
89  shower->set_parent_shower_id(0);
90  }
91  else
92  {
93  // get shower
94  if (G4VUserTrackInformation* p = track->GetUserInformation())
95  {
96  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
97  {
98  if (pp->GetShower())
99  {
100  pp->GetShower()->add_g4particle_id(trackid);
101  pp->GetShower()->add_g4vertex_id(vtxindex);
102  }
103  }
104  }
105  }
106 
107  // tell the primary particle copy in G4 where this output will be stored
108  if (!track->GetParentID())
109  {
110  PHG4UserPrimaryParticleInformation* userdata = static_cast<PHG4UserPrimaryParticleInformation*>(track->GetDynamicParticle()->GetPrimaryParticle()->GetUserInformation());
111  if (userdata)
112  {
113  userdata->set_user_track_id(trackid);
114  userdata->set_user_vtx_id(vtxindex);
115  }
116  }
117 
118  return;
119 }
120 
122 {
123  if (fpTrackingManager)
124  {
125  int trackid = track->GetTrackID();
126  int primaryid = 0;
127  PHG4Shower* shower = nullptr;
128  if (PHG4TrackUserInfoV1* p = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation()))
129  {
130  trackid = p->GetUserTrackId();
131  primaryid = p->GetUserPrimaryId();
132  shower = p->GetShower();
133  }
134 
135  G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries();
136  if (secondaries)
137  {
138  for (size_t i = 0; i < secondaries->size(); ++i)
139  {
140  G4Track* secondary = (*secondaries)[i];
141  PHG4TrackUserInfo::SetUserParentId(const_cast<G4Track*>(secondary), trackid);
142  PHG4TrackUserInfo::SetUserPrimaryId(const_cast<G4Track*>(secondary), primaryid);
143  PHG4TrackUserInfo::SetShower(const_cast<G4Track*>(secondary), shower);
144  }
145  }
146  }
147 
148  if (PHG4TrackUserInfoV1* p = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation()))
149  {
150  if (p->GetKeep())
151  {
152  int trackid = p->GetUserTrackId();
154  }
155  }
156 
157  UpdateG4ParticleStack(track);
158 }
159 
165 {
166  while (!m_G4ParticleStack.empty())
167  {
168  if ( std::abs(m_G4ParticleStack.back().g4track_id) == track->GetParentID() )
169  {
170  break;
171  }
172  else
173  {
174  if (m_G4ParticleStack.back().g4track_id < 0)
175  {
176  m_TruthInfoList->delete_particle( m_G4ParticleStack.back().particle_id );
177  }
178  m_G4ParticleStack.pop_back();
179  }
180  }
181 
183 
184  // Change sign of G4 track id of all upstream tracks in the stack to positive
185  // in order to keep the track
186  PHG4TrackUserInfoV1* p = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation());
187  bool keep_curr_track = p && p->GetKeep() ? true : false;
188 
189  auto stack_iter = m_G4ParticleStack.rbegin();
190  while (keep_curr_track && stack_iter != m_G4ParticleStack.rend() && stack_iter->g4track_id < 0)
191  {
192  stack_iter->g4track_id *= -1;
193  ++stack_iter;
194  }
195 }
196 
198 {
199  //now look for the map and grab a pointer to it.
200  m_TruthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
201 
202  // if we do not find the node we need to make it.
203  if (!m_TruthInfoList)
204  {
205  std::cout << "PHG4TruthEventAction::SetInterfacePointers - unable to find G4TruthInfo" << std::endl;
206  }
207 }
208 
210 {
211  m_VertexMap.clear();
212 
213  while (!m_G4ParticleStack.empty())
214  {
215  if ( m_G4ParticleStack.back().g4track_id < 0)
216  {
217  m_TruthInfoList->delete_particle( m_G4ParticleStack.back().particle_id );
218  }
219  m_G4ParticleStack.pop_back();
220  }
221 
222  return 0;
223 }
224 
226 {
227  int trackid = 0;
228  if (track.GetParentID())
229  {
230  // secondaries get negative user ids and increment downward between geant subevents
231  trackid = truth.mintrkindex() - 1;
232  }
233  else
234  {
235  // primaries get positive user ids and increment upward between geant subevents
236  trackid = truth.maxtrkindex() + 1;
237  }
238 
239  // determine the momentum vector
240  G4ParticleDefinition* def = track.GetDefinition();
241  int pdgid = def->GetPDGEncoding();
242  double mass = def->GetPDGMass();
243  double ke = track.GetVertexKineticEnergy();
244  double ptot = sqrt(ke * ke + 2.0 * mass * ke);
245  G4ThreeVector pdir = track.GetVertexMomentumDirection();
246  pdir *= ptot;
247  PHG4Particle* ti = nullptr;
248  // create a new particle -----------------------------------------------------
249  if (def->IsGeneralIon()) // for ions save a and z in v3 of phg4particle
250  {
251  ti = new PHG4Particlev3();
252  ti->set_A(def->GetAtomicMass());
253  ti->set_Z(def->GetAtomicNumber());
254  }
255  else
256  {
257  ti = new PHG4Particlev2;
258  }
259  ti->set_px(pdir[0] / GeV);
260  ti->set_py(pdir[1] / GeV);
261  ti->set_pz(pdir[2] / GeV);
262  ti->set_track_id(trackid);
263 
264  ti->set_parent_id(track.GetParentID());
265  ti->set_primary_id(trackid);
266 
267  ti->set_pid(pdgid);
268  ti->set_name(def->GetParticleName());
269  ti->set_e(track.GetTotalEnergy() / GeV);
270 
271  // Add new or reuse a vertex and let the track know about it
272  PHG4VtxPoint* vtx = AddVertex(truth, track);
273  ti->set_vtx_id(vtx->get_id());
274 
275  return truth.AddParticle(trackid, ti)->second;
276 }
277 
279 {
280  G4ThreeVector v = track.GetVertexPosition();
281  int vtxindex = (track.GetParentID() == 0 ? truth.maxvtxindex() + 1 : truth.minvtxindex() - 1);
282 
283  auto [iter, inserted] = m_VertexMap.insert(std::make_pair(v, vtxindex));
284 
285  // If could not add a unique vertex => return the existing one
286  if (!inserted)
287  {
288  return truth.GetVtxMap().find(iter->second)->second;
289  }
290  // otherwise, create and add a new one
291  PHG4VtxPoint* vtxpt = new PHG4VtxPointv1(v[0]/cm, v[1]/cm, v[2]/cm, track.GetGlobalTime()/ns, vtxindex);
292 
293  return truth.AddVertex(vtxindex, vtxpt)->second;
294 }