ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4AllDstPileupMerger.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4AllDstPileupMerger.cc
1 
7 
8 #include "PHG4Hit.h" // for PHG4Hit
9 #include "PHG4HitContainer.h"
10 #include "PHG4Hitv1.h"
11 #include "PHG4Particle.h" // for PHG4Particle
12 #include "PHG4Particlev3.h"
13 #include "PHG4TruthInfoContainer.h"
14 #include "PHG4VtxPoint.h" // for PHG4VtxPoint
15 #include "PHG4VtxPointv1.h"
16 
17 #include <phhepmc/PHHepMCGenEvent.h> // for PHHepMCGenEvent
19 
20 #include <phool/PHCompositeNode.h>
21 #include <phool/PHIODataNode.h> // for PHIODataNode
22 #include <phool/PHNode.h> // for PHNode
23 #include <phool/PHNodeIterator.h> // for PHNodeIterator
24 #include <phool/PHNodeOperation.h>
25 #include <phool/PHObject.h> // for PHObject
26 #include <phool/getClass.h>
27 
28 #include <TObject.h>
29 
30 #include <HepMC/GenEvent.h>
31 
32 #include <climits>
33 #include <iostream>
34 #include <iterator>
35 #include <utility>
36 
37 // convenient aliases for deep copying nodes
38 namespace
39 {
40  using PHG4Particle_t = PHG4Particlev3;
41  using PHG4VtxPoint_t = PHG4VtxPointv1;
42  using PHG4Hit_t = PHG4Hitv1;
43 
45  class FindG4HitContainer : public PHNodeOperation
46  {
47  public:
49  using ContainerMap = std::map<std::string, PHG4HitContainer *>;
50 
52  const ContainerMap &containers() const
53  {
54  return m_containers;
55  }
56 
57  protected:
59  void perform(PHNode *node) override
60  {
61  // check type name. Only load PHIODataNode
62  if (node->getType() != "PHIODataNode") return;
63 
64  // cast to IODataNode and check data
65  auto ionode = static_cast<PHIODataNode<TObject> *>(node);
66  auto data = dynamic_cast<PHG4HitContainer *>(ionode->getData());
67  if (data)
68  {
69  m_containers.insert(std::make_pair(node->getName(), data));
70  }
71  }
72 
73  private:
75  ContainerMap m_containers;
76  };
77 
78 } // namespace
79 
80 //_____________________________________________________________________________
82 {
83  // hep mc
84  m_geneventmap = findNode::getClass<PHHepMCGenEventMap>(dstNode, "PHHepMCGenEventMap");
85  if (!m_geneventmap)
86  {
87  std::cout << "Fun4AllDstPileupMerger::load_nodes - creating PHHepMCGenEventMap" << std::endl;
89  dstNode->addNode(new PHIODataNode<PHObject>(m_geneventmap, "PHHepMCGenEventMap", "PHObject"));
90  }
91 
92  // find all G4Hit containers under dstNode
93  FindG4HitContainer nodeFinder;
94  PHNodeIterator(dstNode).forEach(nodeFinder);
95  m_g4hitscontainers = nodeFinder.containers();
96 
97  // g4 truth info
98  m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(dstNode, "G4TruthInfo");
99  if (!m_g4truthinfo)
100  {
101  std::cout << "Fun4AllDstPileupMerger::load_nodes - creating node G4TruthInfo" << std::endl;
103  dstNode->addNode(new PHIODataNode<PHObject>(m_g4truthinfo, "G4TruthInfo", "PHObject"));
104  }
105 }
106 
107 //_____________________________________________________________________________
109 {
110  // copy PHHepMCGenEventMap
111  const auto map = findNode::getClass<PHHepMCGenEventMap>(dstNode, "PHHepMCGenEventMap");
112 
113  // keep track of new embed id, after insertion as background event
114  int new_embed_id = -1;
115 
116  if (map && m_geneventmap)
117  {
118  if (map->size() != 1)
119  {
120  std::cout << "Fun4AllDstPileupMerger::copy_background_event - cannot merge events that contain more than one PHHepMCGenEventMap" << std::endl;
121  return;
122  }
123 
124  // get event and insert in new map
125  auto genevent = map->get_map().begin()->second;
127 
128  /*
129  * this hack prevents a crash when writting out
130  * it boils down to root trying to write deleted items from the HepMC::GenEvent copy if the source has been deleted
131  * it does not happen if the source gets written while the copy is deleted
132  */
133  newevent->getEvent()->swap(*genevent->getEvent());
134 
135  // shift vertex time and store new embed id
136  newevent->moveVertex(0, 0, 0, delta_t);
137  new_embed_id = newevent->get_embedding_id();
138  }
139 
140  // copy truth container
141  // keep track of the correspondance between source index and destination index for vertices, tracks and showers
142  using ConversionMap = std::map<int, int>;
143  ConversionMap vtxid_map;
144  ConversionMap trkid_map;
145 
146  const auto container_truth = findNode::getClass<PHG4TruthInfoContainer>(dstNode, "G4TruthInfo");
147  if (container_truth && m_g4truthinfo)
148  {
149  {
150  // primary vertices
151  auto key = m_g4truthinfo->maxvtxindex();
152  const auto range = container_truth->GetPrimaryVtxRange();
153  for (auto iter = range.first; iter != range.second; ++iter)
154  {
155  // clone vertex, insert in map, and add index conversion
156  const auto &sourceVertex = iter->second;
157  auto newVertex = new PHG4VtxPoint_t(sourceVertex);
158  newVertex->set_t(sourceVertex->get_t() + delta_t);
159  m_g4truthinfo->AddVertex(++key, newVertex);
160  vtxid_map.insert(std::make_pair(sourceVertex->get_id(), key));
161  }
162  }
163 
164  {
165  // secondary vertices
166  auto key = m_g4truthinfo->minvtxindex();
167  const auto range = container_truth->GetSecondaryVtxRange();
168 
169  // loop from last to first to preserve order with respect to the original event
170  for (
171  auto iter = std::reverse_iterator<PHG4TruthInfoContainer::ConstVtxIterator>(range.second);
172  iter != std::reverse_iterator<PHG4TruthInfoContainer::ConstVtxIterator>(range.first);
173  ++iter)
174  {
175  // clone vertex, shift time, insert in map, and add index conversion
176  const auto &sourceVertex = iter->second;
177  auto newVertex = new PHG4VtxPoint_t(sourceVertex);
178  newVertex->set_t(sourceVertex->get_t() + delta_t);
179  m_g4truthinfo->AddVertex(--key, newVertex);
180  vtxid_map.insert(std::make_pair(sourceVertex->get_id(), key));
181  }
182  }
183 
184  {
185  // primary particles
186  auto key = m_g4truthinfo->maxtrkindex();
187  const auto range = container_truth->GetPrimaryParticleRange();
188  for (auto iter = range.first; iter != range.second; ++iter)
189  {
190  const auto &source = iter->second;
191  auto dest = new PHG4Particle_t(source);
192  m_g4truthinfo->AddParticle(++key, dest);
193  dest->set_track_id(key);
194 
195  // set parent to zero
196  dest->set_parent_id(0);
197 
198  // set primary to itself
199  dest->set_primary_id(dest->get_track_id());
200 
201  // update vertex
202  const auto keyiter = vtxid_map.find(source->get_vtx_id());
203  if (keyiter != vtxid_map.end())
204  dest->set_vtx_id(keyiter->second);
205  else
206  std::cout << "Fun4AllDstPileupMerger::copy_background_event - vertex id " << source->get_vtx_id() << " not found in map" << std::endl;
207 
208  // insert in map
209  trkid_map.insert(std::make_pair(source->get_track_id(), dest->get_track_id()));
210  }
211  }
212 
213  {
214  // secondary particles
215  auto key = m_g4truthinfo->mintrkindex();
216  const auto range = container_truth->GetSecondaryParticleRange();
217 
218  /*
219  * loop from last to first to preserve order with respect to the original event
220  * also this ensures that for a given particle its parent has already been converted and thus found in the map
221  */
222  for (
223  auto iter = std::reverse_iterator<PHG4TruthInfoContainer::ConstIterator>(range.second);
224  iter != std::reverse_iterator<PHG4TruthInfoContainer::ConstIterator>(range.first);
225  ++iter)
226  {
227  const auto &source = iter->second;
228  auto dest = new PHG4Particle_t(source);
229  m_g4truthinfo->AddParticle(--key, dest);
230  dest->set_track_id(key);
231 
232  // update parent id
233  auto keyiter = trkid_map.find(source->get_parent_id());
234  if (keyiter != trkid_map.end())
235  dest->set_parent_id(keyiter->second);
236  else
237  std::cout << "Fun4AllDstPileupMerger::copy_background_event - track id " << source->get_parent_id() << " not found in map" << std::endl;
238 
239  // update primary id
240  keyiter = trkid_map.find(source->get_primary_id());
241  if (keyiter != trkid_map.end())
242  dest->set_primary_id(keyiter->second);
243  else
244  std::cout << "Fun4AllDstPileupMerger::copy_background_event - track id " << source->get_primary_id() << " not found in map" << std::endl;
245 
246  // update vertex
247  keyiter = vtxid_map.find(source->get_vtx_id());
248  if (keyiter != vtxid_map.end())
249  dest->set_vtx_id(keyiter->second);
250  else
251  std::cout << "Fun4AllDstPileupMerger::copy_background_event - vertex id " << source->get_vtx_id() << " not found in map" << std::endl;
252 
253  // insert in map
254  trkid_map.insert(std::make_pair(source->get_track_id(), dest->get_track_id()));
255  }
256  }
257 
258  // vertex embed flags
259  /* embed flag is stored only for primary vertices, consistently with PHG4TruthEventAction */
260  for (const auto &pair : vtxid_map)
261  {
262  if (pair.first > 0) m_g4truthinfo->AddEmbededVtxId(pair.second, new_embed_id);
263  }
264 
265  // track embed flags
266  /* embed flag is stored only for primary tracks, consistently with PHG4TruthEventAction */
267  for (const auto &pair : trkid_map)
268  {
269  if (pair.first > 0) m_g4truthinfo->AddEmbededTrkId(pair.second, new_embed_id);
270  }
271  }
272 
273  // copy g4hits
274  // loop over registered maps
275  for (const auto &pair : m_g4hitscontainers)
276  {
277  // check destination node
278  if (!pair.second)
279  {
280  std::cout << "Fun4AllDstPileupMerger::copy_background_event - invalid destination container " << pair.first << std::endl;
281  continue;
282  }
283 
284  // find source node
285  auto container_hit = findNode::getClass<PHG4HitContainer>(dstNode, pair.first);
286  if (!container_hit)
287  {
288  std::cout << "Fun4AllDstPileupMerger::copy_background_event - invalid source container " << pair.first << std::endl;
289  continue;
290  }
291 
292  {
293  // hits
294  const auto range = container_hit->getHits();
295  for (auto iter = range.first; iter != range.second; ++iter)
296  {
297  // clone hit
298  const auto &sourceHit = iter->second;
299  auto newHit = new PHG4Hit_t(sourceHit);
300 
301  // shift time
302  newHit->set_t(0, sourceHit->get_t(0) + delta_t);
303  newHit->set_t(1, sourceHit->get_t(1) + delta_t);
304 
305  // update track id
306  const auto keyiter = trkid_map.find(sourceHit->get_trkid());
307  if (keyiter != trkid_map.end())
308  newHit->set_trkid(keyiter->second);
309  else
310  std::cout << "Fun4AllDstPileupMerger::copy_background_event - track id " << sourceHit->get_trkid() << " not found in map" << std::endl;
311 
312  /*
313  * reset shower ids
314  * it was decided that showers from the background events will not be copied to the merged event
315  * as such we just reset the hits shower id
316  */
317  newHit->set_shower_id(INT_MIN);
318 
319  /*
320  * this will generate a new key for the hit and assign it to the hit
321  * this ensures that there is no conflict with the hits from the 'main' event
322  */
323  pair.second->AddHit(newHit->get_detid(), newHit);
324  }
325  }
326 
327  {
328  // layers
329  const auto range = container_hit->getLayers();
330  for (auto iter = range.first; iter != range.second; ++iter)
331  {
332  pair.second->AddLayer(*iter);
333  }
334  }
335  }
336 }