Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TruthEventAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TruthEventAction.cc
1 #include "PHG4TruthEventAction.h"
2 
3 #include "PHG4Hit.h"
4 #include "PHG4HitContainer.h"
5 #include "PHG4HitDefs.h" // for keytype
6 #include "PHG4Particle.h" // for PHG4Particle
7 #include "PHG4Shower.h"
10 
11 #include <phool/PHCompositeNode.h>
12 #include <phool/PHIODataNode.h> // for PHIODataNode
13 #include <phool/PHNode.h>
14 #include <phool/PHNodeIterator.h> // for PHNodeIterator
16 #include <phool/getClass.h>
17 #include <phool/phool.h> // for PHWHERE
18 
19 #include <Geant4/G4Event.hh>
20 #include <Geant4/G4PrimaryParticle.hh> // for G4PrimaryPart...
21 #include <Geant4/G4PrimaryVertex.hh> // for G4PrimaryVertex
22 #include <Geant4/G4VUserPrimaryParticleInformation.hh> // for G4VUserPrimar...
23 
24 // eigen has some shadowed variables
25 #pragma GCC diagnostic push
26 #pragma GCC diagnostic ignored "-Wshadow"
27 #include <Eigen/Dense>
28 #pragma GCC diagnostic pop
29 
30 #include <cmath> // for isnan
31 #include <cstdlib> // for abs
32 #include <iostream> // for operator<<, endl
33 #include <iterator> // for reverse_iterator
34 #include <map>
35 #include <string> // for operator==
36 #include <utility> // for swap, pair
37 #include <vector> // for vector, vecto...
38 
39 class PHG4VtxPoint;
40 
41 using namespace std;
42 
43 //___________________________________________________
45  : m_TruthInfoContainer(nullptr)
46  , m_LowerKeyPrevExist(0)
47  , m_UpperKeyPrevExist(0)
48 {
49 }
50 
51 //___________________________________________________
52 void PHG4TruthEventAction::BeginOfEventAction(const G4Event* /*evt*/)
53 {
54  // if we do not find the node we need to make it.
56  {
57  std::cout << "PHG4TruthEventAction::EndOfEventAction - unable to find G4TruthInfo node" << std::endl;
58  return;
59  }
60 
62  if (!map.empty())
63  {
64  m_LowerKeyPrevExist = map.begin()->first;
65  m_UpperKeyPrevExist = map.rbegin()->first;
66  }
67 }
68 
69 //___________________________________________________
71 {
72  // if we do not find the node we need to make it.
74  {
75  std::cout << "PHG4TruthEventAction::EndOfEventAction - unable to find G4TruthInfo node" << std::endl;
76  return;
77  }
78  // First deal with the showers - they do need the info which
79  // is removed from the maps in the subsequent cleanup to reduce the
80  // output file size
81  PruneShowers();
83  // construct a list of track ids to preserve in the the output that includes any
84  // track designated in the m_WriteSet during processing or its ancestry chain
85 
86  std::set<int> savelist;
87  std::set<int> savevtxlist;
88 
89  for (std::set<int>::const_iterator write_iter = m_WriteSet.begin();
90  write_iter != m_WriteSet.end();
91  ++write_iter)
92  {
93  std::vector<int> wrttracks;
94  std::vector<int> wrtvtx;
95 
96  // usertrackid
97  int mytrkid = *write_iter;
99 
100  // if track is already in save list, nothing needs to be done
101  if (savelist.find(mytrkid) != savelist.end())
102  {
103  continue;
104  }
105  else
106  {
107  wrttracks.push_back(mytrkid);
108  wrtvtx.push_back(particle->get_vtx_id());
109  }
110 
111  // now crawl up the truth info and add parents until we hit
112  // a track which is already being saved
113  int parentid = particle->get_parent_id();
114  while (savelist.find(parentid) == savelist.end() && parentid != 0)
115  {
116  particle = m_TruthInfoContainer->GetParticle(parentid);
117  wrttracks.push_back(parentid);
118  wrtvtx.push_back(particle->get_vtx_id());
119  parentid = particle->get_parent_id();
120  }
121 
122  // now fill the new tracks into the save list
123  // while emptying the write lists
124  while (wrttracks.begin() != wrttracks.end())
125  {
126  savelist.insert(wrttracks.back());
127  wrttracks.pop_back();
128  }
129 
130  while (wrtvtx.begin() != wrtvtx.end())
131  {
132  savevtxlist.insert(wrtvtx.back());
133  wrtvtx.pop_back();
134  }
135  }
136 
137  // the save lists are filled now, except primary track which never
138  // made it into any active volume and their vertex
139 
140  // loop over particles in truth list and remove them if they are not
141  // in the save list and are not primary particles (parent_id == 0)
142 
143  int removed[4] = {0};
145  PHG4TruthInfoContainer::Iterator truthiter = truth_range.first;
146  while (truthiter != truth_range.second)
147  {
148  removed[0]++;
149  int trackid = truthiter->first;
150  if (savelist.find(trackid) == savelist.end())
151  {
152  // track not in save list
153 
154  // check if particle below offset
155  // primary particles had parentid = 0
156  // for embedding: particles from initial sim do not have their keep flag set, we want to keep particles with trkid <= trackidoffset
157  // tracks from a previous geant pass will not be recorded as leaving
158  // hit in the sim, so we exclude this range from the removal
159  // for regular sims, that range is zero to zero
160  if (((trackid < m_LowerKeyPrevExist) || (trackid > m_UpperKeyPrevExist)) && ((truthiter->second)->get_parent_id() != 0))
161  {
163  removed[1]++;
164  }
165  else
166  {
167  // save vertex id for primary particle which leaves no hit
168  // in active area
169  savevtxlist.insert((truthiter->second)->get_vtx_id());
170  ++truthiter;
171  }
172  }
173  else
174  {
175  // track was in save list, move on
176  ++truthiter;
177  }
178  }
179 
181  PHG4TruthInfoContainer::VtxIterator vtxiter = vtxrange.first;
182  while (vtxiter != vtxrange.second)
183  {
184  removed[2]++;
185  if (savevtxlist.find(vtxiter->first) == savevtxlist.end())
186  {
187  m_TruthInfoContainer->delete_vtx(vtxiter++);
188  removed[3]++;
189  }
190  else
191  {
192  ++vtxiter;
193  }
194  }
195 
196  // loop over all input particles and fish out the ones which have the embed flag set
197  // and store their geant track ids in truthinfo container
198  G4PrimaryVertex* pvtx = evt->GetPrimaryVertex();
199  while (pvtx)
200  {
201  G4PrimaryParticle* part = pvtx->GetPrimary();
202  while (part)
203  {
204  PHG4UserPrimaryParticleInformation* userdata = dynamic_cast<PHG4UserPrimaryParticleInformation*>(part->GetUserInformation());
205  if (userdata)
206  {
207  if (userdata->get_embed())
208  {
211  }
212  }
213  part = part->GetNext();
214  }
215  pvtx = pvtx->GetNext();
216  }
217 
218  return;
219 }
220 
221 //___________________________________________________
223 {
224  m_WriteSet.insert(trackid);
225 }
226 
227 //___________________________________________________
229 {
230  //now look for the map and grab a pointer to it.
231  m_TruthInfoContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
232 
233  // if we do not find the node we need to make it.
235  {
236  std::cout << "PHG4TruthEventAction::SetInterfacePointers - unable to find G4TruthInfo" << std::endl;
237  }
238 
239  SearchNode(topNode);
240 }
241 
243 {
244  // fill a lookup map between the g4hit container ids and the containers themselves
245  // without knowing what the container names are in advance, only that they
246  // begin G4HIT_*
247 
248  PHNodeIterator nodeiter(top);
249  PHPointerListIterator<PHNode> iter(nodeiter.ls());
250  PHNode* thisNode;
251  while ((thisNode = iter()))
252  {
253  if (thisNode->getType() == "PHCompositeNode")
254  {
255  SearchNode(static_cast<PHCompositeNode*>(thisNode));
256  }
257  else if (thisNode->getType() == "PHIODataNode")
258  {
259  if (thisNode->getName().find("G4HIT_") == 0)
260  {
261  PHIODataNode<PHG4HitContainer>* DNode = static_cast<PHIODataNode<PHG4HitContainer>*>(thisNode);
262  if (DNode)
263  {
264  PHG4HitContainer* object = dynamic_cast<PHG4HitContainer*>(DNode->getData());
265  if (object)
266  {
267  m_HitContainerMap[object->GetID()] = object;
268  }
269  }
270  }
271  }
272  }
273 }
274 
276 {
277  m_WriteSet.clear();
278  return 0;
279 }
280 
282 {
284  for (PHG4TruthInfoContainer::ShowerIterator iter = range.first;
285  iter != range.second;
286  ++iter)
287  {
288  PHG4Shower* shower = iter->second;
289 
290  std::set<int> remove_ids;
291  for (PHG4Shower::ParticleIdIter jter = shower->begin_g4particle_id();
292  jter != shower->end_g4particle_id();
293  ++jter)
294  {
295  int g4particle_id = *jter;
297  if (!particle)
298  {
299  remove_ids.insert(g4particle_id);
300  continue;
301  }
302  }
303 
304  for (std::set<int>::iterator jter = remove_ids.begin();
305  jter != remove_ids.end();
306  ++jter)
307  {
308  shower->remove_g4particle_id(*jter);
309  }
310 
311  std::set<int> remove_more_ids;
312  for (std::map<int, std::set<PHG4HitDefs::keytype> >::iterator jter = shower->begin_g4hit_id();
313  jter != shower->end_g4hit_id();
314  ++jter)
315  {
316  int g4hitmap_id = jter->first;
317  std::map<int, PHG4HitContainer*>::iterator mapiter = m_HitContainerMap.find(g4hitmap_id);
318  if (mapiter == m_HitContainerMap.end())
319  {
320  continue;
321  }
322 
323  // get the g4hits from this particle in this volume
324  for (std::set<PHG4HitDefs::keytype>::iterator kter = jter->second.begin();
325  kter != jter->second.end();)
326  {
327  PHG4HitDefs::keytype g4hit_id = *kter;
328 
329  PHG4Hit* g4hit = mapiter->second->findHit(g4hit_id);
330  if (!g4hit)
331  {
332  // some zero edep g4hits have been removed already
333  jter->second.erase(kter++);
334  continue;
335  }
336  else
337  {
338  ++kter;
339  }
340  }
341 
342  if (jter->second.empty())
343  {
344  remove_more_ids.insert(g4hitmap_id);
345  }
346  }
347 
348  for (std::set<int>::iterator jter = remove_more_ids.begin();
349  jter != remove_more_ids.end();
350  ++jter)
351  {
352  shower->remove_g4hit_volume(*jter);
353  }
354  }
355 
357  for (PHG4TruthInfoContainer::ShowerIterator iter = range.first;
358  iter != range.second;)
359  {
360  PHG4Shower* shower = iter->second;
361 
362  if (shower->empty_g4particle_id() && shower->empty_g4hit_id())
363  {
364  if (shower->get_edep() == 0) // check whether this shower has already been processed in the previous simulation cycles
365  {
367  continue;
368  }
369  }
370 
371  ++iter;
372  }
373 }
374 
376 {
378  for (PHG4TruthInfoContainer::ShowerIterator shwiter = range.first;
379  shwiter != range.second;
380  ++shwiter)
381  {
382  PHG4Shower* shower = shwiter->second;
383 
384  // Data structures to hold weighted pca
385  std::vector<std::vector<float> > points;
386  std::vector<float> weights;
387  float sumw = 0.0;
388  float sumw2 = 0.0;
389 
390  for (std::map<int, std::set<PHG4HitDefs::keytype> >::iterator iter = shower->begin_g4hit_id();
391  iter != shower->end_g4hit_id();
392  ++iter)
393  {
394  int g4hitmap_id = iter->first;
395  std::map<int, PHG4HitContainer*>::iterator mapiter = m_HitContainerMap.find(g4hitmap_id);
396  if (mapiter == m_HitContainerMap.end())
397  {
398  continue;
399  }
400 
401  PHG4HitContainer* hits = mapiter->second;
402 
403  unsigned int nhits = 0;
404  float edep = 0.0;
405  float eion = 0.0;
406  float light_yield = 0.0;
407  float edep_e = 0.0;
408  float edep_h = 0.0;
409 
410  // get the g4hits from this particle in this volume
411  for (std::set<PHG4HitDefs::keytype>::iterator kter = iter->second.begin();
412  kter != iter->second.end();
413  ++kter)
414  {
415  PHG4HitDefs::keytype g4hit_id = *kter;
416 
417  PHG4Hit* g4hit = hits->findHit(g4hit_id);
418  if (!g4hit)
419  {
420  cout << PHWHERE << " missing g4hit" << endl;
421  continue;
422  }
423 
425  if (!particle)
426  {
427  cout << PHWHERE << " missing g4particle for track "
428  << g4hit->get_trkid() << endl;
429  continue;
430  }
431 
432  PHG4VtxPoint* vtx = m_TruthInfoContainer->GetVtx(particle->get_vtx_id());
433  if (!vtx)
434  {
435  cout << PHWHERE << " missing g4vertex" << endl;
436  continue;
437  }
438 
439  // shower location and shape info
440 
441  if (!isnan(g4hit->get_x(0)) &&
442  !isnan(g4hit->get_y(0)) &&
443  !isnan(g4hit->get_z(0)))
444  {
445  std::vector<float> entry(3);
446  entry[0] = g4hit->get_x(0);
447  entry[1] = g4hit->get_y(0);
448  entry[2] = g4hit->get_z(0);
449 
450  points.push_back(entry);
451  float w = g4hit->get_edep();
452  weights.push_back(w);
453  sumw += w;
454  sumw2 += w * w;
455  }
456 
457  if (!isnan(g4hit->get_x(1)) &&
458  !isnan(g4hit->get_y(1)) &&
459  !isnan(g4hit->get_z(1)))
460  {
461  std::vector<float> entry(3);
462  entry[0] = g4hit->get_x(1);
463  entry[1] = g4hit->get_y(1);
464  entry[2] = g4hit->get_z(1);
465 
466  points.push_back(entry);
467  float w = g4hit->get_edep();
468  weights.push_back(w);
469  sumw += w;
470  sumw2 += w * w;
471  }
472 
473  // e/h ratio
474 
475  if (!isnan(g4hit->get_edep()))
476  {
477  if (abs(particle->get_pid()) == 11)
478  {
479  edep_e += g4hit->get_edep();
480  }
481  else
482  {
483  edep_h += g4hit->get_edep();
484  }
485  }
486 
487  // summary info
488 
489  ++nhits;
490  if (!isnan(g4hit->get_edep())) edep += g4hit->get_edep();
491  if (!isnan(g4hit->get_eion())) eion += g4hit->get_eion();
492  if (!isnan(g4hit->get_light_yield())) light_yield += g4hit->get_light_yield();
493  } // g4hit loop
494 
495  // summary info
496 
497  if (nhits) shower->set_nhits(g4hitmap_id, nhits);
498  if (edep != 0.0) shower->set_edep(g4hitmap_id, edep);
499  if (eion != 0.0) shower->set_eion(g4hitmap_id, eion);
500  if (light_yield != 0.0) shower->set_light_yield(g4hitmap_id, light_yield);
501  if (edep_h != 0.0) shower->set_eh_ratio(g4hitmap_id, edep_e / edep_h);
502  } // volume loop
503 
504  // fill Eigen matrices to compute wPCA
505  // resizing these non-destructively is expensive
506  // so I fill vectors and then copy
507  Eigen::Matrix<double, Eigen::Dynamic, 3> X(points.size(), 3);
508  Eigen::Matrix<double, Eigen::Dynamic, 1> W(weights.size(), 1);
509 
510  for (unsigned int i = 0; i < points.size(); ++i)
511  {
512  for (unsigned int j = 0; j < 3; ++j)
513  {
514  X(i, j) = points[i][j];
515  }
516  W(i, 0) = weights[i];
517  }
518 
519  // mean value of shower
520  double prefactor = 1.0 / sumw;
521  Eigen::Matrix<double, 1, 3> mean = prefactor * W.transpose() * X;
522 
523  // compute residual relative to the mean
524  for (unsigned int i = 0; i < points.size(); ++i)
525  {
526  for (unsigned int j = 0; j < 3; ++j) X(i, j) = points[i][j] - mean(0, j);
527  }
528 
529  // weighted covariance matrix
530  prefactor = sumw / (sumw*sumw - sumw2); // effectivelly 1/(N-1) when w_i = 1.0
531  Eigen::Matrix<double, 3, 3> covar = prefactor * (X.transpose() * W.asDiagonal() * X);
532 
533  shower->set_x(mean(0, 0));
534  shower->set_y(mean(0, 1));
535  shower->set_z(mean(0, 2));
536 
537  for (unsigned int i = 0; i < 3; ++i)
538  {
539  for (unsigned int j = 0; j <= i; ++j)
540  {
541  shower->set_covar(i, j, covar(i, j));
542  }
543  }
544 
545  // shower->identify();
546  }
547 }