Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2026-08-31 08:14:23

0001 //____________________________________________________________________________..
0002 //
0003 // This is a template for a Fun4All SubsysReco module with all methods from the
0004 // $OFFLINE_MAIN/include/fun4all/SubsysReco.h baseclass
0005 // You do not have to implement all of them, you can just remove unused methods
0006 // here and in VertexCompare.h.
0007 //
0008 // VertexCompare(const std::string &name = "VertexCompare")
0009 // everything is keyed to VertexCompare, duplicate names do work but it makes
0010 // e.g. finding culprits in logs difficult or getting a pointer to the module
0011 // from the command line
0012 //
0013 // VertexCompare::~VertexCompare()
0014 // this is called when the Fun4AllServer is deleted at the end of running. Be
0015 // mindful what you delete - you do loose ownership of object you put on the node tree
0016 //
0017 // int VertexCompare::Init(PHCompositeNode *topNode)
0018 // This method is called when the module is registered with the Fun4AllServer. You
0019 // can create historgrams here or put objects on the node tree but be aware that
0020 // modules which haven't been registered yet did not put antyhing on the node tree
0021 //
0022 // int VertexCompare::InitRun(PHCompositeNode *topNode)
0023 // This method is called when the first event is read (or generated). At
0024 // this point the run number is known (which is mainly interesting for raw data
0025 // processing). Also all objects are on the node tree in case your module's action
0026 // depends on what else is around. Last chance to put nodes under the DST Node
0027 // We mix events during readback if branches are added after the first event
0028 //
0029 // int VertexCompare::process_event(PHCompositeNode *topNode)
0030 // called for every event. Return codes trigger actions, you find them in
0031 // $OFFLINE_MAIN/include/fun4all/Fun4AllReturnCodes.h
0032 //   everything is good:
0033 //     return Fun4AllReturnCodes::EVENT_OK
0034 //   abort event reconstruction, clear everything and process next event:
0035 //     return Fun4AllReturnCodes::ABORT_EVENT;
0036 //   proceed but do not save this event in output (needs output manager setting):
0037 //     return Fun4AllReturnCodes::DISCARD_EVENT;
0038 //   abort processing:
0039 //     return Fun4AllReturnCodes::ABORT_RUN
0040 // all other integers will lead to an error and abort of processing
0041 //
0042 // int VertexCompare::ResetEvent(PHCompositeNode *topNode)
0043 // If you have internal data structures (arrays, stl containers) which needs clearing
0044 // after each event, this is the place to do that. The nodes under the DST node are cleared
0045 // by the framework
0046 //
0047 // int VertexCompare::EndRun(const int runnumber)
0048 // This method is called at the end of a run when an event from a new run is
0049 // encountered. Useful when analyzing multiple runs (raw data). Also called at
0050 // the end of processing (before the End() method)
0051 //
0052 // int VertexCompare::End(PHCompositeNode *topNode)
0053 // This is called at the end of processing. It needs to be called by the macro
0054 // by Fun4AllServer::End(), so do not forget this in your macro
0055 //
0056 // int VertexCompare::Reset(PHCompositeNode *topNode)
0057 // not really used - it is called before the dtor is called
0058 //
0059 // void VertexCompare::Print(const std::string &what) const
0060 // Called from the command line - useful to print information when you need it
0061 //
0062 //____________________________________________________________________________..
0063 
0064 #include "VertexCompare.h"
0065 
0066 #include <fun4all/Fun4AllReturnCodes.h>
0067 
0068 #include <phool/PHCompositeNode.h>
0069 #include <phool/getClass.h>
0070 #include <phool/sphenix_constants.h>
0071 
0072 #include <phhepmc/PHHepMCGenEvent.h>
0073 #include <phhepmc/PHHepMCGenEventMap.h>
0074 
0075 #include <trackbase/ActsGeometry.h>
0076 #include <trackbase/TrkrCluster.h>
0077 #include <trackbase/TrkrClusterContainer.h>
0078 #include <trackbase/TrkrClusterCrossingAssoc.h>
0079 #include <trackbase/TrkrClusterHitAssoc.h>
0080 #include <trackbase_historic/SvtxTrack.h>
0081 #include <trackbase_historic/SvtxTrackMap.h>
0082 #include <trackbase_historic/SvtxPHG4ParticleMap.h>
0083 #include <trackbase_historic/TrackAnalysisUtils.h>
0084 #include <trackbase_historic/TrackSeed.h>
0085 #include <trackbase_historic/TrackSeedContainer.h>
0086 #include <trackbase_historic/TrackSeedHelper.h>
0087 
0088 #include <g4detectors/PHG4TpcGeom.h>
0089 #include <g4detectors/PHG4TpcGeomContainer.h>
0090 
0091 #include <mvtx/SegmentationAlpide.h>
0092 #include <mvtx/MvtxPixelDefs.h>
0093 
0094 #include <trackbase/InttDefs.h>
0095 #include <trackbase/MvtxDefs.h>
0096 #include <trackbase/TrkrDefs.h>
0097 
0098 #include <ffarawobjects/Gl1Packet.h>
0099 #include <ffarawobjects/Gl1RawHit.h>
0100 #include <ffaobjects/EventHeader.h>
0101 
0102 #include <globalvertex/GlobalVertex.h>
0103 #include <globalvertex/GlobalVertexMap.h>
0104 #include <globalvertex/MbdVertex.h>
0105 #include <globalvertex/MbdVertexMap.h>
0106 #include <globalvertex/SvtxVertex.h>
0107 #include <globalvertex/SvtxVertexMap.h>
0108 
0109 #include <mbd/MbdOut.h>
0110 #include <mbd/MbdDefs.h>
0111 #include <mbd/MbdPmtContainer.h>
0112 #include <mbd/MbdPmtHit.h>
0113 
0114 #include <calotrigger/MinimumBiasInfo.h>
0115 #include <centrality/CentralityInfo.h>
0116 
0117 #include "g4eval/SvtxClusterEval.h"
0118 #include <g4eval/SvtxEvalStack.h>
0119 #include <g4eval/SvtxHitEval.h>
0120 #include <g4eval/SvtxTruthEval.h>
0121 #include <g4main/PHG4Hit.h>
0122 #include <g4main/PHG4HitContainer.h>
0123 #include <g4main/PHG4Particle.h>
0124 #include <g4main/PHG4MCProcessDefs.h>
0125 #include <g4main/PHG4TruthInfoContainer.h>
0126 #include <g4main/PHG4VtxPoint.h>
0127 
0128 #include <HepMC/GenEvent.h>
0129 
0130 #include <algorithm>
0131 #include <climits>
0132 #include <cmath>
0133 #include <map>
0134 
0135 #include <TDatabasePDG.h>
0136 #include <TParticlePDG.h>
0137 
0138 namespace
0139 {
0140 constexpr int kMbdSouthArm = 0;
0141 constexpr int kMbdNorthArm = 1;
0142 constexpr int kMbdPmtsPerArm = MbdDefs::MBD_N_PMT / MbdDefs::MBD_N_ARMS;
0143 constexpr float kMbdPmtChargeThreshold = 0.25F;
0144 
0145 template <class Container> void Clean(Container &c) { Container().swap(c); }
0146 } // namespace
0147 
0148 GlobalVertex::VTXTYPE trkType = GlobalVertex::SVTX;
0149 GlobalVertex::VTXTYPE mbdType = GlobalVertex::MBD;
0150 //____________________________________________________________________________..
0151 VertexCompare::VertexCompare(const std::string &name)
0152     : SubsysReco(name)
0153 {
0154     std::cout << "VertexCompare::VertexCompare(const std::string &name) Calling ctor" << std::endl;
0155 }
0156 
0157 //____________________________________________________________________________..
0158 VertexCompare::~VertexCompare()
0159 {
0160     delete svtx_evalstack;
0161     svtx_evalstack = nullptr;
0162     clustereval = nullptr;
0163     hiteval = nullptr;
0164     truth_eval = nullptr;
0165 
0166     if (outFile)
0167     {
0168         if (outFile->IsOpen())
0169         {
0170             outFile->Close();
0171         }
0172         delete outFile;
0173         outFile = nullptr;
0174         outTree = nullptr;
0175     }
0176 
0177     std::cout << "VertexCompare::~VertexCompare() Calling dtor" << std::endl;
0178 }
0179 
0180 //____________________________________________________________________________..
0181 int VertexCompare::Init(PHCompositeNode *topNode)
0182 {
0183     outFile = new TFile(outFileName.c_str(), "RECREATE");
0184     outTree = new TTree("VTX", "VTX");
0185     outTree->SetAutoFlush(-10 * 1024 * 1024);
0186     outTree->SetAutoSave(-100 * 1024 * 1024);
0187     outTree->SetMaxVirtualSize(64 * 1024 * 1024);
0188 
0189     outTree->Branch("counter", &counter, "counter/I");
0190     outTree->Branch("is_min_bias", &is_min_bias);
0191     outTree->Branch("firedTriggers", &firedTriggers);
0192     outTree->Branch("gl1bco", &gl1bco);
0193     outTree->Branch("gl1BunchCrossing", &gl1BunchCrossing);
0194     outTree->Branch("bcotr", &bcotr);
0195     outTree->Branch("centrality_mbd", &centrality_mbd_);
0196     outTree->Branch("n_MBDVertex", &n_MBDVertex, "n_MBDVertex/i");
0197     outTree->Branch("mbdVertex", &mbdVertex);
0198     outTree->Branch("mbdVertexId", &mbdVertexId);
0199     outTree->Branch("mbdVertexCrossing", &mbdVertexCrossing);
0200     outTree->Branch("MBD_charge_sum", &MBD_charge_sum);
0201     outTree->Branch("mbd_north_npmt", &mbd_north_npmt);
0202     outTree->Branch("mbd_south_npmt", &mbd_south_npmt);
0203     outTree->Branch("mbd_south_charge_sum", &mbd_south_charge_sum);
0204     outTree->Branch("mbd_north_charge_sum", &mbd_north_charge_sum);
0205     outTree->Branch("mbd_nhitsoverths_south", &mbd_nhitsoverths_south);
0206     outTree->Branch("mbd_nhitsoverths_north", &mbd_nhitsoverths_north);
0207     outTree->Branch("nSvtxVertices", &nSvtxVertices);
0208     outTree->Branch("nSvtxVertices_validCrossing", &nSvtxVertices_validCrossing);
0209     outTree->Branch("trackerVertexId", &trackerVertexId);
0210     outTree->Branch("trackerVertexX", &trackerVertexX);
0211     outTree->Branch("trackerVertexY", &trackerVertexY);
0212     outTree->Branch("trackerVertexZ", &trackerVertexZ);
0213     outTree->Branch("trackerVertexChisq", &trackerVertexChisq);
0214     outTree->Branch("trackerVertexNdof", &trackerVertexNdof);
0215     outTree->Branch("trackerVertexNTracks", &trackerVertexNTracks);
0216     outTree->Branch("trackerVertexCrossing", &trackerVertexCrossing);
0217     outTree->Branch("trackerVertexTrackIDs", &trackerVertexTrackIDs);
0218     outTree->Branch("nTracks", &nTracks, "nTracks/i");
0219     // outTree->Branch("n_TRKVertex", &n_TRKVertex, "n_TRKVertex/i");
0220     outTree->Branch("hasMBD", &hasMBD, "hasMBD/O");
0221     outTree->Branch("hasTRK", &hasTRK, "hasTRK/O");
0222 
0223     if (writeTrackBranches_)
0224     {
0225         outTree->Branch("nRecoTracks", &nRecoTracks);
0226         outTree->Branch("track_deltapt", &track_deltapt);
0227         outTree->Branch("track_deltaeta", &track_deltaeta);
0228         outTree->Branch("track_deltaphi", &track_deltaphi);
0229         outTree->Branch("track_nhits", &track_nhits);
0230         outTree->Branch("track_nmaps", &track_nmaps);
0231         outTree->Branch("track_nintt", &track_nintt);
0232         outTree->Branch("track_ntpc", &track_ntpc);
0233         outTree->Branch("track_nmms", &track_nmms);
0234         outTree->Branch("track_ntpc1", &track_ntpc1);
0235         outTree->Branch("track_ntpc11", &track_ntpc11);
0236         outTree->Branch("track_ntpc2", &track_ntpc2);
0237         outTree->Branch("track_ntpc3", &track_ntpc3);
0238         outTree->Branch("track_pidedx", &track_pidedx);
0239         outTree->Branch("track_kdedx", &track_kdedx);
0240         outTree->Branch("track_prdedx", &track_prdedx);
0241         outTree->Branch("track_vx", &track_vx);
0242         outTree->Branch("track_vy", &track_vy);
0243         outTree->Branch("track_vz", &track_vz);
0244         outTree->Branch("track_dca2d", &track_dca2d);
0245         outTree->Branch("track_dca2dsigma", &track_dca2dsigma);
0246         outTree->Branch("track_dca3dxy", &track_dca3dxy);
0247         outTree->Branch("track_dca3dxysigma", &track_dca3dxysigma);
0248         outTree->Branch("track_dca3dz", &track_dca3dz);
0249         outTree->Branch("track_dca3dzsigma", &track_dca3dzsigma);
0250         outTree->Branch("track_hlxpt", &track_hlxpt);
0251         outTree->Branch("track_hlxeta", &track_hlxeta);
0252         outTree->Branch("track_hlxphi", &track_hlxphi);
0253         outTree->Branch("track_hlxX0", &track_hlxX0);
0254         outTree->Branch("track_hlxY0", &track_hlxY0);
0255         outTree->Branch("track_hlxZ0", &track_hlxZ0);
0256         outTree->Branch("track_hlxcharge", &track_hlxcharge);
0257         outTree->Branch("track_id", &track_id);
0258         outTree->Branch("track_x", &track_x);
0259         outTree->Branch("track_y", &track_y);
0260         outTree->Branch("track_z", &track_z);
0261         outTree->Branch("track_px", &track_px);
0262         outTree->Branch("track_py", &track_py);
0263         outTree->Branch("track_pz", &track_pz);
0264         outTree->Branch("track_pt", &track_pt);
0265         outTree->Branch("track_eta", &track_eta);
0266         outTree->Branch("track_phi", &track_phi);
0267         outTree->Branch("track_dedx", &track_dedx);
0268         outTree->Branch("track_charge", &track_charge);
0269         outTree->Branch("track_crossing", &track_crossing);
0270         outTree->Branch("track_vertex_id", &track_vertex_id);
0271         outTree->Branch("track_chisq", &track_chisq);
0272         outTree->Branch("track_ndf", &track_ndf);
0273         outTree->Branch("track_quality", &track_quality);
0274         outTree->Branch("track_silseed_id", &track_silseed_id);
0275         outTree->Branch("track_silseed_x", &track_silseed_x);
0276         outTree->Branch("track_silseed_y", &track_silseed_y);
0277         outTree->Branch("track_silseed_z", &track_silseed_z);
0278         outTree->Branch("track_silseed_pt", &track_silseed_pt);
0279         outTree->Branch("track_silseed_eta", &track_silseed_eta);
0280         outTree->Branch("track_silseed_phi", &track_silseed_phi);
0281         outTree->Branch("track_silseed_crossing", &track_silseed_crossing);
0282         outTree->Branch("track_silseed_charge", &track_silseed_charge);
0283         outTree->Branch("track_silseed_nMvtx", &track_silseed_nMvtx);
0284         outTree->Branch("track_silseed_nIntt", &track_silseed_nIntt);
0285         outTree->Branch("track_silseed_clusterKeys", &track_silseed_clusterKeys);
0286         outTree->Branch("track_cluster_layer", &track_cluster_layer);
0287         outTree->Branch("track_cluster_globalX", &track_cluster_globalX);
0288         outTree->Branch("track_cluster_globalY", &track_cluster_globalY);
0289         outTree->Branch("track_cluster_globalZ", &track_cluster_globalZ);
0290         outTree->Branch("track_cluster_phi", &track_cluster_phi);
0291         outTree->Branch("track_cluster_eta", &track_cluster_eta);
0292         outTree->Branch("track_cluster_r", &track_cluster_r);
0293     }
0294 
0295     // silicon seed information
0296     outTree->Branch("nTotalSilSeeds", &nTotalSilSeeds);
0297     outTree->Branch("nSilSeedsValidCrossing", &nSilSeedsValidCrossing);
0298     outTree->Branch("silseed_id", &silseed_id);
0299     outTree->Branch("silseed_assocVtxId", &silseed_assocVtxId);
0300     outTree->Branch("silseed_x", &silseed_x);
0301     outTree->Branch("silseed_y", &silseed_y);
0302     outTree->Branch("silseed_z", &silseed_z);
0303     outTree->Branch("silseed_pt", &silseed_pt);
0304     outTree->Branch("silseed_eta", &silseed_eta);
0305     outTree->Branch("silseed_phi", &silseed_phi);
0306     outTree->Branch("silseed_eta_vtx", &silseed_eta_vtx);
0307     outTree->Branch("silseed_phi_vtx", &silseed_phi_vtx);
0308     outTree->Branch("silseed_crossing", &silseed_crossing);
0309     outTree->Branch("silseed_charge", &silseed_charge);
0310     outTree->Branch("silseed_nMvtx", &silseed_nMvtx);
0311     outTree->Branch("silseed_nIntt", &silseed_nIntt);
0312     outTree->Branch("silseed_clusterKeys", &silseed_clusterKeys);
0313     outTree->Branch("silseed_cluster_layer", &silseed_cluster_layer);
0314     outTree->Branch("silseed_cluster_globalX", &silseed_cluster_globalX);
0315     outTree->Branch("silseed_cluster_globalY", &silseed_cluster_globalY);
0316     outTree->Branch("silseed_cluster_globalZ", &silseed_cluster_globalZ);
0317     outTree->Branch("silseed_cluster_phi", &silseed_cluster_phi);
0318     outTree->Branch("silseed_cluster_eta", &silseed_cluster_eta);
0319     outTree->Branch("silseed_cluster_r", &silseed_cluster_r);
0320     outTree->Branch("silseed_cluster_phiSize", &silseed_cluster_phiSize);
0321     outTree->Branch("silseed_cluster_zSize", &silseed_cluster_zSize);
0322     outTree->Branch("silseed_cluster_strobeID", &silseed_cluster_strobeID);
0323     outTree->Branch("silseed_cluster_timeBucketID", &silseed_cluster_timeBucketID); // only for INTT
0324 
0325     if (writeTpcSeedBranches_)
0326     {
0327         outTree->Branch("nTotalTpcSeeds", &nTotalTpcSeeds);
0328         outTree->Branch("tpcseed_id", &tpcseed_id);
0329         outTree->Branch("tpcseed_x", &tpcseed_x);
0330         outTree->Branch("tpcseed_y", &tpcseed_y);
0331         outTree->Branch("tpcseed_z", &tpcseed_z);
0332         outTree->Branch("tpcseed_pt", &tpcseed_pt);
0333         outTree->Branch("tpcseed_eta", &tpcseed_eta);
0334         outTree->Branch("tpcseed_phi", &tpcseed_phi);
0335         outTree->Branch("tpcseed_crossing", &tpcseed_crossing);
0336         outTree->Branch("tpcseed_crossing_estimate", &tpcseed_crossing_estimate);
0337         outTree->Branch("tpcseed_charge", &tpcseed_charge);
0338         outTree->Branch("tpcseed_nTpc", &tpcseed_nTpc);
0339         outTree->Branch("tpcseed_nMms", &tpcseed_nMms);
0340         outTree->Branch("tpcseed_dedx", &tpcseed_dedx);
0341         outTree->Branch("tpcseed_clusterKeys", &tpcseed_clusterKeys);
0342         outTree->Branch("tpcseed_cluster_layer", &tpcseed_cluster_layer);
0343         outTree->Branch("tpcseed_cluster_globalX", &tpcseed_cluster_globalX);
0344         outTree->Branch("tpcseed_cluster_globalY", &tpcseed_cluster_globalY);
0345         outTree->Branch("tpcseed_cluster_globalZ", &tpcseed_cluster_globalZ);
0346         outTree->Branch("tpcseed_cluster_phi", &tpcseed_cluster_phi);
0347         outTree->Branch("tpcseed_cluster_eta", &tpcseed_cluster_eta);
0348         outTree->Branch("tpcseed_cluster_r", &tpcseed_cluster_r);
0349     }
0350 
0351     // (intt) cluster information
0352     outTree->Branch("clusterKey", &clusterKey);
0353     outTree->Branch("cluster_layer", &cluster_layer);
0354     outTree->Branch("cluster_chip", &cluster_chip);   // for mvtx chip id
0355     outTree->Branch("cluster_stave", &cluster_stave); // for mvtx stave id
0356     outTree->Branch("cluster_globalX", &cluster_globalX);
0357     outTree->Branch("cluster_globalY", &cluster_globalY);
0358     outTree->Branch("cluster_globalZ", &cluster_globalZ);
0359     outTree->Branch("cluster_phi", &cluster_phi);
0360     outTree->Branch("cluster_eta", &cluster_eta);
0361     outTree->Branch("cluster_r", &cluster_r);
0362     outTree->Branch("cluster_phiSize", &cluster_phiSize);
0363     outTree->Branch("cluster_zSize", &cluster_zSize);
0364     outTree->Branch("cluster_adc", &cluster_adc);
0365     outTree->Branch("cluster_timeBucketID", &cluster_timeBucketID);
0366     outTree->Branch("cluster_crossing", &cluster_crossing);
0367     outTree->Branch("cluster_ladderZId", &cluster_ladderZId);     // for intt ladder z id
0368     outTree->Branch("cluster_ladderPhiId", &cluster_ladderPhiId); // for intt ladder phi id
0369     outTree->Branch("cluster_LocalX", &cluster_LocalX);
0370     outTree->Branch("cluster_LocalY", &cluster_LocalY);
0371     // cluster truth matching by max_truth_particle_by_energy
0372     outTree->Branch("cluster_matchedG4P_trackID", &cluster_matchedG4P_trackID);
0373     outTree->Branch("cluster_matchedG4P_PID", &cluster_matchedG4P_PID);
0374     outTree->Branch("cluster_matchedG4P_E", &cluster_matchedG4P_E);
0375     outTree->Branch("cluster_matchedG4P_pT", &cluster_matchedG4P_pT);
0376     outTree->Branch("cluster_matchedG4P_eta", &cluster_matchedG4P_eta);
0377     outTree->Branch("cluster_matchedG4P_phi", &cluster_matchedG4P_phi);
0378 
0379     outTree->Branch("mvtx_seedcluster_key", &mvtx_seedcluster_key);
0380     outTree->Branch("mvtx_seedcluster_layer", &mvtx_seedcluster_layer);
0381     outTree->Branch("mvtx_seedcluster_chip", &mvtx_seedcluster_chip);
0382     outTree->Branch("mvtx_seedcluster_stave", &mvtx_seedcluster_stave);
0383     outTree->Branch("mvtx_seedcluster_globalX", &mvtx_seedcluster_globalX);
0384     outTree->Branch("mvtx_seedcluster_globalY", &mvtx_seedcluster_globalY);
0385     outTree->Branch("mvtx_seedcluster_globalZ", &mvtx_seedcluster_globalZ);
0386     outTree->Branch("mvtx_seedcluster_phi", &mvtx_seedcluster_phi);
0387     outTree->Branch("mvtx_seedcluster_eta", &mvtx_seedcluster_eta);
0388     outTree->Branch("mvtx_seedcluster_r", &mvtx_seedcluster_r);
0389     outTree->Branch("mvtx_seedcluster_phiSize", &mvtx_seedcluster_phiSize);
0390     outTree->Branch("mvtx_seedcluster_zSize", &mvtx_seedcluster_zSize);
0391     outTree->Branch("mvtx_seedcluster_strobeID", &mvtx_seedcluster_strobeID);
0392     outTree->Branch("mvtx_seedcluster_matchedcrossing", &mvtx_seedcluster_matchedcrossing);
0393     outTree->Branch("mvtx_seedcluster_hitX", &mvtx_seedcluster_hitX);
0394     outTree->Branch("mvtx_seedcluster_hitY", &mvtx_seedcluster_hitY);
0395     outTree->Branch("mvtx_seedcluster_hitZ", &mvtx_seedcluster_hitZ);
0396     outTree->Branch("mvtx_seedcluster_hitrow", &mvtx_seedcluster_hitrow);
0397     outTree->Branch("mvtx_seedcluster_hitcol", &mvtx_seedcluster_hitcol);
0398 
0399     if (isSimulation)
0400     {
0401         outTree->Branch("ncoll", &ncoll_);
0402         outTree->Branch("npart", &npart_);
0403         outTree->Branch("N_HepMCGenEvent", &N_HepMCGenEvent);
0404         outTree->Branch("HepMCGenEvent_processID", &HepMCGenEvent_processID);
0405         outTree->Branch("HepMCGenEvent_embeddingID", &HepMCGenEvent_embeddingID);
0406         outTree->Branch("HepMCGenEvent_crossing", &HepMCGenEvent_crossing);
0407         outTree->Branch("nTruthVertex", &nTruthVertex);
0408         outTree->Branch("TruthVertex_isEmbeded", &TruthVertex_isEmbeded);
0409         outTree->Branch("TruthVertexX", &TruthVertexX);
0410         outTree->Branch("TruthVertexY", &TruthVertexY);
0411         outTree->Branch("TruthVertexZ", &TruthVertexZ);
0412         outTree->Branch("TruthVertexT", &TruthVertexT);
0413         outTree->Branch("TruthVertex_crossing", &TruthVertex_crossing);
0414 
0415         outTree->Branch("silseed_ngmvtx", &silseed_ngmvtx);
0416         outTree->Branch("silseed_ngintt", &silseed_ngintt);
0417         outTree->Branch("silseed_cluster_gcluster_key", &silseed_cluster_gcluster_key);
0418         outTree->Branch("silseed_cluster_gcluster_layer", &silseed_cluster_gcluster_layer);
0419         outTree->Branch("silseed_cluster_gcluster_X", &silseed_cluster_gcluster_X);
0420         outTree->Branch("silseed_cluster_gcluster_Y", &silseed_cluster_gcluster_Y);
0421         outTree->Branch("silseed_cluster_gcluster_Z", &silseed_cluster_gcluster_Z);
0422         outTree->Branch("silseed_cluster_gcluster_r", &silseed_cluster_gcluster_r);
0423         outTree->Branch("silseed_cluster_gcluster_phi", &silseed_cluster_gcluster_phi);
0424         outTree->Branch("silseed_cluster_gcluster_eta", &silseed_cluster_gcluster_eta);
0425         outTree->Branch("silseed_cluster_gcluster_edep", &silseed_cluster_gcluster_edep);
0426         outTree->Branch("silseed_cluster_gcluster_adc", &silseed_cluster_gcluster_adc);
0427         outTree->Branch("silseed_cluster_gcluster_phiSize", &silseed_cluster_gcluster_phiSize);
0428         outTree->Branch("silseed_cluster_gcluster_zSize", &silseed_cluster_gcluster_zSize);
0429         outTree->Branch("hasSvtxPHG4ParticleMap", &hasSvtxPHG4ParticleMap);
0430         outTree->Branch("svtxPHG4ParticleMapProcessed", &svtxPHG4ParticleMapProcessed);
0431         outTree->Branch("silseed_f4a_nMatched", &silseed_f4a_nMatched);
0432         outTree->Branch("silseed_f4a_truthTrackID", &silseed_f4a_truthTrackID);
0433         outTree->Branch("silseed_f4a_truthWeight", &silseed_f4a_truthWeight);
0434         outTree->Branch("silseed_f4a_bestTrackID", &silseed_f4a_bestTrackID);
0435         outTree->Branch("silseed_f4a_bestWeight", &silseed_f4a_bestWeight);
0436         outTree->Branch("silseed_f4a_bestG4P_PID", &silseed_f4a_bestG4P_PID);
0437         outTree->Branch("silseed_f4a_bestG4P_E", &silseed_f4a_bestG4P_E);
0438         outTree->Branch("silseed_f4a_bestG4P_pT", &silseed_f4a_bestG4P_pT);
0439         outTree->Branch("silseed_f4a_bestG4P_eta", &silseed_f4a_bestG4P_eta);
0440         outTree->Branch("silseed_f4a_bestG4P_phi", &silseed_f4a_bestG4P_phi);
0441         outTree->Branch("silseed_f4a_bestG4P_ancestor_trackID", &silseed_f4a_bestG4P_ancestor_trackID);
0442         outTree->Branch("silseed_f4a_bestG4P_ancestor_PID", &silseed_f4a_bestG4P_ancestor_PID);
0443 
0444         outTree->Branch("mvtx_seedcluster_matchedG4P_trackID", &mvtx_seedcluster_matchedG4P_trackID);
0445         outTree->Branch("mvtx_seedcluster_matchedG4P_PID", &mvtx_seedcluster_matchedG4P_PID);
0446         outTree->Branch("mvtx_seedcluster_matchedG4P_E", &mvtx_seedcluster_matchedG4P_E);
0447         outTree->Branch("mvtx_seedcluster_matchedG4P_pT", &mvtx_seedcluster_matchedG4P_pT);
0448         outTree->Branch("mvtx_seedcluster_matchedG4P_eta", &mvtx_seedcluster_matchedG4P_eta);
0449         outTree->Branch("mvtx_seedcluster_matchedG4P_phi", &mvtx_seedcluster_matchedG4P_phi);
0450         outTree->Branch("mvtx_seedcluster_matchedG4P_ancestor_trackID", &mvtx_seedcluster_matchedG4P_ancestor_trackID);
0451         outTree->Branch("mvtx_seedcluster_matchedG4P_ancestor_PID", &mvtx_seedcluster_matchedG4P_ancestor_PID);
0452 
0453         outTree->Branch("N_PrimaryPHG4Ptcl", &N_PrimaryPHG4Ptcl);
0454         outTree->Branch("N_sPHENIXPrimary", &N_sPHENIXPrimary);
0455         outTree->Branch("N_AllPHG4Ptcl", &N_AllPHG4Ptcl);
0456 
0457         outTree->Branch("PrimaryPHG4Ptcl_pT", &PrimaryPHG4Ptcl_pT);
0458         outTree->Branch("PrimaryPHG4Ptcl_eta", &PrimaryPHG4Ptcl_eta);
0459         outTree->Branch("PrimaryPHG4Ptcl_phi", &PrimaryPHG4Ptcl_phi);
0460         outTree->Branch("PrimaryPHG4Ptcl_E", &PrimaryPHG4Ptcl_E);
0461         outTree->Branch("PrimaryPHG4Ptcl_PID", &PrimaryPHG4Ptcl_PID);
0462         outTree->Branch("PrimaryPHG4Ptcl_trackID", &PrimaryPHG4Ptcl_trackID);
0463         outTree->Branch("PrimaryPHG4Ptcl_originTrackID", &PrimaryPHG4Ptcl_originTrackID);
0464         outTree->Branch("PrimaryPHG4Ptcl_originVtxID", &PrimaryPHG4Ptcl_originVtxID);
0465         outTree->Branch("PrimaryPHG4Ptcl_originVtxT", &PrimaryPHG4Ptcl_originVtxT);
0466         outTree->Branch("PrimaryPHG4Ptcl_originCrossing", &PrimaryPHG4Ptcl_originCrossing);
0467         outTree->Branch("PrimaryPHG4Ptcl_originIsEmbeded", &PrimaryPHG4Ptcl_originIsEmbeded);
0468         outTree->Branch("PrimaryPHG4Ptcl_truthevalIsPrimary", &PrimaryPHG4Ptcl_truthevalIsPrimary);
0469         outTree->Branch("PrimaryPHG4Ptcl_embedID", &PrimaryPHG4Ptcl_embedID);
0470         outTree->Branch("PrimaryPHG4Ptcl_ParticleClass", &PrimaryPHG4Ptcl_ParticleClass);
0471         outTree->Branch("PrimaryPHG4Ptcl_isStable", &PrimaryPHG4Ptcl_isStable);
0472         outTree->Branch("PrimaryPHG4Ptcl_charge", &PrimaryPHG4Ptcl_charge);
0473         outTree->Branch("PrimaryPHG4Ptcl_isChargedHadron", &PrimaryPHG4Ptcl_isChargedHadron);
0474         outTree->Branch("PrimaryPHG4Ptcl_ancestor_trackID", &PrimaryPHG4Ptcl_ancestor_trackID);
0475         outTree->Branch("PrimaryPHG4Ptcl_ancestor_PID", &PrimaryPHG4Ptcl_ancestor_PID);
0476         outTree->Branch("PrimaryPHG4Ptcl_truthcluster_X", &PrimaryPHG4Ptcl_truthcluster_X);
0477         outTree->Branch("PrimaryPHG4Ptcl_truthcluster_Y", &PrimaryPHG4Ptcl_truthcluster_Y);
0478         outTree->Branch("PrimaryPHG4Ptcl_truthcluster_Z", &PrimaryPHG4Ptcl_truthcluster_Z);
0479         outTree->Branch("PrimaryPHG4Ptcl_truthcluster_edep", &PrimaryPHG4Ptcl_truthcluster_edep);
0480         outTree->Branch("PrimaryPHG4Ptcl_truthcluster_adc", &PrimaryPHG4Ptcl_truthcluster_adc);
0481         outTree->Branch("PrimaryPHG4Ptcl_truthcluster_r", &PrimaryPHG4Ptcl_truthcluster_r);
0482         outTree->Branch("PrimaryPHG4Ptcl_truthcluster_phi", &PrimaryPHG4Ptcl_truthcluster_phi);
0483         outTree->Branch("PrimaryPHG4Ptcl_truthcluster_eta", &PrimaryPHG4Ptcl_truthcluster_eta);
0484         outTree->Branch("PrimaryPHG4Ptcl_truthcluster_phisize", &PrimaryPHG4Ptcl_truthcluster_phisize);
0485         outTree->Branch("PrimaryPHG4Ptcl_truthcluster_zsize", &PrimaryPHG4Ptcl_truthcluster_zsize);
0486         outTree->Branch("PrimaryPHG4Ptcl_recocluster_globalX", &PrimaryPHG4Ptcl_recocluster_globalX);
0487         outTree->Branch("PrimaryPHG4Ptcl_recocluster_globalY", &PrimaryPHG4Ptcl_recocluster_globalY);
0488         outTree->Branch("PrimaryPHG4Ptcl_recocluster_globalZ", &PrimaryPHG4Ptcl_recocluster_globalZ);
0489         outTree->Branch("PrimaryPHG4Ptcl_recocluster_r", &PrimaryPHG4Ptcl_recocluster_r);
0490         outTree->Branch("PrimaryPHG4Ptcl_recocluster_phi", &PrimaryPHG4Ptcl_recocluster_phi);
0491         outTree->Branch("PrimaryPHG4Ptcl_recocluster_eta", &PrimaryPHG4Ptcl_recocluster_eta);
0492         outTree->Branch("PrimaryPHG4Ptcl_recocluster_phisize", &PrimaryPHG4Ptcl_recocluster_phisize);
0493         outTree->Branch("PrimaryPHG4Ptcl_recocluster_zsize", &PrimaryPHG4Ptcl_recocluster_zsize);
0494         outTree->Branch("PrimaryPHG4Ptcl_recocluster_adc", &PrimaryPHG4Ptcl_recocluster_adc);
0495 
0496         outTree->Branch("sPHENIXPrimary_pT", &sPHENIXPrimary_pT);
0497         outTree->Branch("sPHENIXPrimary_eta", &sPHENIXPrimary_eta);
0498         outTree->Branch("sPHENIXPrimary_phi", &sPHENIXPrimary_phi);
0499         outTree->Branch("sPHENIXPrimary_E", &sPHENIXPrimary_E);
0500         outTree->Branch("sPHENIXPrimary_PID", &sPHENIXPrimary_PID);
0501         outTree->Branch("sPHENIXPrimary_trackID", &sPHENIXPrimary_trackID);
0502         outTree->Branch("sPHENIXPrimary_originTrackID", &sPHENIXPrimary_originTrackID);
0503         outTree->Branch("sPHENIXPrimary_originVtxID", &sPHENIXPrimary_originVtxID);
0504         outTree->Branch("sPHENIXPrimary_originVtxT", &sPHENIXPrimary_originVtxT);
0505         outTree->Branch("sPHENIXPrimary_originCrossing", &sPHENIXPrimary_originCrossing);
0506         outTree->Branch("sPHENIXPrimary_originIsEmbeded", &sPHENIXPrimary_originIsEmbeded);
0507         outTree->Branch("sPHENIXPrimary_truthevalIsPrimary", &sPHENIXPrimary_truthevalIsPrimary);
0508         outTree->Branch("sPHENIXPrimary_embedID", &sPHENIXPrimary_embedID);
0509         outTree->Branch("sPHENIXPrimary_ParticleClass", &sPHENIXPrimary_ParticleClass);
0510         outTree->Branch("sPHENIXPrimary_isStable", &sPHENIXPrimary_isStable);
0511         outTree->Branch("sPHENIXPrimary_charge", &sPHENIXPrimary_charge);
0512         outTree->Branch("sPHENIXPrimary_isChargedHadron", &sPHENIXPrimary_isChargedHadron);
0513         outTree->Branch("sPHENIXPrimary_ancestor_trackID", &sPHENIXPrimary_ancestor_trackID);
0514         outTree->Branch("sPHENIXPrimary_ancestor_PID", &sPHENIXPrimary_ancestor_PID);
0515         outTree->Branch("sPHENIXPrimary_truthcluster_X", &sPHENIXPrimary_truthcluster_X);
0516         outTree->Branch("sPHENIXPrimary_truthcluster_Y", &sPHENIXPrimary_truthcluster_Y);
0517         outTree->Branch("sPHENIXPrimary_truthcluster_Z", &sPHENIXPrimary_truthcluster_Z);
0518         outTree->Branch("sPHENIXPrimary_truthcluster_edep", &sPHENIXPrimary_truthcluster_edep);
0519         outTree->Branch("sPHENIXPrimary_truthcluster_adc", &sPHENIXPrimary_truthcluster_adc);
0520         outTree->Branch("sPHENIXPrimary_truthcluster_r", &sPHENIXPrimary_truthcluster_r);
0521         outTree->Branch("sPHENIXPrimary_truthcluster_phi", &sPHENIXPrimary_truthcluster_phi);
0522         outTree->Branch("sPHENIXPrimary_truthcluster_eta", &sPHENIXPrimary_truthcluster_eta);
0523         outTree->Branch("sPHENIXPrimary_truthcluster_phisize", &sPHENIXPrimary_truthcluster_phisize);
0524         outTree->Branch("sPHENIXPrimary_truthcluster_zsize", &sPHENIXPrimary_truthcluster_zsize);
0525         outTree->Branch("sPHENIXPrimary_recocluster_globalX", &sPHENIXPrimary_recocluster_globalX);
0526         outTree->Branch("sPHENIXPrimary_recocluster_globalY", &sPHENIXPrimary_recocluster_globalY);
0527         outTree->Branch("sPHENIXPrimary_recocluster_globalZ", &sPHENIXPrimary_recocluster_globalZ);
0528         outTree->Branch("sPHENIXPrimary_recocluster_r", &sPHENIXPrimary_recocluster_r);
0529         outTree->Branch("sPHENIXPrimary_recocluster_phi", &sPHENIXPrimary_recocluster_phi);
0530         outTree->Branch("sPHENIXPrimary_recocluster_eta", &sPHENIXPrimary_recocluster_eta);
0531         outTree->Branch("sPHENIXPrimary_recocluster_phisize", &sPHENIXPrimary_recocluster_phisize);
0532         outTree->Branch("sPHENIXPrimary_recocluster_zsize", &sPHENIXPrimary_recocluster_zsize);
0533         outTree->Branch("sPHENIXPrimary_recocluster_adc", &sPHENIXPrimary_recocluster_adc);
0534 
0535         outTree->Branch("AllPHG4Ptcl_pT", &AllPHG4Ptcl_pT);
0536         outTree->Branch("AllPHG4Ptcl_eta", &AllPHG4Ptcl_eta);
0537         outTree->Branch("AllPHG4Ptcl_phi", &AllPHG4Ptcl_phi);
0538         outTree->Branch("AllPHG4Ptcl_E", &AllPHG4Ptcl_E);
0539         outTree->Branch("AllPHG4Ptcl_PID", &AllPHG4Ptcl_PID);
0540         outTree->Branch("AllPHG4Ptcl_trackID", &AllPHG4Ptcl_trackID);
0541         outTree->Branch("AllPHG4Ptcl_originTrackID", &AllPHG4Ptcl_originTrackID);
0542         outTree->Branch("AllPHG4Ptcl_originVtxID", &AllPHG4Ptcl_originVtxID);
0543         outTree->Branch("AllPHG4Ptcl_originVtxT", &AllPHG4Ptcl_originVtxT);
0544         outTree->Branch("AllPHG4Ptcl_originCrossing", &AllPHG4Ptcl_originCrossing);
0545         outTree->Branch("AllPHG4Ptcl_originIsEmbeded", &AllPHG4Ptcl_originIsEmbeded);
0546         outTree->Branch("AllPHG4Ptcl_truthevalIsPrimary", &AllPHG4Ptcl_truthevalIsPrimary);
0547         outTree->Branch("AllPHG4Ptcl_embedID", &AllPHG4Ptcl_embedID);
0548         outTree->Branch("AllPHG4Ptcl_ancestor_trackID", &AllPHG4Ptcl_ancestor_trackID);
0549         outTree->Branch("AllPHG4Ptcl_ancestor_PID", &AllPHG4Ptcl_ancestor_PID);
0550         outTree->Branch("AllPHG4Ptcl_truthcluster_X", &AllPHG4Ptcl_truthcluster_X);
0551         outTree->Branch("AllPHG4Ptcl_truthcluster_Y", &AllPHG4Ptcl_truthcluster_Y);
0552         outTree->Branch("AllPHG4Ptcl_truthcluster_Z", &AllPHG4Ptcl_truthcluster_Z);
0553         outTree->Branch("AllPHG4Ptcl_truthcluster_edep", &AllPHG4Ptcl_truthcluster_edep);
0554         outTree->Branch("AllPHG4Ptcl_truthcluster_adc", &AllPHG4Ptcl_truthcluster_adc);
0555         outTree->Branch("AllPHG4Ptcl_truthcluster_r", &AllPHG4Ptcl_truthcluster_r);
0556         outTree->Branch("AllPHG4Ptcl_truthcluster_phi", &AllPHG4Ptcl_truthcluster_phi);
0557         outTree->Branch("AllPHG4Ptcl_truthcluster_eta", &AllPHG4Ptcl_truthcluster_eta);
0558         outTree->Branch("AllPHG4Ptcl_truthcluster_phisize", &AllPHG4Ptcl_truthcluster_phisize);
0559         outTree->Branch("AllPHG4Ptcl_truthcluster_zsize", &AllPHG4Ptcl_truthcluster_zsize);
0560         outTree->Branch("AllPHG4Ptcl_recocluster_globalX", &AllPHG4Ptcl_recocluster_globalX);
0561         outTree->Branch("AllPHG4Ptcl_recocluster_globalY", &AllPHG4Ptcl_recocluster_globalY);
0562         outTree->Branch("AllPHG4Ptcl_recocluster_globalZ", &AllPHG4Ptcl_recocluster_globalZ);
0563         outTree->Branch("AllPHG4Ptcl_recocluster_r", &AllPHG4Ptcl_recocluster_r);
0564         outTree->Branch("AllPHG4Ptcl_recocluster_phi", &AllPHG4Ptcl_recocluster_phi);
0565         outTree->Branch("AllPHG4Ptcl_recocluster_eta", &AllPHG4Ptcl_recocluster_eta);
0566         outTree->Branch("AllPHG4Ptcl_recocluster_phisize", &AllPHG4Ptcl_recocluster_phisize);
0567         outTree->Branch("AllPHG4Ptcl_recocluster_zsize", &AllPHG4Ptcl_recocluster_zsize);
0568         outTree->Branch("AllPHG4Ptcl_recocluster_adc", &AllPHG4Ptcl_recocluster_adc);
0569     }
0570 
0571     return Fun4AllReturnCodes::EVENT_OK;
0572 }
0573 
0574 //____________________________________________________________________________..
0575 int VertexCompare::InitRun(PHCompositeNode *topNode) { return Fun4AllReturnCodes::EVENT_OK; }
0576 
0577 //____________________________________________________________________________..
0578 int VertexCompare::process_event(PHCompositeNode *topNode)
0579 {
0580     std::cout << "VertexCompare::process_event - Processing event " << counter << std::endl;
0581 
0582     if (isSimulation)
0583     {
0584         auto eventheader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0585         if (!eventheader)
0586         {
0587             std::cout << "VertexCompare::process_event - [WARNING] - cannot find EventHeader node EventHeader" << std::endl;
0588         }
0589         else
0590         {
0591             ncoll_ = eventheader->get_ncoll();
0592             npart_ = eventheader->get_npart();
0593         }
0594 
0595         auto geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0596         if (!geneventmap)
0597         {
0598             std::cout << "VertexCompare::process_event - [WARNING] - cannot find PHHepMCGenEventMap node PHHepMCGenEventMap" << std::endl;
0599         }
0600         else
0601         {
0602             for (PHHepMCGenEventMap::ConstIter iter = geneventmap->begin(); iter != geneventmap->end(); ++iter)
0603             {
0604                 PHHepMCGenEvent *genevent = iter->second;
0605                 if (!genevent || !genevent->getEvent())
0606                 {
0607                     continue;
0608                 }
0609 
0610                 HepMCGenEvent_processID.push_back(genevent->getEvent()->signal_process_id());
0611                 HepMCGenEvent_embeddingID.push_back(genevent->get_embedding_id());
0612                 HepMCGenEvent_crossing.push_back(static_cast<int>(std::lround(genevent->get_collision_vertex().t() / sphenix_constants::time_between_crossings)));
0613             }
0614             N_HepMCGenEvent = static_cast<int>(HepMCGenEvent_processID.size());
0615         }
0616     }
0617 
0618     if (truthOnlyOutput_)
0619     {
0620         m_truth_info = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0621         if (!m_truth_info)
0622         {
0623             std::cout << "VertexCompare::process_event - [ERROR/WARNING] - can't find G4TruthInfoContainer node G4TruthInfo" << std::endl;
0624             return Fun4AllReturnCodes::EVENT_OK;
0625         }
0626 
0627         FillTruthParticleTree();
0628         outTree->Fill();
0629         ++counter;
0630         Cleanup();
0631         return Fun4AllReturnCodes::EVENT_OK;
0632     }
0633 
0634     PHNodeIterator dstiter(topNode);
0635     PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "DST"));
0636 
0637     MbdVertexMap *m_dst_mbdvertexmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0638     SvtxVertexMap *m_dst_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0639 
0640     auto globalvertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0641 
0642     clustermap = findNode::getClass<TrkrClusterContainer>(dstNode, clusterContainerName);
0643     clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(dstNode, clusterHitAssocName);
0644     clustercrossingassoc = findNode::getClass<TrkrClusterCrossingAssoc>(dstNode, "TRKR_CLUSTERCROSSINGASSOC");
0645     geometry = findNode::getClass<ActsGeometry>(topNode, geometryNodeName);
0646     tpcgeom = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0647     silseedmap = findNode::getClass<TrackSeedContainer>(topNode, seedContainerName);
0648     tpcseedmap = findNode::getClass<TrackSeedContainer>(topNode, tpcSeedContainerName);
0649     svtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode, svtxTrackMapName);
0650     gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, gl1NodeName);
0651     m_mbdout = findNode::getClass<MbdOut>(topNode, mbdOutNodeName);
0652     m_mbdpmtcontainer = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0653     minimumbiasinfo = findNode::getClass<MinimumBiasInfo>(topNode, "MinimumBiasInfo");
0654     m_CentInfo = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0655 
0656     if (!m_dst_mbdvertexmap)
0657     {
0658         std::cout << "SiliconSeedAnalyzer::process_event - [ERROR/WARNING] - can't find MbdVertexMap node " << "MbdVertexMap" << std::endl;
0659         // return Fun4AllReturnCodes::ABORTEVENT;
0660     }
0661     if (!m_dst_vertexmap)
0662     {
0663         std::cout << "SiliconSeedAnalyzer::process_event - [ERROR/WARNING] - can't find SvtxVertexMap node " << "SvtxVertexMap" << std::endl;
0664         // return Fun4AllReturnCodes::ABORTEVENT;
0665     }
0666     if (!globalvertexmap)
0667     {
0668         std::cout << "SiliconSeedAnalyzer::process_event - [ERROR/WARNING] - can't find GlobalVertexMap node " << "GlobalVertexMap" << std::endl;
0669         // return Fun4AllReturnCodes::ABORTEVENT;
0670     }
0671     if (!clustermap)
0672     {
0673         std::cout << "SiliconSeedAnalyzer::process_event - [ERROR/WARNING] - can't find cluster map node " << clusterContainerName << std::endl;
0674         // return Fun4AllReturnCodes::ABORTEVENT;
0675     }
0676     if (!clusterhitassoc)
0677     {
0678         std::cout << "SiliconSeedAnalyzer::process_event - [ERROR/WARNING] - can't find cluster hit association node " << clusterHitAssocName << std::endl;
0679         // return Fun4AllReturnCodes::ABORTEVENT;
0680     }
0681     if (!clustercrossingassoc)
0682     {
0683         std::cout << "SiliconSeedAnalyzer::process_event - [ERROR/WARNING] - can't find cluster crossing association node TRKR_CLUSTERCROSSINGASSOC" << std::endl;
0684         // return Fun4AllReturnCodes::ABORTEVENT;
0685     }
0686     if (!geometry)
0687     {
0688         std::cout << "SiliconSeedAnalyzer::process_event - [ERROR/WARNING] - can't find ActsGeometry node " << geometryNodeName << std::endl;
0689         // return Fun4AllReturnCodes::ABORTEVENT;
0690     }
0691     if (writeTrackBranches_ && !tpcgeom)
0692     {
0693         std::cout << "SiliconSeedAnalyzer::process_event - [WARNING] - cannot find TPC geometry node TPCGEOMCONTAINER for track dEdx" << std::endl;
0694     }
0695     if (!silseedmap)
0696     {
0697         std::cout << "SiliconSeedAnalyzer::process_event - [ERROR/WARNING] - can't find silicon seed map node " << seedContainerName << std::endl;
0698         // return Fun4AllReturnCodes::ABORTEVENT;
0699     }
0700     if (writeTpcSeedBranches_ && !tpcseedmap)
0701     {
0702         std::cout << "SiliconSeedAnalyzer::process_event - [WARNING] - can't find TPC seed map node " << tpcSeedContainerName << "; will use TPC seeds from SvtxTrackMap" << std::endl;
0703     }
0704     if (!svtxTrackMap)
0705     {
0706         std::cout << "SiliconSeedAnalyzer::process_event - [ERROR/WARNING] - can't find SvtxTrackMap node " << svtxTrackMapName << std::endl;
0707         // return Fun4AllReturnCodes::ABORTEVENT;
0708     }
0709     if (!gl1PacketInfo)
0710     {
0711         std::cout << "SiliconSeedAnalyzer::process_event - [WARNING] - can't find GL1 node " << gl1NodeName << std::endl;
0712         // return Fun4AllReturnCodes::EVENT_OK;
0713     }
0714     else
0715     {
0716         gl1BunchCrossing = gl1PacketInfo->getBunchNumber();
0717         uint64_t triggervec = gl1PacketInfo->getScaledVector();
0718         gl1bco = gl1PacketInfo->getBCO();
0719         auto lbshift = gl1bco << 24U;
0720         bcotr = lbshift >> 24U;
0721         for (int i = 0; i < 64; ++i)
0722         {
0723             bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0724             if (trig_decision)
0725             {
0726                 firedTriggers.push_back(i);
0727             }
0728             triggervec = (triggervec >> 1U) & 0xffffffffU;
0729         }
0730     }
0731     if (!m_mbdout)
0732     {
0733         std::cout << "SiliconSeedAnalyzer::process_event - [WARNING] - can't find MbdOut node " << mbdOutNodeName << std::endl;
0734         // return Fun4AllReturnCodes::EVENT_OK;
0735     }
0736     else
0737     {
0738         mbd_south_npmt = m_mbdout->get_npmt(kMbdSouthArm);
0739         mbd_north_npmt = m_mbdout->get_npmt(kMbdNorthArm);
0740         mbd_south_charge_sum = m_mbdout->get_q(kMbdSouthArm);
0741         mbd_north_charge_sum = m_mbdout->get_q(kMbdNorthArm);
0742         MBD_charge_sum = mbd_south_charge_sum + mbd_north_charge_sum;
0743     }
0744 
0745     if (!m_mbdpmtcontainer)
0746     {
0747         std::cout << "SiliconSeedAnalyzer::process_event - [WARNING] - cannot find MbdPmtContainer node MbdPmtContainer" << std::endl;
0748     }
0749     else
0750     {
0751         for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ++ipmt)
0752         {
0753             const float charge = m_mbdpmtcontainer->get_pmt(ipmt)->get_q();
0754             if (charge <= kMbdPmtChargeThreshold)
0755             {
0756                 continue;
0757             }
0758 
0759             if (ipmt < kMbdPmtsPerArm)
0760             {
0761                 ++mbd_nhitsoverths_south;
0762             }
0763             else
0764             {
0765                 ++mbd_nhitsoverths_north;
0766             }
0767         }
0768     }
0769 
0770     is_min_bias = (minimumbiasinfo) ? minimumbiasinfo->isAuAuMinimumBias() : false;
0771 
0772     // centrality information
0773     if (!m_CentInfo)
0774     {
0775         std::cout << "SiliconSeedAnalyzer::process_event - [WARNING] - can't find CentralityInfo node " << "CentralityInfo" << std::endl;
0776         centrality_mbd_ = -1.;
0777         // return Fun4AllReturnCodes::EVENT_OK;
0778     }
0779     else
0780     {
0781         if (m_CentInfo->has_centrality_bin(CentralityInfo::PROP::mbd_NS))
0782         {
0783             centrality_mbd_ = m_CentInfo->get_centrality_bin(CentralityInfo::PROP::mbd_NS);
0784         }
0785         else
0786         {
0787             std::cout << "[WARNING/ERROR] No centrality information found in CentralityInfo. Setting centrality_mbd_ to -2. Please check!" << std::endl;
0788             m_CentInfo->identify();
0789             centrality_mbd_ = -2.;
0790         }
0791     }
0792 
0793     // mbdVertex = trackerVertexX = trackerVertexY = trackerVertexZ = trackerVertexChisq = std::numeric_limits<float>::quiet_NaN();
0794     // mbdVertex = std::numeric_limits<float>::quiet_NaN();
0795     // nTracks = n_MBDVertex = n_TRKVertex = std::numeric_limits<unsigned int>::quiet_NaN();
0796     nTracks = n_MBDVertex = n_TRKVertex = 0;
0797 
0798     hasMBD = false;
0799     hasTRK = false;
0800 
0801     // std::cout << "Number of vertices in GlobalVertexMap: " << globalvertexmap->size() << std::endl;
0802     for (GlobalVertexMap::ConstIter iter = globalvertexmap->begin(); iter != globalvertexmap->end(); ++iter)
0803     {
0804         GlobalVertex *gvertex = iter->second;
0805 
0806         if (gvertex->count_vtxs(mbdType) != 0)
0807         {
0808             // std::cout << "Found MBD vertex in GlobalVertexMap." << std::endl;
0809             hasMBD = true;
0810 
0811             auto mbditer = gvertex->find_vertexes(mbdType);
0812             auto mbdvertexvector = mbditer->second;
0813 
0814             n_MBDVertex += mbdvertexvector.size();
0815             for (auto &vertex : mbdvertexvector)
0816             {
0817                 MbdVertex *m_dst_vertex = m_dst_mbdvertexmap->find(vertex->get_id())->second;
0818 
0819                 mbdVertexId.push_back(m_dst_vertex->get_id());
0820                 mbdVertex.push_back(m_dst_vertex->get_z());
0821                 mbdVertexCrossing.push_back(m_dst_vertex->get_beam_crossing());
0822 
0823                 // std::cout << "MBD vertex z: " << m_dst_vertex->get_z() << std::endl;
0824             }
0825         }
0826 
0827         if (gvertex->count_vtxs(trkType) != 0)
0828         {
0829             hasTRK = true;
0830 
0831             auto trkiter = gvertex->find_vertexes(trkType);
0832             auto trkvertexvector = trkiter->second;
0833 
0834             n_TRKVertex += trkvertexvector.size();
0835             for (auto &vertex : trkvertexvector)
0836             {
0837                 SvtxVertex *m_dst_vertex = m_dst_vertexmap->find(vertex->get_id())->second;
0838                 // std::cout << "Tracker vertex z: " << m_dst_vertex->get_z() << ", crossing: " << m_dst_vertex->get_beam_crossing() << std::endl;
0839 
0840                 trackerVertexId.push_back(m_dst_vertex->get_id());
0841                 trackerVertexX.push_back(m_dst_vertex->get_x());
0842                 trackerVertexY.push_back(m_dst_vertex->get_y());
0843                 trackerVertexZ.push_back(m_dst_vertex->get_z());
0844                 trackerVertexChisq.push_back(m_dst_vertex->get_chisq());
0845                 trackerVertexNdof.push_back(m_dst_vertex->get_ndof());
0846                 trackerVertexNTracks.push_back(m_dst_vertex->size_tracks());
0847                 trackerVertexCrossing.push_back(m_dst_vertex->get_beam_crossing());
0848 
0849                 // loop over tracks associated with this vertex and save their track IDs
0850                 std::vector<int> trackIDs;
0851                 trackIDs.clear();
0852                 for (auto trackiter = m_dst_vertex->begin_tracks(); trackiter != m_dst_vertex->end_tracks(); ++trackiter)
0853                 {
0854                     trackIDs.push_back(*trackiter);
0855                 }
0856 
0857                 // print out for debugging
0858                 // std::cout << "Tracker vertex ID " << m_dst_vertex->get_id() << " with crossing " << m_dst_vertex->get_beam_crossing() << " m_dst_vertex->size_tracks() = " << m_dst_vertex->size_tracks() << " trackIDs.size() = " << trackIDs.size() << ": [";
0859                 // for (const auto &trackID : trackIDs)
0860                 // {
0861                 //     std::cout << trackID << " ";
0862                 // }
0863                 // std::cout << "]" << std::endl;
0864 
0865                 trackerVertexTrackIDs.push_back(trackIDs);
0866 
0867                 // if (m_dst_vertex->get_beam_crossing() != 0)
0868                 //     continue;
0869                 // if (m_dst_vertex->size_tracks() > nTracks)
0870                 // {
0871                 //     trackerVertexX = m_dst_vertex->get_x();
0872                 //     trackerVertexY = m_dst_vertex->get_y();
0873                 //     trackerVertexZ = m_dst_vertex->get_z();
0874                 //     trackerVertexChisq = m_dst_vertex->get_chisq();
0875                 //     trackerVertexNdof = m_dst_vertex->get_ndof();
0876                 //     nTracks = m_dst_vertex->size_tracks();
0877                 // }
0878                 // if (nTracks == 0)
0879                 //     hasTRK = false;
0880             }
0881         }
0882     }
0883 
0884     nSvtxVertices = trackerVertexId.size();
0885     // count vertices with valid crossing from trackerVertexCrossing (not SHORT_MAX)
0886     nSvtxVertices_validCrossing = std::count_if(trackerVertexCrossing.begin(), trackerVertexCrossing.end(), [](short crossing) { return crossing != std::numeric_limits<short>::max(); });
0887 
0888     // loop over all vertices in MbdVertexMap and fill all vertices to ntuple
0889     // n_MBDVertex = m_dst_mbdvertexmap->size();
0890     // for (const auto& [key, vertex] : *m_dst_mbdvertexmap)
0891     // {
0892     //     mbdVertexId.push_back(vertex->get_id());
0893     //     mbdVertex.push_back(vertex->get_z());
0894     //     mbdVertexCrossing.push_back(vertex->get_beam_crossing());
0895     // }
0896 
0897     // loop over all vertices in SvtxVertexMap and fill all vertices to ntuple
0898     // nSvtxVertices = m_dst_vertexmap->size();
0899     // for (const auto& [key, vertex] : *m_dst_vertexmap)
0900     // {
0901     //     trackerVertexId.push_back(vertex->get_id());
0902     //     trackerVertexX.push_back(vertex->get_x());
0903     //     trackerVertexY.push_back(vertex->get_y());
0904     //     trackerVertexZ.push_back(vertex->get_z());
0905     //     trackerVertexChisq.push_back(vertex->get_chisq()); //! For vertex from PHSimpleVertexFinder, the chisq and ndof are alway 0, need to check if this is intended
0906     //     trackerVertexNdof.push_back(vertex->get_ndof());
0907     //     trackerVertexNTracks.push_back(vertex->size_tracks());
0908     //     trackerVertexCrossing.push_back(vertex->get_beam_crossing());
0909     // }
0910 
0911     // simulation setup
0912     if (isSimulation)
0913     {
0914         svtxPHG4ParticleMap = findNode::getClass<SvtxPHG4ParticleMap>(topNode, svtxPHG4ParticleMapName);
0915         if (!svtxPHG4ParticleMap)
0916         {
0917             std::cout << "[WARNING/ERROR] VertexCompare::process_event - [ERROR/WARNING] - can't find SvtxPHG4ParticleMap node " << svtxPHG4ParticleMapName << std::endl;
0918         }
0919         hasSvtxPHG4ParticleMap = (svtxPHG4ParticleMap != nullptr);
0920         svtxPHG4ParticleMapProcessed = (svtxPHG4ParticleMap != nullptr && svtxPHG4ParticleMap->processed());
0921         // svtxPHG4ParticleMap->identify();
0922 
0923         m_truth_info = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0924         if (!m_truth_info)
0925         {
0926             std::cout << "[WARNING/ERROR] VertexCompare::process_event - [ERROR/WARNING] - can't find G4TruthInfoContainer node " << "G4TruthInfo" << std::endl;
0927             // return Fun4AllReturnCodes::ABORTEVENT;
0928         }
0929 
0930         if (!svtx_evalstack)
0931         {
0932             svtx_evalstack = new SvtxEvalStack(topNode);
0933             svtx_evalstack->set_strict(false);
0934             // svtx_evalstack->do_caching(false);
0935             clustereval = svtx_evalstack->get_cluster_eval();
0936             hiteval = svtx_evalstack->get_hit_eval();
0937             truth_eval = svtx_evalstack->get_truth_eval();
0938         }
0939 
0940         svtx_evalstack->next_event(topNode);
0941     }
0942 
0943     FillSiliconSeedTree();
0944     if (writeTpcSeedBranches_)
0945     {
0946         FillTpcSeedTree();
0947     }
0948     if (writeTrackBranches_)
0949     {
0950         FillTrackTree();
0951     }
0952 
0953     FillClusterTree();
0954 
0955     if (isSimulation)
0956     {
0957         FillTruthParticleTree();
0958     }
0959 
0960     outTree->Fill();
0961 
0962     ++counter;
0963 
0964     Cleanup();
0965 
0966     return Fun4AllReturnCodes::EVENT_OK;
0967 }
0968 
0969 //____________________________________________________________________________..
0970 void VertexCompare::FillTrackTree()
0971 {
0972     const int kMissingInt = std::numeric_limits<int>::max();
0973     const float kMissingFloat = std::numeric_limits<float>::quiet_NaN();
0974 
0975     if (!svtxTrackMap)
0976     {
0977         return;
0978     }
0979 
0980     float layerThicknesses[4] = {0.0, 0.0, 0.0, 0.0};
0981     bool canCalculateDedx = false;
0982     if (clustermap && geometry && tpcgeom)
0983     {
0984         auto *inner1 = tpcgeom->GetLayerCellGeom(7);
0985         auto *inner2 = tpcgeom->GetLayerCellGeom(8);
0986         auto *middle = tpcgeom->GetLayerCellGeom(27);
0987         auto *outer = tpcgeom->GetLayerCellGeom(50);
0988         if (inner1 && inner2 && middle && outer)
0989         {
0990             layerThicknesses[0] = inner1->get_thickness();
0991             layerThicknesses[1] = inner2->get_thickness();
0992             layerThicknesses[2] = middle->get_thickness();
0993             layerThicknesses[3] = outer->get_thickness();
0994             canCalculateDedx = true;
0995         }
0996     }
0997 
0998     const auto count_seed_hits = [](const TrackSeed *seed, int &nhits, int &nmaps, int &nintt, int &ntpc, int &nmms, int &ntpc1, int &ntpc11, int &ntpc2, int &ntpc3)
0999     {
1000         if (!seed)
1001         {
1002             return;
1003         }
1004 
1005         nhits += seed->size_cluster_keys();
1006         for (auto iter = seed->begin_cluster_keys(); iter != seed->end_cluster_keys(); ++iter)
1007         {
1008             const auto clusterKey = *iter;
1009             const unsigned int layer = TrkrDefs::getLayer(clusterKey);
1010             switch (TrkrDefs::getTrkrId(clusterKey))
1011             {
1012             case TrkrDefs::TrkrId::mvtxId:
1013                 ++nmaps;
1014                 break;
1015             case TrkrDefs::TrkrId::inttId:
1016                 ++nintt;
1017                 break;
1018             case TrkrDefs::TrkrId::tpcId:
1019                 ++ntpc;
1020                 if ((layer - 7U) < 8U)
1021                 {
1022                     ++ntpc11;
1023                 }
1024                 if ((layer - 7U) < 16U)
1025                 {
1026                     ++ntpc1;
1027                 }
1028                 else if ((layer - 7U) < 32U)
1029                 {
1030                     ++ntpc2;
1031                 }
1032                 else if ((layer - 7U) < 48U)
1033                 {
1034                     ++ntpc3;
1035                 }
1036                 break;
1037             case TrkrDefs::TrkrId::micromegasId:
1038                 ++nmms;
1039                 break;
1040             default:
1041                 break;
1042             }
1043         }
1044     };
1045 
1046     const auto fill_seed_clusters = [this](const TrackSeed *seed,
1047                                            std::vector<unsigned int> &layers,
1048                                            std::vector<float> &globalX,
1049                                            std::vector<float> &globalY,
1050                                            std::vector<float> &globalZ,
1051                                            std::vector<float> &phis,
1052                                            std::vector<float> &etas,
1053                                            std::vector<float> &radii)
1054     {
1055         if (!seed || !clustermap || !geometry)
1056         {
1057             return;
1058         }
1059 
1060         for (auto iter = seed->begin_cluster_keys(); iter != seed->end_cluster_keys(); ++iter)
1061         {
1062             const auto clusterKey = *iter;
1063             auto *cluster = clustermap->findCluster(clusterKey);
1064             if (!cluster)
1065             {
1066                 continue;
1067             }
1068 
1069             const auto globalpos = geometry->getGlobalPosition(clusterKey, cluster);
1070             layers.push_back(TrkrDefs::getLayer(clusterKey));
1071             globalX.push_back(globalpos.x());
1072             globalY.push_back(globalpos.y());
1073             globalZ.push_back(globalpos.z());
1074             phis.push_back(std::atan2(globalpos.y(), globalpos.x()));
1075             TVector3 posvec(globalpos.x(), globalpos.y(), globalpos.z());
1076             etas.push_back(posvec.Eta());
1077             radii.push_back((globalpos.y() > 0) ? std::sqrt(globalpos.x() * globalpos.x() + globalpos.y() * globalpos.y()) : -std::sqrt(globalpos.x() * globalpos.x() + globalpos.y() * globalpos.y()));
1078         }
1079     };
1080 
1081     for (auto trackIter = svtxTrackMap->begin(); trackIter != svtxTrackMap->end(); ++trackIter)
1082     {
1083         SvtxTrack *track = trackIter->second;
1084         if (!track)
1085         {
1086             continue;
1087         }
1088 
1089         TrackSeed *tpc_seed = track->get_tpc_seed();
1090         TrackSeed *silicon_seed = track->get_silicon_seed();
1091 
1092         int nhits = 0;
1093         int nmaps = 0;
1094         int nintt = 0;
1095         int ntpc = 0;
1096         int nmms = 0;
1097         int ntpc1 = 0;
1098         int ntpc11 = 0;
1099         int ntpc2 = 0;
1100         int ntpc3 = 0;
1101         count_seed_hits(tpc_seed, nhits, nmaps, nintt, ntpc, nmms, ntpc1, ntpc11, ntpc2, ntpc3);
1102         count_seed_hits(silicon_seed, nhits, nmaps, nintt, ntpc, nmms, ntpc1, ntpc11, ntpc2, ntpc3);
1103 
1104         std::vector<unsigned int> cluster_layers;
1105         std::vector<float> cluster_globalX;
1106         std::vector<float> cluster_globalY;
1107         std::vector<float> cluster_globalZ;
1108         std::vector<float> cluster_phi;
1109         std::vector<float> cluster_eta;
1110         std::vector<float> cluster_r;
1111         fill_seed_clusters(tpc_seed, cluster_layers, cluster_globalX, cluster_globalY, cluster_globalZ, cluster_phi, cluster_eta, cluster_r);
1112         fill_seed_clusters(silicon_seed, cluster_layers, cluster_globalX, cluster_globalY, cluster_globalZ, cluster_phi, cluster_eta, cluster_r);
1113 
1114         const float px = track->get_px();
1115         const float py = track->get_py();
1116         const float pz = track->get_pz();
1117         const TVector3 momentum(px, py, pz);
1118         const float pt = momentum.Pt();
1119         const float eta = momentum.Eta();
1120         const float phi = momentum.Phi();
1121 
1122         const float cvxx = track->get_error(3, 3);
1123         const float cvxy = track->get_error(3, 4);
1124         const float cvxz = track->get_error(3, 5);
1125         const float cvyy = track->get_error(4, 4);
1126         const float cvyz = track->get_error(4, 5);
1127         const float cvzz = track->get_error(5, 5);
1128 
1129         const double pt2 = static_cast<double>(px) * px + static_cast<double>(py) * py;
1130         const double p2 = pt2 + static_cast<double>(pz) * pz;
1131 
1132         float deltapt = kMissingFloat;
1133         if (pt2 > 0.)
1134         {
1135             const double arg = (static_cast<double>(cvxx) * px * px + 2. * static_cast<double>(cvxy) * px * py + static_cast<double>(cvyy) * py * py) / pt2;
1136             if (std::isfinite(arg) && arg >= 0.)
1137             {
1138                 deltapt = std::sqrt(arg);
1139             }
1140         }
1141 
1142         float deltaeta = kMissingFloat;
1143         if (pt2 > 0. && p2 > 0.)
1144         {
1145             const double numerator =
1146                 static_cast<double>(cvzz) * pt2 * pt2 +
1147                 static_cast<double>(pz) *
1148                     (-2. * (static_cast<double>(cvxz) * px + static_cast<double>(cvyz) * py) * pt2 +
1149                      static_cast<double>(cvxx) * px * px * pz +
1150                      static_cast<double>(cvyy) * py * py * pz +
1151                      2. * static_cast<double>(cvxy) * px * py * pz);
1152             const double arg = numerator / (pt2 * pt2 * p2);
1153             if (std::isfinite(arg) && arg >= 0.)
1154             {
1155                 deltaeta = std::sqrt(arg);
1156             }
1157         }
1158 
1159         float deltaphi = kMissingFloat;
1160         if (pt2 > 0.)
1161         {
1162             const double arg =
1163                 (static_cast<double>(cvyy) * px * px - 2. * static_cast<double>(cvxy) * px * py + static_cast<double>(cvxx) * py * py) /
1164                 (pt2 * pt2);
1165             if (std::isfinite(arg) && arg >= 0.)
1166             {
1167                 deltaphi = std::sqrt(arg);
1168             }
1169         }
1170 
1171         const int vertexID = track->get_vertex_id();
1172         float vx = kMissingFloat;
1173         float vy = kMissingFloat;
1174         float vz = kMissingFloat;
1175         auto vertexIter = std::find(trackerVertexId.begin(), trackerVertexId.end(), vertexID);
1176         if (vertexIter != trackerVertexId.end())
1177         {
1178             const auto vertexIndex = std::distance(trackerVertexId.begin(), vertexIter);
1179             vx = trackerVertexX[vertexIndex];
1180             vy = trackerVertexY[vertexIndex];
1181             vz = trackerVertexZ[vertexIndex];
1182         }
1183 
1184         float dedx = kMissingFloat;
1185         if (canCalculateDedx && tpc_seed)
1186         {
1187             dedx = TrackAnalysisUtils::calc_dedx(tpc_seed, clustermap, geometry, layerThicknesses);
1188         }
1189 
1190         float hlxpt = kMissingFloat;
1191         float hlxeta = kMissingFloat;
1192         float hlxphi = kMissingFloat;
1193         float hlxX0 = kMissingFloat;
1194         float hlxY0 = kMissingFloat;
1195         float hlxZ0 = kMissingFloat;
1196         int hlxcharge = 0;
1197         if (tpc_seed)
1198         {
1199             hlxpt = tpc_seed->get_pt();
1200             hlxeta = tpc_seed->get_eta();
1201             hlxphi = tpc_seed->get_phi();
1202             hlxX0 = tpc_seed->get_X0();
1203             hlxY0 = tpc_seed->get_Y0();
1204             hlxZ0 = tpc_seed->get_Z0();
1205             if (std::isfinite(tpc_seed->get_qOverR()) && tpc_seed->get_qOverR() != 0.)
1206             {
1207                 hlxcharge = (tpc_seed->get_qOverR() > 0.) ? 1 : -1;
1208             }
1209         }
1210 
1211         track_deltapt.push_back(deltapt);
1212         track_deltaeta.push_back(deltaeta);
1213         track_deltaphi.push_back(deltaphi);
1214         track_nhits.push_back(nhits);
1215         track_nmaps.push_back(nmaps);
1216         track_nintt.push_back(nintt);
1217         track_ntpc.push_back(ntpc);
1218         track_nmms.push_back(nmms);
1219         track_ntpc1.push_back(ntpc1);
1220         track_ntpc11.push_back(ntpc11);
1221         track_ntpc2.push_back(ntpc2);
1222         track_ntpc3.push_back(ntpc3);
1223         track_pidedx.push_back(kMissingFloat);
1224         track_kdedx.push_back(kMissingFloat);
1225         track_prdedx.push_back(kMissingFloat);
1226         track_vx.push_back(vx);
1227         track_vy.push_back(vy);
1228         track_vz.push_back(vz);
1229         track_dca2d.push_back(track->get_dca2d());
1230         track_dca2dsigma.push_back(track->get_dca2d_error());
1231         track_dca3dxy.push_back(track->get_dca3d_xy());
1232         track_dca3dxysigma.push_back(track->get_dca3d_xy_error());
1233         track_dca3dz.push_back(track->get_dca3d_z());
1234         track_dca3dzsigma.push_back(track->get_dca3d_z_error());
1235         track_hlxpt.push_back(hlxpt);
1236         track_hlxeta.push_back(hlxeta);
1237         track_hlxphi.push_back(hlxphi);
1238         track_hlxX0.push_back(hlxX0);
1239         track_hlxY0.push_back(hlxY0);
1240         track_hlxZ0.push_back(hlxZ0);
1241         track_hlxcharge.push_back(hlxcharge);
1242 
1243         track_id.push_back(track->get_id());
1244         track_x.push_back(track->get_x());
1245         track_y.push_back(track->get_y());
1246         track_z.push_back(track->get_z());
1247         track_px.push_back(px);
1248         track_py.push_back(py);
1249         track_pz.push_back(pz);
1250         track_pt.push_back(pt);
1251         track_eta.push_back(eta);
1252         track_phi.push_back(phi);
1253         track_dedx.push_back(dedx);
1254         track_charge.push_back(track->get_charge());
1255         track_crossing.push_back(track->get_crossing());
1256         track_vertex_id.push_back((track->get_vertex_id() == std::numeric_limits<unsigned int>::max()) ? kMissingInt : static_cast<int>(track->get_vertex_id()));
1257         track_chisq.push_back(track->get_chisq());
1258         track_ndf.push_back(static_cast<int>(track->get_ndf()));
1259         track_quality.push_back(track->get_quality());
1260         track_cluster_layer.push_back(cluster_layers);
1261         track_cluster_globalX.push_back(cluster_globalX);
1262         track_cluster_globalY.push_back(cluster_globalY);
1263         track_cluster_globalZ.push_back(cluster_globalZ);
1264         track_cluster_phi.push_back(cluster_phi);
1265         track_cluster_eta.push_back(cluster_eta);
1266         track_cluster_r.push_back(cluster_r);
1267 
1268         if (!silicon_seed)
1269         {
1270             track_silseed_id.push_back(kMissingInt);
1271             track_silseed_x.push_back(kMissingFloat);
1272             track_silseed_y.push_back(kMissingFloat);
1273             track_silseed_z.push_back(kMissingFloat);
1274             track_silseed_pt.push_back(kMissingFloat);
1275             track_silseed_eta.push_back(kMissingFloat);
1276             track_silseed_phi.push_back(kMissingFloat);
1277             track_silseed_crossing.push_back(kMissingInt);
1278             track_silseed_charge.push_back(kMissingInt);
1279             track_silseed_nMvtx.push_back(0);
1280             track_silseed_nIntt.push_back(0);
1281             track_silseed_clusterKeys.push_back(std::vector<uint64_t>());
1282             continue;
1283         }
1284 
1285         const auto seed_pos = TrackSeedHelper::get_xyz(silicon_seed);
1286         const std::vector<uint64_t> seed_cluskeys(silicon_seed->begin_cluster_keys(), silicon_seed->end_cluster_keys());
1287 
1288         int nMvtx = 0;
1289         int nIntt = 0;
1290         for (const auto cluskey : seed_cluskeys)
1291         {
1292             if (TrkrDefs::getTrkrId(cluskey) == TrkrDefs::mvtxId)
1293             {
1294                 ++nMvtx;
1295             }
1296             else if (TrkrDefs::getTrkrId(cluskey) == TrkrDefs::inttId)
1297             {
1298                 ++nIntt;
1299             }
1300         }
1301 
1302         track_silseed_id.push_back(silseedmap ? static_cast<int>(silseedmap->find(silicon_seed)) : kMissingInt);
1303         track_silseed_x.push_back(seed_pos.x());
1304         track_silseed_y.push_back(seed_pos.y());
1305         track_silseed_z.push_back(seed_pos.z());
1306         track_silseed_pt.push_back(silicon_seed->get_pt());
1307         track_silseed_eta.push_back(silicon_seed->get_eta());
1308         track_silseed_phi.push_back(silicon_seed->get_phi());
1309         track_silseed_crossing.push_back(silicon_seed->get_crossing());
1310         track_silseed_charge.push_back((silicon_seed->get_qOverR() > 0) ? 1 : -1);
1311         track_silseed_nMvtx.push_back(nMvtx);
1312         track_silseed_nIntt.push_back(nIntt);
1313         track_silseed_clusterKeys.push_back(seed_cluskeys);
1314     }
1315 
1316     nRecoTracks = track_id.size();
1317     nTracks = nRecoTracks;
1318 
1319     std::cout << "Total number of tracks in this event(trigger frame): " << nRecoTracks << std::endl;
1320 }
1321 
1322 void VertexCompare::FillSiliconSeedTree()
1323 {
1324     constexpr int kUnassociatedVertexId = -std::numeric_limits<int>::max();
1325 
1326     std::vector<float> hit_x;
1327     std::vector<float> hit_y;
1328     std::vector<float> hit_z;
1329     std::vector<int> hit_rows;
1330     std::vector<int> hit_cols;
1331 
1332     std::vector<int> ancestor_trackIDs;
1333     std::vector<int> ancestor_PIDs;
1334 
1335     std::map<unsigned int, unsigned int> svtxTrackIdBySiliconSeed; // silicon seed ID -> SvtxTrack ID
1336     if (svtxTrackMap && doTruthMatching_)
1337     {
1338         for (auto trackIter = svtxTrackMap->begin(); trackIter != svtxTrackMap->end(); ++trackIter)
1339         {
1340             SvtxTrack *track = trackIter->second;
1341             if (!track || !track->get_silicon_seed())
1342             {
1343                 continue;
1344             }
1345 
1346             const auto silicon_seed_id = silseedmap->find(track->get_silicon_seed());
1347 
1348             // print out the svtx track properties and its associated silicon seed properties for debugging --> they are the same, verified
1349             // std::cout << "SvtxTrack ID: " << track->get_id() << ", associated silicon seed ID: " << silicon_seed_id << ", pt: " << track->get_pt() << ", eta: " << track->get_eta() << ", phi: " << track->get_phi() << std::endl;
1350             // std::cout << "Associated silicon seed properties, seed ID: " << silicon_seed_id << ", pt: " << track->get_silicon_seed()->get_pt() << ", eta: " << track->get_silicon_seed()->get_eta() << ", phi: " << track->get_silicon_seed()->get_phi() << std::endl;
1351 
1352             svtxTrackIdBySiliconSeed.emplace(static_cast<unsigned int>(silicon_seed_id), track->get_id());
1353         }
1354     }
1355 
1356     // helper function to find the vertex index based on the seed crossing and index
1357     // look up trackerVertexCrossing and trackerVertexTrackIDs, the seed crossing should match the trackerVertexCrossing
1358     // then from trackerVertexTrackIDs see which of the vertex is the seed associated with, return the vertex index
1359     // If none of the vertex's trackID list contains the seed track ID, it means the seed is not associated with any vertex
1360     // In that case, return -1 and the seed eta/phi w.r.t vertex will be set to some large or small value to indicate they are not associated with any vertex
1361     const auto find_vertex_index_for_crossing = [this](int seed_id, int crossing) -> int
1362     {
1363         // first get all vertex indices with the same crossing
1364         std::vector<int> candidate_vertex_indices;
1365         for (size_t ivtx = 0; ivtx < trackerVertexCrossing.size(); ++ivtx)
1366         {
1367             if (trackerVertexCrossing[ivtx] == crossing)
1368             {
1369                 candidate_vertex_indices.push_back(static_cast<int>(ivtx));
1370             }
1371         }
1372 
1373         // then check if the seed track ID is in any of the candidate vertex track ID list
1374         for (const auto &vtx_index : candidate_vertex_indices)
1375         {
1376             const auto &trackIDs = trackerVertexTrackIDs[vtx_index];
1377             if (std::find(trackIDs.begin(), trackIDs.end(), seed_id) != trackIDs.end())
1378             {
1379                 return vtx_index;
1380             }
1381         }
1382         return -1;
1383     };
1384 
1385     std::vector<std::pair<TrackSeed *, int>> seeds;
1386     std::map<int, int> first_associated_vertex_index_by_crossing;
1387     for (auto iter = silseedmap->begin(); iter != silseedmap->end(); ++iter)
1388     {
1389         auto *seed = *iter;
1390         if (!seed)
1391         {
1392             continue;
1393         }
1394 
1395         const int seed_id = silseedmap->find(seed);
1396         seeds.emplace_back(seed, seed_id);
1397 
1398         const int crossing = seed->get_crossing();
1399         if (first_associated_vertex_index_by_crossing.find(crossing) != first_associated_vertex_index_by_crossing.end())
1400         {
1401             continue;
1402         }
1403 
1404         const int vtx_index = find_vertex_index_for_crossing(seed_id, crossing);
1405         if (vtx_index >= 0)
1406         {
1407             first_associated_vertex_index_by_crossing.emplace(crossing, vtx_index);
1408         }
1409     }
1410 
1411     // loop over all silicon seeds
1412     // make a map of vector of pair of (seed id, associated vertex id) for each unique crossing for debugging purpose
1413     std::map<int, std::vector<std::pair<int, int>>> crossing_SeedIdVertexId_map;
1414     for (const auto &[seed, seed_id] : seeds)
1415     {
1416         if (!seed)
1417         {
1418             // std::cout << "VertexCompare::FillSiliconSeedTree - [ERROR] - silicon seed pointer is null, skip" << std::endl;
1419             continue;
1420         }
1421         silseed_id.push_back(seed_id);
1422         const auto si_pos = TrackSeedHelper::get_xyz(seed);
1423         silseed_x.push_back(si_pos.x());
1424         silseed_y.push_back(si_pos.y());
1425         silseed_z.push_back(si_pos.z());
1426         silseed_crossing.push_back(seed->get_crossing());
1427         if (seed->get_crossing() != SHRT_MAX)
1428         {
1429             ++nSilSeedsValidCrossing;
1430         }
1431 
1432         silseed_pt.push_back(seed->get_pt());
1433         silseed_eta.push_back(seed->get_eta());
1434         silseed_phi.push_back(seed->get_phi());
1435         {
1436             const int associated_vtx_index = find_vertex_index_for_crossing(seed_id, seed->get_crossing());
1437             int eta_phi_vtx_index = associated_vtx_index;
1438             if (eta_phi_vtx_index < 0)
1439             {
1440                 const auto fallback_it = first_associated_vertex_index_by_crossing.find(seed->get_crossing());
1441                 if (fallback_it != first_associated_vertex_index_by_crossing.end())
1442                 {
1443                     eta_phi_vtx_index = fallback_it->second;
1444                 }
1445             }
1446 
1447             if (associated_vtx_index >= 0)
1448             {
1449                 silseed_assocVtxId.push_back(trackerVertexId[associated_vtx_index]);
1450             }
1451             else
1452             {
1453                 silseed_assocVtxId.push_back(kUnassociatedVertexId);
1454             }
1455 
1456             if (eta_phi_vtx_index >= 0)
1457             {
1458                 TVector3 seed_from_vtx(si_pos.x() - trackerVertexX[eta_phi_vtx_index], si_pos.y() - trackerVertexY[eta_phi_vtx_index], si_pos.z() - trackerVertexZ[eta_phi_vtx_index]);
1459                 silseed_eta_vtx.push_back(seed_from_vtx.Eta());
1460                 silseed_phi_vtx.push_back(seed_from_vtx.Phi());
1461             }
1462             else
1463             {
1464                 // set to some minimum value to indicate they are not associated with any vertex
1465                 silseed_eta_vtx.push_back(-1 * std::numeric_limits<float>::max());
1466                 silseed_phi_vtx.push_back(-1 * std::numeric_limits<float>::max());
1467             }
1468         }
1469         silseed_charge.push_back((seed->get_qOverR() > 0) ? 1 : -1);
1470         // filling crossing seed id map for debugging
1471         crossing_SeedIdVertexId_map[seed->get_crossing()].push_back(std::make_pair(seed_id, silseed_assocVtxId.back()));
1472 
1473         std::vector<uint64_t> seed_cluskeys(seed->begin_cluster_keys(), seed->end_cluster_keys());
1474         silseed_clusterKeys.push_back(seed_cluskeys);
1475         silseed_cluster_layer.push_back(std::vector<unsigned int>());
1476         silseed_cluster_globalX.push_back(std::vector<float>());
1477         silseed_cluster_globalY.push_back(std::vector<float>());
1478         silseed_cluster_globalZ.push_back(std::vector<float>());
1479         silseed_cluster_phi.push_back(std::vector<float>());
1480         silseed_cluster_eta.push_back(std::vector<float>());
1481         silseed_cluster_r.push_back(std::vector<float>());
1482         silseed_cluster_phiSize.push_back(std::vector<int>());
1483         silseed_cluster_zSize.push_back(std::vector<int>());
1484         silseed_cluster_strobeID.push_back(std::vector<int>());
1485         silseed_cluster_timeBucketID.push_back(std::vector<int>());
1486         if (isSimulation)
1487         {
1488             silseed_cluster_gcluster_key.push_back(std::vector<uint64_t>());
1489             silseed_cluster_gcluster_layer.push_back(std::vector<unsigned int>());
1490             silseed_cluster_gcluster_X.push_back(std::vector<float>());
1491             silseed_cluster_gcluster_Y.push_back(std::vector<float>());
1492             silseed_cluster_gcluster_Z.push_back(std::vector<float>());
1493             silseed_cluster_gcluster_r.push_back(std::vector<float>());
1494             silseed_cluster_gcluster_phi.push_back(std::vector<float>());
1495             silseed_cluster_gcluster_eta.push_back(std::vector<float>());
1496             silseed_cluster_gcluster_edep.push_back(std::vector<float>());
1497             silseed_cluster_gcluster_adc.push_back(std::vector<int>());
1498             silseed_cluster_gcluster_phiSize.push_back(std::vector<float>());
1499             silseed_cluster_gcluster_zSize.push_back(std::vector<float>());
1500         }
1501         int nMvtx = 0, nIntt = 0;
1502         int ngMvtx = 0, ngIntt = 0;
1503         for (auto cluskey : seed_cluskeys)
1504         {
1505             auto *cluster = clustermap->findCluster(cluskey);
1506 
1507             if (!cluster)
1508             {
1509                 // std::cout << "SiliconSeedAnalyzer::process_event - [ERROR] - cluster pointer is null, skip" << std::endl;
1510                 continue;
1511             }
1512 
1513             unsigned int layer = TrkrDefs::getLayer(cluskey);
1514             silseed_cluster_layer[silseed_cluster_layer.size() - 1].push_back(layer);
1515 
1516             auto globalpos = geometry->getGlobalPosition(cluskey, cluster);
1517             silseed_cluster_globalX[silseed_cluster_globalX.size() - 1].push_back(globalpos.x());
1518             silseed_cluster_globalY[silseed_cluster_globalY.size() - 1].push_back(globalpos.y());
1519             silseed_cluster_globalZ[silseed_cluster_globalZ.size() - 1].push_back(globalpos.z());
1520             float phi = std::atan2(globalpos.y(), globalpos.x());
1521             silseed_cluster_phi[silseed_cluster_phi.size() - 1].push_back(phi);
1522             float r = (globalpos.y() > 0) ? std::sqrt(globalpos.x() * globalpos.x() + globalpos.y() * globalpos.y()) : -std::sqrt(globalpos.x() * globalpos.x() + globalpos.y() * globalpos.y());
1523             silseed_cluster_r[silseed_cluster_r.size() - 1].push_back(r);
1524             TVector3 posvec(globalpos.x(), globalpos.y(), globalpos.z());
1525             silseed_cluster_eta[silseed_cluster_eta.size() - 1].push_back(posvec.Eta());
1526             float phisize = cluster->getPhiSize();
1527             if (phisize <= 0)
1528             {
1529                 phisize += 256;
1530             }
1531             silseed_cluster_phiSize[silseed_cluster_phiSize.size() - 1].push_back(phisize);
1532             float zsize = cluster->getZSize();
1533             silseed_cluster_zSize[silseed_cluster_zSize.size() - 1].push_back(zsize);
1534             if (TrkrDefs::getTrkrId(cluskey) == TrkrDefs::inttId)
1535             {
1536                 ++nIntt;
1537                 silseed_cluster_strobeID[silseed_cluster_strobeID.size() - 1].push_back(std::numeric_limits<int>::min());
1538                 silseed_cluster_timeBucketID[silseed_cluster_timeBucketID.size() - 1].push_back(InttDefs::getTimeBucketId(cluskey));
1539             }
1540             else if (TrkrDefs::getTrkrId(cluskey) == TrkrDefs::mvtxId)
1541             {
1542                 ++nMvtx;
1543                 silseed_cluster_strobeID[silseed_cluster_strobeID.size() - 1].push_back(MvtxDefs::getStrobeId(cluskey));
1544                 silseed_cluster_timeBucketID[silseed_cluster_timeBucketID.size() - 1].push_back(std::numeric_limits<int>::min());
1545 
1546                 // save mvtx seed cluster info
1547                 mvtx_seedcluster_key.push_back(cluskey);
1548                 mvtx_seedcluster_layer.push_back(layer);
1549                 mvtx_seedcluster_chip.push_back(MvtxDefs::getChipId(cluskey));
1550                 mvtx_seedcluster_stave.push_back(MvtxDefs::getStaveId(cluskey));
1551                 mvtx_seedcluster_globalX.push_back(globalpos.x());
1552                 mvtx_seedcluster_globalY.push_back(globalpos.y());
1553                 mvtx_seedcluster_globalZ.push_back(globalpos.z());
1554                 mvtx_seedcluster_phi.push_back(phi);
1555                 mvtx_seedcluster_eta.push_back(posvec.Eta());
1556                 mvtx_seedcluster_r.push_back(r);
1557                 mvtx_seedcluster_phiSize.push_back(phisize);
1558                 mvtx_seedcluster_zSize.push_back(zsize);
1559                 mvtx_seedcluster_strobeID.push_back(MvtxDefs::getStrobeId(cluskey));
1560                 mvtx_seedcluster_matchedcrossing.push_back(seed->get_crossing());
1561                 // get hitrow and hitcol from cluster hit assoc
1562                 if (clusterhitassoc)
1563                 {
1564                     Clean(hit_x);
1565                     Clean(hit_y);
1566                     Clean(hit_z);
1567                     Clean(hit_rows);
1568                     Clean(hit_cols);
1569 
1570                     TrkrClusterHitAssoc::ConstRange hitrange = clusterhitassoc->getHits(cluskey);
1571                     for (TrkrClusterHitAssoc::ConstIterator hititer = hitrange.first; hititer != hitrange.second; ++hititer)
1572                     {
1573                         TrkrDefs::hitkey hitkey = hititer->second;
1574 
1575                         int hitrow = MvtxDefs::getRow(hitkey);
1576                         int hitcol = MvtxDefs::getCol(hitkey);
1577                         hit_rows.push_back(hitrow);
1578                         hit_cols.push_back(hitcol);
1579 
1580                         MvtxPixelDefs::pixelkey pixelkey = MvtxPixelDefs::gen_pixelkey_from_coors(layer, MvtxDefs::getStaveId(cluskey), MvtxDefs::getChipId(cluskey), hitrow, hitcol);
1581                         uint32_t hitsetkey = MvtxPixelDefs::get_hitsetkey(pixelkey);
1582 
1583                         float localX = std::numeric_limits<float>::quiet_NaN();
1584                         float localZ = std::numeric_limits<float>::quiet_NaN();
1585                         SegmentationAlpide::detectorToLocal(hitrow, hitcol, localX, localZ);
1586                         Acts::Vector2 local(localX * Acts::UnitConstants::cm, localZ * Acts::UnitConstants::cm);
1587 
1588                         const auto &surface = geometry->maps().getSiliconSurface(hitsetkey);
1589                         auto glob = surface->localToGlobal(geometry->geometry().getGeoContext(), local, Acts::Vector3());
1590 
1591                         hit_x.push_back(glob.x() / Acts::UnitConstants::cm);
1592                         hit_y.push_back(glob.y() / Acts::UnitConstants::cm);
1593                         hit_z.push_back(glob.z() / Acts::UnitConstants::cm);
1594                     }
1595 
1596                     mvtx_seedcluster_hitX.push_back(hit_x);
1597                     mvtx_seedcluster_hitY.push_back(hit_y);
1598                     mvtx_seedcluster_hitZ.push_back(hit_z);
1599                     mvtx_seedcluster_hitrow.push_back(hit_rows);
1600                     mvtx_seedcluster_hitcol.push_back(hit_cols);
1601                 }
1602                 else
1603                 {
1604                     mvtx_seedcluster_hitX.push_back(std::vector<float>());
1605                     mvtx_seedcluster_hitY.push_back(std::vector<float>());
1606                     mvtx_seedcluster_hitZ.push_back(std::vector<float>());
1607                     mvtx_seedcluster_hitrow.push_back(std::vector<int>());
1608                     mvtx_seedcluster_hitcol.push_back(std::vector<int>());
1609                 }
1610 
1611                 // if simulation, get matched G4 particle info
1612                 if (isSimulation)
1613                 {
1614                     Clean(ancestor_trackIDs);
1615                     Clean(ancestor_PIDs);
1616 
1617                     PHG4Particle *ptcl = clustereval->max_truth_particle_by_energy(cluskey); // alternatively, can also try max_truth_particle_by_cluster_energy
1618                     if (ptcl)
1619                     {
1620                         mvtx_seedcluster_matchedG4P_trackID.push_back(ptcl->get_track_id());
1621                         mvtx_seedcluster_matchedG4P_PID.push_back(ptcl->get_pid());
1622                         mvtx_seedcluster_matchedG4P_E.push_back(ptcl->get_e());
1623                         ROOT::Math::PtEtaPhiEVector g4p4(ptcl->get_px(), ptcl->get_py(), ptcl->get_pz(), ptcl->get_e());
1624                         mvtx_seedcluster_matchedG4P_pT.push_back(g4p4.Pt());
1625                         mvtx_seedcluster_matchedG4P_eta.push_back(g4p4.Eta());
1626                         mvtx_seedcluster_matchedG4P_phi.push_back(g4p4.Phi());
1627 
1628                         // get ancestor info
1629                         PHG4Particle *ancestor = m_truth_info->GetParticle(ptcl->get_parent_id());
1630                         while (ancestor != nullptr)
1631                         {
1632                             ancestor_trackIDs.push_back(ancestor->get_track_id());
1633                             ancestor_PIDs.push_back(ancestor->get_pid());
1634                             ancestor = m_truth_info->GetParticle(ancestor->get_parent_id());
1635                         }
1636                         mvtx_seedcluster_matchedG4P_ancestor_trackID.push_back(ancestor_trackIDs);
1637                         mvtx_seedcluster_matchedG4P_ancestor_PID.push_back(ancestor_PIDs);
1638                     }
1639                     else
1640                     {
1641                         std::cout << "VertexCompare::FillSiliconSeedTree - [WARNING] - no matched G4 particle found for cluster " << cluskey << std::endl;
1642                         mvtx_seedcluster_matchedG4P_trackID.push_back(std::numeric_limits<int>::max());
1643                         mvtx_seedcluster_matchedG4P_PID.push_back(std::numeric_limits<int>::max());
1644                         mvtx_seedcluster_matchedG4P_E.push_back(std::numeric_limits<float>::quiet_NaN());
1645                         mvtx_seedcluster_matchedG4P_pT.push_back(std::numeric_limits<float>::quiet_NaN());
1646                         mvtx_seedcluster_matchedG4P_eta.push_back(std::numeric_limits<float>::quiet_NaN());
1647                         mvtx_seedcluster_matchedG4P_phi.push_back(std::numeric_limits<float>::quiet_NaN());
1648                         mvtx_seedcluster_matchedG4P_ancestor_trackID.push_back(std::vector<int>());
1649                         mvtx_seedcluster_matchedG4P_ancestor_PID.push_back(std::vector<int>());
1650                     }
1651                 }
1652             }
1653             else
1654             {
1655                 std::cout << "SiliconSeedAnalyzer::process_event - [WARNING] - cluster is neither INTT nor MVTX. Please check." << std::endl;
1656                 silseed_cluster_strobeID[silseed_cluster_strobeID.size() - 1].push_back(std::numeric_limits<int>::max());
1657                 silseed_cluster_timeBucketID[silseed_cluster_timeBucketID.size() - 1].push_back(std::numeric_limits<int>::max());
1658             }
1659 
1660             // Fill branch for simulation
1661             if (isSimulation)
1662             {
1663                 std::pair<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster>> truthclus = clustereval->max_truth_cluster_by_energy(cluskey);
1664                 const auto truth_key = truthclus.first;
1665                 const auto &truth_cluster = truthclus.second;
1666                 if (truth_cluster)
1667                 {
1668                     const float gx = truth_cluster->getX();
1669                     const float gy = truth_cluster->getY();
1670                     const float gz = truth_cluster->getZ();
1671                     TVector3 gpos(gx, gy, gz);
1672 
1673                     silseed_cluster_gcluster_key.back().push_back(truth_key);
1674                     silseed_cluster_gcluster_layer.back().push_back(TrkrDefs::getLayer(truth_key));
1675                     silseed_cluster_gcluster_X.back().push_back(gx);
1676                     silseed_cluster_gcluster_Y.back().push_back(gy);
1677                     silseed_cluster_gcluster_Z.back().push_back(gz);
1678                     silseed_cluster_gcluster_r.back().push_back(std::sqrt(gx * gx + gy * gy));
1679                     silseed_cluster_gcluster_phi.back().push_back(gpos.Phi());
1680                     silseed_cluster_gcluster_eta.back().push_back(gpos.Eta());
1681                     silseed_cluster_gcluster_edep.back().push_back(truth_cluster->getError(0, 0));
1682                     silseed_cluster_gcluster_adc.back().push_back(truth_cluster->getAdc());
1683                     silseed_cluster_gcluster_phiSize.back().push_back(truth_cluster->getSize(1, 1));
1684                     silseed_cluster_gcluster_zSize.back().push_back(truth_cluster->getSize(2, 2));
1685 
1686                     if (TrkrDefs::getTrkrId(cluskey) == TrkrDefs::mvtxId)
1687                     {
1688                         ++ngMvtx;
1689                     }
1690                     else if (TrkrDefs::getTrkrId(cluskey) == TrkrDefs::inttId)
1691                     {
1692                         ++ngIntt;
1693                     }
1694                 }
1695                 else
1696                 {
1697                     silseed_cluster_gcluster_key.back().push_back(0);
1698                     silseed_cluster_gcluster_layer.back().push_back(std::numeric_limits<unsigned int>::max());
1699                     silseed_cluster_gcluster_X.back().push_back(-1 * std::numeric_limits<float>::max());
1700                     silseed_cluster_gcluster_Y.back().push_back(-1 * std::numeric_limits<float>::max());
1701                     silseed_cluster_gcluster_Z.back().push_back(-1 * std::numeric_limits<float>::max());
1702                     silseed_cluster_gcluster_r.back().push_back(-1 * std::numeric_limits<float>::max());
1703                     silseed_cluster_gcluster_phi.back().push_back(-1 * std::numeric_limits<float>::max());
1704                     silseed_cluster_gcluster_eta.back().push_back(-1 * std::numeric_limits<float>::max());
1705                     silseed_cluster_gcluster_edep.back().push_back(-1 * std::numeric_limits<float>::max());
1706                     silseed_cluster_gcluster_adc.back().push_back(-1 * std::numeric_limits<int>::max());
1707                     silseed_cluster_gcluster_phiSize.back().push_back(-1 * std::numeric_limits<float>::max());
1708                     silseed_cluster_gcluster_zSize.back().push_back(-1 * std::numeric_limits<float>::max());
1709                 }
1710             } //
1711         } // end loop over clusters associated with the seed
1712         silseed_nMvtx.push_back(nMvtx);
1713         silseed_nIntt.push_back(nIntt);
1714         if (isSimulation)
1715         {
1716             silseed_ngmvtx.push_back(ngMvtx);
1717             silseed_ngintt.push_back(ngIntt);
1718         }
1719 
1720         // get F4A truth matching information if it is simulation
1721         if (isSimulation && doTruthMatching_)
1722         {
1723             std::vector<int> truth_track_ids;
1724             std::vector<float> truth_weights;
1725             std::vector<int> best_ancestor_track_ids;
1726             std::vector<int> best_ancestor_pids;
1727 
1728             int best_track_id = std::numeric_limits<int>::max();
1729             float best_weight = -1 * std::numeric_limits<float>::max();
1730             int best_pid = std::numeric_limits<int>::max();
1731             float best_e = -1 * std::numeric_limits<float>::max();
1732             float best_pt = -1 * std::numeric_limits<float>::max();
1733             float best_eta = -1 * std::numeric_limits<float>::max();
1734             float best_phi = -1 * std::numeric_limits<float>::max();
1735 
1736             const auto svtx_track_id_iter = svtxTrackIdBySiliconSeed.find(static_cast<unsigned int>(seed_id));
1737             if (svtx_track_id_iter != svtxTrackIdBySiliconSeed.end() && svtxPHG4ParticleMap)
1738             {
1739                 const unsigned int svtx_track_id = svtx_track_id_iter->second;
1740                 // std::cout << "Seed ID " << seed_id << " is associated with SvtxTrack ID " << svtx_track_id << std::endl;
1741                 const auto &truth_set = svtxPHG4ParticleMap->get(svtx_track_id);
1742 
1743                 for (auto weight_iter = truth_set.rbegin(); weight_iter != truth_set.rend(); ++weight_iter)
1744                 {
1745                     const float weight = weight_iter->first;
1746                     const auto &truth_ids = weight_iter->second;
1747                     for (auto truth_id_iter = truth_ids.rbegin(); truth_id_iter != truth_ids.rend(); ++truth_id_iter)
1748                     {
1749                         truth_track_ids.push_back(*truth_id_iter);
1750                         truth_weights.push_back(weight);
1751                     }
1752                 }
1753 
1754                 if (!truth_track_ids.empty())
1755                 {
1756                     best_weight = truth_weights.front();
1757                     best_track_id = truth_track_ids.front();
1758 
1759                     PHG4Particle *best_particle = (m_truth_info) ? m_truth_info->GetParticle(best_track_id) : nullptr;
1760                     if (best_particle)
1761                     {
1762                         best_pid = best_particle->get_pid();
1763                         best_e = best_particle->get_e();
1764 
1765                         ROOT::Math::PxPyPzEVector best_p4(best_particle->get_px(), best_particle->get_py(), best_particle->get_pz(), best_particle->get_e());
1766                         best_pt = best_p4.Pt();
1767                         best_eta = best_p4.Eta();
1768                         best_phi = best_p4.Phi();
1769 
1770                         PHG4Particle *ancestor = m_truth_info->GetParticle(best_particle->get_parent_id());
1771                         while (ancestor != nullptr)
1772                         {
1773                             best_ancestor_track_ids.push_back(ancestor->get_track_id());
1774                             best_ancestor_pids.push_back(ancestor->get_pid());
1775                             ancestor = m_truth_info->GetParticle(ancestor->get_parent_id());
1776                         }
1777                     }
1778                 }
1779                 // else
1780                 // {
1781                 //     std::cout << "VertexCompare::FillSiliconSeedTree - [WARNING] - no truth track matched for seed ID " << seed_id << std::endl;
1782                 // }
1783             }
1784 
1785             silseed_f4a_nMatched.push_back(static_cast<int>(truth_track_ids.size()));
1786             silseed_f4a_truthTrackID.push_back(truth_track_ids);
1787             silseed_f4a_truthWeight.push_back(truth_weights);
1788             silseed_f4a_bestTrackID.push_back(best_track_id);
1789             silseed_f4a_bestWeight.push_back(best_weight);
1790             silseed_f4a_bestG4P_PID.push_back(best_pid);
1791             silseed_f4a_bestG4P_E.push_back(best_e);
1792             silseed_f4a_bestG4P_pT.push_back(best_pt);
1793             silseed_f4a_bestG4P_eta.push_back(best_eta);
1794             silseed_f4a_bestG4P_phi.push_back(best_phi);
1795             silseed_f4a_bestG4P_ancestor_trackID.push_back(best_ancestor_track_ids);
1796             silseed_f4a_bestG4P_ancestor_PID.push_back(best_ancestor_pids);
1797         }
1798 
1799     } // end loop over seeds
1800     nTotalSilSeeds = silseed_id.size();
1801     if (isSimulation && doTruthMatching_)
1802     {
1803         int nSilSeedsValidCrossing_noTruthMatch = 0;
1804         for (size_t i = 0; i < silseed_id.size(); ++i)
1805         {
1806             if (silseed_crossing[i] != SHRT_MAX && silseed_f4a_nMatched[i] < 1)
1807             {
1808                 ++nSilSeedsValidCrossing_noTruthMatch;
1809             }
1810         }
1811 
1812         std::cout << "Total silicon seeds in this event: " << nTotalSilSeeds << " with valid crossing: " << nSilSeedsValidCrossing << " and valid crossing but no truth match: " << nSilSeedsValidCrossing_noTruthMatch << std::endl;
1813     }
1814     else if (isSimulation)
1815     {
1816         std::cout << "Total silicon seeds in this event: " << nTotalSilSeeds << " with valid crossing: " << nSilSeedsValidCrossing << std::endl;
1817     }
1818     else
1819     {
1820         std::cout << "Total silicon seeds in this event: " << nTotalSilSeeds << " with valid crossing: " << nSilSeedsValidCrossing << std::endl;
1821     }
1822 
1823     // now print out the crossing seed id map for debugging
1824     // std::cout << "Crossing to seed ID mapping:" << std::endl;
1825     // for (const auto &[crossing, seed_id_vertex_id_pairs] : crossing_SeedIdVertexId_map)
1826     // {
1827     //     std::cout << "  Crossing " << crossing << ": Seed IDs [";
1828     //     for (size_t i = 0; i < seed_id_vertex_id_pairs.size(); ++i)
1829     //     {
1830     //         std::cout << "(" << seed_id_vertex_id_pairs[i].first << ", " << seed_id_vertex_id_pairs[i].second << ")";
1831     //         if (i != seed_id_vertex_id_pairs.size() - 1)
1832     //         {
1833     //             std::cout << ", ";
1834     //         }
1835     //     }
1836     //     std::cout << "]" << std::endl;
1837     // }
1838 }
1839 
1840 //____________________________________________________________________________..
1841 void VertexCompare::FillTpcSeedTree()
1842 {
1843     const float kMissingFloat = std::numeric_limits<float>::quiet_NaN();
1844 
1845     std::vector<std::pair<TrackSeed *, unsigned int>> seeds;
1846 
1847     if (tpcseedmap)
1848     {
1849         for (auto iter = tpcseedmap->begin(); iter != tpcseedmap->end(); ++iter)
1850         {
1851             auto *seed = *iter;
1852             if (seed)
1853             {
1854                 seeds.emplace_back(seed, static_cast<unsigned int>(tpcseedmap->find(seed)));
1855             }
1856         }
1857     }
1858     else if (svtxTrackMap)
1859     {
1860         std::vector<TrackSeed *> seen_seeds;
1861         for (auto trackIter = svtxTrackMap->begin(); trackIter != svtxTrackMap->end(); ++trackIter)
1862         {
1863             SvtxTrack *track = trackIter->second;
1864             TrackSeed *seed = track ? track->get_tpc_seed() : nullptr;
1865             if (!seed || std::find(seen_seeds.begin(), seen_seeds.end(), seed) != seen_seeds.end())
1866             {
1867                 continue;
1868             }
1869 
1870             seen_seeds.push_back(seed);
1871             seeds.emplace_back(seed, static_cast<unsigned int>(seeds.size()));
1872         }
1873     }
1874 
1875     float layerThicknesses[4] = {0.0, 0.0, 0.0, 0.0};
1876     bool canCalculateDedx = false;
1877     if (clustermap && geometry && tpcgeom)
1878     {
1879         auto *inner1 = tpcgeom->GetLayerCellGeom(7);
1880         auto *inner2 = tpcgeom->GetLayerCellGeom(8);
1881         auto *middle = tpcgeom->GetLayerCellGeom(27);
1882         auto *outer = tpcgeom->GetLayerCellGeom(50);
1883         if (inner1 && inner2 && middle && outer)
1884         {
1885             layerThicknesses[0] = inner1->get_thickness();
1886             layerThicknesses[1] = inner2->get_thickness();
1887             layerThicknesses[2] = middle->get_thickness();
1888             layerThicknesses[3] = outer->get_thickness();
1889             canCalculateDedx = true;
1890         }
1891     }
1892 
1893     for (const auto &[seed, seed_id] : seeds)
1894     {
1895         tpcseed_id.push_back(seed_id);
1896 
1897         const auto seed_pos = TrackSeedHelper::get_xyz(seed);
1898         tpcseed_x.push_back(seed_pos.x());
1899         tpcseed_y.push_back(seed_pos.y());
1900         tpcseed_z.push_back(seed_pos.z());
1901         tpcseed_pt.push_back(seed->get_pt());
1902         tpcseed_eta.push_back(seed->get_eta());
1903         tpcseed_phi.push_back(seed->get_phi());
1904         tpcseed_crossing.push_back(seed->get_crossing());
1905         tpcseed_crossing_estimate.push_back(seed->get_crossing_estimate());
1906         tpcseed_charge.push_back((seed->get_qOverR() > 0) ? 1 : -1);
1907         tpcseed_dedx.push_back(canCalculateDedx ? TrackAnalysisUtils::calc_dedx(seed, clustermap, geometry, layerThicknesses) : kMissingFloat);
1908 
1909         std::vector<uint64_t> seed_cluskeys(seed->begin_cluster_keys(), seed->end_cluster_keys());
1910         tpcseed_clusterKeys.push_back(seed_cluskeys);
1911         tpcseed_cluster_layer.push_back(std::vector<unsigned int>());
1912         tpcseed_cluster_globalX.push_back(std::vector<float>());
1913         tpcseed_cluster_globalY.push_back(std::vector<float>());
1914         tpcseed_cluster_globalZ.push_back(std::vector<float>());
1915         tpcseed_cluster_phi.push_back(std::vector<float>());
1916         tpcseed_cluster_eta.push_back(std::vector<float>());
1917         tpcseed_cluster_r.push_back(std::vector<float>());
1918 
1919         int nTpc = 0;
1920         int nMms = 0;
1921         for (const auto cluskey : seed_cluskeys)
1922         {
1923             if (TrkrDefs::getTrkrId(cluskey) == TrkrDefs::tpcId)
1924             {
1925                 ++nTpc;
1926             }
1927             else if (TrkrDefs::getTrkrId(cluskey) == TrkrDefs::micromegasId)
1928             {
1929                 ++nMms;
1930             }
1931 
1932             auto *cluster = clustermap ? clustermap->findCluster(cluskey) : nullptr;
1933             if (!cluster || !geometry)
1934             {
1935                 continue;
1936             }
1937 
1938             const auto globalpos = geometry->getGlobalPosition(cluskey, cluster);
1939             tpcseed_cluster_layer.back().push_back(TrkrDefs::getLayer(cluskey));
1940             tpcseed_cluster_globalX.back().push_back(globalpos.x());
1941             tpcseed_cluster_globalY.back().push_back(globalpos.y());
1942             tpcseed_cluster_globalZ.back().push_back(globalpos.z());
1943             tpcseed_cluster_phi.back().push_back(std::atan2(globalpos.y(), globalpos.x()));
1944             TVector3 posvec(globalpos.x(), globalpos.y(), globalpos.z());
1945             tpcseed_cluster_eta.back().push_back(posvec.Eta());
1946             tpcseed_cluster_r.back().push_back((globalpos.y() > 0) ? std::sqrt(globalpos.x() * globalpos.x() + globalpos.y() * globalpos.y()) : -std::sqrt(globalpos.x() * globalpos.x() + globalpos.y() * globalpos.y()));
1947         }
1948 
1949         tpcseed_nTpc.push_back(nTpc);
1950         tpcseed_nMms.push_back(nMms);
1951     }
1952 
1953     nTotalTpcSeeds = tpcseed_id.size();
1954     std::cout << "Total TPC seeds in this event: " << nTotalTpcSeeds << std::endl;
1955 }
1956 
1957 //____________________________________________________________________________..
1958 void VertexCompare::FillClusterTree()
1959 {
1960     for (const auto &det : {TrkrDefs::TrkrId::inttId})
1961     {
1962         for (const auto &hitsetkey : clustermap->getHitSetKeys(det))
1963         {
1964             auto range = clustermap->getClusters(hitsetkey);
1965             for (auto iter = range.first; iter != range.second; ++iter)
1966             {
1967                 uint64_t key = iter->first;
1968                 clusterKey.push_back(key);
1969                 unsigned int layer = TrkrDefs::getLayer(key);
1970                 cluster_layer.push_back(layer);
1971                 auto *cluster = iter->second;
1972                 auto globalpos = geometry->getGlobalPosition(key, cluster);
1973                 cluster_globalX.push_back(globalpos.x());
1974                 cluster_globalY.push_back(globalpos.y());
1975                 cluster_globalZ.push_back(globalpos.z());
1976                 cluster_phi.push_back(std::atan2(globalpos.y(), globalpos.x()));
1977                 TVector3 posvec(globalpos.x(), globalpos.y(), globalpos.z());
1978                 cluster_eta.push_back(posvec.Eta());
1979                 cluster_r.push_back((globalpos.y() > 0) ? std::sqrt(globalpos.x() * globalpos.x() + globalpos.y() * globalpos.y()) : -std::sqrt(globalpos.x() * globalpos.x() + globalpos.y() * globalpos.y()));
1980                 int phiSize = cluster->getPhiSize();
1981                 if (phiSize <= 0)
1982                 {
1983                     phiSize += 256;
1984                 }
1985                 cluster_phiSize.push_back(phiSize);
1986                 int zSize = cluster->getZSize();
1987                 cluster_zSize.push_back(zSize);
1988                 cluster_adc.push_back(cluster->getAdc());
1989                 switch (det)
1990                 {
1991                 case TrkrDefs::TrkrId::mvtxId:
1992                     cluster_chip.push_back(MvtxDefs::getChipId(key));
1993                     cluster_stave.push_back(MvtxDefs::getStaveId(key));
1994                     cluster_timeBucketID.push_back(std::numeric_limits<int>::min());
1995                     cluster_LocalX.push_back(cluster->getLocalX());
1996                     cluster_LocalY.push_back(cluster->getLocalY());
1997                     break;
1998                 case TrkrDefs::TrkrId::inttId:
1999                     cluster_ladderZId.push_back(InttDefs::getLadderZId(key));
2000                     cluster_ladderPhiId.push_back(InttDefs::getLadderPhiId(key));
2001                     cluster_timeBucketID.push_back(InttDefs::getTimeBucketId(key));
2002                     {
2003                         int crossing = -std::numeric_limits<int>::max();
2004                         if (clustercrossingassoc)
2005                         {
2006                             auto crossing_range = clustercrossingassoc->getCrossings(key);
2007                             if (crossing_range.first != crossing_range.second)
2008                             {
2009                                 crossing = crossing_range.first->second;
2010                             }
2011                         }
2012                         cluster_crossing.push_back(crossing);
2013                     }
2014                     // std::cout << "INTT cluster key: " << key << ", timeBucketID: " << InttDefs::getTimeBucketId(key) << ", crossing: " << cluster_crossing.back() << std::endl;
2015                     cluster_LocalX.push_back(cluster->getLocalX());
2016                     cluster_LocalY.push_back(cluster->getLocalY());
2017                     break;
2018                 default:
2019                     std::cout << "SiliconSeedAnalyzer::process_event - [WARNING] - cluster is neither INTT nor MVTX. Please check." << std::endl;
2020                     cluster_chip.push_back(std::numeric_limits<int>::max());
2021                     cluster_stave.push_back(std::numeric_limits<int>::max());
2022                     cluster_ladderZId.push_back(std::numeric_limits<int>::max());
2023                     cluster_ladderPhiId.push_back(std::numeric_limits<int>::max());
2024                     cluster_timeBucketID.push_back(std::numeric_limits<int>::max());
2025                     cluster_LocalX.push_back(std::numeric_limits<float>::max());
2026                     cluster_LocalY.push_back(std::numeric_limits<float>::max());
2027                     break;
2028                 }
2029 
2030                 if (isSimulation)
2031                 {
2032                     PHG4Particle *ptcl_maxE = (clustereval) ? clustereval->max_truth_particle_by_energy(key) : nullptr;
2033                     if (ptcl_maxE)
2034                     {
2035                         cluster_matchedG4P_trackID.push_back(ptcl_maxE->get_track_id());
2036                         cluster_matchedG4P_PID.push_back(ptcl_maxE->get_pid());
2037                         cluster_matchedG4P_E.push_back(ptcl_maxE->get_e());
2038 
2039                         ROOT::Math::PxPyPzEVector g4p4(ptcl_maxE->get_px(), ptcl_maxE->get_py(), ptcl_maxE->get_pz(), ptcl_maxE->get_e());
2040                         cluster_matchedG4P_pT.push_back(g4p4.Pt());
2041                         cluster_matchedG4P_eta.push_back(g4p4.Eta());
2042                         cluster_matchedG4P_phi.push_back(g4p4.Phi());
2043                     }
2044                     else
2045                     {
2046                         cluster_matchedG4P_trackID.push_back(std::numeric_limits<int>::max());
2047                         cluster_matchedG4P_PID.push_back(std::numeric_limits<int>::max());
2048                         cluster_matchedG4P_E.push_back(-1 * std::numeric_limits<float>::max());
2049                         cluster_matchedG4P_pT.push_back(-1 * std::numeric_limits<float>::max());
2050                         cluster_matchedG4P_eta.push_back(-1 * std::numeric_limits<float>::max());
2051                         cluster_matchedG4P_phi.push_back(-1 * std::numeric_limits<float>::max());
2052                     }
2053                 }
2054             }
2055         }
2056     }
2057 }
2058 
2059 //____________________________________________________________________________..
2060 void VertexCompare::FillTruthParticleTree()
2061 {
2062     if (!m_truth_info)
2063     {
2064         std::cout << "VertexCompare::FillTruthParticleTree - [WARNING] - missing G4TruthInfo, skip filling truth particle branches" << std::endl;
2065         return;
2066     }
2067 
2068     // truth primary vertices
2069     const auto vtx_range = m_truth_info->GetPrimaryVtxRange();
2070     for (auto iter = vtx_range.first; iter != vtx_range.second; ++iter)
2071     {
2072         const int point_id = iter->first;
2073         PHG4VtxPoint *point = iter->second;
2074         if (!point)
2075         {
2076             continue;
2077         }
2078 
2079         ++nTruthVertex;
2080         // if (m_truth_info->isEmbededVtx(point_id) == 0)
2081         {
2082             TruthVertex_isEmbeded.push_back(m_truth_info->isEmbededVtx(point_id));
2083             TruthVertexX.push_back(point->get_x());
2084             TruthVertexY.push_back(point->get_y());
2085             TruthVertexZ.push_back(point->get_z());
2086             TruthVertexT.push_back(point->get_t());
2087             TruthVertex_crossing.push_back(std::isfinite(point->get_t()) ? static_cast<int>(std::round(point->get_t() / sphenix_constants::time_between_crossings)) : std::numeric_limits<int>::max());
2088         }
2089     }
2090 
2091     auto fill_ancestor_info = [this](PHG4Particle *ptcl, std::vector<int> &ancestor_trackIDs, std::vector<int> &ancestor_PIDs)
2092     {
2093         PHG4Particle *ancestor = m_truth_info->GetParticle(ptcl->get_parent_id());
2094         while (ancestor != nullptr)
2095         {
2096             ancestor_trackIDs.push_back(ancestor->get_track_id());
2097             ancestor_PIDs.push_back(ancestor->get_pid());
2098             ancestor = m_truth_info->GetParticle(ancestor->get_parent_id());
2099         }
2100     };
2101 
2102     auto fill_particle_kinematics = [](PHG4Particle *ptcl, std::vector<float> &out_pT, std::vector<float> &out_eta, std::vector<float> &out_phi, std::vector<float> &out_E, std::vector<int> &out_PID, std::vector<int> &out_trackID)
2103     {
2104         ROOT::Math::PxPyPzEVector p4(ptcl->get_px(), ptcl->get_py(), ptcl->get_pz(), ptcl->get_e());
2105         out_pT.push_back(p4.Pt());
2106         out_eta.push_back(p4.Eta());
2107         out_phi.push_back(p4.Phi());
2108         out_E.push_back(ptcl->get_e());
2109         out_PID.push_back(ptcl->get_pid());
2110         out_trackID.push_back(ptcl->get_track_id());
2111     };
2112 
2113     auto fill_origin_info = [this](PHG4Particle *ptcl, std::vector<int> &out_originTrackID, std::vector<int> &out_originVtxID, std::vector<float> &out_originVtxT, std::vector<int> &out_originCrossing, std::vector<int> &out_originIsEmbeded)
2114     {
2115         const int kMissingInt = std::numeric_limits<int>::max();
2116         const float kMissingFloat = std::numeric_limits<float>::quiet_NaN();
2117 
2118         PHG4Particle *origin = ptcl;
2119         PHG4Particle *parent = origin ? m_truth_info->GetParticle(origin->get_parent_id()) : nullptr;
2120         while (parent)
2121         {
2122             origin = parent;
2123             parent = m_truth_info->GetParticle(origin->get_parent_id());
2124         }
2125 
2126         if (!origin)
2127         {
2128             out_originTrackID.push_back(kMissingInt);
2129             out_originVtxID.push_back(kMissingInt);
2130             out_originVtxT.push_back(kMissingFloat);
2131             out_originCrossing.push_back(kMissingInt);
2132             out_originIsEmbeded.push_back(kMissingInt);
2133             return;
2134         }
2135 
2136         const int origin_track_id = origin->get_track_id();
2137         const int origin_vtx_id = origin->get_vtx_id();
2138         PHG4VtxPoint *origin_vtx = m_truth_info->GetVtx(origin_vtx_id);
2139         const float origin_vtx_t = origin_vtx ? origin_vtx->get_t() : kMissingFloat;
2140         const int origin_crossing = std::isfinite(origin_vtx_t) ? static_cast<int>(std::round(origin_vtx_t / sphenix_constants::time_between_crossings)) : kMissingInt;
2141 
2142         out_originTrackID.push_back(origin_track_id);
2143         out_originVtxID.push_back(origin_vtx_id);
2144         out_originVtxT.push_back(origin_vtx_t);
2145         out_originCrossing.push_back(origin_crossing);
2146         out_originIsEmbeded.push_back(m_truth_info->isEmbeded(origin_track_id));
2147     };
2148 
2149     auto fill_primary_particle_info = [](PHG4Particle *ptcl, std::vector<TString> &out_particleClass, std::vector<bool> &out_isStable, std::vector<double> &out_charge, std::vector<bool> &out_isChargedHadron)
2150     {
2151         TString particleClass = "Unknown";
2152         bool isStable = false;
2153         double charge = 0;
2154 
2155         TParticlePDG *pdg = TDatabasePDG::Instance()->GetParticle(ptcl->get_pid());
2156         if (pdg)
2157         {
2158             particleClass = TString(pdg->ParticleClass());
2159             isStable = (pdg->Stable() == 1);
2160             charge = pdg->Charge();
2161         }
2162 
2163         bool isHadron = (particleClass.Contains("Baryon") || particleClass.Contains("Meson"));
2164         bool isChargedHadron = (isStable && (charge != 0) && isHadron);
2165 
2166         out_particleClass.push_back(particleClass);
2167         out_isStable.push_back(isStable);
2168         out_charge.push_back(charge);
2169         out_isChargedHadron.push_back(isChargedHadron);
2170     };
2171 
2172     auto fill_truthcluster_match_info = [this](PHG4Particle *ptcl,                                        //
2173                                                std::vector<std::vector<float>> &out_truthcluster_X,       //
2174                                                std::vector<std::vector<float>> &out_truthcluster_Y,       //
2175                                                std::vector<std::vector<float>> &out_truthcluster_Z,       //
2176                                                std::vector<std::vector<float>> &out_truthcluster_edep,    //
2177                                                std::vector<std::vector<float>> &out_truthcluster_adc,     //
2178                                                std::vector<std::vector<float>> &out_truthcluster_r,       //
2179                                                std::vector<std::vector<float>> &out_truthcluster_phi,     //
2180                                                std::vector<std::vector<float>> &out_truthcluster_eta,     //
2181                                                std::vector<std::vector<float>> &out_truthcluster_phisize, //
2182                                                std::vector<std::vector<float>> &out_truthcluster_zsize,   //
2183                                                std::vector<std::vector<float>> &out_recocluster_globalX,  //
2184                                                std::vector<std::vector<float>> &out_recocluster_globalY,  //
2185                                                std::vector<std::vector<float>> &out_recocluster_globalZ,  //
2186                                                std::vector<std::vector<float>> &out_recocluster_r,        //
2187                                                std::vector<std::vector<float>> &out_recocluster_phi,      //
2188                                                std::vector<std::vector<float>> &out_recocluster_eta,      //
2189                                                std::vector<std::vector<float>> &out_recocluster_phisize,  //
2190                                                std::vector<std::vector<float>> &out_recocluster_zsize,    //
2191                                                std::vector<std::vector<float>> &out_recocluster_adc       //
2192                                         )
2193     {
2194         std::vector<float> truthcluster_X;
2195         std::vector<float> truthcluster_Y;
2196         std::vector<float> truthcluster_Z;
2197         std::vector<float> truthcluster_edep;
2198         std::vector<float> truthcluster_adc;
2199         std::vector<float> truthcluster_r;
2200         std::vector<float> truthcluster_phi;
2201         std::vector<float> truthcluster_eta;
2202         std::vector<float> truthcluster_phisize;
2203         std::vector<float> truthcluster_zsize;
2204         std::vector<float> recocluster_globalX;
2205         std::vector<float> recocluster_globalY;
2206         std::vector<float> recocluster_globalZ;
2207         std::vector<float> recocluster_r;
2208         std::vector<float> recocluster_phi;
2209         std::vector<float> recocluster_eta;
2210         std::vector<float> recocluster_phisize;
2211         std::vector<float> recocluster_zsize;
2212         std::vector<float> recocluster_adc;
2213 
2214         if (truth_eval && clustereval)
2215         {
2216             const auto truth_clusters = truth_eval->all_truth_clusters(ptcl);
2217 
2218             for (const auto &[ckey, gclus] : truth_clusters)
2219             {
2220                 if (!gclus)
2221                 {
2222                     continue;
2223                 }
2224 
2225                 // keep tracker truth clusters that have analysis branches below
2226                 if (TrkrDefs::getTrkrId(ckey) != TrkrDefs::mvtxId && TrkrDefs::getTrkrId(ckey) != TrkrDefs::inttId && TrkrDefs::getTrkrId(ckey) != TrkrDefs::tpcId && TrkrDefs::getTrkrId(ckey) != TrkrDefs::micromegasId)
2227                 {
2228                     continue;
2229                 }
2230 
2231                 const float gx = gclus->getX();
2232                 const float gy = gclus->getY();
2233                 const float gz = gclus->getZ();
2234                 const float gedep = gclus->getError(0, 0);
2235                 const float gadc = static_cast<float>(gclus->getAdc());
2236 
2237                 TVector3 gpos(gx, gy, gz);
2238                 truthcluster_X.push_back(gx);
2239                 truthcluster_Y.push_back(gy);
2240                 truthcluster_Z.push_back(gz);
2241                 truthcluster_edep.push_back(gedep);
2242                 truthcluster_adc.push_back(gadc);
2243                 truthcluster_r.push_back(std::sqrt(gx * gx + gy * gy));
2244                 truthcluster_phi.push_back(gpos.Phi());
2245                 truthcluster_eta.push_back(gpos.Eta());
2246                 truthcluster_phisize.push_back(gclus->getSize(1, 1));
2247                 truthcluster_zsize.push_back(gclus->getSize(2, 2));
2248 
2249                 float x = std::numeric_limits<float>::quiet_NaN();
2250                 float y = std::numeric_limits<float>::quiet_NaN();
2251                 float z = std::numeric_limits<float>::quiet_NaN();
2252                 float r = std::numeric_limits<float>::quiet_NaN();
2253                 float phi = std::numeric_limits<float>::quiet_NaN();
2254                 float eta = std::numeric_limits<float>::quiet_NaN();
2255                 float phisize = std::numeric_limits<float>::quiet_NaN();
2256                 float zsize = std::numeric_limits<float>::quiet_NaN();
2257                 float adc = std::numeric_limits<float>::quiet_NaN();
2258 
2259                 const auto [reco_ckey, reco_cluster] = clustereval->reco_cluster_from_truth_cluster(ckey, gclus);
2260                 if (reco_cluster)
2261                 {
2262                     const auto global = geometry->getGlobalPosition(reco_ckey, reco_cluster);
2263                     x = global.x();
2264                     y = global.y();
2265                     z = global.z();
2266 
2267                     TVector3 reco_pos(x, y, z);
2268                     r = std::sqrt(x * x + y * y);
2269                     phi = reco_pos.Phi();
2270                     eta = reco_pos.Eta();
2271                     phisize = reco_cluster->getPhiSize();
2272                     zsize = reco_cluster->getZSize();
2273                     adc = static_cast<float>(reco_cluster->getAdc());
2274                 }
2275 
2276                 recocluster_globalX.push_back(x);
2277                 recocluster_globalY.push_back(y);
2278                 recocluster_globalZ.push_back(z);
2279                 recocluster_r.push_back(r);
2280                 recocluster_phi.push_back(phi);
2281                 recocluster_eta.push_back(eta);
2282                 recocluster_phisize.push_back(phisize);
2283                 recocluster_zsize.push_back(zsize);
2284                 recocluster_adc.push_back(adc);
2285             }
2286         }
2287 
2288         out_truthcluster_X.push_back(truthcluster_X);
2289         out_truthcluster_Y.push_back(truthcluster_Y);
2290         out_truthcluster_Z.push_back(truthcluster_Z);
2291         out_truthcluster_edep.push_back(truthcluster_edep);
2292         out_truthcluster_adc.push_back(truthcluster_adc);
2293         out_truthcluster_r.push_back(truthcluster_r);
2294         out_truthcluster_phi.push_back(truthcluster_phi);
2295         out_truthcluster_eta.push_back(truthcluster_eta);
2296         out_truthcluster_phisize.push_back(truthcluster_phisize);
2297         out_truthcluster_zsize.push_back(truthcluster_zsize);
2298         out_recocluster_globalX.push_back(recocluster_globalX);
2299         out_recocluster_globalY.push_back(recocluster_globalY);
2300         out_recocluster_globalZ.push_back(recocluster_globalZ);
2301         out_recocluster_r.push_back(recocluster_r);
2302         out_recocluster_phi.push_back(recocluster_phi);
2303         out_recocluster_eta.push_back(recocluster_eta);
2304         out_recocluster_phisize.push_back(recocluster_phisize);
2305         out_recocluster_zsize.push_back(recocluster_zsize);
2306         out_recocluster_adc.push_back(recocluster_adc);
2307     };
2308 
2309     std::vector<int> ancestor_trackIDs;
2310     std::vector<int> ancestor_PIDs;
2311 
2312     // primary PHG4 particles
2313     const auto primary_range = m_truth_info->GetPrimaryParticleRange();
2314     for (auto iter = primary_range.first; iter != primary_range.second; ++iter)
2315     {
2316         PHG4Particle *ptcl = iter->second;
2317         if (!ptcl)
2318         {
2319             continue;
2320         }
2321 
2322         Clean(ancestor_trackIDs);
2323         Clean(ancestor_PIDs);
2324         fill_particle_kinematics(ptcl, PrimaryPHG4Ptcl_pT, PrimaryPHG4Ptcl_eta, PrimaryPHG4Ptcl_phi, PrimaryPHG4Ptcl_E, PrimaryPHG4Ptcl_PID, PrimaryPHG4Ptcl_trackID);
2325         fill_origin_info(ptcl, PrimaryPHG4Ptcl_originTrackID, PrimaryPHG4Ptcl_originVtxID, PrimaryPHG4Ptcl_originVtxT, PrimaryPHG4Ptcl_originCrossing, PrimaryPHG4Ptcl_originIsEmbeded);
2326         PrimaryPHG4Ptcl_truthevalIsPrimary.push_back(truth_eval->is_primary(ptcl));
2327         PrimaryPHG4Ptcl_embedID.push_back(truth_eval->get_embed(ptcl));
2328         fill_primary_particle_info(ptcl, PrimaryPHG4Ptcl_ParticleClass, PrimaryPHG4Ptcl_isStable, PrimaryPHG4Ptcl_charge, PrimaryPHG4Ptcl_isChargedHadron);
2329         fill_ancestor_info(ptcl, ancestor_trackIDs, ancestor_PIDs);
2330         PrimaryPHG4Ptcl_ancestor_trackID.push_back(ancestor_trackIDs);
2331         PrimaryPHG4Ptcl_ancestor_PID.push_back(ancestor_PIDs);
2332         fill_truthcluster_match_info(ptcl, PrimaryPHG4Ptcl_truthcluster_X, PrimaryPHG4Ptcl_truthcluster_Y, PrimaryPHG4Ptcl_truthcluster_Z, PrimaryPHG4Ptcl_truthcluster_edep, PrimaryPHG4Ptcl_truthcluster_adc, PrimaryPHG4Ptcl_truthcluster_r, PrimaryPHG4Ptcl_truthcluster_phi, PrimaryPHG4Ptcl_truthcluster_eta, PrimaryPHG4Ptcl_truthcluster_phisize, PrimaryPHG4Ptcl_truthcluster_zsize,
2333                                      PrimaryPHG4Ptcl_recocluster_globalX, PrimaryPHG4Ptcl_recocluster_globalY, PrimaryPHG4Ptcl_recocluster_globalZ, PrimaryPHG4Ptcl_recocluster_r, PrimaryPHG4Ptcl_recocluster_phi, PrimaryPHG4Ptcl_recocluster_eta, PrimaryPHG4Ptcl_recocluster_phisize, PrimaryPHG4Ptcl_recocluster_zsize, PrimaryPHG4Ptcl_recocluster_adc);
2334     }
2335     // clean ancestor info vectors before reusing for sPHENIX primary particles
2336     Clean(ancestor_trackIDs);
2337     Clean(ancestor_PIDs);
2338 
2339     const auto sPHENIXprimary_particle_range = m_truth_info->GetSPHENIXPrimaryParticleRange();
2340     if (VertexCompareVerbosity::fillTruthParticle > 5)
2341     {
2342         std::cout << "Number of sPHENIX primary particles: " << std::distance(sPHENIXprimary_particle_range.first, sPHENIXprimary_particle_range.second) << std::endl;
2343     }
2344     for (auto iter = sPHENIXprimary_particle_range.first; iter != sPHENIXprimary_particle_range.second; ++iter)
2345     {
2346         PHG4Particle *ptcl = iter->second;
2347         if (!ptcl)
2348         {
2349             continue;
2350         }
2351         // get the particle in truth info container with the same track ID
2352         PHG4Particle *real_ptcl = m_truth_info->GetParticle(ptcl->get_track_id());
2353 
2354         Clean(ancestor_trackIDs);
2355         Clean(ancestor_PIDs);
2356         fill_particle_kinematics(real_ptcl, sPHENIXPrimary_pT, sPHENIXPrimary_eta, sPHENIXPrimary_phi, sPHENIXPrimary_E, sPHENIXPrimary_PID, sPHENIXPrimary_trackID);
2357         fill_origin_info(real_ptcl, sPHENIXPrimary_originTrackID, sPHENIXPrimary_originVtxID, sPHENIXPrimary_originVtxT, sPHENIXPrimary_originCrossing, sPHENIXPrimary_originIsEmbeded);
2358         sPHENIXPrimary_truthevalIsPrimary.push_back(truth_eval->is_primary(real_ptcl));
2359         sPHENIXPrimary_embedID.push_back(truth_eval->get_embed(real_ptcl));
2360         fill_primary_particle_info(real_ptcl, sPHENIXPrimary_ParticleClass, sPHENIXPrimary_isStable, sPHENIXPrimary_charge, sPHENIXPrimary_isChargedHadron);
2361         fill_ancestor_info(real_ptcl, ancestor_trackIDs, ancestor_PIDs);
2362         sPHENIXPrimary_ancestor_trackID.push_back(ancestor_trackIDs);
2363         sPHENIXPrimary_ancestor_PID.push_back(ancestor_PIDs);
2364         fill_truthcluster_match_info(real_ptcl, sPHENIXPrimary_truthcluster_X, sPHENIXPrimary_truthcluster_Y, sPHENIXPrimary_truthcluster_Z, sPHENIXPrimary_truthcluster_edep, sPHENIXPrimary_truthcluster_adc, sPHENIXPrimary_truthcluster_r, sPHENIXPrimary_truthcluster_phi, sPHENIXPrimary_truthcluster_eta, sPHENIXPrimary_truthcluster_phisize, sPHENIXPrimary_truthcluster_zsize,
2365                                      sPHENIXPrimary_recocluster_globalX, sPHENIXPrimary_recocluster_globalY, sPHENIXPrimary_recocluster_globalZ, sPHENIXPrimary_recocluster_r, sPHENIXPrimary_recocluster_phi, sPHENIXPrimary_recocluster_eta, sPHENIXPrimary_recocluster_phisize, sPHENIXPrimary_recocluster_zsize, sPHENIXPrimary_recocluster_adc);
2366 
2367         // print out the truth particle info and how many truth and reco clusters are matched for debugging
2368         if (VertexCompareVerbosity::fillTruthParticle > 5)
2369         {
2370             std::cout << "sPHENIX Primary Particle - trackID: " << real_ptcl->get_track_id() << ", PID: " << real_ptcl->get_pid() << ", pT: " << sPHENIXPrimary_pT.back() << ", eta: " << sPHENIXPrimary_eta.back() << ", phi: " << sPHENIXPrimary_phi.back() << ", E: " << sPHENIXPrimary_E.back() << std::endl;
2371             std::cout << "  Matched truth clusters: " << sPHENIXPrimary_truthcluster_X.back().size() << std::endl;
2372             std::cout << "  Matched reco clusters: " << sPHENIXPrimary_recocluster_globalX.back().size() << std::endl;
2373             // flag a particle if it has more reco clusters than truth clusters
2374             if (sPHENIXPrimary_recocluster_globalX.back().size() > sPHENIXPrimary_truthcluster_X.back().size())
2375             {
2376                 std::cout << "  *** Particle has more reco clusters than truth clusters ***" << std::endl;
2377             }
2378         }
2379     }
2380     // again clean it before reusing for all PHG4 particles
2381     Clean(ancestor_trackIDs);
2382     Clean(ancestor_PIDs);
2383 
2384     // all PHG4 particles
2385     const auto all_particle_range = m_truth_info->GetParticleRange();
2386     for (auto iter = all_particle_range.first; iter != all_particle_range.second; ++iter)
2387     {
2388         PHG4Particle *ptcl = iter->second;
2389         if (!ptcl)
2390         {
2391             continue;
2392         }
2393 
2394         Clean(ancestor_trackIDs);
2395         Clean(ancestor_PIDs);
2396         fill_particle_kinematics(ptcl, AllPHG4Ptcl_pT, AllPHG4Ptcl_eta, AllPHG4Ptcl_phi, AllPHG4Ptcl_E, AllPHG4Ptcl_PID, AllPHG4Ptcl_trackID);
2397         fill_origin_info(ptcl, AllPHG4Ptcl_originTrackID, AllPHG4Ptcl_originVtxID, AllPHG4Ptcl_originVtxT, AllPHG4Ptcl_originCrossing, AllPHG4Ptcl_originIsEmbeded);
2398         AllPHG4Ptcl_truthevalIsPrimary.push_back(truth_eval->is_primary(ptcl));
2399         AllPHG4Ptcl_embedID.push_back(truth_eval->get_embed(ptcl));
2400         fill_ancestor_info(ptcl, ancestor_trackIDs, ancestor_PIDs);
2401         AllPHG4Ptcl_ancestor_trackID.push_back(ancestor_trackIDs);
2402         AllPHG4Ptcl_ancestor_PID.push_back(ancestor_PIDs);
2403         fill_truthcluster_match_info(ptcl, AllPHG4Ptcl_truthcluster_X, AllPHG4Ptcl_truthcluster_Y, AllPHG4Ptcl_truthcluster_Z, AllPHG4Ptcl_truthcluster_edep, AllPHG4Ptcl_truthcluster_adc, AllPHG4Ptcl_truthcluster_r, AllPHG4Ptcl_truthcluster_phi, AllPHG4Ptcl_truthcluster_eta, AllPHG4Ptcl_truthcluster_phisize, AllPHG4Ptcl_truthcluster_zsize, AllPHG4Ptcl_recocluster_globalX,
2404                                      AllPHG4Ptcl_recocluster_globalY, AllPHG4Ptcl_recocluster_globalZ, AllPHG4Ptcl_recocluster_r, AllPHG4Ptcl_recocluster_phi, AllPHG4Ptcl_recocluster_eta, AllPHG4Ptcl_recocluster_phisize, AllPHG4Ptcl_recocluster_zsize, AllPHG4Ptcl_recocluster_adc);
2405     }
2406 
2407     N_PrimaryPHG4Ptcl = PrimaryPHG4Ptcl_trackID.size();
2408     N_sPHENIXPrimary = sPHENIXPrimary_trackID.size();
2409     N_AllPHG4Ptcl = AllPHG4Ptcl_trackID.size();
2410     // final clean up of ancestor info vectors
2411     Clean(ancestor_trackIDs);
2412     Clean(ancestor_PIDs);
2413 }
2414 
2415 //____________________________________________________________________________..
2416 void VertexCompare::Cleanup()
2417 {
2418     ncoll_ = 0;
2419     npart_ = 0;
2420     N_HepMCGenEvent = 0;
2421     Clean(HepMCGenEvent_processID);
2422     Clean(HepMCGenEvent_embeddingID);
2423     Clean(HepMCGenEvent_crossing);
2424     nTotalSilSeeds = 0;
2425     nSilSeedsValidCrossing = 0;
2426     MBD_charge_sum = 0;
2427     mbd_north_npmt = 0;
2428     mbd_south_npmt = 0;
2429     mbd_south_charge_sum = 0;
2430     mbd_north_charge_sum = 0;
2431     mbd_nhitsoverths_south = 0;
2432     mbd_nhitsoverths_north = 0;
2433     N_PrimaryPHG4Ptcl = 0;
2434     N_sPHENIXPrimary = 0;
2435     N_AllPHG4Ptcl = 0;
2436     nTruthVertex = 0;
2437     hasSvtxPHG4ParticleMap = false;
2438     svtxPHG4ParticleMapProcessed = false;
2439     nRecoTracks = 0;
2440     Clean(track_deltapt);
2441     Clean(track_deltaeta);
2442     Clean(track_deltaphi);
2443     Clean(track_nhits);
2444     Clean(track_nmaps);
2445     Clean(track_nintt);
2446     Clean(track_ntpc);
2447     Clean(track_nmms);
2448     Clean(track_ntpc1);
2449     Clean(track_ntpc11);
2450     Clean(track_ntpc2);
2451     Clean(track_ntpc3);
2452     Clean(track_pidedx);
2453     Clean(track_kdedx);
2454     Clean(track_prdedx);
2455     Clean(track_vx);
2456     Clean(track_vy);
2457     Clean(track_vz);
2458     Clean(track_dca2d);
2459     Clean(track_dca2dsigma);
2460     Clean(track_dca3dxy);
2461     Clean(track_dca3dxysigma);
2462     Clean(track_dca3dz);
2463     Clean(track_dca3dzsigma);
2464     Clean(track_hlxpt);
2465     Clean(track_hlxeta);
2466     Clean(track_hlxphi);
2467     Clean(track_hlxX0);
2468     Clean(track_hlxY0);
2469     Clean(track_hlxZ0);
2470     Clean(track_hlxcharge);
2471     Clean(track_id);
2472     Clean(track_x);
2473     Clean(track_y);
2474     Clean(track_z);
2475     Clean(track_px);
2476     Clean(track_py);
2477     Clean(track_pz);
2478     Clean(track_pt);
2479     Clean(track_eta);
2480     Clean(track_phi);
2481     Clean(track_dedx);
2482     Clean(track_charge);
2483     Clean(track_crossing);
2484     Clean(track_vertex_id);
2485     Clean(track_chisq);
2486     Clean(track_ndf);
2487     Clean(track_quality);
2488     Clean(track_silseed_id);
2489     Clean(track_silseed_x);
2490     Clean(track_silseed_y);
2491     Clean(track_silseed_z);
2492     Clean(track_silseed_pt);
2493     Clean(track_silseed_eta);
2494     Clean(track_silseed_phi);
2495     Clean(track_silseed_crossing);
2496     Clean(track_silseed_charge);
2497     Clean(track_silseed_nMvtx);
2498     Clean(track_silseed_nIntt);
2499     Clean(track_silseed_clusterKeys);
2500     Clean(track_cluster_layer);
2501     Clean(track_cluster_globalX);
2502     Clean(track_cluster_globalY);
2503     Clean(track_cluster_globalZ);
2504     Clean(track_cluster_phi);
2505     Clean(track_cluster_eta);
2506     Clean(track_cluster_r);
2507 
2508     Clean(firedTriggers);
2509     Clean(TruthVertex_isEmbeded);
2510     Clean(TruthVertexX);
2511     Clean(TruthVertexY);
2512     Clean(TruthVertexZ);
2513     Clean(TruthVertexT);
2514     Clean(TruthVertex_crossing);
2515 
2516     Clean(mbdVertex);
2517     Clean(mbdVertexId);
2518     Clean(mbdVertexCrossing);
2519 
2520     Clean(trackerVertexId);
2521     Clean(trackerVertexX);
2522     Clean(trackerVertexY);
2523     Clean(trackerVertexZ);
2524     Clean(trackerVertexChisq);
2525     Clean(trackerVertexNdof);
2526     Clean(trackerVertexNTracks);
2527     Clean(trackerVertexCrossing);
2528     Clean(trackerVertexTrackIDs);
2529 
2530     Clean(silseed_id);
2531     Clean(silseed_assocVtxId);
2532     Clean(silseed_x);
2533     Clean(silseed_y);
2534     Clean(silseed_z);
2535     Clean(silseed_pt);
2536     Clean(silseed_eta);
2537     Clean(silseed_phi);
2538     Clean(silseed_eta_vtx);
2539     Clean(silseed_phi_vtx);
2540     Clean(silseed_crossing);
2541     Clean(silseed_charge);
2542     Clean(silseed_nMvtx);
2543     Clean(silseed_nIntt);
2544     Clean(silseed_clusterKeys);
2545     Clean(silseed_cluster_layer);
2546     Clean(silseed_cluster_globalX);
2547     Clean(silseed_cluster_globalY);
2548     Clean(silseed_cluster_globalZ);
2549     Clean(silseed_cluster_phi);
2550     Clean(silseed_cluster_eta);
2551     Clean(silseed_cluster_r);
2552     Clean(silseed_cluster_phiSize);
2553     Clean(silseed_cluster_zSize);
2554     Clean(silseed_cluster_strobeID);
2555     Clean(silseed_cluster_timeBucketID);
2556     Clean(silseed_ngmvtx);
2557     Clean(silseed_ngintt);
2558     Clean(silseed_cluster_gcluster_key);
2559     Clean(silseed_cluster_gcluster_layer);
2560     Clean(silseed_cluster_gcluster_X);
2561     Clean(silseed_cluster_gcluster_Y);
2562     Clean(silseed_cluster_gcluster_Z);
2563     Clean(silseed_cluster_gcluster_r);
2564     Clean(silseed_cluster_gcluster_phi);
2565     Clean(silseed_cluster_gcluster_eta);
2566     Clean(silseed_cluster_gcluster_edep);
2567     Clean(silseed_cluster_gcluster_adc);
2568     Clean(silseed_cluster_gcluster_phiSize);
2569     Clean(silseed_cluster_gcluster_zSize);
2570     Clean(silseed_f4a_nMatched);
2571     Clean(silseed_f4a_truthTrackID);
2572     Clean(silseed_f4a_truthWeight);
2573     Clean(silseed_f4a_bestTrackID);
2574     Clean(silseed_f4a_bestWeight);
2575     Clean(silseed_f4a_bestG4P_PID);
2576     Clean(silseed_f4a_bestG4P_E);
2577     Clean(silseed_f4a_bestG4P_pT);
2578     Clean(silseed_f4a_bestG4P_eta);
2579     Clean(silseed_f4a_bestG4P_phi);
2580     Clean(silseed_f4a_bestG4P_ancestor_trackID);
2581     Clean(silseed_f4a_bestG4P_ancestor_PID);
2582 
2583     nTotalTpcSeeds = 0;
2584     Clean(tpcseed_id);
2585     Clean(tpcseed_x);
2586     Clean(tpcseed_y);
2587     Clean(tpcseed_z);
2588     Clean(tpcseed_pt);
2589     Clean(tpcseed_eta);
2590     Clean(tpcseed_phi);
2591     Clean(tpcseed_crossing);
2592     Clean(tpcseed_crossing_estimate);
2593     Clean(tpcseed_charge);
2594     Clean(tpcseed_nTpc);
2595     Clean(tpcseed_nMms);
2596     Clean(tpcseed_dedx);
2597     Clean(tpcseed_clusterKeys);
2598     Clean(tpcseed_cluster_layer);
2599     Clean(tpcseed_cluster_globalX);
2600     Clean(tpcseed_cluster_globalY);
2601     Clean(tpcseed_cluster_globalZ);
2602     Clean(tpcseed_cluster_phi);
2603     Clean(tpcseed_cluster_eta);
2604     Clean(tpcseed_cluster_r);
2605 
2606     Clean(clusterKey);
2607     Clean(cluster_layer);
2608     Clean(cluster_chip);
2609     Clean(cluster_stave);
2610     Clean(cluster_globalX);
2611     Clean(cluster_globalY);
2612     Clean(cluster_globalZ);
2613     Clean(cluster_phi);
2614     Clean(cluster_eta);
2615     Clean(cluster_r);
2616     Clean(cluster_phiSize);
2617     Clean(cluster_zSize);
2618     Clean(cluster_adc);
2619     Clean(cluster_timeBucketID);
2620     Clean(cluster_crossing);
2621     Clean(cluster_ladderZId);
2622     Clean(cluster_ladderPhiId);
2623     Clean(cluster_LocalX);
2624     Clean(cluster_LocalY);
2625     Clean(cluster_matchedG4P_trackID);
2626     Clean(cluster_matchedG4P_PID);
2627     Clean(cluster_matchedG4P_E);
2628     Clean(cluster_matchedG4P_pT);
2629     Clean(cluster_matchedG4P_eta);
2630     Clean(cluster_matchedG4P_phi);
2631 
2632     Clean(mvtx_seedcluster_key);
2633     Clean(mvtx_seedcluster_layer);
2634     Clean(mvtx_seedcluster_chip);
2635     Clean(mvtx_seedcluster_stave);
2636     Clean(mvtx_seedcluster_globalX);
2637     Clean(mvtx_seedcluster_globalY);
2638     Clean(mvtx_seedcluster_globalZ);
2639     Clean(mvtx_seedcluster_phi);
2640     Clean(mvtx_seedcluster_eta);
2641     Clean(mvtx_seedcluster_r);
2642     Clean(mvtx_seedcluster_phiSize);
2643     Clean(mvtx_seedcluster_zSize);
2644     Clean(mvtx_seedcluster_strobeID);
2645     Clean(mvtx_seedcluster_matchedcrossing);
2646     Clean(mvtx_seedcluster_hitX);
2647     Clean(mvtx_seedcluster_hitY);
2648     Clean(mvtx_seedcluster_hitZ);
2649     Clean(mvtx_seedcluster_hitrow);
2650     Clean(mvtx_seedcluster_hitcol);
2651 
2652     Clean(mvtx_seedcluster_matchedG4P_trackID);
2653     Clean(mvtx_seedcluster_matchedG4P_PID);
2654     Clean(mvtx_seedcluster_matchedG4P_E);
2655     Clean(mvtx_seedcluster_matchedG4P_pT);
2656     Clean(mvtx_seedcluster_matchedG4P_eta);
2657     Clean(mvtx_seedcluster_matchedG4P_phi);
2658     Clean(mvtx_seedcluster_matchedG4P_ancestor_trackID);
2659     Clean(mvtx_seedcluster_matchedG4P_ancestor_PID);
2660 
2661     Clean(PrimaryPHG4Ptcl_pT);
2662     Clean(PrimaryPHG4Ptcl_eta);
2663     Clean(PrimaryPHG4Ptcl_phi);
2664     Clean(PrimaryPHG4Ptcl_E);
2665     Clean(PrimaryPHG4Ptcl_PID);
2666     Clean(PrimaryPHG4Ptcl_trackID);
2667     Clean(PrimaryPHG4Ptcl_originTrackID);
2668     Clean(PrimaryPHG4Ptcl_originVtxID);
2669     Clean(PrimaryPHG4Ptcl_originVtxT);
2670     Clean(PrimaryPHG4Ptcl_originCrossing);
2671     Clean(PrimaryPHG4Ptcl_originIsEmbeded);
2672     Clean(PrimaryPHG4Ptcl_truthevalIsPrimary);
2673     Clean(PrimaryPHG4Ptcl_embedID);
2674     Clean(PrimaryPHG4Ptcl_ParticleClass);
2675     Clean(PrimaryPHG4Ptcl_isStable);
2676     Clean(PrimaryPHG4Ptcl_charge);
2677     Clean(PrimaryPHG4Ptcl_isChargedHadron);
2678     Clean(PrimaryPHG4Ptcl_ancestor_trackID);
2679     Clean(PrimaryPHG4Ptcl_ancestor_PID);
2680     Clean(PrimaryPHG4Ptcl_truthcluster_X);
2681     Clean(PrimaryPHG4Ptcl_truthcluster_Y);
2682     Clean(PrimaryPHG4Ptcl_truthcluster_Z);
2683     Clean(PrimaryPHG4Ptcl_truthcluster_edep);
2684     Clean(PrimaryPHG4Ptcl_truthcluster_adc);
2685     Clean(PrimaryPHG4Ptcl_truthcluster_r);
2686     Clean(PrimaryPHG4Ptcl_truthcluster_phi);
2687     Clean(PrimaryPHG4Ptcl_truthcluster_eta);
2688     Clean(PrimaryPHG4Ptcl_truthcluster_phisize);
2689     Clean(PrimaryPHG4Ptcl_truthcluster_zsize);
2690     Clean(PrimaryPHG4Ptcl_recocluster_globalX);
2691     Clean(PrimaryPHG4Ptcl_recocluster_globalY);
2692     Clean(PrimaryPHG4Ptcl_recocluster_globalZ);
2693     Clean(PrimaryPHG4Ptcl_recocluster_r);
2694     Clean(PrimaryPHG4Ptcl_recocluster_phi);
2695     Clean(PrimaryPHG4Ptcl_recocluster_eta);
2696     Clean(PrimaryPHG4Ptcl_recocluster_phisize);
2697     Clean(PrimaryPHG4Ptcl_recocluster_zsize);
2698     Clean(PrimaryPHG4Ptcl_recocluster_adc);
2699 
2700     Clean(sPHENIXPrimary_pT);
2701     Clean(sPHENIXPrimary_eta);
2702     Clean(sPHENIXPrimary_phi);
2703     Clean(sPHENIXPrimary_E);
2704     Clean(sPHENIXPrimary_PID);
2705     Clean(sPHENIXPrimary_trackID);
2706     Clean(sPHENIXPrimary_originTrackID);
2707     Clean(sPHENIXPrimary_originVtxID);
2708     Clean(sPHENIXPrimary_originVtxT);
2709     Clean(sPHENIXPrimary_originCrossing);
2710     Clean(sPHENIXPrimary_originIsEmbeded);
2711     Clean(sPHENIXPrimary_truthevalIsPrimary);
2712     Clean(sPHENIXPrimary_embedID);
2713     Clean(sPHENIXPrimary_ParticleClass);
2714     Clean(sPHENIXPrimary_isStable);
2715     Clean(sPHENIXPrimary_charge);
2716     Clean(sPHENIXPrimary_isChargedHadron);
2717     Clean(sPHENIXPrimary_ancestor_trackID);
2718     Clean(sPHENIXPrimary_ancestor_PID);
2719     Clean(sPHENIXPrimary_truthcluster_X);
2720     Clean(sPHENIXPrimary_truthcluster_Y);
2721     Clean(sPHENIXPrimary_truthcluster_Z);
2722     Clean(sPHENIXPrimary_truthcluster_edep);
2723     Clean(sPHENIXPrimary_truthcluster_adc);
2724     Clean(sPHENIXPrimary_truthcluster_r);
2725     Clean(sPHENIXPrimary_truthcluster_phi);
2726     Clean(sPHENIXPrimary_truthcluster_eta);
2727     Clean(sPHENIXPrimary_truthcluster_phisize);
2728     Clean(sPHENIXPrimary_truthcluster_zsize);
2729     Clean(sPHENIXPrimary_recocluster_globalX);
2730     Clean(sPHENIXPrimary_recocluster_globalY);
2731     Clean(sPHENIXPrimary_recocluster_globalZ);
2732     Clean(sPHENIXPrimary_recocluster_r);
2733     Clean(sPHENIXPrimary_recocluster_phi);
2734     Clean(sPHENIXPrimary_recocluster_eta);
2735     Clean(sPHENIXPrimary_recocluster_phisize);
2736     Clean(sPHENIXPrimary_recocluster_zsize);
2737     Clean(sPHENIXPrimary_recocluster_adc);
2738 
2739     Clean(AllPHG4Ptcl_pT);
2740     Clean(AllPHG4Ptcl_eta);
2741     Clean(AllPHG4Ptcl_phi);
2742     Clean(AllPHG4Ptcl_E);
2743     Clean(AllPHG4Ptcl_PID);
2744     Clean(AllPHG4Ptcl_trackID);
2745     Clean(AllPHG4Ptcl_originTrackID);
2746     Clean(AllPHG4Ptcl_originVtxID);
2747     Clean(AllPHG4Ptcl_originVtxT);
2748     Clean(AllPHG4Ptcl_originCrossing);
2749     Clean(AllPHG4Ptcl_originIsEmbeded);
2750     Clean(AllPHG4Ptcl_truthevalIsPrimary);
2751     Clean(AllPHG4Ptcl_embedID);
2752     Clean(AllPHG4Ptcl_ancestor_trackID);
2753     Clean(AllPHG4Ptcl_ancestor_PID);
2754     Clean(AllPHG4Ptcl_truthcluster_X);
2755     Clean(AllPHG4Ptcl_truthcluster_Y);
2756     Clean(AllPHG4Ptcl_truthcluster_Z);
2757     Clean(AllPHG4Ptcl_truthcluster_edep);
2758     Clean(AllPHG4Ptcl_truthcluster_adc);
2759     Clean(AllPHG4Ptcl_truthcluster_r);
2760     Clean(AllPHG4Ptcl_truthcluster_phi);
2761     Clean(AllPHG4Ptcl_truthcluster_eta);
2762     Clean(AllPHG4Ptcl_truthcluster_phisize);
2763     Clean(AllPHG4Ptcl_truthcluster_zsize);
2764     Clean(AllPHG4Ptcl_recocluster_globalX);
2765     Clean(AllPHG4Ptcl_recocluster_globalY);
2766     Clean(AllPHG4Ptcl_recocluster_globalZ);
2767     Clean(AllPHG4Ptcl_recocluster_r);
2768     Clean(AllPHG4Ptcl_recocluster_phi);
2769     Clean(AllPHG4Ptcl_recocluster_eta);
2770     Clean(AllPHG4Ptcl_recocluster_phisize);
2771     Clean(AllPHG4Ptcl_recocluster_zsize);
2772     Clean(AllPHG4Ptcl_recocluster_adc);
2773 }
2774 
2775 //____________________________________________________________________________..
2776 int VertexCompare::ResetEvent(PHCompositeNode *topNode) { return Fun4AllReturnCodes::EVENT_OK; }
2777 
2778 //____________________________________________________________________________..
2779 int VertexCompare::EndRun(const int runnumber) { return Fun4AllReturnCodes::EVENT_OK; }
2780 
2781 //____________________________________________________________________________..
2782 int VertexCompare::End(PHCompositeNode *topNode)
2783 {
2784     if (outFile)
2785     {
2786         outFile->cd();
2787         outTree->Write("", TObject::kOverwrite);
2788         outFile->Close();
2789         delete outFile;
2790         outFile = nullptr;
2791         outTree = nullptr;
2792     }
2793 
2794     delete svtx_evalstack;
2795     svtx_evalstack = nullptr;
2796     clustereval = nullptr;
2797     hiteval = nullptr;
2798     truth_eval = nullptr;
2799 
2800     return Fun4AllReturnCodes::EVENT_OK;
2801 }
2802 
2803 //____________________________________________________________________________..
2804 int VertexCompare::Reset(PHCompositeNode *topNode) { return Fun4AllReturnCodes::EVENT_OK; }
2805 
2806 //____________________________________________________________________________..
2807 void VertexCompare::Print(const std::string &what) const {}