Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:21:44 

0001 #include "EventEvaluator.h"
0002 
0003 #include "CaloEvalStack.h"
0004 #include "CaloRawClusterEval.h"
0005 #include "CaloRawTowerEval.h"
0006 #include "CaloTruthEval.h"
0007 
0008 #include <g4main/PHG4Hit.h>
0009 #include <g4main/PHG4HitContainer.h>
0010 #include <g4main/PHG4Particle.h>
0011 #include <g4main/PHG4Shower.h>
0012 #include <g4main/PHG4TruthInfoContainer.h>
0013 #include <g4main/PHG4VtxPoint.h>
0014 
0015 #include <globalvertex/GlobalVertex.h>
0016 #include <globalvertex/GlobalVertexMap.h>
0017 
0018 #include <globalvertex/SvtxVertex.h>  // for SvtxVertex
0019 #include <globalvertex/SvtxVertexMap.h>
0020 
0021 #include <trackbase_historic/SvtxTrack.h>
0022 #include <trackbase_historic/SvtxTrackMap.h>
0023 #include <trackbase_historic/SvtxTrack_FastSim.h>
0024 
0025 #include <calobase/RawCluster.h>
0026 #include <calobase/RawClusterContainer.h>
0027 #include <calobase/RawClusterUtility.h>
0028 #include <calobase/RawTower.h>
0029 #include <calobase/RawTowerContainer.h>
0030 #include <calobase/RawTowerGeom.h>
0031 #include <calobase/RawTowerGeomContainer.h>
0032 #include <calobase/RawTowerv2.h>
0033 
0034 #include <phhepmc/PHGenIntegral.h>
0035 #include <phhepmc/PHHepMCGenEvent.h>
0036 #include <phhepmc/PHHepMCGenEventMap.h>
0037 
0038 #include <fun4all/Fun4AllReturnCodes.h>
0039 #include <fun4all/SubsysReco.h>
0040 
0041 #include <phool/PHCompositeNode.h>
0042 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0043 #include <phool/getClass.h>
0044 #include <phool/phool.h>
0045 
0046 #include <TFile.h>
0047 #include <TNtuple.h>
0048 #include <TTree.h>
0049 
0050 #include <CLHEP/Vector/ThreeVector.h>
0051 
0052 #include <HepMC/GenEvent.h>
0053 #include <HepMC/GenVertex.h>
0054 
0055 #include <cmath>
0056 #include <cstdlib>
0057 #include <iostream>
0058 #include <set>
0059 #include <utility>
0060 
0061 EventEvaluator::EventEvaluator(const std::string& name, const std::string& filename)
0062   : SubsysReco(name)
0063   , _geometry_done(new int[20])
0064   , _filename(filename)
0065 {
0066   _hits_layerID = new int[_maxNHits];
0067   _hits_trueID = new int[_maxNHits];
0068   _hits_x = new float[_maxNHits];
0069   _hits_y = new float[_maxNHits];
0070   _hits_z = new float[_maxNHits];
0071   _hits_t = new float[_maxNHits];
0072 
0073   _tower_HCALIN_E = new float[_maxNTowersCentral];
0074   _tower_HCALIN_iEta = new int[_maxNTowersCentral];
0075   _tower_HCALIN_iPhi = new int[_maxNTowersCentral];
0076   _tower_HCALIN_trueID = new int[_maxNTowersCentral];
0077   _cluster_HCALIN_E = new float[_maxNclustersCentral];
0078   _cluster_HCALIN_Eta = new float[_maxNclustersCentral];
0079   _cluster_HCALIN_Phi = new float[_maxNclustersCentral];
0080   _cluster_HCALIN_NTower = new int[_maxNclustersCentral];
0081   _cluster_HCALIN_trueID = new int[_maxNclustersCentral];
0082 
0083   _tower_HCALOUT_E = new float[_maxNTowersCentral];
0084   _tower_HCALOUT_iEta = new int[_maxNTowersCentral];
0085   _tower_HCALOUT_iPhi = new int[_maxNTowersCentral];
0086   _tower_HCALOUT_trueID = new int[_maxNTowersCentral];
0087   _cluster_HCALOUT_E = new float[_maxNclustersCentral];
0088   _cluster_HCALOUT_Eta = new float[_maxNclustersCentral];
0089   _cluster_HCALOUT_Phi = new float[_maxNclustersCentral];
0090   _cluster_HCALOUT_NTower = new int[_maxNclustersCentral];
0091   _cluster_HCALOUT_trueID = new int[_maxNclustersCentral];
0092 
0093   _tower_CEMC_E = new float[_maxNTowersCentral];
0094   _tower_CEMC_iEta = new int[_maxNTowersCentral];
0095   _tower_CEMC_iPhi = new int[_maxNTowersCentral];
0096   _tower_CEMC_trueID = new int[_maxNTowersCentral];
0097   _cluster_CEMC_E = new float[_maxNclustersCentral];
0098   _cluster_CEMC_Eta = new float[_maxNclustersCentral];
0099   _cluster_CEMC_Phi = new float[_maxNclustersCentral];
0100   _cluster_CEMC_NTower = new int[_maxNclustersCentral];
0101   _cluster_CEMC_trueID = new int[_maxNclustersCentral];
0102 
0103   _track_ID = new float[_maxNTracks];
0104   _track_trueID = new float[_maxNTracks];
0105   _track_px = new float[_maxNTracks];
0106   _track_py = new float[_maxNTracks];
0107   _track_pz = new float[_maxNTracks];
0108   _track_dca = new float[_maxNTracks];
0109   _track_dca_2d = new float[_maxNTracks];
0110   _track_source = new unsigned short[_maxNTracks];
0111   _track_ProjTrackID = new float[_maxNProjections];
0112   _track_ProjLayer = new int[_maxNProjections];
0113   _track_TLP_x = new float[_maxNProjections];
0114   _track_TLP_y = new float[_maxNProjections];
0115   _track_TLP_z = new float[_maxNProjections];
0116   _track_TLP_t = new float[_maxNProjections];
0117   _track_TLP_true_x = new float[_maxNProjections];
0118   _track_TLP_true_y = new float[_maxNProjections];
0119   _track_TLP_true_z = new float[_maxNProjections];
0120   _track_TLP_true_t = new float[_maxNProjections];
0121 
0122   _mcpart_ID = new int[_maxNMCPart];
0123   _mcpart_ID_parent = new int[_maxNMCPart];
0124   _mcpart_PDG = new int[_maxNMCPart];
0125   _mcpart_E = new float[_maxNMCPart];
0126   _mcpart_px = new float[_maxNMCPart];
0127   _mcpart_py = new float[_maxNMCPart];
0128   _mcpart_pz = new float[_maxNMCPart];
0129   _mcpart_BCID = new int[_maxNMCPart];
0130 
0131   _hepmcp_BCID = new int[_maxNHepmcp];
0132   //  _hepmcp_ID_parent = new float[_maxNHepmcp];
0133   _hepmcp_status = new int[_maxNHepmcp];
0134   _hepmcp_PDG = new int[_maxNHepmcp];
0135   _hepmcp_E = new float[_maxNHepmcp];
0136   _hepmcp_px = new float[_maxNHepmcp];
0137   _hepmcp_py = new float[_maxNHepmcp];
0138   _hepmcp_pz = new float[_maxNHepmcp];
0139   _hepmcp_m1 = new int[_maxNHepmcp];
0140   _hepmcp_m2 = new int[_maxNHepmcp];
0141 
0142   _calo_towers_iEta = new int[_maxNTowersCalo];
0143   _calo_towers_iPhi = new int[_maxNTowersCalo];
0144   _calo_towers_Eta = new float[_maxNTowersCalo];
0145   _calo_towers_Phi = new float[_maxNTowersCalo];
0146   _calo_towers_x = new float[_maxNTowersCalo];
0147   _calo_towers_y = new float[_maxNTowersCalo];
0148   _calo_towers_z = new float[_maxNTowersCalo];
0149 
0150   for (int igem = 0; igem < 20; igem++)
0151   {
0152     _geometry_done[igem] = 0;
0153   }
0154 }
0155 
0156 int EventEvaluator::Init(PHCompositeNode* /*topNode*/)
0157 {
0158   _ievent = 0;
0159 
0160   _tfile = new TFile(_filename.c_str(), "RECREATE");
0161 
0162   _event_tree = new TTree("event_tree", "event_tree");
0163   if (_do_store_event_info)
0164   {
0165     // Event level info. This isn't the most efficient way to store this info, but it's straightforward
0166     // within the structure of the class, so the size is small compared to the rest of the output.
0167     _event_tree->Branch("cross_section", &_cross_section, "cross_section/F");
0168     _event_tree->Branch("event_weight", &_event_weight, "event_weight/F");
0169     _event_tree->Branch("n_generator_accepted", &_n_generator_accepted, "n_generator_accepted/I");
0170   }
0171   // tracks and hits
0172   if (_do_HITS)
0173   {
0174     _event_tree->Branch("nHits", &_nHitsLayers, "nHits/I");
0175     _event_tree->Branch("hits_layerID", _hits_layerID, "hits_layerID[nHits]/I");
0176     _event_tree->Branch("hits_trueID", _hits_trueID, "hits_trueID[nHits]/I");
0177     _event_tree->Branch("hits_x", _hits_x, "hits_x[nHits]/F");
0178     _event_tree->Branch("hits_y", _hits_y, "hits_y[nHits]/F");
0179     _event_tree->Branch("hits_z", _hits_z, "hits_z[nHits]/F");
0180     _event_tree->Branch("hits_t", _hits_t, "hits_t[nHits]/F");
0181   }
0182   if (_do_TRACKS)
0183   {
0184     _event_tree->Branch("nTracks", &_nTracks, "nTracks/I");
0185     _event_tree->Branch("tracks_ID", _track_ID, "tracks_ID[nTracks]/F");
0186     _event_tree->Branch("tracks_px", _track_px, "tracks_px[nTracks]/F");
0187     _event_tree->Branch("tracks_py", _track_py, "tracks_py[nTracks]/F");
0188     _event_tree->Branch("tracks_pz", _track_pz, "tracks_pz[nTracks]/F");
0189     _event_tree->Branch("tracks_dca", _track_dca, "tracks_dca[nTracks]/F");
0190     _event_tree->Branch("tracks_dca_2d", _track_dca_2d, "tracks_dca_2d[nTracks]/F");
0191     _event_tree->Branch("tracks_trueID", _track_trueID, "tracks_trueID[nTracks]/F");
0192     _event_tree->Branch("tracks_source", _track_source, "tracks_source[nTracks]/s");
0193   }
0194   if (_do_PROJECTIONS)
0195   {
0196     _event_tree->Branch("nProjections", &_nProjections, "nProjections/I");
0197     _event_tree->Branch("track_ProjTrackID", _track_ProjTrackID, "track_ProjTrackID[nProjections]/F");
0198     _event_tree->Branch("track_ProjLayer", _track_ProjLayer, "track_ProjLayer[nProjections]/I");
0199     _event_tree->Branch("track_TLP_x", _track_TLP_x, "track_TLP_x[nProjections]/F");
0200     _event_tree->Branch("track_TLP_y", _track_TLP_y, "track_TLP_y[nProjections]/F");
0201     _event_tree->Branch("track_TLP_z", _track_TLP_z, "track_TLP_z[nProjections]/F");
0202     _event_tree->Branch("track_TLP_t", _track_TLP_t, "track_TLP_t[nProjections]/F");
0203     _event_tree->Branch("track_TLP_true_x", _track_TLP_true_x, "track_TLP_true_x[nProjections]/F");
0204     _event_tree->Branch("track_TLP_true_y", _track_TLP_true_y, "track_TLP_true_y[nProjections]/F");
0205     _event_tree->Branch("track_TLP_true_z", _track_TLP_true_z, "track_TLP_true_z[nProjections]/F");
0206     _event_tree->Branch("track_TLP_true_t", _track_TLP_true_t, "track_TLP_true_t[nProjections]/F");
0207   }
0208   if (_do_HCALIN)
0209   {
0210     // towers HCAL-in
0211     _event_tree->Branch("tower_HCALIN_N", &_nTowers_HCALIN, "tower_HCALIN_N/I");
0212     _event_tree->Branch("tower_HCALIN_E", _tower_HCALIN_E, "tower_HCALIN_E[tower_HCALIN_N]/F");
0213     _event_tree->Branch("tower_HCALIN_iEta", _tower_HCALIN_iEta, "tower_HCALIN_iEta[tower_HCALIN_N]/I");
0214     _event_tree->Branch("tower_HCALIN_iPhi", _tower_HCALIN_iPhi, "tower_HCALIN_iPhi[tower_HCALIN_N]/I");
0215     _event_tree->Branch("tower_HCALIN_trueID", _tower_HCALIN_trueID, "tower_HCALIN_trueID[tower_HCALIN_N]/I");
0216     if (_do_CLUSTERS)
0217     {
0218       // clusters HCAL-in
0219       _event_tree->Branch("cluster_HCALIN_N", &_nclusters_HCALIN, "cluster_HCALIN_N/I");
0220       _event_tree->Branch("cluster_HCALIN_E", _cluster_HCALIN_E, "cluster_HCALIN_E[cluster_HCALIN_N]/F");
0221       _event_tree->Branch("cluster_HCALIN_Eta", _cluster_HCALIN_Eta, "cluster_HCALIN_Eta[cluster_HCALIN_N]/F");
0222       _event_tree->Branch("cluster_HCALIN_Phi", _cluster_HCALIN_Phi, "cluster_HCALIN_Phi[cluster_HCALIN_N]/F");
0223       _event_tree->Branch("cluster_HCALIN_NTower", _cluster_HCALIN_NTower, "cluster_HCALIN_NTower[cluster_HCALIN_N]/I");
0224       _event_tree->Branch("cluster_HCALIN_trueID", _cluster_HCALIN_trueID, "cluster_HCALIN_trueID[cluster_HCALIN_N]/I");
0225     }
0226   }
0227   if (_do_HCALOUT)
0228   {
0229     // towers HCAL-out
0230     _event_tree->Branch("tower_HCALOUT_N", &_nTowers_HCALOUT, "tower_HCALOUT_N/I");
0231     _event_tree->Branch("tower_HCALOUT_E", _tower_HCALOUT_E, "tower_HCALOUT_E[tower_HCALOUT_N]/F");
0232     _event_tree->Branch("tower_HCALOUT_iEta", _tower_HCALOUT_iEta, "tower_HCALOUT_iEta[tower_HCALOUT_N]/I");
0233     _event_tree->Branch("tower_HCALOUT_iPhi", _tower_HCALOUT_iPhi, "tower_HCALOUT_iPhi[tower_HCALOUT_N]/I");
0234     _event_tree->Branch("tower_HCALOUT_trueID", _tower_HCALOUT_trueID, "tower_HCALOUT_trueID[tower_HCALOUT_N]/I");
0235     if (_do_CLUSTERS)
0236     {
0237       // clusters HCAL-out
0238       _event_tree->Branch("cluster_HCALOUT_N", &_nclusters_HCALOUT, "cluster_HCALOUT_N/I");
0239       _event_tree->Branch("cluster_HCALOUT_E", _cluster_HCALOUT_E, "cluster_HCALOUT_E[cluster_HCALOUT_N]/F");
0240       _event_tree->Branch("cluster_HCALOUT_Eta", _cluster_HCALOUT_Eta, "cluster_HCALOUT_Eta[cluster_HCALOUT_N]/F");
0241       _event_tree->Branch("cluster_HCALOUT_Phi", _cluster_HCALOUT_Phi, "cluster_HCALOUT_Phi[cluster_HCALOUT_N]/F");
0242       _event_tree->Branch("cluster_HCALOUT_NTower", _cluster_HCALOUT_NTower, "cluster_HCALOUT_NTower[cluster_HCALOUT_N]/I");
0243       _event_tree->Branch("cluster_HCALOUT_trueID", _cluster_HCALOUT_trueID, "cluster_HCALOUT_trueID[cluster_HCALOUT_N]/I");
0244     }
0245   }
0246   if (_do_CEMC)
0247   {
0248     // towers CEMC
0249     _event_tree->Branch("tower_CEMC_N", &_nTowers_CEMC, "tower_CEMC_N/I");
0250     _event_tree->Branch("tower_CEMC_E", _tower_CEMC_E, "tower_CEMC_E[tower_CEMC_N]/F");
0251     _event_tree->Branch("tower_CEMC_iEta", _tower_CEMC_iEta, "tower_CEMC_iEta[tower_CEMC_N]/I");
0252     _event_tree->Branch("tower_CEMC_iPhi", _tower_CEMC_iPhi, "tower_CEMC_iPhi[tower_CEMC_N]/I");
0253     _event_tree->Branch("tower_CEMC_trueID", _tower_CEMC_trueID, "tower_CEMC_trueID[tower_CEMC_N]/I");
0254     if (_do_CLUSTERS)
0255     {
0256       // clusters CEMC
0257       _event_tree->Branch("cluster_CEMC_N", &_nclusters_CEMC, "cluster_CEMC_N/I");
0258       _event_tree->Branch("cluster_CEMC_E", _cluster_CEMC_E, "cluster_CEMC_E[cluster_CEMC_N]/F");
0259       _event_tree->Branch("cluster_CEMC_Eta", _cluster_CEMC_Eta, "cluster_CEMC_Eta[cluster_CEMC_N]/F");
0260       _event_tree->Branch("cluster_CEMC_Phi", _cluster_CEMC_Phi, "cluster_CEMC_Phi[cluster_CEMC_N]/F");
0261       _event_tree->Branch("cluster_CEMC_NTower", _cluster_CEMC_NTower, "cluster_CEMC_NTower[cluster_CEMC_N]/I");
0262       _event_tree->Branch("cluster_CEMC_trueID", _cluster_CEMC_trueID, "cluster_CEMC_trueID[cluster_CEMC_N]/I");
0263     }
0264   }
0265   if (_do_VERTEX)
0266   {
0267     // vertex
0268     _event_tree->Branch("vertex_x", &_vertex_x, "vertex_x/F");
0269     _event_tree->Branch("vertex_y", &_vertex_y, "vertex_y/F");
0270     _event_tree->Branch("vertex_z", &_vertex_z, "vertex_z/F");
0271     _event_tree->Branch("vertex_NCont", &_vertex_NCont, "vertex_NCont/I");
0272     _event_tree->Branch("vertex_true_x", &_vertex_true_x, "vertex_true_x/F");
0273     _event_tree->Branch("vertex_true_y", &_vertex_true_y, "vertex_true_y/F");
0274     _event_tree->Branch("vertex_true_z", &_vertex_true_z, "vertex_true_z/F");
0275   }
0276   if (_do_MCPARTICLES)
0277   {
0278     // MC particles
0279     _event_tree->Branch("nMCPart", &_nMCPart, "nMCPart/I");
0280     _event_tree->Branch("mcpart_ID", _mcpart_ID, "mcpart_ID[nMCPart]/I");
0281     _event_tree->Branch("mcpart_ID_parent", _mcpart_ID_parent, "mcpart_ID_parent[nMCPart]/I");
0282     _event_tree->Branch("mcpart_PDG", _mcpart_PDG, "mcpart_PDG[nMCPart]/I");
0283     _event_tree->Branch("mcpart_E", _mcpart_E, "mcpart_E[nMCPart]/F");
0284     _event_tree->Branch("mcpart_px", _mcpart_px, "mcpart_px[nMCPart]/F");
0285     _event_tree->Branch("mcpart_py", _mcpart_py, "mcpart_py[nMCPart]/F");
0286     _event_tree->Branch("mcpart_pz", _mcpart_pz, "mcpart_pz[nMCPart]/F");
0287     _event_tree->Branch("mcpart_BCID", _mcpart_BCID, "mcpart_BCID[nMCPart]/I");
0288   }
0289   if (_do_HEPMC)
0290   {
0291     // MC particles
0292     _event_tree->Branch("nHepmcp", &_nHepmcp, "nHepmcp/I");
0293     _event_tree->Branch("hepmcp_procid", &_hepmcp_procid, "hepmcp_procid/I");
0294     _event_tree->Branch("hepmcp_x1", &_hepmcp_x1, "hepmcp_x1/F");
0295     _event_tree->Branch("hepmcp_x2", &_hepmcp_x2, "hepmcp_x2/F");
0296 
0297     //    _event_tree->Branch("hepmcp_ID_parent", _hepmcp_ID_parent, "hepmcp_ID_parent[nHepmcp]/F");
0298     _event_tree->Branch("hepmcp_status", _hepmcp_status, "hepmcp_status[nHepmcp]/I");
0299     _event_tree->Branch("hepmcp_PDG", _hepmcp_PDG, "hepmcp_PDG[nHepmcp]/I");
0300     _event_tree->Branch("hepmcp_E", _hepmcp_E, "hepmcp_E[nHepmcp]/F");
0301     _event_tree->Branch("hepmcp_px", _hepmcp_px, "hepmcp_px[nHepmcp]/F");
0302     _event_tree->Branch("hepmcp_py", _hepmcp_py, "hepmcp_py[nHepmcp]/F");
0303     _event_tree->Branch("hepmcp_pz", _hepmcp_pz, "hepmcp_pz[nHepmcp]/F");
0304     _event_tree->Branch("hepmcp_BCID", _hepmcp_BCID, "hepmcp_BCID[nHepmcp]/I");
0305     _event_tree->Branch("hepmcp_m1", _hepmcp_m1, "hepmcp_m1[nHepmcp]/I");
0306     _event_tree->Branch("hepmcp_m2", _hepmcp_m2, "hepmcp_m2[nHepmcp]/I");
0307   }
0308 
0309   if (_do_GEOMETRY)
0310   {
0311     _tfile_geometry = new TFile("geometry.root", "RECREATE");
0312 
0313     _geometry_tree = new TTree("geometry_tree", "geometry_tree");
0314     // tracks and hits
0315     _geometry_tree->Branch("calo", &_calo_ID, "nHits/I");
0316     _geometry_tree->Branch("calo_towers_N", &_calo_towers_N, "calo_towers_N/I");
0317     _geometry_tree->Branch("calo_towers_iEta", _calo_towers_iEta, "calo_towers_iEta[calo_towers_N]/I");
0318     _geometry_tree->Branch("calo_towers_iPhi", _calo_towers_iPhi, "calo_towers_iPhi[calo_towers_N]/I");
0319     _geometry_tree->Branch("calo_towers_Eta", _calo_towers_Eta, "calo_towers_Eta[calo_towers_N]/F");
0320     _geometry_tree->Branch("calo_towers_Phi", _calo_towers_Phi, "calo_towers_Phi[calo_towers_N]/F");
0321     _geometry_tree->Branch("calo_towers_x", _calo_towers_x, "calo_towers_x[calo_towers_N]/F");
0322     _geometry_tree->Branch("calo_towers_y", _calo_towers_y, "calo_towers_y[calo_towers_N]/F");
0323     _geometry_tree->Branch("calo_towers_z", _calo_towers_z, "calo_towers_z[calo_towers_N]/F");
0324   }
0325 
0326   return Fun4AllReturnCodes::EVENT_OK;
0327 }
0328 
0329 int EventEvaluator::process_event(PHCompositeNode* topNode)
0330 {
0331   if (Verbosity() > 0)
0332   {
0333     std::cout << "entered process_event" << std::endl;
0334   }
0335   if (_do_HCALIN)
0336   {
0337     if (!_caloevalstackHCALIN)
0338     {
0339       _caloevalstackHCALIN = new CaloEvalStack(topNode, "HCALIN");
0340       _caloevalstackHCALIN->set_strict(_strict);
0341       _caloevalstackHCALIN->set_verbosity(Verbosity() + 1);
0342     }
0343     else
0344     {
0345       _caloevalstackHCALIN->next_event(topNode);
0346     }
0347   }
0348   if (_do_HCALOUT)
0349   {
0350     if (!_caloevalstackHCALOUT)
0351     {
0352       _caloevalstackHCALOUT = new CaloEvalStack(topNode, "HCALOUT");
0353       _caloevalstackHCALOUT->set_strict(_strict);
0354       _caloevalstackHCALOUT->set_verbosity(Verbosity() + 1);
0355     }
0356     else
0357     {
0358       _caloevalstackHCALOUT->next_event(topNode);
0359     }
0360   }
0361   if (_do_CEMC)
0362   {
0363     if (!_caloevalstackCEMC)
0364     {
0365       _caloevalstackCEMC = new CaloEvalStack(topNode, "CEMC");
0366       _caloevalstackCEMC->set_strict(_strict);
0367       _caloevalstackCEMC->set_verbosity(Verbosity() + 1);
0368     }
0369     else
0370     {
0371       _caloevalstackCEMC->next_event(topNode);
0372     }
0373   }
0374 
0375   if (Verbosity() > 0)
0376   {
0377     std::cout << "loaded evalstack" << std::endl;
0378   }
0379 
0380   // fill the Evaluator Tree
0381   fillOutputNtuples(topNode);
0382 
0383   ++_ievent;
0384 
0385   return Fun4AllReturnCodes::EVENT_OK;
0386 }
0387 
0388 void EventEvaluator::fillOutputNtuples(PHCompositeNode* topNode)
0389 {
0390   if (Verbosity() > 2)
0391   {
0392     std::cout << "EventEvaluator::fillOutputNtuples() entered" << std::endl;
0393   }
0394 
0395   //----------------------
0396   // fill the Event Tree
0397   //----------------------
0398 
0399   //----------------------
0400   // Event level info
0401   //---------------------
0402   // Extract weight info from the stored HepMC event.
0403   if (_do_store_event_info)
0404   {
0405     PHHepMCGenEventMap* hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0406     if (hepmceventmap)
0407     {
0408       if (Verbosity() > 0)
0409       {
0410         std::cout << "saving event level info" << std::endl;
0411       }
0412 
0413       for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
0414            eventIter != hepmceventmap->end();
0415            ++eventIter)
0416       {
0417         PHHepMCGenEvent* hepmcevent = eventIter->second;
0418 
0419         if (hepmcevent)
0420         {
0421           HepMC::GenEvent* truthevent = hepmcevent->getEvent();
0422           if (!truthevent)
0423           {
0424             std::cout << PHWHERE
0425                       << "no evt pointer under phhepmvgeneventmap found "
0426                       << std::endl;
0427             return;
0428           }
0429 
0430           auto* xsec = truthevent->cross_section();
0431           if (xsec)
0432           {
0433             _cross_section = xsec->cross_section();
0434           }
0435           // Only fill the event weight if available.
0436           // The overall event weight will be stored in the last entry in the vector.
0437           auto weights = truthevent->weights();
0438           if (!weights.empty())
0439           {
0440             _event_weight = weights[weights.size() - 1];
0441           }
0442         }
0443       }
0444     }
0445     else
0446     {
0447       if (Verbosity() > 0)
0448       {
0449         std::cout << PHWHERE << " PHHepMCGenEventMap node (for event level info) not found on node tree" << std::endl;
0450       }
0451       return;
0452     }
0453 
0454     // Retrieve the number of generator accepted events
0455     // Following how this was implemented in PHPythia8
0456     PHNodeIterator iter(topNode);
0457     PHCompositeNode* sumNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
0458     if (!sumNode)
0459     {
0460       std::cout << PHWHERE << "RUN Node missing doing nothing" << std::endl;
0461       return;
0462     }
0463     auto* integralNode = findNode::getClass<PHGenIntegral>(sumNode, "PHGenIntegral");
0464     if (integralNode)
0465     {
0466       _n_generator_accepted = integralNode->get_N_Generator_Accepted_Event();
0467     }
0468     else
0469     {
0470       if (Verbosity() > 0)
0471       {
0472         std::cout << PHWHERE << " PHGenIntegral node (for n generator accepted) not found on node tree. Continuing" << std::endl;
0473       }
0474     }
0475   }
0476   //----------------------
0477   //    VERTEX
0478   //----------------------
0479   SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0480   if (_do_VERTEX)
0481   {
0482     if (vertexmap)
0483     {
0484       if (!vertexmap->empty())
0485       {
0486         if (Verbosity() > 0)
0487         {
0488           std::cout << "saving vertex" << std::endl;
0489         }
0490         SvtxVertex* vertex = (vertexmap->begin())->second;
0491 
0492         _vertex_x = vertex->get_x();
0493         _vertex_y = vertex->get_y();
0494         _vertex_z = vertex->get_z();
0495         _vertex_NCont = vertex->size_tracks();
0496       }
0497       else
0498       {
0499         _vertex_x = 0.;
0500         _vertex_y = 0.;
0501         _vertex_z = 0.;
0502         _vertex_NCont = -1;
0503       }
0504     }
0505   }
0506   //----------------------
0507   //    HITS
0508   //----------------------
0509   if (_do_HITS)
0510   {
0511     if (Verbosity() > 0)
0512     {
0513       std::cout << "saving hits" << std::endl;
0514     }
0515     _nHitsLayers = 0;
0516     PHG4TruthInfoContainer* truthinfocontainerHits = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0517     for (int iIndex = 0; iIndex < 60; ++iIndex)
0518     {
0519       // you need to add your layer name here to be saved! This has to be done
0520       // as we do not want to save thousands of calorimeter hits!
0521       if ((GetProjectionNameFromIndex(iIndex).find("MVTX") != std::string::npos) ||
0522           (GetProjectionNameFromIndex(iIndex).find("INTT") != std::string::npos))
0523       {
0524         std::string nodename = "G4HIT_" + GetProjectionNameFromIndex(iIndex);
0525         PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
0526         if (hits)
0527         {
0528           if (Verbosity() > 1)
0529           {
0530             std::cout << __PRETTY_FUNCTION__ << " number of hits: " << hits->size() << std::endl;
0531           }
0532           PHG4HitContainer::ConstRange hit_range = hits->getHits();
0533           for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0534           {
0535             if (Verbosity() > 1)
0536             {
0537               std::cout << __PRETTY_FUNCTION__ << " found hit with id " << hit_iter->second->get_trkid() << std::endl;
0538             }
0539             if (_nHitsLayers > _maxNHits)
0540             {
0541               std::cout << __PRETTY_FUNCTION__ << " exceededed maximum hit array size! Please check where these hits come from!" << std::endl;
0542               break;
0543             }
0544             _hits_x[_nHitsLayers] = hit_iter->second->get_x(0);
0545             _hits_y[_nHitsLayers] = hit_iter->second->get_y(0);
0546             _hits_z[_nHitsLayers] = hit_iter->second->get_z(0);
0547             _hits_t[_nHitsLayers] = hit_iter->second->get_t(0);
0548             _hits_layerID[_nHitsLayers] = iIndex;
0549             if (truthinfocontainerHits)
0550             {
0551               PHG4Particle* particle = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
0552 
0553               if (particle->get_parent_id() != 0)
0554               {
0555                 PHG4Particle* g4particleMother = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
0556                 int mcSteps = 0;
0557                 while (g4particleMother->get_parent_id() != 0)
0558                 {
0559                   g4particleMother = truthinfocontainerHits->GetParticle(g4particleMother->get_parent_id());
0560                   if (g4particleMother == nullptr)
0561                   {
0562                     break;
0563                   }
0564                   mcSteps += 1;
0565                 }
0566                 if (mcSteps <= _depth_MCstack)
0567                 {
0568                   _hits_trueID[_nHitsLayers] = hit_iter->second->get_trkid();
0569                 }
0570                 else
0571                 {
0572                   PHG4Particle* g4particleMother2 = truthinfocontainerHits->GetParticle(hit_iter->second->get_trkid());
0573                   int mcSteps2 = 0;
0574                   while (g4particleMother2->get_parent_id() != 0 && (mcSteps2 < (mcSteps - _depth_MCstack + 1)))
0575                   {
0576                     g4particleMother2 = truthinfocontainerHits->GetParticle(g4particleMother2->get_parent_id());
0577                     if (g4particleMother2 == nullptr)
0578                     {
0579                       break;
0580                     }
0581 
0582                     _hits_trueID[_nHitsLayers] = g4particleMother2->get_parent_id();
0583                     mcSteps2 += 1;
0584                   }
0585                 }
0586               }
0587               else
0588               {
0589                 _hits_trueID[_nHitsLayers] = hit_iter->second->get_trkid();
0590               }
0591             }
0592             _nHitsLayers++;
0593           }
0594           if (Verbosity() > 0)
0595           {
0596             std::cout << "saved\t" << _nHitsLayers << "\thits for " << GetProjectionNameFromIndex(iIndex) << std::endl;
0597           }
0598         }
0599         else
0600         {
0601           if (Verbosity() > 0)
0602           {
0603             std::cout << __PRETTY_FUNCTION__ << " could not find " << nodename << std::endl;
0604           }
0605           continue;
0606         }
0607       }
0608     }
0609   }
0610   //----------------------
0611   //    TOWERS HCALIN
0612   //----------------------
0613   if (_do_HCALIN)
0614   {
0615     CaloRawTowerEval* towerevalHCALIN = _caloevalstackHCALIN->get_rawtower_eval();
0616     _nTowers_HCALIN = 0;
0617     std::string towernodeHCALIN = "TOWER_CALIB_HCALIN";
0618     RawTowerContainer* towersHCALIN = findNode::getClass<RawTowerContainer>(topNode, towernodeHCALIN);
0619     if (towersHCALIN)
0620     {
0621       if (Verbosity() > 0)
0622       {
0623         std::cout << "saving HCAL towers" << std::endl;
0624       }
0625       std::string towergeomnodeHCALIN = "TOWERGEOM_HCALIN";
0626       RawTowerGeomContainer* towergeomHCALIN = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeHCALIN);
0627       if (towergeomHCALIN)
0628       {
0629         if (_do_GEOMETRY && !_geometry_done[kHCALIN])
0630         {
0631           RawTowerGeomContainer::ConstRange all_towers = towergeomHCALIN->get_tower_geometries();
0632           for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
0633                it != all_towers.second; ++it)
0634           {
0635             _calo_ID = kHCALIN;
0636             _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
0637             _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
0638             _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
0639             _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
0640             _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
0641             _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
0642             _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
0643             _calo_towers_N++;
0644           }
0645           _geometry_done[kHCALIN] = 1;
0646           _geometry_tree->Fill();
0647           resetGeometryArrays();
0648         }
0649         RawTowerContainer::ConstRange begin_end = towersHCALIN->getTowers();
0650         RawTowerContainer::ConstIterator rtiter;
0651         for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0652         {
0653           RawTower* tower = rtiter->second;
0654           if (tower)
0655           {
0656             // min energy cut
0657             if (tower->get_energy() < _reco_e_threshold)
0658             {
0659               continue;
0660             }
0661             _tower_HCALIN_iEta[_nTowers_HCALIN] = tower->get_bineta();
0662             _tower_HCALIN_iPhi[_nTowers_HCALIN] = tower->get_binphi();
0663             _tower_HCALIN_E[_nTowers_HCALIN] = tower->get_energy();
0664 
0665             PHG4Particle* primary = towerevalHCALIN->max_truth_primary_particle_by_energy(tower);
0666             if (primary)
0667             {
0668               _tower_HCALIN_trueID[_nTowers_HCALIN] = primary->get_track_id();
0669               // gflavor = primary->get_pid();
0670               // efromtruth = towerevalHCALIN->get_energy_contribution(tower, primary);
0671             }
0672             else
0673             {
0674               _tower_HCALIN_trueID[_nTowers_HCALIN] = -10;
0675             }
0676             _nTowers_HCALIN++;
0677           }
0678         }
0679       }
0680       else
0681       {
0682         if (Verbosity() > 0)
0683         {
0684           std::cout << PHWHERE << " ERROR: Can't find " << towergeomnodeHCALIN << std::endl;
0685         }
0686         // return;
0687       }
0688       if (Verbosity() > 0)
0689       {
0690         std::cout << "saved\t" << _nTowers_HCALIN << "\tHCALIN towers" << std::endl;
0691       }
0692     }
0693     else
0694     {
0695       if (Verbosity() > 0)
0696       {
0697         std::cout << PHWHERE << " ERROR: Can't find " << towernodeHCALIN << std::endl;
0698       }
0699       // return;
0700     }
0701   }
0702   //----------------------
0703   //    TOWERS HCALOUT
0704   //----------------------
0705   if (_do_HCALOUT)
0706   {
0707     CaloRawTowerEval* towerevalHCALOUT = _caloevalstackHCALOUT->get_rawtower_eval();
0708     _nTowers_HCALOUT = 0;
0709     std::string towernodeHCALOUT = "TOWER_CALIB_HCALOUT";
0710     RawTowerContainer* towersHCALOUT = findNode::getClass<RawTowerContainer>(topNode, towernodeHCALOUT);
0711     if (towersHCALOUT)
0712     {
0713       if (Verbosity() > 0)
0714       {
0715         std::cout << "saving HCAL towers" << std::endl;
0716       }
0717       std::string towergeomnodeHCALOUT = "TOWERGEOM_HCALOUT";
0718       RawTowerGeomContainer* towergeomHCALOUT = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeHCALOUT);
0719       if (towergeomHCALOUT)
0720       {
0721         if (_do_GEOMETRY && !_geometry_done[kHCALOUT])
0722         {
0723           RawTowerGeomContainer::ConstRange all_towers = towergeomHCALOUT->get_tower_geometries();
0724           for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
0725                it != all_towers.second; ++it)
0726           {
0727             _calo_ID = kHCALOUT;
0728             _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
0729             _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
0730             _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
0731             _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
0732             _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
0733             _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
0734             _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
0735             _calo_towers_N++;
0736           }
0737           _geometry_done[kHCALOUT] = 1;
0738           _geometry_tree->Fill();
0739           resetGeometryArrays();
0740         }
0741         RawTowerContainer::ConstRange begin_end = towersHCALOUT->getTowers();
0742         RawTowerContainer::ConstIterator rtiter;
0743         for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0744         {
0745           RawTower* tower = rtiter->second;
0746           if (tower)
0747           {
0748             // min energy cut
0749             if (tower->get_energy() < _reco_e_threshold)
0750             {
0751               continue;
0752             }
0753             _tower_HCALOUT_iEta[_nTowers_HCALOUT] = tower->get_bineta();
0754             _tower_HCALOUT_iPhi[_nTowers_HCALOUT] = tower->get_binphi();
0755             _tower_HCALOUT_E[_nTowers_HCALOUT] = tower->get_energy();
0756 
0757             PHG4Particle* primary = towerevalHCALOUT->max_truth_primary_particle_by_energy(tower);
0758             if (primary)
0759             {
0760               _tower_HCALOUT_trueID[_nTowers_HCALOUT] = primary->get_track_id();
0761               // gflavor = primary->get_pid();
0762               // efromtruth = towerevalHCALOUT->get_energy_contribution(tower, primary);
0763             }
0764             else
0765             {
0766               _tower_HCALOUT_trueID[_nTowers_HCALOUT] = -10;
0767             }
0768             _nTowers_HCALOUT++;
0769           }
0770         }
0771       }
0772       else
0773       {
0774         if (Verbosity() > 0)
0775         {
0776           std::cout << PHWHERE << " ERROR: Can't find " << towergeomnodeHCALOUT << std::endl;
0777         }
0778         // return;
0779       }
0780       if (Verbosity() > 0)
0781       {
0782         std::cout << "saved\t" << _nTowers_HCALOUT << "\tHCALOUT towers" << std::endl;
0783       }
0784     }
0785     else
0786     {
0787       if (Verbosity() > 0)
0788       {
0789         std::cout << PHWHERE << " ERROR: Can't find " << towernodeHCALOUT << std::endl;
0790       }
0791       // return;
0792     }
0793   }
0794   //----------------------
0795   //    TOWERS CEMC
0796   //----------------------
0797   if (_do_CEMC)
0798   {
0799     CaloRawTowerEval* towerevalCEMC = _caloevalstackCEMC->get_rawtower_eval();
0800     _nTowers_CEMC = 0;
0801     std::string towernodeCEMC = "TOWER_CALIB_CEMC";
0802     RawTowerContainer* towersCEMC = findNode::getClass<RawTowerContainer>(topNode, towernodeCEMC);
0803     if (towersCEMC)
0804     {
0805       if (Verbosity() > 0)
0806       {
0807         std::cout << "saving EMC towers" << std::endl;
0808       }
0809       std::string towergeomnodeCEMC = "TOWERGEOM_CEMC";
0810       RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodeCEMC);
0811       if (towergeom)
0812       {
0813         if (_do_GEOMETRY && !_geometry_done[kCEMC])
0814         {
0815           RawTowerGeomContainer::ConstRange all_towers = towergeom->get_tower_geometries();
0816           for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
0817                it != all_towers.second; ++it)
0818           {
0819             _calo_ID = kCEMC;
0820             _calo_towers_iEta[_calo_towers_N] = it->second->get_bineta();
0821             _calo_towers_iPhi[_calo_towers_N] = it->second->get_binphi();
0822             _calo_towers_Eta[_calo_towers_N] = it->second->get_eta();
0823             _calo_towers_Phi[_calo_towers_N] = it->second->get_phi();
0824             _calo_towers_x[_calo_towers_N] = it->second->get_center_x();
0825             _calo_towers_y[_calo_towers_N] = it->second->get_center_y();
0826             _calo_towers_z[_calo_towers_N] = it->second->get_center_z();
0827             _calo_towers_N++;
0828           }
0829           _geometry_done[kCEMC] = 1;
0830           _geometry_tree->Fill();
0831           resetGeometryArrays();
0832         }
0833         RawTowerContainer::ConstRange begin_end = towersCEMC->getTowers();
0834         RawTowerContainer::ConstIterator rtiter;
0835         for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0836         {
0837           RawTower* tower = rtiter->second;
0838           if (tower)
0839           {
0840             // min energy cut
0841             if (tower->get_energy() < _reco_e_threshold)
0842             {
0843               continue;
0844             }
0845 
0846             _tower_CEMC_iEta[_nTowers_CEMC] = tower->get_bineta();
0847             _tower_CEMC_iPhi[_nTowers_CEMC] = tower->get_binphi();
0848             _tower_CEMC_E[_nTowers_CEMC] = tower->get_energy();
0849 
0850             PHG4Particle* primary = towerevalCEMC->max_truth_primary_particle_by_energy(tower);
0851             if (primary)
0852             {
0853               _tower_CEMC_trueID[_nTowers_CEMC] = primary->get_track_id();
0854               // gflavor = primary->get_pid();
0855               // efromtruth = towerevalCEMC->get_energy_contribution(tower, primary);
0856             }
0857             else
0858             {
0859               _tower_CEMC_trueID[_nTowers_CEMC] = -10;
0860             }
0861             _nTowers_CEMC++;
0862           }
0863         }
0864       }
0865       else
0866       {
0867         if (Verbosity() > 0)
0868         {
0869           std::cout << PHWHERE << " ERROR: Can't find " << towergeomnodeCEMC << std::endl;
0870         }
0871         // return;
0872       }
0873       if (Verbosity() > 0)
0874       {
0875         std::cout << "saved\t" << _nTowers_CEMC << "\tCEMC towers" << std::endl;
0876       }
0877     }
0878     else
0879     {
0880       if (Verbosity() > 0)
0881       {
0882         std::cout << PHWHERE << " ERROR: Can't find " << towernodeCEMC << std::endl;
0883       }
0884       // return;
0885     }
0886   }
0887 
0888   //------------------------
0889   // CLUSTERS HCALIN
0890   //------------------------
0891   if (_do_HCALIN && _do_CLUSTERS)
0892   {
0893     CaloRawClusterEval* clusterevalHCALIN = _caloevalstackHCALIN->get_rawcluster_eval();
0894     _nclusters_HCALIN = 0;
0895     if (Verbosity() > 1)
0896     {
0897       std::cout << "CaloEvaluator::filling gcluster ntuple..." << std::endl;
0898     }
0899 
0900     std::string clusternodeHCALIN = "CLUSTER_HCALIN";
0901     RawClusterContainer* clustersHCALIN = findNode::getClass<RawClusterContainer>(topNode, clusternodeHCALIN);
0902     if (clustersHCALIN)
0903     {
0904       // for every cluster
0905       for (const auto& iterator : clustersHCALIN->getClustersMap())
0906       {
0907         RawCluster* cluster = iterator.second;
0908 
0909         if (cluster->get_energy() < _reco_e_threshold)
0910         {
0911           continue;
0912         }
0913 
0914         _cluster_HCALIN_E[_nclusters_HCALIN] = cluster->get_energy();
0915         _cluster_HCALIN_NTower[_nclusters_HCALIN] = cluster->getNTowers();
0916         _cluster_HCALIN_Phi[_nclusters_HCALIN] = cluster->get_phi();
0917 
0918         // require vertex for cluster eta calculation
0919         if (vertexmap)
0920         {
0921           if (!vertexmap->empty())
0922           {
0923             SvtxVertex* vertex = (vertexmap->begin()->second);
0924             _cluster_HCALIN_Eta[_nclusters_HCALIN] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
0925           }
0926           else
0927           {
0928             _cluster_HCALIN_Eta[_nclusters_HCALIN] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
0929           };
0930         }
0931         else
0932         {
0933           _cluster_HCALIN_Eta[_nclusters_HCALIN] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
0934         };
0935 
0936         PHG4Particle* primary = clusterevalHCALIN->max_truth_primary_particle_by_energy(cluster);
0937 
0938         if (primary)
0939         {
0940           _cluster_HCALIN_trueID[_nclusters_HCALIN] = primary->get_track_id();
0941         }
0942         else
0943         {
0944           _cluster_HCALIN_trueID[_nclusters_HCALIN] = -10;
0945         }
0946 
0947         _nclusters_HCALIN++;
0948       }
0949     }
0950     else
0951     {
0952       std::cout << PHWHERE << " ERROR: Can't find " << clusternodeHCALIN << std::endl;
0953       // return;
0954     }
0955     if (Verbosity() > 0)
0956     {
0957       std::cout << "saved\t" << _nclusters_HCALIN << "\tHCALIN clusters" << std::endl;
0958     }
0959   }
0960   //------------------------
0961   // CLUSTERS HCALOUT
0962   //------------------------
0963   if (_do_HCALOUT && _do_CLUSTERS)
0964   {
0965     CaloRawClusterEval* clusterevalHCALOUT = _caloevalstackHCALOUT->get_rawcluster_eval();
0966     _nclusters_HCALOUT = 0;
0967     if (Verbosity() > 1)
0968     {
0969       std::cout << "CaloEvaluator::filling gcluster ntuple..." << std::endl;
0970     }
0971 
0972     std::string clusternodeHCALOUT = "CLUSTER_HCALOUT";
0973     RawClusterContainer* clustersHCALOUT = findNode::getClass<RawClusterContainer>(topNode, clusternodeHCALOUT);
0974     if (clustersHCALOUT)
0975     {
0976       // for every cluster
0977       for (const auto& iterator : clustersHCALOUT->getClustersMap())
0978       {
0979         RawCluster* cluster = iterator.second;
0980 
0981         if (cluster->get_energy() < _reco_e_threshold)
0982         {
0983           continue;
0984         }
0985 
0986         _cluster_HCALOUT_E[_nclusters_HCALOUT] = cluster->get_energy();
0987         _cluster_HCALOUT_NTower[_nclusters_HCALOUT] = cluster->getNTowers();
0988         _cluster_HCALOUT_Phi[_nclusters_HCALOUT] = cluster->get_phi();
0989 
0990         // require vertex for cluster eta calculation
0991         if (vertexmap)
0992         {
0993           if (!vertexmap->empty())
0994           {
0995             SvtxVertex* vertex = (vertexmap->begin()->second);
0996             _cluster_HCALOUT_Eta[_nclusters_HCALOUT] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
0997           }
0998           else
0999           {
1000             _cluster_HCALOUT_Eta[_nclusters_HCALOUT] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
1001           };
1002         }
1003         else
1004         {
1005           _cluster_HCALOUT_Eta[_nclusters_HCALOUT] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
1006         };
1007 
1008         PHG4Particle* primary = clusterevalHCALOUT->max_truth_primary_particle_by_energy(cluster);
1009 
1010         if (primary)
1011         {
1012           _cluster_HCALOUT_trueID[_nclusters_HCALOUT] = primary->get_track_id();
1013         }
1014         else
1015         {
1016           _cluster_HCALOUT_trueID[_nclusters_HCALOUT] = -10;
1017         }
1018 
1019         _nclusters_HCALOUT++;
1020       }
1021     }
1022     else
1023     {
1024       std::cout << PHWHERE << " ERROR: Can't find " << clusternodeHCALOUT << std::endl;
1025       // return;
1026     }
1027     if (Verbosity() > 0)
1028     {
1029       std::cout << "saved\t" << _nclusters_HCALOUT << "\tHCALOUT clusters" << std::endl;
1030     }
1031   }
1032   //------------------------
1033   // CLUSTERS CEMC
1034   //------------------------
1035   if (_do_CEMC && _do_CLUSTERS)
1036   {
1037     CaloRawClusterEval* clusterevalCEMC = _caloevalstackCEMC->get_rawcluster_eval();
1038     _nclusters_CEMC = 0;
1039     if (Verbosity() > 1)
1040     {
1041       std::cout << "CaloEvaluator::filling gcluster ntuple..." << std::endl;
1042     }
1043 
1044     std::string clusternodeCEMC = "CLUSTER_CEMC";
1045     RawClusterContainer* clustersCEMC = findNode::getClass<RawClusterContainer>(topNode, clusternodeCEMC);
1046     if (clustersCEMC)
1047     {
1048       // for every cluster
1049       for (const auto& iterator : clustersCEMC->getClustersMap())
1050       {
1051         RawCluster* cluster = iterator.second;
1052 
1053         if (cluster->get_energy() < _reco_e_threshold)
1054         {
1055           continue;
1056         }
1057 
1058         _cluster_CEMC_E[_nclusters_CEMC] = cluster->get_energy();
1059         _cluster_CEMC_NTower[_nclusters_CEMC] = cluster->getNTowers();
1060         _cluster_CEMC_Phi[_nclusters_CEMC] = cluster->get_phi();
1061 
1062         // require vertex for cluster eta calculation
1063         if (vertexmap)
1064         {
1065           if (!vertexmap->empty())
1066           {
1067             SvtxVertex* vertex = (vertexmap->begin()->second);
1068             _cluster_CEMC_Eta[_nclusters_CEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
1069           }
1070           else
1071           {
1072             _cluster_CEMC_Eta[_nclusters_CEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
1073           };
1074         }
1075         else
1076         {
1077           _cluster_CEMC_Eta[_nclusters_CEMC] = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0, 0, 0));
1078         };
1079 
1080         PHG4Particle* primary = clusterevalCEMC->max_truth_primary_particle_by_energy(cluster);
1081 
1082         if (primary)
1083         {
1084           _cluster_CEMC_trueID[_nclusters_CEMC] = primary->get_track_id();
1085         }
1086         else
1087         {
1088           _cluster_CEMC_trueID[_nclusters_CEMC] = -10;
1089         }
1090 
1091         _nclusters_CEMC++;
1092       }
1093     }
1094     else
1095     {
1096       std::cout << PHWHERE << " ERROR: Can't find " << clusternodeCEMC << std::endl;
1097       // return;
1098     }
1099     if (Verbosity() > 0)
1100     {
1101       std::cout << "saved\t" << _nclusters_CEMC << "\tCEMC clusters" << std::endl;
1102     }
1103   }
1104 
1105   //------------------------
1106   // TRACKS
1107   //------------------------
1108   if (_do_TRACKS)
1109   {
1110     _nTracks = 0;
1111     _nProjections = 0;
1112     // Loop over track maps, identifiy each source.
1113     // Although this configuration is fixed here, it doesn't require multiple sources.
1114     // It will only store them if they're available.
1115     std::vector<std::pair<std::string, TrackSource_t>> trackMapInfovec = {
1116         {"TrackMap", TrackSource_t::all},
1117         {"TrackMapInner", TrackSource_t::inner}};
1118     bool foundAtLeastOneTrackSource = false;
1119     for (const auto& trackMapInfo : trackMapInfovec)
1120     {
1121       SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, trackMapInfo.first);
1122       if (trackmap)
1123       {
1124         foundAtLeastOneTrackSource = true;
1125         int nTracksInASource = 0;
1126         if (Verbosity() > 0)
1127         {
1128           std::cout << "saving tracks for track map: " << trackMapInfo.first << std::endl;
1129         }
1130         for (SvtxTrackMap::ConstIter track_itr = trackmap->begin(); track_itr != trackmap->end(); track_itr++)
1131         {
1132           SvtxTrack_FastSim* track = dynamic_cast<SvtxTrack_FastSim*>(track_itr->second);
1133           if (track)
1134           {
1135             _track_ID[_nTracks] = track->get_id();
1136             _track_px[_nTracks] = track->get_px();
1137             _track_py[_nTracks] = track->get_py();
1138             _track_pz[_nTracks] = track->get_pz();
1139             // Ideally, would be dca3d_xy and dca3d_z, but these don't seem to be calculated properly in the
1140             // current (June 2021) simulations (they return NaN). So we take dca (seems to be ~ the 3d distance)
1141             // and dca_2d (seems to be ~ the distance in the transverse plane).
1142             // The names of the branches are based on the method names.
1143             _track_dca[_nTracks] = track->get_dca();
1144             _track_dca_2d[_nTracks] = track->get_dca2d();
1145             _track_trueID[_nTracks] = track->get_truth_track_id();
1146             _track_source[_nTracks] = static_cast<unsigned short>(trackMapInfo.second);
1147             if (_do_PROJECTIONS)
1148             {
1149               // find projections
1150               for (SvtxTrack::ConstStateIter trkstates = track->begin_states(); trkstates != track->end_states(); ++trkstates)
1151               {
1152                 if (Verbosity() > 1)
1153                 {
1154                   std::cout << __PRETTY_FUNCTION__ << " processing " << trkstates->second->get_name() << std::endl;
1155                 }
1156                 std::string trackStateName = trkstates->second->get_name();
1157                 if (Verbosity() > 1)
1158                 {
1159                   std::cout << __PRETTY_FUNCTION__ << " found " << trkstates->second->get_name() << std::endl;
1160                 }
1161                 int trackStateIndex = GetProjectionIndex(trackStateName);
1162                 if (trackStateIndex > -1)
1163                 {
1164                   // save true projection info to given branch
1165                   _track_TLP_true_x[_nProjections] = trkstates->second->get_pos(0);
1166                   _track_TLP_true_y[_nProjections] = trkstates->second->get_pos(1);
1167                   _track_TLP_true_z[_nProjections] = trkstates->second->get_pos(2);
1168                   _track_TLP_true_t[_nProjections] = trkstates->first;
1169                   _track_ProjLayer[_nProjections] = trackStateIndex;
1170                   _track_ProjTrackID[_nProjections] = _nTracks;
1171 
1172                   std::string nodename = "G4HIT_" + trkstates->second->get_name();
1173                   PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
1174                   if (hits)
1175                   {
1176                     if (Verbosity() > 1)
1177                     {
1178                       std::cout << __PRETTY_FUNCTION__ << " number of hits: " << hits->size() << std::endl;
1179                     }
1180                     PHG4HitContainer::ConstRange hit_range = hits->getHits();
1181                     for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
1182                     {
1183                       if (Verbosity() > 1)
1184                       {
1185                         std::cout << __PRETTY_FUNCTION__ << " checking hit id " << hit_iter->second->get_trkid() << " against " << track->get_truth_track_id() << std::endl;
1186                       }
1187                       if (hit_iter->second->get_trkid() - track->get_truth_track_id() == 0)
1188                       {
1189                         if (Verbosity() > 1)
1190                         {
1191                           std::cout << __PRETTY_FUNCTION__ << " found hit with id " << hit_iter->second->get_trkid() << std::endl;
1192                         }
1193                         // save reco projection info to given branch
1194                         _track_TLP_x[_nProjections] = hit_iter->second->get_x(0);
1195                         _track_TLP_y[_nProjections] = hit_iter->second->get_y(0);
1196                         _track_TLP_z[_nProjections] = hit_iter->second->get_z(0);
1197                         _track_TLP_t[_nProjections] = hit_iter->second->get_t(0);
1198                       }
1199                     }
1200                   }
1201                   else
1202                   {
1203                     if (Verbosity() > 1)
1204                     {
1205                       std::cout << __PRETTY_FUNCTION__ << " could not find " << nodename << std::endl;
1206                     }
1207                     continue;
1208                   }
1209                   _nProjections++;
1210                 }
1211               }
1212             }
1213             _nTracks++;
1214             nTracksInASource++;
1215           }
1216           else
1217           {
1218             if (Verbosity() > 0)  // Verbosity()
1219             {
1220               std::cout << "PHG4TrackFastSimEval::fill_track_tree - ignore track that is not a SvtxTrack_FastSim:";
1221               track_itr->second->identify();
1222             }
1223             continue;
1224           }
1225         }
1226         if (Verbosity() > 0)
1227         {
1228           std::cout << "saved\t" << nTracksInASource << "\ttracks from track map " << trackMapInfo.first << ". Total saved tracks: " << _nTracks << std::endl;
1229         }
1230       }
1231       else
1232       {
1233         if (Verbosity() > 0)
1234         {
1235           std::cout << PHWHERE << "SvtxTrackMap node with name '" << trackMapInfo.first << "' not found on node tree" << std::endl;
1236         }
1237       }
1238     }
1239     if (foundAtLeastOneTrackSource == false)
1240     {
1241       std::cout << PHWHERE << "Requested tracks, but found no sources on node tree. Returning" << std::endl;
1242       return;
1243     }
1244   }
1245   //------------------------
1246   // MC PARTICLES
1247   //------------------------
1248   _nMCPart = 0;
1249   if (_do_MCPARTICLES)
1250   {
1251     PHG4TruthInfoContainer* truthinfocontainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1252     if (truthinfocontainer)
1253     {
1254       if (Verbosity() > 0)
1255       {
1256         std::cout << "saving MC particles" << std::endl;
1257       }
1258       // GetParticleRange for all particles
1259       // GetPrimaryParticleRange for primary particles
1260       PHG4TruthInfoContainer::ConstRange range = truthinfocontainer->GetParticleRange();
1261       for (PHG4TruthInfoContainer::ConstIterator truth_itr = range.first; truth_itr != range.second; ++truth_itr)
1262       {
1263         PHG4Particle* g4particle = truth_itr->second;
1264         if (!g4particle)
1265         {
1266           continue;
1267         }
1268 
1269         int mcSteps = 0;
1270         PHG4Particle* g4particleMother = truth_itr->second;
1271         if (g4particle->get_parent_id() != 0)
1272         {
1273           while (g4particleMother->get_parent_id() != 0)
1274           {
1275             g4particleMother = truthinfocontainer->GetParticle(g4particleMother->get_parent_id());
1276             if (g4particleMother == nullptr)
1277             {
1278               break;
1279             }
1280             mcSteps += 1;
1281           }
1282         }
1283         if (mcSteps > _depth_MCstack)
1284         {
1285           continue;
1286         }
1287 
1288         // evaluating true primary vertex
1289         if (_do_VERTEX && _nMCPart == 0)
1290         {
1291           PHG4VtxPoint* vtx = truthinfocontainer->GetVtx(g4particle->get_vtx_id());
1292           if (vtx)
1293           {
1294             _vertex_true_x = vtx->get_x();
1295             _vertex_true_y = vtx->get_y();
1296             _vertex_true_z = vtx->get_z();
1297           }
1298         }
1299 
1300         // in case of all MC particles, make restrictions on the secondary selection
1301         // if(g4particle->get_track_id()<0 && g4particle->get_e()<0.5) continue;
1302         // primary (g4particle->get_parent_id() == 0) selection via:
1303         // if(gtrackID < 0) continue;
1304 
1305         // using the e threshold also for the truth particles gets rid of all the low energy secondary particles
1306         if (g4particle->get_e() < _reco_e_threshold)
1307         {
1308           continue;
1309         }
1310 
1311         _mcpart_ID[_nMCPart] = g4particle->get_track_id();
1312         _mcpart_ID_parent[_nMCPart] = g4particle->get_parent_id();
1313         _mcpart_PDG[_nMCPart] = g4particle->get_pid();
1314         _mcpart_E[_nMCPart] = g4particle->get_e();
1315         _mcpart_px[_nMCPart] = g4particle->get_px();
1316         _mcpart_py[_nMCPart] = g4particle->get_py();
1317         _mcpart_pz[_nMCPart] = g4particle->get_pz();
1318         // BCID added for G4Particle --  HEPMC particle matching
1319         _mcpart_BCID[_nMCPart] = g4particle->get_barcode();
1320         // TVector3 projvec(_mcpart_px[0],_mcpart_py[0],_mcpart_pz[0]);
1321         // float projeta = projvec.Eta();
1322         _nMCPart++;
1323       }
1324       if (Verbosity() > 0)
1325       {
1326         std::cout << "saved\t" << _nMCPart << "\tMC particles" << std::endl;
1327       }
1328     }
1329     else
1330     {
1331       if (Verbosity() > 0)
1332       {
1333         std::cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree" << std::endl;
1334       }
1335       return;
1336     }
1337   }
1338 
1339   //------------------------
1340   // HEPMC
1341   //------------------------
1342   _nHepmcp = 0;
1343   if (_do_HEPMC)
1344   {
1345     PHHepMCGenEventMap* hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
1346     if (hepmceventmap)
1347     {
1348       if (Verbosity() > 0)
1349       {
1350         std::cout << "saving HepMC output" << std::endl;
1351       }
1352       if (Verbosity() > 0)
1353       {
1354         hepmceventmap->Print();
1355       }
1356 
1357       for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
1358            eventIter != hepmceventmap->end();
1359            ++eventIter)
1360       {
1361         PHHepMCGenEvent* hepmcevent = eventIter->second;
1362 
1363         if (hepmcevent)
1364         {
1365           HepMC::GenEvent* truthevent = hepmcevent->getEvent();
1366           if (!truthevent)
1367           {
1368             std::cout << PHWHERE
1369                       << "no evt pointer under phhepmvgeneventmap found "
1370                       << std::endl;
1371             return;
1372           }
1373 
1374           HepMC::PdfInfo* pdfinfo = truthevent->pdf_info();
1375 
1376           //     m_partid1 = pdfinfo->id1();
1377           // m_partid2 = pdfinfo->id2();
1378           _hepmcp_x1 = pdfinfo->x1();
1379           _hepmcp_x2 = pdfinfo->x2();
1380 
1381           // m_mpi = truthevent->mpi();
1382 
1383           _hepmcp_procid = truthevent->signal_process_id();
1384 
1385           if (Verbosity() > 2)
1386           {
1387             std::cout << " Iterating over an event" << std::endl;
1388           }
1389           for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
1390                iter != truthevent->particles_end();
1391                ++iter)
1392           {
1393             _hepmcp_E[_nHepmcp] = (*iter)->momentum().e();
1394             _hepmcp_PDG[_nHepmcp] = (*iter)->pdg_id();
1395             _hepmcp_px[_nHepmcp] = (*iter)->momentum().px();
1396             _hepmcp_py[_nHepmcp] = (*iter)->momentum().py();
1397             _hepmcp_pz[_nHepmcp] = (*iter)->momentum().pz();
1398             _hepmcp_status[_nHepmcp] = (*iter)->status();
1399             _hepmcp_BCID[_nHepmcp] = (*iter)->barcode();
1400             _hepmcp_m2[_nHepmcp] = 0;
1401             _hepmcp_m1[_nHepmcp] = 0;
1402             if ((*iter)->production_vertex())
1403             {
1404               for (HepMC::GenVertex::particle_iterator mother = (*iter)->production_vertex()->particles_begin(HepMC::parents);
1405                    mother != (*iter)->production_vertex()->particles_end(HepMC::parents);
1406                    ++mother)
1407               {
1408                 _hepmcp_m2[_nHepmcp] = (*mother)->barcode();
1409                 if (_hepmcp_m1[_nHepmcp] == 0)
1410                 {
1411                   _hepmcp_m1[_nHepmcp] = (*mother)->barcode();
1412                 }
1413               }
1414             }
1415             if (Verbosity() > 2)
1416             {
1417               std::cout << "nHepmcp " << _nHepmcp << "\tPDG " << _hepmcp_PDG[_nHepmcp] << "\tEnergy " << _hepmcp_E[_nHepmcp] << "\tbarcode " << _hepmcp_BCID[_nHepmcp] << "\tMother1 " << _hepmcp_m1[_nHepmcp] << "\tMother2 " << _hepmcp_m2[_nHepmcp] << std::endl;
1418             }
1419             _nHepmcp++;
1420           }
1421         }
1422       }
1423     }
1424     else
1425     {
1426       if (Verbosity() > 0)
1427       {
1428         std::cout << PHWHERE << " PHHepMCGenEventMap node not found on node tree" << std::endl;
1429       }
1430       return;
1431     }
1432   }  // hepmc
1433 
1434   _event_tree->Fill();
1435 
1436   if (Verbosity() > 0)
1437   {
1438     std::cout << "Resetting buffer ..." << std::endl;
1439   }
1440   resetBuffer();
1441   if (Verbosity() > 0)
1442   {
1443     std::cout << "EventEvaluator buffer reset" << std::endl;
1444   }
1445   return;
1446 }
1447 
1448 int EventEvaluator::End(PHCompositeNode* /*topNode*/)
1449 {
1450   _tfile->cd();
1451 
1452   _event_tree->Write();
1453 
1454   _tfile->Close();
1455 
1456   delete _tfile;
1457 
1458   if (_do_GEOMETRY)
1459   {
1460     _tfile_geometry->cd();
1461 
1462     _geometry_tree->Write();
1463 
1464     _tfile_geometry->Close();
1465 
1466     delete _tfile_geometry;
1467   }
1468   if (Verbosity() > 0)
1469   {
1470     std::cout << "========================= " << Name() << "::End() ============================" << std::endl;
1471     std::cout << " " << _ievent << " events of output written to: " << _filename << std::endl;
1472     std::cout << "===========================================================================" << std::endl;
1473   }
1474 
1475   delete _caloevalstackHCALIN;
1476   delete _caloevalstackHCALOUT;
1477   delete _caloevalstackCEMC;
1478 
1479   return Fun4AllReturnCodes::EVENT_OK;
1480 }
1481 
1482 int EventEvaluator::GetProjectionIndex(const std::string& projname)
1483 {
1484   if (projname.find("HCALIN") != std::string::npos)
1485   {
1486     return 1;
1487   }
1488   if (projname.find("HCALOUT") != std::string::npos)
1489   {
1490     return 2;
1491   }
1492   if (projname.find("CEMC") != std::string::npos)
1493   {
1494     return 3;
1495   }
1496 
1497   return -1;
1498 
1499   return -1;
1500 }
1501 
1502 std::string EventEvaluator::GetProjectionNameFromIndex(int projindex)
1503 {
1504   switch (projindex)
1505   {
1506   case 1:
1507     return "HCALIN";
1508   case 2:
1509     return "HCALOUT";
1510   case 3:
1511     return "CEMC";
1512   default:
1513     return "NOTHING";
1514   }
1515 }
1516 
1517 void EventEvaluator::resetGeometryArrays()
1518 {
1519   for (Int_t igeo = 0; igeo < _calo_towers_N; igeo++)
1520   {
1521     _calo_towers_iEta[_calo_towers_N] = -10000;
1522     _calo_towers_iPhi[_calo_towers_N] = -10000;
1523     _calo_towers_Eta[_calo_towers_N] = -10000;
1524     _calo_towers_Phi[_calo_towers_N] = -10000;
1525     _calo_towers_x[_calo_towers_N] = -10000;
1526     _calo_towers_y[_calo_towers_N] = -10000;
1527     _calo_towers_z[_calo_towers_N] = -10000;
1528   }
1529   _calo_ID = -1;
1530   _calo_towers_N = 0;
1531 }
1532 void EventEvaluator::resetBuffer()
1533 {
1534   if (_do_store_event_info)
1535   {
1536     _cross_section = 0;
1537     _event_weight = 0;
1538     _n_generator_accepted = 0;
1539     if (Verbosity() > 0)
1540     {
1541       std::cout << "\t... event info variables reset" << std::endl;
1542     }
1543   }
1544   if (_do_VERTEX)
1545   {
1546     _vertex_x = -1000;
1547     _vertex_y = -1000;
1548     _vertex_z = -1000;
1549     _vertex_NCont = 0;
1550     _vertex_true_x = -1000;
1551     _vertex_true_y = -1000;
1552     _vertex_true_z = -1000;
1553     if (Verbosity() > 0)
1554     {
1555       std::cout << "\t... vertex variables reset" << std::endl;
1556     }
1557   }
1558   if (_do_HITS)
1559   {
1560     _nHitsLayers = 0;
1561     for (Int_t ihit = 0; ihit < _maxNHits; ihit++)
1562     {
1563       _hits_layerID[ihit] = 0;
1564       _hits_trueID[ihit] = 0;
1565       _hits_x[ihit] = 0;
1566       _hits_y[ihit] = 0;
1567       _hits_z[ihit] = 0;
1568       _hits_t[ihit] = 0;
1569     }
1570     if (Verbosity() > 0)
1571     {
1572       std::cout << "\t... hit variables reset" << std::endl;
1573     }
1574   }
1575   if (_do_CEMC)
1576   {
1577     _nTowers_CEMC = 0;
1578     for (Int_t itow = 0; itow < _maxNTowersCentral; itow++)
1579     {
1580       _tower_CEMC_E[itow] = 0;
1581       _tower_CEMC_iEta[itow] = 0;
1582       _tower_CEMC_iPhi[itow] = 0;
1583       _tower_CEMC_trueID[itow] = 0;
1584     }
1585     if (_do_CLUSTERS)
1586     {
1587       _nclusters_CEMC = 0;
1588       for (Int_t itow = 0; itow < _maxNclustersCentral; itow++)
1589       {
1590         _cluster_CEMC_E[itow] = 0;
1591         _cluster_CEMC_Eta[itow] = 0;
1592         _cluster_CEMC_Phi[itow] = 0;
1593         _cluster_CEMC_NTower[itow] = 0;
1594         _cluster_CEMC_trueID[itow] = 0;
1595       }
1596     }
1597     if (Verbosity() > 0)
1598     {
1599       std::cout << "\t... CEMC variables reset" << std::endl;
1600     }
1601   }
1602   if (_do_HCALIN)
1603   {
1604     _nTowers_HCALIN = 0;
1605     for (Int_t itow = 0; itow < _maxNTowersCentral; itow++)
1606     {
1607       _tower_HCALIN_E[itow] = 0;
1608       _tower_HCALIN_iEta[itow] = 0;
1609       _tower_HCALIN_iPhi[itow] = 0;
1610       _tower_HCALIN_trueID[itow] = 0;
1611     }
1612     if (_do_CLUSTERS)
1613     {
1614       _nclusters_HCALIN = 0;
1615       for (Int_t itow = 0; itow < _maxNclustersCentral; itow++)
1616       {
1617         _cluster_HCALIN_E[itow] = 0;
1618         _cluster_HCALIN_Eta[itow] = 0;
1619         _cluster_HCALIN_Phi[itow] = 0;
1620         _cluster_HCALIN_NTower[itow] = 0;
1621         _cluster_HCALIN_trueID[itow] = 0;
1622       }
1623     }
1624     if (Verbosity() > 0)
1625     {
1626       std::cout << "\t... HCALIN variables reset" << std::endl;
1627     }
1628   }
1629   if (_do_HCALOUT)
1630   {
1631     if (Verbosity() > 0)
1632     {
1633       std::cout << "\t... resetting HCALOUT variables" << std::endl;
1634     }
1635     _nTowers_HCALOUT = 0;
1636     for (Int_t itow = 0; itow < _maxNTowersCentral; itow++)
1637     {
1638       _tower_HCALOUT_E[itow] = 0;
1639       _tower_HCALOUT_iEta[itow] = 0;
1640       _tower_HCALOUT_iPhi[itow] = 0;
1641       _tower_HCALOUT_trueID[itow] = 0;
1642     }
1643     if (_do_CLUSTERS)
1644     {
1645       _nclusters_HCALOUT = 0;
1646       for (Int_t itow = 0; itow < _maxNclustersCentral; itow++)
1647       {
1648         _cluster_HCALOUT_E[itow] = 0;
1649         _cluster_HCALOUT_Eta[itow] = 0;
1650         _cluster_HCALOUT_Phi[itow] = 0;
1651         _cluster_HCALOUT_NTower[itow] = 0;
1652         _cluster_HCALOUT_trueID[itow] = 0;
1653       }
1654     }
1655     if (Verbosity() > 0)
1656     {
1657       std::cout << "\t... HCALOUT variables reset" << std::endl;
1658     }
1659   }
1660   if (_do_TRACKS)
1661   {
1662     if (Verbosity() > 0)
1663     {
1664       std::cout << "\t... resetting Track variables" << std::endl;
1665     }
1666     _nTracks = 0;
1667     for (Int_t itrk = 0; itrk < _maxNTracks; itrk++)
1668     {
1669       _track_ID[itrk] = 0;
1670       _track_trueID[itrk] = 0;
1671       _track_px[itrk] = 0;
1672       _track_py[itrk] = 0;
1673       _track_pz[itrk] = 0;
1674       _track_dca[itrk] = 0;
1675       _track_dca_2d[itrk] = 0;
1676       _track_source[itrk] = 0;
1677     }
1678     if (_do_PROJECTIONS)
1679     {
1680       _nProjections = 0;
1681       for (Int_t iproj = 0; iproj < _maxNProjections; iproj++)
1682       {
1683         _track_ProjLayer[iproj] = -1;
1684         _track_ProjTrackID[iproj] = 0;
1685         _track_TLP_x[iproj] = 0;
1686         _track_TLP_y[iproj] = 0;
1687         _track_TLP_z[iproj] = 0;
1688         _track_TLP_t[iproj] = 0;
1689         _track_TLP_true_x[iproj] = 0;
1690         _track_TLP_true_y[iproj] = 0;
1691         _track_TLP_true_z[iproj] = 0;
1692         _track_TLP_true_t[iproj] = 0;
1693       }
1694     }
1695     if (Verbosity() > 0)
1696     {
1697       std::cout << "\t... track variables reset" << std::endl;
1698     }
1699   }
1700   if (_do_MCPARTICLES)
1701   {
1702     _nMCPart = 0;
1703     for (Int_t imcpart = 0; imcpart < _maxNMCPart; imcpart++)
1704     {
1705       _mcpart_ID[imcpart] = 0;
1706       _mcpart_ID_parent[imcpart] = 0;
1707       _mcpart_PDG[imcpart] = 0;
1708       _mcpart_E[imcpart] = 0;
1709       _mcpart_px[imcpart] = 0;
1710       _mcpart_py[imcpart] = 0;
1711       _mcpart_pz[imcpart] = 0;
1712       _mcpart_BCID[imcpart] = -10;
1713     }
1714   }
1715 
1716   if (_do_HEPMC)
1717   {
1718     _nHepmcp = 0;
1719     for (Int_t iHepmcp = 0; iHepmcp < _maxNHepmcp; iHepmcp++)
1720     {
1721       _hepmcp_E[iHepmcp] = 0;
1722       _hepmcp_PDG[iHepmcp] = 0;
1723       _hepmcp_px[iHepmcp] = 0;
1724       _hepmcp_py[iHepmcp] = 0;
1725       _hepmcp_pz[iHepmcp] = 0;
1726       _hepmcp_status[iHepmcp] = -10;
1727       _hepmcp_BCID[iHepmcp] = 0;
1728       _hepmcp_m2[iHepmcp] = 0;
1729       _hepmcp_m1[iHepmcp] = 0;
1730     }
1731     if (Verbosity() > 0)
1732     {
1733       std::cout << "\t... MC variables reset" << std::endl;
1734     }
1735   }
1736 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!