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