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