EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EventEvaluatorEIC.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EventEvaluatorEIC.cc
1 #include "EventEvaluatorEIC.h"
2 
3 #include "g4eval/CaloEvalStack.h"
4 #include "g4eval/CaloRawClusterEval.h"
5 #include "g4eval/CaloRawTowerEval.h"
6 #include "g4eval/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 <eicpidbase/EICPIDParticle.h>
34 #include <eicpidbase/EICPIDParticleContainer.h>
35 
36 #include <HepMC/GenEvent.h>
37 #include <HepMC/GenVertex.h>
38 #include <phhepmc/PHGenIntegral.h>
39 #include <phhepmc/PHHepMCGenEvent.h>
40 #include <phhepmc/PHHepMCGenEventMap.h>
41 
42 #include <fun4all/Fun4AllReturnCodes.h>
43 #include <fun4all/SubsysReco.h>
44 
45 #include <phool/PHCompositeNode.h>
46 #include <phool/PHNodeIterator.h> // for PHNodeIterator
47 #include <phool/getClass.h>
48 #include <phool/phool.h>
49 
50 #include <TFile.h>
51 #include <TNtuple.h>
52 #include <TTree.h>
53 
54 #include <CLHEP/Vector/ThreeVector.h>
55 
56 #include <cmath>
57 #include <cstdlib>
58 #include <iostream>
59 #include <set>
60 #include <utility>
61 // #include <fstream>
62 
63 using namespace std;
64 
65 EventEvaluatorEIC::EventEvaluatorEIC(const string& name, const string& filename)
66  : SubsysReco(name)
67  , _do_store_event_info(false)
68  , _do_FHCAL(false)
69  , _do_BECAL(false)
70  , _do_HCALIN(false)
71  , _do_HCALOUT(false)
72  , _do_EHCAL(false)
73  , _do_FEMC(false)
74  , _do_CEMC(false)
75  , _do_EEMC(false)
76  , _do_EEMCG(false)
77  , _do_DRCALO(false)
78  , _do_FOCAL(false)
79  , _do_LFHCAL(false)
80  , _do_HITS(false)
81  , _do_HITS_ABSORBER(false)
82  , _do_HITS_CALO(false)
83  , _do_TRACKS(false)
84  , _do_CLUSTERS(false)
85  , _do_VERTEX(false)
86  , _do_PROJECTIONS(false)
87  , _do_MCPARTICLES(false)
88  , _do_HEPMC(false)
89  , _do_GEOMETRY(false)
90  , _do_BLACKHOLE(false)
91  , _ievent(0)
92  , _cross_section(0)
93  , _event_weight(0)
94  , _n_generator_accepted(0)
95  , _nHitsLayers(0)
96  , _hits_layerID(0)
97  , _hits_trueID(0)
98  , _hits_x(0)
99  , _hits_y(0)
100  , _hits_z(0)
101  , _hits_x2(0)
102  , _hits_y2(0)
103  , _hits_z2(0)
104  , _hits_t(0)
105  , _hits_edep(0)
106  , _hits_lightyield(0)
107  , _hits_isAbsorber(0)
108 
109  , _nTowers_FHCAL(0)
110  , _tower_FHCAL_E(0)
111  , _tower_FHCAL_iEta(0)
112  , _tower_FHCAL_iPhi(0)
113  , _tower_FHCAL_trueID(0)
114 
115  , _nTowers_BECAL(0)
116  , _tower_BECAL_E(0)
117  , _tower_BECAL_iEta(0)
118  , _tower_BECAL_iPhi(0)
119  , _tower_BECAL_trueID(0)
120 
121  , _nTowers_HCALIN(0)
122  , _tower_HCALIN_E(0)
123  , _tower_HCALIN_iEta(0)
124  , _tower_HCALIN_iPhi(0)
125  , _tower_HCALIN_trueID(0)
126 
127  , _nTowers_HCALOUT(0)
128  , _tower_HCALOUT_E(0)
129  , _tower_HCALOUT_iEta(0)
130  , _tower_HCALOUT_iPhi(0)
131  , _tower_HCALOUT_trueID(0)
132 
133  , _nTowers_EHCAL(0)
134  , _tower_EHCAL_E(0)
135  , _tower_EHCAL_iEta(0)
136  , _tower_EHCAL_iPhi(0)
137  , _tower_EHCAL_trueID(0)
138 
139  , _nTowers_DRCALO(0)
140  , _tower_DRCALO_E(0)
141  , _tower_DRCALO_NScint(0)
142  , _tower_DRCALO_NCerenkov(0)
143  , _tower_DRCALO_iEta(0)
144  , _tower_DRCALO_iPhi(0)
145  , _tower_DRCALO_trueID(0)
146 
147  , _nTowers_FOCAL(0)
148  , _tower_FOCAL_E(0)
149  , _tower_FOCAL_NScint(0)
150  , _tower_FOCAL_NCerenkov(0)
151  , _tower_FOCAL_iEta(0)
152  , _tower_FOCAL_iPhi(0)
153  , _tower_FOCAL_trueID(0)
154 
155  , _nTowers_LFHCAL(0)
156  , _tower_LFHCAL_E(0)
157  , _tower_LFHCAL_iEta(0)
158  , _tower_LFHCAL_iPhi(0)
159  , _tower_LFHCAL_iL(0)
160  , _tower_LFHCAL_trueID(0)
161 
162  , _nTowers_FEMC(0)
163  , _tower_FEMC_E(0)
164  , _tower_FEMC_iEta(0)
165  , _tower_FEMC_iPhi(0)
166  , _tower_FEMC_trueID(0)
167 
168  , _nTowers_EEMC(0)
169  , _tower_EEMC_E(0)
170  , _tower_EEMC_iEta(0)
171  , _tower_EEMC_iPhi(0)
172  , _tower_EEMC_trueID(0)
173 
174  , _nTowers_EEMCG(0)
175  , _tower_EEMCG_E(0)
176  , _tower_EEMCG_iEta(0)
177  , _tower_EEMCG_iPhi(0)
178  , _tower_EEMCG_trueID(0)
179 
180  , _nTowers_CEMC(0)
181  , _tower_CEMC_E(0)
182  , _tower_CEMC_iEta(0)
183  , _tower_CEMC_iPhi(0)
184  , _tower_CEMC_trueID(0)
185 
186  , _nclusters_FHCAL(0)
187  , _cluster_FHCAL_E(0)
188  , _cluster_FHCAL_Eta(0)
189  , _cluster_FHCAL_Phi(0)
190  , _cluster_FHCAL_NTower(0)
191  , _cluster_FHCAL_trueID(0)
192 
193  , _nclusters_HCALIN(0)
194  , _cluster_HCALIN_E(0)
195  , _cluster_HCALIN_Eta(0)
196  , _cluster_HCALIN_Phi(0)
197  , _cluster_HCALIN_NTower(0)
198  , _cluster_HCALIN_trueID(0)
199 
200  , _nclusters_HCALOUT(0)
201  , _cluster_HCALOUT_E(0)
202  , _cluster_HCALOUT_Eta(0)
203  , _cluster_HCALOUT_Phi(0)
204  , _cluster_HCALOUT_NTower(0)
205  , _cluster_HCALOUT_trueID(0)
206 
207  , _nclusters_EHCAL(0)
208  , _cluster_EHCAL_E(0)
209  , _cluster_EHCAL_Eta(0)
210  , _cluster_EHCAL_Phi(0)
211  , _cluster_EHCAL_NTower(0)
212  , _cluster_EHCAL_trueID(0)
213 
214  , _nclusters_FEMC(0)
215  , _cluster_FEMC_E(0)
216  , _cluster_FEMC_Eta(0)
217  , _cluster_FEMC_Phi(0)
218  , _cluster_FEMC_NTower(0)
219  , _cluster_FEMC_trueID(0)
220 
221  , _nclusters_CEMC(0)
222  , _cluster_CEMC_E(0)
223  , _cluster_CEMC_Eta(0)
224  , _cluster_CEMC_Phi(0)
225  , _cluster_CEMC_NTower(0)
226  , _cluster_CEMC_trueID(0)
227 
228  , _nclusters_EEMC(0)
229  , _cluster_EEMC_E(0)
230  , _cluster_EEMC_Eta(0)
231  , _cluster_EEMC_Phi(0)
232  , _cluster_EEMC_NTower(0)
233  , _cluster_EEMC_trueID(0)
234 
235  , _nclusters_EEMCG(0)
236  , _cluster_EEMCG_E(0)
237  , _cluster_EEMCG_Eta(0)
238  , _cluster_EEMCG_Phi(0)
239  , _cluster_EEMCG_NTower(0)
240  , _cluster_EEMCG_trueID(0)
241 
242  , _vertex_x(0)
243  , _vertex_y(0)
244  , _vertex_z(0)
245  , _vertex_NCont(0)
246  , _vertex_true_x(0)
247  , _vertex_true_y(0)
248  , _vertex_true_z(0)
249 
250  , _nTracks(0)
251  , _track_ID(0)
252  , _track_charge(0)
253  , _track_px(0)
254  , _track_py(0)
255  , _track_pz(0)
256  , _track_x(0)
257  , _track_y(0)
258  , _track_z(0)
259  , _track_ndf(0)
260  , _track_chi2(0)
261  , _track_dca(0)
262  , _track_dca_2d(0)
263  , _track_trueID(0)
264  , _track_source(0)
265  , _nProjections(0)
266  , _track_ProjTrackID(0)
267  , _track_ProjLayer(0)
268  , _track_TLP_x(0)
269  , _track_TLP_y(0)
270  , _track_TLP_z(0)
271  , _track_TLP_t(0)
272  , _track_TLP_px(0)
273  , _track_TLP_py(0)
274  , _track_TLP_pz(0)
275  , _track_TLP_true_x(0)
276  , _track_TLP_true_y(0)
277  , _track_TLP_true_z(0)
278  , _track_TLP_true_t(0)
279 
280  , _nMCPart(0)
281  , _mcpart_ID(0)
282  , _mcpart_ID_parent(0)
283  , _mcpart_PDG(0)
284  , _mcpart_E(0)
285  , _mcpart_px(0)
286  , _mcpart_py(0)
287  , _mcpart_pz(0)
288  , _mcpart_x(0)
289  , _mcpart_y(0)
290  , _mcpart_z(0)
291  , _mcpart_BCID(0)
292 
293  , _nHepmcp(0)
294  , _hepmcp_procid(0)
295  , _hepmcp_x1(NAN)
296  , _hepmcp_x2(NAN)
297  , _hepmcp_Q2(NAN)
298  , _hepmcp_vtx_x(NAN)
299  , _hepmcp_vtx_y(NAN)
300  , _hepmcp_vtx_z(NAN)
301  , _hepmcp_vtx_t(NAN)
302 
303  // , _hepmcp_ID_parent(0)
304  , _hepmcp_status(0)
305  , _hepmcp_PDG(0)
306  , _hepmcp_E(0)
307  , _hepmcp_px(0)
308  , _hepmcp_py(0)
309  , _hepmcp_pz(0)
310  , _hepmcp_m1(0)
311  , _hepmcp_m2(0)
312  , _hepmcp_BCID(0)
313 
314  , _calo_ID(0)
315  , _calo_towers_N(0)
316  , _calo_towers_iEta(0)
317  , _calo_towers_iPhi(0)
318  , _calo_towers_iL(0)
319  , _calo_towers_Eta(0)
320  , _calo_towers_Phi(0)
321  , _calo_towers_x(0)
322  , _calo_towers_y(0)
323  , _calo_towers_z(0)
324  , _geometry_done(0)
325 
326  , _reco_e_threshold(0)
327  , _reco_e_thresholdMC(0.001)
328  , _depth_MCstack(0)
329  , _caloevalstackFHCAL(nullptr)
330  , _caloevalstackBECAL(nullptr)
331  , _caloevalstackHCALIN(nullptr)
332  , _caloevalstackHCALOUT(nullptr)
333  , _caloevalstackEHCAL(nullptr)
334  , _caloevalstackDRCALO(nullptr)
335  , _caloevalstackFOCAL(nullptr)
336  , _caloevalstackLFHCAL(nullptr)
337  , _caloevalstackFEMC(nullptr)
338  , _caloevalstackCEMC(nullptr)
339  , _caloevalstackEEMC(nullptr)
340  , _caloevalstackEEMCG(nullptr)
341  , _strict(false)
342  , _event_tree(nullptr)
343  , _geometry_tree(nullptr)
344  , _filename(filename)
345  , _tfile(nullptr)
346  , _tfile_geometry(nullptr)
347 {
348  _reco_e_threshold = new float[_maxNCalo];
349  _reco_e_threshold[kFHCAL] = 0.05;
350  _reco_e_threshold[kFEMC] = 0.005;
351  _reco_e_threshold[kDRCALO] = 0.0;
352  _reco_e_threshold[kFOCAL] = 0.0;
353  _reco_e_threshold[kEEMC] = 0.005;
354  _reco_e_threshold[kCEMC] = 0.01;
355  _reco_e_threshold[kEHCAL] = 0.05;
356  _reco_e_threshold[kHCALIN] = 0.01;
357  _reco_e_threshold[kHCALOUT] = 0.05;
358  _reco_e_threshold[kLFHCAL] = 0.001;
359  _reco_e_threshold[kEEMCG] = 0.005;
360  _reco_e_threshold[kBECAL] = 0.001;
361 
362  _hits_layerID = new int[_maxNHits];
363  _hits_trueID = new int[_maxNHits];
364  _hits_x = new float[_maxNHits];
365  _hits_y = new float[_maxNHits];
366  _hits_z = new float[_maxNHits];
367  _hits_x2 = new float[_maxNHits];
368  _hits_y2 = new float[_maxNHits];
369  _hits_z2 = new float[_maxNHits];
370  _hits_t = new float[_maxNHits];
371  _hits_edep = new float[_maxNHits];
372  _hits_lightyield = new float[_maxNHits];
373  _hits_isAbsorber = new int[_maxNHits];
374 
375  _tower_FHCAL_E = new float[_maxNTowers];
376  _tower_FHCAL_iEta = new int[_maxNTowers];
377  _tower_FHCAL_iPhi = new int[_maxNTowers];
378  _tower_FHCAL_trueID = new int[_maxNTowers];
379  _cluster_FHCAL_E = new float[_maxNclusters];
380  _cluster_FHCAL_Eta = new float[_maxNclusters];
381  _cluster_FHCAL_Phi = new float[_maxNclusters];
382  _cluster_FHCAL_NTower = new int[_maxNclusters];
383  _cluster_FHCAL_trueID = new int[_maxNclusters];
384 
385  _tower_BECAL_E = new float[_maxNTowers];
386  _tower_BECAL_iEta = new int[_maxNTowers];
387  _tower_BECAL_iPhi = new int[_maxNTowers];
388  _tower_BECAL_trueID = new int[_maxNTowers];
389 
390  _tower_HCALIN_E = new float[_maxNTowersCentral];
391  _tower_HCALIN_iEta = new int[_maxNTowersCentral];
392  _tower_HCALIN_iPhi = new int[_maxNTowersCentral];
393  _tower_HCALIN_trueID = new int[_maxNTowersCentral];
394  _cluster_HCALIN_E = new float[_maxNclusters];
395  _cluster_HCALIN_Eta = new float[_maxNclusters];
396  _cluster_HCALIN_Phi = new float[_maxNclusters];
397  _cluster_HCALIN_NTower = new int[_maxNclusters];
398  _cluster_HCALIN_trueID = new int[_maxNclusters];
399 
400  _tower_HCALOUT_E = new float[_maxNTowersCentral];
401  _tower_HCALOUT_iEta = new int[_maxNTowersCentral];
402  _tower_HCALOUT_iPhi = new int[_maxNTowersCentral];
403  _tower_HCALOUT_trueID = new int[_maxNTowersCentral];
404  _cluster_HCALOUT_E = new float[_maxNclusters];
405  _cluster_HCALOUT_Eta = new float[_maxNclusters];
406  _cluster_HCALOUT_Phi = new float[_maxNclusters];
407  _cluster_HCALOUT_NTower = new int[_maxNclusters];
408  _cluster_HCALOUT_trueID = new int[_maxNclusters];
409 
410  _tower_EHCAL_E = new float[_maxNTowers];
411  _tower_EHCAL_iEta = new int[_maxNTowers];
412  _tower_EHCAL_iPhi = new int[_maxNTowers];
413  _tower_EHCAL_trueID = new int[_maxNTowers];
414  _cluster_EHCAL_E = new float[_maxNclusters];
415  _cluster_EHCAL_Eta = new float[_maxNclusters];
416  _cluster_EHCAL_Phi = new float[_maxNclusters];
417  _cluster_EHCAL_NTower = new int[_maxNclusters];
418  _cluster_EHCAL_trueID = new int[_maxNclusters];
419 
420  _tower_DRCALO_E = new float[_maxNTowersDR];
421  _tower_DRCALO_NScint = new int[_maxNTowersDR];
422  _tower_DRCALO_NCerenkov = new int[_maxNTowersDR];
423  _tower_DRCALO_iEta = new int[_maxNTowersDR];
424  _tower_DRCALO_iPhi = new int[_maxNTowersDR];
425  _tower_DRCALO_trueID = new int[_maxNTowersDR];
426 
427  _tower_FOCAL_E = new float[_maxNTowersDR];
428  _tower_FOCAL_NScint = new int[_maxNTowersDR];
429  _tower_FOCAL_NCerenkov = new int[_maxNTowersDR];
430  _tower_FOCAL_iEta = new int[_maxNTowersDR];
431  _tower_FOCAL_iPhi = new int[_maxNTowersDR];
432  _tower_FOCAL_trueID = new int[_maxNTowersDR];
433 
434  _tower_LFHCAL_E = new float[_maxNTowers];
435  _tower_LFHCAL_iEta = new int[_maxNTowers];
436  _tower_LFHCAL_iPhi = new int[_maxNTowers];
437  _tower_LFHCAL_iL = new int[_maxNTowers];
438  _tower_LFHCAL_trueID = new int[_maxNTowers];
439 
440  _tower_FEMC_E = new float[_maxNTowers];
441  _tower_FEMC_iEta = new int[_maxNTowers];
442  _tower_FEMC_iPhi = new int[_maxNTowers];
443  _tower_FEMC_trueID = new int[_maxNTowers];
444  _cluster_FEMC_E = new float[_maxNclusters];
445  _cluster_FEMC_Eta = new float[_maxNclusters];
446  _cluster_FEMC_Phi = new float[_maxNclusters];
447  _cluster_FEMC_NTower = new int[_maxNclusters];
448  _cluster_FEMC_trueID = new int[_maxNclusters];
449 
450  _tower_CEMC_E = new float[_maxNTowersCentral];
451  _tower_CEMC_iEta = new int[_maxNTowersCentral];
452  _tower_CEMC_iPhi = new int[_maxNTowersCentral];
453  _tower_CEMC_trueID = new int[_maxNTowersCentral];
454  _cluster_CEMC_E = new float[_maxNclustersCentral];
455  _cluster_CEMC_Eta = new float[_maxNclustersCentral];
456  _cluster_CEMC_Phi = new float[_maxNclustersCentral];
457  _cluster_CEMC_NTower = new int[_maxNclustersCentral];
458  _cluster_CEMC_trueID = new int[_maxNclustersCentral];
459 
460  _tower_EEMC_E = new float[_maxNTowers];
461  _tower_EEMC_iEta = new int[_maxNTowers];
462  _tower_EEMC_iPhi = new int[_maxNTowers];
463  _tower_EEMC_trueID = new int[_maxNTowers];
464  _cluster_EEMC_E = new float[_maxNclusters];
465  _cluster_EEMC_Eta = new float[_maxNclusters];
466  _cluster_EEMC_Phi = new float[_maxNclusters];
467  _cluster_EEMC_NTower = new int[_maxNclusters];
468  _cluster_EEMC_trueID = new int[_maxNclusters];
469 
470  _tower_EEMCG_E = new float[_maxNTowers];
471  _tower_EEMCG_iEta = new int[_maxNTowers];
472  _tower_EEMCG_iPhi = new int[_maxNTowers];
473  _tower_EEMCG_trueID = new int[_maxNTowers];
474  _cluster_EEMCG_E = new float[_maxNclusters];
475  _cluster_EEMCG_Eta = new float[_maxNclusters];
476  _cluster_EEMCG_Phi = new float[_maxNclusters];
477  _cluster_EEMCG_NTower = new int[_maxNclusters];
478  _cluster_EEMCG_trueID = new int[_maxNclusters];
479 
480  _track_ID = new float[_maxNTracks];
481  _track_charge = new short[_maxNTracks];
482  _track_trueID = new float[_maxNTracks];
483  _track_px = new float[_maxNTracks];
484  _track_py = new float[_maxNTracks];
485  _track_pz = new float[_maxNTracks];
486  _track_x = new float[_maxNTracks];
487  _track_y = new float[_maxNTracks];
488  _track_z = new float[_maxNTracks];
489  _track_ndf = new float[_maxNTracks];
490  _track_chi2 = new float[_maxNTracks];
491  _track_dca = new float[_maxNTracks];
492  _track_dca_2d = new float[_maxNTracks];
493  _track_source = new unsigned short[_maxNTracks];
494  _track_ProjTrackID = new float[_maxNProjections];
495  _track_ProjLayer = new int[_maxNProjections];
496  _track_TLP_x = new float[_maxNProjections];
497  _track_TLP_y = new float[_maxNProjections];
498  _track_TLP_z = new float[_maxNProjections];
499  _track_TLP_t = new float[_maxNProjections];
500  _track_TLP_px = new float[_maxNProjections];
501  _track_TLP_py = new float[_maxNProjections];
502  _track_TLP_pz = new float[_maxNProjections];
503  _track_TLP_true_x = new float[_maxNProjections];
504  _track_TLP_true_y = new float[_maxNProjections];
505  _track_TLP_true_z = new float[_maxNProjections];
506  _track_TLP_true_t = new float[_maxNProjections];
507  _track_pion_LL.resize(_maxNTracks, -100);
508  _track_kaon_LL.resize(_maxNTracks, -100);
509  _track_proton_LL.resize(_maxNTracks, -100);
510 
511  _mcpart_ID = new int[_maxNMCPart];
512  _mcpart_ID_parent = new int[_maxNMCPart];
513  _mcpart_PDG = new int[_maxNMCPart];
514  _mcpart_E = new float[_maxNMCPart];
515  _mcpart_px = new float[_maxNMCPart];
516  _mcpart_py = new float[_maxNMCPart];
517  _mcpart_pz = new float[_maxNMCPart];
518  _mcpart_x = new float[_maxNMCPart];
519  _mcpart_y = new float[_maxNMCPart];
520  _mcpart_z = new float[_maxNMCPart];
521  _mcpart_BCID = new int[_maxNMCPart];
522 
523  _hepmcp_BCID = new int[_maxNHepmcp];
524  // _hepmcp_ID_parent = new float[_maxNHepmcp];
525  _hepmcp_status = new int[_maxNHepmcp];
526  _hepmcp_PDG = new int[_maxNHepmcp];
527  _hepmcp_E = new float[_maxNHepmcp];
528  _hepmcp_px = new float[_maxNHepmcp];
529  _hepmcp_py = new float[_maxNHepmcp];
530  _hepmcp_pz = new float[_maxNHepmcp];
531  _hepmcp_m1 = new int[_maxNHepmcp];
532  _hepmcp_m2 = new int[_maxNHepmcp];
533 
534  _calo_towers_iEta = new int[_maxNTowersCalo];
535  _calo_towers_iPhi = new int[_maxNTowersCalo];
536  _calo_towers_iL = new int[_maxNTowersCalo];
537  _calo_towers_Eta = new float[_maxNTowersCalo];
538  _calo_towers_Phi = new float[_maxNTowersCalo];
539  _calo_towers_x = new float[_maxNTowersCalo];
540  _calo_towers_y = new float[_maxNTowersCalo];
541  _calo_towers_z = new float[_maxNTowersCalo];
542  _geometry_done = new int[20];
543  for (int igem = 0; igem < 20; igem++) _geometry_done[igem] = 0;
544 }
545 
546 int EventEvaluatorEIC::Init(PHCompositeNode* topNode)
547 {
548  _ievent = 0;
549 
550  _tfile = new TFile(_filename.c_str(), "RECREATE");
551 
552  _event_tree = new TTree("event_tree", "event_tree");
553  if (_do_store_event_info)
554  {
555  // Event level info. This isn't the most efficient way to store this info, but it's straightforward
556  // within the structure of the class, so the size is small compared to the rest of the output.
557  _event_tree->Branch("cross_section", &_cross_section, "cross_section/F");
558  _event_tree->Branch("event_weight", &_event_weight, "event_weight/F");
559  _event_tree->Branch("n_generator_accepted", &_n_generator_accepted, "n_generator_accepted/I");
560  }
561  // tracks and hits
562  if (_do_HITS)
563  {
564  _event_tree->Branch("nHits", &_nHitsLayers, "nHits/I");
565  _event_tree->Branch("hits_layerID", _hits_layerID, "hits_layerID[nHits]/I");
566  _event_tree->Branch("hits_trueID", _hits_trueID, "hits_trueID[nHits]/I");
567  _event_tree->Branch("hits_x", _hits_x, "hits_x[nHits]/F");
568  _event_tree->Branch("hits_y", _hits_y, "hits_y[nHits]/F");
569  _event_tree->Branch("hits_z", _hits_z, "hits_z[nHits]/F");
570  _event_tree->Branch("hits_x2", _hits_x2, "hits_x2[nHits]/F");
571  _event_tree->Branch("hits_y2", _hits_y2, "hits_y2[nHits]/F");
572  _event_tree->Branch("hits_z2", _hits_z2, "hits_z2[nHits]/F");
573  _event_tree->Branch("hits_t", _hits_t, "hits_t[nHits]/F");
574  _event_tree->Branch("hits_edep", _hits_edep, "hits_edep[nHits]/F");
575  _event_tree->Branch("hits_lightyield", _hits_lightyield, "hits_lightyield[nHits]/F");
576  _event_tree->Branch("hits_isAbsorber", _hits_isAbsorber, "hits_isAbsorber[nHits]/I");
577  }
578 
579  if (_do_TRACKS)
580  {
581  _event_tree->Branch("nTracks", &_nTracks, "nTracks/I");
582  _event_tree->Branch("tracks_ID", _track_ID, "tracks_ID[nTracks]/F");
583  _event_tree->Branch("tracks_charge", _track_charge, "tracks_charge[nTracks]/S");
584  _event_tree->Branch("tracks_px", _track_px, "tracks_px[nTracks]/F");
585  _event_tree->Branch("tracks_py", _track_py, "tracks_py[nTracks]/F");
586  _event_tree->Branch("tracks_pz", _track_pz, "tracks_pz[nTracks]/F");
587  _event_tree->Branch("tracks_x", _track_x, "tracks_x[nTracks]/F");
588  _event_tree->Branch("tracks_y", _track_y, "tracks_y[nTracks]/F");
589  _event_tree->Branch("tracks_z", _track_z, "tracks_z[nTracks]/F");
590  _event_tree->Branch("tracks_ndf", _track_ndf, "tracks_ndf[nTracks]/F");
591  _event_tree->Branch("tracks_chi2", _track_chi2, "tracks_chi2[nTracks]/F");
592  _event_tree->Branch("tracks_dca", _track_dca, "tracks_dca[nTracks]/F");
593  _event_tree->Branch("tracks_dca_2d", _track_dca_2d, "tracks_dca_2d[nTracks]/F");
594  _event_tree->Branch("tracks_trueID", _track_trueID, "tracks_trueID[nTracks]/F");
595  _event_tree->Branch("tracks_source", _track_source, "tracks_source[nTracks]/s");
596 
597  if (_do_PID_LogLikelihood)
598  {
599  // save hadron PID log likelihood
600  _event_tree->Branch("track_pion_LL", _track_pion_LL.data(), "track_pion_LL[nTracks]/F");
601  _event_tree->Branch("track_kaon_LL", _track_kaon_LL.data(), "track_kaon_LL[nTracks]/F");
602  _event_tree->Branch("track_proton_LL", _track_proton_LL.data(), "track_proton_LL[nTracks]/F");
603  }
604  }
605  if (_do_PROJECTIONS)
606  {
607  _event_tree->Branch("nProjections", &_nProjections, "nProjections/I");
608  _event_tree->Branch("track_ProjTrackID", _track_ProjTrackID, "track_ProjTrackID[nProjections]/F");
609  _event_tree->Branch("track_ProjLayer", _track_ProjLayer, "track_ProjLayer[nProjections]/I");
610  _event_tree->Branch("track_TLP_x", _track_TLP_x, "track_TLP_x[nProjections]/F");
611  _event_tree->Branch("track_TLP_y", _track_TLP_y, "track_TLP_y[nProjections]/F");
612  _event_tree->Branch("track_TLP_z", _track_TLP_z, "track_TLP_z[nProjections]/F");
613  _event_tree->Branch("track_TLP_t", _track_TLP_t, "track_TLP_t[nProjections]/F");
614  _event_tree->Branch("track_TLP_px", _track_TLP_px, "track_TLP_px[nProjections]/F");
615  _event_tree->Branch("track_TLP_py", _track_TLP_py, "track_TLP_py[nProjections]/F");
616  _event_tree->Branch("track_TLP_pz", _track_TLP_pz, "track_TLP_pz[nProjections]/F");
617  _event_tree->Branch("track_TLP_true_x", _track_TLP_true_x, "track_TLP_true_x[nProjections]/F");
618  _event_tree->Branch("track_TLP_true_y", _track_TLP_true_y, "track_TLP_true_y[nProjections]/F");
619  _event_tree->Branch("track_TLP_true_z", _track_TLP_true_z, "track_TLP_true_z[nProjections]/F");
620  _event_tree->Branch("track_TLP_true_t", _track_TLP_true_t, "track_TLP_true_t[nProjections]/F");
621  }
622  if (_do_FHCAL)
623  {
624  // towers FHCAL
625  _event_tree->Branch("tower_FHCAL_N", &_nTowers_FHCAL, "tower_FHCAL_N/I");
626  _event_tree->Branch("tower_FHCAL_E", _tower_FHCAL_E, "tower_FHCAL_E[tower_FHCAL_N]/F");
627  _event_tree->Branch("tower_FHCAL_iEta", _tower_FHCAL_iEta, "tower_FHCAL_iEta[tower_FHCAL_N]/I");
628  _event_tree->Branch("tower_FHCAL_iPhi", _tower_FHCAL_iPhi, "tower_FHCAL_iPhi[tower_FHCAL_N]/I");
629  _event_tree->Branch("tower_FHCAL_trueID", _tower_FHCAL_trueID, "tower_FHCAL_trueID[tower_FHCAL_N]/I");
630  if (_do_CLUSTERS)
631  {
632  // clusters FHCAL
633  _event_tree->Branch("cluster_FHCAL_N", &_nclusters_FHCAL, "cluster_FHCAL_N/I");
634  _event_tree->Branch("cluster_FHCAL_E", _cluster_FHCAL_E, "cluster_FHCAL_E[cluster_FHCAL_N]/F");
635  _event_tree->Branch("cluster_FHCAL_Eta", _cluster_FHCAL_Eta, "cluster_FHCAL_Eta[cluster_FHCAL_N]/F");
636  _event_tree->Branch("cluster_FHCAL_Phi", _cluster_FHCAL_Phi, "cluster_FHCAL_Phi[cluster_FHCAL_N]/F");
637  _event_tree->Branch("cluster_FHCAL_NTower", _cluster_FHCAL_NTower, "cluster_FHCAL_NTower[cluster_FHCAL_N]/I");
638  _event_tree->Branch("cluster_FHCAL_trueID", _cluster_FHCAL_trueID, "cluster_FHCAL_trueID[cluster_FHCAL_N]/I");
639  }
640  }
641  if (_do_BECAL)
642  {
643  // towers BECAL
644  _event_tree->Branch("tower_BECAL_N", &_nTowers_BECAL, "tower_BECAL_N/I");
645  _event_tree->Branch("tower_BECAL_E", _tower_BECAL_E, "tower_BECAL_E[tower_BECAL_N]/F");
646  _event_tree->Branch("tower_BECAL_iEta", _tower_BECAL_iEta, "tower_BECAL_iEta[tower_BECAL_N]/I");
647  _event_tree->Branch("tower_BECAL_iPhi", _tower_BECAL_iPhi, "tower_BECAL_iPhi[tower_BECAL_N]/I");
648  _event_tree->Branch("tower_BECAL_trueID", _tower_BECAL_trueID, "tower_BECAL_trueID[tower_BECAL_N]/I");
649  }
650  if (_do_HCALIN)
651  {
652  // towers HCAL-in
653  _event_tree->Branch("tower_HCALIN_N", &_nTowers_HCALIN, "tower_HCALIN_N/I");
654  _event_tree->Branch("tower_HCALIN_E", _tower_HCALIN_E, "tower_HCALIN_E[tower_HCALIN_N]/F");
655  _event_tree->Branch("tower_HCALIN_iEta", _tower_HCALIN_iEta, "tower_HCALIN_iEta[tower_HCALIN_N]/I");
656  _event_tree->Branch("tower_HCALIN_iPhi", _tower_HCALIN_iPhi, "tower_HCALIN_iPhi[tower_HCALIN_N]/I");
657  _event_tree->Branch("tower_HCALIN_trueID", _tower_HCALIN_trueID, "tower_HCALIN_trueID[tower_HCALIN_N]/I");
658  if (_do_CLUSTERS)
659  {
660  // clusters HCAL-in
661  _event_tree->Branch("cluster_HCALIN_N", &_nclusters_HCALIN, "cluster_HCALIN_N/I");
662  _event_tree->Branch("cluster_HCALIN_E", _cluster_HCALIN_E, "cluster_HCALIN_E[cluster_HCALIN_N]/F");
663  _event_tree->Branch("cluster_HCALIN_Eta", _cluster_HCALIN_Eta, "cluster_HCALIN_Eta[cluster_HCALIN_N]/F");
664  _event_tree->Branch("cluster_HCALIN_Phi", _cluster_HCALIN_Phi, "cluster_HCALIN_Phi[cluster_HCALIN_N]/F");
665  _event_tree->Branch("cluster_HCALIN_NTower", _cluster_HCALIN_NTower, "cluster_HCALIN_NTower[cluster_HCALIN_N]/I");
666  _event_tree->Branch("cluster_HCALIN_trueID", _cluster_HCALIN_trueID, "cluster_HCALIN_trueID[cluster_HCALIN_N]/I");
667  }
668  }
669  if (_do_HCALOUT)
670  {
671  // towers HCAL-out
672  _event_tree->Branch("tower_HCALOUT_N", &_nTowers_HCALOUT, "tower_HCALOUT_N/I");
673  _event_tree->Branch("tower_HCALOUT_E", _tower_HCALOUT_E, "tower_HCALOUT_E[tower_HCALOUT_N]/F");
674  _event_tree->Branch("tower_HCALOUT_iEta", _tower_HCALOUT_iEta, "tower_HCALOUT_iEta[tower_HCALOUT_N]/I");
675  _event_tree->Branch("tower_HCALOUT_iPhi", _tower_HCALOUT_iPhi, "tower_HCALOUT_iPhi[tower_HCALOUT_N]/I");
676  _event_tree->Branch("tower_HCALOUT_trueID", _tower_HCALOUT_trueID, "tower_HCALOUT_trueID[tower_HCALOUT_N]/I");
677  if (_do_CLUSTERS)
678  {
679  // clusters HCAL-out
680  _event_tree->Branch("cluster_HCALOUT_N", &_nclusters_HCALOUT, "cluster_HCALOUT_N/I");
681  _event_tree->Branch("cluster_HCALOUT_E", _cluster_HCALOUT_E, "cluster_HCALOUT_E[cluster_HCALOUT_N]/F");
682  _event_tree->Branch("cluster_HCALOUT_Eta", _cluster_HCALOUT_Eta, "cluster_HCALOUT_Eta[cluster_HCALOUT_N]/F");
683  _event_tree->Branch("cluster_HCALOUT_Phi", _cluster_HCALOUT_Phi, "cluster_HCALOUT_Phi[cluster_HCALOUT_N]/F");
684  _event_tree->Branch("cluster_HCALOUT_NTower", _cluster_HCALOUT_NTower, "cluster_HCALOUT_NTower[cluster_HCALOUT_N]/I");
685  _event_tree->Branch("cluster_HCALOUT_trueID", _cluster_HCALOUT_trueID, "cluster_HCALOUT_trueID[cluster_HCALOUT_N]/I");
686  }
687  }
688  if (_do_EHCAL)
689  {
690  // towers EHCAL
691  _event_tree->Branch("tower_EHCAL_N", &_nTowers_EHCAL, "tower_EHCAL_N/I");
692  _event_tree->Branch("tower_EHCAL_E", _tower_EHCAL_E, "tower_EHCAL_E[tower_EHCAL_N]/F");
693  _event_tree->Branch("tower_EHCAL_iEta", _tower_EHCAL_iEta, "tower_EHCAL_iEta[tower_EHCAL_N]/I");
694  _event_tree->Branch("tower_EHCAL_iPhi", _tower_EHCAL_iPhi, "tower_EHCAL_iPhi[tower_EHCAL_N]/I");
695  _event_tree->Branch("tower_EHCAL_trueID", _tower_EHCAL_trueID, "tower_EHCAL_trueID[tower_EHCAL_N]/I");
696  if (_do_CLUSTERS)
697  {
698  // clusters EHCAL
699  _event_tree->Branch("cluster_EHCAL_N", &_nclusters_EHCAL, "cluster_EHCAL_N/I");
700  _event_tree->Branch("cluster_EHCAL_E", _cluster_EHCAL_E, "cluster_EHCAL_E[cluster_EHCAL_N]/F");
701  _event_tree->Branch("cluster_EHCAL_Eta", _cluster_EHCAL_Eta, "cluster_EHCAL_Eta[cluster_EHCAL_N]/F");
702  _event_tree->Branch("cluster_EHCAL_Phi", _cluster_EHCAL_Phi, "cluster_EHCAL_Phi[cluster_EHCAL_N]/F");
703  _event_tree->Branch("cluster_EHCAL_NTower", _cluster_EHCAL_NTower, "cluster_EHCAL_NTower[cluster_EHCAL_N]/I");
704  _event_tree->Branch("cluster_EHCAL_trueID", _cluster_EHCAL_trueID, "cluster_EHCAL_trueID[cluster_EHCAL_N]/I");
705  }
706  }
707  if (_do_DRCALO)
708  {
709  // towers DRCALO
710  _event_tree->Branch("tower_DRCALO_N", &_nTowers_DRCALO, "tower_DRCALO_N/I");
711  _event_tree->Branch("tower_DRCALO_E", _tower_DRCALO_E, "tower_DRCALO_E[tower_DRCALO_N]/F");
712  _event_tree->Branch("tower_DRCALO_NScint", _tower_DRCALO_NScint, "tower_DRCALO_NScint[tower_DRCALO_N]/I");
713  _event_tree->Branch("tower_DRCALO_NCerenkov", _tower_DRCALO_NCerenkov, "tower_DRCALO_NCerenkov[tower_DRCALO_N]/I");
714  _event_tree->Branch("tower_DRCALO_iEta", _tower_DRCALO_iEta, "tower_DRCALO_iEta[tower_DRCALO_N]/I");
715  _event_tree->Branch("tower_DRCALO_iPhi", _tower_DRCALO_iPhi, "tower_DRCALO_iPhi[tower_DRCALO_N]/I");
716  _event_tree->Branch("tower_DRCALO_trueID", _tower_DRCALO_trueID, "tower_DRCALO_trueID[tower_DRCALO_N]/I");
717  }
718  if (_do_FOCAL)
719  {
720  // towers FOCAL
721  _event_tree->Branch("tower_FOCAL_N", &_nTowers_FOCAL, "tower_FOCAL_N/I");
722  _event_tree->Branch("tower_FOCAL_E", _tower_FOCAL_E, "tower_FOCAL_E[tower_FOCAL_N]/F");
723  _event_tree->Branch("tower_FOCAL_NScint", _tower_FOCAL_NScint, "tower_FOCAL_NScint[tower_FOCAL_N]/I");
724  _event_tree->Branch("tower_FOCAL_NCerenkov", _tower_FOCAL_NCerenkov, "tower_FOCAL_NCerenkov[tower_FOCAL_N]/I");
725  _event_tree->Branch("tower_FOCAL_iEta", _tower_FOCAL_iEta, "tower_FOCAL_iEta[tower_FOCAL_N]/I");
726  _event_tree->Branch("tower_FOCAL_iPhi", _tower_FOCAL_iPhi, "tower_FOCAL_iPhi[tower_FOCAL_N]/I");
727  _event_tree->Branch("tower_FOCAL_trueID", _tower_FOCAL_trueID, "tower_FOCAL_trueID[tower_FOCAL_N]/I");
728  }
729  if (_do_LFHCAL)
730  {
731  // towers LFHCALO
732  _event_tree->Branch("tower_LFHCAL_N", &_nTowers_LFHCAL, "tower_LFHCAL_N/I");
733  _event_tree->Branch("tower_LFHCAL_E", _tower_LFHCAL_E, "tower_LFHCAL_E[tower_LFHCAL_N]/F");
734  _event_tree->Branch("tower_LFHCAL_iEta", _tower_LFHCAL_iEta, "tower_LFHCAL_iEta[tower_LFHCAL_N]/I");
735  _event_tree->Branch("tower_LFHCAL_iPhi", _tower_LFHCAL_iPhi, "tower_LFHCAL_iPhi[tower_LFHCAL_N]/I");
736  _event_tree->Branch("tower_LFHCAL_iL", _tower_LFHCAL_iL, "tower_LFHCAL_iL[tower_LFHCAL_N]/I");
737  _event_tree->Branch("tower_LFHCAL_trueID", _tower_LFHCAL_trueID, "tower_LFHCAL_trueID[tower_LFHCAL_N]/I");
738  }
739  if (_do_FEMC)
740  {
741  // towers FEMC
742  _event_tree->Branch("tower_FEMC_N", &_nTowers_FEMC, "tower_FEMC_N/I");
743  _event_tree->Branch("tower_FEMC_E", _tower_FEMC_E, "tower_FEMC_E[tower_FEMC_N]/F");
744  _event_tree->Branch("tower_FEMC_iEta", _tower_FEMC_iEta, "tower_FEMC_iEta[tower_FEMC_N]/I");
745  _event_tree->Branch("tower_FEMC_iPhi", _tower_FEMC_iPhi, "tower_FEMC_iPhi[tower_FEMC_N]/I");
746  _event_tree->Branch("tower_FEMC_trueID", _tower_FEMC_trueID, "tower_FEMC_trueID[tower_FEMC_N]/I");
747  if (_do_CLUSTERS)
748  {
749  // clusters FEMC
750  _event_tree->Branch("cluster_FEMC_N", &_nclusters_FEMC, "cluster_FEMC_N/I");
751  _event_tree->Branch("cluster_FEMC_E", _cluster_FEMC_E, "cluster_FEMC_E[cluster_FEMC_N]/F");
752  _event_tree->Branch("cluster_FEMC_Eta", _cluster_FEMC_Eta, "cluster_FEMC_Eta[cluster_FEMC_N]/F");
753  _event_tree->Branch("cluster_FEMC_Phi", _cluster_FEMC_Phi, "cluster_FEMC_Phi[cluster_FEMC_N]/F");
754  _event_tree->Branch("cluster_FEMC_NTower", _cluster_FEMC_NTower, "cluster_FEMC_NTower[cluster_FEMC_N]/I");
755  _event_tree->Branch("cluster_FEMC_trueID", _cluster_FEMC_trueID, "cluster_FEMC_trueID[cluster_FEMC_N]/I");
756  }
757  }
758  if (_do_CEMC)
759  {
760  // towers CEMC
761  _event_tree->Branch("tower_CEMC_N", &_nTowers_CEMC, "tower_CEMC_N/I");
762  _event_tree->Branch("tower_CEMC_E", _tower_CEMC_E, "tower_CEMC_E[tower_CEMC_N]/F");
763  _event_tree->Branch("tower_CEMC_iEta", _tower_CEMC_iEta, "tower_CEMC_iEta[tower_CEMC_N]/I");
764  _event_tree->Branch("tower_CEMC_iPhi", _tower_CEMC_iPhi, "tower_CEMC_iPhi[tower_CEMC_N]/I");
765  _event_tree->Branch("tower_CEMC_trueID", _tower_CEMC_trueID, "tower_CEMC_trueID[tower_CEMC_N]/I");
766  if (_do_CLUSTERS)
767  {
768  // clusters CEMC
769  _event_tree->Branch("cluster_CEMC_N", &_nclusters_CEMC, "cluster_CEMC_N/I");
770  _event_tree->Branch("cluster_CEMC_E", _cluster_CEMC_E, "cluster_CEMC_E[cluster_CEMC_N]/F");
771  _event_tree->Branch("cluster_CEMC_Eta", _cluster_CEMC_Eta, "cluster_CEMC_Eta[cluster_CEMC_N]/F");
772  _event_tree->Branch("cluster_CEMC_Phi", _cluster_CEMC_Phi, "cluster_CEMC_Phi[cluster_CEMC_N]/F");
773  _event_tree->Branch("cluster_CEMC_NTower", _cluster_CEMC_NTower, "cluster_CEMC_NTower[cluster_CEMC_N]/I");
774  _event_tree->Branch("cluster_CEMC_trueID", _cluster_CEMC_trueID, "cluster_CEMC_trueID[cluster_CEMC_N]/I");
775  }
776  }
777  if (_do_EEMC)
778  {
779  // towers EEMC
780  _event_tree->Branch("tower_EEMC_N", &_nTowers_EEMC, "tower_EEMC_N/I");
781  _event_tree->Branch("tower_EEMC_E", _tower_EEMC_E, "tower_EEMC_E[tower_EEMC_N]/F");
782  _event_tree->Branch("tower_EEMC_iEta", _tower_EEMC_iEta, "tower_EEMC_iEta[tower_EEMC_N]/I");
783  _event_tree->Branch("tower_EEMC_iPhi", _tower_EEMC_iPhi, "tower_EEMC_iPhi[tower_EEMC_N]/I");
784  _event_tree->Branch("tower_EEMC_trueID", _tower_EEMC_trueID, "tower_EEMC_trueID[tower_EEMC_N]/I");
785  if (_do_CLUSTERS)
786  {
787  // clusters EEMC
788  _event_tree->Branch("cluster_EEMC_N", &_nclusters_EEMC, "cluster_EEMC_N/I");
789  _event_tree->Branch("cluster_EEMC_E", _cluster_EEMC_E, "cluster_EEMC_E[cluster_EEMC_N]/F");
790  _event_tree->Branch("cluster_EEMC_Eta", _cluster_EEMC_Eta, "cluster_EEMC_Eta[cluster_EEMC_N]/F");
791  _event_tree->Branch("cluster_EEMC_Phi", _cluster_EEMC_Phi, "cluster_EEMC_Phi[cluster_EEMC_N]/F");
792  _event_tree->Branch("cluster_EEMC_NTower", _cluster_EEMC_NTower, "cluster_EEMC_NTower[cluster_EEMC_N]/I");
793  _event_tree->Branch("cluster_EEMC_trueID", _cluster_EEMC_trueID, "cluster_EEMC_trueID[cluster_EEMC_N]/I");
794  }
795  }
796  if (_do_EEMCG)
797  {
798  // towers EEMCG
799  _event_tree->Branch("tower_EEMCG_N", &_nTowers_EEMCG, "tower_EEMCG_N/I");
800  _event_tree->Branch("tower_EEMCG_E", _tower_EEMCG_E, "tower_EEMCG_E[tower_EEMCG_N]/F");
801  _event_tree->Branch("tower_EEMCG_iEta", _tower_EEMCG_iEta, "tower_EEMCG_iEta[tower_EEMCG_N]/I");
802  _event_tree->Branch("tower_EEMCG_iPhi", _tower_EEMCG_iPhi, "tower_EEMCG_iPhi[tower_EEMCG_N]/I");
803  _event_tree->Branch("tower_EEMCG_trueID", _tower_EEMCG_trueID, "tower_EEMCG_trueID[tower_EEMCG_N]/I");
804  if (_do_CLUSTERS)
805  {
806  // clusters EEMCG
807  _event_tree->Branch("cluster_EEMCG_N", &_nclusters_EEMCG, "cluster_EEMCG_N/I");
808  _event_tree->Branch("cluster_EEMCG_E", _cluster_EEMCG_E, "cluster_EEMCG_E[cluster_EEMCG_N]/F");
809  _event_tree->Branch("cluster_EEMCG_Eta", _cluster_EEMCG_Eta, "cluster_EEMCG_Eta[cluster_EEMCG_N]/F");
810  _event_tree->Branch("cluster_EEMCG_Phi", _cluster_EEMCG_Phi, "cluster_EEMCG_Phi[cluster_EEMCG_N]/F");
811  _event_tree->Branch("cluster_EEMCG_NTower", _cluster_EEMCG_NTower, "cluster_EEMCG_NTower[cluster_EEMCG_N]/I");
812  _event_tree->Branch("cluster_EEMCG_trueID", _cluster_EEMCG_trueID, "cluster_EEMCG_trueID[cluster_EEMCG_N]/I");
813  }
814  }
815  if (_do_VERTEX)
816  {
817  // vertex
818  _event_tree->Branch("vertex_x", &_vertex_x, "vertex_x/F");
819  _event_tree->Branch("vertex_y", &_vertex_y, "vertex_y/F");
820  _event_tree->Branch("vertex_z", &_vertex_z, "vertex_z/F");
821  _event_tree->Branch("vertex_NCont", &_vertex_NCont, "vertex_NCont/I");
822  _event_tree->Branch("vertex_true_x", &_vertex_true_x, "vertex_true_x/F");
823  _event_tree->Branch("vertex_true_y", &_vertex_true_y, "vertex_true_y/F");
824  _event_tree->Branch("vertex_true_z", &_vertex_true_z, "vertex_true_z/F");
825  }
826  if (_do_MCPARTICLES)
827  {
828  // MC particles
829  _event_tree->Branch("nMCPart", &_nMCPart, "nMCPart/I");
830  _event_tree->Branch("mcpart_ID", _mcpart_ID, "mcpart_ID[nMCPart]/I");
831  _event_tree->Branch("mcpart_ID_parent", _mcpart_ID_parent, "mcpart_ID_parent[nMCPart]/I");
832  _event_tree->Branch("mcpart_PDG", _mcpart_PDG, "mcpart_PDG[nMCPart]/I");
833  _event_tree->Branch("mcpart_E", _mcpart_E, "mcpart_E[nMCPart]/F");
834  _event_tree->Branch("mcpart_px", _mcpart_px, "mcpart_px[nMCPart]/F");
835  _event_tree->Branch("mcpart_py", _mcpart_py, "mcpart_py[nMCPart]/F");
836  _event_tree->Branch("mcpart_pz", _mcpart_pz, "mcpart_pz[nMCPart]/F");
837  _event_tree->Branch("mcpart_x", _mcpart_x, "mcpart_x[nMCPart]/F");
838  _event_tree->Branch("mcpart_y", _mcpart_y, "mcpart_y[nMCPart]/F");
839  _event_tree->Branch("mcpart_z", _mcpart_z, "mcpart_z[nMCPart]/F");
840  _event_tree->Branch("mcpart_BCID", _mcpart_BCID, "mcpart_BCID[nMCPart]/I");
841  }
842  if (_do_HEPMC)
843  {
844  // MC particles
845  _event_tree->Branch("nHepmcp", &_nHepmcp, "nHepmcp/I");
846  _event_tree->Branch("hepmcp_procid", &_hepmcp_procid, "hepmcp_procid/I");
847  _event_tree->Branch("hepmcp_x1", &_hepmcp_x1, "hepmcp_x1/F");
848  _event_tree->Branch("hepmcp_x2", &_hepmcp_x2, "hepmcp_x2/F");
849  _event_tree->Branch("hepmcp_Q2", &_hepmcp_Q2, "hepmcp_Q2/F");
850  _event_tree->Branch("hepmcp_vtx_x", &_hepmcp_vtx_x, "hepmcp_vtx_x/F");
851  _event_tree->Branch("hepmcp_vtx_y", &_hepmcp_vtx_y, "hepmcp_vtx_y/F");
852  _event_tree->Branch("hepmcp_vtx_z", &_hepmcp_vtx_z, "hepmcp_vtx_z/F");
853  _event_tree->Branch("hepmcp_vtx_t", &_hepmcp_vtx_t, "hepmcp_vtx_t/F");
854 
855  // _event_tree->Branch("hepmcp_ID_parent", _hepmcp_ID_parent, "hepmcp_ID_parent[nHepmcp]/F");
856  _event_tree->Branch("hepmcp_status", _hepmcp_status, "hepmcp_status[nHepmcp]/I");
857  _event_tree->Branch("hepmcp_PDG", _hepmcp_PDG, "hepmcp_PDG[nHepmcp]/I");
858  _event_tree->Branch("hepmcp_E", _hepmcp_E, "hepmcp_E[nHepmcp]/F");
859  _event_tree->Branch("hepmcp_px", _hepmcp_px, "hepmcp_px[nHepmcp]/F");
860  _event_tree->Branch("hepmcp_py", _hepmcp_py, "hepmcp_py[nHepmcp]/F");
861  _event_tree->Branch("hepmcp_pz", _hepmcp_pz, "hepmcp_pz[nHepmcp]/F");
862  _event_tree->Branch("hepmcp_BCID", _hepmcp_BCID, "hepmcp_BCID[nHepmcp]/I");
863  _event_tree->Branch("hepmcp_m1", _hepmcp_m1, "hepmcp_m1[nHepmcp]/I");
864  _event_tree->Branch("hepmcp_m2", _hepmcp_m2, "hepmcp_m2[nHepmcp]/I");
865  }
866 
867  if (_do_GEOMETRY)
868  {
869  _tfile_geometry = new TFile("geometry.root", "RECREATE");
870 
871  _geometry_tree = new TTree("geometry_tree", "geometry_tree");
872  // tracks and hits
873  _geometry_tree->Branch("calo", &_calo_ID, "nHits/I");
874  _geometry_tree->Branch("calo_towers_N", &_calo_towers_N, "calo_towers_N/I");
875  _geometry_tree->Branch("calo_towers_iEta", _calo_towers_iEta, "calo_towers_iEta[calo_towers_N]/I");
876  _geometry_tree->Branch("calo_towers_iPhi", _calo_towers_iPhi, "calo_towers_iPhi[calo_towers_N]/I");
877  _geometry_tree->Branch("calo_towers_iL", _calo_towers_iL, "calo_towers_iL[calo_towers_N]/I");
878  _geometry_tree->Branch("calo_towers_Eta", _calo_towers_Eta, "calo_towers_Eta[calo_towers_N]/F");
879  _geometry_tree->Branch("calo_towers_Phi", _calo_towers_Phi, "calo_towers_Phi[calo_towers_N]/F");
880  _geometry_tree->Branch("calo_towers_x", _calo_towers_x, "calo_towers_x[calo_towers_N]/F");
881  _geometry_tree->Branch("calo_towers_y", _calo_towers_y, "calo_towers_y[calo_towers_N]/F");
882  _geometry_tree->Branch("calo_towers_z", _calo_towers_z, "calo_towers_z[calo_towers_N]/F");
883  }
884 
885  return Fun4AllReturnCodes::EVENT_OK;
886 }
887 
888 int EventEvaluatorEIC::process_event(PHCompositeNode* topNode)
889 {
890  if (Verbosity() > 0)
891  {
892  cout << "entered process_event" << endl;
893  }
894  if (_do_FHCAL)
895  {
896  if (!_caloevalstackFHCAL)
897  {
898  _caloevalstackFHCAL = new CaloEvalStack(topNode, "FHCAL");
899  _caloevalstackFHCAL->set_strict(_strict);
900  _caloevalstackFHCAL->set_verbosity(Verbosity() + 1);
901  }
902  else
903  {
904  _caloevalstackFHCAL->next_event(topNode);
905  }
906  }
907  if (_do_BECAL)
908  {
909  if (!_caloevalstackBECAL)
910  {
911  _caloevalstackBECAL = new CaloEvalStack(topNode, "BECAL");
912  _caloevalstackBECAL->set_strict(_strict);
913  _caloevalstackBECAL->set_verbosity(Verbosity() + 1);
914  }
915  else
916  {
917  _caloevalstackBECAL->next_event(topNode);
918  }
919  }
920  if (_do_HCALIN)
921  {
922  if (!_caloevalstackHCALIN)
923  {
924  _caloevalstackHCALIN = new CaloEvalStack(topNode, "HCALIN");
925  _caloevalstackHCALIN->set_strict(_strict);
926  _caloevalstackHCALIN->set_verbosity(Verbosity() + 1);
927  }
928  else
929  {
930  _caloevalstackHCALIN->next_event(topNode);
931  }
932  }
933  if (_do_HCALOUT)
934  {
935  if (!_caloevalstackHCALOUT)
936  {
937  _caloevalstackHCALOUT = new CaloEvalStack(topNode, "HCALOUT");
938  _caloevalstackHCALOUT->set_strict(_strict);
939  _caloevalstackHCALOUT->set_verbosity(Verbosity() + 1);
940  }
941  else
942  {
943  _caloevalstackHCALOUT->next_event(topNode);
944  }
945  }
946  if (_do_EHCAL)
947  {
948  if (!_caloevalstackEHCAL)
949  {
950  _caloevalstackEHCAL = new CaloEvalStack(topNode, "EHCAL");
951  _caloevalstackEHCAL->set_strict(_strict);
952  _caloevalstackEHCAL->set_verbosity(Verbosity() + 1);
953  }
954  else
955  {
956  _caloevalstackEHCAL->next_event(topNode);
957  }
958  }
959  if (_do_DRCALO)
960  {
961  if (!_caloevalstackDRCALO)
962  {
963  _caloevalstackDRCALO = new CaloEvalStack(topNode, "DRCALO");
964  _caloevalstackDRCALO->set_strict(_strict);
965  _caloevalstackDRCALO->set_verbosity(Verbosity() + 1);
966  }
967  else
968  {
969  _caloevalstackDRCALO->next_event(topNode);
970  }
971  }
972  if (_do_FOCAL)
973  {
974  if (!_caloevalstackFOCAL)
975  {
976  _caloevalstackFOCAL = new CaloEvalStack(topNode, "FOCAL");
977  _caloevalstackFOCAL->set_strict(_strict);
978  _caloevalstackFOCAL->set_verbosity(Verbosity() + 1);
979  }
980  else
981  {
982  _caloevalstackFOCAL->next_event(topNode);
983  }
984  }
985  if (_do_LFHCAL)
986  {
987  if (!_caloevalstackLFHCAL)
988  {
989  _caloevalstackLFHCAL = new CaloEvalStack(topNode, "LFHCAL");
990  _caloevalstackLFHCAL->set_strict(_strict);
991  _caloevalstackLFHCAL->set_verbosity(Verbosity() + 1);
992  }
993  else
994  {
995  _caloevalstackLFHCAL->next_event(topNode);
996  }
997  }
998 
999  if (_do_FEMC)
1000  {
1001  if (!_caloevalstackFEMC)
1002  {
1003  _caloevalstackFEMC = new CaloEvalStack(topNode, "FEMC");
1004  _caloevalstackFEMC->set_strict(_strict);
1005  _caloevalstackFEMC->set_verbosity(Verbosity() + 1);
1006  }
1007  else
1008  {
1009  _caloevalstackFEMC->next_event(topNode);
1010  }
1011  }
1012  if (_do_CEMC)
1013  {
1014  if (!_caloevalstackCEMC)
1015  {
1016  _caloevalstackCEMC = new CaloEvalStack(topNode, "CEMC");
1017  _caloevalstackCEMC->set_strict(_strict);
1018  _caloevalstackCEMC->set_verbosity(Verbosity() + 1);
1019  }
1020  else
1021  {
1022  _caloevalstackCEMC->next_event(topNode);
1023  }
1024  }
1025  if (_do_EEMC)
1026  {
1027  if (!_caloevalstackEEMC)
1028  {
1029  _caloevalstackEEMC = new CaloEvalStack(topNode, "EEMC");
1030  _caloevalstackEEMC->set_strict(_strict);
1031  _caloevalstackEEMC->set_verbosity(Verbosity() + 1);
1032  }
1033  else
1034  {
1035  _caloevalstackEEMC->next_event(topNode);
1036  }
1037  }
1038 
1039  if (_do_EEMCG)
1040  {
1041  if (!_caloevalstackEEMCG)
1042  {
1043  _caloevalstackEEMCG = new CaloEvalStack(topNode, "EEMC_glass");
1044  _caloevalstackEEMCG->set_strict(_strict);
1045  _caloevalstackEEMCG->set_verbosity(Verbosity() + 1);
1046  }
1047  else
1048  {
1049  _caloevalstackEEMCG->next_event(topNode);
1050  }
1051  }
1052 
1053  if (Verbosity() > 0)
1054  {
1055  cout << "loaded evalstack" << endl;
1056  }
1057 
1058  // fill the Evaluator Tree
1059  fillOutputNtuples(topNode);
1060 
1061  ++_ievent;
1062 
1063  return Fun4AllReturnCodes::EVENT_OK;
1064 }
1065 
1066 void EventEvaluatorEIC::fillOutputNtuples(PHCompositeNode* topNode)
1067 {
1068  if (Verbosity() > 2)
1069  {
1070  cout << "EventEvaluatorEIC::fillOutputNtuples() entered" << endl;
1071  }
1072 
1073  //----------------------
1074  // fill the Event Tree
1075  //----------------------
1076 
1077  //----------------------
1078  // Event level info
1079  //---------------------
1080  // Extract weight info from the stored HepMC event.
1081  if (_do_store_event_info)
1082  {
1083  PHHepMCGenEventMap* hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
1084  if (hepmceventmap)
1085  {
1086  if (Verbosity() > 0)
1087  {
1088  cout << "saving event level info" << endl;
1089  }
1090 
1091  for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
1092  eventIter != hepmceventmap->end();
1093  ++eventIter)
1094  {
1095  PHHepMCGenEvent* hepmcevent = eventIter->second;
1096 
1097  if (hepmcevent)
1098  {
1099  HepMC::GenEvent* truthevent = hepmcevent->getEvent();
1100  if (!truthevent)
1101  {
1102  cout << PHWHERE
1103  << "no evt pointer under phhepmvgeneventmap found "
1104  << endl;
1105  return;
1106  }
1107 
1108  auto xsec = truthevent->cross_section();
1109  if (xsec)
1110  {
1111  _cross_section = xsec->cross_section();
1112  }
1113  // Only fill the event weight if available.
1114  // The overall event weight will be stored in the last entry in the vector.
1115  auto weights = truthevent->weights();
1116  if (weights.size() > 0)
1117  {
1118  _event_weight = weights[weights.size() - 1];
1119  }
1120  }
1121  }
1122  }
1123  else
1124  {
1125  if (Verbosity() > 0)
1126  {
1127  cout << PHWHERE << " PHHepMCGenEventMap node (for event level info) not found on node tree" << endl;
1128  }
1129  return;
1130  }
1131 
1132  // Retrieve the number of generator accepted events
1133  // Following how this was implemented in PHPythia8
1134  PHNodeIterator iter(topNode);
1135  PHCompositeNode* sumNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
1136  if (!sumNode)
1137  {
1138  cout << PHWHERE << "RUN Node missing doing nothing" << endl;
1139  return;
1140  }
1141  auto* integralNode = findNode::getClass<PHGenIntegral>(sumNode, "PHGenIntegral");
1142  if (integralNode)
1143  {
1144  _n_generator_accepted = integralNode->get_N_Generator_Accepted_Event();
1145  }
1146  else
1147  {
1148  if (Verbosity() > 0)
1149  {
1150  cout << PHWHERE << " PHGenIntegral node (for n generator accepted) not found on node tree. Continuing" << endl;
1151  }
1152  }
1153  }
1154  //----------------------
1155  // VERTEX
1156  //----------------------
1157  SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
1158  if (_do_VERTEX)
1159  {
1160  if (vertexmap)
1161  {
1162  if (!vertexmap->empty())
1163  {
1164  if (Verbosity() > 0)
1165  {
1166  cout << "saving vertex" << endl;
1167  }
1168  SvtxVertex* vertex = (vertexmap->begin())->second;
1169 
1170  _vertex_x = vertex->get_x();
1171  _vertex_y = vertex->get_y();
1172  _vertex_z = vertex->get_z();
1173  _vertex_NCont = vertex->size_tracks();
1174  }
1175  else
1176  {
1177  _vertex_x = 0.;
1178  _vertex_y = 0.;
1179  _vertex_z = 0.;
1180  _vertex_NCont = -1;
1181  }
1182  }
1183  }
1184  //----------------------
1185  // HITS
1186  //----------------------
1187  if (_do_HITS)
1188  {
1189  if (Verbosity() > 0)
1190  {
1191  cout << "saving hits" << endl;
1192  }
1193  _nHitsLayers = 0;
1194  PHG4TruthInfoContainer* truthinfocontainerHits = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1195  for (int iIndex = 0; iIndex < 200; ++iIndex)
1196  {
1197  // you need to add your layer name here to be saved! This has to be done
1198  // as we do not want to save thousands of calorimeter hits!
1199  if (
1200  (_do_HITS_CALO && ( (GetProjectionNameFromIndex(iIndex).find("LFHCAL") != std::string::npos) ||
1201  (GetProjectionNameFromIndex(iIndex).find("HCALIN") != std::string::npos) ||
1202  (GetProjectionNameFromIndex(iIndex).find("HCALOUT") != std::string::npos) ||
1203  (GetProjectionNameFromIndex(iIndex).find("BECAL") != std::string::npos) ||
1204  (GetProjectionNameFromIndex(iIndex).find("FEMC") != std::string::npos) ||
1205  (GetProjectionNameFromIndex(iIndex).find("EHCAL") != std::string::npos) ||
1206  (GetProjectionNameFromIndex(iIndex).find("EEMCH") != std::string::npos) ||
1207  (GetProjectionNameFromIndex(iIndex).find("EEMC") != std::string::npos))) ||
1208  (GetProjectionNameFromIndex(iIndex).find("TTL") != std::string::npos) ||
1209  (GetProjectionNameFromIndex(iIndex).find("LBLVTX") != std::string::npos) ||
1210  (GetProjectionNameFromIndex(iIndex).find("BARR") != std::string::npos) ||
1211  (GetProjectionNameFromIndex(iIndex).find("FST") != std::string::npos) ||
1212  (GetProjectionNameFromIndex(iIndex).find("ZDCsurrogate") != std::string::npos) ||
1213  (GetProjectionNameFromIndex(iIndex).find("rpTruth") != std::string::npos) ||
1214  (GetProjectionNameFromIndex(iIndex).find("rpTruth2") != std::string::npos) || // needed for IP8
1215  (GetProjectionNameFromIndex(iIndex).find("offMomTruth") != std::string::npos) ||
1216  (GetProjectionNameFromIndex(iIndex).find("b0Truth") != std::string::npos) ||
1217  (GetProjectionNameFromIndex(iIndex).find("EGEM") != std::string::npos) ||
1218  (GetProjectionNameFromIndex(iIndex).find("FGEM") != std::string::npos) ||
1219  (GetProjectionNameFromIndex(iIndex).find("RWELL") != std::string::npos) ||
1220  (((GetProjectionNameFromIndex(iIndex).find("BH_1") != std::string::npos) || (GetProjectionNameFromIndex(iIndex).find("BH_FORWARD_PLUS") != std::string::npos) || (GetProjectionNameFromIndex(iIndex).find("BH_FORWARD_NEG") != std::string::npos)) && _do_BLACKHOLE))
1221  {
1222 
1223  string nodename = "G4HIT_" + GetProjectionNameFromIndex(iIndex);
1224  PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
1225  if (hits)
1226  {
1227  if (Verbosity() > 1)
1228  {
1229  cout << __PRETTY_FUNCTION__ << " number of hits: " << hits->size() << endl;
1230  }
1231  PHG4HitContainer::ConstRange hit_range = hits->getHits();
1232  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
1233  {
1234  if (Verbosity() > 1)
1235  {
1236  cout << __PRETTY_FUNCTION__ << " found hit with id " << hit_iter->second->get_trkid() << endl;
1237  }
1238  if (_nHitsLayers > _maxNHits)
1239  {
1240  cout << __PRETTY_FUNCTION__ << " exceededed maximum hit array size! Please check where these hits come from!" << endl;
1241  break;
1242  }
1243 
1244  _hits_x[_nHitsLayers] = hit_iter->second->get_x(0);
1245  _hits_y[_nHitsLayers] = hit_iter->second->get_y(0);
1246  _hits_z[_nHitsLayers] = hit_iter->second->get_z(0);
1247  _hits_x2[_nHitsLayers] = hit_iter->second->get_x(1);
1248  _hits_y2[_nHitsLayers] = hit_iter->second->get_y(1);
1249  _hits_z2[_nHitsLayers] = hit_iter->second->get_z(1);
1250  _hits_t[_nHitsLayers] = hit_iter->second->get_t(0);
1251  _hits_edep[_nHitsLayers] = hit_iter->second->get_edep();
1252  _hits_lightyield[_nHitsLayers] = hit_iter->second->get_light_yield();
1253  _hits_isAbsorber[_nHitsLayers] = 0;
1254  _hits_layerID[_nHitsLayers] = iIndex;
1255  // cout << "i " << hit_iter->second->get_index_i() << "\tj " <<hit_iter->second->get_index_j() << "\tk " <<hit_iter->second->get_index_k() << "\tl " << hit_iter->second->get_index_l() << "\tsens_x "<< hit_iter->second->get_strip_z_index()<< "\tsens_y "<< hit_iter->second->get_strip_y_index() << endl;
1256  if (truthinfocontainerHits)
1257  {
1258  PHG4Particle* particle = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
1259 
1260  if (particle->get_parent_id() != 0)
1261  {
1262  PHG4Particle* g4particleMother = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
1263  int mcSteps = 0;
1264  while (g4particleMother->get_parent_id() != 0)
1265  {
1266  g4particleMother = truthinfocontainerHits->GetParticle(g4particleMother->get_parent_id());
1267  if (g4particleMother == NULL) break;
1268  mcSteps += 1;
1269  }
1270  if (mcSteps <= _depth_MCstack)
1271  {
1272  _hits_trueID[_nHitsLayers] = hit_iter->second->get_trkid();
1273  }
1274  else
1275  {
1276  PHG4Particle* g4particleMother2 = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
1277  int mcSteps2 = 0;
1278  while (g4particleMother2->get_parent_id() != 0 && (mcSteps2 < (mcSteps - _depth_MCstack + 1)))
1279  {
1280  g4particleMother2 = truthinfocontainerHits->GetParticle(g4particleMother2->get_parent_id());
1281  if (g4particleMother2 == NULL)
1282  {
1283  break;
1284  }
1285  else
1286  {
1287  _hits_trueID[_nHitsLayers] = g4particleMother2->get_parent_id();
1288  mcSteps2 += 1;
1289  }
1290  }
1291  }
1292  }
1293  else
1294  {
1295  _hits_trueID[_nHitsLayers] = hit_iter->second->get_trkid();
1296  }
1297  }
1298  _nHitsLayers++;
1299  }
1300  if (Verbosity() > 0)
1301  {
1302  cout << "saved\t" << _nHitsLayers << "\thits for " << GetProjectionNameFromIndex(iIndex) << endl;
1303  }
1304  }
1305  else
1306  {
1307  if (Verbosity() > 0)
1308  {
1309  cout << __PRETTY_FUNCTION__ << " could not find " << nodename << endl;
1310  }
1311  continue;
1312  }
1313  if(_do_HITS_ABSORBER){
1314  string nodenameAbs = "G4HIT_ABSORBER_" + GetProjectionNameFromIndex(iIndex);
1315  PHG4HitContainer* hitsAbs = findNode::getClass<PHG4HitContainer>(topNode, nodenameAbs);
1316  if (hitsAbs)
1317  {
1318  if (Verbosity() > 1)
1319  {
1320  cout << __PRETTY_FUNCTION__ << " absorber number of hits: " << hitsAbs->size() << endl;
1321  }
1322  PHG4HitContainer::ConstRange hit_range = hitsAbs->getHits();
1323  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
1324  {
1325  if (Verbosity() > 1)
1326  {
1327  cout << __PRETTY_FUNCTION__ << " found hit with id " << hit_iter->second->get_trkid() << endl;
1328  }
1329  if (_nHitsLayers > _maxNHits)
1330  {
1331  cout << __PRETTY_FUNCTION__ << " exceededed maximum hit array size! Please check where these hits come from!" << endl;
1332  break;
1333  }
1334 
1335  _hits_x[_nHitsLayers] = hit_iter->second->get_x(0);
1336  _hits_y[_nHitsLayers] = hit_iter->second->get_y(0);
1337  _hits_z[_nHitsLayers] = hit_iter->second->get_z(0);
1338  _hits_x2[_nHitsLayers] = hit_iter->second->get_x(1);
1339  _hits_y2[_nHitsLayers] = hit_iter->second->get_y(1);
1340  _hits_z2[_nHitsLayers] = hit_iter->second->get_z(1);
1341  _hits_t[_nHitsLayers] = hit_iter->second->get_t(0);
1342  _hits_edep[_nHitsLayers] = hit_iter->second->get_edep();
1343  _hits_lightyield[_nHitsLayers] = hit_iter->second->get_light_yield();
1344  _hits_isAbsorber[_nHitsLayers] = 1;
1345  _hits_layerID[_nHitsLayers] = iIndex;
1346  // cout << "i " << hit_iter->second->get_index_i() << "\tj " <<hit_iter->second->get_index_j() << "\tk " <<hit_iter->second->get_index_k() << "\tl " << hit_iter->second->get_index_l() << "\tsens_x "<< hit_iter->second->get_strip_z_index()<< "\tsens_y "<< hit_iter->second->get_strip_y_index() << endl;
1347  if (truthinfocontainerHits)
1348  {
1349  PHG4Particle* particle = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
1350 
1351  if (particle->get_parent_id() != 0)
1352  {
1353  PHG4Particle* g4particleMother = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
1354  int mcSteps = 0;
1355  while (g4particleMother->get_parent_id() != 0)
1356  {
1357  g4particleMother = truthinfocontainerHits->GetParticle(g4particleMother->get_parent_id());
1358  if (g4particleMother == NULL) break;
1359  mcSteps += 1;
1360  }
1361  if (mcSteps <= _depth_MCstack)
1362  {
1363  _hits_trueID[_nHitsLayers] = hit_iter->second->get_trkid();
1364  }
1365  else
1366  {
1367  PHG4Particle* g4particleMother2 = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
1368  int mcSteps2 = 0;
1369  while (g4particleMother2->get_parent_id() != 0 && (mcSteps2 < (mcSteps - _depth_MCstack + 1)))
1370  {
1371  g4particleMother2 = truthinfocontainerHits->GetParticle(g4particleMother2->get_parent_id());
1372  if (g4particleMother2 == NULL)
1373  {
1374  break;
1375  }
1376  else
1377  {
1378  _hits_trueID[_nHitsLayers] = g4particleMother2->get_parent_id();
1379  mcSteps2 += 1;
1380  }
1381  }
1382  }
1383  }
1384  else
1385  {
1386  _hits_trueID[_nHitsLayers] = hit_iter->second->get_trkid();
1387  }
1388  }
1389  _nHitsLayers++;
1390  }
1391  if (Verbosity() > 0)
1392  {
1393  cout << "saved\t" << _nHitsLayers << "\thits for absorber of " << GetProjectionNameFromIndex(iIndex) << endl;
1394  }
1395  }
1396  else
1397  {
1398  if (Verbosity() > 0)
1399  {
1400  cout << __PRETTY_FUNCTION__ << " could not find " << nodename << endl;
1401  }
1402  continue;
1403  }
1404  }
1405  }
1406  }
1407  }
1408  //----------------------
1409  // TOWERS FHCAL
1410  //----------------------
1411  if (_do_FHCAL)
1412  {
1413  CaloRawTowerEval* towerevalFHCAL = _caloevalstackFHCAL->get_rawtower_eval();
1414  _nTowers_FHCAL = 0;
1415  string towernodeFHCAL = "TOWER_CALIB_FHCAL";
1416  RawTowerContainer* towersFHCAL = findNode::getClass<RawTowerContainer>(topNode, towernodeFHCAL.c_str());
1417  if (towersFHCAL)
1418  {
1419  if (Verbosity() > 0)
1420  {
1421  cout << "saving HCAL towers" << endl;
1422  }
1423  string towergeomnodeFHCAL = "TOWERGEOM_FHCAL";
1424  RawTowerGeomContainer* towergeomFHCAL = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeFHCAL.c_str());
1425  if (towergeomFHCAL)
1426  {
1427  if (_do_GEOMETRY && !_geometry_done[kFHCAL])
1428  {
1429  RawTowerGeomContainer::ConstRange all_towers = towergeomFHCAL->get_tower_geometries();
1430  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
1431  it != all_towers.second; ++it)
1432  {
1433  _calo_ID = kFHCAL;
1434  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
1435  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
1436  _calo_towers_iL[_calo_towers_N] = -1;
1437  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
1438  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
1439  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
1440  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
1441  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
1442  _calo_towers_N++;
1443  }
1444  _geometry_done[kFHCAL] = 1;
1445  _geometry_tree->Fill();
1446  resetGeometryArrays();
1447  }
1448 
1449  RawTowerContainer::ConstRange begin_end = towersFHCAL->getTowers();
1450  RawTowerContainer::ConstIterator rtiter;
1451  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
1452  {
1453  RawTower* tower = rtiter->second;
1454  if (tower)
1455  {
1456  // min energy cut
1457  if (tower->get_energy() < _reco_e_threshold[kFHCAL]) continue;
1458  // cout << "\tnew FHCAL tower" << endl;
1459  _tower_FHCAL_iEta[_nTowers_FHCAL] = tower->get_bineta();
1460  _tower_FHCAL_iPhi[_nTowers_FHCAL] = tower->get_binphi();
1461  _tower_FHCAL_E[_nTowers_FHCAL] = tower->get_energy();
1462 
1463  // PHG4TruthInfoContainer* truthinfocontainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1464  // RawTower::ShowerConstRange shower_range = tower->get_g4showers();
1465  // for (RawTower::ShowerConstIterator iter = shower_range.first;
1466  // iter != shower_range.second;
1467  // ++iter)
1468  // {
1469  // PHG4Shower* shower = truthinfocontainer->GetShower(iter->first);
1470  // // PHG4Particle* particleParent = truthinfocontainer->GetParticle(shower->get_parent_particle_id());
1471  // // if (particleParent)
1472  // // {
1473  // // cout << "\t\tcurr shower parent id: " << shower->get_parent_particle_id()<< "\tPDG " << particleParent->get_pid() << "\tedep: " << shower->get_edep() << endl;
1474  // // } else {
1475  // // cout << "\t\tcurr shower parent id: " << shower->get_parent_particle_id() << endl;
1476  // // }
1477  // // for (PHG4Shower::ParticleIdIter jter = shower->begin_g4particle_id();jter != shower->end_g4particle_id(); ++jter)
1478  // // {
1479  // // int g4particle_id = *jter;
1480  // // PHG4Particle* particle = truthinfocontainer->GetParticle(g4particle_id);
1481  // // if (particle)
1482  // // {
1483  // // if(particle->get_e()>0.01*shower->get_edep()){
1484  // // // cout << "\t\t\tparticle ID in shower: " << particle->get_track_id() << "\tPDG " << particle->get_pid() << "\tenergy: " << particle->get_e()<< endl;
1485  // // }
1486  // // }
1487  // // }
1488  // }
1489 
1490  PHG4Particle* primary = towerevalFHCAL->max_truth_primary_particle_by_energy(tower);
1491  if (primary)
1492  {
1493  _tower_FHCAL_trueID[_nTowers_FHCAL] = primary->get_track_id();
1494  // gflavor = primary->get_pid();
1495  // efromtruth = towerevalFHCAL->get_energy_contribution(tower, primary);
1496  }
1497  else
1498  {
1499  _tower_FHCAL_trueID[_nTowers_FHCAL] = -10;
1500  }
1501  _nTowers_FHCAL++;
1502  }
1503  }
1504  }
1505  else
1506  {
1507  if (Verbosity() > 0)
1508  {
1509  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeFHCAL << endl;
1510  }
1511  // return;
1512  }
1513  if (Verbosity() > 0)
1514  {
1515  cout << "saved\t" << _nTowers_FHCAL << "\tFHCAL towers" << endl;
1516  }
1517  }
1518  else
1519  {
1520  if (Verbosity() > 0)
1521  {
1522  cout << PHWHERE << " ERROR: Can't find " << towernodeFHCAL << endl;
1523  }
1524  // return;
1525  }
1526  }
1527  //----------------------
1528  // TOWERS BECAL
1529  //----------------------
1530  if (_do_BECAL)
1531  {
1532  CaloRawTowerEval* towerevalBECAL = _caloevalstackBECAL->get_rawtower_eval();
1533  _nTowers_BECAL = 0;
1534  string towernodeBECAL = "TOWER_CALIB_BECAL";
1535  RawTowerContainer* towersBECAL = findNode::getClass<RawTowerContainer>(topNode, towernodeBECAL.c_str());
1536 
1537  if (Verbosity() > 0)
1538  {
1539  RawTowerContainer* towersBECAL1 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_RAW_BECAL");
1540  RawTowerContainer* towersBECAL2 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_SIM_BECAL");
1541  cout << "BECAL sim: " << towersBECAL2->size() << endl;
1542  cout << "BECAL raw: " << towersBECAL1->size() << endl;
1543  cout << "BECAL calib: " << towersBECAL->size() << endl;
1544  }
1545  if (towersBECAL)
1546  {
1547  if (Verbosity() > 0)
1548  {
1549  cout << "saving BECAL towers" << endl;
1550  }
1551  string towergeomnodeBECAL = "TOWERGEOM_BECAL";
1552  RawTowerGeomContainer* towergeomBECAL = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeBECAL.c_str());
1553  if (towergeomBECAL)
1554  {
1555  if (_do_GEOMETRY && !_geometry_done[kBECAL])
1556  {
1557  // std::ostream *fout= new ofstream("test.list");
1558  RawTowerGeomContainer::ConstRange all_towers = towergeomBECAL->get_tower_geometries();
1559  // *fout << "Calorimeter ID: " << towergeomBECAL->get_calorimeter_id() << endl;
1560  // *fout << "size: " << towergeomBECAL->size() << endl;
1561  // towergeomBECAL->identify(*fout);
1562  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
1563  it != all_towers.second; ++it)
1564  {
1565  _calo_ID = kBECAL;
1566  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
1567  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
1568  _calo_towers_iL[_calo_towers_N] = -1;
1569  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
1570  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
1571  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
1572  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
1573  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
1574  // it->second->identify(*fout);
1575  _calo_towers_N++;
1576  }
1577  _geometry_done[kBECAL] = 1;
1578  _geometry_tree->Fill();
1579  resetGeometryArrays();
1580  // _calo_ID = kBECAL;
1581  // _calo_ID = kBECAL;
1582  }
1583 
1584  RawTowerContainer::ConstRange begin_end = towersBECAL->getTowers();
1585  RawTowerContainer::ConstIterator rtiter;
1586  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
1587  {
1588  RawTower* tower = rtiter->second;
1589  if (tower)
1590  {
1591  // min energy cut
1592  if (tower->get_energy() < _reco_e_threshold[kBECAL]) continue;
1593  _tower_BECAL_iEta[_nTowers_BECAL] = tower->get_bineta();
1594  _tower_BECAL_iPhi[_nTowers_BECAL] = tower->get_binphi();
1595  _tower_BECAL_E[_nTowers_BECAL] = tower->get_energy();
1596 
1597  PHG4Particle* primary = towerevalBECAL->max_truth_primary_particle_by_energy(tower);
1598  if (primary)
1599  {
1600  _tower_BECAL_trueID[_nTowers_BECAL] = primary->get_track_id();
1601  // gflavor = primary->get_pid();
1602  // efromtruth = towerevalBECAL->get_energy_contribution(tower, primary);
1603  }
1604  else
1605  {
1606  _tower_BECAL_trueID[_nTowers_BECAL] = -10;
1607  }
1608  _nTowers_BECAL++;
1609  }
1610  }
1611  }
1612  else
1613  {
1614  if (Verbosity() > 0)
1615  {
1616  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeBECAL << endl;
1617  }
1618  // return;
1619  }
1620  if (Verbosity() > 0)
1621  {
1622  cout << "saved\t" << _nTowers_BECAL << "\tBECAL towers" << endl;
1623  }
1624  }
1625  else
1626  {
1627  if (Verbosity() > 0)
1628  {
1629  cout << PHWHERE << " ERROR: Can't find " << towernodeBECAL << endl;
1630  }
1631  // return;
1632  }
1633  }
1634  //----------------------
1635  // TOWERS HCALIN
1636  //----------------------
1637  if (_do_HCALIN)
1638  {
1639  CaloRawTowerEval* towerevalHCALIN = _caloevalstackHCALIN->get_rawtower_eval();
1640  _nTowers_HCALIN = 0;
1641  string towernodeHCALIN = "TOWER_CALIB_HCALIN";
1642  RawTowerContainer* towersHCALIN = findNode::getClass<RawTowerContainer>(topNode, towernodeHCALIN.c_str());
1643  if (towersHCALIN)
1644  {
1645  if (Verbosity() > 0)
1646  {
1647  cout << "saving HCAL towers" << endl;
1648  }
1649  string towergeomnodeHCALIN = "TOWERGEOM_HCALIN";
1650  RawTowerGeomContainer* towergeomHCALIN = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeHCALIN.c_str());
1651  if (towergeomHCALIN)
1652  {
1653  if (_do_GEOMETRY && !_geometry_done[kHCALIN])
1654  {
1655  RawTowerGeomContainer::ConstRange all_towers = towergeomHCALIN->get_tower_geometries();
1656  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
1657  it != all_towers.second; ++it)
1658  {
1659  _calo_ID = kHCALIN;
1660  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
1661  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
1662  _calo_towers_iL[_calo_towers_N] = -1;
1663  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
1664  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
1665  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
1666  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
1667  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
1668  _calo_towers_N++;
1669  }
1670  _geometry_done[kHCALIN] = 1;
1671  _geometry_tree->Fill();
1672  resetGeometryArrays();
1673  }
1674  RawTowerContainer::ConstRange begin_end = towersHCALIN->getTowers();
1675  RawTowerContainer::ConstIterator rtiter;
1676  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
1677  {
1678  RawTower* tower = rtiter->second;
1679  if (tower)
1680  {
1681  // min energy cut
1682  if (tower->get_energy() < _reco_e_threshold[kHCALIN]) continue;
1683  _tower_HCALIN_iEta[_nTowers_HCALIN] = tower->get_bineta();
1684  _tower_HCALIN_iPhi[_nTowers_HCALIN] = tower->get_binphi();
1685  _tower_HCALIN_E[_nTowers_HCALIN] = tower->get_energy();
1686 
1687  PHG4Particle* primary = towerevalHCALIN->max_truth_primary_particle_by_energy(tower);
1688  if (primary)
1689  {
1690  _tower_HCALIN_trueID[_nTowers_HCALIN] = primary->get_track_id();
1691  // gflavor = primary->get_pid();
1692  // efromtruth = towerevalHCALIN->get_energy_contribution(tower, primary);
1693  }
1694  else
1695  {
1696  _tower_HCALIN_trueID[_nTowers_HCALIN] = -10;
1697  }
1698  _nTowers_HCALIN++;
1699  }
1700  }
1701  }
1702  else
1703  {
1704  if (Verbosity() > 0)
1705  {
1706  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeHCALIN << endl;
1707  }
1708  // return;
1709  }
1710  if (Verbosity() > 0)
1711  {
1712  cout << "saved\t" << _nTowers_HCALIN << "\tHCALIN towers" << endl;
1713  }
1714  }
1715  else
1716  {
1717  if (Verbosity() > 0)
1718  {
1719  cout << PHWHERE << " ERROR: Can't find " << towernodeHCALIN << endl;
1720  }
1721  // return;
1722  }
1723  }
1724  //----------------------
1725  // TOWERS HCALOUT
1726  //----------------------
1727  if (_do_HCALOUT)
1728  {
1729  CaloRawTowerEval* towerevalHCALOUT = _caloevalstackHCALOUT->get_rawtower_eval();
1730  _nTowers_HCALOUT = 0;
1731  string towernodeHCALOUT = "TOWER_CALIB_HCALOUT";
1732  RawTowerContainer* towersHCALOUT = findNode::getClass<RawTowerContainer>(topNode, towernodeHCALOUT.c_str());
1733  if (towersHCALOUT)
1734  {
1735  if (Verbosity() > 0)
1736  {
1737  cout << "saving HCAL towers" << endl;
1738  }
1739  string towergeomnodeHCALOUT = "TOWERGEOM_HCALOUT";
1740  RawTowerGeomContainer* towergeomHCALOUT = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeHCALOUT.c_str());
1741  if (towergeomHCALOUT)
1742  {
1743  if (_do_GEOMETRY && !_geometry_done[kHCALOUT])
1744  {
1745  RawTowerGeomContainer::ConstRange all_towers = towergeomHCALOUT->get_tower_geometries();
1746  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
1747  it != all_towers.second; ++it)
1748  {
1749  _calo_ID = kHCALOUT;
1750  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
1751  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
1752  _calo_towers_iL[_calo_towers_N] = -1;
1753  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
1754  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
1755  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
1756  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
1757  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
1758  _calo_towers_N++;
1759  }
1760  _geometry_done[kHCALOUT] = 1;
1761  _geometry_tree->Fill();
1762  resetGeometryArrays();
1763  }
1764  RawTowerContainer::ConstRange begin_end = towersHCALOUT->getTowers();
1765  RawTowerContainer::ConstIterator rtiter;
1766  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
1767  {
1768  RawTower* tower = rtiter->second;
1769  if (tower)
1770  {
1771  // min energy cut
1772  if (tower->get_energy() < _reco_e_threshold[kHCALOUT]) continue;
1773  _tower_HCALOUT_iEta[_nTowers_HCALOUT] = tower->get_bineta();
1774  _tower_HCALOUT_iPhi[_nTowers_HCALOUT] = tower->get_binphi();
1775  _tower_HCALOUT_E[_nTowers_HCALOUT] = tower->get_energy();
1776 
1777  PHG4Particle* primary = towerevalHCALOUT->max_truth_primary_particle_by_energy(tower);
1778  if (primary)
1779  {
1780  _tower_HCALOUT_trueID[_nTowers_HCALOUT] = primary->get_track_id();
1781  // gflavor = primary->get_pid();
1782  // efromtruth = towerevalHCALOUT->get_energy_contribution(tower, primary);
1783  }
1784  else
1785  {
1786  _tower_HCALOUT_trueID[_nTowers_HCALOUT] = -10;
1787  }
1788  _nTowers_HCALOUT++;
1789  }
1790  }
1791  }
1792  else
1793  {
1794  if (Verbosity() > 0)
1795  {
1796  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeHCALOUT << endl;
1797  }
1798  // return;
1799  }
1800  if (Verbosity() > 0)
1801  {
1802  cout << "saved\t" << _nTowers_HCALOUT << "\tHCALOUT towers" << endl;
1803  }
1804  }
1805  else
1806  {
1807  if (Verbosity() > 0)
1808  {
1809  cout << PHWHERE << " ERROR: Can't find " << towernodeHCALOUT << endl;
1810  }
1811  // return;
1812  }
1813  }
1814  //----------------------
1815  // TOWERS EHCAL
1816  //----------------------
1817  if (_do_EHCAL)
1818  {
1819  CaloRawTowerEval* towerevalEHCAL = _caloevalstackEHCAL->get_rawtower_eval();
1820  _nTowers_EHCAL = 0;
1821  string towernodeEHCAL = "TOWER_CALIB_EHCAL";
1822  RawTowerContainer* towersEHCAL = findNode::getClass<RawTowerContainer>(topNode, towernodeEHCAL.c_str());
1823  if (towersEHCAL)
1824  {
1825  if (Verbosity() > 0)
1826  {
1827  cout << "saving HCAL towers" << endl;
1828  }
1829  string towergeomnodeEHCAL = "TOWERGEOM_EHCAL";
1830  RawTowerGeomContainer* towergeomEHCAL = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeEHCAL.c_str());
1831  if (towergeomEHCAL)
1832  {
1833  if (_do_GEOMETRY && !_geometry_done[kEHCAL])
1834  {
1835  RawTowerGeomContainer::ConstRange all_towers = towergeomEHCAL->get_tower_geometries();
1836  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
1837  it != all_towers.second; ++it)
1838  {
1839  _calo_ID = kEHCAL;
1840  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
1841  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
1842  _calo_towers_iL[_calo_towers_N] = -1;
1843  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
1844  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
1845  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
1846  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
1847  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
1848  _calo_towers_N++;
1849  }
1850  _geometry_done[kEHCAL] = 1;
1851  _geometry_tree->Fill();
1852  resetGeometryArrays();
1853  }
1854  RawTowerContainer::ConstRange begin_end = towersEHCAL->getTowers();
1855  RawTowerContainer::ConstIterator rtiter;
1856  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
1857  {
1858  RawTower* tower = rtiter->second;
1859  if (tower)
1860  {
1861  // min energy cut
1862  if (tower->get_energy() < _reco_e_threshold[kEHCAL]) continue;
1863  _tower_EHCAL_iEta[_nTowers_EHCAL] = tower->get_bineta();
1864  _tower_EHCAL_iPhi[_nTowers_EHCAL] = tower->get_binphi();
1865  _tower_EHCAL_E[_nTowers_EHCAL] = tower->get_energy();
1866 
1867  PHG4Particle* primary = towerevalEHCAL->max_truth_primary_particle_by_energy(tower);
1868  if (primary)
1869  {
1870  _tower_EHCAL_trueID[_nTowers_EHCAL] = primary->get_track_id();
1871  // gflavor = primary->get_pid();
1872  // efromtruth = towerevalEHCAL->get_energy_contribution(tower, primary);
1873  }
1874  else
1875  {
1876  _tower_EHCAL_trueID[_nTowers_EHCAL] = -10;
1877  }
1878  _nTowers_EHCAL++;
1879  }
1880  }
1881  }
1882  else
1883  {
1884  if (Verbosity() > 0)
1885  {
1886  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeEHCAL << endl;
1887  }
1888  // return;
1889  }
1890  if (Verbosity() > 0)
1891  {
1892  cout << "saved\t" << _nTowers_EHCAL << "\tEHCAL towers" << endl;
1893  }
1894  }
1895  else
1896  {
1897  if (Verbosity() > 0)
1898  {
1899  cout << PHWHERE << " ERROR: Can't find " << towernodeEHCAL << endl;
1900  }
1901  // return;
1902  }
1903  }
1904  //----------------------
1905  // TOWERS DRCALO
1906  //----------------------
1907  if (_do_DRCALO)
1908  {
1909  CaloRawTowerEval* towerevalDRCALO = _caloevalstackDRCALO->get_rawtower_eval();
1910  _nTowers_DRCALO = 0;
1911  string towernodeDRCALO = "TOWER_CALIB_DRCALO";
1912  RawTowerContainer* towersDRCALO = findNode::getClass<RawTowerContainer>(topNode, towernodeDRCALO.c_str());
1913  if (towersDRCALO)
1914  {
1915  if (Verbosity() > 0)
1916  {
1917  cout << "saving DRCALO towers" << endl;
1918  }
1919  string towergeomnodeDRCALO = "TOWERGEOM_DRCALO";
1920  RawTowerGeomContainer* towergeomDRCALO = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeDRCALO.c_str());
1921  if (towergeomDRCALO)
1922  {
1923  if (_do_GEOMETRY && !_geometry_done[kDRCALO])
1924  {
1925  RawTowerGeomContainer::ConstRange all_towers = towergeomDRCALO->get_tower_geometries();
1926  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
1927  it != all_towers.second; ++it)
1928  {
1929  _calo_ID = kDRCALO;
1930  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
1931  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
1932  _calo_towers_iL[_calo_towers_N] = -1;
1933  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
1934  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
1935  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
1936  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
1937  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
1938  _calo_towers_N++;
1939  }
1940  _geometry_done[kDRCALO] = 1;
1941  _geometry_tree->Fill();
1942  resetGeometryArrays();
1943  }
1944  if (Verbosity() > 0)
1945  {
1946  cout << "found DRCALO geom" << endl;
1947  }
1948  RawTowerContainer::ConstRange begin_end = towersDRCALO->getTowers();
1949  RawTowerContainer::ConstIterator rtiter;
1950  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
1951  {
1952  RawTower* tower = rtiter->second;
1953  // if (dynamic_cast<RawTowerv2 *>(tower) == nullptr)
1954  // {
1955  // cout << __PRETTY_FUNCTION__ << " : Fatal Error! "
1956  // << "Expect RawTowerv2, but found this tower:";
1957  // tower->identify();
1958  // // exit(1);
1959  // }
1960  if (tower)
1961  {
1962  // min energy cut
1963  if (tower->get_energy() < _reco_e_threshold[kDRCALO]) continue;
1964 
1965  _tower_DRCALO_iEta[_nTowers_DRCALO] = tower->get_bineta();
1966  _tower_DRCALO_iPhi[_nTowers_DRCALO] = tower->get_binphi();
1967  _tower_DRCALO_E[_nTowers_DRCALO] = tower->get_energy();
1968  // cout << __LINE__ << endl;
1969  _tower_DRCALO_NScint[_nTowers_DRCALO] = tower->get_scint_gammas();
1970  _tower_DRCALO_NCerenkov[_nTowers_DRCALO] = tower->get_cerenkov_gammas();
1971  // cout << "sci gammas: " << tower->get_scint_gammas() << "\tcerenk gammas: " << tower->get_cerenkov_gammas() << endl;
1972 
1973  PHG4Particle* primary = towerevalDRCALO->max_truth_primary_particle_by_energy(tower);
1974  if (primary)
1975  {
1976  _tower_DRCALO_trueID[_nTowers_DRCALO] = primary->get_track_id();
1977  // gflavor = primary->get_pid();
1978  // efromtruth = towerevalDRCALO->get_energy_contribution(tower, primary);
1979  }
1980  else
1981  {
1982  _tower_DRCALO_trueID[_nTowers_DRCALO] = -10;
1983  }
1984  _nTowers_DRCALO++;
1985  }
1986  }
1987  if (Verbosity() > 0)
1988  {
1989  cout << "finished DRCALO twr loop" << endl;
1990  }
1991  }
1992  else
1993  {
1994  if (Verbosity() > 0)
1995  {
1996  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeDRCALO << endl;
1997  }
1998  // return;
1999  }
2000  if (Verbosity() > 0)
2001  {
2002  cout << "saved\t" << _nTowers_DRCALO << "\tDRCALO towers" << endl;
2003  }
2004  }
2005  else
2006  {
2007  if (Verbosity() > 0)
2008  {
2009  cout << PHWHERE << " ERROR: Can't find " << towernodeDRCALO << endl;
2010  }
2011  // return;
2012  }
2013  }
2014  //----------------------
2015  // TOWERS FOCAL
2016  //----------------------
2017  if (_do_FOCAL)
2018  {
2019  CaloRawTowerEval* towerevalFOCAL = _caloevalstackFOCAL->get_rawtower_eval();
2020  _nTowers_FOCAL = 0;
2021  string towernodeFOCAL = "TOWER_CALIB_FOCAL";
2022  RawTowerContainer* towersFOCAL = findNode::getClass<RawTowerContainer>(topNode, towernodeFOCAL.c_str());
2023  if (towersFOCAL)
2024  {
2025  if (Verbosity() > 0)
2026  {
2027  cout << "saving FOCAL towers" << endl;
2028  }
2029  string towergeomnodeFOCAL = "TOWERGEOM_FOCAL";
2030  RawTowerGeomContainer* towergeomFOCAL = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeFOCAL.c_str());
2031  if (towergeomFOCAL)
2032  {
2033  if (_do_GEOMETRY && !_geometry_done[kFOCAL])
2034  {
2035  RawTowerGeomContainer::ConstRange all_towers = towergeomFOCAL->get_tower_geometries();
2036  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
2037  it != all_towers.second; ++it)
2038  {
2039  _calo_ID = kFOCAL;
2040  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
2041  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
2042  _calo_towers_iL[_calo_towers_N] = -1;
2043  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
2044  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
2045  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
2046  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
2047  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
2048  _calo_towers_N++;
2049  }
2050  _geometry_done[kFOCAL] = 1;
2051  _geometry_tree->Fill();
2052  resetGeometryArrays();
2053  }
2054  if (Verbosity() > 0)
2055  {
2056  cout << "found FOCAL geom" << endl;
2057  }
2058  RawTowerContainer::ConstRange begin_end = towersFOCAL->getTowers();
2059  RawTowerContainer::ConstIterator rtiter;
2060  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
2061  {
2062  RawTower* tower = rtiter->second;
2063  if (tower)
2064  {
2065  // min energy cut
2066  if (tower->get_energy() < _reco_e_threshold[kFOCAL]) continue;
2067 
2068  _tower_FOCAL_iEta[_nTowers_FOCAL] = tower->get_bineta();
2069  _tower_FOCAL_iPhi[_nTowers_FOCAL] = tower->get_binphi();
2070  _tower_FOCAL_E[_nTowers_FOCAL] = tower->get_energy();
2071  _tower_FOCAL_NScint[_nTowers_FOCAL] = tower->get_scint_gammas();
2072  _tower_FOCAL_NCerenkov[_nTowers_FOCAL] = tower->get_cerenkov_gammas();
2073 
2074  PHG4Particle* primary = towerevalFOCAL->max_truth_primary_particle_by_energy(tower);
2075  if (primary)
2076  {
2077  _tower_FOCAL_trueID[_nTowers_FOCAL] = primary->get_track_id();
2078  }
2079  else
2080  {
2081  _tower_FOCAL_trueID[_nTowers_FOCAL] = -10;
2082  }
2083  _nTowers_FOCAL++;
2084  }
2085  }
2086  if (Verbosity() > 0)
2087  {
2088  cout << "finished FOCAL twr loop" << endl;
2089  }
2090  }
2091  else
2092  {
2093  if (Verbosity() > 0)
2094  {
2095  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeFOCAL << endl;
2096  }
2097  // return;
2098  }
2099  if (Verbosity() > 0)
2100  {
2101  cout << "saved\t" << _nTowers_FOCAL << "\tFOCAL towers" << endl;
2102  }
2103  }
2104  else
2105  {
2106  if (Verbosity() > 0)
2107  {
2108  cout << PHWHERE << " ERROR: Can't find " << towernodeFOCAL << endl;
2109  }
2110  // return;
2111  }
2112  }
2113  //----------------------
2114  // TOWERS LFHCAL
2115  //----------------------
2116  if (_do_LFHCAL)
2117  {
2118  CaloRawTowerEval* towerevalLFHCAL = _caloevalstackLFHCAL->get_rawtower_eval();
2119  _nTowers_LFHCAL = 0;
2120  string towernodeLFHCAL = "TOWER_CALIB_LFHCAL";
2121  RawTowerContainer* towersLFHCAL = findNode::getClass<RawTowerContainer>(topNode, towernodeLFHCAL.c_str());
2122  if (Verbosity() > 1)
2123  std::cout << "reading towers: " << towersLFHCAL->size() << std::endl;
2124  if (towersLFHCAL)
2125  {
2126  if (Verbosity() > 0)
2127  {
2128  cout << "saving LFHCAL towers" << endl;
2129  }
2130  string towergeomnodeLFHCAL = "TOWERGEOM_LFHCAL";
2131  RawTowerGeomContainer* towergeomLFHCAL = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeLFHCAL.c_str());
2132  if (towergeomLFHCAL)
2133  {
2134  if (_do_GEOMETRY && !_geometry_done[kLFHCAL])
2135  {
2136  RawTowerGeomContainer::ConstRange all_towers = towergeomLFHCAL->get_tower_geometries();
2137  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
2138  it != all_towers.second; ++it)
2139  {
2140  _calo_ID = kLFHCAL;
2141  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
2142  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
2143  _calo_towers_iL[_calo_towers_N] = it->second->get_binl();
2144  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
2145  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
2146  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
2147  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
2148  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
2149  _calo_towers_N++;
2150  }
2151  _geometry_done[kLFHCAL] = 1;
2152  _geometry_tree->Fill();
2153  resetGeometryArrays();
2154  }
2155  if (Verbosity() > 0)
2156  {
2157  cout << "found LFHCAL geom" << endl;
2158  }
2159 
2160  RawTowerContainer::ConstRange begin_end = towersLFHCAL->getTowers();
2161  RawTowerContainer::ConstIterator rtiter;
2162 
2163  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
2164  {
2165  RawTower* tower = rtiter->second;
2166  if (tower)
2167  {
2168  // min energy cut
2169  if (tower->get_energy() < _reco_e_threshold[kLFHCAL]) continue;
2170  if (Verbosity() > 1) cout << "\n event eval: \t" << tower->get_energy() << "\t ieta: " << tower->get_bineta() << "\t iphi: " << tower->get_binphi() << "\t iZ: " << tower->get_binl() << endl;
2171  _tower_LFHCAL_iEta[_nTowers_LFHCAL] = tower->get_bineta();
2172  _tower_LFHCAL_iPhi[_nTowers_LFHCAL] = tower->get_binphi();
2173  _tower_LFHCAL_iL[_nTowers_LFHCAL] = tower->get_binl();
2174  _tower_LFHCAL_E[_nTowers_LFHCAL] = tower->get_energy();
2175  PHG4Particle* primary = towerevalLFHCAL->max_truth_primary_particle_by_energy(tower);
2176  if (primary)
2177  {
2178  _tower_LFHCAL_trueID[_nTowers_LFHCAL] = primary->get_track_id();
2179  // gflavor = primary->get_pid();
2180  // efromtruth = towerevalLFHCAL->get_energy_contribution(tower, primary);
2181  }
2182  else
2183  {
2184  _tower_LFHCAL_trueID[_nTowers_LFHCAL] = -10;
2185  }
2186  _nTowers_LFHCAL++;
2187  }
2188  }
2189  }
2190  else
2191  {
2192  if (Verbosity() > 0)
2193  {
2194  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeLFHCAL << endl;
2195  }
2196  // return;
2197  }
2198  if (Verbosity() > 0)
2199  {
2200  cout << "saved\t" << _nTowers_LFHCAL << "\tLFHCAL towers" << endl;
2201  }
2202  }
2203  else
2204  {
2205  if (Verbosity() > 0)
2206  {
2207  cout << PHWHERE << " ERROR: Can't find " << towernodeLFHCAL << endl;
2208  }
2209  // return;
2210  }
2211  }
2212 
2213  //----------------------
2214  // TOWERS FEMC
2215  //----------------------
2216  if (_do_FEMC)
2217  {
2218  CaloRawTowerEval* towerevalFEMC = _caloevalstackFEMC->get_rawtower_eval();
2219  _nTowers_FEMC = 0;
2220  string towernodeFEMC = "TOWER_CALIB_FEMC";
2221  RawTowerContainer* towersFEMC = findNode::getClass<RawTowerContainer>(topNode, towernodeFEMC.c_str());
2222  if (towersFEMC)
2223  {
2224  if (Verbosity() > 0)
2225  {
2226  cout << "saving EMC towers" << endl;
2227  }
2228  string towergeomnodeFEMC = "TOWERGEOM_FEMC";
2229  RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeFEMC.c_str());
2230  if (towergeom)
2231  {
2232  if (_do_GEOMETRY && !_geometry_done[kFEMC])
2233  {
2234  RawTowerGeomContainer::ConstRange all_towers = towergeom->get_tower_geometries();
2235  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
2236  it != all_towers.second; ++it)
2237  {
2238  _calo_ID = kFEMC;
2239  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
2240  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
2241  _calo_towers_iL[_calo_towers_N] = -1;
2242  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
2243  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
2244  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
2245  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
2246  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
2247  _calo_towers_N++;
2248  }
2249  _geometry_done[kFEMC] = 1;
2250  _geometry_tree->Fill();
2251  resetGeometryArrays();
2252  }
2253  RawTowerContainer::ConstRange begin_end = towersFEMC->getTowers();
2254  RawTowerContainer::ConstIterator rtiter;
2255  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
2256  {
2257  RawTower* tower = rtiter->second;
2258  if (tower)
2259  {
2260  // min energy cut
2261  if (tower->get_energy() < _reco_e_threshold[kFEMC]) continue;
2262 
2263  _tower_FEMC_iEta[_nTowers_FEMC] = tower->get_bineta();
2264  _tower_FEMC_iPhi[_nTowers_FEMC] = tower->get_binphi();
2265  _tower_FEMC_E[_nTowers_FEMC] = tower->get_energy();
2266 
2267  // cout << "\tnew FEMC tower " << tower->get_energy() << endl;
2268 
2269  PHG4Particle* primary = towerevalFEMC->max_truth_primary_particle_by_energy(tower);
2270  if (primary)
2271  {
2272  _tower_FEMC_trueID[_nTowers_FEMC] = primary->get_track_id();
2273  // gflavor = primary->get_pid();
2274  // efromtruth = towerevalFEMC->get_energy_contribution(tower, primary);
2275  }
2276  else
2277  {
2278  _tower_FEMC_trueID[_nTowers_FEMC] = -10;
2279  }
2280 
2281  // PHG4TruthInfoContainer* truthinfocontainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
2282  // RawTower::ShowerConstRange shower_range = tower->get_g4showers();
2283  // for (RawTower::ShowerConstIterator iter = shower_range.first;
2284  // iter != shower_range.second;
2285  // ++iter)
2286  // {
2287  // PHG4Shower* shower = truthinfocontainer->GetShower(iter->first);
2288  // PHG4Particle* particleParent = truthinfocontainer->GetParticle(shower->get_parent_particle_id());
2289  // if (particleParent)
2290  // {
2291  // cout << "\t\tcurr shower parent id: " << shower->get_parent_particle_id() << " ("<< _tower_FEMC_trueID[_nTowers_FEMC] << ")\tPDG " << particleParent->get_pid() << "\tedep: " << shower->get_edep() << endl;
2292  // } else {
2293  // cout << "\t\tcurr shower parent id: " << shower->get_parent_particle_id() << endl;
2294  // }
2295  // for (PHG4Shower::ParticleIdIter jter = shower->begin_g4particle_id();jter != shower->end_g4particle_id(); ++jter)
2296  // {
2297  // int g4particle_id = *jter;
2298  // PHG4Particle* particle = truthinfocontainer->GetParticle(g4particle_id);
2299  // if (particle)
2300  // {
2301  // // if(particle->get_e()>0.1*shower->get_edep()){
2302  // // cout << "\t\t\tparticle ID in shower: " << particle->get_track_id() << "\tPDG " << particle->get_pid() << "\tenergy: " << particle->get_e()<< endl;
2303  // // }
2304  // }
2305  // }
2306  // }
2307  _nTowers_FEMC++;
2308  }
2309  }
2310  }
2311  else
2312  {
2313  if (Verbosity() > 0)
2314  {
2315  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeFEMC << endl;
2316  }
2317  // return;
2318  }
2319  if (Verbosity() > 0)
2320  {
2321  cout << "saved\t" << _nTowers_FEMC << "\tFEMC towers" << endl;
2322  }
2323  }
2324  else
2325  {
2326  if (Verbosity() > 0)
2327  {
2328  cout << PHWHERE << " ERROR: Can't find " << towernodeFEMC << endl;
2329  }
2330  // return;
2331  }
2332  }
2333  //----------------------
2334  // TOWERS CEMC
2335  //----------------------
2336  if (_do_CEMC)
2337  {
2338  CaloRawTowerEval* towerevalCEMC = _caloevalstackCEMC->get_rawtower_eval();
2339  _nTowers_CEMC = 0;
2340  string towernodeCEMC = "TOWER_CALIB_CEMC";
2341  RawTowerContainer* towersCEMC = findNode::getClass<RawTowerContainer>(topNode, towernodeCEMC.c_str());
2342  if (towersCEMC)
2343  {
2344  if (Verbosity() > 0)
2345  {
2346  cout << "saving EMC towers" << endl;
2347  }
2348  string towergeomnodeCEMC = "TOWERGEOM_CEMC";
2349  RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeCEMC.c_str());
2350  if (towergeom)
2351  {
2352  if (_do_GEOMETRY && !_geometry_done[kCEMC])
2353  {
2354  RawTowerGeomContainer::ConstRange all_towers = towergeom->get_tower_geometries();
2355  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
2356  it != all_towers.second; ++it)
2357  {
2358  _calo_ID = kCEMC;
2359  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
2360  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
2361  _calo_towers_iL[_calo_towers_N] = -1;
2362  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
2363  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
2364  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
2365  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
2366  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
2367  _calo_towers_N++;
2368  }
2369  _geometry_done[kCEMC] = 1;
2370  _geometry_tree->Fill();
2371  resetGeometryArrays();
2372  }
2373  RawTowerContainer::ConstRange begin_end = towersCEMC->getTowers();
2374  RawTowerContainer::ConstIterator rtiter;
2375  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
2376  {
2377  RawTower* tower = rtiter->second;
2378  if (tower)
2379  {
2380  // min energy cut
2381  if (tower->get_energy() < _reco_e_threshold[kCEMC]) continue;
2382 
2383  _tower_CEMC_iEta[_nTowers_CEMC] = tower->get_bineta();
2384  _tower_CEMC_iPhi[_nTowers_CEMC] = tower->get_binphi();
2385  _tower_CEMC_E[_nTowers_CEMC] = tower->get_energy();
2386 
2387  PHG4Particle* primary = towerevalCEMC->max_truth_primary_particle_by_energy(tower);
2388  if (primary)
2389  {
2390  _tower_CEMC_trueID[_nTowers_CEMC] = primary->get_track_id();
2391  // gflavor = primary->get_pid();
2392  // efromtruth = towerevalCEMC->get_energy_contribution(tower, primary);
2393  }
2394  else
2395  {
2396  _tower_CEMC_trueID[_nTowers_CEMC] = -10;
2397  }
2398  _nTowers_CEMC++;
2399  }
2400  }
2401  }
2402  else
2403  {
2404  if (Verbosity() > 0)
2405  {
2406  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeCEMC << endl;
2407  }
2408  // return;
2409  }
2410  if (Verbosity() > 0)
2411  {
2412  cout << "saved\t" << _nTowers_CEMC << "\tCEMC towers" << endl;
2413  }
2414  }
2415  else
2416  {
2417  if (Verbosity() > 0)
2418  {
2419  cout << PHWHERE << " ERROR: Can't find " << towernodeCEMC << endl;
2420  }
2421  // return;
2422  }
2423  }
2424  //----------------------
2425  // TOWERS EEMC
2426  //----------------------
2427  if (_do_EEMC)
2428  {
2429  CaloRawTowerEval* towerevalEEMC = _caloevalstackEEMC->get_rawtower_eval();
2430  _nTowers_EEMC = 0;
2431  string towernodeEEMC = "TOWER_CALIB_EEMC";
2432  RawTowerContainer* towersEEMC = findNode::getClass<RawTowerContainer>(topNode, towernodeEEMC.c_str());
2433  if (towersEEMC)
2434  {
2435  if (Verbosity() > 0)
2436  {
2437  cout << "saving EMC towers" << endl;
2438  }
2439  string towergeomnodeEEMC = "TOWERGEOM_EEMC";
2440  RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeEEMC.c_str());
2441  if (towergeom)
2442  {
2443  if (_do_GEOMETRY && !_geometry_done[kEEMC])
2444  {
2445  RawTowerGeomContainer::ConstRange all_towers = towergeom->get_tower_geometries();
2446  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
2447  it != all_towers.second; ++it)
2448  {
2449  _calo_ID = kEEMC;
2450  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
2451  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
2452  _calo_towers_iL[_calo_towers_N] = -1;
2453  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
2454  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
2455  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
2456  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
2457  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
2458  _calo_towers_N++;
2459  }
2460  _geometry_done[kEEMC] = 1;
2461  _geometry_tree->Fill();
2462  resetGeometryArrays();
2463  }
2464  RawTowerContainer::ConstRange begin_end = towersEEMC->getTowers();
2465  RawTowerContainer::ConstIterator rtiter;
2466  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
2467  {
2468  RawTower* tower = rtiter->second;
2469  if (tower)
2470  {
2471  // min energy cut
2472  if (tower->get_energy() < _reco_e_threshold[kEEMC]) continue;
2473 
2474  _tower_EEMC_iEta[_nTowers_EEMC] = tower->get_bineta();
2475  _tower_EEMC_iPhi[_nTowers_EEMC] = tower->get_binphi();
2476  _tower_EEMC_E[_nTowers_EEMC] = tower->get_energy();
2477 
2478  PHG4Particle* primary = towerevalEEMC->max_truth_primary_particle_by_energy(tower);
2479  if (primary)
2480  {
2481  _tower_EEMC_trueID[_nTowers_EEMC] = primary->get_track_id();
2482  // gflavor = primary->get_pid();
2483  // efromtruth = towerevalEEMC->get_energy_contribution(tower, primary);
2484  }
2485  else
2486  {
2487  _tower_EEMC_trueID[_nTowers_EEMC] = -10;
2488  }
2489  _nTowers_EEMC++;
2490  }
2491  }
2492  }
2493  else
2494  {
2495  if (Verbosity() > 0)
2496  {
2497  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeEEMC << endl;
2498  }
2499  // return;
2500  }
2501  if (Verbosity() > 0)
2502  {
2503  cout << "saved\t" << _nTowers_EEMC << "\tEEMC towers" << endl;
2504  }
2505  }
2506  else
2507  {
2508  if (Verbosity() > 0)
2509  {
2510  cout << PHWHERE << " ERROR: Can't find " << towernodeEEMC << endl;
2511  }
2512  // return;
2513  }
2514  }
2515 
2516  //----------------------
2517  // TOWERS EEMC
2518  //----------------------
2519  if (_do_EEMCG)
2520  {
2521  CaloRawTowerEval* towerevalEEMCG = _caloevalstackEEMCG->get_rawtower_eval();
2522  _nTowers_EEMCG = 0;
2523  string towernodeEEMCG = "TOWER_CALIB_EEMC_glass";
2524  RawTowerContainer* towersEEMCG = findNode::getClass<RawTowerContainer>(topNode, towernodeEEMCG.c_str());
2525  if (towersEEMCG)
2526  {
2527  if (Verbosity() > 0)
2528  {
2529  cout << "saving EMC towers" << endl;
2530  }
2531  string towergeomnodeEEMCG = "TOWERGEOM_EEMC_glass";
2532  RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeEEMCG.c_str());
2533  if (towergeom)
2534  {
2535  if (_do_GEOMETRY && !_geometry_done[kEEMCG])
2536  {
2537  RawTowerGeomContainer::ConstRange all_towers = towergeom->get_tower_geometries();
2538  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
2539  it != all_towers.second; ++it)
2540  {
2541  _calo_ID = kEEMCG;
2542  _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
2543  _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
2544  _calo_towers_iL[_calo_towers_N] = -1;
2545  _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
2546  _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
2547  _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
2548  _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
2549  _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
2550  _calo_towers_N++;
2551  }
2552  _geometry_done[kEEMCG] = 1;
2553  _geometry_tree->Fill();
2554  resetGeometryArrays();
2555  }
2556  RawTowerContainer::ConstRange begin_end = towersEEMCG->getTowers();
2557  RawTowerContainer::ConstIterator rtiter;
2558  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
2559  {
2560  RawTower* tower = rtiter->second;
2561  if (tower)
2562  {
2563  // min energy cut
2564  if (tower->get_energy() < _reco_e_threshold[kEEMCG]) continue;
2565 
2566  _tower_EEMCG_iEta[_nTowers_EEMCG] = tower->get_bineta();
2567  _tower_EEMCG_iPhi[_nTowers_EEMCG] = tower->get_binphi();
2568  _tower_EEMCG_E[_nTowers_EEMCG] = tower->get_energy();
2569 
2570  PHG4Particle* primary = towerevalEEMCG->max_truth_primary_particle_by_energy(tower);
2571  if (primary)
2572  {
2573  _tower_EEMCG_trueID[_nTowers_EEMCG] = primary->get_track_id();
2574  // gflavor = primary->get_pid();
2575  // efromtruth = towerevalEEMCG->get_energy_contribution(tower, primary);
2576  }
2577  else
2578  {
2579  _tower_EEMCG_trueID[_nTowers_EEMCG] = -10;
2580  }
2581  _nTowers_EEMCG++;
2582  }
2583  }
2584  }
2585  else
2586  {
2587  if (Verbosity() > 0)
2588  {
2589  cout << PHWHERE << " ERROR: Can't find " << towergeomnodeEEMCG << endl;
2590  }
2591  // return;
2592  }
2593  if (Verbosity() > 0)
2594  {
2595  cout << "saved\t" << _nTowers_EEMCG << "\tEEMCG towers" << endl;
2596  }
2597  }
2598  else
2599  {
2600  if (Verbosity() > 0)
2601  {
2602  cout << PHWHERE << " ERROR: Can't find " << towernodeEEMCG << endl;
2603  }
2604  // return;
2605  }
2606  }
2607 
2608  //------------------------
2609  // CLUSTERS FHCAL
2610  //------------------------
2611  if (Verbosity() > 0)
2612  {
2613  cout << "saving clusters" << endl;
2614  }
2615  if (_do_FHCAL && _do_CLUSTERS)
2616  {
2617  CaloRawClusterEval* clusterevalFHCAL = _caloevalstackFHCAL->get_rawcluster_eval();
2618  _nclusters_FHCAL = 0;
2619  if (Verbosity() > 1)
2620  {
2621  cout << "CaloEvaluator::filling gcluster ntuple..." << endl;
2622  }
2623 
2624  string clusternodeFHCAL = "CLUSTER_FHCAL";
2625  RawClusterContainer* clustersFHCAL = findNode::getClass<RawClusterContainer>(topNode, clusternodeFHCAL.c_str());
2626  if (clustersFHCAL)
2627  {
2628  // for every cluster
2629  for (const auto& iterator : clustersFHCAL->getClustersMap())
2630  {
2631  RawCluster* cluster = iterator.second;
2632 
2633  if (cluster->get_energy() < _reco_e_threshold[kFHCAL]) continue;
2634 
2635  _cluster_FHCAL_E[_nclusters_FHCAL] = cluster->get_energy();
2636  _cluster_FHCAL_NTower[_nclusters_FHCAL] = cluster->getNTowers();
2637  _cluster_FHCAL_Phi[_nclusters_FHCAL] = cluster->get_phi();
2638 
2639  // require vertex for cluster eta calculation
2640  if (vertexmap)
2641  {
2642  if (!vertexmap->empty())
2643  {
2644  SvtxVertex* vertex = (vertexmap->begin()->second);
2645  _cluster_FHCAL_Eta[_nclusters_FHCAL] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
2646  }
2647  else
2648  _cluster_FHCAL_Eta[_nclusters_FHCAL] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2649  }
2650  else
2651  _cluster_FHCAL_Eta[_nclusters_FHCAL] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2652  ;
2653 
2654  PHG4Particle* primary = clusterevalFHCAL->max_truth_primary_particle_by_energy(cluster);
2655 
2656  if (primary)
2657  {
2658  _cluster_FHCAL_trueID[_nclusters_FHCAL] = primary->get_track_id();
2659  }
2660  else
2661  {
2662  _cluster_FHCAL_trueID[_nclusters_FHCAL] = -10;
2663  }
2664 
2665  _nclusters_FHCAL++;
2666  }
2667  }
2668  else
2669  {
2670  cerr << PHWHERE << " ERROR: Can't find " << clusternodeFHCAL << endl;
2671  // return;
2672  }
2673  if (Verbosity() > 0)
2674  {
2675  cout << "saved\t" << _nclusters_FHCAL << "\tFHCAL clusters" << endl;
2676  }
2677  }
2678  //------------------------
2679  // CLUSTERS HCALIN
2680  //------------------------
2681  if (_do_HCALIN && _do_CLUSTERS)
2682  {
2683  CaloRawClusterEval* clusterevalHCALIN = _caloevalstackHCALIN->get_rawcluster_eval();
2684  _nclusters_HCALIN = 0;
2685  if (Verbosity() > 1)
2686  {
2687  cout << "CaloEvaluator::filling gcluster ntuple..." << endl;
2688  }
2689 
2690  string clusternodeHCALIN = "CLUSTER_HCALIN";
2691  RawClusterContainer* clustersHCALIN = findNode::getClass<RawClusterContainer>(topNode, clusternodeHCALIN.c_str());
2692  if (clustersHCALIN)
2693  {
2694  // for every cluster
2695  for (const auto& iterator : clustersHCALIN->getClustersMap())
2696  {
2697  RawCluster* cluster = iterator.second;
2698 
2699  if (cluster->get_energy() < _reco_e_threshold[kHCALIN]) continue;
2700 
2701  _cluster_HCALIN_E[_nclusters_HCALIN] = cluster->get_energy();
2702  _cluster_HCALIN_NTower[_nclusters_HCALIN] = cluster->getNTowers();
2703  _cluster_HCALIN_Phi[_nclusters_HCALIN] = cluster->get_phi();
2704 
2705  // require vertex for cluster eta calculation
2706  if (vertexmap)
2707  {
2708  if (!vertexmap->empty())
2709  {
2710  SvtxVertex* vertex = (vertexmap->begin()->second);
2711  _cluster_HCALIN_Eta[_nclusters_HCALIN] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
2712  }
2713  else
2714  _cluster_HCALIN_Eta[_nclusters_HCALIN] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2715  ;
2716  }
2717  else
2718  _cluster_HCALIN_Eta[_nclusters_HCALIN] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2719  ;
2720 
2721  PHG4Particle* primary = clusterevalHCALIN->max_truth_primary_particle_by_energy(cluster);
2722 
2723  if (primary)
2724  {
2725  _cluster_HCALIN_trueID[_nclusters_HCALIN] = primary->get_track_id();
2726  }
2727  else
2728  {
2729  _cluster_HCALIN_trueID[_nclusters_HCALIN] = -10;
2730  }
2731 
2732  _nclusters_HCALIN++;
2733  }
2734  }
2735  else
2736  {
2737  cerr << PHWHERE << " ERROR: Can't find " << clusternodeHCALIN << endl;
2738  // return;
2739  }
2740  if (Verbosity() > 0)
2741  {
2742  cout << "saved\t" << _nclusters_HCALIN << "\tHCALIN clusters" << endl;
2743  }
2744  }
2745  //------------------------
2746  // CLUSTERS HCALOUT
2747  //------------------------
2748  if (_do_HCALOUT && _do_CLUSTERS)
2749  {
2750  CaloRawClusterEval* clusterevalHCALOUT = _caloevalstackHCALOUT->get_rawcluster_eval();
2751  _nclusters_HCALOUT = 0;
2752  if (Verbosity() > 1)
2753  {
2754  cout << "CaloEvaluator::filling gcluster ntuple..." << endl;
2755  }
2756 
2757  string clusternodeHCALOUT = "CLUSTER_HCALOUT";
2758  RawClusterContainer* clustersHCALOUT = findNode::getClass<RawClusterContainer>(topNode, clusternodeHCALOUT.c_str());
2759  if (clustersHCALOUT)
2760  {
2761  // for every cluster
2762  for (const auto& iterator : clustersHCALOUT->getClustersMap())
2763  {
2764  RawCluster* cluster = iterator.second;
2765 
2766  if (cluster->get_energy() < _reco_e_threshold[kHCALOUT]) continue;
2767 
2768  _cluster_HCALOUT_E[_nclusters_HCALOUT] = cluster->get_energy();
2769  _cluster_HCALOUT_NTower[_nclusters_HCALOUT] = cluster->getNTowers();
2770  _cluster_HCALOUT_Phi[_nclusters_HCALOUT] = cluster->get_phi();
2771 
2772  // require vertex for cluster eta calculation
2773  if (vertexmap)
2774  {
2775  if (!vertexmap->empty())
2776  {
2777  SvtxVertex* vertex = (vertexmap->begin()->second);
2778  _cluster_HCALOUT_Eta[_nclusters_HCALOUT] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
2779  }
2780  else
2781  _cluster_HCALOUT_Eta[_nclusters_HCALOUT] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2782  ;
2783  }
2784  else
2785  _cluster_HCALOUT_Eta[_nclusters_HCALOUT] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2786  ;
2787 
2788  PHG4Particle* primary = clusterevalHCALOUT->max_truth_primary_particle_by_energy(cluster);
2789 
2790  if (primary)
2791  {
2792  _cluster_HCALOUT_trueID[_nclusters_HCALOUT] = primary->get_track_id();
2793  }
2794  else
2795  {
2796  _cluster_HCALOUT_trueID[_nclusters_HCALOUT] = -10;
2797  }
2798 
2799  _nclusters_HCALOUT++;
2800  }
2801  }
2802  else
2803  {
2804  cerr << PHWHERE << " ERROR: Can't find " << clusternodeHCALOUT << endl;
2805  // return;
2806  }
2807  if (Verbosity() > 0)
2808  {
2809  cout << "saved\t" << _nclusters_HCALOUT << "\tHCALOUT clusters" << endl;
2810  }
2811  }
2812  //------------------------
2813  // CLUSTERS EHCAL
2814  //------------------------
2815  if (_do_EHCAL && _do_CLUSTERS)
2816  {
2817  CaloRawClusterEval* clusterevalEHCAL = _caloevalstackEHCAL->get_rawcluster_eval();
2818  _nclusters_EHCAL = 0;
2819  if (Verbosity() > 1)
2820  {
2821  cout << "CaloEvaluator::filling gcluster ntuple..." << endl;
2822  }
2823 
2824  string clusternodeEHCAL = "CLUSTER_EHCAL";
2825  RawClusterContainer* clustersEHCAL = findNode::getClass<RawClusterContainer>(topNode, clusternodeEHCAL.c_str());
2826  if (clustersEHCAL)
2827  {
2828  // for every cluster
2829  for (const auto& iterator : clustersEHCAL->getClustersMap())
2830  {
2831  RawCluster* cluster = iterator.second;
2832 
2833  if (cluster->get_energy() < _reco_e_threshold[kEHCAL]) continue;
2834 
2835  _cluster_EHCAL_E[_nclusters_EHCAL] = cluster->get_energy();
2836  _cluster_EHCAL_NTower[_nclusters_EHCAL] = cluster->getNTowers();
2837  _cluster_EHCAL_Phi[_nclusters_EHCAL] = cluster->get_phi();
2838 
2839  // require vertex for cluster eta calculation
2840  if (vertexmap)
2841  {
2842  if (!vertexmap->empty())
2843  {
2844  SvtxVertex* vertex = (vertexmap->begin()->second);
2845  _cluster_EHCAL_Eta[_nclusters_EHCAL] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
2846  }
2847  else
2848  _cluster_EHCAL_Eta[_nclusters_EHCAL] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2849  ;
2850  }
2851  else
2852  _cluster_EHCAL_Eta[_nclusters_EHCAL] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2853  ;
2854 
2855  PHG4Particle* primary = clusterevalEHCAL->max_truth_primary_particle_by_energy(cluster);
2856 
2857  if (primary)
2858  {
2859  _cluster_EHCAL_trueID[_nclusters_EHCAL] = primary->get_track_id();
2860  }
2861  else
2862  {
2863  _cluster_EHCAL_trueID[_nclusters_EHCAL] = -10;
2864  }
2865 
2866  _nclusters_EHCAL++;
2867  }
2868  }
2869  else
2870  {
2871  cerr << PHWHERE << " ERROR: Can't find " << clusternodeEHCAL << endl;
2872  // return;
2873  }
2874  if (Verbosity() > 0)
2875  {
2876  cout << "saved\t" << _nclusters_EHCAL << "\tEHCAL clusters" << endl;
2877  }
2878  }
2879 
2880  //------------------------
2881  // CLUSTERS FEMC
2882  //------------------------
2883  if (_do_FEMC && _do_CLUSTERS)
2884  {
2885  CaloRawClusterEval* clusterevalFEMC = _caloevalstackFEMC->get_rawcluster_eval();
2886  _nclusters_FEMC = 0;
2887  if (Verbosity() > 1)
2888  {
2889  cout << "CaloEvaluator::filling gcluster ntuple..." << endl;
2890  }
2891 
2892  string clusternodeFEMC = "CLUSTER_FEMC";
2893  RawClusterContainer* clustersFEMC = findNode::getClass<RawClusterContainer>(topNode, clusternodeFEMC.c_str());
2894  if (clustersFEMC)
2895  {
2896  // for every cluster
2897  for (const auto& iterator : clustersFEMC->getClustersMap())
2898  {
2899  RawCluster* cluster = iterator.second;
2900 
2901  if (cluster->get_energy() < _reco_e_threshold[kFEMC]) continue;
2902 
2903  _cluster_FEMC_E[_nclusters_FEMC] = cluster->get_energy();
2904  _cluster_FEMC_NTower[_nclusters_FEMC] = cluster->getNTowers();
2905  _cluster_FEMC_Phi[_nclusters_FEMC] = cluster->get_phi();
2906 
2907  // require vertex for cluster eta calculation
2908  if (vertexmap)
2909  {
2910  if (!vertexmap->empty())
2911  {
2912  SvtxVertex* vertex = (vertexmap->begin()->second);
2913  _cluster_FEMC_Eta[_nclusters_FEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
2914  }
2915  else
2916  _cluster_FEMC_Eta[_nclusters_FEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2917  ;
2918  }
2919  else
2920  _cluster_FEMC_Eta[_nclusters_FEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2921  ;
2922 
2923  PHG4Particle* primary = clusterevalFEMC->max_truth_primary_particle_by_energy(cluster);
2924 
2925  if (primary)
2926  {
2927  _cluster_FEMC_trueID[_nclusters_FEMC] = primary->get_track_id();
2928  }
2929  else
2930  {
2931  _cluster_FEMC_trueID[_nclusters_FEMC] = -10;
2932  }
2933 
2934  _nclusters_FEMC++;
2935  }
2936  }
2937  else
2938  {
2939  cerr << PHWHERE << " ERROR: Can't find " << clusternodeFEMC << endl;
2940  // return;
2941  }
2942  if (Verbosity() > 0)
2943  {
2944  cout << "saved\t" << _nclusters_FEMC << "\tFEMC clusters" << endl;
2945  }
2946  }
2947  //------------------------
2948  // CLUSTERS CEMC
2949  //------------------------
2950  if (_do_CEMC && _do_CLUSTERS)
2951  {
2952  CaloRawClusterEval* clusterevalCEMC = _caloevalstackCEMC->get_rawcluster_eval();
2953  _nclusters_CEMC = 0;
2954  if (Verbosity() > 1)
2955  {
2956  cout << "CaloEvaluator::filling gcluster ntuple..." << endl;
2957  }
2958 
2959  string clusternodeCEMC = "CLUSTER_CEMC";
2960  RawClusterContainer* clustersCEMC = findNode::getClass<RawClusterContainer>(topNode, clusternodeCEMC.c_str());
2961  if (clustersCEMC)
2962  {
2963  // for every cluster
2964  for (const auto& iterator : clustersCEMC->getClustersMap())
2965  {
2966  RawCluster* cluster = iterator.second;
2967 
2968  if (cluster->get_energy() < _reco_e_threshold[kCEMC]) continue;
2969 
2970  _cluster_CEMC_E[_nclusters_CEMC] = cluster->get_energy();
2971  _cluster_CEMC_NTower[_nclusters_CEMC] = cluster->getNTowers();
2972  _cluster_CEMC_Phi[_nclusters_CEMC] = cluster->get_phi();
2973 
2974  // require vertex for cluster eta calculation
2975  if (vertexmap)
2976  {
2977  if (!vertexmap->empty())
2978  {
2979  SvtxVertex* vertex = (vertexmap->begin()->second);
2980  _cluster_CEMC_Eta[_nclusters_CEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
2981  }
2982  else
2983  _cluster_CEMC_Eta[_nclusters_CEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2984  ;
2985  }
2986  else
2987  _cluster_CEMC_Eta[_nclusters_CEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
2988  ;
2989 
2990  PHG4Particle* primary = clusterevalCEMC->max_truth_primary_particle_by_energy(cluster);
2991 
2992  if (primary)
2993  {
2994  _cluster_CEMC_trueID[_nclusters_CEMC] = primary->get_track_id();
2995  }
2996  else
2997  {
2998  _cluster_CEMC_trueID[_nclusters_CEMC] = -10;
2999  }
3000 
3001  _nclusters_CEMC++;
3002  }
3003  }
3004  else
3005  {
3006  cerr << PHWHERE << " ERROR: Can't find " << clusternodeCEMC << endl;
3007  // return;
3008  }
3009  if (Verbosity() > 0)
3010  {
3011  cout << "saved\t" << _nclusters_CEMC << "\tCEMC clusters" << endl;
3012  }
3013  }
3014  //------------------------
3015  // CLUSTERS EEMC
3016  //------------------------
3017  if (_do_EEMC && _do_CLUSTERS)
3018  {
3019  CaloRawClusterEval* clusterevalEEMC = _caloevalstackEEMC->get_rawcluster_eval();
3020  _nclusters_EEMC = 0;
3021  if (Verbosity() > 1)
3022  {
3023  cout << "CaloEvaluator::filling gcluster ntuple..." << endl;
3024  }
3025 
3026  string clusternodeEEMC = "CLUSTER_EEMC";
3027  RawClusterContainer* clustersEEMC = findNode::getClass<RawClusterContainer>(topNode, clusternodeEEMC.c_str());
3028  if (clustersEEMC)
3029  {
3030  // for every cluster
3031  for (const auto& iterator : clustersEEMC->getClustersMap())
3032  {
3033  RawCluster* cluster = iterator.second;
3034 
3035  if (cluster->get_energy() < _reco_e_threshold[kEEMC]) continue;
3036 
3037  _cluster_EEMC_E[_nclusters_EEMC] = cluster->get_energy();
3038  _cluster_EEMC_NTower[_nclusters_EEMC] = cluster->getNTowers();
3039  _cluster_EEMC_Phi[_nclusters_EEMC] = cluster->get_phi();
3040 
3041  // require vertex for cluster eta calculation
3042  if (vertexmap)
3043  {
3044  if (!vertexmap->empty())
3045  {
3046  SvtxVertex* vertex = (vertexmap->begin()->second);
3047  _cluster_EEMC_Eta[_nclusters_EEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
3048  }
3049  else
3050  _cluster_EEMC_Eta[_nclusters_EEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
3051  ;
3052  }
3053  else
3054  _cluster_EEMC_Eta[_nclusters_EEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
3055  ;
3056 
3057  PHG4Particle* primary = clusterevalEEMC->max_truth_primary_particle_by_energy(cluster);
3058 
3059  if (primary)
3060  {
3061  _cluster_EEMC_trueID[_nclusters_EEMC] = primary->get_track_id();
3062  }
3063  else
3064  {
3065  _cluster_EEMC_trueID[_nclusters_EEMC] = -10;
3066  }
3067 
3068  _nclusters_EEMC++;
3069  }
3070  }
3071  else
3072  {
3073  cerr << PHWHERE << " ERROR: Can't find " << clusternodeEEMC << endl;
3074  // return;
3075  }
3076  if (Verbosity() > 0)
3077  {
3078  cout << "saved\t" << _nclusters_EEMC << "\tEEMC clusters" << endl;
3079  }
3080  }
3081  //------------------------
3082  // CLUSTERS EEMCG
3083  //------------------------
3084  if (_do_EEMCG && _do_CLUSTERS)
3085  {
3086  CaloRawClusterEval* clusterevalEEMCG = _caloevalstackEEMCG->get_rawcluster_eval();
3087  _nclusters_EEMCG = 0;
3088  if (Verbosity() > 1)
3089  {
3090  cout << "CaloEvaluator::filling gcluster ntuple..." << endl;
3091  }
3092 
3093  string clusternodeEEMCG = "CLUSTER_EEMC_glass";
3094  RawClusterContainer* clustersEEMCG = findNode::getClass<RawClusterContainer>(topNode, clusternodeEEMCG.c_str());
3095  if (clustersEEMCG)
3096  {
3097  // for every cluster
3098  for (const auto& iterator : clustersEEMCG->getClustersMap())
3099  {
3100  RawCluster* cluster = iterator.second;
3101 
3102  if (cluster->get_energy() < _reco_e_threshold[kEEMCG]) continue;
3103 
3104  _cluster_EEMCG_E[_nclusters_EEMCG] = cluster->get_energy();
3105  _cluster_EEMCG_NTower[_nclusters_EEMCG] = cluster->getNTowers();
3106  _cluster_EEMCG_Phi[_nclusters_EEMCG] = cluster->get_phi();
3107 
3108  // require vertex for cluster eta calculation
3109  if (vertexmap)
3110  {
3111  if (!vertexmap->empty())
3112  {
3113  SvtxVertex* vertex = (vertexmap->begin()->second);
3114  _cluster_EEMCG_Eta[_nclusters_EEMCG] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
3115  }
3116  else
3117  _cluster_EEMCG_Eta[_nclusters_EEMCG] = -10000;
3118  }
3119  else
3120  _cluster_EEMCG_Eta[_nclusters_EEMCG] = -10000;
3121 
3122  PHG4Particle* primary = clusterevalEEMCG->max_truth_primary_particle_by_energy(cluster);
3123 
3124  if (primary)
3125  {
3126  _cluster_EEMCG_trueID[_nclusters_EEMCG] = primary->get_track_id();
3127  }
3128  else
3129  {
3130  _cluster_EEMCG_trueID[_nclusters_EEMCG] = -10;
3131  }
3132 
3133  _nclusters_EEMCG++;
3134  }
3135  }
3136  else
3137  {
3138  cerr << PHWHERE << " ERROR: Can't find " << clusternodeEEMCG << endl;
3139  // return;
3140  }
3141  if (Verbosity() > 0)
3142  {
3143  cout << "saved\t" << _nclusters_EEMCG << "\tEEMCG clusters" << endl;
3144  }
3145  }
3146 
3147  //------------------------
3148  // TRACKS
3149  //------------------------
3150  if (_do_TRACKS)
3151  {
3152  _nTracks = 0;
3153  _nProjections = 0;
3154 
3155  EICPIDParticleContainer* pidcontainer(nullptr);
3156 
3157  if (_do_PID_LogLikelihood)
3158  {
3159  pidcontainer = findNode::getClass<EICPIDParticleContainer>(topNode, "EICPIDParticleMap");
3160  if (pidcontainer == nullptr)
3161  {
3162  cout << __PRETTY_FUNCTION__ << " Error: missing EICPIDParticleMap while _do_PID_LogLikelihood = "
3163  << _do_PID_LogLikelihood << endl;
3164  }
3165  }
3166 
3167  // Loop over track maps, identifiy each source.
3168  // Although this configuration is fixed here, it doesn't require multiple sources.
3169  // It will only store them if they're available.
3170  std::vector<std::pair<std::string, TrackSource_t>> trackMapPairs = {
3171  {"TrackMap", TrackSource_t::all},
3172  {"InnerTrackMap", TrackSource_t::inner},
3173  {"SiliconTrackMap", TrackSource_t::silicon},
3174  {"TTLTrackMap", TrackSource_t::ttl}
3175  };
3176  bool foundAtLeastOneTrackSource = false;
3177  for (const auto& trackMapInfo : trackMapPairs)
3178  {
3179  if (_nTracks >= _maxNTracks) break;
3180 
3181  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, trackMapInfo.first);
3182  if (trackmap)
3183  {
3184  foundAtLeastOneTrackSource = true;
3185  int nTracksInASource = 0;
3186  if (Verbosity() > 0)
3187  {
3188  cout << "saving tracks for track map: " << trackMapInfo.first << endl;
3189  }
3190  for (SvtxTrackMap::ConstIter track_itr = trackmap->begin(); track_itr != trackmap->end(); track_itr++)
3191  {
3192  if (_nTracks >= _maxNTracks) break;
3193 
3194  SvtxTrack_FastSim* track = dynamic_cast<SvtxTrack_FastSim*>(track_itr->second);
3195  if (track)
3196  {
3197  _track_ID[_nTracks] = track->get_id();
3198  _track_charge[_nTracks] = track->get_charge();
3199  _track_px[_nTracks] = track->get_px();
3200  _track_py[_nTracks] = track->get_py();
3201  _track_pz[_nTracks] = track->get_pz();
3202  _track_x[_nTracks] = track->get_x();
3203  _track_y[_nTracks] = track->get_y();
3204  _track_z[_nTracks] = track->get_z();
3205  _track_ndf[_nTracks] = track->get_ndf();
3206  _track_chi2[_nTracks] = track->get_chisq();
3207  // Ideally, would be dca3d_xy and dca3d_z, but these don't seem to be calculated properly in the
3208  // current (June 2021) simulations (they return NaN). So we take dca (seems to be ~ the 3d distance)
3209  // and dca_2d (seems to be ~ the distance in the transverse plane).
3210  // The names of the branches are based on the method names.
3211  _track_dca[_nTracks] = static_cast<float>(track->get_dca());
3212  _track_dca_2d[_nTracks] = static_cast<float>(track->get_dca2d());
3213  _track_trueID[_nTracks] = track->get_truth_track_id();
3214  _track_source[_nTracks] = static_cast<unsigned short>(trackMapInfo.second);
3215  if (_do_PROJECTIONS)
3216  {
3217  // find projections
3218  for (SvtxTrack::ConstStateIter trkstates = track->begin_states(); trkstates != track->end_states(); ++trkstates)
3219  {
3220  if (Verbosity() > 1)
3221  {
3222  cout << __PRETTY_FUNCTION__ << " processing " << trkstates->second->get_name() << endl;
3223  }
3224  string trackStateName = trkstates->second->get_name();
3225  if (Verbosity() > 1)
3226  {
3227  cout << __PRETTY_FUNCTION__ << " found " << trkstates->second->get_name() << endl;
3228  }
3229  int trackStateIndex = GetProjectionIndex(trackStateName);
3230  if (trackStateIndex > -1)
3231  {
3232  // save true projection info to given branch
3233  _track_TLP_x[_nProjections] = trkstates->second->get_pos(0);
3234  _track_TLP_y[_nProjections] = trkstates->second->get_pos(1);
3235  _track_TLP_z[_nProjections] = trkstates->second->get_pos(2);
3236  _track_TLP_t[_nProjections] = trkstates->first;
3237  _track_TLP_px[_nProjections] = trkstates->second->get_px();
3238  _track_TLP_py[_nProjections] = trkstates->second->get_py();
3239  _track_TLP_pz[_nProjections] = trkstates->second->get_pz();
3240  _track_ProjLayer[_nProjections] = trackStateIndex;
3241  _track_ProjTrackID[_nProjections] = _nTracks;
3242 
3243  string nodename = "G4HIT_" + trkstates->second->get_name();
3244  PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
3245  if (hits)
3246  {
3247  if (Verbosity() > 1)
3248  {
3249  cout << __PRETTY_FUNCTION__ << " number of hits: " << hits->size() << endl;
3250  }
3251  PHG4HitContainer::ConstRange hit_range = hits->getHits();
3252  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
3253  {
3254  if (Verbosity() > 1)
3255  {
3256  cout << __PRETTY_FUNCTION__ << " checking hit id " << hit_iter->second->get_trkid() << " against " << track->get_truth_track_id() << endl;
3257  }
3258  if (hit_iter->second->get_trkid() - track->get_truth_track_id() == 0)
3259  {
3260  if (Verbosity() > 1)
3261  {
3262  cout << __PRETTY_FUNCTION__ << " found hit with id " << hit_iter->second->get_trkid() << endl;
3263  }
3264  // save reco projection info to given branch
3265  _track_TLP_true_x[_nProjections] = hit_iter->second->get_x(0);
3266  _track_TLP_true_y[_nProjections] = hit_iter->second->get_y(0);
3267  _track_TLP_true_z[_nProjections] = hit_iter->second->get_z(0);
3268  _track_TLP_true_t[_nProjections] = hit_iter->second->get_t(0);
3269  }
3270  }
3271  }
3272  else
3273  {
3274  if (Verbosity() > 1)
3275  {
3276  cout << __PRETTY_FUNCTION__ << " could not find " << nodename << endl;
3277  }
3278  _track_TLP_true_x[_nProjections] = -10000;
3279  _track_TLP_true_y[_nProjections] = -10000;
3280  _track_TLP_true_z[_nProjections] = -10000;
3281  _track_TLP_true_t[_nProjections] = -10000;
3282  }
3283  _nProjections++;
3284  }
3285  }
3286  } // if (_do_PROJECTIONS)
3287 
3288  if (_do_PID_LogLikelihood)
3289  {
3290  // perform PID matching
3291  if (trackMapInfo.second == TrackSource_t::all and pidcontainer != nullptr)
3292  {
3293  // only do so for the TrackSource_t::all
3294  const EICPIDParticle* pid_particle =
3295  pidcontainer->findEICPIDParticle(track->get_id());
3296 
3297  if (pid_particle)
3298  {
3299  if ((unsigned int) _nTracks >= _track_pion_LL.size())
3300  {
3301  cout << __PRETTY_FUNCTION__
3302  << " logical error _nTracks = " << _nTracks
3303  << " logical error _track_pion_LL.size() = " << _track_pion_LL.size()
3304  << endl;
3305  exit(1);
3306  }
3307 
3308  _track_pion_LL[_nTracks] = pid_particle->get_SumLogLikelyhood(EICPIDDefs::PionCandiate);
3309  _track_kaon_LL[_nTracks] = pid_particle->get_SumLogLikelyhood(EICPIDDefs::KaonCandiate);
3310  _track_proton_LL[_nTracks] = pid_particle->get_SumLogLikelyhood(EICPIDDefs::ProtonCandiate);
3311  }
3312  }
3313  } // if (_do_PID_LogLikelihood )
3314 
3315  _nTracks++;
3316  nTracksInASource++;
3317  }
3318  else
3319  {
3320  if (Verbosity() > 0) //Verbosity()
3321  {
3322  cout << "PHG4TrackFastSimEval::fill_track_tree - ignore track that is not a SvtxTrack_FastSim:";
3323  track_itr->second->identify();
3324  }
3325  continue;
3326  }
3327  }
3328  if (Verbosity() > 0)
3329  {
3330  cout << "saved\t" << nTracksInASource << "\ttracks from track map " << trackMapInfo.first << ". Total saved tracks: " << _nTracks << endl;
3331  }
3332  }
3333  else
3334  {
3335  if (Verbosity() > 0)
3336  {
3337  cout << PHWHERE << "SvtxTrackMap node with name '" << trackMapInfo.first << "' not found on node tree" << endl;
3338  }
3339  }
3340  }
3341  if (foundAtLeastOneTrackSource == false)
3342  {
3343  cout << PHWHERE << "Requested tracks, but found no sources on node tree. Returning" << endl;
3344  return;
3345  }
3346  }
3347  //------------------------
3348  // MC PARTICLES
3349  //------------------------
3350  _nMCPart = 0;
3351  if (_do_MCPARTICLES)
3352  {
3353  PHG4TruthInfoContainer* truthinfocontainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
3354  if (truthinfocontainer)
3355  {
3356  if (Verbosity() > 0)
3357  {
3358  cout << "saving MC particles" << endl;
3359  }
3360  //GetParticleRange for all particles
3361  //GetPrimaryParticleRange for primary particles
3362  PHG4TruthInfoContainer::ConstRange range = truthinfocontainer->GetParticleRange();
3363  for (PHG4TruthInfoContainer::ConstIterator truth_itr = range.first; truth_itr != range.second; ++truth_itr)
3364  {
3365  PHG4Particle* g4particle = truth_itr->second;
3366  if (!g4particle) continue;
3367 
3368  int mcSteps = 0;
3369  PHG4Particle* g4particleMother = truth_itr->second;
3370  if (g4particle->get_parent_id() != 0)
3371  {
3372  while (g4particleMother->get_parent_id() != 0)
3373  {
3374  g4particleMother = truthinfocontainer->GetParticle(g4particleMother->get_parent_id());
3375  if (g4particleMother == NULL) break;
3376  mcSteps += 1;
3377  }
3378  }
3379  if (mcSteps > _depth_MCstack) continue;
3380 
3381  // evaluating true primary vertex
3382  if (_do_VERTEX && _nMCPart == 0)
3383  {
3384  PHG4VtxPoint* vtx = truthinfocontainer->GetVtx(g4particle->get_vtx_id());
3385  if (vtx)
3386  {
3387  _vertex_true_x = vtx->get_x();
3388  _vertex_true_y = vtx->get_y();
3389  _vertex_true_z = vtx->get_z();
3390  }
3391  }
3392 
3393  // in case of all MC particles, make restrictions on the secondary selection
3394  // if(g4particle->get_track_id()<0 && g4particle->get_e()<0.5) continue;
3395  // primary (g4particle->get_parent_id() == 0) selection via:
3396  // if(gtrackID < 0) continue;
3397 
3398  //using the e threshold also for the truth particles gets rid of all the low energy secondary particles
3399  if (g4particle->get_e() < _reco_e_thresholdMC) continue;
3400 
3401  _mcpart_ID[_nMCPart] = g4particle->get_track_id();
3402  _mcpart_ID_parent[_nMCPart] = g4particle->get_parent_id();
3403  _mcpart_PDG[_nMCPart] = g4particle->get_pid();
3404  _mcpart_E[_nMCPart] = g4particle->get_e();
3405  _mcpart_px[_nMCPart] = g4particle->get_px();
3406  _mcpart_py[_nMCPart] = g4particle->get_py();
3407  _mcpart_pz[_nMCPart] = g4particle->get_pz();
3408  PHG4VtxPoint* vtxtmp = truthinfocontainer->GetVtx(g4particle->get_vtx_id());
3409  if (vtxtmp)
3410  {
3411  _mcpart_x[_nMCPart] = vtxtmp->get_x();
3412  _mcpart_y[_nMCPart] = vtxtmp->get_y();
3413  _mcpart_z[_nMCPart] = vtxtmp->get_z();
3414  }
3415  //BCID added for G4Particle -- HEPMC particle matching
3416  _mcpart_BCID[_nMCPart] = g4particle->get_barcode();
3417  // TVector3 projvec(_mcpart_px[0],_mcpart_py[0],_mcpart_pz[0]);
3418  // float projeta = projvec.Eta();
3419  _nMCPart++;
3420  }
3421  if (Verbosity() > 0)
3422  {
3423  cout << "saved\t" << _nMCPart << "\tMC particles" << endl;
3424  }
3425  }
3426  else
3427  {
3428  if (Verbosity() > 0)
3429  {
3430  cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree" << endl;
3431  }
3432  return;
3433  }
3434  }
3435 
3436  //------------------------
3437  // HEPMC
3438  //------------------------
3439  _nHepmcp = 0;
3440  if (_do_HEPMC)
3441  {
3442  PHHepMCGenEventMap* hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
3443  if (hepmceventmap)
3444  {
3445  if (Verbosity() > 0)
3446  {
3447  cout << "saving HepMC output" << endl;
3448  }
3449  if (Verbosity() > 0)
3450  {
3451  hepmceventmap->Print();
3452  }
3453 
3454  for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
3455  eventIter != hepmceventmap->end();
3456  ++eventIter)
3457  {
3458  PHHepMCGenEvent* hepmcevent = eventIter->second;
3459 
3460  if (hepmcevent)
3461  {
3462  HepMC::GenEvent* truthevent = hepmcevent->getEvent();
3463  if (!truthevent)
3464  {
3465  cout << PHWHERE
3466  << "no evt pointer under phhepmvgeneventmap found "
3467  << endl;
3468  return;
3469  }
3470 
3471  _hepmcp_vtx_x = hepmcevent->get_collision_vertex().x();
3472  _hepmcp_vtx_y = hepmcevent->get_collision_vertex().y();
3473  _hepmcp_vtx_z = hepmcevent->get_collision_vertex().z();
3474  _hepmcp_vtx_t = hepmcevent->get_collision_vertex().t();
3475 
3476  HepMC::PdfInfo* pdfinfo = truthevent->pdf_info();
3477 
3478  // m_partid1 = pdfinfo->id1();
3479  // m_partid2 = pdfinfo->id2();
3480  _hepmcp_x1 = pdfinfo->x1();
3481  _hepmcp_x2 = pdfinfo->x2();
3482  _hepmcp_Q2 = pdfinfo->scalePDF();
3483 
3484  // m_mpi = truthevent->mpi();
3485 
3486  _hepmcp_procid = truthevent->signal_process_id();
3487 
3488  if (Verbosity() > 2)
3489  {
3490  cout << " Iterating over an event" << endl;
3491  }
3492  for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
3493  iter != truthevent->particles_end();
3494  ++iter)
3495  {
3496  _hepmcp_E[_nHepmcp] = (*iter)->momentum().e();
3497  _hepmcp_PDG[_nHepmcp] = (*iter)->pdg_id();
3498  _hepmcp_px[_nHepmcp] = (*iter)->momentum().px();
3499  _hepmcp_py[_nHepmcp] = (*iter)->momentum().py();
3500  _hepmcp_pz[_nHepmcp] = (*iter)->momentum().pz();
3501  _hepmcp_status[_nHepmcp] = (*iter)->status();
3502  _hepmcp_BCID[_nHepmcp] = (*iter)->barcode();
3503  _hepmcp_m2[_nHepmcp] = 0;
3504  _hepmcp_m1[_nHepmcp] = 0;
3505  if ((*iter)->production_vertex())
3506  {
3507  for (HepMC::GenVertex::particle_iterator mother = (*iter)->production_vertex()->particles_begin(HepMC::parents);
3508  mother != (*iter)->production_vertex()->particles_end(HepMC::parents);
3509  ++mother)
3510  {
3511  _hepmcp_m2[_nHepmcp] = (*mother)->barcode();
3512  if (_hepmcp_m1[_nHepmcp] == 0)
3513  _hepmcp_m1[_nHepmcp] = (*mother)->barcode();
3514  }
3515  }
3516  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;
3517  _nHepmcp++;
3518  }
3519  }
3520  }
3521  }
3522  else
3523  {
3524  if (Verbosity() > 0)
3525  {
3526  cout << PHWHERE << " PHHepMCGenEventMap node not found on node tree" << endl;
3527  }
3528  return;
3529  }
3530  } //hepmc
3531 
3532  _event_tree->Fill();
3533 
3534  if (Verbosity() > 0)
3535  {
3536  cout << "Resetting buffer ..." << endl;
3537  }
3538  resetBuffer();
3539  if (Verbosity() > 0)
3540  {
3541  cout << "EventEvaluatorEIC buffer reset" << endl;
3542  }
3543  return;
3544 }
3545 
3546 int EventEvaluatorEIC::End(PHCompositeNode* topNode)
3547 {
3548  _tfile->cd();
3549 
3550  _event_tree->Write();
3551 
3552  _tfile->Close();
3553 
3554  delete _tfile;
3555 
3556  if (_do_GEOMETRY)
3557  {
3558  _tfile_geometry->cd();
3559 
3560  _geometry_tree->Write();
3561 
3562  _tfile_geometry->Close();
3563 
3564  delete _tfile_geometry;
3565  }
3566  if (Verbosity() > 0)
3567  {
3568  cout << "========================= " << Name() << "::End() ============================" << endl;
3569  cout << " " << _ievent << " events of output written to: " << _filename << endl;
3570  cout << "===========================================================================" << endl;
3571  }
3572 
3573  if (_caloevalstackFHCAL) delete _caloevalstackFHCAL;
3574  if (_caloevalstackBECAL) delete _caloevalstackBECAL;
3575  if (_caloevalstackHCALIN) delete _caloevalstackHCALIN;
3576  if (_caloevalstackHCALOUT) delete _caloevalstackHCALOUT;
3577  if (_caloevalstackEHCAL) delete _caloevalstackEHCAL;
3578  if (_caloevalstackDRCALO) delete _caloevalstackDRCALO;
3579  if (_caloevalstackFOCAL) delete _caloevalstackFOCAL;
3580  if (_caloevalstackLFHCAL) delete _caloevalstackLFHCAL;
3581  if (_caloevalstackFEMC) delete _caloevalstackFEMC;
3582  if (_caloevalstackCEMC) delete _caloevalstackCEMC;
3583  if (_caloevalstackEEMC) delete _caloevalstackEEMC;
3584  if (_caloevalstackEEMCG) delete _caloevalstackEEMCG;
3585 
3586  return Fun4AllReturnCodes::EVENT_OK;
3587 }
3588 
3589 int EventEvaluatorEIC::GetProjectionIndex(std::string projname)
3590 {
3591  if (projname.find("FTTL_0") != std::string::npos)
3592  return 0;
3593  else if (projname.find("FTTL_1") != std::string::npos)
3594  return 1;
3595  else if (projname.find("FTTL_2") != std::string::npos)
3596  return 2;
3597  else if (projname.find("ETTL_0") != std::string::npos)
3598  return 3;
3599  else if (projname.find("ETTL_1") != std::string::npos)
3600  return 4;
3601  else if (projname.find("LFHCAL") != std::string::npos)
3602  return 5;
3603  else if (projname.find("FEMC") != std::string::npos)
3604  return 6;
3605  else if (projname.find("CTTL_0") != std::string::npos)
3606  return 7;
3607  else if (projname.find("CTTL_1") != std::string::npos)
3608  return 8;
3609 
3610  else if (projname.find("LBLVTX_CENTRAL_10") != std::string::npos)
3611  return 10;
3612  else if (projname.find("LBLVTX_CENTRAL_11") != std::string::npos)
3613  return 11;
3614  else if (projname.find("LBLVTX_CENTRAL_12") != std::string::npos)
3615  return 12;
3616  else if (projname.find("LBLVTX_CENTRAL_13") != std::string::npos)
3617  return 13;
3618  else if (projname.find("LBLVTX_CENTRAL_14") != std::string::npos)
3619  return 14;
3620  else if (projname.find("LBLVTX_CENTRAL_15") != std::string::npos)
3621  return 15;
3622 
3623  else if (projname.find("LBLVTX_FORWARD_20") != std::string::npos)
3624  return 20;
3625  else if (projname.find("LBLVTX_FORWARD_21") != std::string::npos)
3626  return 21;
3627  else if (projname.find("LBLVTX_FORWARD_22") != std::string::npos)
3628  return 22;
3629  else if (projname.find("LBLVTX_FORWARD_23") != std::string::npos)
3630  return 23;
3631  else if (projname.find("LBLVTX_FORWARD_24") != std::string::npos)
3632  return 24;
3633 
3634  else if (projname.find("LBLVTX_BACKWARD_30") != std::string::npos)
3635  return 30;
3636  else if (projname.find("LBLVTX_BACKWARD_31") != std::string::npos)
3637  return 31;
3638  else if (projname.find("LBLVTX_BACKWARD_32") != std::string::npos)
3639  return 32;
3640  else if (projname.find("LBLVTX_BACKWARD_33") != std::string::npos)
3641  return 33;
3642  else if (projname.find("LBLVTX_BACKWARD_34") != std::string::npos)
3643  return 34;
3644 
3645  else if (projname.find("BARREL_0") != std::string::npos)
3646  return 40;
3647  else if (projname.find("BARREL_1") != std::string::npos)
3648  return 41;
3649  else if (projname.find("BARREL_2") != std::string::npos)
3650  return 42;
3651  else if (projname.find("BARREL_3") != std::string::npos)
3652  return 43;
3653  else if (projname.find("BARREL_4") != std::string::npos)
3654  return 44;
3655  else if (projname.find("BARREL_5") != std::string::npos)
3656  return 45;
3657 
3658  else if (projname.find("EFST_0") != std::string::npos)
3659  return 100;
3660  else if (projname.find("EFST_1") != std::string::npos)
3661  return 101;
3662  else if (projname.find("EFST_2") != std::string::npos)
3663  return 102;
3664  else if (projname.find("EFST_3") != std::string::npos)
3665  return 103;
3666  else if (projname.find("EFST_4") != std::string::npos)
3667  return 104;
3668  else if (projname.find("EFST_5") != std::string::npos)
3669  return 105;
3670  else if (projname.find("EFST_6") != std::string::npos)
3671  return 106;
3672 
3673  else if (projname.find("FST_0") != std::string::npos)
3674  return 50;
3675  else if (projname.find("FST_1") != std::string::npos)
3676  return 51;
3677  else if (projname.find("FST_2") != std::string::npos)
3678  return 52;
3679  else if (projname.find("FST_3") != std::string::npos)
3680  return 53;
3681  else if (projname.find("FST_4") != std::string::npos)
3682  return 54;
3683  else if (projname.find("FST_5") != std::string::npos)
3684  return 55;
3685  else if (projname.find("FST_6") != std::string::npos)
3686  return 56;
3687 
3688  else if (projname.find("EHCAL") != std::string::npos)
3689  return 60;
3690  else if (projname.find("EEMC") != std::string::npos)
3691  return 61;
3692  else if (projname.find("HCALIN") != std::string::npos)
3693  return 62;
3694  else if (projname.find("HCALOUT") != std::string::npos)
3695  return 63;
3696  else if (projname.find("CEMC") != std::string::npos)
3697  return 64;
3698  else if (projname.find("EEMC_glass_0") != std::string::npos)
3699  return 65;
3700  else if (projname.find("BECAL") != std::string::npos)
3701  return 66;
3702  // else if (projname.find("FHCAL") != std::string::npos)
3703  // return 67;
3704 
3705  else if (projname.find("ZDCsurrogate") != std::string::npos)
3706  return 70;
3707  else if (projname.find("rpTruth") != std::string::npos)
3708  return 71;
3709  else if (projname.find("b0Truth") != std::string::npos)
3710  return 72;
3711  else if (projname.find("rpTruth2") != std::string::npos)
3712  return 73;
3713  else if (projname.find("offMomTruth") != std::string::npos)
3714  return 74;
3715 
3716  else if (projname.find("BH_1") != std::string::npos)
3717  return 90;
3718  else if (projname.find("BH_FORWARD_PLUS") != std::string::npos)
3719  return 91;
3720  else if (projname.find("BH_FORWARD_NEG") != std::string::npos)
3721  return 92;
3722 
3723  else if (projname.find("RWELL_0") != std::string::npos)
3724  return 110;
3725  else if (projname.find("RWELL_1") != std::string::npos)
3726  return 111;
3727  else if (projname.find("RWELL_2") != std::string::npos)
3728  return 112;
3729 
3730 
3731  else if (projname.find("EST_0") != std::string::npos)
3732  return 113;
3733  else if (projname.find("EST_1") != std::string::npos)
3734  return 114;
3735  else if (projname.find("EST_2") != std::string::npos)
3736  return 115;
3737  else if (projname.find("EST_3") != std::string::npos)
3738  return 116;
3739  else if (projname.find("EST_4") != std::string::npos)
3740  return 117;
3741  else if (projname.find("EST_5") != std::string::npos)
3742  return 118;
3743  else if (projname.find("EST_6") != std::string::npos)
3744  return 119;
3745 
3746  else if (projname.find("FGEM_0") != std::string::npos)
3747  return 120;
3748  else if (projname.find("FGEM_1") != std::string::npos)
3749  return 121;
3750 
3751  else if (projname.find("EGEM_0") != std::string::npos)
3752  return 130;
3753  else if (projname.find("EGEM_1") != std::string::npos)
3754  return 131;
3755 
3756 
3757  else if (projname.find("LFHCAL_0") != std::string::npos)
3758  return 140;
3759  else if (projname.find("LFHCAL_1") != std::string::npos)
3760  return 141;
3761  else if (projname.find("LFHCAL_2") != std::string::npos)
3762  return 142;
3763  else if (projname.find("LFHCAL_3") != std::string::npos)
3764  return 143;
3765  else if (projname.find("LFHCAL_4") != std::string::npos)
3766  return 144;
3767  else if (projname.find("LFHCAL_5") != std::string::npos)
3768  return 145;
3769  else if (projname.find("LFHCAL_6") != std::string::npos)
3770  return 146;
3771  else if (projname.find("LFHCAL_7") != std::string::npos)
3772  return 147;
3773 
3774  else if (projname.find("BARR_0") != std::string::npos)
3775  return 150;
3776  else if (projname.find("BARR_1") != std::string::npos)
3777  return 151;
3778  else if (projname.find("BARR_2") != std::string::npos)
3779  return 152;
3780  else if (projname.find("BARR_3") != std::string::npos)
3781  return 153;
3782  else if (projname.find("BARR_4") != std::string::npos)
3783  return 154;
3784  else if (projname.find("SVTX_0") != std::string::npos)
3785  return 155;
3786  else if (projname.find("SVTX_1") != std::string::npos)
3787  return 156;
3788  else if (projname.find("SVTX_2") != std::string::npos)
3789  return 157;
3790  else if (projname.find("SVTX_3") != std::string::npos)
3791  return 158;
3792  else if (projname.find("SVTX_4") != std::string::npos)
3793  return 159;
3794 
3795  else if (projname.find("BARR") != std::string::npos)
3796  return 160;
3797  else if (projname.find("SVTX") != std::string::npos)
3798  return 161;
3799  else
3800  return -1;
3801  return -1;
3802 }
3803 
3804 std::string EventEvaluatorEIC::GetProjectionNameFromIndex(int projindex)
3805 {
3806  switch (projindex)
3807  {
3808  case 0:
3809  return "FTTL_0";
3810  case 1:
3811  return "FTTL_1";
3812  case 2:
3813  return "FTTL_2";
3814  case 3:
3815  return "ETTL_0";
3816  case 4:
3817  return "ETTL_1";
3818  case 5:
3819  return "LFHCAL";
3820  case 6:
3821  return "FEMC";
3822  case 7:
3823  return "CTTL_0";
3824  case 8:
3825  return "CTTL_1";
3826 
3827  case 10:
3828  return "LBLVTX_CENTRAL_10";
3829  case 11:
3830  return "LBLVTX_CENTRAL_11";
3831  case 12:
3832  return "LBLVTX_CENTRAL_12";
3833  case 13:
3834  return "LBLVTX_CENTRAL_13";
3835  case 14:
3836  return "LBLVTX_CENTRAL_14";
3837  case 15:
3838  return "LBLVTX_CENTRAL_15";
3839 
3840  case 20:
3841  return "LBLVTX_FORWARD_20";
3842  case 21:
3843  return "LBLVTX_FORWARD_21";
3844  case 22:
3845  return "LBLVTX_FORWARD_22";
3846  case 23:
3847  return "LBLVTX_FORWARD_23";
3848  case 24:
3849  return "LBLVTX_FORWARD_24";
3850 
3851  case 30:
3852  return "LBLVTX_BACKWARD_30";
3853  case 31:
3854  return "LBLVTX_BACKWARD_31";
3855  case 32:
3856  return "LBLVTX_BACKWARD_32";
3857  case 33:
3858  return "LBLVTX_BACKWARD_33";
3859  case 34:
3860  return "LBLVTX_BACKWARD_34";
3861 
3862  case 40:
3863  return "BARREL_0";
3864  case 41:
3865  return "BARREL_1";
3866  case 42:
3867  return "BARREL_2";
3868  case 43:
3869  return "BARREL_3";
3870  case 44:
3871  return "BARREL_4";
3872  case 45:
3873  return "BARREL_5";
3874 
3875  case 50:
3876  return "FST_0";
3877  case 51:
3878  return "FST_1";
3879  case 52:
3880  return "FST_2";
3881  case 53:
3882  return "FST_3";
3883  case 54:
3884  return "FST_4";
3885  case 55:
3886  return "FST_5";
3887  case 56 :
3888  return "FST_6";
3889 
3890  case 60:
3891  return "EHCAL";
3892  case 61:
3893  return "EEMC";
3894  case 62:
3895  return "HCALIN";
3896  case 63:
3897  return "HCALOUT";
3898  case 64:
3899  return "CEMC";
3900  case 65:
3901  return "EEMC_glass";
3902  case 66:
3903  return "BECAL";
3904  // case 67:
3905  // return "LFHCAL";
3906 
3907  case 70:
3908  return "ZDCsurrogate";
3909  case 71:
3910  return "rpTruth";
3911  case 72:
3912  return "b0Truth";
3913  case 73:
3914  return "rpTruth2";
3915  case 74:
3916  return "offMomTruth";
3917 
3918  case 90:
3919  return "BH_1";
3920  case 91:
3921  return "BH_FORWARD_PLUS";
3922  case 92:
3923  return "BH_FORWARD_NEG";
3924 
3925  case 100:
3926  return "EFST_0";
3927  case 101:
3928  return "EFST_1";
3929  case 102:
3930  return "EFST_2";
3931  case 103:
3932  return "EFST_3";
3933  case 104:
3934  return "EFST_4";
3935  case 105:
3936  return "EFST_5";
3937 
3938  case 110:
3939  return "RWELL_0";
3940  case 111:
3941  return "RWELL_1";
3942  case 112:
3943  return "RWELL_2";
3944 
3945  case 113:
3946  return "EST_0";
3947  case 114:
3948  return "EST_1";
3949  case 115:
3950  return "EST_2";
3951  case 116:
3952  return "EST_3";
3953  case 117:
3954  return "EST_4";
3955  case 118:
3956  return "EST_5";
3957 
3958 
3959  case 120:
3960  return "FGEM_0";
3961  case 121:
3962  return "FGEM_1";
3963 
3964  case 130:
3965  return "EGEM_0";
3966  case 131:
3967  return "EGEM_1";
3968 
3969 
3970  case 140:
3971  return "LFHCAL_0";
3972  case 141:
3973  return "LFHCAL_1";
3974  case 142:
3975  return "LFHCAL_2";
3976  case 143:
3977  return "LFHCAL_3";
3978  case 144:
3979  return "LFHCAL_4";
3980  case 145:
3981  return "LFHCAL_5";
3982  case 146:
3983  return "LFHCAL_6";
3984  case 147:
3985  return "LFHCAL_7";
3986 
3987  case 150:
3988  return "BARR_0";
3989  case 151:
3990  return "BARR_1";
3991  case 152:
3992  return "BARR_2";
3993  case 153:
3994  return "BARR_3";
3995  case 154:
3996  return "BARR_4";
3997  case 155:
3998  return "SVTX_0";
3999  case 156:
4000  return "SVTX_1";
4001  case 157:
4002  return "SVTX_2";
4003  case 158:
4004  return "SVTX_3";
4005  case 159:
4006  return "SVTX_4";
4007 
4008  case 160:
4009  return "BARR";
4010  case 161:
4011  return "SVTX";
4012  default:
4013  return "NOTHING";
4014  }
4015 }
4016 
4017 void EventEvaluatorEIC::resetGeometryArrays()
4018 {
4019  for (Int_t igeo = 0; igeo < _calo_towers_N; igeo++)
4020  {
4021  _calo_towers_iEta[_calo_towers_N] = -10000;
4022  _calo_towers_iPhi[_calo_towers_N] = -10000;
4023  _calo_towers_iL[_calo_towers_N] = -10000;
4024  _calo_towers_Eta[_calo_towers_N] = -10000;
4025  _calo_towers_Phi[_calo_towers_N] = -10000;
4026  _calo_towers_x[_calo_towers_N] = -10000;
4027  _calo_towers_y[_calo_towers_N] = -10000;
4028  _calo_towers_z[_calo_towers_N] = -10000;
4029  }
4030  _calo_ID = -1;
4031  _calo_towers_N = 0;
4032 }
4033 void EventEvaluatorEIC::resetBuffer()
4034 {
4035  if (_do_store_event_info)
4036  {
4037  _cross_section = 0;
4038  _event_weight = 0;
4039  _n_generator_accepted = 0;
4040  if (Verbosity() > 0)
4041  {
4042  cout << "\t... event info variables reset" << endl;
4043  }
4044  }
4045  if (_do_VERTEX)
4046  {
4047  _vertex_x = -1000;
4048  _vertex_y = -1000;
4049  _vertex_z = -1000;
4050  _vertex_NCont = 0;
4051  _vertex_true_x = -1000;
4052  _vertex_true_y = -1000;
4053  _vertex_true_z = -1000;
4054  if (Verbosity() > 0)
4055  {
4056  cout << "\t... vertex variables reset" << endl;
4057  }
4058  }
4059  if (_do_HITS)
4060  {
4061  _nHitsLayers = 0;
4062  for (Int_t ihit = 0; ihit < _maxNHits; ihit++)
4063  {
4064  _hits_layerID[ihit] = 0;
4065  _hits_trueID[ihit] = 0;
4066  _hits_x[ihit] = 0;
4067  _hits_y[ihit] = 0;
4068  _hits_z[ihit] = 0;
4069  _hits_x2[ihit] = 0;
4070  _hits_y2[ihit] = 0;
4071  _hits_z2[ihit] = 0;
4072  _hits_t[ihit] = 0;
4073  _hits_edep[ihit] = 0;
4074  _hits_lightyield[ihit] = 0;
4075  _hits_isAbsorber[ihit] = 0;
4076  }
4077  if (Verbosity() > 0)
4078  {
4079  cout << "\t... hit variables reset" << endl;
4080  }
4081  }
4082  if (_do_FHCAL)
4083  {
4084  _nTowers_FHCAL = 0;
4085  for (Int_t itow = 0; itow < _maxNTowers; itow++)
4086  {
4087  _tower_FHCAL_E[itow] = 0;
4088  _tower_FHCAL_iEta[itow] = 0;
4089  _tower_FHCAL_iPhi[itow] = 0;
4090  _tower_FHCAL_trueID[itow] = 0;
4091  }
4092  if (_do_CLUSTERS)
4093  {
4094  _nclusters_FHCAL = 0;
4095  for (Int_t itow = 0; itow < _maxNclusters; itow++)
4096  {
4097  _cluster_FHCAL_E[itow] = 0;
4098  _cluster_FHCAL_Eta[itow] = 0;
4099  _cluster_FHCAL_Phi[itow] = 0;
4100  _cluster_FHCAL_NTower[itow] = 0;
4101  _cluster_FHCAL_trueID[itow] = 0;
4102  }
4103  }
4104  if (Verbosity() > 0)
4105  {
4106  cout << "\t... FHCAL variables reset" << endl;
4107  }
4108  }
4109  if (_do_BECAL)
4110  {
4111  _nTowers_BECAL = 0;
4112  for (Int_t itow = 0; itow < _maxNTowers; itow++)
4113  {
4114  _tower_BECAL_E[itow] = 0;
4115  _tower_BECAL_iEta[itow] = 0;
4116  _tower_BECAL_iPhi[itow] = 0;
4117  _tower_BECAL_trueID[itow] = 0;
4118  }
4119  if (Verbosity() > 0)
4120  {
4121  cout << "\t... BECAL variables reset" << endl;
4122  }
4123  }
4124  if (_do_FEMC)
4125  {
4126  _nTowers_FEMC = 0;
4127  for (Int_t itow = 0; itow < _maxNTowers; itow++)
4128  {
4129  _tower_FEMC_E[itow] = 0;
4130  _tower_FEMC_iEta[itow] = 0;
4131  _tower_FEMC_iPhi[itow] = 0;
4132  _tower_FEMC_trueID[itow] = 0;
4133  }
4134  if (_do_CLUSTERS)
4135  {
4136  _nclusters_FEMC = 0;
4137  for (Int_t itow = 0; itow < _maxNclusters; itow++)
4138  {
4139  _cluster_FEMC_E[itow] = 0;
4140  _cluster_FEMC_Eta[itow] = 0;
4141  _cluster_FEMC_Phi[itow] = 0;
4142  _cluster_FEMC_NTower[itow] = 0;
4143  _cluster_FEMC_trueID[itow] = 0;
4144  }
4145  }
4146  if (Verbosity() > 0)
4147  {
4148  cout << "\t... FEMC variables reset" << endl;
4149  }
4150  }
4151  if (_do_CEMC)
4152  {
4153  _nTowers_CEMC = 0;
4154  for (Int_t itow = 0; itow < _maxNTowersCentral; itow++)
4155  {
4156  _tower_CEMC_E[itow] = 0;
4157  _tower_CEMC_iEta[itow] = 0;
4158  _tower_CEMC_iPhi[itow] = 0;
4159  _tower_CEMC_trueID[itow] = 0;
4160  }
4161  if (_do_CLUSTERS)
4162  {
4163  _nclusters_CEMC = 0;
4164  for (Int_t itow = 0; itow < _maxNclustersCentral; itow++)
4165  {
4166  _cluster_CEMC_E[itow] = 0;
4167  _cluster_CEMC_Eta[itow] = 0;
4168  _cluster_CEMC_Phi[itow] = 0;
4169  _cluster_CEMC_NTower[itow] = 0;
4170  _cluster_CEMC_trueID[itow] = 0;
4171  }
4172  }
4173  if (Verbosity() > 0)
4174  {
4175  cout << "\t... CEMC variables reset" << endl;
4176  }
4177  }
4178  if (_do_HCALIN)
4179  {
4180  _nTowers_HCALIN = 0;
4181  for (Int_t itow = 0; itow < _maxNTowersCentral; itow++)
4182  {
4183  _tower_HCALIN_E[itow] = 0;
4184  _tower_HCALIN_iEta[itow] = 0;
4185  _tower_HCALIN_iPhi[itow] = 0;
4186  _tower_HCALIN_trueID[itow] = 0;
4187  }
4188  if (_do_CLUSTERS)
4189  {
4190  _nclusters_HCALIN = 0;
4191  for (Int_t itow = 0; itow < _maxNclusters; itow++)
4192  {
4193  _cluster_HCALIN_E[itow] = 0;
4194  _cluster_HCALIN_Eta[itow] = 0;
4195  _cluster_HCALIN_Phi[itow] = 0;
4196  _cluster_HCALIN_NTower[itow] = 0;
4197  _cluster_HCALIN_trueID[itow] = 0;
4198  }
4199  }
4200  if (Verbosity() > 0)
4201  {
4202  cout << "\t... HCALIN variables reset" << endl;
4203  }
4204  }
4205  if (_do_HCALOUT)
4206  {
4207  if (Verbosity() > 0)
4208  {
4209  cout << "\t... resetting HCALOUT variables" << endl;
4210  }
4211  _nTowers_HCALOUT = 0;
4212  for (Int_t itow = 0; itow < _maxNTowersCentral; itow++)
4213  {
4214  _tower_HCALOUT_E[itow] = 0;
4215  _tower_HCALOUT_iEta[itow] = 0;
4216  _tower_HCALOUT_iPhi[itow] = 0;
4217  _tower_HCALOUT_trueID[itow] = 0;
4218  }
4219  if (_do_CLUSTERS)
4220  {
4221  _nclusters_HCALOUT = 0;
4222  for (Int_t itow = 0; itow < _maxNclusters; itow++)
4223  {
4224  _cluster_HCALOUT_E[itow] = 0;
4225  _cluster_HCALOUT_Eta[itow] = 0;
4226  _cluster_HCALOUT_Phi[itow] = 0;
4227  _cluster_HCALOUT_NTower[itow] = 0;
4228  _cluster_HCALOUT_trueID[itow] = 0;
4229  }
4230  }
4231  if (Verbosity() > 0)
4232  {
4233  cout << "\t... HCALOUT variables reset" << endl;
4234  }
4235  }
4236  if (_do_EEMC)
4237  {
4238  if (Verbosity() > 0)
4239  {
4240  cout << "\t... resetting EEMC variables" << endl;
4241  }
4242  _nTowers_EEMC = 0;
4243  for (Int_t itow = 0; itow < _maxNTowers; itow++)
4244  {
4245  _tower_EEMC_E[itow] = 0;
4246  _tower_EEMC_iEta[itow] = 0;
4247  _tower_EEMC_iPhi[itow] = 0;
4248  _tower_EEMC_trueID[itow] = 0;
4249  }
4250  if (_do_CLUSTERS)
4251  {
4252  _nclusters_EEMC = 0;
4253  for (Int_t itow = 0; itow < _maxNclusters; itow++)
4254  {
4255  _cluster_EEMC_E[itow] = 0;
4256  _cluster_EEMC_Eta[itow] = 0;
4257  _cluster_EEMC_Phi[itow] = 0;
4258  _cluster_EEMC_NTower[itow] = 0;
4259  _cluster_EEMC_trueID[itow] = 0;
4260  }
4261  }
4262  if (Verbosity() > 0)
4263  {
4264  cout << "\t... EEMC variables reset" << endl;
4265  }
4266  }
4267  if (_do_EEMCG)
4268  {
4269  if (Verbosity() > 0)
4270  {
4271  cout << "\t... resetting EEMCG variables" << endl;
4272  }
4273  _nTowers_EEMCG = 0;
4274  for (Int_t itow = 0; itow < _maxNTowers; itow++)
4275  {
4276  _tower_EEMCG_E[itow] = 0;
4277  _tower_EEMCG_iEta[itow] = 0;
4278  _tower_EEMCG_iPhi[itow] = 0;
4279  _tower_EEMCG_trueID[itow] = 0;
4280  }
4281  if (_do_CLUSTERS)
4282  {
4283  _nclusters_EEMCG = 0;
4284  for (Int_t itow = 0; itow < _maxNclusters; itow++)
4285  {
4286  _cluster_EEMCG_E[itow] = 0;
4287  _cluster_EEMCG_Eta[itow] = 0;
4288  _cluster_EEMCG_Phi[itow] = 0;
4289  _cluster_EEMCG_NTower[itow] = 0;
4290  _cluster_EEMCG_trueID[itow] = 0;
4291  }
4292  }
4293  if (Verbosity() > 0)
4294  {
4295  cout << "\t... EEMCG variables reset" << endl;
4296  }
4297  }
4298  if (_do_DRCALO)
4299  {
4300  if (Verbosity() > 0)
4301  {
4302  cout << "\t... resetting DRCALO variables" << endl;
4303  }
4304  _nTowers_DRCALO = 0;
4305  for (Int_t itow = 0; itow < _maxNTowersDR; itow++)
4306  {
4307  _tower_DRCALO_E[itow] = 0;
4308  _tower_DRCALO_NScint[itow] = 0;
4309  _tower_DRCALO_NCerenkov[itow] = 0;
4310  _tower_DRCALO_iEta[itow] = 0;
4311  _tower_DRCALO_iPhi[itow] = 0;
4312  _tower_DRCALO_trueID[itow] = 0;
4313  }
4314  if (Verbosity() > 0)
4315  {
4316  cout << "\t... DRCALO variables reset" << endl;
4317  }
4318  }
4319  if (_do_FOCAL)
4320  {
4321  if (Verbosity() > 0)
4322  {
4323  cout << "\t... resetting FOCAL variables" << endl;
4324  }
4325  _nTowers_FOCAL = 0;
4326  for (Int_t itow = 0; itow < _maxNTowersDR; itow++)
4327  {
4328  _tower_FOCAL_E[itow] = 0;
4329  _tower_FOCAL_NScint[itow] = 0;
4330  _tower_FOCAL_NCerenkov[itow] = 0;
4331  _tower_FOCAL_iEta[itow] = 0;
4332  _tower_FOCAL_iPhi[itow] = 0;
4333  _tower_FOCAL_trueID[itow] = 0;
4334  }
4335  if (Verbosity() > 0)
4336  {
4337  cout << "\t... FOCAL variables reset" << endl;
4338  }
4339  }
4340  if (_do_LFHCAL)
4341  {
4342  if (Verbosity() > 0)
4343  {
4344  cout << "\t... resetting LFHCAL variables" << endl;
4345  }
4346  _nTowers_LFHCAL = 0;
4347  for (Int_t itow = 0; itow < _maxNTowers; itow++)
4348  {
4349  _tower_LFHCAL_E[itow] = 0;
4350  _tower_LFHCAL_iEta[itow] = 0;
4351  _tower_LFHCAL_iPhi[itow] = 0;
4352  _tower_LFHCAL_iL[itow] = 0;
4353  _tower_LFHCAL_trueID[itow] = 0;
4354  }
4355  if (Verbosity() > 0)
4356  {
4357  cout << "\t... LFHCAL variables reset" << endl;
4358  }
4359  }
4360  if (_do_TRACKS)
4361  {
4362  if (Verbosity() > 0)
4363  {
4364  cout << "\t... resetting Track variables" << endl;
4365  }
4366  _nTracks = 0;
4367  for (Int_t itrk = 0; itrk < _maxNTracks; itrk++)
4368  {
4369  _track_ID[itrk] = 0;
4370  _track_charge[itrk] = 0;
4371  _track_trueID[itrk] = 0;
4372  _track_px[itrk] = 0;
4373  _track_py[itrk] = 0;
4374  _track_pz[itrk] = 0;
4375  _track_x[itrk] = 0;
4376  _track_y[itrk] = 0;
4377  _track_z[itrk] = 0;
4378  _track_ndf[itrk] = 0;
4379  _track_chi2[itrk] = 0;
4380  _track_dca[itrk] = 0;
4381  _track_dca_2d[itrk] = 0;
4382  _track_source[itrk] = 0;
4383 
4384  // Default to not performing PID matching. Set default LL to same value for all PIDs
4385  _track_pion_LL[itrk] = -100;
4386  _track_kaon_LL[itrk] = -100;
4387  _track_proton_LL[itrk] = -100;
4388  }
4389  if (_do_PROJECTIONS)
4390  {
4391  _nProjections = 0;
4392  for (Int_t iproj = 0; iproj < _maxNProjections; iproj++)
4393  {
4394  _track_ProjLayer[iproj] = -1;
4395  _track_ProjTrackID[iproj] = 0;
4396  _track_TLP_x[iproj] = 0;
4397  _track_TLP_y[iproj] = 0;
4398  _track_TLP_z[iproj] = 0;
4399  _track_TLP_t[iproj] = 0;
4400  _track_TLP_px[iproj] = 0;
4401  _track_TLP_py[iproj] = 0;
4402  _track_TLP_pz[iproj] = 0;
4403  _track_TLP_true_x[iproj] = 0;
4404  _track_TLP_true_y[iproj] = 0;
4405  _track_TLP_true_z[iproj] = 0;
4406  _track_TLP_true_t[iproj] = 0;
4407  }
4408  }
4409  if (Verbosity() > 0)
4410  {
4411  cout << "\t... track variables reset" << endl;
4412  }
4413  }
4414  if (_do_MCPARTICLES)
4415  {
4416  _nMCPart = 0;
4417  for (Int_t imcpart = 0; imcpart < _maxNMCPart; imcpart++)
4418  {
4419  _mcpart_ID[imcpart] = 0;
4420  _mcpart_ID_parent[imcpart] = 0;
4421  _mcpart_PDG[imcpart] = 0;
4422  _mcpart_E[imcpart] = 0;
4423  _mcpart_px[imcpart] = 0;
4424  _mcpart_py[imcpart] = 0;
4425  _mcpart_pz[imcpart] = 0;
4426  _mcpart_x[imcpart] = 0;
4427  _mcpart_y[imcpart] = 0;
4428  _mcpart_z[imcpart] = 0;
4429  _mcpart_BCID[imcpart] = -10;
4430  }
4431  }
4432 
4433  if (_do_HEPMC)
4434  {
4435  _nHepmcp = 0;
4436  _hepmcp_vtx_x = 0;
4437  _hepmcp_vtx_y = 0;
4438  _hepmcp_vtx_z = 0;
4439  _hepmcp_vtx_t = 0;
4440  for (Int_t iHepmcp = 0; iHepmcp < _maxNHepmcp; iHepmcp++)
4441  {
4442  _hepmcp_E[iHepmcp] = 0;
4443  _hepmcp_PDG[iHepmcp] = 0;
4444  _hepmcp_px[iHepmcp] = 0;
4445  _hepmcp_py[iHepmcp] = 0;
4446  _hepmcp_pz[iHepmcp] = 0;
4447  _hepmcp_status[iHepmcp] = -10;
4448  _hepmcp_BCID[iHepmcp] = 0;
4449  _hepmcp_m2[iHepmcp] = 0;
4450  _hepmcp_m1[iHepmcp] = 0;
4451  }
4452  if (Verbosity() > 0)
4453  {
4454  cout << "\t... MC variables reset" << endl;
4455  }
4456  }
4457 }