Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EventEvaluator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EventEvaluator.cc
1 #include "EventEvaluator.h"
2 
3 #include "CaloEvalStack.h"
4 #include "CaloRawClusterEval.h"
5 #include "CaloRawTowerEval.h"
6 #include "CaloTruthEval.h"
7 
8 #include <g4main/PHG4Hit.h>
9 #include <g4main/PHG4HitContainer.h>
10 #include <g4main/PHG4Particle.h>
11 #include <g4main/PHG4Shower.h>
12 #include <g4main/PHG4TruthInfoContainer.h>
13 #include <g4main/PHG4VtxPoint.h>
14 
15 #include <globalvertex/GlobalVertex.h>
16 #include <globalvertex/GlobalVertexMap.h>
17 
18 #include <trackbase_historic/SvtxTrack.h>
19 #include <trackbase_historic/SvtxTrackMap.h>
20 #include <trackbase_historic/SvtxTrack_FastSim.h>
21 #include <globalvertex/SvtxVertex.h> // for SvtxVertex
22 #include <globalvertex/SvtxVertexMap.h>
23 
24 #include <calobase/RawCluster.h>
25 #include <calobase/RawClusterContainer.h>
26 #include <calobase/RawClusterUtility.h>
27 #include <calobase/RawTower.h>
28 #include <calobase/RawTowerContainer.h>
29 #include <calobase/RawTowerGeom.h>
30 #include <calobase/RawTowerGeomContainer.h>
31 #include <calobase/RawTowerv2.h>
32 
33 #include <phhepmc/PHGenIntegral.h>
34 #include <phhepmc/PHHepMCGenEvent.h>
35 #include <phhepmc/PHHepMCGenEventMap.h>
36 
37 #include <fun4all/Fun4AllReturnCodes.h>
38 #include <fun4all/SubsysReco.h>
39 
40 #include <phool/PHCompositeNode.h>
41 #include <phool/PHNodeIterator.h> // for PHNodeIterator
42 #include <phool/getClass.h>
43 #include <phool/phool.h>
44 
45 #include <TFile.h>
46 #include <TNtuple.h>
47 #include <TTree.h>
48 
49 #include <CLHEP/Vector/ThreeVector.h>
50 
51 #pragma GCC diagnostic push
52 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
53 #include <HepMC/GenEvent.h>
54 #include <HepMC/GenVertex.h>
55 #pragma GCC diagnostic pop
56 
57 #include <cmath>
58 #include <cstdlib>
59 #include <iostream>
60 #include <set>
61 #include <utility>
62 
63 EventEvaluator::EventEvaluator(const std::string& name, const std::string& filename)
64  : SubsysReco(name)
65  , _filename(filename)
66 {
67  _hits_layerID = new int[_maxNHits];
68  _hits_trueID = new int[_maxNHits];
69  _hits_x = new float[_maxNHits];
70  _hits_y = new float[_maxNHits];
71  _hits_z = new float[_maxNHits];
72  _hits_t = new float[_maxNHits];
73 
74  _tower_HCALIN_E = new float[_maxNTowersCentral];
75  _tower_HCALIN_iEta = new int[_maxNTowersCentral];
76  _tower_HCALIN_iPhi = new int[_maxNTowersCentral];
77  _tower_HCALIN_trueID = new int[_maxNTowersCentral];
78  _cluster_HCALIN_E = new float[_maxNclustersCentral];
79  _cluster_HCALIN_Eta = new float[_maxNclustersCentral];
80  _cluster_HCALIN_Phi = new float[_maxNclustersCentral];
81  _cluster_HCALIN_NTower = new int[_maxNclustersCentral];
82  _cluster_HCALIN_trueID = new int[_maxNclustersCentral];
83 
84  _tower_HCALOUT_E = new float[_maxNTowersCentral];
85  _tower_HCALOUT_iEta = new int[_maxNTowersCentral];
86  _tower_HCALOUT_iPhi = new int[_maxNTowersCentral];
87  _tower_HCALOUT_trueID = new int[_maxNTowersCentral];
88  _cluster_HCALOUT_E = new float[_maxNclustersCentral];
89  _cluster_HCALOUT_Eta = new float[_maxNclustersCentral];
90  _cluster_HCALOUT_Phi = new float[_maxNclustersCentral];
91  _cluster_HCALOUT_NTower = new int[_maxNclustersCentral];
92  _cluster_HCALOUT_trueID = new int[_maxNclustersCentral];
93 
94  _tower_CEMC_E = new float[_maxNTowersCentral];
95  _tower_CEMC_iEta = new int[_maxNTowersCentral];
96  _tower_CEMC_iPhi = new int[_maxNTowersCentral];
97  _tower_CEMC_trueID = new int[_maxNTowersCentral];
98  _cluster_CEMC_E = new float[_maxNclustersCentral];
99  _cluster_CEMC_Eta = new float[_maxNclustersCentral];
100  _cluster_CEMC_Phi = new float[_maxNclustersCentral];
101  _cluster_CEMC_NTower = new int[_maxNclustersCentral];
102  _cluster_CEMC_trueID = new int[_maxNclustersCentral];
103 
104  _track_ID = new float[_maxNTracks];
105  _track_trueID = new float[_maxNTracks];
106  _track_px = new float[_maxNTracks];
107  _track_py = new float[_maxNTracks];
108  _track_pz = new float[_maxNTracks];
109  _track_dca = new float[_maxNTracks];
110  _track_dca_2d = new float[_maxNTracks];
111  _track_source = new unsigned short[_maxNTracks];
112  _track_ProjTrackID = new float[_maxNProjections];
113  _track_ProjLayer = new int[_maxNProjections];
114  _track_TLP_x = new float[_maxNProjections];
115  _track_TLP_y = new float[_maxNProjections];
116  _track_TLP_z = new float[_maxNProjections];
117  _track_TLP_t = new float[_maxNProjections];
118  _track_TLP_true_x = new float[_maxNProjections];
119  _track_TLP_true_y = new float[_maxNProjections];
120  _track_TLP_true_z = new float[_maxNProjections];
121  _track_TLP_true_t = new float[_maxNProjections];
122 
123  _mcpart_ID = new int[_maxNMCPart];
124  _mcpart_ID_parent = new int[_maxNMCPart];
125  _mcpart_PDG = new int[_maxNMCPart];
126  _mcpart_E = new float[_maxNMCPart];
127  _mcpart_px = new float[_maxNMCPart];
128  _mcpart_py = new float[_maxNMCPart];
129  _mcpart_pz = new float[_maxNMCPart];
130  _mcpart_BCID = new int[_maxNMCPart];
131 
132  _hepmcp_BCID = new int[_maxNHepmcp];
133  // _hepmcp_ID_parent = new float[_maxNHepmcp];
134  _hepmcp_status = new int[_maxNHepmcp];
135  _hepmcp_PDG = new int[_maxNHepmcp];
136  _hepmcp_E = new float[_maxNHepmcp];
137  _hepmcp_px = new float[_maxNHepmcp];
138  _hepmcp_py = new float[_maxNHepmcp];
139  _hepmcp_pz = new float[_maxNHepmcp];
140  _hepmcp_m1 = new int[_maxNHepmcp];
141  _hepmcp_m2 = new int[_maxNHepmcp];
142 
143  _calo_towers_iEta = new int[_maxNTowersCalo];
144  _calo_towers_iPhi = new int[_maxNTowersCalo];
145  _calo_towers_Eta = new float[_maxNTowersCalo];
146  _calo_towers_Phi = new float[_maxNTowersCalo];
147  _calo_towers_x = new float[_maxNTowersCalo];
148  _calo_towers_y = new float[_maxNTowersCalo];
149  _calo_towers_z = new float[_maxNTowersCalo];
150  _geometry_done = new int[20];
151  for (int igem = 0; igem < 20; igem++)
152  {
153  _geometry_done[igem] = 0;
154  }
155 }
156 
157 int EventEvaluator::Init(PHCompositeNode* /*topNode*/)
158 {
159  _ievent = 0;
160 
161  _tfile = new TFile(_filename.c_str(), "RECREATE");
162 
163  _event_tree = new TTree("event_tree", "event_tree");
164  if (_do_store_event_info)
165  {
166  // Event level info. This isn't the most efficient way to store this info, but it's straightforward
167  // within the structure of the class, so the size is small compared to the rest of the output.
168  _event_tree->Branch("cross_section", &_cross_section, "cross_section/F");
169  _event_tree->Branch("event_weight", &_event_weight, "event_weight/F");
170  _event_tree->Branch("n_generator_accepted", &_n_generator_accepted, "n_generator_accepted/I");
171  }
172  // tracks and hits
173  if (_do_HITS)
174  {
175  _event_tree->Branch("nHits", &_nHitsLayers, "nHits/I");
176  _event_tree->Branch("hits_layerID", _hits_layerID, "hits_layerID[nHits]/I");
177  _event_tree->Branch("hits_trueID", _hits_trueID, "hits_trueID[nHits]/I");
178  _event_tree->Branch("hits_x", _hits_x, "hits_x[nHits]/F");
179  _event_tree->Branch("hits_y", _hits_y, "hits_y[nHits]/F");
180  _event_tree->Branch("hits_z", _hits_z, "hits_z[nHits]/F");
181  _event_tree->Branch("hits_t", _hits_t, "hits_t[nHits]/F");
182  }
183  if (_do_TRACKS)
184  {
185  _event_tree->Branch("nTracks", &_nTracks, "nTracks/I");
186  _event_tree->Branch("tracks_ID", _track_ID, "tracks_ID[nTracks]/F");
187  _event_tree->Branch("tracks_px", _track_px, "tracks_px[nTracks]/F");
188  _event_tree->Branch("tracks_py", _track_py, "tracks_py[nTracks]/F");
189  _event_tree->Branch("tracks_pz", _track_pz, "tracks_pz[nTracks]/F");
190  _event_tree->Branch("tracks_dca", _track_dca, "tracks_dca[nTracks]/F");
191  _event_tree->Branch("tracks_dca_2d", _track_dca_2d, "tracks_dca_2d[nTracks]/F");
192  _event_tree->Branch("tracks_trueID", _track_trueID, "tracks_trueID[nTracks]/F");
193  _event_tree->Branch("tracks_source", _track_source, "tracks_source[nTracks]/s");
194  }
195  if (_do_PROJECTIONS)
196  {
197  _event_tree->Branch("nProjections", &_nProjections, "nProjections/I");
198  _event_tree->Branch("track_ProjTrackID", _track_ProjTrackID, "track_ProjTrackID[nProjections]/F");
199  _event_tree->Branch("track_ProjLayer", _track_ProjLayer, "track_ProjLayer[nProjections]/I");
200  _event_tree->Branch("track_TLP_x", _track_TLP_x, "track_TLP_x[nProjections]/F");
201  _event_tree->Branch("track_TLP_y", _track_TLP_y, "track_TLP_y[nProjections]/F");
202  _event_tree->Branch("track_TLP_z", _track_TLP_z, "track_TLP_z[nProjections]/F");
203  _event_tree->Branch("track_TLP_t", _track_TLP_t, "track_TLP_t[nProjections]/F");
204  _event_tree->Branch("track_TLP_true_x", _track_TLP_true_x, "track_TLP_true_x[nProjections]/F");
205  _event_tree->Branch("track_TLP_true_y", _track_TLP_true_y, "track_TLP_true_y[nProjections]/F");
206  _event_tree->Branch("track_TLP_true_z", _track_TLP_true_z, "track_TLP_true_z[nProjections]/F");
207  _event_tree->Branch("track_TLP_true_t", _track_TLP_true_t, "track_TLP_true_t[nProjections]/F");
208  }
209  if (_do_HCALIN)
210  {
211  // towers HCAL-in
212  _event_tree->Branch("tower_HCALIN_N", &_nTowers_HCALIN, "tower_HCALIN_N/I");
213  _event_tree->Branch("tower_HCALIN_E", _tower_HCALIN_E, "tower_HCALIN_E[tower_HCALIN_N]/F");
214  _event_tree->Branch("tower_HCALIN_iEta", _tower_HCALIN_iEta, "tower_HCALIN_iEta[tower_HCALIN_N]/I");
215  _event_tree->Branch("tower_HCALIN_iPhi", _tower_HCALIN_iPhi, "tower_HCALIN_iPhi[tower_HCALIN_N]/I");
216  _event_tree->Branch("tower_HCALIN_trueID", _tower_HCALIN_trueID, "tower_HCALIN_trueID[tower_HCALIN_N]/I");
217  if (_do_CLUSTERS)
218  {
219  // clusters HCAL-in
220  _event_tree->Branch("cluster_HCALIN_N", &_nclusters_HCALIN, "cluster_HCALIN_N/I");
221  _event_tree->Branch("cluster_HCALIN_E", _cluster_HCALIN_E, "cluster_HCALIN_E[cluster_HCALIN_N]/F");
222  _event_tree->Branch("cluster_HCALIN_Eta", _cluster_HCALIN_Eta, "cluster_HCALIN_Eta[cluster_HCALIN_N]/F");
223  _event_tree->Branch("cluster_HCALIN_Phi", _cluster_HCALIN_Phi, "cluster_HCALIN_Phi[cluster_HCALIN_N]/F");
224  _event_tree->Branch("cluster_HCALIN_NTower", _cluster_HCALIN_NTower, "cluster_HCALIN_NTower[cluster_HCALIN_N]/I");
225  _event_tree->Branch("cluster_HCALIN_trueID", _cluster_HCALIN_trueID, "cluster_HCALIN_trueID[cluster_HCALIN_N]/I");
226  }
227  }
228  if (_do_HCALOUT)
229  {
230  // towers HCAL-out
231  _event_tree->Branch("tower_HCALOUT_N", &_nTowers_HCALOUT, "tower_HCALOUT_N/I");
232  _event_tree->Branch("tower_HCALOUT_E", _tower_HCALOUT_E, "tower_HCALOUT_E[tower_HCALOUT_N]/F");
233  _event_tree->Branch("tower_HCALOUT_iEta", _tower_HCALOUT_iEta, "tower_HCALOUT_iEta[tower_HCALOUT_N]/I");
234  _event_tree->Branch("tower_HCALOUT_iPhi", _tower_HCALOUT_iPhi, "tower_HCALOUT_iPhi[tower_HCALOUT_N]/I");
235  _event_tree->Branch("tower_HCALOUT_trueID", _tower_HCALOUT_trueID, "tower_HCALOUT_trueID[tower_HCALOUT_N]/I");
236  if (_do_CLUSTERS)
237  {
238  // clusters HCAL-out
239  _event_tree->Branch("cluster_HCALOUT_N", &_nclusters_HCALOUT, "cluster_HCALOUT_N/I");
240  _event_tree->Branch("cluster_HCALOUT_E", _cluster_HCALOUT_E, "cluster_HCALOUT_E[cluster_HCALOUT_N]/F");
241  _event_tree->Branch("cluster_HCALOUT_Eta", _cluster_HCALOUT_Eta, "cluster_HCALOUT_Eta[cluster_HCALOUT_N]/F");
242  _event_tree->Branch("cluster_HCALOUT_Phi", _cluster_HCALOUT_Phi, "cluster_HCALOUT_Phi[cluster_HCALOUT_N]/F");
243  _event_tree->Branch("cluster_HCALOUT_NTower", _cluster_HCALOUT_NTower, "cluster_HCALOUT_NTower[cluster_HCALOUT_N]/I");
244  _event_tree->Branch("cluster_HCALOUT_trueID", _cluster_HCALOUT_trueID, "cluster_HCALOUT_trueID[cluster_HCALOUT_N]/I");
245  }
246  }
247  if (_do_CEMC)
248  {
249  // towers CEMC
250  _event_tree->Branch("tower_CEMC_N", &_nTowers_CEMC, "tower_CEMC_N/I");
251  _event_tree->Branch("tower_CEMC_E", _tower_CEMC_E, "tower_CEMC_E[tower_CEMC_N]/F");
252  _event_tree->Branch("tower_CEMC_iEta", _tower_CEMC_iEta, "tower_CEMC_iEta[tower_CEMC_N]/I");
253  _event_tree->Branch("tower_CEMC_iPhi", _tower_CEMC_iPhi, "tower_CEMC_iPhi[tower_CEMC_N]/I");
254  _event_tree->Branch("tower_CEMC_trueID", _tower_CEMC_trueID, "tower_CEMC_trueID[tower_CEMC_N]/I");
255  if (_do_CLUSTERS)
256  {
257  // clusters CEMC
258  _event_tree->Branch("cluster_CEMC_N", &_nclusters_CEMC, "cluster_CEMC_N/I");
259  _event_tree->Branch("cluster_CEMC_E", _cluster_CEMC_E, "cluster_CEMC_E[cluster_CEMC_N]/F");
260  _event_tree->Branch("cluster_CEMC_Eta", _cluster_CEMC_Eta, "cluster_CEMC_Eta[cluster_CEMC_N]/F");
261  _event_tree->Branch("cluster_CEMC_Phi", _cluster_CEMC_Phi, "cluster_CEMC_Phi[cluster_CEMC_N]/F");
262  _event_tree->Branch("cluster_CEMC_NTower", _cluster_CEMC_NTower, "cluster_CEMC_NTower[cluster_CEMC_N]/I");
263  _event_tree->Branch("cluster_CEMC_trueID", _cluster_CEMC_trueID, "cluster_CEMC_trueID[cluster_CEMC_N]/I");
264  }
265  }
266  if (_do_VERTEX)
267  {
268  // vertex
269  _event_tree->Branch("vertex_x", &_vertex_x, "vertex_x/F");
270  _event_tree->Branch("vertex_y", &_vertex_y, "vertex_y/F");
271  _event_tree->Branch("vertex_z", &_vertex_z, "vertex_z/F");
272  _event_tree->Branch("vertex_NCont", &_vertex_NCont, "vertex_NCont/I");
273  _event_tree->Branch("vertex_true_x", &_vertex_true_x, "vertex_true_x/F");
274  _event_tree->Branch("vertex_true_y", &_vertex_true_y, "vertex_true_y/F");
275  _event_tree->Branch("vertex_true_z", &_vertex_true_z, "vertex_true_z/F");
276  }
277  if (_do_MCPARTICLES)
278  {
279  // MC particles
280  _event_tree->Branch("nMCPart", &_nMCPart, "nMCPart/I");
281  _event_tree->Branch("mcpart_ID", _mcpart_ID, "mcpart_ID[nMCPart]/I");
282  _event_tree->Branch("mcpart_ID_parent", _mcpart_ID_parent, "mcpart_ID_parent[nMCPart]/I");
283  _event_tree->Branch("mcpart_PDG", _mcpart_PDG, "mcpart_PDG[nMCPart]/I");
284  _event_tree->Branch("mcpart_E", _mcpart_E, "mcpart_E[nMCPart]/F");
285  _event_tree->Branch("mcpart_px", _mcpart_px, "mcpart_px[nMCPart]/F");
286  _event_tree->Branch("mcpart_py", _mcpart_py, "mcpart_py[nMCPart]/F");
287  _event_tree->Branch("mcpart_pz", _mcpart_pz, "mcpart_pz[nMCPart]/F");
288  _event_tree->Branch("mcpart_BCID", _mcpart_BCID, "mcpart_BCID[nMCPart]/I");
289  }
290  if (_do_HEPMC)
291  {
292  // MC particles
293  _event_tree->Branch("nHepmcp", &_nHepmcp, "nHepmcp/I");
294  _event_tree->Branch("hepmcp_procid", &_hepmcp_procid, "hepmcp_procid/I");
295  _event_tree->Branch("hepmcp_x1", &_hepmcp_x1, "hepmcp_x1/F");
296  _event_tree->Branch("hepmcp_x2", &_hepmcp_x2, "hepmcp_x2/F");
297 
298  // _event_tree->Branch("hepmcp_ID_parent", _hepmcp_ID_parent, "hepmcp_ID_parent[nHepmcp]/F");
299  _event_tree->Branch("hepmcp_status", _hepmcp_status, "hepmcp_status[nHepmcp]/I");
300  _event_tree->Branch("hepmcp_PDG", _hepmcp_PDG, "hepmcp_PDG[nHepmcp]/I");
301  _event_tree->Branch("hepmcp_E", _hepmcp_E, "hepmcp_E[nHepmcp]/F");
302  _event_tree->Branch("hepmcp_px", _hepmcp_px, "hepmcp_px[nHepmcp]/F");
303  _event_tree->Branch("hepmcp_py", _hepmcp_py, "hepmcp_py[nHepmcp]/F");
304  _event_tree->Branch("hepmcp_pz", _hepmcp_pz, "hepmcp_pz[nHepmcp]/F");
305  _event_tree->Branch("hepmcp_BCID", _hepmcp_BCID, "hepmcp_BCID[nHepmcp]/I");
306  _event_tree->Branch("hepmcp_m1", _hepmcp_m1, "hepmcp_m1[nHepmcp]/I");
307  _event_tree->Branch("hepmcp_m2", _hepmcp_m2, "hepmcp_m2[nHepmcp]/I");
308  }
309 
310  if (_do_GEOMETRY)
311  {
312  _tfile_geometry = new TFile("geometry.root", "RECREATE");
313 
314  _geometry_tree = new TTree("geometry_tree", "geometry_tree");
315  // tracks and hits
316  _geometry_tree->Branch("calo", &_calo_ID, "nHits/I");
317  _geometry_tree->Branch("calo_towers_N", &_calo_towers_N, "calo_towers_N/I");
318  _geometry_tree->Branch("calo_towers_iEta", _calo_towers_iEta, "calo_towers_iEta[calo_towers_N]/I");
319  _geometry_tree->Branch("calo_towers_iPhi", _calo_towers_iPhi, "calo_towers_iPhi[calo_towers_N]/I");
320  _geometry_tree->Branch("calo_towers_Eta", _calo_towers_Eta, "calo_towers_Eta[calo_towers_N]/F");
321  _geometry_tree->Branch("calo_towers_Phi", _calo_towers_Phi, "calo_towers_Phi[calo_towers_N]/F");
322  _geometry_tree->Branch("calo_towers_x", _calo_towers_x, "calo_towers_x[calo_towers_N]/F");
323  _geometry_tree->Branch("calo_towers_y", _calo_towers_y, "calo_towers_y[calo_towers_N]/F");
324  _geometry_tree->Branch("calo_towers_z", _calo_towers_z, "calo_towers_z[calo_towers_N]/F");
325  }
326 
327  return Fun4AllReturnCodes::EVENT_OK;
328 }
329 
330 int EventEvaluator::process_event(PHCompositeNode* topNode)
331 {
332  if (Verbosity() > 0)
333  {
334  std::cout << "entered process_event" << std::endl;
335  }
336  if (_do_HCALIN)
337  {
338  if (!_caloevalstackHCALIN)
339  {
340  _caloevalstackHCALIN = new CaloEvalStack(topNode, "HCALIN");
341  _caloevalstackHCALIN->set_strict(_strict);
342  _caloevalstackHCALIN->set_verbosity(Verbosity() + 1);
343  }
344  else
345  {
346  _caloevalstackHCALIN->next_event(topNode);
347  }
348  }
349  if (_do_HCALOUT)
350  {
351  if (!_caloevalstackHCALOUT)
352  {
353  _caloevalstackHCALOUT = new CaloEvalStack(topNode, "HCALOUT");
354  _caloevalstackHCALOUT->set_strict(_strict);
355  _caloevalstackHCALOUT->set_verbosity(Verbosity() + 1);
356  }
357  else
358  {
359  _caloevalstackHCALOUT->next_event(topNode);
360  }
361  }
362  if (_do_CEMC)
363  {
364  if (!_caloevalstackCEMC)
365  {
366  _caloevalstackCEMC = new CaloEvalStack(topNode, "CEMC");
367  _caloevalstackCEMC->set_strict(_strict);
368  _caloevalstackCEMC->set_verbosity(Verbosity() + 1);
369  }
370  else
371  {
372  _caloevalstackCEMC->next_event(topNode);
373  }
374  }
375 
376  if (Verbosity() > 0)
377  {
378  std::cout << "loaded evalstack" << std::endl;
379  }
380 
381  // fill the Evaluator Tree
382  fillOutputNtuples(topNode);
383 
384  ++_ievent;
385 
386  return Fun4AllReturnCodes::EVENT_OK;
387 }
388 
389 void EventEvaluator::fillOutputNtuples(PHCompositeNode* topNode)
390 {
391  if (Verbosity() > 2)
392  {
393  std::cout << "EventEvaluator::fillOutputNtuples() entered" << std::endl;
394  }
395 
396  //----------------------
397  // fill the Event Tree
398  //----------------------
399 
400  //----------------------
401  // Event level info
402  //---------------------
403  // Extract weight info from the stored HepMC event.
404  if (_do_store_event_info)
405  {
406  PHHepMCGenEventMap* hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
407  if (hepmceventmap)
408  {
409  if (Verbosity() > 0)
410  {
411  std::cout << "saving event level info" << std::endl;
412  }
413 
414  for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
415  eventIter != hepmceventmap->end();
416  ++eventIter)
417  {
418  PHHepMCGenEvent* hepmcevent = eventIter->second;
419 
420  if (hepmcevent)
421  {
422  HepMC::GenEvent* truthevent = hepmcevent->getEvent();
423  if (!truthevent)
424  {
425  std::cout << PHWHERE
426  << "no evt pointer under phhepmvgeneventmap found "
427  << std::endl;
428  return;
429  }
430 
431  auto xsec = truthevent->cross_section();
432  if (xsec)
433  {
434  _cross_section = xsec->cross_section();
435  }
436  // Only fill the event weight if available.
437  // The overall event weight will be stored in the last entry in the vector.
438  auto weights = truthevent->weights();
439  if (weights.size() > 0)
440  {
441  _event_weight = weights[weights.size() - 1];
442  }
443  }
444  }
445  }
446  else
447  {
448  if (Verbosity() > 0)
449  {
450  std::cout << PHWHERE << " PHHepMCGenEventMap node (for event level info) not found on node tree" << std::endl;
451  }
452  return;
453  }
454 
455  // Retrieve the number of generator accepted events
456  // Following how this was implemented in PHPythia8
457  PHNodeIterator iter(topNode);
458  PHCompositeNode* sumNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
459  if (!sumNode)
460  {
461  std::cout << PHWHERE << "RUN Node missing doing nothing" << std::endl;
462  return;
463  }
464  auto* integralNode = findNode::getClass<PHGenIntegral>(sumNode, "PHGenIntegral");
465  if (integralNode)
466  {
467  _n_generator_accepted = integralNode->get_N_Generator_Accepted_Event();
468  }
469  else
470  {
471  if (Verbosity() > 0)
472  {
473  std::cout << PHWHERE << " PHGenIntegral node (for n generator accepted) not found on node tree. Continuing" << std::endl;
474  }
475  }
476  }
477  //----------------------
478  // VERTEX
479  //----------------------
480  SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
481  if (_do_VERTEX)
482  {
483  if (vertexmap)
484  {
485  if (!vertexmap->empty())
486  {
487  if (Verbosity() > 0)
488  {
489  std::cout << "saving vertex" << std::endl;
490  }
491  SvtxVertex* vertex = (vertexmap->begin())->second;
492 
493  _vertex_x = vertex->get_x();
494  _vertex_y = vertex->get_y();
495  _vertex_z = vertex->get_z();
496  _vertex_NCont = vertex->size_tracks();
497  }
498  else
499  {
500  _vertex_x = 0.;
501  _vertex_y = 0.;
502  _vertex_z = 0.;
503  _vertex_NCont = -1;
504  }
505  }
506  }
507  //----------------------
508  // HITS
509  //----------------------
510  if (_do_HITS)
511  {
512  if (Verbosity() > 0)
513  {
514  std::cout << "saving hits" << std::endl;
515  }
516  _nHitsLayers = 0;
517  PHG4TruthInfoContainer* truthinfocontainerHits = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
518  for (int iIndex = 0; iIndex < 60; ++iIndex)
519  {
520  // you need to add your layer name here to be saved! This has to be done
521  // as we do not want to save thousands of calorimeter hits!
522  if ((GetProjectionNameFromIndex(iIndex).find("MVTX") != std::string::npos) ||
523  (GetProjectionNameFromIndex(iIndex).find("INTT") != std::string::npos))
524  {
525  std::string nodename = "G4HIT_" + GetProjectionNameFromIndex(iIndex);
526  PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
527  if (hits)
528  {
529  if (Verbosity() > 1)
530  {
531  std::cout << __PRETTY_FUNCTION__ << " number of hits: " << hits->size() << std::endl;
532  }
533  PHG4HitContainer::ConstRange hit_range = hits->getHits();
534  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
535  {
536  if (Verbosity() > 1)
537  {
538  std::cout << __PRETTY_FUNCTION__ << " found hit with id " << hit_iter->second->get_trkid() << std::endl;
539  }
540  if (_nHitsLayers > _maxNHits)
541  {
542  std::cout << __PRETTY_FUNCTION__ << " exceededed maximum hit array size! Please check where these hits come from!" << std::endl;
543  break;
544  }
545  _hits_x[_nHitsLayers] = hit_iter->second->get_x(0);
546  _hits_y[_nHitsLayers] = hit_iter->second->get_y(0);
547  _hits_z[_nHitsLayers] = hit_iter->second->get_z(0);
548  _hits_t[_nHitsLayers] = hit_iter->second->get_t(0);
549  _hits_layerID[_nHitsLayers] = iIndex;
550  if (truthinfocontainerHits)
551  {
552  PHG4Particle* particle = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
553 
554  if (particle->get_parent_id() != 0)
555  {
556  PHG4Particle* g4particleMother = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
557  int mcSteps = 0;
558  while (g4particleMother->get_parent_id() != 0)
559  {
560  g4particleMother = truthinfocontainerHits->GetParticle(g4particleMother->get_parent_id());
561  if (g4particleMother == nullptr)
562  {
563  break;
564  }
565  mcSteps += 1;
566  }
567  if (mcSteps <= _depth_MCstack)
568  {
569  _hits_trueID[_nHitsLayers] = hit_iter->second->get_trkid();
570  }
571  else
572  {
573  PHG4Particle* g4particleMother2 = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
574  int mcSteps2 = 0;
575  while (g4particleMother2->get_parent_id() != 0 && (mcSteps2 < (mcSteps - _depth_MCstack + 1)))
576  {
577  g4particleMother2 = truthinfocontainerHits->GetParticle(g4particleMother2->get_parent_id());
578  if (g4particleMother2 == nullptr)
579  {
580  break;
581  }
582  else
583  {
584  _hits_trueID[_nHitsLayers] = g4particleMother2->get_parent_id();
585  mcSteps2 += 1;
586  }
587  }
588  }
589  }
590  else
591  {
592  _hits_trueID[_nHitsLayers] = hit_iter->second->get_trkid();
593  }
594  }
595  _nHitsLayers++;
596  }
597  if (Verbosity() > 0)
598  {
599  std::cout << "saved\t" << _nHitsLayers << "\thits for " << GetProjectionNameFromIndex(iIndex) << std::endl;
600  }
601  }
602  else
603  {
604  if (Verbosity() > 0)
605  {
606  std::cout << __PRETTY_FUNCTION__ << " could not find " << nodename << std::endl;
607  }
608  continue;
609  }
610  }
611  }
612  }
613  //----------------------
614  // TOWERS HCALIN
615  //----------------------
616  if (_do_HCALIN)
617  {
618  CaloRawTowerEval* towerevalHCALIN = _caloevalstackHCALIN->get_rawtower_eval();
619  _nTowers_HCALIN = 0;
620  std::string towernodeHCALIN = "TOWER_CALIB_HCALIN";
621  RawTowerContainer* towersHCALIN = findNode::getClass<RawTowerContainer>(topNode, towernodeHCALIN.c_str());
622  if (towersHCALIN)
623  {
624  if (Verbosity() > 0)
625  {
626  std::cout << "saving HCAL towers" << std::endl;
627  }
628  std::string towergeomnodeHCALIN = "TOWERGEOM_HCALIN";
629  RawTowerGeomContainer* towergeomHCALIN = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeHCALIN.c_str());
630  if (towergeomHCALIN)
631  {
632  if (_do_GEOMETRY && !_geometry_done[kHCALIN])
633  {
634  RawTowerGeomContainer::ConstRange all_towers = towergeomHCALIN->get_tower_geometries();
635  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
636  it != all_towers.second; ++it)
637  {
638  _calo_ID = kHCALIN;
639  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
640  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
641  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
642  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
643  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
644  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
645  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
646  _calo_towers_N++;
647  }
648  _geometry_done[kHCALIN] = 1;
649  _geometry_tree->Fill();
650  resetGeometryArrays();
651  }
652  RawTowerContainer::ConstRange begin_end = towersHCALIN->getTowers();
653  RawTowerContainer::ConstIterator rtiter;
654  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
655  {
656  RawTower* tower = rtiter->second;
657  if (tower)
658  {
659  // min energy cut
660  if (tower->get_energy() < _reco_e_threshold)
661  {
662  continue;
663  }
664  _tower_HCALIN_iEta[_nTowers_HCALIN] = tower->get_bineta();
665  _tower_HCALIN_iPhi[_nTowers_HCALIN] = tower->get_binphi();
666  _tower_HCALIN_E[_nTowers_HCALIN] = tower->get_energy();
667 
668  PHG4Particle* primary = towerevalHCALIN->max_truth_primary_particle_by_energy(tower);
669  if (primary)
670  {
671  _tower_HCALIN_trueID[_nTowers_HCALIN] = primary->get_track_id();
672  // gflavor = primary->get_pid();
673  // efromtruth = towerevalHCALIN->get_energy_contribution(tower, primary);
674  }
675  else
676  {
677  _tower_HCALIN_trueID[_nTowers_HCALIN] = -10;
678  }
679  _nTowers_HCALIN++;
680  }
681  }
682  }
683  else
684  {
685  if (Verbosity() > 0)
686  {
687  std::cout << PHWHERE << " ERROR: Can't find " << towergeomnodeHCALIN << std::endl;
688  }
689  // return;
690  }
691  if (Verbosity() > 0)
692  {
693  std::cout << "saved\t" << _nTowers_HCALIN << "\tHCALIN towers" << std::endl;
694  }
695  }
696  else
697  {
698  if (Verbosity() > 0)
699  {
700  std::cout << PHWHERE << " ERROR: Can't find " << towernodeHCALIN << std::endl;
701  }
702  // return;
703  }
704  }
705  //----------------------
706  // TOWERS HCALOUT
707  //----------------------
708  if (_do_HCALOUT)
709  {
710  CaloRawTowerEval* towerevalHCALOUT = _caloevalstackHCALOUT->get_rawtower_eval();
711  _nTowers_HCALOUT = 0;
712  std::string towernodeHCALOUT = "TOWER_CALIB_HCALOUT";
713  RawTowerContainer* towersHCALOUT = findNode::getClass<RawTowerContainer>(topNode, towernodeHCALOUT.c_str());
714  if (towersHCALOUT)
715  {
716  if (Verbosity() > 0)
717  {
718  std::cout << "saving HCAL towers" << std::endl;
719  }
720  std::string towergeomnodeHCALOUT = "TOWERGEOM_HCALOUT";
721  RawTowerGeomContainer* towergeomHCALOUT = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeHCALOUT.c_str());
722  if (towergeomHCALOUT)
723  {
724  if (_do_GEOMETRY && !_geometry_done[kHCALOUT])
725  {
726  RawTowerGeomContainer::ConstRange all_towers = towergeomHCALOUT->get_tower_geometries();
727  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
728  it != all_towers.second; ++it)
729  {
730  _calo_ID = kHCALOUT;
731  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
732  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
733  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
734  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
735  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
736  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
737  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
738  _calo_towers_N++;
739  }
740  _geometry_done[kHCALOUT] = 1;
741  _geometry_tree->Fill();
742  resetGeometryArrays();
743  }
744  RawTowerContainer::ConstRange begin_end = towersHCALOUT->getTowers();
745  RawTowerContainer::ConstIterator rtiter;
746  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
747  {
748  RawTower* tower = rtiter->second;
749  if (tower)
750  {
751  // min energy cut
752  if (tower->get_energy() < _reco_e_threshold)
753  {
754  continue;
755  }
756  _tower_HCALOUT_iEta[_nTowers_HCALOUT] = tower->get_bineta();
757  _tower_HCALOUT_iPhi[_nTowers_HCALOUT] = tower->get_binphi();
758  _tower_HCALOUT_E[_nTowers_HCALOUT] = tower->get_energy();
759 
760  PHG4Particle* primary = towerevalHCALOUT->max_truth_primary_particle_by_energy(tower);
761  if (primary)
762  {
763  _tower_HCALOUT_trueID[_nTowers_HCALOUT] = primary->get_track_id();
764  // gflavor = primary->get_pid();
765  // efromtruth = towerevalHCALOUT->get_energy_contribution(tower, primary);
766  }
767  else
768  {
769  _tower_HCALOUT_trueID[_nTowers_HCALOUT] = -10;
770  }
771  _nTowers_HCALOUT++;
772  }
773  }
774  }
775  else
776  {
777  if (Verbosity() > 0)
778  {
779  std::cout << PHWHERE << " ERROR: Can't find " << towergeomnodeHCALOUT << std::endl;
780  }
781  // return;
782  }
783  if (Verbosity() > 0)
784  {
785  std::cout << "saved\t" << _nTowers_HCALOUT << "\tHCALOUT towers" << std::endl;
786  }
787  }
788  else
789  {
790  if (Verbosity() > 0)
791  {
792  std::cout << PHWHERE << " ERROR: Can't find " << towernodeHCALOUT << std::endl;
793  }
794  // return;
795  }
796  }
797  //----------------------
798  // TOWERS CEMC
799  //----------------------
800  if (_do_CEMC)
801  {
802  CaloRawTowerEval* towerevalCEMC = _caloevalstackCEMC->get_rawtower_eval();
803  _nTowers_CEMC = 0;
804  std::string towernodeCEMC = "TOWER_CALIB_CEMC";
805  RawTowerContainer* towersCEMC = findNode::getClass<RawTowerContainer>(topNode, towernodeCEMC.c_str());
806  if (towersCEMC)
807  {
808  if (Verbosity() > 0)
809  {
810  std::cout << "saving EMC towers" << std::endl;
811  }
812  std::string towergeomnodeCEMC = "TOWERGEOM_CEMC";
813  RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeCEMC.c_str());
814  if (towergeom)
815  {
816  if (_do_GEOMETRY && !_geometry_done[kCEMC])
817  {
818  RawTowerGeomContainer::ConstRange all_towers = towergeom->get_tower_geometries();
819  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
820  it != all_towers.second; ++it)
821  {
822  _calo_ID = kCEMC;
823  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
824  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
825  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
826  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
827  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
828  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
829  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
830  _calo_towers_N++;
831  }
832  _geometry_done[kCEMC] = 1;
833  _geometry_tree->Fill();
834  resetGeometryArrays();
835  }
836  RawTowerContainer::ConstRange begin_end = towersCEMC->getTowers();
837  RawTowerContainer::ConstIterator rtiter;
838  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
839  {
840  RawTower* tower = rtiter->second;
841  if (tower)
842  {
843  // min energy cut
844  if (tower->get_energy() < _reco_e_threshold)
845  {
846  continue;
847  }
848 
849  _tower_CEMC_iEta[_nTowers_CEMC] = tower->get_bineta();
850  _tower_CEMC_iPhi[_nTowers_CEMC] = tower->get_binphi();
851  _tower_CEMC_E[_nTowers_CEMC] = tower->get_energy();
852 
853  PHG4Particle* primary = towerevalCEMC->max_truth_primary_particle_by_energy(tower);
854  if (primary)
855  {
856  _tower_CEMC_trueID[_nTowers_CEMC] = primary->get_track_id();
857  // gflavor = primary->get_pid();
858  // efromtruth = towerevalCEMC->get_energy_contribution(tower, primary);
859  }
860  else
861  {
862  _tower_CEMC_trueID[_nTowers_CEMC] = -10;
863  }
864  _nTowers_CEMC++;
865  }
866  }
867  }
868  else
869  {
870  if (Verbosity() > 0)
871  {
872  std::cout << PHWHERE << " ERROR: Can't find " << towergeomnodeCEMC << std::endl;
873  }
874  // return;
875  }
876  if (Verbosity() > 0)
877  {
878  std::cout << "saved\t" << _nTowers_CEMC << "\tCEMC towers" << std::endl;
879  }
880  }
881  else
882  {
883  if (Verbosity() > 0)
884  {
885  std::cout << PHWHERE << " ERROR: Can't find " << towernodeCEMC << std::endl;
886  }
887  // return;
888  }
889  }
890 
891  //------------------------
892  // CLUSTERS HCALIN
893  //------------------------
894  if (_do_HCALIN && _do_CLUSTERS)
895  {
896  CaloRawClusterEval* clusterevalHCALIN = _caloevalstackHCALIN->get_rawcluster_eval();
897  _nclusters_HCALIN = 0;
898  if (Verbosity() > 1)
899  {
900  std::cout << "CaloEvaluator::filling gcluster ntuple..." << std::endl;
901  }
902 
903  std::string clusternodeHCALIN = "CLUSTER_HCALIN";
904  RawClusterContainer* clustersHCALIN = findNode::getClass<RawClusterContainer>(topNode, clusternodeHCALIN.c_str());
905  if (clustersHCALIN)
906  {
907  // for every cluster
908  for (const auto& iterator : clustersHCALIN->getClustersMap())
909  {
910  RawCluster* cluster = iterator.second;
911 
912  if (cluster->get_energy() < _reco_e_threshold)
913  {
914  continue;
915  }
916 
917  _cluster_HCALIN_E[_nclusters_HCALIN] = cluster->get_energy();
918  _cluster_HCALIN_NTower[_nclusters_HCALIN] = cluster->getNTowers();
919  _cluster_HCALIN_Phi[_nclusters_HCALIN] = cluster->get_phi();
920 
921  // require vertex for cluster eta calculation
922  if (vertexmap)
923  {
924  if (!vertexmap->empty())
925  {
926  SvtxVertex* vertex = (vertexmap->begin()->second);
927  _cluster_HCALIN_Eta[_nclusters_HCALIN] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
928  }
929  else
930  {
931  _cluster_HCALIN_Eta[_nclusters_HCALIN] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
932  };
933  }
934  else
935  {
936  _cluster_HCALIN_Eta[_nclusters_HCALIN] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
937  };
938 
939  PHG4Particle* primary = clusterevalHCALIN->max_truth_primary_particle_by_energy(cluster);
940 
941  if (primary)
942  {
943  _cluster_HCALIN_trueID[_nclusters_HCALIN] = primary->get_track_id();
944  }
945  else
946  {
947  _cluster_HCALIN_trueID[_nclusters_HCALIN] = -10;
948  }
949 
950  _nclusters_HCALIN++;
951  }
952  }
953  else
954  {
955  std::cout << PHWHERE << " ERROR: Can't find " << clusternodeHCALIN << std::endl;
956  // return;
957  }
958  if (Verbosity() > 0)
959  {
960  std::cout << "saved\t" << _nclusters_HCALIN << "\tHCALIN clusters" << std::endl;
961  }
962  }
963  //------------------------
964  // CLUSTERS HCALOUT
965  //------------------------
966  if (_do_HCALOUT && _do_CLUSTERS)
967  {
968  CaloRawClusterEval* clusterevalHCALOUT = _caloevalstackHCALOUT->get_rawcluster_eval();
969  _nclusters_HCALOUT = 0;
970  if (Verbosity() > 1)
971  {
972  std::cout << "CaloEvaluator::filling gcluster ntuple..." << std::endl;
973  }
974 
975  std::string clusternodeHCALOUT = "CLUSTER_HCALOUT";
976  RawClusterContainer* clustersHCALOUT = findNode::getClass<RawClusterContainer>(topNode, clusternodeHCALOUT.c_str());
977  if (clustersHCALOUT)
978  {
979  // for every cluster
980  for (const auto& iterator : clustersHCALOUT->getClustersMap())
981  {
982  RawCluster* cluster = iterator.second;
983 
984  if (cluster->get_energy() < _reco_e_threshold)
985  {
986  continue;
987  }
988 
989  _cluster_HCALOUT_E[_nclusters_HCALOUT] = cluster->get_energy();
990  _cluster_HCALOUT_NTower[_nclusters_HCALOUT] = cluster->getNTowers();
991  _cluster_HCALOUT_Phi[_nclusters_HCALOUT] = cluster->get_phi();
992 
993  // require vertex for cluster eta calculation
994  if (vertexmap)
995  {
996  if (!vertexmap->empty())
997  {
998  SvtxVertex* vertex = (vertexmap->begin()->second);
999  _cluster_HCALOUT_Eta[_nclusters_HCALOUT] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
1000  }
1001  else
1002  {
1003  _cluster_HCALOUT_Eta[_nclusters_HCALOUT] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
1004  };
1005  }
1006  else
1007  {
1008  _cluster_HCALOUT_Eta[_nclusters_HCALOUT] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
1009  };
1010 
1011  PHG4Particle* primary = clusterevalHCALOUT->max_truth_primary_particle_by_energy(cluster);
1012 
1013  if (primary)
1014  {
1015  _cluster_HCALOUT_trueID[_nclusters_HCALOUT] = primary->get_track_id();
1016  }
1017  else
1018  {
1019  _cluster_HCALOUT_trueID[_nclusters_HCALOUT] = -10;
1020  }
1021 
1022  _nclusters_HCALOUT++;
1023  }
1024  }
1025  else
1026  {
1027  std::cout << PHWHERE << " ERROR: Can't find " << clusternodeHCALOUT << std::endl;
1028  // return;
1029  }
1030  if (Verbosity() > 0)
1031  {
1032  std::cout << "saved\t" << _nclusters_HCALOUT << "\tHCALOUT clusters" << std::endl;
1033  }
1034  }
1035  //------------------------
1036  // CLUSTERS CEMC
1037  //------------------------
1038  if (_do_CEMC && _do_CLUSTERS)
1039  {
1040  CaloRawClusterEval* clusterevalCEMC = _caloevalstackCEMC->get_rawcluster_eval();
1041  _nclusters_CEMC = 0;
1042  if (Verbosity() > 1)
1043  {
1044  std::cout << "CaloEvaluator::filling gcluster ntuple..." << std::endl;
1045  }
1046 
1047  std::string clusternodeCEMC = "CLUSTER_CEMC";
1048  RawClusterContainer* clustersCEMC = findNode::getClass<RawClusterContainer>(topNode, clusternodeCEMC.c_str());
1049  if (clustersCEMC)
1050  {
1051  // for every cluster
1052  for (const auto& iterator : clustersCEMC->getClustersMap())
1053  {
1054  RawCluster* cluster = iterator.second;
1055 
1056  if (cluster->get_energy() < _reco_e_threshold)
1057  {
1058  continue;
1059  }
1060 
1061  _cluster_CEMC_E[_nclusters_CEMC] = cluster->get_energy();
1062  _cluster_CEMC_NTower[_nclusters_CEMC] = cluster->getNTowers();
1063  _cluster_CEMC_Phi[_nclusters_CEMC] = cluster->get_phi();
1064 
1065  // require vertex for cluster eta calculation
1066  if (vertexmap)
1067  {
1068  if (!vertexmap->empty())
1069  {
1070  SvtxVertex* vertex = (vertexmap->begin()->second);
1071  _cluster_CEMC_Eta[_nclusters_CEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
1072  }
1073  else
1074  {
1075  _cluster_CEMC_Eta[_nclusters_CEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
1076  };
1077  }
1078  else
1079  {
1080  _cluster_CEMC_Eta[_nclusters_CEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
1081  };
1082 
1083  PHG4Particle* primary = clusterevalCEMC->max_truth_primary_particle_by_energy(cluster);
1084 
1085  if (primary)
1086  {
1087  _cluster_CEMC_trueID[_nclusters_CEMC] = primary->get_track_id();
1088  }
1089  else
1090  {
1091  _cluster_CEMC_trueID[_nclusters_CEMC] = -10;
1092  }
1093 
1094  _nclusters_CEMC++;
1095  }
1096  }
1097  else
1098  {
1099  std::cout << PHWHERE << " ERROR: Can't find " << clusternodeCEMC << std::endl;
1100  // return;
1101  }
1102  if (Verbosity() > 0)
1103  {
1104  std::cout << "saved\t" << _nclusters_CEMC << "\tCEMC clusters" << std::endl;
1105  }
1106  }
1107 
1108  //------------------------
1109  // TRACKS
1110  //------------------------
1111  if (_do_TRACKS)
1112  {
1113  _nTracks = 0;
1114  _nProjections = 0;
1115  // Loop over track maps, identifiy each source.
1116  // Although this configuration is fixed here, it doesn't require multiple sources.
1117  // It will only store them if they're available.
1118  std::vector<std::pair<std::string, TrackSource_t>> trackMapInfovec = {
1119  {"TrackMap", TrackSource_t::all},
1120  {"TrackMapInner", TrackSource_t::inner}};
1121  bool foundAtLeastOneTrackSource = false;
1122  for (const auto& trackMapInfo : trackMapInfovec)
1123  {
1124  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, trackMapInfo.first);
1125  if (trackmap)
1126  {
1127  foundAtLeastOneTrackSource = true;
1128  int nTracksInASource = 0;
1129  if (Verbosity() > 0)
1130  {
1131  std::cout << "saving tracks for track map: " << trackMapInfo.first << std::endl;
1132  }
1133  for (SvtxTrackMap::ConstIter track_itr = trackmap->begin(); track_itr != trackmap->end(); track_itr++)
1134  {
1135  SvtxTrack_FastSim* track = dynamic_cast<SvtxTrack_FastSim*>(track_itr->second);
1136  if (track)
1137  {
1138  _track_ID[_nTracks] = track->get_id();
1139  _track_px[_nTracks] = track->get_px();
1140  _track_py[_nTracks] = track->get_py();
1141  _track_pz[_nTracks] = track->get_pz();
1142  // Ideally, would be dca3d_xy and dca3d_z, but these don't seem to be calculated properly in the
1143  // current (June 2021) simulations (they return NaN). So we take dca (seems to be ~ the 3d distance)
1144  // and dca_2d (seems to be ~ the distance in the transverse plane).
1145  // The names of the branches are based on the method names.
1146  _track_dca[_nTracks] = static_cast<float>(track->get_dca());
1147  _track_dca_2d[_nTracks] = static_cast<float>(track->get_dca2d());
1148  _track_trueID[_nTracks] = track->get_truth_track_id();
1149  _track_source[_nTracks] = static_cast<unsigned short>(trackMapInfo.second);
1150  if (_do_PROJECTIONS)
1151  {
1152  // find projections
1153  for (SvtxTrack::ConstStateIter trkstates = track->begin_states(); trkstates != track->end_states(); ++trkstates)
1154  {
1155  if (Verbosity() > 1)
1156  {
1157  std::cout << __PRETTY_FUNCTION__ << " processing " << trkstates->second->get_name() << std::endl;
1158  }
1159  std::string trackStateName = trkstates->second->get_name();
1160  if (Verbosity() > 1)
1161  {
1162  std::cout << __PRETTY_FUNCTION__ << " found " << trkstates->second->get_name() << std::endl;
1163  }
1164  int trackStateIndex = GetProjectionIndex(trackStateName);
1165  if (trackStateIndex > -1)
1166  {
1167  // save true projection info to given branch
1168  _track_TLP_true_x[_nProjections] = trkstates->second->get_pos(0);
1169  _track_TLP_true_y[_nProjections] = trkstates->second->get_pos(1);
1170  _track_TLP_true_z[_nProjections] = trkstates->second->get_pos(2);
1171  _track_TLP_true_t[_nProjections] = trkstates->first;
1172  _track_ProjLayer[_nProjections] = trackStateIndex;
1173  _track_ProjTrackID[_nProjections] = _nTracks;
1174 
1175  std::string nodename = "G4HIT_" + trkstates->second->get_name();
1176  PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
1177  if (hits)
1178  {
1179  if (Verbosity() > 1)
1180  {
1181  std::cout << __PRETTY_FUNCTION__ << " number of hits: " << hits->size() << std::endl;
1182  }
1183  PHG4HitContainer::ConstRange hit_range = hits->getHits();
1184  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
1185  {
1186  if (Verbosity() > 1)
1187  {
1188  std::cout << __PRETTY_FUNCTION__ << " checking hit id " << hit_iter->second->get_trkid() << " against " << track->get_truth_track_id() << std::endl;
1189  }
1190  if (hit_iter->second->get_trkid() - track->get_truth_track_id() == 0)
1191  {
1192  if (Verbosity() > 1)
1193  {
1194  std::cout << __PRETTY_FUNCTION__ << " found hit with id " << hit_iter->second->get_trkid() << std::endl;
1195  }
1196  // save reco projection info to given branch
1197  _track_TLP_x[_nProjections] = hit_iter->second->get_x(0);
1198  _track_TLP_y[_nProjections] = hit_iter->second->get_y(0);
1199  _track_TLP_z[_nProjections] = hit_iter->second->get_z(0);
1200  _track_TLP_t[_nProjections] = hit_iter->second->get_t(0);
1201  }
1202  }
1203  }
1204  else
1205  {
1206  if (Verbosity() > 1)
1207  {
1208  std::cout << __PRETTY_FUNCTION__ << " could not find " << nodename << std::endl;
1209  }
1210  continue;
1211  }
1212  _nProjections++;
1213  }
1214  }
1215  }
1216  _nTracks++;
1217  nTracksInASource++;
1218  }
1219  else
1220  {
1221  if (Verbosity() > 0) // Verbosity()
1222  {
1223  std::cout << "PHG4TrackFastSimEval::fill_track_tree - ignore track that is not a SvtxTrack_FastSim:";
1224  track_itr->second->identify();
1225  }
1226  continue;
1227  }
1228  }
1229  if (Verbosity() > 0)
1230  {
1231  std::cout << "saved\t" << nTracksInASource << "\ttracks from track map " << trackMapInfo.first << ". Total saved tracks: " << _nTracks << std::endl;
1232  }
1233  }
1234  else
1235  {
1236  if (Verbosity() > 0)
1237  {
1238  std::cout << PHWHERE << "SvtxTrackMap node with name '" << trackMapInfo.first << "' not found on node tree" << std::endl;
1239  }
1240  }
1241  }
1242  if (foundAtLeastOneTrackSource == false)
1243  {
1244  std::cout << PHWHERE << "Requested tracks, but found no sources on node tree. Returning" << std::endl;
1245  return;
1246  }
1247  }
1248  //------------------------
1249  // MC PARTICLES
1250  //------------------------
1251  _nMCPart = 0;
1252  if (_do_MCPARTICLES)
1253  {
1254  PHG4TruthInfoContainer* truthinfocontainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1255  if (truthinfocontainer)
1256  {
1257  if (Verbosity() > 0)
1258  {
1259  std::cout << "saving MC particles" << std::endl;
1260  }
1261  // GetParticleRange for all particles
1262  // GetPrimaryParticleRange for primary particles
1263  PHG4TruthInfoContainer::ConstRange range = truthinfocontainer->GetParticleRange();
1264  for (PHG4TruthInfoContainer::ConstIterator truth_itr = range.first; truth_itr != range.second; ++truth_itr)
1265  {
1266  PHG4Particle* g4particle = truth_itr->second;
1267  if (!g4particle)
1268  {
1269  continue;
1270  }
1271 
1272  int mcSteps = 0;
1273  PHG4Particle* g4particleMother = truth_itr->second;
1274  if (g4particle->get_parent_id() != 0)
1275  {
1276  while (g4particleMother->get_parent_id() != 0)
1277  {
1278  g4particleMother = truthinfocontainer->GetParticle(g4particleMother->get_parent_id());
1279  if (g4particleMother == nullptr)
1280  {
1281  break;
1282  }
1283  mcSteps += 1;
1284  }
1285  }
1286  if (mcSteps > _depth_MCstack)
1287  {
1288  continue;
1289  }
1290 
1291  // evaluating true primary vertex
1292  if (_do_VERTEX && _nMCPart == 0)
1293  {
1294  PHG4VtxPoint* vtx = truthinfocontainer->GetVtx(g4particle->get_vtx_id());
1295  if (vtx)
1296  {
1297  _vertex_true_x = vtx->get_x();
1298  _vertex_true_y = vtx->get_y();
1299  _vertex_true_z = vtx->get_z();
1300  }
1301  }
1302 
1303  // in case of all MC particles, make restrictions on the secondary selection
1304  // if(g4particle->get_track_id()<0 && g4particle->get_e()<0.5) continue;
1305  // primary (g4particle->get_parent_id() == 0) selection via:
1306  // if(gtrackID < 0) continue;
1307 
1308  // using the e threshold also for the truth particles gets rid of all the low energy secondary particles
1309  if (g4particle->get_e() < _reco_e_threshold)
1310  {
1311  continue;
1312  }
1313 
1314  _mcpart_ID[_nMCPart] = g4particle->get_track_id();
1315  _mcpart_ID_parent[_nMCPart] = g4particle->get_parent_id();
1316  _mcpart_PDG[_nMCPart] = g4particle->get_pid();
1317  _mcpart_E[_nMCPart] = g4particle->get_e();
1318  _mcpart_px[_nMCPart] = g4particle->get_px();
1319  _mcpart_py[_nMCPart] = g4particle->get_py();
1320  _mcpart_pz[_nMCPart] = g4particle->get_pz();
1321  // BCID added for G4Particle -- HEPMC particle matching
1322  _mcpart_BCID[_nMCPart] = g4particle->get_barcode();
1323  // TVector3 projvec(_mcpart_px[0],_mcpart_py[0],_mcpart_pz[0]);
1324  // float projeta = projvec.Eta();
1325  _nMCPart++;
1326  }
1327  if (Verbosity() > 0)
1328  {
1329  std::cout << "saved\t" << _nMCPart << "\tMC particles" << std::endl;
1330  }
1331  }
1332  else
1333  {
1334  if (Verbosity() > 0)
1335  {
1336  std::cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree" << std::endl;
1337  }
1338  return;
1339  }
1340  }
1341 
1342  //------------------------
1343  // HEPMC
1344  //------------------------
1345  _nHepmcp = 0;
1346  if (_do_HEPMC)
1347  {
1348  PHHepMCGenEventMap* hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
1349  if (hepmceventmap)
1350  {
1351  if (Verbosity() > 0)
1352  {
1353  std::cout << "saving HepMC output" << std::endl;
1354  }
1355  if (Verbosity() > 0)
1356  {
1357  hepmceventmap->Print();
1358  }
1359 
1360  for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
1361  eventIter != hepmceventmap->end();
1362  ++eventIter)
1363  {
1364  PHHepMCGenEvent* hepmcevent = eventIter->second;
1365 
1366  if (hepmcevent)
1367  {
1368  HepMC::GenEvent* truthevent = hepmcevent->getEvent();
1369  if (!truthevent)
1370  {
1371  std::cout << PHWHERE
1372  << "no evt pointer under phhepmvgeneventmap found "
1373  << std::endl;
1374  return;
1375  }
1376 
1377  HepMC::PdfInfo* pdfinfo = truthevent->pdf_info();
1378 
1379  // m_partid1 = pdfinfo->id1();
1380  // m_partid2 = pdfinfo->id2();
1381  _hepmcp_x1 = pdfinfo->x1();
1382  _hepmcp_x2 = pdfinfo->x2();
1383 
1384  // m_mpi = truthevent->mpi();
1385 
1386  _hepmcp_procid = truthevent->signal_process_id();
1387 
1388  if (Verbosity() > 2)
1389  {
1390  std::cout << " Iterating over an event" << std::endl;
1391  }
1392  for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
1393  iter != truthevent->particles_end();
1394  ++iter)
1395  {
1396  _hepmcp_E[_nHepmcp] = (*iter)->momentum().e();
1397  _hepmcp_PDG[_nHepmcp] = (*iter)->pdg_id();
1398  _hepmcp_px[_nHepmcp] = (*iter)->momentum().px();
1399  _hepmcp_py[_nHepmcp] = (*iter)->momentum().py();
1400  _hepmcp_pz[_nHepmcp] = (*iter)->momentum().pz();
1401  _hepmcp_status[_nHepmcp] = (*iter)->status();
1402  _hepmcp_BCID[_nHepmcp] = (*iter)->barcode();
1403  _hepmcp_m2[_nHepmcp] = 0;
1404  _hepmcp_m1[_nHepmcp] = 0;
1405  if ((*iter)->production_vertex())
1406  {
1407  for (HepMC::GenVertex::particle_iterator mother = (*iter)->production_vertex()->particles_begin(HepMC::parents);
1408  mother != (*iter)->production_vertex()->particles_end(HepMC::parents);
1409  ++mother)
1410  {
1411  _hepmcp_m2[_nHepmcp] = (*mother)->barcode();
1412  if (_hepmcp_m1[_nHepmcp] == 0)
1413  {
1414  _hepmcp_m1[_nHepmcp] = (*mother)->barcode();
1415  }
1416  }
1417  }
1418  if (Verbosity() > 2)
1419  {
1420  std::cout << "nHepmcp " << _nHepmcp << "\tPDG " << _hepmcp_PDG[_nHepmcp] << "\tEnergy " << _hepmcp_E[_nHepmcp] << "\tbarcode " << _hepmcp_BCID[_nHepmcp] << "\tMother1 " << _hepmcp_m1[_nHepmcp] << "\tMother2 " << _hepmcp_m2[_nHepmcp] << std::endl;
1421  }
1422  _nHepmcp++;
1423  }
1424  }
1425  }
1426  }
1427  else
1428  {
1429  if (Verbosity() > 0)
1430  {
1431  std::cout << PHWHERE << " PHHepMCGenEventMap node not found on node tree" << std::endl;
1432  }
1433  return;
1434  }
1435  } // hepmc
1436 
1437  _event_tree->Fill();
1438 
1439  if (Verbosity() > 0)
1440  {
1441  std::cout << "Resetting buffer ..." << std::endl;
1442  }
1443  resetBuffer();
1444  if (Verbosity() > 0)
1445  {
1446  std::cout << "EventEvaluator buffer reset" << std::endl;
1447  }
1448  return;
1449 }
1450 
1451 int EventEvaluator::End(PHCompositeNode* /*topNode*/)
1452 {
1453  _tfile->cd();
1454 
1455  _event_tree->Write();
1456 
1457  _tfile->Close();
1458 
1459  delete _tfile;
1460 
1461  if (_do_GEOMETRY)
1462  {
1463  _tfile_geometry->cd();
1464 
1465  _geometry_tree->Write();
1466 
1467  _tfile_geometry->Close();
1468 
1469  delete _tfile_geometry;
1470  }
1471  if (Verbosity() > 0)
1472  {
1473  std::cout << "========================= " << Name() << "::End() ============================" << std::endl;
1474  std::cout << " " << _ievent << " events of output written to: " << _filename << std::endl;
1475  std::cout << "===========================================================================" << std::endl;
1476  }
1477 
1478  if (_caloevalstackHCALIN)
1479  {
1480  delete _caloevalstackHCALIN;
1481  }
1482  if (_caloevalstackHCALOUT)
1483  {
1484  delete _caloevalstackHCALOUT;
1485  }
1486  if (_caloevalstackCEMC)
1487  {
1488  delete _caloevalstackCEMC;
1489  }
1490 
1491  return Fun4AllReturnCodes::EVENT_OK;
1492 }
1493 
1494 int EventEvaluator::GetProjectionIndex(const std::string& projname)
1495 {
1496  if (projname.find("HCALIN") != std::string::npos)
1497  {
1498  return 1;
1499  }
1500  else if (projname.find("HCALOUT") != std::string::npos)
1501  {
1502  return 2;
1503  }
1504  else if (projname.find("CEMC") != std::string::npos)
1505  {
1506  return 3;
1507  }
1508  else
1509  {
1510  return -1;
1511  }
1512  return -1;
1513 }
1514 
1515 std::string EventEvaluator::GetProjectionNameFromIndex(int projindex)
1516 {
1517  switch (projindex)
1518  {
1519  case 1:
1520  return "HCALIN";
1521  case 2:
1522  return "HCALOUT";
1523  case 3:
1524  return "CEMC";
1525  default:
1526  return "NOTHING";
1527  }
1528 }
1529 
1530 void EventEvaluator::resetGeometryArrays()
1531 {
1532  for (Int_t igeo = 0; igeo < _calo_towers_N; igeo++)
1533  {
1534  _calo_towers_iEta[_calo_towers_N] = -10000;
1535  _calo_towers_iPhi[_calo_towers_N] = -10000;
1536  _calo_towers_Eta[_calo_towers_N] = -10000;
1537  _calo_towers_Phi[_calo_towers_N] = -10000;
1538  _calo_towers_x[_calo_towers_N] = -10000;
1539  _calo_towers_y[_calo_towers_N] = -10000;
1540  _calo_towers_z[_calo_towers_N] = -10000;
1541  }
1542  _calo_ID = -1;
1543  _calo_towers_N = 0;
1544 }
1545 void EventEvaluator::resetBuffer()
1546 {
1547  if (_do_store_event_info)
1548  {
1549  _cross_section = 0;
1550  _event_weight = 0;
1551  _n_generator_accepted = 0;
1552  if (Verbosity() > 0)
1553  {
1554  std::cout << "\t... event info variables reset" << std::endl;
1555  }
1556  }
1557  if (_do_VERTEX)
1558  {
1559  _vertex_x = -1000;
1560  _vertex_y = -1000;
1561  _vertex_z = -1000;
1562  _vertex_NCont = 0;
1563  _vertex_true_x = -1000;
1564  _vertex_true_y = -1000;
1565  _vertex_true_z = -1000;
1566  if (Verbosity() > 0)
1567  {
1568  std::cout << "\t... vertex variables reset" << std::endl;
1569  }
1570  }
1571  if (_do_HITS)
1572  {
1573  _nHitsLayers = 0;
1574  for (Int_t ihit = 0; ihit < _maxNHits; ihit++)
1575  {
1576  _hits_layerID[ihit] = 0;
1577  _hits_trueID[ihit] = 0;
1578  _hits_x[ihit] = 0;
1579  _hits_y[ihit] = 0;
1580  _hits_z[ihit] = 0;
1581  _hits_t[ihit] = 0;
1582  }
1583  if (Verbosity() > 0)
1584  {
1585  std::cout << "\t... hit variables reset" << std::endl;
1586  }
1587  }
1588  if (_do_CEMC)
1589  {
1590  _nTowers_CEMC = 0;
1591  for (Int_t itow = 0; itow < _maxNTowersCentral; itow++)
1592  {
1593  _tower_CEMC_E[itow] = 0;
1594  _tower_CEMC_iEta[itow] = 0;
1595  _tower_CEMC_iPhi[itow] = 0;
1596  _tower_CEMC_trueID[itow] = 0;
1597  }
1598  if (_do_CLUSTERS)
1599  {
1600  _nclusters_CEMC = 0;
1601  for (Int_t itow = 0; itow < _maxNclustersCentral; itow++)
1602  {
1603  _cluster_CEMC_E[itow] = 0;
1604  _cluster_CEMC_Eta[itow] = 0;
1605  _cluster_CEMC_Phi[itow] = 0;
1606  _cluster_CEMC_NTower[itow] = 0;
1607  _cluster_CEMC_trueID[itow] = 0;
1608  }
1609  }
1610  if (Verbosity() > 0)
1611  {
1612  std::cout << "\t... CEMC variables reset" << std::endl;
1613  }
1614  }
1615  if (_do_HCALIN)
1616  {
1617  _nTowers_HCALIN = 0;
1618  for (Int_t itow = 0; itow < _maxNTowersCentral; itow++)
1619  {
1620  _tower_HCALIN_E[itow] = 0;
1621  _tower_HCALIN_iEta[itow] = 0;
1622  _tower_HCALIN_iPhi[itow] = 0;
1623  _tower_HCALIN_trueID[itow] = 0;
1624  }
1625  if (_do_CLUSTERS)
1626  {
1627  _nclusters_HCALIN = 0;
1628  for (Int_t itow = 0; itow < _maxNclustersCentral; itow++)
1629  {
1630  _cluster_HCALIN_E[itow] = 0;
1631  _cluster_HCALIN_Eta[itow] = 0;
1632  _cluster_HCALIN_Phi[itow] = 0;
1633  _cluster_HCALIN_NTower[itow] = 0;
1634  _cluster_HCALIN_trueID[itow] = 0;
1635  }
1636  }
1637  if (Verbosity() > 0)
1638  {
1639  std::cout << "\t... HCALIN variables reset" << std::endl;
1640  }
1641  }
1642  if (_do_HCALOUT)
1643  {
1644  if (Verbosity() > 0)
1645  {
1646  std::cout << "\t... resetting HCALOUT variables" << std::endl;
1647  }
1648  _nTowers_HCALOUT = 0;
1649  for (Int_t itow = 0; itow < _maxNTowersCentral; itow++)
1650  {
1651  _tower_HCALOUT_E[itow] = 0;
1652  _tower_HCALOUT_iEta[itow] = 0;
1653  _tower_HCALOUT_iPhi[itow] = 0;
1654  _tower_HCALOUT_trueID[itow] = 0;
1655  }
1656  if (_do_CLUSTERS)
1657  {
1658  _nclusters_HCALOUT = 0;
1659  for (Int_t itow = 0; itow < _maxNclustersCentral; itow++)
1660  {
1661  _cluster_HCALOUT_E[itow] = 0;
1662  _cluster_HCALOUT_Eta[itow] = 0;
1663  _cluster_HCALOUT_Phi[itow] = 0;
1664  _cluster_HCALOUT_NTower[itow] = 0;
1665  _cluster_HCALOUT_trueID[itow] = 0;
1666  }
1667  }
1668  if (Verbosity() > 0)
1669  {
1670  std::cout << "\t... HCALOUT variables reset" << std::endl;
1671  }
1672  }
1673  if (_do_TRACKS)
1674  {
1675  if (Verbosity() > 0)
1676  {
1677  std::cout << "\t... resetting Track variables" << std::endl;
1678  }
1679  _nTracks = 0;
1680  for (Int_t itrk = 0; itrk < _maxNTracks; itrk++)
1681  {
1682  _track_ID[itrk] = 0;
1683  _track_trueID[itrk] = 0;
1684  _track_px[itrk] = 0;
1685  _track_py[itrk] = 0;
1686  _track_pz[itrk] = 0;
1687  _track_dca[itrk] = 0;
1688  _track_dca_2d[itrk] = 0;
1689  _track_source[itrk] = 0;
1690  }
1691  if (_do_PROJECTIONS)
1692  {
1693  _nProjections = 0;
1694  for (Int_t iproj = 0; iproj < _maxNProjections; iproj++)
1695  {
1696  _track_ProjLayer[iproj] = -1;
1697  _track_ProjTrackID[iproj] = 0;
1698  _track_TLP_x[iproj] = 0;
1699  _track_TLP_y[iproj] = 0;
1700  _track_TLP_z[iproj] = 0;
1701  _track_TLP_t[iproj] = 0;
1702  _track_TLP_true_x[iproj] = 0;
1703  _track_TLP_true_y[iproj] = 0;
1704  _track_TLP_true_z[iproj] = 0;
1705  _track_TLP_true_t[iproj] = 0;
1706  }
1707  }
1708  if (Verbosity() > 0)
1709  {
1710  std::cout << "\t... track variables reset" << std::endl;
1711  }
1712  }
1713  if (_do_MCPARTICLES)
1714  {
1715  _nMCPart = 0;
1716  for (Int_t imcpart = 0; imcpart < _maxNMCPart; imcpart++)
1717  {
1718  _mcpart_ID[imcpart] = 0;
1719  _mcpart_ID_parent[imcpart] = 0;
1720  _mcpart_PDG[imcpart] = 0;
1721  _mcpart_E[imcpart] = 0;
1722  _mcpart_px[imcpart] = 0;
1723  _mcpart_py[imcpart] = 0;
1724  _mcpart_pz[imcpart] = 0;
1725  _mcpart_BCID[imcpart] = -10;
1726  }
1727  }
1728 
1729  if (_do_HEPMC)
1730  {
1731  _nHepmcp = 0;
1732  for (Int_t iHepmcp = 0; iHepmcp < _maxNHepmcp; iHepmcp++)
1733  {
1734  _hepmcp_E[iHepmcp] = 0;
1735  _hepmcp_PDG[iHepmcp] = 0;
1736  _hepmcp_px[iHepmcp] = 0;
1737  _hepmcp_py[iHepmcp] = 0;
1738  _hepmcp_pz[iHepmcp] = 0;
1739  _hepmcp_status[iHepmcp] = -10;
1740  _hepmcp_BCID[iHepmcp] = 0;
1741  _hepmcp_m2[iHepmcp] = 0;
1742  _hepmcp_m1[iHepmcp] = 0;
1743  }
1744  if (Verbosity() > 0)
1745  {
1746  std::cout << "\t... MC variables reset" << std::endl;
1747  }
1748  }
1749 }