Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-10-15 08:17:17 

0001 //-------------------------------------------------------
0002 //  w.j. llope, wjllope@wayne.edu, run-2 2024
0003 //
0004 #include "Collect.h"
0005 
0006 //---- for tracking...
0007 #include <limits>
0008 #include <map>
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010 #include <fun4all/Fun4AllInputManager.h>
0011 #include <fun4all/Fun4AllServer.h>
0012 #include <fun4all/Fun4AllUtils.h>
0013 #include <phool/phool.h>
0014 #include <phool/getClass.h>
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/PHNodeIterator.h>
0017 #include <trackbase/ClusterErrorPara.h>
0018 #include <trackbase/InttDefs.h>
0019 #include <trackbase/MvtxDefs.h>
0020 #include <trackbase/TpcDefs.h>
0021 #include <trackbase/TrackFitUtils.h>
0022 #include <trackbase/TrkrCluster.h>
0023 #include <trackbase/TrkrClusterContainer.h>
0024 #include <trackbase/TrkrHit.h>
0025 #include <trackbase/TrkrHitSet.h>
0026 #include <trackbase/TrkrHitSetContainer.h>
0027 #include <g4detectors/PHG4CylinderGeomContainer.h>
0028 #include <g4detectors/PHG4TpcGeom.h>
0029 #include <g4detectors/PHG4TpcGeomContainer.h>
0030 #include <globalvertex/GlobalVertex.h>
0031 #include <globalvertex/GlobalVertexMap.h>
0032 #include <globalvertex/SvtxVertex.h>
0033 #include <globalvertex/SvtxVertexMap.h>
0034 // #include <globalvertex/SvtxVertexMap_v1.h>
0035 // #include <trackbase_historic/ActsTransformations.h>
0036 // #include <trackbase_historic/SvtxAlignmentState.h>
0037 // #include <trackbase_historic/SvtxAlignmentStateMap.h>
0038 // #include <trackbase_historic/SvtxTrackMap.h>
0039 #include <trackbase_historic/TrackSeed.h>
0040 #include <trackbase_historic/TrackSeedContainer.h>
0041 #include <tpc/TpcGlobalPositionWrapper.h>
0042 #include <tpc/TpcDistortionCorrectionContainer.h>
0043 
0044 #include <trackbase_historic/SvtxTrack.h>
0045 #include <trackbase_historic/SvtxTrackMap.h>
0046 #include <trackbase_historic/SvtxTrackState.h>
0047 // #include <trackbase_historic/SvtxTrack_v4.h>
0048 // #include <trackbase_historic/SvtxTrackMap_v2.h>
0049 // #include <trackbase_historic/SvtxTrackState_v3.h>
0050 #include <trackbase_historic/SvtxTrackSeed_v1.h>
0051 #include <trackbase_historic/TrackAnalysisUtils.h>
0052 
0053 // Calo includes
0054 #include <calobase/RawCluster.h>
0055 #include <calobase/RawClusterContainer.h>
0056 #include <calobase/RawClusterUtility.h>
0057 #include <calobase/RawTower.h>
0058 #include <calobase/RawTowerContainer.h>
0059 #include <calobase/RawTowerGeom.h>
0060 #include <calobase/RawTowerGeomContainer.h>
0061 #include <calotrigger/CaloTriggerInfo.h>
0062 #include <calobase/TowerInfo.h>
0063 #include <calobase/TowerInfoContainer.h>
0064 #include <calobase/TowerInfoDefs.h>
0065 
0066 //#include <KFParticle.h>
0067 //#include <trackbase/ActsGeometry.h>
0068 #include <Acts/Definitions/Algebra.hpp>
0069 
0070 #include <ffaobjects/RunHeader.h>
0071 #include <ffaobjects/EventHeader.h>
0072 #include <ffarawobjects/Gl1Packet.h>
0073 #include <ffarawobjects/Gl1RawHit.h>
0074 
0075 #include <epd/EpdGeom.h>
0076 #include "eventplaneinfo/Eventplaneinfo.h"
0077 #include "eventplaneinfo/EventplaneinfoMap.h"
0078 #include "eventplaneinfo/EventplaneinfoMapv1.h"
0079 #include "eventplaneinfo/Eventplaneinfov1.h"
0080 
0081 #include <TDatabasePDG.h>
0082 
0083 #include <kfparticle_sphenix/KFParticle_Container.h>
0084 #include <kfparticle_sphenix/KFParticle_Tools.h>
0085 #include "KFPTrack.h"
0086 //KFParticle_Tools kf_tools;
0087 //#include <KFParticle.h>
0088 
0089 #include <vector>
0090 using std::cout;
0091 using std::endl;
0092 #include <stdexcept>
0093 
0094 // #include <odbc++/connection.h>
0095 // #include <odbc++/drivermanager.h>
0096 // #include <odbc++/resultset.h>
0097 // #include <odbc++/statement.h>
0098 // #include <odbc++/types.h>
0099 
0100 double EdgeEMeta[97] = {-1.13381,-1.11254,-1.09107,-1.06949,-1.04781,-1.02602,-1.00413,-0.98214,-0.960045,-0.937848,-0.915549,-0.89315,-0.870648,-0.848046,-0.825343,
0101             -0.802541,-0.779637,-0.756634,-0.733531,-0.710327,-0.687024,-0.66362,-0.640116,-0.616512,-0.592806,-0.569,-0.545091,-0.52108,-0.496965,-0.472747,
0102             -0.448423,-0.423994,-0.399457,-0.374812,-0.350056,-0.325188,-0.300207,-0.27511,-0.249896,-0.224561,-0.199104,-0.173523,-0.148255,-0.123759,-0.0996162,
0103             -0.0754161,-0.0511715,-0.0268971,0.,0.0268971,0.0511715,0.0754161,0.0996162,0.123759,0.148255,0.173523,0.199104,0.224561,0.249896,0.27511,
0104             0.300207,0.325188,0.350056,0.374812,0.399457,0.423994,0.448423,0.472747,0.496965,0.52108,0.545091,0.569,0.592806,0.616512,0.640116,
0105             0.66362,0.687024,0.710327,0.733531,0.756634,0.779637,0.802541,0.825343,0.848046,0.870648,0.89315,0.915549,0.937848,0.960045,0.98214,
0106             1.00413,1.02602,1.04781,1.06949,1.09107,1.11254,1.13381};
0107 double EdgeIHeta[25] = {-1.1,-1.00833,-0.916667,-0.825,-0.733333,-0.641667,-0.55,-0.458333,-0.366667,-0.275,-0.183333,-0.0916667,0.,0.0916667,0.183333,
0108             0.275,0.366667,0.458333,0.55,0.641667,0.733333,0.825,0.916667,1.00833,1.1};
0109 double EdgeOHeta[25] = {-1.1,-1.00833,-0.916667,-0.825,-0.733333,-0.641667,-0.55,-0.458333,-0.366667,-0.275,-0.183333,-0.0916667,0.,0.0916667,0.183333,
0110             0.275,0.366667,0.458333,0.55,0.641667,0.733333,0.825,0.916667,1.00833,1.1};
0111 
0112 #define FIXFLOAT(x) std::fixed<<std::setprecision(3)<<(x)
0113 
0114 // Optional "bco" branch (README_SplitTracks.md sec.12a.6): GL1 BCO of the
0115 // trigger frame (one per row, 8 bytes). ON: bco + crossing is the absolute
0116 // crossing, the key Corral uses to remove collisions that appear in two
0117 // overlapping TFs (README_ana573.md). (The silocx/silocy branches added with
0118 // it on 2026-09-23 were removed 2026-09-26 - never used by Corral.)
0119 static constexpr bool ADD_BCO     = true;
0120 
0121 //  enum CAL_LAYER
0122 //   {
0123 //     PRES = 0,
0124 //     CEMC = 1,
0125 //     HCALIN = 2,
0126 //     HCALOUT = 3,
0127 //     OUTER_CEMC = 4,
0128 //     OUTER_HCALIN = 5,
0129 //     OUTER_HCALOUT = 6
0130 //   };
0131 
0132 //----- take arbitrary phi value in sPH coord system and get phi index of calorimeter tower
0133 int Collect::GetIPHI(int layer, double state_phi){
0134     double anphi=-1.;               // number of towers azimuthally
0135     if (layer==SvtxTrack::CEMC){
0136         anphi   = 256.; 
0137     } else if ( layer==SvtxTrack::HCALIN || layer==SvtxTrack::HCALOUT ){    
0138         anphi   =  64.; 
0139     }
0140     double phib = state_phi + 2.0*M_PI/64.;
0141     if (phib < 0.     ){ phib += 2.0*M_PI; } 
0142     if (phib >=2.*M_PI){ phib -= 2.0*M_PI; }
0143     int iphi    = phib / (2.0*M_PI/anphi);
0144     return iphi;
0145 }
0146 //----- take arbitrary eta value in sPH coord system and get eta index of calorimeter tower
0147 int Collect::GetIETA(int layer, double state_eta){
0148     double aneta=-1;
0149     double limlo=-1,limup=-1;
0150     if ( layer==SvtxTrack::CEMC ){
0151         aneta   = 96.;
0152         limlo   = EdgeEMeta[0];
0153         limup   = EdgeEMeta[int(aneta)];    // note there are NETA+1 values in this array!
0154     } else if ( layer==SvtxTrack::HCALIN || layer==SvtxTrack::HCALOUT ){
0155         aneta   = 24.; 
0156         limlo   = EdgeIHeta[0];
0157         limup   = EdgeIHeta[int(aneta)];    // note there are NETA+1 values in this array!
0158     }
0159     int ieta    = -1;
0160     if ( state_eta>=limlo && state_eta<limup ){
0161         for (int ie=0;ie<int(aneta);ie++){ 
0162             if (layer==SvtxTrack::CEMC    && state_eta<EdgeEMeta[ie+1]){ ieta=ie; break; } 
0163             if (layer==SvtxTrack::HCALIN  && state_eta<EdgeIHeta[ie+1]){ ieta=ie; break; } 
0164             if (layer==SvtxTrack::HCALOUT && state_eta<EdgeOHeta[ie+1]){ ieta=ie; break; }
0165         }
0166     }
0167     return ieta;
0168 }
0169 
0170 //____________________________________________________________________________..
0171 Collect::Collect(const std::string &name):SubsysReco(name){
0172     //
0173     fEvtSeqFirst    = -1;
0174     //
0175     kVerbosity      = 0;
0176     DEBUG           = false;
0177     DEBUG2          = false;
0178     foutfileName    = TString("new.root");
0179     for (int i=0;i<MAXNTRACKS;i++){
0180         fTprimary[i]        = 0;
0181     }
0182     fKillSplitTracks    = false;    // off until validated on real data - see README_SplitTracks.md
0183     ievtSeen        = -1;           // incremented at TOP of process_event!
0184     fTtriggerVector.clear();        //
0185     fTbco       = ~0ULL;
0186     fNTracks    = 0;
0187     fTnv0       = 0;
0188     sumNfoundKS = 0;
0189     sumNfoundLA = 0;
0190     sumNfoundAL = 0;
0191     fHasCaloGeo = false;
0192     //
0193     THRESH_em   = 0.150;
0194     THRESH_ih   = 0.012;
0195     THRESH_oh   = 0.06;
0196     THRESH_ep   = 0.30;
0197     Kshort_PDGID    =   310;
0198     Lambda_PDGID    =  3122;
0199     ALambda_PDGID   = -3122;
0200     //
0201     //std::cout << "Collect::Collect(const std::string &name) Calling ctor" << std::endl;
0202 }
0203 
0204 //____________________________________________________________________________..
0205 Collect::~Collect(){
0206     //std::cout << "Collect::~Collect() Calling dtor" << std::endl;
0207 }
0208 
0209 //____________________________________________________________________________..
0210 int Collect::Init(PHCompositeNode *topNode){
0211     //std::cout << "Collect::Init -- Initializing" << std::endl;
0212     //
0213     bool TREETRACKS = true;
0214     bool TREEV0     = true;
0215     //
0216     //---- book
0217     std::cout<<"Collect::Init -- opening "<<foutfileName<<std::endl;
0218     fout    = new TFile(foutfileName,"RECREATE");
0219     //
0220     hnvtx_frame = new TH1D("hnvtx_frame","No. Vertices/TriggerFrame vs triggerframe",10003,-1.5,10001.5);
0221     hdXing      = new TH1D("hdXing","Vertexed Crossings Interval (crossings)",500,-0.5,499.5);
0222     hech_em     = new TH1D("hech_em","channel energy em ",2000,-0.5,1.5);
0223     hech_ih     = new TH1D("hech_ih","channel energy ih ",2000,-0.5,1.5);
0224     hech_oh     = new TH1D("hech_oh","channel energy oh ",2000,-0.5,1.5);
0225     hech_ep     = new TH1D("hech_ep","channel energy epd",2000,-0.5,4.5);
0226     //
0227 //  MatchCounter    = 0;
0228 //  for (int im=0;im<100;im++){
0229 //      hmatch_EMetaphiE[im]    = new TH2D(Form("hmatch_EMetaphiE_%d",im),Form("hmatch_EMetaphiE_%d",im),96,-1.1,1.1,256,0.,2.*M_PI);
0230 //      hmatch_HCetaphiE[im]    = new TH2D(Form("hmatch_HCetaphiE_%d",im),Form("hmatch_HCetaphiE_%d",im),24,-1.1,1.1, 64,0.,2.*M_PI);
0231 //      gmatch_EMetaphi[im]     = new TGraph();     // position of track in eta,phi at EM layer
0232 //      gmatch_HCetaphi[im]     = new TGraph();     // position of track in eta,phi at HC layer
0233 //      gmatch_EMetaphi[im]     ->SetMarkerStyle(20);
0234 //      gmatch_HCetaphi[im]     ->SetMarkerStyle(20);
0235 //      gmatch_EMetaphi[im]     ->SetMarkerColor(418);  // kGreen+2
0236 //      gmatch_HCetaphi[im]     ->SetMarkerColor(600);  // kBlue
0237 //  }
0238     //
0239     if (CALODETAIL){
0240         //
0241         hmapEM_eta  = new TProfile("hmapEM_eta","hmapEM_eta"  ,10.* 96,-1.152,1.152,-1,300.);
0242         hmapEM_phi  = new TProfile("hmapEM_phi","hmapEM_phi"  ,10.*256,0.,2.*M_PI,-1,300.);
0243         hmapIH_eta  = new TProfile("hmapIH_eta","hmapIH_eta"  ,10.* 24,-1.2,  1.2,-1,300.);
0244         hmapIH_phi  = new TProfile("hmapIH_phi","hmapIH_phi"  ,10.* 64,0.,2.*M_PI,-1,300.);
0245         hmapOH_eta  = new TProfile("hmapOH_eta","hmapOH_eta"  ,10.* 24,-1.2,  1.2,-1,300.);
0246         hmapOH_phi  = new TProfile("hmapOH_phi","hmapOH_phi"  ,10.* 64,0.,2.*M_PI,-1,300.);
0247         hmapEM_eta2 = new TProfile("hmapEM_eta2","hmapEM_eta2",10.* 96,-1.152,1.152,-1,300.);
0248         hmapEM_phi2 = new TProfile("hmapEM_phi2","hmapEM_phi2",10.*256,0.,2.*M_PI,-1,300.);
0249         hmapIH_eta2 = new TProfile("hmapIH_eta2","hmapIH_eta2",10.* 24,-1.2,  1.2,-1,300.);
0250         hmapIH_phi2 = new TProfile("hmapIH_phi2","hmapIH_phi2",10.* 64,0.,2.*M_PI,-1,300.);
0251         hmapOH_eta2 = new TProfile("hmapOH_eta2","hmapOH_eta2",10.* 24,-1.2,  1.2,-1,300.);
0252         hmapOH_phi2 = new TProfile("hmapOH_phi2","hmapOH_phi2",10.* 64,0.,2.*M_PI,-1,300.);
0253         hmapEM_ieta = new TProfile("hmapEM_ieta","hmapEM_ieta", 96,0, 96,    -10,10);
0254         hmapEM_iphi = new TProfile("hmapEM_iphi","hmapEM_iphi",256,0,256,    -10,10);
0255         hmapIH_ieta = new TProfile("hmapIH_ieta","hmapIH_ieta", 24,0, 24,    -10,10);
0256         hmapIH_iphi = new TProfile("hmapIH_iphi","hmapIH_iphi", 64,0, 64,    -10,10);
0257         hmapOH_ieta = new TProfile("hmapOH_ieta","hmapOH_ieta", 24,0, 24,    -10,10);
0258         hmapOH_iphi = new TProfile("hmapOH_iphi","hmapOH_iphi", 64,0, 64,    -10,10);
0259         //
0260         for (int iev=0;iev<MAXEVTDETAIL;iev++){
0261             for (int ilayer=0;ilayer<3;ilayer++){
0262                 int nbeta=24,nbphi=64;
0263                 if (ilayer==0){ nbeta=96; nbphi=256;}
0264                 hlayer[ilayer][iev]         = new TH2D(Form("hlayer_%d_%d"      ,ilayer,iev),Form("hlayers_%d_%d"      ,ilayer,iev),nbeta,-0.5,((double)nbeta)-0.5,nbphi,-0.5,((double)nbphi)-0.5);
0265                 if (ilayer==0){ 
0266                     hlayer_efDET[ilayer][iev]   = new TH2D(Form("hlayer_efDET_%d_%d",ilayer,iev),Form("hlayers_efDET_%d_%d",ilayer,iev),nbeta,EdgeEMeta,nbphi,-M_PI,M_PI);
0267                     hPTOT_etaDETem[iev]         = new TH1D(Form("hPTOT_etaDETem_%d",iev),Form("hPTOT_etaDETem_%d",iev),nbeta,EdgeEMeta );
0268                     hPTOT_phiDETem[iev]         = new TH1D(Form("hPTOT_phiDETem_%d",iev),Form("hPTOT_phiDETem_%d",iev),nbphi,-M_PI,M_PI);
0269                 } else if (ilayer==1){
0270                     hlayer_efDET[ilayer][iev]   = new TH2D(Form("hlayer_efDET_%d_%d",ilayer,iev),Form("hlayers_efDET_%d_%d",ilayer,iev),nbeta,EdgeIHeta,nbphi,-M_PI,M_PI);
0271                     hPTOT_etaDETih[iev]         = new TH1D(Form("hPTOT_etaDETih_%d",iev),Form("hPTOT_etaDETih_%d",iev),nbeta,EdgeIHeta );
0272                     hPTOT_phiDETih[iev]         = new TH1D(Form("hPTOT_phiDETih_%d",iev),Form("hPTOT_phiDETih_%d",iev),4*nbphi,-M_PI,M_PI);
0273                 } else if (ilayer==2){
0274                     hlayer_efDET[ilayer][iev]   = new TH2D(Form("hlayer_efDET_%d_%d",ilayer,iev),Form("hlayers_efDET_%d_%d",ilayer,iev),nbeta,EdgeOHeta,nbphi,-M_PI,M_PI);
0275                     hPTOT_etaDEToh[iev]         = new TH1D(Form("hPTOT_etaDEToh_%d",iev),Form("hPTOT_etaDEToh_%d",iev),nbeta,EdgeOHeta );
0276                     hPTOT_phiDEToh[iev]         = new TH1D(Form("hPTOT_phiDEToh_%d",iev),Form("hPTOT_phiDEToh_%d",iev),4*nbphi,-M_PI,M_PI);
0277                 }
0278             }
0279         }
0280         hPedestalEM = new TProfile2D("hPedestalEM","Pedestal EM vs (etaDET,phiDET)",96,EdgeEMeta,256,-M_PI,M_PI,"s");   // ERRORS ARE STD DEV!!!
0281         hPedestalIH = new TProfile2D("hPedestalIH","Pedestal IH vs (etaDET,phiDET)",24,EdgeIHeta, 64,-M_PI,M_PI,"s");   // ERRORS ARE STD DEV!!!
0282         hPedestalOH = new TProfile2D("hPedestalOH","Pedestal OH vs (etaDET,phiDET)",24,EdgeOHeta, 64,-M_PI,M_PI,"s");   // ERRORS ARE STD DEV!!!
0283     //  hPedestalCorrEM = new TH2D("hPedestalCorrEM","Pedestal Corr EM ",100,-THRESH_em,THRESH_em,100,-THRESH_em,THRESH_em);
0284     //  hPedestalCorrIH = new TH2D("hPedestalCorrIH","Pedestal Corr IH ",100,-THRESH_ih,THRESH_ih,100,-THRESH_ih,THRESH_ih);
0285     //  hPedestalCorrOH = new TH2D("hPedestalCorrOH","Pedestal Corr OH ",100,-THRESH_oh,THRESH_oh,100,-THRESH_oh,THRESH_oh);
0286         for (int iet=0;iet<96;iet++){
0287             for (int ifi=0;ifi<256;ifi++){
0288                     hEevt_em[iet][ifi]  = new TH1D(Form("hEevt_em_%d_%d",iet,ifi),Form("hEevt_em ieta=%d iphi=%d",iet,ifi),10003,-1.5,10001.5);
0289                 if (iet<24&&ifi<64){
0290                     hEevt_ih[iet][ifi]  = new TH1D(Form("hEevt_ih_%d_%d",iet,ifi),Form("hEevt_ih ieta=%d iphi=%d",iet,ifi),10003,-1.5,10001.5);
0291                     hEevt_oh[iet][ifi]  = new TH1D(Form("hEevt_oh_%d_%d",iet,ifi),Form("hEevt_oh ieta=%d iphi=%d",iet,ifi),10003,-1.5,10001.5);
0292                 }
0293             }
0294         }
0295     }   // end CALODETAIL
0296     //
0297     outTree = new TTree("outTree","Collect Tree");
0298     //
0299     outTree->Branch("run"     ,&fRunNum     ,"run/I");
0300     outTree->Branch("segment" ,&fSegment    ,"segment/I");
0301     outTree->Branch("evt"     ,&fEvtSeq     ,"evt/I");
0302     outTree->Branch("trigVec" ,&fTtriggerVector);
0303     outTree->Branch("nvtx"    ,&fTnvtx      ,"nvtx/I");
0304     outTree->Branch("vtxntr"  ,&fTvtxntr    ,"vtxntr/I");
0305     outTree->Branch("nvtxXing",&fTnvtxXing  ,"nvtxXing/I"); // vertices in this crossing (>1 = pileup or split vertex); only the hottest one gets a row
0306     outTree->Branch("vtxntr2" ,&fTvtxntr2   ,"vtxntr2/I");  // runner-up vertex attached-track count in this crossing (0 if nvtxXing==1)
0307     outTree->Branch("vtxx"    ,&fTvtxx      ,"vtxx/D");
0308     outTree->Branch("vtxy"    ,&fTvtxy      ,"vtxy/D");
0309     outTree->Branch("vtxz"    ,&fTvtxz      ,"vtxz/D");
0310 //? outTree->Branch("vtxzGlob",&fTvtxzGlob  ,"vtxzGlob/D");
0311     if (ADD_BCO) outTree->Branch("bco"     ,&fTbco       ,"bco/l"); // GL1 BCO of this trigger frame; ~0ULL if no GL1RAWHIT (README_SplitTracks.md sec.12a.6)
0312     outTree->Branch("crossing",&fTcrossing  ,"crossing/I");
0313     outTree->Branch("etotem"  ,&fTTotE_em   ,"etotem/D");
0314     outTree->Branch("etotih"  ,&fTTotE_ih   ,"etotih/D");
0315     outTree->Branch("etotoh"  ,&fTTotE_oh   ,"etotoh/D");
0316     outTree->Branch("etotioh" ,&fTTotE_ioh  ,"etotioh/D");
0317     outTree->Branch("etotepd" ,&fTTotE_epd  ,"etotepd/D");
0318     //
0319     if (TREETRACKS){
0320         cout<<"Collect::Init .. Enabled TREETRACKS"<<endl;
0321         outTree->Branch("ntr"     ,&fNTracks    ,"ntr/I");
0322         outTree->Branch("pid"     ,fTpid        ,"pid[ntr]/B");     // =-1 undefined, =0 pion, =1 kaon, =2 proton (ID'd as V0 daughters) - sec.12c: was /I
0323         outTree->Branch("indv0"   ,fTindv0      ,"indv0[ntr]/S");   // if >=0, this is the index of the v0 mother of this track in v0 tree! - sec.12c: was /I
0324         outTree->Branch("primary" ,fTprimary    ,"primary[ntr]/B"); // if ==1, vertex finder attached this track to the row's vertex (no DCA test since 2026-10-13) - sec.12c: was /I
0325         outTree->Branch("quality" ,fTquality    ,"quality[ntr]/F"); // sec.12c: was /D (this block down through dedxKFP)
0326         outTree->Branch("chisq"   ,fTchisq      ,"chisq[ntr]/F");
0327         outTree->Branch("chg"     ,fTchg        ,"chg[ntr]/F");
0328 //      outTree->Branch("px"      ,fTpx         ,"px[ntr]/D");
0329 //      outTree->Branch("py"      ,fTpy         ,"py[ntr]/D");
0330 //      outTree->Branch("pz"      ,fTpz         ,"pz[ntr]/D");
0331         outTree->Branch("ptot"    ,fTptot       ,"ptot[ntr]/F");
0332         outTree->Branch("eta"     ,fTeta        ,"eta[ntr]/F");
0333         outTree->Branch("phi"     ,fTphi        ,"phi[ntr]/F");
0334         outTree->Branch("pt"      ,fTpt         ,"pt[ntr]/F");
0335         outTree->Branch("seedeta" ,fTseedeta    ,"seedeta[ntr]/F");
0336         outTree->Branch("seedphi" ,fTseedphi    ,"seedphi[ntr]/F");
0337         outTree->Branch("seedpt"  ,fTseedpt     ,"seedpt[ntr]/F");
0338         outTree->Branch("ntpc"    ,fTnhitstpc   ,"ntpc[ntr]/b");    // sec.12c: was /I
0339         outTree->Branch("nmvtx"   ,fTnhitsmvt   ,"nmvtx[ntr]/b");
0340         outTree->Branch("nintt"   ,fTnhitsint   ,"nintt[ntr]/b");
0341         outTree->Branch("dcaxy"   ,fTdcaxy      ,"dcaxy[ntr]/F");
0342         outTree->Branch("dcaz"    ,fTdcaz       ,"dcaz[ntr]/F");
0343         outTree->Branch("dedx70n" ,fTdedx70n    ,"dedx70n[ntr]/b");     //  70% of the hits! - sec.12c: was /I
0344         outTree->Branch("dedx70s" ,fTdedx70s    ,"dedx70s[ntr]/F");     //  70% of the hits! - sec.12c: was /D
0345         outTree->Branch("dedxKFP" ,fTdedxKFP    ,"dedxKFP[ntr]/F");     //  TrackAnalysisUtils::calc_dedx() - see Collect.h - sec.12c: was /D
0346         outTree->Branch("layermask",fTlayermask ,"layermask[ntr]/l");   // bit per absolute TrkrDefs layer (0-2 MVTX, 3-6 INTT, 7-54 TPC) with a cluster on this track - lets outread_stream compute STAR-style SL itself for any pair, sibling or mixed (see README_SplitTracks.md sec. 8)
0347         outTree->Branch("siclukey" ,fTsiclukey   ,"siclukey[ntr][7]/l");    // raw silicon cluster key per absolute Si layer 0-6 (slot == layermask bit index); sentinel ~0ULL means no cluster on that layer (see README_SplitTracks.md sec. 12)
0348         outTree->Branch("tpcsector",fTtpcsector  ,"tpcsector[ntr][48]/b");  // TpcDefs sector*2+side per TPC layer (slot = layer-7, 0-47); 255 = no cluster (see README_SplitTracks.md sec. 12b)
0349         outTree->Branch("tpcarclen",fTtpcarclen  ,"tpcarclen[ntr][48]/B");  // local arc-length within that sector, 2.0 mm/count; valid only where tpcsector!=255
0350         outTree->Branch("tpcz"     ,fTtpcz       ,"tpcz[ntr][48]/S");   // global z, 0.1 mm/count; valid only where tpcsector!=255
0351         if (fAddTowersToTrackTree){
0352             cout<<"Collect::Init .. Adding calorimeter energy grids!"<<endl;
0353             outTree->Branch("stateTI_x" , fTstateTI_x ,"stateTI_x[ntr]/D");
0354             outTree->Branch("stateTI_y" , fTstateTI_y ,"stateTI_y[ntr]/D");
0355             outTree->Branch("stateTI_z" , fTstateTI_z ,"stateTI_z[ntr]/D");
0356             outTree->Branch("stateTO_x" , fTstateTO_x ,"stateTO_x[ntr]/D");
0357             outTree->Branch("stateTO_y" , fTstateTO_y ,"stateTO_y[ntr]/D");
0358             outTree->Branch("stateTO_z" , fTstateTO_z ,"stateTO_z[ntr]/D");
0359             outTree->Branch("stateEM_x" , fTstateEM_x ,"stateEM_x[ntr]/D");
0360             outTree->Branch("stateEM_y" , fTstateEM_y ,"stateEM_y[ntr]/D");
0361             outTree->Branch("stateEM_z" , fTstateEM_z ,"stateEM_z[ntr]/D");
0362             outTree->Branch("stateIH_x" , fTstateIH_x ,"stateIH_x[ntr]/D");
0363             outTree->Branch("stateIH_y" , fTstateIH_y ,"stateIH_y[ntr]/D");
0364             outTree->Branch("stateIH_z" , fTstateIH_z ,"stateIH_z[ntr]/D");
0365             outTree->Branch("stateOH_x" , fTstateOH_x ,"stateOH_x[ntr]/D");
0366             outTree->Branch("stateOH_y" , fTstateOH_y ,"stateOH_y[ntr]/D");
0367             outTree->Branch("stateOH_z" , fTstateOH_z ,"stateOH_z[ntr]/D");
0368             outTree->Branch("egridEM"   , fTegridEM   ,"egridEM[ntr][7][7]/D");
0369             outTree->Branch("egridIH"   , fTegridIH   ,"egridIH[ntr][7][7]/D");
0370             outTree->Branch("egridOH"   , fTegridOH   ,"egridOH[ntr][7][7]/D");
0371             outTree->Branch("eopEM"     , fTeopEM     ,"eopEM[ntr]/D"); // EMCal 7x7-sum / track p ("E/p"), -9 if not available
0372             outTree->Branch("eopHC"     , fTeopHC     ,"eopHC[ntr]/D"); // (IHCal+OHCal) 7x7-sum / track p, -9 if neither available
0373             outTree->Branch("eopCA"     , fTeopCA     ,"eopCA[ntr]/D"); // (EMCal+IHCal+OHCal) 7x7-sum / track p, -9 if none available
0374         }   // end addtowers check...
0375     }
0376     //
0377     if (TREEV0){
0378         cout<<"Collect::Init .. Enabled TREEVO"<<endl;
0379         outTree->Branch("nv0"       ,&fTnv0         ,"nv0/I");
0380         outTree->Branch("v0pid"     ,fTv0pid        ,"v0pid[nv0]/I");
0381         outTree->Branch("v0x"       ,fTv0x          ,"v0x[nv0]/D");
0382         outTree->Branch("v0y"       ,fTv0y          ,"v0y[nv0]/D");
0383         outTree->Branch("v0z"       ,fTv0z          ,"v0z[nv0]/D");
0384         outTree->Branch("v0px"      ,fTv0px         ,"v0px[nv0]/D");
0385         outTree->Branch("v0py"      ,fTv0py         ,"v0py[nv0]/D");
0386         outTree->Branch("v0pz"      ,fTv0pz         ,"v0pz[nv0]/D");
0387         outTree->Branch("v0pt"      ,fTv0pt         ,"v0pt[nv0]/D");
0388         outTree->Branch("v0ptot"    ,fTv0ptot       ,"v0ptot[nv0]/D");
0389         outTree->Branch("v0eta"     ,fTv0eta        ,"v0eta[nv0]/D");
0390         outTree->Branch("v0phi"     ,fTv0phi        ,"v0phi[nv0]/D");
0391         outTree->Branch("v0ene"     ,fTv0ene        ,"v0ene[nv0]/D");
0392         outTree->Branch("v0mass"    ,fTv0mass       ,"v0mass[nv0]/D");
0393         outTree->Branch("v0ctau"    ,fTv0ctau       ,"v0ctau[nv0]/D");
0394         outTree->Branch("v0decaylen",fTv0decaylen   ,"v0decaylen[nv0]/D");
0395         outTree->Branch("v0chi2ndf" ,fTv0chi2ndf    ,"v0chi2ndf[nv0]/D");
0396         outTree->Branch("v0dira"    ,fTv0dira       ,"v0dira[nv0]/D");  // cos(pointing angle) wrt PV
0397         outTree->Branch("v0pvdca"   ,fTv0pvdca      ,"v0pvdca[nv0]/D"); // V0 trajectory DCA to PV (cm)
0398         outTree->Branch("v0decayx"  ,fTv0decayx     ,"v0decayx[nv0]/D");    // daughter-pair decay vtx (what dira/pvdca use, NOT v0x/y/z)
0399         outTree->Branch("v0decayy"  ,fTv0decayy     ,"v0decayy[nv0]/D");
0400         outTree->Branch("v0decayz"  ,fTv0decayz     ,"v0decayz[nv0]/D");
0401         outTree->Branch("v0indtr1"  ,fTv0indtr1     ,"v0indtr1[nv0]/I");    // index to entry for track1 in Tracks OR PFLOW tree!
0402         outTree->Branch("v0indtr2"  ,fTv0indtr2     ,"v0indtr2[nv0]/I");    // index to entry for track2 in Tracks OR PFLOW tree!
0403     }
0404     //
0405     //  to get to EMCal:    pT >~ 0.20
0406     //  to get to IHCal:    pT >~ 0.25
0407     //  to get to OHCal:    pT >~ 0.36
0408     //
0409     PTLL    = 0.1;  // LL -- same track-quality floor regardless of fAddTowersToTrackTree, so track/v0 content doesn't depend on it (0.2 until 2026-09-26: dE/dx vs ptot of kept v0 daughters showed tracking fine below 0.2)
0410     cout<<"Collect::Init -- pT lower limit: PTLL= "<<PTLL<<endl;
0411     //
0412     //
0413     ntpcmin = 20;   // LL
0414     nmvtmin =  2;   // LL
0415     nintmin =  1;   // LL
0416     cout<<"Collect::Init -- Ntpc  lower limit = "<<ntpcmin<<endl;
0417     cout<<"Collect::Init -- Nmvtx lower limit = "<<nmvtmin<<endl;
0418     cout<<"Collect::Init -- Nintt lower limit = "<<nintmin<<endl;
0419     //
0420     return Fun4AllReturnCodes::EVENT_OK;
0421     //
0422     //---- done init
0423 }
0424 
0425 //____________________________________________________________________________..
0426 int Collect::InitRun(PHCompositeNode *topNode){
0427     std::cout << "Collect::InitRun -- Initializing for Run XXX" << std::endl;
0428 
0429     // crossing+distortion-corrected TPC global position (README_SplitTracks.md sec.12b) -
0430     // loads ActsGeometry + the Tpc distortion-correction-container nodes once per run.
0431     m_globalPositionWrapper.loadNodes(topNode);
0432 
0433     if (kVerbosity==0){ DEBUG = false; DEBUG2 = false; } else
0434     if (kVerbosity==1){ DEBUG = true;  DEBUG2 = false; } else
0435     if (kVerbosity>=2){ DEBUG = true;  DEBUG2 = true;  }
0436 
0437     //---- calo tower geometry comes from a Fun4AllRunNodeInputManager-supplied
0438     //---- Run node - not guaranteed populated yet at Init() time, but is by
0439     //---- InitRun() (called once the framework has synced to an actual run).
0440     RawTowerGeomContainer*  EMCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0441     RawTowerGeomContainer*  IHCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0442     RawTowerGeomContainer*  OHCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0443     fHasCaloGeo = (EMCalGeo && IHCalGeo && OHCalGeo);
0444     if (!fHasCaloGeo){
0445         //---- no calo geometry available (e.g. running off a seed-level/tracking-only input) --
0446         //---- track + KFParticle v0 reconstruction still runs, just without the calo-layer-around-track info.
0447         std::cout<<"Collect::InitRun -- calorimeter geometry not found (EMCalGeo="<<(bool)EMCalGeo
0448                  <<" IHCalGeo="<<(bool)IHCalGeo<<" OHCalGeo="<<(bool)OHCalGeo
0449                  <<") -- running without calo-layer-around-track info"<<std::endl;
0450     } else {
0451         caloRadiusInnerEMCal    = EMCalGeo->get_radius();       //  93.5    cm
0452         caloRadiusInnerIHCal    = IHCalGeo->get_radius();       // 115.0    cm
0453         caloRadiusInnerOHCal    = OHCalGeo->get_radius();       // 177.423  cm
0454         caloThicknessEMCal      = EMCalGeo->get_thickness();    //  20.4997 cm
0455         caloThicknessIHCal      = IHCalGeo->get_thickness();    //  25.005  cm
0456         caloThicknessOHCal      = OHCalGeo->get_thickness();    //  96.894  cm
0457         caloRadiusOuterEMCal    = caloRadiusInnerEMCal + caloThicknessEMCal;
0458         caloRadiusOuterIHCal    = caloRadiusInnerIHCal + caloThicknessIHCal;
0459         caloRadiusOuterOHCal    = caloRadiusInnerOHCal + caloThicknessOHCal;
0460         cout<<"Collect::InitRun -- caloRadiusInnerEMCal= "<<caloRadiusInnerEMCal<<endl;
0461         cout<<"Collect::InitRun -- caloRadiusOuterEMCal= "<<caloRadiusOuterEMCal<<endl;
0462         cout<<"Collect::InitRun -- caloRadiusInnerIHCal= "<<caloRadiusInnerIHCal<<endl;
0463         cout<<"Collect::InitRun -- caloRadiusOuterIHCal= "<<caloRadiusOuterIHCal<<endl;
0464         cout<<"Collect::InitRun -- caloRadiusInnerOHCal= "<<caloRadiusInnerOHCal<<endl;
0465         cout<<"Collect::InitRun -- caloRadiusOuterOHCal= "<<caloRadiusOuterOHCal<<endl;
0466         caloRadiusMidEMCal  = caloRadiusInnerEMCal + caloThicknessEMCal/2.;
0467         caloRadiusMidIHCal  = caloRadiusInnerIHCal + caloThicknessIHCal/2.;
0468         caloRadiusMidOHCal  = caloRadiusInnerOHCal + caloThicknessOHCal/2.;
0469         cout<<"Collect::InitRun -- caloRadiusMidEMCal= "<<caloRadiusMidEMCal<<endl;
0470         cout<<"Collect::InitRun -- caloRadiusMidIHCal= "<<caloRadiusMidIHCal<<endl;
0471         cout<<"Collect::InitRun -- caloRadiusMidOHCal= "<<caloRadiusMidOHCal<<endl;
0472     }
0473 
0474     return Fun4AllReturnCodes::EVENT_OK;
0475 }
0476 
0477 // ============================================================================
0478 // DIAGTEMP-20260920: throwaway merged-track exploration (README_SplitTracks.md
0479 // sec.12b discussion, 2026-09-20). NOT for commit - every block tagged
0480 // DIAGTEMP-20260920 below (and in process_event()/End()) should be deleted
0481 // once this exploration is done. Goal: (1) how often does a literal TPC
0482 // TrkrDefs::cluskey get claimed by 2+ tracks in the same event, and (2) does
0483 // per-cluster (untruncated) adc show a second population near ~2x the bulk
0484 // peak for tracks in a broad, charge-agnostic kinematic-proximity window (the
0485 // STAR "double-pion dE/dx band" the user described from 30 years at STAR).
0486 // ============================================================================
0487 struct DiagTempClusterInfo20260920 { int layer; float adc; TrkrDefs::cluskey key; float gx; float gy; float gz; };
0488 static TH1D* sDiagAdcAll20260920        = nullptr;  // every TPC cluster's normalized adc, all tracks
0489 static TH1D* sDiagAdcSharedKey20260920  = nullptr;  // adc of clusters whose key is claimed by 2+ tracks this event
0490 static TH1D* sDiagAdcClosePair20260920  = nullptr;  // adc of clusters on a shared occupied layer between a broad-window close pair
0491 static TH1D* sDiagDist20260920          = nullptr;  // DIAGTEMP-20260920 v2: 3D global-position separation (cm) between same-layer cluster pairs of a close track pair
0492 static TH2D* sDiagAdcVsDist20260920     = nullptr;  // DIAGTEMP-20260920 v2: adc vs that separation - the real STAR-style FMH-adjacent test
0493 static TProfile* sDiagGapRateVsDist20260920 = nullptr;  // DIAGTEMP-20260920 v4: rate of a layermask hit-presence asymmetry on L+-1 vs the L-separation of the pair
0494 static long  sDiagNTracks20260920       = 0;
0495 static long  sDiagNClusters20260920     = 0;
0496 static long  sDiagNSharedKeyInstances20260920 = 0;  // (track,cluster) instances where the key is also claimed elsewhere
0497 static long  sDiagNEventsWithSharedKey20260920 = 0;
0498 static long  sDiagNClosePairsChecked20260920   = 0;
0499 static long  sDiagNClosePairsSharedKeyToo20260920 = 0;
0500 static long  sDiagNGapChecksClose20260920 = 0;  // DIAGTEMP-20260920 v4: neighbor-layer checks where the confirming pair's dist3d<2cm
0501 static long  sDiagNGapAsymClose20260920   = 0;
0502 static long  sDiagNGapChecksFar20260920   = 0;  // dist3d>10cm control
0503 static long  sDiagNGapAsymFar20260920     = 0;
0504 // DIAGTEMP-20260920 v5: is TpcDefs::getSectorId/getSide(key) - decodable from the raw
0505 // cluster key with zero position data, exactly like siclukey's design - a good enough
0506 // proxy for "these two clusters are really close" to be usable downstream without
0507 // persisting full 3D position? Test it against the real dist3d we already compute.
0508 static TH1D* sDiagDistSameSector20260920 = nullptr;
0509 static TH1D* sDiagDistDiffSector20260920 = nullptr;
0510 static TProfile* sDiagGapRateVsDistSameSector20260920 = nullptr;
0511 static TProfile* sDiagGapRateVsDistDiffSector20260920 = nullptr;
0512 
0513 // DIAGTEMP-20260920: per-vertex analysis, called right after FlagSplitTracks()
0514 // so it sees the same fNTracks/fTeta/fTphi/fTchg this vertex just filled.
0515 static void DiagMergedTracksExplore20260920(
0516         int nTracks, const uint64_t* tlayermask,
0517         const std::vector<std::vector<DiagTempClusterInfo20260920>>& clustersByTrack,
0518         const std::map<TrkrDefs::cluskey,std::vector<int>>& keyOwners){
0519     //
0520     if (!sDiagAdcAll20260920){
0521         sDiagAdcAll20260920       = new TH1D("hDiagAdcAll20260920","DIAGTEMP adc, all TPC clusters;normalized adc",300,0,6000);
0522         sDiagAdcSharedKey20260920 = new TH1D("hDiagAdcSharedKey20260920","DIAGTEMP adc, shared-key clusters;normalized adc",300,0,6000);
0523         sDiagAdcClosePair20260920 = new TH1D("hDiagAdcClosePair20260920","DIAGTEMP adc, close-pair shared-layer clusters;normalized adc",300,0,6000);
0524         sDiagDist20260920         = new TH1D("hDiagDist20260920","DIAGTEMP 3D separation, same-layer clusters of a close track pair;cm",300,0,30);
0525         sDiagAdcVsDist20260920    = new TH2D("hDiagAdcVsDist20260920","DIAGTEMP adc vs 3D separation;cm;normalized adc",150,0,30,300,0,6000);
0526         sDiagGapRateVsDist20260920 = new TProfile("hDiagGapRateVsDist20260920","DIAGTEMP L+-1 hit-presence asymmetry rate vs 3D separation;cm;asymmetry rate",150,0,30);
0527         sDiagAdcAll20260920->SetDirectory(0);
0528         sDiagAdcSharedKey20260920->SetDirectory(0);
0529         sDiagAdcClosePair20260920->SetDirectory(0);
0530         sDiagDist20260920->SetDirectory(0);
0531         sDiagAdcVsDist20260920->SetDirectory(0);
0532         sDiagGapRateVsDist20260920->SetDirectory(0);
0533         sDiagDistSameSector20260920 = new TH1D("hDiagDistSameSector20260920","DIAGTEMP dist3d, TpcDefs sector+side match;cm",300,0,30);
0534         sDiagDistDiffSector20260920 = new TH1D("hDiagDistDiffSector20260920","DIAGTEMP dist3d, TpcDefs sector+side differ;cm",300,0,30);
0535         sDiagGapRateVsDistSameSector20260920 = new TProfile("hDiagGapRateVsDistSameSector20260920","DIAGTEMP gap rate vs dist, sector+side match;cm;rate",150,0,30);
0536         sDiagGapRateVsDistDiffSector20260920 = new TProfile("hDiagGapRateVsDistDiffSector20260920","DIAGTEMP gap rate vs dist, sector+side differ;cm;rate",150,0,30);
0537         sDiagDistSameSector20260920->SetDirectory(0);
0538         sDiagDistDiffSector20260920->SetDirectory(0);
0539         sDiagGapRateVsDistSameSector20260920->SetDirectory(0);
0540         sDiagGapRateVsDistDiffSector20260920->SetDirectory(0);
0541     }
0542     //
0543     bool eventHasSharedKey = false;
0544     for (int it=0; it<nTracks; it++){
0545         sDiagNTracks20260920++;
0546         for (const auto &ci : clustersByTrack[it]){
0547             sDiagNClusters20260920++;
0548             sDiagAdcAll20260920->Fill(ci.adc);
0549             auto ownIt = keyOwners.find(ci.key);
0550             if (ownIt!=keyOwners.end() && ownIt->second.size()>=2){
0551                 sDiagNSharedKeyInstances20260920++;
0552                 sDiagAdcSharedKey20260920->Fill(ci.adc);
0553                 eventHasSharedKey = true;
0554             }
0555         }
0556     }
0557     if (eventHasSharedKey) sDiagNEventsWithSharedKey20260920++;
0558     //
0559     // DIAGTEMP-20260920 v3: NO kinematic (deta,dphi) pre-gate here anymore - v1/v2
0560     // copied the split-track gate, which assumes one real particle refit twice
0561     // (so its fragments share nearly the same trajectory everywhere). Two merged
0562     // tracks are two different particles that can be close at one radius and
0563     // diverge elsewhere in eta/phi - a global kinematic gate would (and, per the
0564     // v2 distance histogram, apparently did) throw away most of the real
0565     // candidate population before it ever reached the per-layer position check.
0566     // Scan every same-event pair directly instead; with ~4-5 tracks/vertex here
0567     // this is trivially cheap (O(N^2 x 48 layers)), so no shortcut is needed.
0568     for (int i=0;i<nTracks;i++){
0569         for (int j=i+1;j<nTracks;j++){
0570             sDiagNClosePairsChecked20260920++;
0571             bool pairSharesKeyToo = false;
0572             // index clusters of j by layer for a quick lookup
0573             std::map<int,const DiagTempClusterInfo20260920*> jByLayer;
0574             for (const auto &cj : clustersByTrack[j]) jByLayer[cj.layer] = &cj;
0575             for (const auto &ci : clustersByTrack[i]){
0576                 auto itj = jByLayer.find(ci.layer);
0577                 if (itj==jByLayer.end()) continue;
0578                 sDiagAdcClosePair20260920->Fill(ci.adc);
0579                 sDiagAdcClosePair20260920->Fill(itj->second->adc);
0580                 if (ci.key==itj->second->key) pairSharesKeyToo = true;
0581                     // DIAGTEMP-20260920 v2: real 3D separation between the two same-layer clusters
0582                     double ddx = ci.gx - itj->second->gx;
0583                     double ddy = ci.gy - itj->second->gy;
0584                     double ddz = ci.gz - itj->second->gz;
0585     double dist3d = sqrt(ddx*ddx + ddy*ddy + ddz*ddz);
0586                     sDiagDist20260920->Fill(dist3d);
0587                     sDiagAdcVsDist20260920->Fill(dist3d, ci.adc);
0588                     sDiagAdcVsDist20260920->Fill(dist3d, itj->second->adc);
0589                     // DIAGTEMP-20260920 v5: TpcDefs::getSectorId/getSide decode from the raw key
0590                     // alone (no position needed) - is "same sector+side" a good enough proxy for
0591                     // "really close" to be usable downstream from a persisted key, instead of
0592                     // needing full 3D position (much bigger branch, see README sec.12c size debate)?
0593                     bool sameSector = (TpcDefs::getSectorId(ci.key)==TpcDefs::getSectorId(itj->second->key)) &&
0594                                        (TpcDefs::getSide(ci.key)==TpcDefs::getSide(itj->second->key));
0595                     if (sameSector){ sDiagDistSameSector20260920->Fill(dist3d); }
0596                     else            { sDiagDistDiffSector20260920->Fill(dist3d); }
0597                     // DIAGTEMP-20260920 v4: does a neighboring TPC layer show a hit-presence
0598                     // asymmetry (one track has it, the other doesn't) right next to this
0599                     // confirmed close encounter? User's hypothesis: with hit assignment
0600                     // apparently globally exclusive (zero shared keys, section 12 chat
0601                     // 2026-09-20), two genuinely close real tracks can't double-count a
0602                     // contested cluster the way STAR's tracker tolerated - the "loser"
0603                     // just goes missing at that padrow instead. Test on L-1/L+1 (not L
0604                     // itself, which by construction has a hit for both here).
0605                     for (int dl : {-1, 1}){
0606                         int nl = ci.layer + dl;
0607                         if (nl<7 || nl>54) continue;
0608                         bool hitI = (tlayermask[i] >> nl) & 1ULL;
0609                         bool hitJ = (tlayermask[j] >> nl) & 1ULL;
0610                         bool asym = (hitI != hitJ);
0611                         sDiagGapRateVsDist20260920->Fill(dist3d, asym?1.0:0.0);
0612                         if (dist3d<2.0){ sDiagNGapChecksClose20260920++; if (asym) sDiagNGapAsymClose20260920++; }
0613                         if (dist3d>10.0){ sDiagNGapChecksFar20260920++;  if (asym) sDiagNGapAsymFar20260920++;  }
0614                         if (sameSector){ sDiagGapRateVsDistSameSector20260920->Fill(dist3d, asym?1.0:0.0); }
0615                         else            { sDiagGapRateVsDistDiffSector20260920->Fill(dist3d, asym?1.0:0.0); }
0616                     }
0617             }
0618             if (pairSharesKeyToo) sDiagNClosePairsSharedKeyToo20260920++;
0619         }
0620     }
0621 }
0622 // ============================================================================
0623 // end DIAGTEMP-20260920 declarations
0624 // ============================================================================
0625 
0626 //____________________________________________________________________________..
0627 int Collect::process_event(PHCompositeNode *topNode){
0628 
0629     bool SHOWGRIDS = false;     // show calorimeter energy grids if DEBUG=true
0630 
0631     ++ievtSeen;
0632 
0633     //---- get some collections...
0634     //
0635     auto trackmap       = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0636     auto clustermap     = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0637     if (!clustermap) clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER_SEED");   // ana573 seed/track DSTs - same fallback as KFParticle_Tools.cc
0638 //? auto globvertexmap  = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0639     auto svtxvertexmap  = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0640     auto tpcGeom        = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0641     auto geometry       = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0642     //auto mvtxGeom     = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
0643     //auto inttGeom     = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
0644     //auto mmGeom       = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
0645     //
0646     TowerInfoContainer* towersEM    = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");  //
0647     TowerInfoContainer* towersIH    = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0648     TowerInfoContainer* towersOH    = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0649     RawTowerGeomContainer* EMCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0650     RawTowerGeomContainer* IHCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0651     RawTowerGeomContainer* OHCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");  
0652 
0653     
0654     //---- header info...
0655     EventHeader *evtHeader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0656     //int bunch     = -1;
0657     //float epang   = 0;
0658     if (!evtHeader){
0659         if (DEBUG) std::cout << "Collect::process_event -- no EventHeader... " << std::endl;
0660         fRunNum = -1;   // flag rather than silently keep the previous event's run/evt numbers
0661         fEvtSeq = -1;
0662     } else {
0663         fRunNum = evtHeader->get_RunNumber();
0664         fEvtSeq = evtHeader->get_EvtSequence();
0665         //if (DEBUG) cout<<fRunNum<<" "<<fEvtSeq<<endl;
0666     }
0667     //---- per-event segment from the file the named input manager is reading right now
0668     //     (a job may read several segments of one run - see setSegmentFromInputManager())
0669     if (!fSegmentInputManager.empty()){
0670         Fun4AllInputManager *im = Fun4AllServer::instance()->getInputManager(fSegmentInputManager);
0671         if (im && !im->FileName().empty()){
0672             fSegment = Fun4AllUtils::GetRunSegment(im->FileName()).second;
0673         } else {
0674             fSegment = -1;  // flag rather than silently keep the previous segment
0675         }
0676     }
0677     if (fEvtSeqFirst<0){
0678         fEvtSeqFirst    = fEvtSeq;
0679         cout<<"Collect::process_event -- fEvtSeqFirst= "<<fEvtSeqFirst<<"\t ievtSeen= "<<ievtSeen<<endl;    
0680     }
0681 
0682 //  if (fEvtSeq==11||fEvtSeq==13||fEvtSeq==15){ DEBUG = true; } else { DEBUG = false; }
0683     if (DEBUG) cout<<"COLLECT START ------------- \t fEvtSeq-0= "<<fEvtSeq-fEvtSeqFirst<<"\t ievtSeen= "<<ievtSeen<<endl;
0684 
0685     //---- Grab the GL1 data
0686     Gl1Packet* gl1  = findNode::getClass<Gl1Packet>(topNode, "GL1RAWHIT");
0687     if (!gl1){
0688         if (DEBUG) std::cout<<"Collect::process_event -- no GL1RAWHIT..."<<std::endl;
0689         fTtriggerVector.clear();    // flag rather than silently keep the previous event's trigger vector
0690         fTbco = ~0ULL;
0691     } else {
0692         fTbco = gl1->getBCO();
0693         //uint64_t evtseq     = gl1->getEvtSequence();
0694         uint64_t triggervec = gl1->getScaledVector();
0695         fTtriggerVector.clear();
0696         for (int i = 0; i < 64; i++){
0697             bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0698             if (DEBUG){ if (trig_decision){ cout<<"1"; }else{ cout<<"0"; } }
0699             fTtriggerVector.push_back(trig_decision);
0700             triggervec = (triggervec >> 1U) & 0xffffffffU;
0701         }   if (DEBUG){ cout<<endl; }
0702     }
0703 
0704     //---- KFParticle area...
0705     //
0706     nv0tf   = 0;    // number of V0's total in this TRIGGER FRAME across all three V0 species
0707     int NfoundKS    = 0;
0708     int NfoundLA    = 0;
0709     int NfoundAL    = 0;
0710     //
0711     NfoundKS = GetKFP( topNode, Kshort_PDGID,  "KFPrecoKS_KFParticle_Container", "KFPrecoKS_SvtxTrackMap" );
0712     NfoundLA = GetKFP( topNode, Lambda_PDGID,  "KFPrecoLA_KFParticle_Container", "KFPrecoLA_SvtxTrackMap" );
0713     NfoundAL = GetKFP( topNode, ALambda_PDGID, "KFPrecoAL_KFParticle_Container", "KFPrecoAL_SvtxTrackMap" );
0714     //
0715     sumNfoundKS+=NfoundKS;  sumNfoundLA+=NfoundLA;  sumNfoundAL+=NfoundAL;
0716     if ((fEvtSeq-fEvtSeqFirst)%10==0) std::cout << "processing event " << fEvtSeq << " (" << fEvtSeq-fEvtSeqFirst << ") from run " << fRunNum
0717                     << ".  SUM Nk,l,a: " << sumNfoundKS << " " << sumNfoundLA << " " << sumNfoundAL << std::endl;
0718     bool v0found    = false;
0719     if (NfoundKS||NfoundLA||NfoundAL) v0found = true;
0720     if (v0found && DEBUG){
0721         std::cout   <<"KFP Done... nv0 = "<<fTnv0
0722                     //<<"\t Nk,l,a= "<<NfoundKS<<" "<<NfoundLA<<" "<<NfoundAL
0723                     <<"\t SUM Nk,l,a= "<<sumNfoundKS<<" "<<sumNfoundLA<<" "<<sumNfoundAL
0724                     <<endl;
0725     }
0726     //
0727     //---- reject any v0's that share the same daughter tracks! (in the same crossing obvs) 
0728     //----  pick the v0 with the lower chi2/ndf...
0729     if (nv0tf>1){
0730         for (int iv0=0;iv0<nv0tf;iv0++){
0731             if (fTv0TFuseit[iv0]){
0732                 unsigned int trkid1 = V0Track1id[iv0];
0733                 unsigned int trkid2 = V0Track2id[iv0];
0734                 for (int jv0=0;jv0<iv0;jv0++){
0735                     if (fTv0TFuseit[jv0]){
0736                         if ( V0Crossing[iv0]==V0Crossing[jv0]   // make sure these two V0s are even in the same crossing!
0737                          && ( trkid1 == V0Track1id[jv0] ||      //    is trk1 duplicated in 2 different V0s? 
0738                               trkid2 == V0Track2id[jv0]) ){     // OR is trk2 duplicated in 2 different V0s?
0739                             double chi2i    = fTv0TFchi2ndf[iv0];
0740                             double chi2j    = fTv0TFchi2ndf[jv0];
0741                             if (chi2i>=chi2j){ fTv0TFuseit[iv0] = false; } else
0742                             if (chi2i< chi2j){ fTv0TFuseit[jv0] = false; }
0743                             if (DEBUG) cout<<"V0 DAUGHTER TRACK DUPLICATED .. chi2: "<<chi2i<<" "<<chi2j
0744                                 <<"\t iv0: "<<trkid1<<" "<<trkid2<<"\t jv0: "<<V0Track1id[jv0]<<" "<<V0Track2id[jv0]
0745                                 <<endl;
0746                         }   // end check of track ids used in these two V0s...
0747                     }   // end check use jv0...
0748                 }   // end loop jv0...
0749             }   // end check use iv0...
0750         }   // end loop iv0...
0751     }   // end check of >1 V0
0752     //
0753     //---- end KFP....
0754 
0755 
0756     //---- calorimeter info...
0757     //
0758     fTTotE_em       = 0.;
0759     fTTotE_ih       = 0.;
0760     fTTotE_oh       = 0.;
0761     fTTotE_ioh      = 0.;
0762     fTTotE_epd      = 0.;
0763     double TOTE_em  = 0.;
0764     double TOTE_ih  = 0.;
0765     double TOTE_oh  = 0.;
0766     double TOTE_ioh = 0.;
0767     double TOTE_epd = 0.;
0768     //
0769     double calophishift = 2.0*M_PI/64.;
0770     //
0771     if(towersEM){
0772         for (unsigned int i=0; i<towersEM->size(); i++){            // i is towerIndex...
0773             TowerInfo* towerInfo    = towersEM->get_tower_at_channel(i);
0774             double etow             = towerInfo->get_energy();
0775             if (CALODETAIL){
0776                 unsigned int towkey     = towersEM->encode_key(i);
0777                 int ieta     = towersEM->getTowerEtaBin(towkey);
0778                 int iphi     = towersEM->getTowerPhiBin(towkey);
0779                 double eta=-9,phi=-9;
0780                 if (ievtSeen<MAXEVTDETAIL){
0781                     RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
0782                     RawTowerGeom* towergeom = EMCalGeo->get_tower_geometry(key);
0783                     int iphit=-9;
0784                     if (towergeom){
0785                         eta     = towergeom->get_eta();             // eta-DET
0786                         phi     = towergeom->get_phi();             // phi-DET
0787                         double phib     = phi + calophishift;       // shift is 4*bw
0788                         if (phib< 0.     ) phib += 2.*M_PI;
0789                         if (phib>=2.*M_PI) phib -= 2.*M_PI;
0790                         double phibw    = 2.0*M_PI/256.;
0791                         iphit           = phib/phibw;
0792                         hmapEM_eta  ->Fill(eta,ieta);
0793                         hmapEM_phi  ->Fill(phi,iphi);
0794         //              hmapEM_eta2 ->Fill(eta,ietat-ieta);
0795                         hmapEM_phi2 ->Fill(phi,iphit-iphi);
0796                         hmapEM_ieta ->Fill(ieta,eta);
0797                         hmapEM_iphi ->Fill(iphi,phib);
0798                     } 
0799                 }
0800                 hEevt_em[ieta][iphi]            ->Fill(fEvtSeq,etow);
0801                 if (ievtSeen<MAXEVTDETAIL){
0802                     hlayer[0][ievtSeen]         ->Fill(ieta,iphi,etow);
0803                     hlayer_efDET[0][ievtSeen]   ->Fill( eta, phi,etow);
0804                 }
0805                 if (etow>=-THRESH_em&&etow<THRESH_em){ 
0806                     if (CALODETAIL) hPedestalEM ->Fill(eta,phi,etow);   
0807                 }
0808             }   // end CALODETAIL
0809             if (etow>=THRESH_em){ 
0810                 TOTE_em     += etow;
0811             }
0812             hech_em     ->Fill(etow);
0813 //          delete towergeom;
0814         }   // end tower loop
0815     }   // end towersEM
0816     //
0817     if(towersIH){
0818         for(unsigned int i = 0; i < towersIH->size(); i++){
0819             TowerInfo* towerInfo    = towersIH->get_tower_at_channel(i);
0820             double etow             = towerInfo->get_energy();
0821             if (CALODETAIL){
0822                 unsigned int towkey     = towersIH->encode_key(i);
0823                 int ieta     = towersIH->getTowerEtaBin(towkey);
0824                 int iphi     = towersIH->getTowerPhiBin(towkey);
0825                 double eta=-9,phi=-9;
0826                 if (ievtSeen<MAXEVTDETAIL){
0827                     RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0828                     RawTowerGeom* towergeom = IHCalGeo->get_tower_geometry(key);
0829                     int iphit=-9;
0830                     if (towergeom){
0831                         eta     = towergeom->get_eta();             // eta-DET
0832                         phi     = towergeom->get_phi();             // phi-DET
0833                         double phib     = phi + calophishift;       // shift is 1*bw
0834                         if (phib< 0.     ) phib += 2.*M_PI;
0835                         if (phib>=2.*M_PI) phib -= 2.*M_PI;
0836                         double phibw    = 2.0*M_PI/64.;
0837                         iphit           = phib/phibw;
0838                         //cout<<eta<<" "<<ieta<<" "<<ietat<<"\t "<<phi<<" "<<iphi<<" "<<iphit<<endl;
0839                         hmapIH_eta  ->Fill(eta,ieta);
0840                         hmapIH_phi  ->Fill(phi,iphi);
0841         //              hmapIH_eta2 ->Fill(eta,ietat-ieta);
0842                         hmapIH_phi2 ->Fill(phi,iphit-iphi);
0843                         hmapIH_ieta ->Fill(ieta,eta);
0844                         hmapIH_iphi ->Fill(iphi,phib);
0845                     } 
0846                 }
0847                 hEevt_ih[ieta][iphi]            ->Fill(fEvtSeq,etow);
0848                 if (ievtSeen<MAXEVTDETAIL){
0849                     hlayer[1][ievtSeen]         ->Fill(ieta,iphi,etow);
0850                     hlayer_efDET[1][ievtSeen]   ->Fill( eta, phi,etow);
0851                 }
0852                 if (etow>=-THRESH_ih&&etow<THRESH_ih){ 
0853                     if (CALODETAIL) hPedestalIH ->Fill(eta,phi,etow);   
0854                 }
0855             }   // end CALODETAIL
0856             if (etow>=THRESH_ih){ 
0857                 TOTE_ih     += etow;
0858             }
0859             hech_ih     ->Fill(etow);
0860             //delete towergeom;
0861         }   // end tower loop
0862     }   // end towersEM
0863     //
0864     if(towersOH){
0865         for(unsigned int i = 0; i < towersOH->size(); i++){
0866             TowerInfo* towerInfo    = towersOH->get_tower_at_channel(i);
0867             double etow             = towerInfo->get_energy();
0868             if (CALODETAIL){
0869                 unsigned int towkey     = towersOH->encode_key(i);
0870                 int ieta     = towersOH->getTowerEtaBin(towkey);
0871                 int iphi     = towersOH->getTowerPhiBin(towkey);
0872                 double eta=-9,phi=-9;
0873                 if (ievtSeen<MAXEVTDETAIL){
0874                     RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0875                     RawTowerGeom* towergeom = OHCalGeo->get_tower_geometry(key);
0876                     int iphit=-9;
0877                     if (towergeom){
0878                         eta     = towergeom->get_eta();         // eta-DET
0879                         phi     = towergeom->get_phi();         // phi-DET
0880                         double phib     = phi + calophishift;   // shift is 1*bw
0881                         if (phib< 0.     ) phib += 2.*M_PI;
0882                         if (phib>=2.*M_PI) phib -= 2.*M_PI;
0883                         double phibw    = 2.0*M_PI/64.;
0884                         iphit           = phib/phibw;
0885                         //cout<<eta<<" "<<ieta<<" "<<ietat<<"\t "<<phi<<" "<<iphi<<" "<<iphit<<endl;
0886                         hmapOH_eta  ->Fill(eta,ieta);
0887                         hmapOH_phi  ->Fill(phi,iphi);
0888         //              hmapOH_eta2 ->Fill(eta,ietat-ieta);
0889                         hmapOH_phi2 ->Fill(phi,iphit-iphi);
0890                         hmapOH_ieta ->Fill(ieta,eta);
0891                         hmapOH_iphi ->Fill(iphi,phib);
0892                     } 
0893                 }
0894                 hEevt_oh[ieta][iphi]            ->Fill(fEvtSeq,etow);
0895                 if (ievtSeen<MAXEVTDETAIL){
0896                     hlayer[2][ievtSeen]         ->Fill(ieta,iphi,etow);
0897                     hlayer_efDET[2][ievtSeen]   ->Fill( eta, phi,etow);
0898                 }
0899                 if (etow>=-THRESH_oh&&etow<THRESH_oh){
0900                     if (CALODETAIL) hPedestalOH ->Fill(eta,phi,etow);
0901                 }
0902             }   // end CALODETAIL
0903             if (etow>=THRESH_oh){
0904                 TOTE_oh     += etow;
0905             }
0906             hech_oh     ->Fill(etow);
0907             //delete towergeom;
0908         }   // end tower loop
0909     }   // end towersEM
0910     TOTE_ioh    = TOTE_ih + TOTE_oh;
0911     //
0912     //---- sepd info...
0913     TowerInfoContainer* epdtowerinfo    = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_SEPD");
0914     EpdGeom* epdgeom                    = findNode::getClass<EpdGeom>(topNode, "TOWERGEOM_EPD");
0915     if ( !epdtowerinfo || !epdgeom ){   // require both!
0916         // EPD not present in this dataset - silently skip, not an error.
0917     } else {
0918         TowerInfo *towerInfo = nullptr;
0919         for(unsigned int i=0; i<epdtowerinfo->size(); i++){
0920             towerInfo = epdtowerinfo->get_tower_at_channel(i);
0921             double etow  = towerInfo->get_energy();
0922             if (etow>=THRESH_ep) TOTE_epd   += etow;
0923             hech_ep     ->Fill(etow);
0924         }
0925         //
0926         unsigned int ntowers = epdtowerinfo->size();
0927         if (DEBUG){ cout<<"epd ntowers = "<<ntowers<<" tote_epd = "<<TOTE_epd<<endl; }
0928         //
0929     }
0930     //
0931     if (DEBUG) cout<<"Eemc="<<TOTE_em<<"\t Eih="<<TOTE_ih<<"\t Eoh="<<TOTE_oh<<"\t Eepd="<<TOTE_epd<<endl;
0932 
0933     //
0934     std::vector<int> Xings;
0935     Xings.clear();              // is this needed?????
0936     //
0937 
0938     //---- loop over tracking vertices (from tpc/silicon)
0939     //
0940     int nverticeskept   = 0;
0941     if (svtxvertexmap){
0942         fTnvtx      = svtxvertexmap->size();        // number of vertices found in this event
0943         hnvtx_frame ->Fill(fEvtSeq-fEvtSeqFirst,fTnvtx);
0944         if (DEBUG) cout<<"----- Vertex NVERT = "<<fTnvtx<<endl;
0945         //
0946         //---- Protect against multiple vertices sharing the same crossing (real
0947         //     pileup, or a split-vertex reconstruction artifact - either way,
0948         //     giving both their own row double-counts every track/v0 whose
0949         //     crossing matches, since the track/v0 inclusion below is
0950         //     crossing-based, not vertex-based - see the collect_v0_track_index_bug
0951         //     memory, "vertex-row duplication" finding, for how this was found).
0952         //     Keep only the vertex with the most attached tracks per crossing.
0953         //     Independent of fRequireXingZero, which is a separate decision
0954         //     (whether to look at crossing zero only) - this generalizes what
0955         //     used to be a crossing-zero-only protection to every crossing.
0956         std::map<int, unsigned int> vertexIdKeepByXing; // crossing -> vertex id with the most attached tracks
0957         std::map<int, unsigned int> maxNtrkByXing;          // crossing -> that vertex's track count
0958         std::map<int, unsigned int> secondNtrkByXing;       // crossing -> runner-up vertex's track count (0 if none)
0959         std::map<int, int>          nVtxByXing;             // crossing -> number of vertices in it
0960         for (const auto& [key, vertex] : *svtxvertexmap){   // loop over vertices...
0961             int          m_vertex_crossing  = vertex->get_beam_crossing();
0962             // Hotness metric: count tracks in this crossing with get_vertex_id()==key
0963             // (2026-10-13: the 1.5 cm DCA test was dropped - Collect applies no DCA
0964             // selection anywhere; the choice only matters for nvtxXing>1 rows),
0965             // NOT gated by the ntpc/nmvtx/nintt/pT track-quality
0966             // floor used for tree inclusion. Author request: a genuinely busy
0967             // real vertex can have many real but lower-quality tracks attached,
0968             // and should not lose the "hottest" pick to a vertex with fewer but
0969             // higher-quality tracks just because of that floor. Replaces
0970             // vertex->size_tracks(), an opaque vertex-finder-internal count we
0971             // do not otherwise use or trust anywhere else in this file.
0972             unsigned int ntrk = 0;
0973             if (trackmap){
0974                 for (auto &iter2 : *trackmap){
0975                     SvtxTrack *trk2 = iter2.second;
0976                     if (!trk2) continue;
0977                     if (trk2->get_crossing()  != m_vertex_crossing) continue;
0978                     if (trk2->get_vertex_id() != key) continue;
0979                     ntrk++;
0980                 }
0981             }
0982             nVtxByXing[m_vertex_crossing]++;
0983             auto it = maxNtrkByXing.find(m_vertex_crossing);
0984             if (it==maxNtrkByXing.end() || ntrk>it->second){
0985                 if (it!=maxNtrkByXing.end()) secondNtrkByXing[m_vertex_crossing] = it->second;  // old leader becomes runner-up
0986                 maxNtrkByXing[m_vertex_crossing]        = ntrk;
0987                 vertexIdKeepByXing[m_vertex_crossing]   = key;
0988             } else if (ntrk>secondNtrkByXing[m_vertex_crossing]){
0989                 secondNtrkByXing[m_vertex_crossing]     = ntrk;
0990             }
0991         }
0992         //
0993         for (const auto& [key, vertex] : *svtxvertexmap){   // loop over vertices...
0994             unsigned int m_vertexid     = key;
0995             int    m_vertex_crossing    = vertex->get_beam_crossing();
0996             //
0997             if (fRequireXingZero && m_vertex_crossing!=0) continue;     // keep only crossing zero if requested
0998             if (m_vertexid != vertexIdKeepByXing[m_vertex_crossing]){
0999                 if (DEBUG) cout<<"Skipping vertex ID= "<<m_vertexid<<"\t crossing= "<<m_vertex_crossing
1000                             <<"\t Ntracks= "<<vertex->size_tracks()
1001                             <<" (keeping ID= "<<vertexIdKeepByXing[m_vertex_crossing]<<" instead)"<<endl;
1002                 continue;       // not the hottest vertex at this crossing - skip to avoid double-counting
1003             }
1004             //
1005             Xings.push_back(m_vertex_crossing);
1006             double m_vx                 = vertex->get_x();
1007             double m_vy                 = vertex->get_y();
1008             double m_vz                 = vertex->get_z();
1009             int    m_ntracks            = maxNtrkByXing[m_vertex_crossing]; // attached-track count, not vertex->size_tracks() (see above)
1010             if (DEBUG) cout <<"-----   vtx: "<<m_vertexid<<" crossing="<<m_vertex_crossing
1011                             <<" x="<<m_vx<<" y="<<m_vy<<" z="<<m_vz
1012                             <<"\t ntrkattached = "<<m_ntracks<<" ------ "<<endl;
1013             fTvtxntr            = m_ntracks;            // to outTree... total number of tracks assigned to chosen vertex
1014             fTnvtxXing          = nVtxByXing[m_vertex_crossing];        // to outTree... vertices sharing this crossing (CF analyses skip >1)
1015             fTvtxntr2           = secondNtrkByXing[m_vertex_crossing];  // to outTree... runner-up vertex attached-track count (0 if none)
1016             fTvtxx              = m_vx;                 // to outTree... PV x 
1017             fTvtxy              = m_vy;                 // to outTree... PV y
1018             fTvtxz              = m_vz;                 // to outTree... PV z
1019             fTcrossing          = m_vertex_crossing;
1020             if (m_vertex_crossing == 0){
1021                 fTTotE_em       = TOTE_em;
1022                 fTTotE_ih       = TOTE_ih;
1023                 fTTotE_oh       = TOTE_oh;
1024                 fTTotE_ioh      = TOTE_ioh;
1025                 fTTotE_epd      = TOTE_epd;
1026             } else {
1027                 fTTotE_em=fTTotE_ih=fTTotE_oh=fTTotE_ioh=fTTotE_epd=-1.0;
1028             }
1029 //?         fTvtxzGlob  = Crossing2GlobVtxZ[m_vertex_crossing];
1030             //
1031             fNTracks    = 0;
1032             fTnv0       = 0;
1033             std::vector<std::vector<DiagTempClusterInfo20260920>> diagClustersByTrack20260920;  // DIAGTEMP-20260920
1034             std::map<TrkrDefs::cluskey,std::vector<int>> diagKeyOwners20260920;                 // DIAGTEMP-20260920
1035             for (int it=0;it<MAXNTRACKS;it++){  //---- before track loop, zero arrays for towers near tracks...
1036                 //for (int idphi=-3;idphi<=3;idphi++){
1037                 //  for (int ideta=-3;ideta<=3;ideta++){
1038                 //      fTegridEM[it][3+ideta][3+idphi] = 0.;
1039                 //      fTegridIH[it][3+ideta][3+idphi] = 0.;
1040                 //      fTegridOH[it][3+ideta][3+idphi] = 0.;
1041                 //  }
1042                 //}
1043                 fTstateTI_x[it] = fTstateTO_x[it] = 0.; 
1044                 fTstateTI_y[it] = fTstateTO_y[it] = 0.; 
1045                 fTstateTI_z[it] = fTstateTO_z[it] = 0.; 
1046                 fTstateEM_x[it] = fTstateIH_x[it] = fTstateOH_x[it] = 0.; 
1047                 fTstateEM_y[it] = fTstateIH_y[it] = fTstateOH_y[it] = 0.; 
1048                 fTstateEM_z[it] = fTstateIH_z[it] = fTstateOH_z[it] = 0.; 
1049                 fTegridEMsum7x7[it] = fTegridIHsum7x7[it] = fTegridOHsum7x7[it] = 0.;
1050                 fTeopEM[it] = fTeopIH[it] = fTeopOH[it] = -9.;
1051             }
1052             //
1053             //---- find the V0s in this crossing.
1054             if (DEBUG2) cout<<"finding v0s in this crossing:  nv0tf= "<<nv0tf<<endl;
1055             std::vector<int> trk1inds;
1056             std::vector<int> trk2inds;
1057             for (int iv0=0;iv0<nv0tf;iv0++){
1058                 if (m_vertex_crossing   == V0Crossing[iv0]  // take V0s from /this/ crossing
1059                  && fTv0TFuseit[iv0]    ){                  // & make sure these V0s have distinct daughters!
1060                     fTv0pid[fTnv0]      = fTv0TFpid[iv0];
1061                     fTv0x[fTnv0]        = fTv0TFx[iv0];
1062                     fTv0y[fTnv0]        = fTv0TFy[iv0];
1063                     fTv0z[fTnv0]        = fTv0TFz[iv0];
1064                     fTv0px[fTnv0]       = fTv0TFpx[iv0];
1065                     fTv0py[fTnv0]       = fTv0TFpy[iv0];
1066                     fTv0pz[fTnv0]       = fTv0TFpz[iv0];
1067                     fTv0pt[fTnv0]       = fTv0TFpt[iv0];
1068                     fTv0ptot[fTnv0]     = fTv0TFptot[iv0];
1069                     {   // DIRA & PV_DCA: daughter-pair decay vtx (V0DecayX/Y/Z, NOT fTv0x/y/z - see
1070                         // Collect.h comment: the mother KFParticle stored in fTv0x/y/z has already
1071                         // been TransportToProductionVertex()'d by KFParticle_DST before Collect ever
1072                         // reads it, so it sits ~at the PV, not the decay vertex - confirmed empirically
1073                         // 2026-09-17, |v0pos-PV| was uncorrelated with v0decaylen and dira was scattered
1074                         // from -0.75 to +0.1 on a DIRA>=0.99 sample) vs PV (m_vx/m_vy/m_vz), along fTv0p{x,y,z}
1075                         fTv0decayx[fTnv0]   = V0DecayX[iv0];
1076                         fTv0decayy[fTnv0]   = V0DecayY[iv0];
1077                         fTv0decayz[fTnv0]   = V0DecayZ[iv0];
1078                         double dx   = fTv0decayx[fTnv0] - m_vx;
1079                         double dy   = fTv0decayy[fTnv0] - m_vy;
1080                         double dz   = fTv0decayz[fTnv0] - m_vz;
1081                         double dmag = sqrt(dx*dx + dy*dy + dz*dz);
1082                         double pmag = fTv0ptot[fTnv0];
1083                         fTv0dira[fTnv0]     = (dmag>0. && pmag>0.) ? (fTv0px[fTnv0]*dx + fTv0py[fTnv0]*dy + fTv0pz[fTnv0]*dz)/(dmag*pmag) : -9.;
1084                         double crx  = dy*fTv0pz[fTnv0] - dz*fTv0py[fTnv0];
1085                         double cry  = dz*fTv0px[fTnv0] - dx*fTv0pz[fTnv0];
1086                         double crz  = dx*fTv0py[fTnv0] - dy*fTv0px[fTnv0];
1087                         fTv0pvdca[fTnv0]    = (pmag>0.) ? sqrt(crx*crx + cry*cry + crz*crz)/pmag : -9.;
1088                     }
1089                     fTv0eta[fTnv0]      = fTv0TFeta[iv0];
1090                     fTv0phi[fTnv0]      = fTv0TFphi[iv0];   
1091                     fTv0ene[fTnv0]      = fTv0TFene[iv0];   
1092                     fTv0mass[fTnv0]     = fTv0TFmass[iv0];
1093                     fTv0ctau[fTnv0]     = fTv0TFctau[iv0];      
1094                     fTv0decaylen[fTnv0] = fTv0TFdecaylen[iv0];  
1095                     fTv0chi2ndf[fTnv0]  = fTv0TFchi2ndf[iv0];   
1096                     fTv0indtr1[fTnv0]   = fTv0TFindtr1[iv0];    // still set in vertex&track loops below!
1097                     fTv0indtr2[fTnv0]   = fTv0TFindtr2[iv0];    // still set in vertex&track loops below!
1098                     fV0FinalToTF[fTnv0] = iv0;                  // remember which trigger-frame slot this final slot came from (dedup skips shift these out of alignment!)
1099                     ++fTnv0;
1100                 }
1101             }
1102             if (DEBUG2) cout<<"done finding v0s in this crossing:  nv0tf= "<<nv0tf<<endl;
1103             //
1104             //---- loop over tracks attached to this vertex...
1105             //for (SvtxVertex::TrackIter iter2=vertex->begin_tracks(); iter2!=vertex->end_tracks();++iter2 ){
1106             //  SvtxTrack *track        = trackmap->get(*iter2);
1107             //---- loop over ALL tracks in this trigger frame, and require crossing #s match vertex->track!
1108             //
1109             //
1110             if (DEBUG2) cout<<"looping over tracks... "<<trackmap->size()<<endl;
1111             for (auto &iter : *trackmap){
1112                 SvtxTrack *track    = iter.second;
1113                 if (!track) continue;
1114                 int m_track_crossing    = track->get_crossing();
1115                 if (m_track_crossing != m_vertex_crossing) continue; 
1116                 //
1117                 unsigned int trackid    = track->get_id();
1118                 double chg              = track->get_charge();
1119                 double px               = track->get_px();
1120                 double py               = track->get_py();
1121                 double pz               = track->get_pz();
1122                 double chisq            = track->get_chisq();
1123                 double quality          = track->get_quality();
1124                 double pt               = sqrt(px*px + py*py);
1125                 double ptot             = sqrt(pt*pt + pz*pz);
1126                 double eta              = track->get_eta();
1127                 double phi              = track->get_phi();
1128                 double dca3d_xy         = NAN;
1129                 double dca3d_xy_error   = NAN;
1130                 double dca3d_z          = NAN;
1131                 double dca3d_z_error    = NAN;
1132                 if (DEBUG2) cout<<"TrackLoop: getting dca..."<<endl;
1133                 if (svtxvertexmap){
1134                     calc3DDCA(track, svtxvertexmap, m_vertexid, dca3d_xy, dca3d_xy_error, dca3d_z, dca3d_z_error);
1135                 } else {
1136                     cout<<"no vertex map."<<endl;
1137                 }
1138                 //
1139                 int vertexFlag  = 0;
1140                 if (track->get_vertex_id() == m_vertexid) vertexFlag = 1;   // this track is attached to this vertex!  (otherwise just same crossing)
1141                 // 2026-10-13: no DCA test here any more - "primary" is purely the vertex finder's
1142                 // attachment; dcaxy/dcaz are in the tree, so any DCA selection is made downstream (Corral AcceptTrack).
1143                 //
1144                 if (DEBUG2) cout<<"TrackLoop: doing dedx..."<<endl;
1145                 double  dedx70n     =    0;
1146                 double  dedx70s     = -999;
1147                 double  dedxKFP     = -999;
1148                 int     nmaps   =    0;
1149                 int     nintt   =    0;
1150                 int     ntpc    =    0;
1151                 int     nmms    =    0;
1152                 //double  seedpx    =    0; 
1153                 //double  seedpy    =    0; 
1154                 //double  seedpz    =    0; 
1155                 double  seedpt  =    0;
1156                 double  seedeta =    0;
1157                 double  seedphi =    0;
1158                 uint64_t layermask  = 0;    // one bit per absolute TrkrDefs layer with a cluster on this track - see fTlayermask in Collect.h
1159                 TrkrDefs::cluskey siclukey_local[7] = {~0ULL,~0ULL,~0ULL,~0ULL,~0ULL,~0ULL,~0ULL};  // raw Si cluster key per layer 0-6; ~0ULL sentinel (0 itself is a valid key - see README_SplitTracks.md sec. 12.3) - see fTsiclukey in Collect.h
1160                 UChar_t tpcsector_local[48];    // see fTtpcsector in Collect.h - README_SplitTracks.md sec. 12b
1161                 Char_t  tpcarclen_local[48];
1162                 Short_t tpcz_local[48];
1163                 for (int itl=0;itl<48;itl++){ tpcsector_local[itl]=255; tpcarclen_local[itl]=0; tpcz_local[itl]=0; }
1164                 std::vector<DiagTempClusterInfo20260920> thisTrackClusters20260920; // DIAGTEMP-20260920
1165                 auto tpcseed = track->get_tpc_seed();
1166                 if (tpcseed){
1167                     //
1168                     //seedpx    = tpcseed->get_px();
1169                     //seedpy    = tpcseed->get_py();
1170                     //seedpz    = tpcseed->get_pz();
1171                     //seedpt    = sqrt(seedpx*seedpx + seedpy*seedpy);
1172                     seedpt  = tpcseed->get_pt();
1173                     seedeta = tpcseed->get_eta();
1174                     seedphi = tpcseed->get_phi();
1175                     bool  signsame = std::signbit(eta) == std::signbit(seedeta);
1176                     if ( !signsame ) seedeta = -seedeta;
1177                     //
1178                     std::vector<double> dedxvalues  = calc_dedx(tpcseed, clustermap, tpcGeom);
1179                     //
1180                     if (dedxvalues.size()==2){
1181                         dedx70n     = dedxvalues.at(0);     //     Nclus ,  70% truncation
1182                         dedx70s     = dedxvalues.at(1);     //     <dedx>,  70% truncation
1183                     }
1184                     //
1185                     // same call KFParticle_Tools::get_dEdx() makes internally for the usePID()
1186                     // dE/dx cut - written as its own branch (dedxKFP) so it can be compared
1187                     // directly against our own dedx70s above without reimplementing anything
1188                     // (see README_localKFPchanges.md, 2026-09-18 dE/dx scale hack entry)
1189                     {
1190                         float layerThicknesses[4] = {0.0, 0.0, 0.0, 0.0};
1191                         layerThicknesses[0] = tpcGeom->GetLayerCellGeom(7 )->get_thickness();
1192                         layerThicknesses[1] = tpcGeom->GetLayerCellGeom(8 )->get_thickness();
1193                         layerThicknesses[2] = tpcGeom->GetLayerCellGeom(27)->get_thickness();
1194                         layerThicknesses[3] = tpcGeom->GetLayerCellGeom(50)->get_thickness();
1195                         dedxKFP = TrackAnalysisUtils::calc_dedx(tpcseed, clustermap, geometry, layerThicknesses);
1196                     }
1197                     //
1198                     for (auto it = tpcseed->begin_cluster_keys(); it != tpcseed->end_cluster_keys(); ++it){
1199                         auto ckey = *it;
1200                         if (TrkrDefs::getTrkrId(ckey) == TrkrDefs::TrkrId::tpcId       ){ ntpc++;  } 
1201                         if (TrkrDefs::getTrkrId(ckey) == TrkrDefs::TrkrId::micromegasId){ nmms++;  } 
1202                         if (TrkrDefs::getTrkrId(ckey) == TrkrDefs::TrkrId::tpcId ){
1203                             auto cluster    = clustermap->findCluster(ckey);
1204                             if (!cluster) continue;     //!!!!! PROTECTION DST MISMATCH !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1205                             auto surf       = geometry->maps().getSurface(ckey, cluster);
1206                             //if (surf){
1207                             //  Acts::Vector3 global_in = geometry->getGlobalPosition(key, cluster);
1208                             //}
1209                             unsigned int cluslayer  = TrkrDefs::getLayer(ckey);
1210                             if (cluslayer<=54){ layermask |= (1ULL << cluslayer); } // TPC bits (7-54) of fTlayermask
1211                         if (cluslayer>=7 && cluslayer<=54){
1212                             // Merged-track support (README_SplitTracks.md sec.12b): compact per-layer
1213                             // local position. Real merging can only happen within one TPC sector
1214                             // (sectors are separated by dead-space frames), so we persist a sector ID
1215                             // plus a local (arc-length, z) position rather than full global 3D position -
1216                             // same information content where it matters, far fewer bytes. See Collect.h
1217                             // for the fixed mm/count scales (chosen with margin, not per-layer geometry).
1218                             int tpcidx = (int)cluslayer - 7;
1219                             Acts::Vector3 tpcGlobalPos = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, (short)m_track_crossing);
1220                             double gx = tpcGlobalPos.x();
1221                             double gy = tpcGlobalPos.y();
1222                             double gz = tpcGlobalPos.z();
1223                             double phiGlobal = atan2(gy, gx);
1224                             double radiusCm  = sqrt(gx*gx + gy*gy);
1225                             constexpr double TPC_SECTOR_WIDTH_RAD = 2.0*M_PI / 12.0;    // TpcDefs::NSectors
1226                             double phiWrapped = fmod(phiGlobal, TPC_SECTOR_WIDTH_RAD);
1227                             if (phiWrapped < 0) phiWrapped += TPC_SECTOR_WIDTH_RAD;
1228                             if (phiWrapped > TPC_SECTOR_WIDTH_RAD/2.0) phiWrapped -= TPC_SECTOR_WIDTH_RAD;
1229                             double localArcLenMM = radiusCm * 10.0 * phiWrapped;    // cm -> mm
1230                             long arclenCount = lround(localArcLenMM / 2.0); // 2.0 mm/count, see Collect.h
1231                             if (arclenCount >  127) arclenCount =  127;
1232                             if (arclenCount < -127) arclenCount = -127;
1233                             long zCount = lround(gz * 100.0);   // 0.1 mm/count, see Collect.h
1234                             if (zCount >  32767) zCount =  32767;
1235                             if (zCount < -32767) zCount = -32767;
1236                             int sectorId = (int)TpcDefs::getSectorId(ckey);
1237                             int sideId   = (int)TpcDefs::getSide(ckey);
1238                             tpcsector_local[tpcidx] = (UChar_t)(sectorId*2 + sideId);   // 0-23, never 255 (the sentinel)
1239                             tpcarclen_local[tpcidx] = (Char_t)arclenCount;
1240                             tpcz_local[tpcidx]      = (Short_t)zCount;
1241                         }
1242                             {   // DIAGTEMP-20260920: same normalized-adc formula as calc_dedx(), just tagged with layer+key
1243                                 PHG4TpcGeom* diagGeo20260920 = tpcGeom->GetLayerCellGeom(cluslayer);
1244                                 if (diagGeo20260920){
1245                                     float diagAdc20260920   = cluster->getAdc();
1246                                     float diagThick20260920 = diagGeo20260920->get_thickness();
1247                                     float diagR20260920     = diagGeo20260920->get_radius();
1248                                     float diagAlpha20260920 = (diagR20260920*diagR20260920) / (2*diagR20260920*TMath::Abs(1.0/tpcseed->get_qOverR()));
1249                                     float diagBeta20260920  = atan(tpcseed->get_slope());
1250                                     diagAdc20260920 /= diagThick20260920;
1251                                     diagAdc20260920 *= cos(diagAlpha20260920);
1252                                     diagAdc20260920 *= cos(diagBeta20260920);
1253                                     // DIAGTEMP-20260920 v2: real global position, crossing-corrected (v6, 2026-09-20 - see README sec.12b)
1254                                     Acts::Vector3 diagGlobalPos20260920 = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, (short)m_track_crossing);
1255                                     thisTrackClusters20260920.push_back({(int)cluslayer, diagAdc20260920, ckey,
1256                                         (float)diagGlobalPos20260920.x(), (float)diagGlobalPos20260920.y(), (float)diagGlobalPos20260920.z()});
1257                                 }
1258                             }
1259                         }   // end tpcId check...
1260                     }   // end cluster key loop...
1261                     //
1262                 }   // tpcseed exists...
1263                 //if (DEBUG) cout<<"TrackLoop: done doing dedx..."<<endl;
1264                 //
1265                 if (DEBUG2) cout<<"TrackLoop: count silicon hits..."<<endl;
1266                 auto silseed = track->get_silicon_seed();
1267                 if (silseed){
1268                     for (auto it = silseed->begin_cluster_keys(); it != silseed->end_cluster_keys(); ++it){
1269                         auto ckey = *it;
1270                         unsigned int siliconlayer = TrkrDefs::getLayer(ckey);
1271                         switch (TrkrDefs::getTrkrId(ckey)){
1272                             case TrkrDefs::mvtxId:
1273                                 nmaps++;
1274                                 layermask |= (1ULL << siliconlayer);    // MVTX bits (0-2) of fTlayermask
1275                                 siclukey_local[siliconlayer] = ckey;    // MVTX slots (0-2) of fTsiclukey
1276                                 break;
1277                             case TrkrDefs::inttId:
1278                                 nintt++;
1279                                 layermask |= (1ULL << siliconlayer);    // INTT bits (3-6) of fTlayermask
1280                                 siclukey_local[siliconlayer] = ckey;    // INTT slots (3-6) of fTsiclukey
1281                                 break;
1282                         }   // ensd detId check...
1283                     }   // end cluster key loop...
1284                 }   // silseed exists...
1285                 //if (DEBUG) cout<<"TrackLoop: done counting silicon hits..."<<endl;
1286                 //
1287                 //---- good track selection (daughters from v0s are kept even if they fail this goodtrack cut!)
1288                 bool KEEPTRK                  = true;
1289                 if (ntpc   < ntpcmin  ){ KEEPTRK = false; } else    // open in general, tight if adding towers to track tree!
1290                 if (nmaps  < nmvtmin  ){ KEEPTRK = false; } else
1291                 if (nintt  < nintmin  ){ KEEPTRK = false; } else
1292                 if (pt     < PTLL     ){ KEEPTRK = false; } else
1293                 if (pt    >= PTLL+20. ){ KEEPTRK = false; }
1294                 //
1295                 //
1296                 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1297                 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1298                 // Collect::process_event -- caloRadiusInnerEMCal= 93.5
1299                 // Collect::process_event -- caloRadiusOuterEMCal= 114
1300                 // Collect::process_event -- caloRadiusInnerIHCal= 115
1301                 // Collect::process_event -- caloRadiusOuterIHCal= 140.005
1302                 // Collect::process_event -- caloRadiusInnerOHCal= 177.423
1303                 // Collect::process_event -- caloRadiusOuterOHCal= 274.317
1304                 // Collect::process_event -- caloRadiusMidEMCal= 103.75
1305                 // Collect::process_event -- caloRadiusMidIHCal= 127.502
1306                 // Collect::process_event -- caloRadiusMidOHCal= 225.87
1307                 //
1308                 //---- look at the calorimeters near this track!
1309                 if (DEBUG2) cout<<"TrackLoop: Swim to calos..."<<endl;
1310                 //
1311                 //SvtxTrackState *thisState = nullptr;
1312                 SvtxTrackState *thisStateI  = nullptr;
1313                 SvtxTrackState *thisStateO  = nullptr;
1314                 double stx=0,sty=0,stz=0;
1315                 //
1316                 //---- loop over states on this track to find first & last tpc states...    
1317                 double rhomin       = 999.; 
1318                 double rhomax       =   0.; 
1319                 for (SvtxTrack::ConstStateIter state_itr = track->begin_states();
1320                                                state_itr != track->end_states(); 
1321                                                state_itr++){    
1322                     SvtxTrackState* astate = dynamic_cast<SvtxTrackState*>(state_itr->second);
1323                     stx = astate->get_x();
1324                     sty = astate->get_y();
1325                     stz = astate->get_z();
1326                     double rho  = sqrt( stx*stx + sty*sty );
1327                     //cout<<track->get_id()<<"\t "<<astate->get_pathlength()<<"\t "<<rho<<endl;
1328                     if (rho>30. && rho<rhomin){
1329                         rhomin                  = rho;
1330                         fTstateTI_x[fNTracks]   = stx;  // first TPC state on track
1331                         fTstateTI_y[fNTracks]   = sty;  // first TPC state on track
1332                         fTstateTI_z[fNTracks]   = stz;  // first TPC state on track
1333                     }
1334                     if (rho<76. && rho>rhomax){
1335                         rhomax                  = rho;
1336                         fTstateTO_x[fNTracks]   = stx;  // last TPC state on track
1337                         fTstateTO_y[fNTracks]   = sty;  // last TPC state on track
1338                         fTstateTO_z[fNTracks]   = stz;  // last TPC state on track
1339                     }
1340                 }
1341                 //
1342                 if (fHasCaloGeo && fAddTowersToTrackTree){
1343                 //---- take track to emcal depth
1344                 thisStateI  = track->get_state(caloRadiusInnerEMCal);
1345                 thisStateO  = track->get_state(caloRadiusOuterEMCal);
1346                 if ( thisStateI && !thisStateO){ 
1347                     stx = thisStateI->get_x();
1348                     sty = thisStateI->get_y();
1349                     stz = thisStateI->get_z();
1350                 } else if ( !thisStateI && thisStateO){ 
1351                     stx = thisStateO->get_x();
1352                     sty = thisStateO->get_y();
1353                     stz = thisStateO->get_z();
1354                 } else if ( thisStateI && thisStateO){ 
1355                     stx = (thisStateI->get_x() + thisStateO->get_x())/2.;
1356                     sty = (thisStateI->get_y() + thisStateO->get_y())/2.;
1357                     stz = (thisStateI->get_z() + thisStateO->get_z())/2.;
1358                 } 
1359                 if (thisStateI||thisStateO){
1360                     fTstateEM_x[fNTracks]   = stx;
1361                     fTstateEM_y[fNTracks]   = sty;
1362                     fTstateEM_z[fNTracks]   = stz;
1363                     fTstateEM_eta[fNTracks] = asinh(stz/sqrt(stx*stx + sty*sty));
1364                     fTstateEM_phi[fNTracks] = atan2(sty,stx);
1365                     int ieta                = GetIETA(SvtxTrack::CEMC, fTstateEM_eta[fNTracks]);
1366                     int iphi                = GetIPHI(SvtxTrack::CEMC, fTstateEM_phi[fNTracks]);
1367                     if (ieta>=0&&ieta<96&&iphi>=0&&iphi<256){   // >>>---> ieta,iphi in hand for this state! 
1368                         //RawTowerDefs::keytype keyRaw  = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
1369                         //RawTowerGeom* towergeom       = EMCalGeo->get_tower_geometry(keyRaw);
1370                         //double toweta                 = towergeom->get_eta();
1371                         //double towphi                 = towergeom->get_phi();
1372                         //unsigned int keyInfo  = TowerInfoDefs::encode_emcal(ieta,iphi);
1373                         //unsigned int ichannel = TowerInfoDefs::decode_emcal(keyInfo);
1374                         //TowerInfo* towinfo    = towersEM->get_tower_at_channel(ichannel);
1375                         double esum =  0.;
1376                         for (int idphi=-3;idphi<=3;idphi++){
1377                             for (int ideta=-3;ideta<=3;ideta++){
1378                                 int ietatow = ieta + ideta;
1379                                 int iphitow = iphi + idphi;
1380                                 if (iphitow<   0) iphitow+=256; 
1381                                 if (iphitow>=256) iphitow-=256; 
1382                                 double etow = -9.;
1383                                 if (ietatow>=0&&ietatow<96&&iphitow>=0&&iphitow<256){   // >>>---> ieta,iphi in hand for this state! 
1384                                     unsigned int keytow      = TowerInfoDefs::encode_emcal(ietatow,iphitow);    // ieta,iphi->key
1385                                     unsigned int ichan       = TowerInfoDefs::decode_emcal(keytow);             //       key->chan
1386                                     TowerInfo* towi          = towersEM->get_tower_at_channel(ichan);           //      chan->info
1387                                     etow                     = towi->get_energy();
1388                                     esum                    += etow;
1389                                 }   
1390                                 fTegridEM[fNTracks][3+ideta][3+idphi]    = etow;
1391                             }   // end ideta
1392                         }   // end idphi
1393                         fTegridEMsum7x7[fNTracks]   = esum;
1394                         fTeopEM[fNTracks]           = esum/track->get_p();
1395                         if (SHOWGRIDS && track->get_p()>5.){
1396                             cout<<fEvtSeq<<"\tEM "<<track->get_id()<<" "<<track->get_p()<<" "<<fTegridEMsum7x7[fNTracks]<<endl;
1397                             for (int idphi=-3;idphi<=3;idphi++){
1398                                 for (int ideta=-3;ideta<=3;ideta++){
1399                                     cout<<FIXFLOAT(fTegridEM[fNTracks][3+ideta][3+idphi])<<" ";
1400                                 }   cout<<endl;
1401                             }       cout<<endl;
1402                         }   // end DEBUG check
1403                         //
1404                         if (CALODETAIL && ievtSeen<MAXEVTDETAIL){
1405                             hPTOT_etaDETem[ievtSeen]    ->Fill( fTstateEM_eta[fNTracks], track->get_p() );
1406                             hPTOT_phiDETem[ievtSeen]    ->Fill( fTstateEM_phi[fNTracks], track->get_p() );
1407                         }
1408                         //
1409                     } else {
1410                         fTstateEM_x[fNTracks]   =  0.;
1411                         fTstateEM_y[fNTracks]   =  0.;
1412                         fTstateEM_z[fNTracks]   =  0.;
1413                         fTstateEM_eta[fNTracks] = -9.;
1414                         fTstateEM_phi[fNTracks] = -9.;
1415                         for (int idphi=-3;idphi<=3;idphi++){
1416                             for (int ideta=-3;ideta<=3;ideta++){
1417                                 fTegridEM[fNTracks][3+ideta][3+idphi]   = -9;
1418                             }
1419                         }
1420                         fTeopEM[fNTracks]   = -9.;
1421                     }   //  end ieta,iphi check
1422                 } else {
1423                     fTstateEM_x[fNTracks]   =  0.;
1424                     fTstateEM_y[fNTracks]   =  0.;
1425                     fTstateEM_z[fNTracks]   =  0.;
1426                     fTstateEM_eta[fNTracks] = -9.;
1427                     fTstateEM_phi[fNTracks] = -9.;
1428                     for (int idphi=-3;idphi<=3;idphi++){
1429                         for (int ideta=-3;ideta<=3;ideta++){
1430                             fTegridEM[fNTracks][3+ideta][3+idphi]   = -9;
1431                         }
1432                     }
1433                     fTeopEM[fNTracks]   = -9.;
1434                 }   //  end thisState
1435                 //
1436                 //---- take track to ihcal depth
1437                 thisStateI  = track->get_state(caloRadiusInnerIHCal);
1438                 thisStateO  = track->get_state(caloRadiusOuterIHCal);
1439                 if ( thisStateI && !thisStateO){ 
1440                     stx = thisStateI->get_x();
1441                     sty = thisStateI->get_y();
1442                     stz = thisStateI->get_z();
1443                 } else if ( !thisStateI && thisStateO){ 
1444                     stx = thisStateO->get_x();
1445                     sty = thisStateO->get_y();
1446                     stz = thisStateO->get_z();
1447                 } else if ( thisStateI && thisStateO){ 
1448                     stx = (thisStateI->get_x() + thisStateO->get_x())/2.;
1449                     sty = (thisStateI->get_y() + thisStateO->get_y())/2.;
1450                     stz = (thisStateI->get_z() + thisStateO->get_z())/2.;
1451                 } 
1452                 if (thisStateI||thisStateO){
1453                     fTstateIH_x[fNTracks]   = stx;
1454                     fTstateIH_y[fNTracks]   = sty;
1455                     fTstateIH_z[fNTracks]   = stz;
1456                     fTstateIH_eta[fNTracks] = asinh(stz/sqrt(stx*stx + sty*sty));
1457                     fTstateIH_phi[fNTracks] = atan2(sty,stx);
1458                     int ieta                = GetIETA(SvtxTrack::HCALIN, fTstateIH_eta[fNTracks]);
1459                     int iphi                = GetIPHI(SvtxTrack::HCALIN, fTstateIH_phi[fNTracks]);
1460                     if (ieta>=0&&ieta<24&&iphi>=0&&iphi<64){    // >>>---> ieta,iphi in hand for this state! 
1461                         //RawTowerDefs::keytype keyRaw  = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
1462                         //RawTowerGeom* towergeom       = IHCalGeo->get_tower_geometry(keyRaw);
1463                         //double toweta                 = towergeom->get_eta();
1464                         //double towphi                 = towergeom->get_phi();
1465                         //unsigned int keyInfo  = TowerInfoDefs::encode_hcal(ieta,iphi);
1466                         //unsigned int ichannel = TowerInfoDefs::decode_hcal(keyInfo);
1467                         //TowerInfo* towinfo    = towersIH->get_tower_at_channel(ichannel);
1468                         double esum =  0.;
1469                         for (int idphi=-3;idphi<=3;idphi++){
1470                             for (int ideta=-3;ideta<=3;ideta++){
1471                                 int ietatow = ieta + ideta;
1472                                 int iphitow = iphi + idphi;
1473                                 if (iphitow<   0) iphitow+= 64; 
1474                                 if (iphitow>= 64) iphitow-= 64; 
1475                                 double etow = -9.;
1476                                 if (ietatow>=0&&ietatow<24&&iphitow>=0&&iphitow<64){    // >>>---> ieta,iphi in hand for this state! 
1477                                     unsigned int keytow      = TowerInfoDefs::encode_hcal(ietatow,iphitow);     // ieta,iphi->key
1478                                     unsigned int ichan       = TowerInfoDefs::decode_hcal(keytow);              //       key->chan
1479                                     TowerInfo* towi          = towersIH->get_tower_at_channel(ichan);           //      chan->info
1480                                     etow                     = towi->get_energy();
1481                                     esum                    += etow;
1482                                 }   
1483                                 fTegridIH[fNTracks][3+ideta][3+idphi]    = etow;
1484                             }   // end ideta
1485                         }   // end idphi
1486                         fTegridIHsum7x7[fNTracks]   = esum;
1487                         fTeopIH[fNTracks]           = esum/track->get_p();
1488                         if (SHOWGRIDS && track->get_p()>5.){
1489                             cout<<fEvtSeq<<"\tIH "<<track->get_id()<<" "<<track->get_p()<<" "<<fTegridIHsum7x7[fNTracks]<<endl;
1490                             for (int idphi=-3;idphi<=3;idphi++){
1491                                 for (int ideta=-3;ideta<=3;ideta++){
1492                                     cout<<FIXFLOAT(fTegridIH[fNTracks][3+ideta][3+idphi])<<" ";
1493                                 }   cout<<endl;
1494                             }       cout<<endl;
1495                         }   // end DEBUG check
1496                         //
1497                         if (CALODETAIL && ievtSeen<MAXEVTDETAIL){
1498                             hPTOT_etaDETih[ievtSeen]    ->Fill( fTstateIH_eta[fNTracks], track->get_p() );
1499                             hPTOT_phiDETih[ievtSeen]    ->Fill( fTstateIH_phi[fNTracks], track->get_p() );
1500                         }
1501                         //
1502                     } else {
1503                         fTstateIH_x[fNTracks]   =  0.;
1504                         fTstateIH_y[fNTracks]   =  0.;
1505                         fTstateIH_z[fNTracks]   =  0.;
1506                         fTstateIH_eta[fNTracks] = -9.;
1507                         fTstateIH_phi[fNTracks] = -9.;
1508                         for (int idphi=-3;idphi<=3;idphi++){
1509                             for (int ideta=-3;ideta<=3;ideta++){
1510                                 fTegridIH[fNTracks][3+ideta][3+idphi]   = -9;
1511                             }
1512                         }
1513                         fTeopIH[fNTracks]   = -9.;
1514                     }   //  end ieta,iphi check
1515                 } else {
1516                     fTstateIH_x[fNTracks]   =  0.;
1517                     fTstateIH_y[fNTracks]   =  0.;
1518                     fTstateIH_z[fNTracks]   =  0.;
1519                     fTstateIH_eta[fNTracks] = -9.;
1520                     fTstateIH_phi[fNTracks] = -9.;
1521                     for (int idphi=-3;idphi<=3;idphi++){
1522                         for (int ideta=-3;ideta<=3;ideta++){
1523                             fTegridIH[fNTracks][3+ideta][3+idphi]   = -9;
1524                         }
1525                     }
1526                     fTeopIH[fNTracks]   = -9.;
1527                 }   //  end thisState
1528                 //
1529                 //---- take track to ohcal depth
1530                 thisStateI  = track->get_state(caloRadiusInnerOHCal);
1531                 thisStateO  = track->get_state(caloRadiusOuterOHCal);
1532                 if ( thisStateI && !thisStateO){ 
1533                     stx = thisStateI->get_x();
1534                     sty = thisStateI->get_y();
1535                     stz = thisStateI->get_z();
1536                 } else if ( !thisStateI && thisStateO){ 
1537                     stx = thisStateO->get_x();
1538                     sty = thisStateO->get_y();
1539                     stz = thisStateO->get_z();
1540                 } else if ( thisStateI && thisStateO){ 
1541                     stx = (thisStateI->get_x() + thisStateO->get_x())/2.;
1542                     sty = (thisStateI->get_y() + thisStateO->get_y())/2.;
1543                     stz = (thisStateI->get_z() + thisStateO->get_z())/2.;
1544                 } 
1545                 if (thisStateI||thisStateO){
1546                     fTstateOH_x[fNTracks]   = stx;
1547                     fTstateOH_y[fNTracks]   = sty;
1548                     fTstateOH_z[fNTracks]   = stz;
1549                     fTstateOH_eta[fNTracks] = asinh(stz/sqrt(stx*stx + sty*sty));
1550                     fTstateOH_phi[fNTracks] = atan2(sty,stx);
1551                     int ieta                = GetIETA(SvtxTrack::HCALOUT, fTstateOH_eta[fNTracks]);
1552                     int iphi                = GetIPHI(SvtxTrack::HCALOUT, fTstateOH_phi[fNTracks]);
1553                     if (ieta>=0&&ieta<24&&iphi>=0&&iphi<64){    // >>>---> ieta,iphi in hand for this state! 
1554                         //RawTowerDefs::keytype keyRaw  = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
1555                         //RawTowerGeom* towergeom       = OHCalGeo->get_tower_geometry(keyRaw);
1556                         //double toweta                 = towergeom->get_eta();
1557                         //double towphi                 = towergeom->get_phi();
1558                         //unsigned int keyInfo  = TowerInfoDefs::encode_hcal(ieta,iphi);
1559                         //unsigned int ichannel = TowerInfoDefs::decode_hcal(keyInfo);
1560                         //TowerInfo* towinfo    = towersOH->get_tower_at_channel(ichannel);
1561                         double esum =  0.;
1562                         for (int idphi=-3;idphi<=3;idphi++){
1563                             for (int ideta=-3;ideta<=3;ideta++){
1564                                 int ietatow = ieta + ideta;
1565                                 int iphitow = iphi + idphi;
1566                                 if (iphitow<   0) iphitow+= 64; 
1567                                 if (iphitow>= 64) iphitow-= 64; 
1568                                 double etow = -9.;
1569                                 if (ietatow>=0&&ietatow<24&&iphitow>=0&&iphitow<64){    // >>>---> ieta,iphi in hand for this state! 
1570                                     unsigned int keytow      = TowerInfoDefs::encode_hcal(ietatow,iphitow);     // ieta,iphi->key
1571                                     unsigned int ichan       = TowerInfoDefs::decode_hcal(keytow);              //       key->chan
1572                                     TowerInfo* towi          = towersOH->get_tower_at_channel(ichan);           //      chan->info
1573                                     etow                     = towi->get_energy();
1574                                     esum                    += etow;
1575                                 }   
1576                                 fTegridOH[fNTracks][3+ideta][3+idphi]   = etow;
1577                             }   // end ideta
1578                         }   // end idphi
1579                         fTegridOHsum7x7[fNTracks]   = esum;
1580                         fTeopOH[fNTracks]           = esum/track->get_p();
1581                         if (SHOWGRIDS && track->get_p()>5.){
1582                             cout<<fEvtSeq<<"\tOH "<<track->get_id()<<" "<<track->get_p()<<" "<<fTegridOHsum7x7[fNTracks]<<endl;
1583                             for (int idphi=-3;idphi<=3;idphi++){
1584                                 for (int ideta=-3;ideta<=3;ideta++){
1585                                     cout<<FIXFLOAT(fTegridOH[fNTracks][3+ideta][3+idphi])<<" ";
1586                                 }   cout<<endl;
1587                             }       cout<<endl;
1588                         }   // end DEBUG check
1589                         //
1590                         if (CALODETAIL && ievtSeen<MAXEVTDETAIL){
1591                             hPTOT_etaDEToh[ievtSeen]    ->Fill( fTstateOH_eta[fNTracks], track->get_p() );
1592                             hPTOT_phiDEToh[ievtSeen]    ->Fill( fTstateOH_phi[fNTracks], track->get_p() );
1593                         }
1594                         //
1595                     } else {
1596                         fTstateOH_x[fNTracks]   =  0.;
1597                         fTstateOH_y[fNTracks]   =  0.;
1598                         fTstateOH_z[fNTracks]   =  0.;
1599                         fTstateOH_eta[fNTracks] = -9.;
1600                         fTstateOH_phi[fNTracks] = -9.;
1601                         for (int idphi=-3;idphi<=3;idphi++){
1602                             for (int ideta=-3;ideta<=3;ideta++){
1603                                 fTegridOH[fNTracks][3+ideta][3+idphi]   = -9;
1604                             }
1605                         }
1606                         fTeopOH[fNTracks]   = -9.;
1607                     }   //  end ieta,iphi check
1608                 } else {
1609                     fTstateOH_x[fNTracks]   =  0.;
1610                     fTstateOH_y[fNTracks]   =  0.;
1611                     fTstateOH_z[fNTracks]   =  0.;
1612                     fTstateOH_eta[fNTracks] = -9.;
1613                     fTstateOH_phi[fNTracks] = -9.;
1614                     for (int idphi=-3;idphi<=3;idphi++){
1615                         for (int ideta=-3;ideta<=3;ideta++){
1616                             fTegridOH[fNTracks][3+ideta][3+idphi]   = -9;
1617                         }
1618                     }
1619                     fTeopOH[fNTracks]   = -9.;
1620                 }   //  end thisState
1621                 } else {    // !fHasCaloGeo or !fAddTowersToTrackTree -- calo info unavailable or not wanted, leave calo-layer-around-track fields at "no info" defaults
1622                     fTstateEM_eta[fNTracks] = fTstateIH_eta[fNTracks] = fTstateOH_eta[fNTracks] = -9.;
1623                     fTstateEM_phi[fNTracks] = fTstateIH_phi[fNTracks] = fTstateOH_phi[fNTracks] = -9.;
1624                     fTegridEMsum7x7[fNTracks] = fTegridIHsum7x7[fNTracks] = fTegridOHsum7x7[fNTracks] = 0.;
1625                     fTeopEM[fNTracks] = fTeopIH[fNTracks] = fTeopOH[fNTracks] = -9.;
1626                     for (int idphi=-3;idphi<=3;idphi++){
1627                         for (int ideta=-3;ideta<=3;ideta++){
1628                             fTegridEM[fNTracks][3+ideta][3+idphi] = -9;
1629                             fTegridIH[fNTracks][3+ideta][3+idphi] = -9;
1630                             fTegridOH[fNTracks][3+ideta][3+idphi] = -9;
1631                         }
1632                     }
1633                 }   // end fHasCaloGeo
1634                 {   // eopHC/eopCA: sum of per-layer eop values (same denominator
1635                     // track->get_p(), so addition is exact) - a hadronic shower does
1636                     // not respect the IH/OH boundary, so the combined hadronic
1637                     // (HC=IH+OH) and total-calorimeter (CA=EM+IH+OH) ratios are more
1638                     // physically meaningful than the individual IH/OH ratios alone.
1639                     bool hasEM = (fTeopEM[fNTracks] > -8.5);
1640                     bool hasIH = (fTeopIH[fNTracks] > -8.5);
1641                     bool hasOH = (fTeopOH[fNTracks] > -8.5);
1642                     fTeopHC[fNTracks] = (hasIH||hasOH) ? (hasIH?fTeopIH[fNTracks]:0.) + (hasOH?fTeopOH[fNTracks]:0.) : -9.;
1643                     fTeopCA[fNTracks] = (hasEM||hasIH||hasOH) ? (hasEM?fTeopEM[fNTracks]:0.) + (hasIH?fTeopIH[fNTracks]:0.) + (hasOH?fTeopOH[fNTracks]:0.) : -9.;
1644                 }
1645                 if (DEBUG2) cout<<fNTracks<<"\t CALO "<<fTegridEMsum7x7[fNTracks]<<" "<<fTegridIHsum7x7[fNTracks]<<" "<<fTegridOHsum7x7[fNTracks]<<" "<<endl;
1646                 //
1647                 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1648                 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1649                 //              
1650                 //
1651                 //---- loop through v0 crossing info and find any daughters
1652                 if (DEBUG2) cout<<"TrackLoop: loop over v0s and find daughters..."<<endl;
1653                 int thispid     =  0;
1654                 int thisindv0   = -1;
1655                 if (v0found){
1656                     for (int iv0=0;iv0<fTnv0;iv0++){
1657                         int itf         = fV0FinalToTF[iv0];    // fTv0*[iv0] came from trigger-frame slot itf - V0Crossing/V0pid/V0Track1id/V0Track2id below are STILL in trigger-frame space, NOT final space, so they must be read at itf, not iv0!
1658                         int thisv0xing  = V0Crossing[itf];
1659                         if (thisv0xing == m_track_crossing      // this v0 is in the same crossing as this track!
1660                          && thispid    == 0 ){                  // this line prevents assigning this track to multiple V0s!
1661                             if (trackid== V0Track1id[itf]){
1662                                 fTv0indtr1[iv0] = fNTracks;     // set index to this track(1) in v0 tree!
1663                             }
1664                             if (trackid== V0Track2id[itf]){
1665                                 fTv0indtr2[iv0] = fNTracks;     // set index to this track(2) in v0 tree!
1666                             }
1667                             if (trackid==V0Track1id[itf] && V0pid[itf]== Kshort_PDGID){ thispid = 1; }  //trk1 pi+
1668                             if (trackid==V0Track2id[itf] && V0pid[itf]== Kshort_PDGID){ thispid = 1; }  //trk2 pi-
1669                             if (trackid==V0Track1id[itf] && V0pid[itf]== Lambda_PDGID){ thispid = 3; }  //trk1 proton
1670                             if (trackid==V0Track2id[itf] && V0pid[itf]== Lambda_PDGID){ thispid = 1; }  //trk2 pi-
1671                             if (trackid==V0Track1id[itf] && V0pid[itf]==ALambda_PDGID){ thispid = 1; }  //trk1 pi+
1672                             if (trackid==V0Track2id[itf] && V0pid[itf]==ALambda_PDGID){ thispid = 3; }  //trk2 antiproton
1673                             //if (DEBUG) cout   <<"V0 LOOP .. iv0="<<iv0<<"\t v0xing: "<<thisv0xing<<" "<<m_track_crossing
1674                             //              <<"\t trkid= "<<trackid<<"\t fNTracks= "<<fNTracks
1675                             //              <<"\t pid= "<<thispid<<"\t thisindv0= "<<thisindv0
1676                             //              <<endl;
1677                             if (thispid>0){     //---- set index to this V0
1678                                 thisindv0   = iv0;
1679                                 if (DEBUG) cout<<"V0-Track Match: iv0="<<iv0
1680                                     <<"\t trackid= "<<trackid
1681                                     <<"\t V0trackid= "<<V0Track1id[itf]<<" "<<V0Track2id[itf]
1682                                     <<"\t trkxing= "<<m_track_crossing
1683                                     //<<"\t v0xing= "<<thisv0xing
1684                                     <<"\t pid,chg= "<<thispid<<" "<<chg
1685                                     <<"\t indv0= "<<thisindv0
1686                                     <<endl;
1687                             }
1688                         }
1689                     }
1690                 }   // end v0found...
1691                 //
1692                 if (thispid>0 && !KEEPTRK) KEEPTRK = true;  //!!!!! force keep if v0 daughter..
1693                 //
1694                 //if (DEBUG) cout<<"TrackLoop: done loop over v0s and find daughters..."<<endl;
1695                 //
1696                 if (KEEPTRK){
1697                     //
1698                     if (fNTracks>=MAXNTRACKS-1){
1699                         std::cout<<"Collect::process_event -- Arrays for tree are full!  Increase MAXNTRACKS."<<std::endl;
1700                         //exit(0);
1701                     }
1702                     if (fNTracks<MAXNTRACKS-1){
1703                         fTpid[fNTracks]         = thispid;
1704                         fTindv0[fNTracks]       = thisindv0;
1705                         fTprimary[fNTracks]     = vertexFlag;
1706                         fTquality[fNTracks]     = quality;
1707                         fTchisq[fNTracks]       = chisq;
1708                         fTchg[fNTracks]         = chg;
1709                         fTpx[fNTracks]          = px;
1710                         fTpy[fNTracks]          = py;
1711                         fTpz[fNTracks]          = pz;
1712                         fTptot[fNTracks]        = ptot;
1713                         fTeta[fNTracks]         = eta;
1714                         fTphi[fNTracks]         = phi;
1715                         fTpt[fNTracks]          = pt;
1716                         fTseedeta[fNTracks]     = seedeta;
1717                         fTseedphi[fNTracks]     = seedphi;
1718                         fTseedpt[fNTracks]      = seedpt;
1719                         fTnhitstpc[fNTracks]    = ntpc;
1720                         fTnhitsmvt[fNTracks]    = nmaps;
1721                         fTnhitsint[fNTracks]    = nintt;
1722                         fTdcaxy[fNTracks]       = dca3d_xy;
1723                         fTdcaz[fNTracks]        = dca3d_z;
1724                         fTdedx70n[fNTracks]     = (int)dedx70n;
1725                         fTdedx70s[fNTracks]     = dedx70s;
1726                         fTdedxKFP[fNTracks]     = dedxKFP;
1727                         fTlayermask[fNTracks]   = layermask;
1728                         for (int isl=0;isl<7;isl++){ fTsiclukey[fNTracks][isl] = siclukey_local[isl]; }
1729                         for (int itl=0;itl<48;itl++){
1730                             fTtpcsector[fNTracks][itl] = tpcsector_local[itl];
1731                             fTtpcarclen[fNTracks][itl] = tpcarclen_local[itl];
1732                             fTtpcz[fNTracks][itl]      = tpcz_local[itl];
1733                         }
1734                         if (DEBUG){
1735                             std::cout<<trackid<<"\t chg="<<chg<<"\t"
1736                                     //<<deltapt<<"\t"
1737                                     <<vertexFlag<<" zv="<<fTvtxz<<" "<<fTvtxzGlob<<"\t "
1738                                     //<<px<<" "<<py<<" "<<pz<<" "<<ptot<<"\t"
1739                                     <<"\t eta="<<eta<<" ("<<seedeta<<")"
1740                                     <<"\t phi="<<phi<<" ("<<seedphi<<")"
1741                                     <<"\t pt=" <<pt <<" ("<<seedpt <<")"
1742                                     //<<vx<<" "<<vy<<" "<<vz<<"\t"
1743                                     //<<pcax<<" "<<pcay<<" "<<pcaz<<"\t"
1744                                     <<"\t dedx="<<dedx70n<<" "<<dedx70s
1745                                     <<"\t Ncl="<<nmaps<<" "<<nintt<<" "<<ntpc
1746                                     <<"\t CALO: "<<fTegridEMsum7x7[fNTracks]<<" "<<fTegridIHsum7x7[fNTracks]<<" "<<fTegridOHsum7x7[fNTracks]
1747                                     <<std::endl;
1748                         }
1749                         // DIAGTEMP-20260920: commit this track's cluster list under its fNTracks index
1750                         diagClustersByTrack20260920.push_back(thisTrackClusters20260920);
1751                         for (const auto &ci20260920 : thisTrackClusters20260920){
1752                             diagKeyOwners20260920[ci20260920.key].push_back(fNTracks);
1753                         }
1754                         ++fNTracks;
1755                         //
1756                     } else {
1757                         std::cout<<"Collect::process_event -- Track "<<trackid<<" could not be saved... "<<std::endl;
1758                     }   // end array protection on fNTracks...
1759                 }   // end KEEPTRK
1760                 //
1761             }   // end loop over tracks...
1762             //
1763             if (fNTracks>0) FlagSplitTracks();
1764             if (fNTracks>0) DiagMergedTracksExplore20260920(fNTracks, fTlayermask, diagClustersByTrack20260920, diagKeyOwners20260920); // DIAGTEMP-20260920
1765             //
1766             if (fNTracks){
1767                 ++nverticeskept;
1768                 outTree->Fill();
1769             }
1770             //
1771             if (fNTracks>0){
1772                 for (int i=0;i<=fNTracks;i++){
1773                     fTprimary[i]    = 0;        // reset track/vertex flag
1774                 }
1775             }
1776             //
1777         }   // end loop over vertices...
1778     }   // end svtxvertexmap check...
1779     //
1780     if (Xings.size()>=2){
1781         std::stable_sort(Xings.begin(), Xings.end());
1782         for (unsigned int ix=0;ix<Xings.size()-1;ix++){
1783             double dXing    = Xings[ix+1] - Xings[ix];
1784             hdXing          ->Fill(dXing);
1785         }
1786     }
1787     //
1788     //
1789     //std::cout<<"Event "<<ievtSeen-1<<"\t nvertices="<<nverticeskept<<endl;
1790     //
1791     //
1792     return Fun4AllReturnCodes::EVENT_OK;
1793     //
1794 }
1795 
1796 
1797 //------------------------------------------------------------------------
1798 //
1799 int Collect::GetKFP( PHCompositeNode *topNode, const int parentpdgid, const std::string &containername, const std::string &trackmapname ){
1800     //
1801     int linecount       = 0;
1802     int initialnv0tf    = nv0tf;            // V0 counter INCLUDING v0's found in previous calls to this routine!!!
1803     int thisv0id        = initialnv0tf;     // V0 counter INCLUDING v0's found in previous calls to this routine!!!
1804     //if (DEBUG) cout<<"Collect::GetKFP START ..  nv0tf= "<<initialnv0tf<<endl;
1805     //
1806     KFParticle_Container* kfpContainerKS = findNode::getClass<KFParticle_Container>(topNode, containername);
1807     if (kfpContainerKS){
1808         //
1809         int pdgidA=0,pdgidB=0,trkidA=0,trkidB=0;
1810         double posxA=0,posyA=0,poszA=0,posxB=0,posyB=0,poszB=0;
1811         //
1812         for (KFParticle_Container::Iter kfp_iter = kfpContainerKS->begin(); kfp_iter != kfpContainerKS->end(); ++kfp_iter){
1813             KFParticle *thisKFpart      = kfp_iter->second;
1814             if (DEBUG) cout<<"KScontainer .. "<<thisKFpart->GetPDG()<<" "<<thisKFpart->Id()<<" "<<thisv0id<<endl;
1815             //
1816             if        (linecount%3==0){
1817                 pdgidA = thisKFpart->GetPDG();
1818                 trkidA = thisKFpart->Id();
1819                 posxA  = thisKFpart->GetX();
1820                 posyA  = thisKFpart->GetY();
1821                 poszA  = thisKFpart->GetZ();
1822             } else if (linecount%3==1){
1823                 pdgidB = thisKFpart->GetPDG();
1824                 trkidB = thisKFpart->Id();
1825                 posxB  = thisKFpart->GetX();
1826                 posyB  = thisKFpart->GetY();
1827                 poszB  = thisKFpart->GetZ();
1828             } else if (linecount%3==2){
1829                 if (thisv0id < MAXNV0S){
1830                     V0pid[thisv0id]         = parentpdgid;
1831                     V0DecayX[thisv0id]      = 0.5*(posxA+posxB);
1832                     V0DecayY[thisv0id]      = 0.5*(posyA+posyB);
1833                     V0DecayZ[thisv0id]      = 0.5*(poszA+poszB);
1834                     TParticlePDG* particleA = TDatabasePDG::Instance()->GetParticle(pdgidA);
1835                     if (particleA->Charge()>0.){            // particleA is positive so it's particle 1
1836                         V0Track1id[thisv0id]    = trkidA;   //      set A as 1 because it's positive
1837                         V0Track2id[thisv0id]    = trkidB;   //      set B as 2 because it's negative
1838                     } else if (particleA->Charge()<0.){     // particleB is positive so it's particle 1
1839                         V0Track1id[thisv0id]    = trkidB;   //      set B as 1 because it's positive
1840                         V0Track2id[thisv0id]    = trkidA;   //      set A as 2 because it's negative
1841                     }                                       //
1842                     if (DEBUG) cout<<"KScont view .. "<<pdgidA<<" "<<pdgidB<<"\t "<<trkidA<<" "<<trkidB
1843                                 <<"\t v0id: "<<thisv0id<<" trkids: "<<V0Track1id[thisv0id]<<" "<<V0Track2id[thisv0id]<<endl;
1844                 } else {
1845                     std::cout<<"Collect::GetKFP -- Pass 1: too many V0s for MAXNV0S="<<MAXNV0S<<", dropping v0id="<<thisv0id<<std::endl;
1846                 }
1847                 ++thisv0id;
1848             }
1849             ++linecount;
1850         }
1851         //
1852         if (linecount==0){
1853             return Fun4AllReturnCodes::EVENT_OK;    // nothing to do, so leave. 
1854         }
1855         //
1856         //----- get v0s from containers (one container, many crossings!) TF=TriggerFrame
1857         for (auto &iter : *kfpContainerKS){
1858             if (iter.second->GetPDG() == parentpdgid){
1859                 KFParticle* thisKFpart  = iter.second;
1860                 double tmass    = thisKFpart->GetMass();
1861                 double tpx      = thisKFpart->GetPx();
1862                 double tpy      = thisKFpart->GetPy();
1863                 double tpz      = thisKFpart->GetPz();
1864                 if (nv0tf < MAXNV0S-1 ){
1865                     fTv0TFuseit[nv0tf]      = true;
1866                     fTv0TFpid[nv0tf]        = parentpdgid;
1867                     fTv0TFx[nv0tf]          = thisKFpart->GetX();
1868                     fTv0TFy[nv0tf]          = thisKFpart->GetY();
1869                     fTv0TFz[nv0tf]          = thisKFpart->GetZ();
1870                     fTv0TFpx[nv0tf]         = tpx;
1871                     fTv0TFpy[nv0tf]         = tpy;
1872                     fTv0TFpz[nv0tf]         = tpz;
1873                     fTv0TFpt[nv0tf]         = sqrt(tpx*tpx + tpy*tpy);
1874                     fTv0TFptot[nv0tf]       = sqrt(tpz*tpz + fTv0TFpt[nv0tf]*fTv0TFpt[nv0tf]);
1875                     fTv0TFeta[nv0tf]        = thisKFpart->GetEta();
1876                     fTv0TFphi[nv0tf]        = thisKFpart->GetPhi();
1877                     fTv0TFene[nv0tf]        = thisKFpart->GetE();
1878                     fTv0TFmass[nv0tf]       = tmass;
1879                     fTv0TFctau[nv0tf]       = thisKFpart->GetLifeTime();
1880                     fTv0TFdecaylen[nv0tf]   = thisKFpart->GetDecayLength();
1881                     double chi2ndf          = thisKFpart->GetChi2();
1882                     if (thisKFpart->GetNDF()){ chi2ndf /= thisKFpart->GetNDF(); }else{ chi2ndf = -1; }
1883                     fTv0TFchi2ndf[nv0tf]    = chi2ndf;
1884                     fTv0TFindtr1[nv0tf]     = -1;       // set in later code! 
1885                     fTv0TFindtr2[nv0tf]     = -1;       // set in later code! 
1886                     ++nv0tf;
1887                 } else {
1888                     std::cout<<"Collect .. TOO MANY V0s!!!!!!!   fix this. "<<std::endl;
1889                     return Fun4AllReturnCodes::ABORTPROCESSING;
1890                 }
1891                 if (DEBUG)
1892                 cout<<"KFP v0 .. "<<nv0tf           <<"\t"
1893                     <<"mass=" <<thisKFpart->GetMass()           <<"\t"
1894                     <<"pdgid="<<thisKFpart->GetPDG()            <<"\t"
1895                     <<"eta="  <<thisKFpart->GetEta()            <<" "
1896                     <<"phi="  <<thisKFpart->GetPhi()            <<" "
1897                     <<"pt="   <<thisKFpart->GetPt()             <<"\t"
1898                     <<"X="    <<thisKFpart->GetX()              <<" "
1899                     <<"Y="    <<thisKFpart->GetY()              <<" "
1900                     <<"Z="    <<thisKFpart->GetZ()              <<"\t"
1901                     <<"ctau=" <<thisKFpart->GetLifeTime()       <<"\t"
1902                     <<"DL="   <<thisKFpart->GetDecayLength()    <<"\t"
1903                     <<"chi2=" <<thisKFpart->GetChi2()<<" "<<thisKFpart->GetNDF() <<"\t"
1904                     <<endl;
1905             }
1906         }
1907 //!!        for (KFParticle_Container::Iter kfp_iter = kfpContainerKS->begin(); kfp_iter != kfpContainerKS->end(); ++kfp_iter){
1908 //!!            KFParticle *thisKFpart      = kfp_iter->second;
1909 //!!            unsigned int thisKFtrackid  = thisKFpart->Id();
1910 //!!            int pdgid                   = thisKFpart->GetPDG();
1911 //!!            //cout<<"TEST "<<thisKFtrackid<<" "<<pdgid<<" "<<thisKFpart->GetEta()<<endl;
1912 //!!        }       
1913     }   // end KFPcontainer KS......
1914     //
1915     SvtxTrackMap* KFP_trackMapKS = findNode::getClass<SvtxTrackMap>(topNode, trackmapname);
1916     if(KFP_trackMapKS){
1917         // Look up each v0's crossing directly by its daughter track ID (already
1918         // correctly assigned above, in the KFParticle_Container pass) instead of
1919         // scanning this map positionally in assumed [daughter,daughter,mother]
1920         // triples. KFP_trackMapKS is sorted by key, and KFParticle_DST::
1921         // fillParticleNode_Track gives mother/intermediate pseudo-tracks huge
1922         // keys near UINT_MAX (a decrementing resonanceCounter) so they don't
1923         // collide with real track IDs - so with 2+ candidates of this species
1924         // in one event, every real daughter (small key) sorts before every
1925         // mother placeholder (huge key), breaking the triple-per-candidate
1926         // assumption the old position-based scan relied on. That silently
1927         // left V0Crossing[] holding a mother-placeholder's bogus crossing
1928         // (typically 32767, a sentinel) for every candidate past the first,
1929         // which then could never match any real vertex's crossing and was
1930         // dropped from Collect's v0 tree with no trace - see the
1931         // collect_v0_track_index_bug memory for how this was found.
1932         for (int iv0 = initialnv0tf; iv0 < nv0tf && iv0 < MAXNV0S; ++iv0){
1933             SvtxTrack* daughter1 = KFP_trackMapKS->get(V0Track1id[iv0]);
1934             if (daughter1){
1935                 V0Crossing[iv0] = daughter1->get_crossing();
1936             }
1937             if (DEBUG){
1938                 std::cout<<"KFP trk .. indexV0: "<<iv0
1939                         <<"\t trackids: "<<V0Track1id[iv0]<<" "<<V0Track2id[iv0]
1940                         <<"\t xing="<<V0Crossing[iv0]
1941                         <<std::endl;
1942             }
1943         }
1944     }
1945     //
1946     if (DEBUG) cout<<"Collect::GetKFP END ..  running nv0tf= "<<initialnv0tf<<"\t this nv0tf= "<<nv0tf-initialnv0tf<<"\t ?: "<<thisv0id<<endl;
1947     return nv0tf-initialnv0tf;
1948     //
1949 }   // end GetKFP
1950 
1951 
1952 
1953 //----------------------------------------------------------------------
1954 //
1955 //  modified WJL 
1956 //  now gives dca values even if the given track is /not attached/ to the given vertex!
1957 //
1958 void Collect::calc3DDCA(SvtxTrack* track, SvtxVertexMap* vertexmap, unsigned int vtxid,
1959                double& dca3d_xy, double& dca3d_xy_error,
1960                double& dca3d_z , double& dca3d_z_error )
1961 {
1962   Acts::Vector3 pos(track->get_x(),
1963             track->get_y(),
1964             track->get_z());
1965   Acts::Vector3 mom(track->get_px(),
1966             track->get_py(),
1967             track->get_pz());
1968 //WJL  auto vtxid = track->get_vertex_id();
1969   auto svtxVertex = vertexmap->get(vtxid);
1970   Acts::Vector3 vertex(svtxVertex->get_x(),
1971                svtxVertex->get_y(),
1972                svtxVertex->get_z());
1973   pos -= vertex;
1974   Acts::ActsMatrix<3,3> posCov;
1975   for(int i = 0; i < 3; ++i){
1976       for(int j = 0; j < 3; ++j){
1977           posCov(i, j) = track->get_error(i, j);
1978       } 
1979   }
1980   Acts::Vector3 r = mom.cross(Acts::Vector3(0.,0.,1.));
1981   double phi = atan2(r(1), r(0));
1982   Acts::RotationMatrix3 rot;
1983   Acts::RotationMatrix3 rot_T;
1984   rot(0,0) = cos(phi);
1985   rot(0,1) = -sin(phi);
1986   rot(0,2) = 0;
1987   rot(1,0) = sin(phi);
1988   rot(1,1) = cos(phi);
1989   rot(1,2) = 0;
1990   rot(2,0) = 0;
1991   rot(2,1) = 0;
1992   rot(2,2) = 1;
1993   rot_T = rot.transpose();
1994   Acts::Vector3 pos_R = rot * pos;
1995   Acts::ActsMatrix<3,3> rotCov = rot * posCov * rot_T;
1996   dca3d_xy= pos_R(0);
1997   dca3d_z = pos_R(2);
1998   dca3d_xy_error = sqrt(rotCov(0,0));
1999   dca3d_z_error = sqrt(rotCov(2,2));
2000 }
2001 
2002 //____________________________________________________________________________..
2003 //
2004 //  this method produces several different "dedx" values
2005 //      70% truncated mean 
2006 //      85% truncated mean 
2007 //      100% mean
2008 //      mean-0.2*sig
2009 //  also saved available is first 3 moments: mean, variance (rms^2), skewness
2010 //
2011 //double Collect::calc_dedx(TrackSeed* tpcseed, TrkrClusterContainer* clustermap, PHG4TpcGeomContainer* tpcGeom){
2012 //
2013 std::vector<double> Collect::calc_dedx(TrackSeed* tpcseed, TrkrClusterContainer* clustermap, PHG4TpcGeomContainer* tpcGeom){
2014 //
2015     //
2016     if (!clustermap){
2017         cout<<"no cluster map... cannot continue"<<endl;
2018         return {};
2019     }
2020     if (!tpcGeom){
2021         cout<<"no tpc geom container (TPCGEOMCONTAINER)... cannot continue"<<endl;
2022         return {};
2023     }
2024     //
2025     std::vector<TrkrDefs::cluskey> clusterKeys;
2026     clusterKeys.insert(clusterKeys.end(), tpcseed->begin_cluster_keys(), tpcseed->end_cluster_keys());
2027     std::vector<float> dedxlist;
2028     for (unsigned long cluster_key : clusterKeys){
2029         auto detid = TrkrDefs::getTrkrId(cluster_key);
2030         if (detid != TrkrDefs::TrkrId::tpcId){ continue; }   // the micromegas clusters are added to the TPC seeds
2031         TrkrCluster* cluster                = clustermap->findCluster(cluster_key);
2032         if (!cluster) continue;     //!!!!! PROTECTION DST MISMATCH !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2033         unsigned int layer_local            = TrkrDefs::getLayer(cluster_key);
2034         PHG4TpcGeom* GeoLayer_local = tpcGeom->GetLayerCellGeom(layer_local);
2035         if (!GeoLayer_local) continue;  //!!!!! PROTECTION: layer not in this geometry container !!!!!
2036         float adc       = cluster->getAdc();
2037         float thick     = GeoLayer_local->get_thickness();
2038         float r         = GeoLayer_local->get_radius();
2039         float alpha     = (r * r) / (2 * r * TMath::Abs(1.0 / tpcseed->get_qOverR()));
2040         float beta      = atan(tpcseed->get_slope());
2041         float alphacorr = cos(alpha);
2042         float betacorr  = cos(beta);        // betacorr is exactly the same as 1./cosh(tpcseed->get_eta())      
2043 //?     if( alphacorr<0 || alphacorr >4 ){ alphacorr=4; }
2044 //?     if( betacorr <0 || betacorr  >4 ){ betacorr =4; }
2045         adc /= thick;
2046         adc *= alphacorr;
2047         adc *= betacorr;
2048         dedxlist.push_back(adc);
2049 //if (DEBUG) cout<<"DEDX: got a cluster? adc= "<<adc<<endl;
2050     }
2051 //if (DEBUG) cout<<"DEDX: done clus loop. dedxsize= "<<dedxlist.size()<<endl;
2052     if (dedxlist.size() < 1){
2053         return {};
2054     }
2055 //if (DEBUG) cout<<"DEDX: sorting... "<<endl;
2056     sort(dedxlist.begin(), dedxlist.end());
2057 //if (DEBUG) cout<<"DEDX: sorted size= "<<dedxlist.size()<<endl;
2058     //
2059     int trunc_max   = (int)(dedxlist.size() * 0.70);
2060     double dedx70n  = 0.;
2061     double dedx70s  = 0.;
2062     for (int i = 0; i < trunc_max; ++i){
2063         dedx70s += dedxlist.at(i);
2064         dedx70n++;
2065     }
2066     if (dedx70n>0) dedx70s /= dedx70n;
2067     //
2068     return { dedx70n, dedx70s };
2069     //
2070 }
2071 
2072 //____________________________________________________________________________..
2073 // STAR-style Splitting Level (SL) split-track logic - Y. Adams et al. (STAR),
2074 // nucl-ex/0411036 sec. 2.4.1 - using real TPC layer (pad-row) occupancy from
2075 // fTlayermask (persisted; see Collect.h), masked here to just the TPC bits
2076 // (7-54): MVTX/INTT bits are deliberately excluded from the SL arithmetic
2077 // because a real split pair shares one common silicon seed by construction
2078 // (README_SplitTracks.md sec. 3), so those bits are identical for both
2079 // fragments and would only dilute/bias the TPC-specific SL value, not
2080 // discriminate anything.
2081 //
2082 // This function exists SOLELY to support fKillSplitTracks's own production-
2083 // time culling of SIBLING (same trigger-frame) pairs - nothing it computes is
2084 // written to the tree. Mixed-event pairs are structurally invisible to
2085 // Collect (two different trigger frames, processed at two different times),
2086 // so outread_stream/CalcRm.cxx must implement this identical SL(i,j)
2087 // calculation itself from fTlayermask to cover mixed pairs anyway (see README
2088 // sec. 8) - once that general function exists there, it also covers the
2089 // sibling case, making a separately-persisted sibling-only answer here purely
2090 // redundant (and a second, possibly-diverging "which pairs are split"
2091 // bookkeeping to keep in sync). If fKillSplitTracks stays false (the
2092 // default), this function is a no-op past the pairing loop.
2093 //
2094 // Pairing is restricted to same reconstructed charge: splitting is one real
2095 // particle refit as two tracks, so it can never flip charge (matches STAR's
2096 // own LS-only scope for this cut). The (deta,dphi) gate below is purely a
2097 // computational shortcut to skip obviously-unrelated pairs before the
2098 // popcount - SL itself is what actually discriminates split tracks, and
2099 // Collect never thresholds on SL's value to decide anything.
2100 void Collect::FlagSplitTracks(){
2101     //
2102     if (!fKillSplitTracks) return;
2103     //
2104     constexpr uint64_t TPC_LAYER_MASK   = ((1ULL<<48)-1) << 7;  // bits 7-54 set
2105     const double DETA_GATE  = 0.01;     // = outread_loop.cxx's already-validated window (README_SplitTracks.md sec. 3/5) -
2106     const double DPHI_GATE  = 0.03;     // the CF spike this whole investigation started from is ~single-bin sharp, so there is
2107                                         // no evidence real splits extend beyond this; deliberately NOT opened up further.
2108     //
2109     std::vector<double> splitSL(fNTracks, 9.);      // local only - sentinel: no close same-charge pair found
2110     std::vector<int>    splitPartner(fNTracks, -1); // local only
2111     std::vector<int>    splitLoser(fNTracks, 0);    // local only
2112     //
2113     for (int i=0;i<fNTracks;i++){
2114         uint64_t maski = fTlayermask[i] & TPC_LAYER_MASK;
2115         if (maski==0) continue;
2116         for (int j=i+1;j<fNTracks;j++){
2117             uint64_t maskj = fTlayermask[j] & TPC_LAYER_MASK;
2118             if (maskj==0) continue;
2119             if (fTchg[i] != fTchg[j]) continue;         // splitting cannot flip charge
2120             double deta = fabs(fTeta[i]-fTeta[j]);
2121             double dphi = fabs(fTphi[i]-fTphi[j]);
2122             if (dphi>M_PI) dphi = 2.*M_PI - dphi;
2123             if (deta>DETA_GATE || dphi>DPHI_GATE) continue;
2124             //
2125             int nboth   = __builtin_popcountll( maski & maskj );
2126             int nxor    = __builtin_popcountll( maski ^ maskj );
2127             int n1      = __builtin_popcountll( maski );
2128             int n2      = __builtin_popcountll( maskj );
2129             if ((n1+n2)==0) continue;
2130             double SL   = ( (double)nxor - (double)nboth ) / (double)(n1+n2);       // STAR eq.: SL=-0.5 -> same track split, SL=1.0 -> independent tracks
2131             //
2132             if (SL < splitSL[i]){ splitSL[i] = SL; splitPartner[i] = j; }       // keep the most split-like (lowest SL) partner per track
2133             if (SL < splitSL[j]){ splitSL[j] = SL; splitPartner[j] = i; }
2134         }   // end j loop...
2135     }   // end i loop...
2136     //
2137     //---- ALWAYS decide winner/loser for every found pair - independent of SL's value, gated only by having a partner at all
2138     for (int it=0;it<fNTracks;it++){
2139         int jt = splitPartner[it];
2140         if (jt<0) continue;
2141         //
2142         if      (fTnhitstpc[it] != fTnhitstpc[jt]){ splitLoser[it] = (fTnhitstpc[it] < fTnhitstpc[jt]) ? 1 : 0; }   // fewer TPC hits loses - directly which fragment got the smaller share of the real hits (README sec. 3)
2143         else if (fTquality[it]  != fTquality[jt] ){ splitLoser[it] = (fTquality[it]  > fTquality[jt])  ? 1 : 0; }   // ...tie-break: worse (higher) quality loses - direct fit-consistency metric (README sec. 5)
2144         // else: both tied on ntpc AND quality - no independently-justified way to pick a loser (pt is only a
2145         // correlated echo of ntpc, not new evidence), so leave both at their default splitLoser=0 - don't touch a
2146         // track without real evidence something is wrong with it.
2147     }   // end winner/loser loop...
2148     //
2149     for (int it=0;it<fNTracks;it++){
2150         if (!splitLoser[it]) continue;
2151         if (fTindv0[it] >= 0) continue;                 // never kill a V0 daughter
2152         fTprimary[it] = -1;                             // sentinel already understood downstream as "skip this track"
2153     }   // end kill loop...
2154     //
2155 }
2156 
2157 //____________________________________________________________________________..
2158 //---- reset any counters/arrays used in/for single events
2159 int Collect::ResetEvent(PHCompositeNode *topNode){
2160     //std::cout << "Collect::ResetEvent -- Resetting internal structures, prepare for next event" << std::endl;
2161     fNTracks    = 0;
2162     return Fun4AllReturnCodes::EVENT_OK;
2163 }
2164 
2165 //____________________________________________________________________________..
2166 //---- called at end of sim run or end of entire job
2167 int Collect::EndRun(const int runnumber){
2168     std::cout << "Collect::EndRun -- Ending Run for Run " << runnumber << std::endl;
2169     return Fun4AllReturnCodes::EVENT_OK;
2170 }
2171 
2172 //____________________________________________________________________________..
2173 //
2174 //  take a specific tprofile and find the calorimeter tower edges in eta
2175 //
2176 void Collect::DoEdgesEta( TProfile* hmap, std::vector<double>& edges, TString layerstr ){
2177     //
2178     edges.clear();
2179     if (!hmap){ cout<<"Collect::DoEdgesEta -- no map data. "<<endl; return; }
2180     //
2181     //---- get some info from map hists...
2182     int nbins       = hmap->GetNbinsX();
2183     double detamin  = fabs(hmap->GetBinContent(2)-hmap->GetBinContent(1));
2184     double detamax  = fabs(hmap->GetBinContent(nbins)-hmap->GetBinContent(nbins-1));
2185     for (int ib=1;ib<=hmap->GetNbinsX();ib++){
2186         double etac = hmap->GetBinContent(ib);
2187         double deta=0.,val=0.;
2188         if (ib==1){ 
2189             val     = etac-detamin/2.;
2190         } else if (ib> 1){ 
2191             deta    = hmap->GetBinContent(ib) - hmap->GetBinContent(ib-1);
2192             val     = etac-deta/2.;
2193         }
2194         if (fabs(val)<1.0E-10) val=0.0;
2195         edges.push_back(val);
2196         if (ib==nbins){             // now add one more (the upper limit)
2197             edges.push_back(etac+detamax/2.);
2198         }
2199     }
2200     //  
2201     cout<<"\tdouble Edge"<<layerstr.Data()<<"eta["<<edges.size()<<"] = {";
2202     int icnt = 0;
2203     for (size_t i=0;i<edges.size();i++){
2204         cout<<edges.at(i);
2205         if (i<edges.size()-1) cout<<",";
2206         icnt++;
2207         if (icnt==15){ cout<<endl<<"\t\t\t\t"; icnt=0; }
2208     }   cout<<"};"<<endl;
2209     //
2210 }
2211 
2212 //____________________________________________________________________________..
2213 //
2214 //  take a specific tprofile and find the calorimeter tower edges in phi
2215 //
2216 void Collect::DoEdgesPhi( TProfile* hmap, std::vector<double>& edges ){
2217     //
2218     edges.clear();
2219     if (!hmap){ cout<<"Collect::DoEdgesEta -- no map data. "<<endl; return; }
2220     //
2221     //---- get some info from map hists...
2222     int nbins       = hmap->GetNbinsX();
2223     double dphi     = fabs(hmap->GetBinContent(12)-hmap->GetBinContent(11));
2224     for (int ib=1;ib<=hmap->GetNbinsX();ib++){
2225         double phic = hmap->GetBinContent(ib);
2226         double val=0.;
2227         val = phic-dphi/2.;
2228         if (fabs(val)<1.0E-10) val=0.0;
2229         edges.push_back(val);
2230         if (ib==nbins){             // now add one more (the upper limit)
2231             edges.push_back(phic+dphi/2.);
2232         }
2233     }
2234     //
2235 }
2236 
2237 //____________________________________________________________________________..
2238 //---- called at end of all processing in this job
2239 int Collect::End(PHCompositeNode *topNode){
2240 
2241     std::cout << "Collect::End -- This is the End... " << std::endl;
2242 
2243     if (CALODETAIL){
2244         std::vector<double> EdgesEM_eta;    
2245         DoEdgesEta( hmapEM_ieta, EdgesEM_eta, TString("EM") );
2246         std::vector<double> EdgesIH_eta;    
2247         DoEdgesEta( hmapIH_ieta, EdgesIH_eta, TString("IH") );
2248         std::vector<double> EdgesOH_eta;    
2249         DoEdgesEta( hmapOH_ieta, EdgesOH_eta, TString("OH") );
2250         std::vector<double> EdgesEM_phi;    
2251         DoEdgesPhi( hmapEM_iphi, EdgesEM_phi );
2252         std::vector<double> EdgesIH_phi;    
2253         DoEdgesPhi( hmapIH_iphi, EdgesIH_phi );
2254         std::vector<double> EdgesOH_phi;    
2255         DoEdgesPhi( hmapOH_iphi, EdgesOH_phi );
2256     }
2257     
2258     //---- write and close...
2259     std::cout << "Collect::End -- This is the End...  Writing root objects... " << std::endl;
2260     fout->cd(); 
2261     TString origtitle   = hnvtx_frame->GetTitle();
2262     TString newtitle    = origtitle + TString(Form(", Offset= %d, RunNum= %d",fEvtSeqFirst,fRunNum));
2263     hnvtx_frame->SetTitle(newtitle);
2264     hnvtx_frame->Write();
2265     hdXing->Write();
2266     hech_em->Write();
2267     hech_ih->Write();
2268     hech_oh->Write();
2269     hech_ep->Write();
2270     //
2271     if (CALODETAIL) {
2272         hmapEM_eta->Write();
2273         hmapEM_phi->Write();
2274         hmapIH_eta->Write();
2275         hmapIH_phi->Write();
2276         hmapOH_eta->Write();
2277         hmapOH_phi->Write();
2278         hmapEM_eta2->Write();
2279         hmapEM_phi2->Write();
2280         hmapIH_eta2->Write();
2281         hmapIH_phi2->Write();
2282         hmapOH_eta2->Write();
2283         hmapOH_phi2->Write();
2284         hmapEM_ieta->Write();
2285         hmapEM_iphi->Write();
2286         hmapIH_ieta->Write();
2287         hmapIH_iphi->Write();
2288         hmapOH_ieta->Write();
2289         hmapOH_iphi->Write();
2290         //for (int im=0;im<MatchCounter;im++){
2291         //  hmatch_EMetaphiE[im]->Write();
2292         //  hmatch_HCetaphiE[im]->Write();
2293         //}
2294         for (int iev=0;iev<MAXEVTDETAIL;iev++){
2295             for (int ilayer=0;ilayer<3;ilayer++){
2296                 hlayer[ilayer][iev]         ->Write();
2297                 hlayer_efDET[ilayer][iev]   ->Write();
2298             }
2299             hPTOT_etaDETem[iev] ->Write();
2300             hPTOT_phiDETem[iev] ->Write();
2301             hPTOT_etaDETih[iev] ->Write();
2302             hPTOT_phiDETih[iev] ->Write();
2303             hPTOT_etaDEToh[iev] ->Write();
2304             hPTOT_phiDEToh[iev] ->Write();
2305         }
2306         hPedestalEM ->Write();
2307         hPedestalIH ->Write();
2308         hPedestalOH ->Write();
2309         for (int iet=0;iet<96;iet++){
2310             for (int ifi=0;ifi<256;ifi++){
2311                     hEevt_em[iet][ifi]->Write();
2312                 if (iet<24&&ifi<64){
2313                     hEevt_ih[iet][ifi]->Write();
2314                     hEevt_oh[iet][ifi]->Write();
2315                 }
2316             }
2317         }
2318     }   // end CALODETAIL
2319     //
2320     // DIAGTEMP-20260920: throwaway merged-track exploration output - see README_SplitTracks.md sec.12b discussion
2321     fout->cd();
2322     std::cout<<"DIAGTEMP-20260920: nTracks="<<sDiagNTracks20260920<<" nClusters="<<sDiagNClusters20260920
2323             <<" nSharedKeyInstances="<<sDiagNSharedKeyInstances20260920<<" nEventsWithSharedKey="<<sDiagNEventsWithSharedKey20260920
2324             <<" nClosePairsChecked="<<sDiagNClosePairsChecked20260920<<" nClosePairsSharedKeyToo="<<sDiagNClosePairsSharedKeyToo20260920
2325             <<std::endl;
2326     std::cout<<"DIAGTEMP-20260920 gap-check: close(dist<2cm) asym="<<sDiagNGapAsymClose20260920<<"/"<<sDiagNGapChecksClose20260920
2327             <<" far(dist>10cm) asym="<<sDiagNGapAsymFar20260920<<"/"<<sDiagNGapChecksFar20260920
2328             <<std::endl;
2329     if (sDiagAdcAll20260920)       sDiagAdcAll20260920->Write();
2330     if (sDiagAdcSharedKey20260920) sDiagAdcSharedKey20260920->Write();
2331     if (sDiagAdcClosePair20260920) sDiagAdcClosePair20260920->Write();
2332     if (sDiagDist20260920)         sDiagDist20260920->Write();
2333     if (sDiagAdcVsDist20260920)    sDiagAdcVsDist20260920->Write();
2334     if (sDiagGapRateVsDist20260920) sDiagGapRateVsDist20260920->Write();
2335     if (sDiagDistSameSector20260920) sDiagDistSameSector20260920->Write();
2336     if (sDiagDistDiffSector20260920) sDiagDistDiffSector20260920->Write();
2337     if (sDiagGapRateVsDistSameSector20260920) sDiagGapRateVsDistSameSector20260920->Write();
2338     if (sDiagGapRateVsDistDiffSector20260920) sDiagGapRateVsDistDiffSector20260920->Write();
2339     //
2340     outTree->Write();
2341     fout->Close();
2342     std::cout << "Collect::End -- This is the End...  rootfile closed!" << std::endl;
2343     return Fun4AllReturnCodes::EVENT_OK;
2344     //
2345 }
2346 
2347 //____________________________________________________________________________..
2348 //---- not really used...
2349 int Collect::Reset(PHCompositeNode *topNode){
2350     std::cout << "Collect::Reset -- resetting..." << std::endl;
2351     return Fun4AllReturnCodes::EVENT_OK;
2352 }
2353 
2354 //____________________________________________________________________________..
2355 void Collect::Print(const std::string &what) const{
2356     std::cout << "Collect::Print -- Printing info for " << what << std::endl;
2357 }
2358
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!