Back to home page

sPhenix code displayed by LXR

 
 

    


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

0001 //  -*- C++ -*-.
0002 #ifndef COLLECT_H
0003 #define COLLECT_H
0004 
0005 #include <trackbase/ActsGeometry.h>
0006 #include <trackbase/ClusterErrorPara.h>
0007 #include <trackbase/TrkrDefs.h>
0008 #include <tpc/TpcGlobalPositionWrapper.h>
0009 #include <fun4all/SubsysReco.h>
0010 #include <TFile.h>
0011 #include <TTree.h>
0012 #include <TH1.h>
0013 #include <TH2.h>
0014 #include <TGraph.h>
0015 #include <TProfile.h>
0016 #include <TProfile2D.h>
0017 #include <cmath>
0018 #include <cstdint>
0019 #include <iostream>
0020 #include <limits>
0021 #include <string>
0022 
0023 class RunHeader;
0024 class EventHeader;
0025 class TrkrCluster;
0026 class PHCompositeNode;
0027 class ActsGeometry;
0028 //class SvtxTrackMap_v2;        //
0029 //class SvtxTrack_v4;           //
0030 class SvtxTrackMap;         //
0031 class SvtxTrack;            //
0032 class SvtxVertexMap;        //
0033 class SvtxVertex;           //
0034 class GlobalVertexMap;      //
0035 class GlobalVertex;         //
0036 class TrackSeed;
0037 class TrkrClusterContainer;
0038 class TrkrHitSetContainer;
0039 class PHG4TpcGeomContainer;
0040 class PHG4CylinderGeomContainer;
0041 class TpcDistortionCorrection;
0042 class TpcDistortionCorrectionContainer;
0043 class TTree;
0044 class TFile;
0045 class KFParticle_Container;
0046 
0047 class Collect : public SubsysReco
0048 {
0049  public:
0050 
0051     Collect(const std::string &name = "Collect");
0052     ~Collect() override;
0053     int Init(PHCompositeNode *topNode) override;
0054     int InitRun(PHCompositeNode *topNode) override;
0055     int process_event(PHCompositeNode *topNode) override;
0056     int ResetEvent(PHCompositeNode *topNode) override;
0057     int EndRun(const int runnumber) override;
0058     int End(PHCompositeNode *topNode) override;
0059     int Reset(PHCompositeNode * /*topNode*/) override;
0060     void Print(const std::string &what = "ALL") const override;
0061     //
0062     void outfileName(const std::string &name) { 
0063                         foutfileName = name; 
0064                         std::cout<<"Collect::outfileName -- "<<foutfileName<<std::endl; 
0065                     }
0066     void setRequireXingZero(bool value){
0067                         fRequireXingZero    = value;
0068                         if ( fRequireXingZero){ std::cout<<"Collect::fRequireXingZero is TRUE "<<std::endl; } else
0069                         if (!fRequireXingZero){ std::cout<<"Collect::fRequireXingZero is FALSE"<<std::endl; }
0070                     }
0071     void setAddTowersToTrackTree(bool value){
0072                         fAddTowersToTrackTree       = value;
0073                         if ( fAddTowersToTrackTree){ std::cout<<"Collect::AddTowersToTrackTree is TRUE "<<std::endl; } else
0074                         if (!fAddTowersToTrackTree){ std::cout<<"Collect::AddTowersToTrackTree is FALSE"<<std::endl; }
0075                     }
0076     void setVerbosity(int iverb){ 
0077                         kVerbosity = iverb; 
0078                              if (kVerbosity==0){ DEBUG=false; DEBUG2=false; }
0079                         else if (kVerbosity==1){ DEBUG=true;  DEBUG2=false; }
0080                         else if (kVerbosity> 1){ DEBUG=true;  DEBUG2=true;  }
0081                     }
0082     void setSegment(int seg){
0083                         fSegment = seg;
0084                     }
0085     // Take "segment" per event from the filename the named input manager is
0086     // currently reading (overrides setSegment) - for jobs reading several
0087     // segments. Empty (default) = keep the fixed setSegment() value.
0088     void setSegmentFromInputManager(const std::string &name){
0089                         fSegmentInputManager = name;
0090                     }
0091     void setKillSplitTracks(bool value){
0092                         fKillSplitTracks    = value;
0093                         if ( fKillSplitTracks){ std::cout<<"Collect::fKillSplitTracks is TRUE "<<std::endl; } else
0094                         if (!fKillSplitTracks){ std::cout<<"Collect::fKillSplitTracks is FALSE"<<std::endl; }
0095                     }
0096     //
0097     void calc3DDCA(SvtxTrack* track, SvtxVertexMap* vertexmap, unsigned int vtxid,
0098          double& dca3d_xy, double& dca3d_xy_error,
0099          double& dca3d_z , double& dca3d_z_error );
0100     // 
0101     std::vector<double> calc_dedx(TrackSeed *tpcseed,
0102                 TrkrClusterContainer *clusters, 
0103                 PHG4TpcGeomContainer *tpcGeom);
0104     //
0105     int GetKFP( PHCompositeNode *topNode,
0106                  const int parentpdgid,
0107                  const std::string &containername,
0108                  const std::string &ctrackmapname );
0109     //
0110     // STAR-style ("SL", nucl-ex/0411036 sec. 2.4.1) split-track logic, using the
0111     // persisted per-track layer-occupancy mask (fTlayermask) - see
0112     // README_SplitTracks.md sec. 5-8. SL/pair-finding/winner-loser are computed
0113     // here purely to support fKillSplitTracks's own production-time culling of
0114     // SIBLING (same trigger-frame) pairs - none of it is written to the tree as
0115     // separate branches, since mixed-event pairs (which this function structurally
0116     // cannot see) also need the identical calculation redone downstream anyway;
0117     // outread_stream must implement its own general SL(i,j) from fTlayermask for
0118     // both sibling and mixed pairs regardless of what Collect does internally, so
0119     // there is no benefit to also persisting Collect's sibling-only, gate-frozen
0120     // answer (see README sec. 8). Called once per tree row, after the track loop
0121     // fills fT*[0..fNTracks-1], before Fill(). Never touches a V0 daughter
0122     // (fTindv0>=0), and does nothing at all unless fKillSplitTracks is true.
0123     void FlagSplitTracks();
0124     //
0125     void DoEdgesEta( TProfile* h, std::vector<double>& edges, TString str );
0126     void DoEdgesPhi( TProfile* h, std::vector<double>& edges );
0127     int GetIETA( int layer, double e );
0128     int GetIPHI( int layer, double f );
0129     //
0130  private:
0131     //
0132     static const bool CALODETAIL = false;   // enables extra plots, BIG
0133     //
0134     static const int MAXNTRACKS     = 500;
0135     static const int MAXNV0S        =  50;
0136     static const int MAXEVTDETAIL   =  20;
0137     double THRESH_em;
0138     double THRESH_ih;
0139     double THRESH_oh;
0140     double THRESH_ep;
0141     int Kshort_PDGID;
0142     int Lambda_PDGID;
0143     int ALambda_PDGID;
0144     //
0145     TH1D*   hnvtx_frame;
0146     TH1D*   hdXing;
0147     TH1D*   hech_em;
0148     TH1D*   hech_ih;
0149     TH1D*   hech_oh;
0150     TH1D*   hech_ep;
0151     TProfile*   hmapEM_eta;
0152     TProfile*   hmapEM_phi;
0153     TProfile*   hmapIH_eta;
0154     TProfile*   hmapIH_phi;
0155     TProfile*   hmapOH_eta;
0156     TProfile*   hmapOH_phi;
0157     TProfile*   hmapEM_eta2;
0158     TProfile*   hmapEM_phi2;
0159     TProfile*   hmapIH_eta2;
0160     TProfile*   hmapIH_phi2;
0161     TProfile*   hmapOH_eta2;
0162     TProfile*   hmapOH_phi2;
0163     TProfile*   hmapEM_ieta;
0164     TProfile*   hmapEM_iphi;
0165     TProfile*   hmapIH_ieta;
0166     TProfile*   hmapIH_iphi;
0167     TProfile*   hmapOH_ieta;
0168     TProfile*   hmapOH_iphi;
0169 //  int         MatchCounter;
0170 //  TH2D*       hmatch_EMetaphiE[100];
0171 //  TH2D*       hmatch_HCetaphiE[100];
0172 //  TGraph*     gmatch_EMetaphi[100];
0173 //  TGraph*     gmatch_HCetaphi[100];
0174     TH2D*   hlayer[3][MAXEVTDETAIL];
0175     TH2D*   hlayer_efDET[3][MAXEVTDETAIL];
0176     TH1D*   hPTOT_etaDETem[MAXEVTDETAIL];
0177     TH1D*   hPTOT_phiDETem[MAXEVTDETAIL];
0178     TH1D*   hPTOT_etaDETih[MAXEVTDETAIL];
0179     TH1D*   hPTOT_phiDETih[MAXEVTDETAIL];
0180     TH1D*   hPTOT_etaDEToh[MAXEVTDETAIL];
0181     TH1D*   hPTOT_phiDEToh[MAXEVTDETAIL];
0182     TProfile2D* hPedestalEM;
0183     TProfile2D* hPedestalIH;
0184     TProfile2D* hPedestalOH;
0185 //  TProfile2D* hPedestalCorrEM;
0186 //  TProfile2D* hPedestalCorrIH;
0187 //  TProfile2D* hPedestalCorrOH;
0188     TH1D* hEevt_em[96][256];
0189     TH1D* hEevt_ih[24][64];
0190     TH1D* hEevt_oh[24][64];
0191     //
0192     TpcGlobalPositionWrapper m_globalPositionWrapper;   // crossing+distortion-corrected TPC global position - README_SplitTracks.md sec.12b
0193     bool    fRequireXingZero;
0194     bool    fAddTowersToTrackTree;
0195     bool    fHasCaloGeo;
0196     bool    fKillSplitTracks;       // see FlagSplitTracks() / README_SplitTracks.md; default off until validated on real data
0197     double  PTLL;       // track pT lower limit (usually 0.1)
0198     //
0199     int     kVerbosity;
0200     bool    DEBUG;      // basic verbosity
0201     bool    DEBUG2;     // enhanced verbosity
0202     TString foutfileName;
0203     TFile*  fout;
0204     TTree*  outTree;
0205     int     fRunNum;
0206     int     fSegment;
0207     std::string fSegmentInputManager;   // see setSegmentFromInputManager()
0208     int     fEvtSeq;
0209     int     fEvtSeqFirst;
0210     int     ievtSeen;
0211     int     ntpcmin;
0212     int     nmvtmin;
0213     int     nintmin;
0214     //
0215     double caloRadiusInnerEMCal;
0216     double caloRadiusInnerOHCal;
0217     double caloRadiusInnerIHCal;
0218     double caloThicknessEMCal;
0219     double caloThicknessOHCal;
0220     double caloThicknessIHCal;
0221     double caloRadiusOuterEMCal;
0222     double caloRadiusOuterOHCal;
0223     double caloRadiusOuterIHCal;
0224     double caloRadiusMidEMCal;
0225     double caloRadiusMidOHCal;
0226     double caloRadiusMidIHCal;
0227     //
0228     int sumNfoundKS;
0229     int sumNfoundLA;
0230     int sumNfoundAL;
0231     //
0232     int             V0Crossing[MAXNV0S] ;
0233     int             V0pid[MAXNV0S]      ;
0234     unsigned int    V0Track1id[MAXNV0S];
0235     unsigned int    V0Track2id[MAXNV0S];
0236     double          V0DecayX[MAXNV0S];  // daughter-pair vertex position (NOT the mother's own GetX/Y/Z! -
0237     double          V0DecayY[MAXNV0S];  // see 2026-09-17 dev-log entry: the mother KFParticle is transported
0238     double          V0DecayZ[MAXNV0S];  // to the production vertex before Collect ever reads it, so its own
0239                                         // position is ~the PV, not the decay vertex - daughters are not
0240                                         // transported and still sit at the decay vertex where they were fit.
0241     //
0242     std::vector<bool> fTtriggerVector;
0243     ULong64_t   fTbco;                  // GL1 BCO of this trigger frame (absolute beam clock) - lets outread_stream tell whether two
0244                                         // rows' trigger frames could overlap in time; ~0ULL if no GL1RAWHIT (README_SplitTracks.md sec.12a.6)
0245     double  fTTotE_em;
0246     double  fTTotE_ih;
0247     double  fTTotE_oh;
0248     double  fTTotE_ioh;
0249     double  fTTotE_epd;
0250     int     fTnvtx;
0251     int     fTvtxntr;
0252     int     fTnvtxXing;     // number of vertices sharing this row's crossing (all sizes; >1 = pileup or split vertex)
0253     int     fTvtxntr2;      // attached-track count (same count as vtxntr) of the runner-up vertex in this crossing, 0 if none
0254     double  fTvtxx;
0255     double  fTvtxy;
0256     double  fTvtxz;
0257     double  fTvtxzGlob;
0258     int     fTcrossing;
0259     //
0260     int     fNTracks;
0261     Char_t  fTpid[MAXNTRACKS];      // tag for v0 daughters!  0=unknown, 1=pi, 2=K, 3=p - README_SplitTracks.md sec.12c: was Int_t, now Char_t (signed, small range)
0262     Short_t fTindv0[MAXNTRACKS];    // index of v0 in tree v0 arrays, -1 if this track not a daughter of a found v0 - sec.12c: was Int_t, now Short_t
0263     Char_t  fTprimary[MAXNTRACKS];  // =1 is matched to hottest vertex, =2 if matched to /a/ vertex, =0 otherwise, =-1 killed by FlagSplitTracks() - sec.12c: was Int_t, now Char_t (signed - -1 is real)
0264     Float_t fTquality[MAXNTRACKS];  // sec.12c: was Double_t, now Float_t (throughout this block)
0265     Float_t fTchisq[MAXNTRACKS];
0266     Float_t fTchg[MAXNTRACKS];
0267     double  fTpx[MAXNTRACKS];
0268     double  fTpy[MAXNTRACKS];
0269     double  fTpz[MAXNTRACKS];
0270     Float_t fTptot[MAXNTRACKS];
0271     Float_t fTeta[MAXNTRACKS];
0272     Float_t fTphi[MAXNTRACKS];
0273     Float_t fTpt[MAXNTRACKS];
0274     Float_t fTseedeta[MAXNTRACKS];
0275     Float_t fTseedphi[MAXNTRACKS];
0276     Float_t fTseedpt[MAXNTRACKS];
0277     UChar_t fTnhitstpc[MAXNTRACKS]; // sec.12c: was Int_t, now UChar_t (always >=0, <=48)
0278     UChar_t fTnhitsmvt[MAXNTRACKS];
0279     UChar_t fTnhitsint[MAXNTRACKS];
0280     Float_t fTdcaxy[MAXNTRACKS];
0281     Float_t fTdcaz[MAXNTRACKS];
0282     //
0283     UChar_t fTdedx70n[MAXNTRACKS];  // sec.12c: was Int_t, now UChar_t
0284     Float_t fTdedx70s[MAXNTRACKS];  // sec.12c: was Double_t, now Float_t
0285     Float_t fTdedxKFP[MAXNTRACKS];  // TrackAnalysisUtils::calc_dedx() - the exact function KFParticle_Tools::get_dEdx()
0286                                     // calls internally for the usePID() dE/dx cut - written alongside our own
0287                                     // dedx70s as a direct apples-to-apples cross-check (2026-09-18)
0288     //
0289     uint64_t    fTlayermask[MAXNTRACKS];    // one bit per absolute TrkrDefs layer (0-2 MVTX, 3-6 INTT, 7-54 TPC) with a cluster on this track -
0290                                         // persisted (branch "layermask") so outread_stream can compute STAR-style SL itself for ANY pair
0291                                         // (sibling or mixed - see README_SplitTracks.md sec. 8); Collect's own internal use (FlagSplitTracks(),
0292                                         // for fKillSplitTracks) masks this to the TPC bits only before running the SL arithmetic.
0293     TrkrDefs::cluskey fTsiclukey[MAXNTRACKS][7];    // raw silicon cluster key per absolute Si layer (slot index == layermask's bit 0-6,
0294                                         // 0-2 MVTX / 3-6 INTT); sentinel ~0ULL for "no cluster on this layer" (0 itself is NOT safe -
0295                                         // see README_SplitTracks.md sec. 12.3). Additive next to layermask: lets outread_stream check
0296                                         // literal shared-cluster identity between two tracks, not just coverage-pattern agreement
0297                                         // (README_SplitTracks.md sec. 12).
0298     // Merged-track ("STAR FMH"-adjacent) support - README_SplitTracks.md sec.12b.
0299     // Compact per-TPC-layer local position, slot index = absolute layer - 7 (0-47,
0300     // TPC layers 7-54). All three arrays share one sentinel: tpcsector==255 means
0301     // "no TPC cluster on this layer" - tpcarclen/tpcz are don't-care in that case,
0302     // they do NOT carry an independent sentinel. Real merging is only physically
0303     // possible within one TPC sector (sectors are separated by dead-space frames -
0304     // section 12b design discussion, 2026-09-20), so two tracks' local positions
0305     // are only meaningfully comparable when their tpcsector values match.
0306     UChar_t fTtpcsector[MAXNTRACKS][48];    // TpcDefs sector*2+side, 0-23; 255 = no cluster on this layer (sentinel)
0307     Char_t  fTtpcarclen[MAXNTRACKS][48];    // local arc-length within that sector, signed, 2.0 mm/count (fixed scale,
0308                                         // sized for the outermost-layer sector half-width with margin - same scale
0309                                         // for every layer, no per-layer geometry lookup needed downstream)
0310     Short_t fTtpcz[MAXNTRACKS][48];     // global z, signed, 0.1 mm/count (huge margin vs the ~110cm TPC half-length)
0311     //
0312     double    fTstateEM_eta[MAXNTRACKS];    
0313     double    fTstateEM_phi[MAXNTRACKS];    
0314     double      fTstateEM_x[MAXNTRACKS];    
0315     double      fTstateEM_y[MAXNTRACKS];    
0316     double      fTstateEM_z[MAXNTRACKS];    
0317     double  fTegridEM[MAXNTRACKS][7][7];
0318     double  fTegridEMsum7x7[MAXNTRACKS];
0319     double  fTeopEM[MAXNTRACKS];    // EMCal 7x7-sum / track momentum ("E/p")
0320     double    fTstateIH_eta[MAXNTRACKS];
0321     double    fTstateIH_phi[MAXNTRACKS];    
0322     double      fTstateIH_x[MAXNTRACKS];    
0323     double      fTstateIH_y[MAXNTRACKS];    
0324     double      fTstateIH_z[MAXNTRACKS];    
0325     double  fTegridIH[MAXNTRACKS][7][7];
0326     double  fTegridIHsum7x7[MAXNTRACKS];
0327     double  fTeopIH[MAXNTRACKS];    // IHCal 7x7-sum / track momentum ("E/p") - internal ingredient for eopHC/eopCA, not itself a branch
0328     double    fTstateOH_eta[MAXNTRACKS];
0329     double    fTstateOH_phi[MAXNTRACKS];    
0330     double      fTstateOH_x[MAXNTRACKS];    
0331     double      fTstateOH_y[MAXNTRACKS];    
0332     double      fTstateOH_z[MAXNTRACKS];    
0333     double  fTegridOH[MAXNTRACKS][7][7];
0334     double  fTegridOHsum7x7[MAXNTRACKS];
0335     double  fTeopOH[MAXNTRACKS];    // OHCal 7x7-sum / track momentum ("E/p") - internal ingredient for eopHC/eopCA, not itself a branch
0336     double  fTeopHC[MAXNTRACKS];    // (IHCal+OHCal) 7x7-sum / track momentum - "E/p" for hadronic calorimetry combined
0337     double  fTeopCA[MAXNTRACKS];    // (EMCal+IHCal+OHCal) 7x7-sum / track momentum - "E/p" for all 3 calorimeters combined
0338     //
0339     double      fTstateTI_x[MAXNTRACKS];    
0340     double      fTstateTI_y[MAXNTRACKS];    
0341     double      fTstateTI_z[MAXNTRACKS];    
0342     double      fTstateTO_x[MAXNTRACKS];    
0343     double      fTstateTO_y[MAXNTRACKS];    
0344     double      fTstateTO_z[MAXNTRACKS];    
0345     //
0346     int     nv0tf;
0347     bool    fTv0TFuseit[MAXNV0S];       // internal only
0348     int     fTv0TFxing[MAXNV0S];        // internal only
0349     int     fTv0TFpid[MAXNV0S];
0350     double  fTv0TFx[MAXNV0S];
0351     double  fTv0TFy[MAXNV0S];
0352     double  fTv0TFz[MAXNV0S];
0353     double  fTv0TFpx[MAXNV0S];
0354     double  fTv0TFpy[MAXNV0S];
0355     double  fTv0TFpz[MAXNV0S];
0356     double  fTv0TFptot[MAXNV0S];
0357     double  fTv0TFeta[MAXNV0S];
0358     double  fTv0TFphi[MAXNV0S];
0359     double  fTv0TFpt[MAXNV0S];
0360     double  fTv0TFene[MAXNV0S];
0361     double  fTv0TFmass[MAXNV0S];
0362     double  fTv0TFctau[MAXNV0S];
0363     double  fTv0TFdecaylen[MAXNV0S];
0364     double  fTv0TFchi2ndf[MAXNV0S];
0365     int     fTv0TFindtr1[MAXNV0S];  // index of this v0's track1 in tracks tree!
0366     int     fTv0TFindtr2[MAXNV0S];  // index of this v0's track2 in tracks tree!
0367     //
0368     int     fTnv0;      // MAXNV0S
0369     int     fTv0pid[MAXNV0S];
0370     double  fTv0x[MAXNV0S];
0371     double  fTv0y[MAXNV0S];
0372     double  fTv0z[MAXNV0S];
0373     double  fTv0px[MAXNV0S];
0374     double  fTv0py[MAXNV0S];
0375     double  fTv0pz[MAXNV0S];
0376     double  fTv0ptot[MAXNV0S];
0377     double  fTv0eta[MAXNV0S];
0378     double  fTv0phi[MAXNV0S];
0379     double  fTv0pt[MAXNV0S];
0380     double  fTv0ene[MAXNV0S];
0381     double  fTv0mass[MAXNV0S];
0382     double  fTv0ctau[MAXNV0S];
0383     double  fTv0decaylen[MAXNV0S];
0384     double  fTv0chi2ndf[MAXNV0S];
0385     double  fTv0dira[MAXNV0S];      // cos(pointing angle): V0 momentum vs (decay vtx - PV)
0386     double  fTv0pvdca[MAXNV0S];     // V0 trajectory DCA to the primary vertex (cm)
0387     double  fTv0decayx[MAXNV0S];    // daughter-pair decay vertex (see V0DecayX comment) - what dira/pvdca are computed from
0388     double  fTv0decayy[MAXNV0S];
0389     double  fTv0decayz[MAXNV0S];
0390     int     fTv0indtr1[MAXNV0S];    // index of this v0's track1 in tracks tree!
0391     int     fTv0indtr2[MAXNV0S];    // index of this v0's track2 in tracks tree!
0392     int     fV0FinalToTF[MAXNV0S];  // internal only: final (fTv0*) index -> trigger-frame (V0Track1id/etc) index, so the daughter-matching loop reads the right trigger-frame slot
0393     // needs v0 chi2, decayLength, ctau
0394     //
0395 //  unsigned int thisndedx;
0396     //
0397 };
0398 
0399 #endif // COLLECT_H