Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:17:55 

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