Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-03 08:17:42

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003 #ifndef TRACKRESIDUALS_H
0004 #define TRACKRESIDUALS_H
0005 
0006 #include <tpc/TpcClusterMover.h>
0007 #include <tpc/TpcGlobalPositionWrapper.h>
0008 
0009 #include <trackbase/ClusterErrorPara.h>
0010 #include <trackbase/TrkrDefs.h>
0011 
0012 #include <fun4all/SubsysReco.h>
0013 
0014 #include <TFile.h>
0015 #include <TH1.h>
0016 #include <TTree.h>
0017 
0018 #include <cmath>
0019 #include <iostream>
0020 #include <limits>
0021 #include <string>
0022 
0023 class TrkrCluster;
0024 class PHCompositeNode;
0025 class ActsGeometry;
0026 class SvtxTrack;
0027 class TrackSeed;
0028 class TrkrClusterContainer;
0029 class TrkrHitSetContainer;
0030 class PHG4TpcCylinderGeomContainer;
0031 class PHG4CylinderGeomContainer;
0032 class TpcDistortionCorrectionContainer;
0033 class TrackResiduals : public SubsysReco
0034 {
0035  public:
0036   TrackResiduals(const std::string &name = "TrackResiduals");
0037 
0038   ~TrackResiduals() override = default;
0039 
0040   int InitRun(PHCompositeNode *topNode) override;
0041   int process_event(PHCompositeNode *topNode) override;
0042   int End(PHCompositeNode *topNode) override;
0043   void outfileName(const std::string &name) { m_outfileName = name; }
0044   void alignment(bool align) { m_doAlignment = align; }
0045   void alignmentmapName(const std::string &name) { m_alignmentMapName = name; }
0046   void trackmapName(const std::string &name) { m_trackMapName = name; }
0047   void clusterTree() { m_doClusters = true; }
0048   void vertexTree() { m_doVertex = true; }
0049   void hitTree() { m_doHits = true; }
0050   void eventTree() { m_doEventTree = true; }
0051   void MatchedTracksOnly() { m_doMatchedOnly = true; }
0052   void ppmode() { m_ppmode = true; }
0053   void convertSeeds(bool flag) { m_convertSeeds = flag; }
0054   void dropClustersNoState(bool flag) { m_dropClustersNoState = flag; }
0055   void zeroField() { m_zeroField = true; }
0056   void runnumber(const int run) { m_runnumber = run; }
0057   void segment(const int seg) { m_segment = seg; }
0058   void linefitAll() { m_linefitTPCOnly = false; }
0059   void setClusterMinSize(unsigned int size) { m_min_cluster_size = size; }
0060   void failedTree() { m_doFailedSeeds = true; }
0061   void setSegment(const int segment) { m_segment = segment; }
0062 
0063   void set_doMicromegasOnly(bool value) { m_doMicromegasOnly = value; }
0064   void setTrkrClusterContainerName(std::string &name) { m_clusterContainerName = name; }
0065 
0066   void set_use_clustermover(bool flag) { m_use_clustermover = flag; }
0067 
0068  private:
0069   void fillStatesWithLineFit(const TrkrDefs::cluskey &ckey,
0070                              TrkrCluster *cluster, ActsGeometry *geometry);
0071   void clearClusterStateVectors();
0072   void createBranches();
0073   static float convertTimeToZ(ActsGeometry *geometry, TrkrDefs::cluskey cluster_key, TrkrCluster *cluster);
0074   void fillEventTree(PHCompositeNode *topNode);
0075   void fillClusterTree(TrkrClusterContainer *clusters, ActsGeometry *geometry);
0076   void fillHitTree(TrkrHitSetContainer *hitmap, ActsGeometry *geometry,
0077                    PHG4TpcCylinderGeomContainer *tpcGeom, PHG4CylinderGeomContainer *mvtxGeom,
0078                    PHG4CylinderGeomContainer *inttGeom, PHG4CylinderGeomContainer *mmGeom);
0079   void fillResidualTreeKF(PHCompositeNode *topNode);
0080   void fillResidualTreeSeeds(PHCompositeNode *topNode);
0081   void fillClusterBranchesKF(TrkrDefs::cluskey ckey, SvtxTrack *track,
0082                              const std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> &global_moved,
0083                              PHCompositeNode *topNode);
0084   void fillClusterBranchesSeeds(TrkrDefs::cluskey ckey,  // SvtxTrack* track,
0085                                 const std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> &global,
0086                                 PHCompositeNode *topNode);
0087   void lineFitClusters(std::vector<TrkrDefs::cluskey> &keys, TrkrClusterContainer *clusters, const short int &crossing);
0088   void circleFitClusters(std::vector<TrkrDefs::cluskey> &keys, TrkrClusterContainer *clusters, const short int &crossing);
0089   void fillStatesWithCircleFit(const TrkrDefs::cluskey &key, TrkrCluster *cluster,
0090                                Acts::Vector3 &glob, ActsGeometry *geometry);
0091   void fillVertexTree(PHCompositeNode *topNode);
0092   void fillFailedSeedTree(PHCompositeNode *topNode, std::set<unsigned int> &tpc_seed_ids);
0093   static float calc_dedx(TrackSeed *tpcseed, TrkrClusterContainer *clustermap, PHG4TpcCylinderGeomContainer *tpcGeom);
0094 
0095   bool m_use_clustermover = true;
0096 
0097   std::string m_outfileName = "";
0098   TFile *m_outfile = nullptr;
0099   TTree *m_tree = nullptr;
0100   TTree *m_clustree = nullptr;
0101   TTree *m_eventtree = nullptr;
0102   TTree *m_hittree = nullptr;
0103   TTree *m_vertextree = nullptr;
0104   TTree *m_failedfits = nullptr;
0105 
0106   bool m_doVertex = false;
0107   bool m_doClusters = false;
0108   bool m_doHits = false;
0109   bool m_doEventTree = false;
0110   bool m_zeroField = false;
0111   bool m_doFailedSeeds = false;
0112   bool m_doMatchedOnly = false;
0113 
0114   TpcClusterMover m_clusterMover;
0115   TpcGlobalPositionWrapper m_globalPositionWrapper;
0116 
0117   ClusterErrorPara m_clusErrPara;
0118   std::string m_alignmentMapName = "SvtxAlignmentStateMap";
0119   std::string m_trackMapName = "SvtxTrackMap";
0120   std::string m_clusterContainerName = "TRKR_CLUSTER";
0121 
0122   bool m_doAlignment = false;
0123   bool m_ppmode = false;
0124   bool m_convertSeeds = false;
0125   bool m_linefitTPCOnly = true;
0126   bool m_dropClustersNoState = false;
0127   unsigned int m_min_cluster_size = 0;
0128 
0129   bool m_doMicromegasOnly = false;
0130 
0131   int m_event = 0;
0132   int m_segment = std::numeric_limits<int>::quiet_NaN();
0133   int m_runnumber = std::numeric_limits<int>::quiet_NaN();
0134   int m_ntpcclus = std::numeric_limits<int>::quiet_NaN();
0135   float m_totalmbd = std::numeric_limits<float>::quiet_NaN();
0136   std::vector<int> m_firedTriggers;
0137   uint64_t m_gl1BunchCrossing = std::numeric_limits<uint64_t>::quiet_NaN();
0138 
0139   //! Event level quantities
0140   int m_nmvtx_all = std::numeric_limits<int>::quiet_NaN();
0141   int m_nintt_all = std::numeric_limits<int>::quiet_NaN();
0142   int m_ntpc_hits0 = std::numeric_limits<int>::quiet_NaN();
0143   int m_ntpc_hits1 = std::numeric_limits<int>::quiet_NaN();
0144   int m_ntpc_clus0 = std::numeric_limits<int>::quiet_NaN();
0145   int m_ntpc_clus1 = std::numeric_limits<int>::quiet_NaN();
0146   int m_nmms_all = std::numeric_limits<int>::quiet_NaN();
0147   int m_nsiseed = std::numeric_limits<int>::quiet_NaN();
0148   int m_ntpcseed = std::numeric_limits<int>::quiet_NaN();
0149   int m_ntracks_all = std::numeric_limits<int>::quiet_NaN();
0150 
0151   //! Track level quantities
0152   uint64_t m_bco = std::numeric_limits<uint64_t>::quiet_NaN();
0153   uint64_t m_bcotr = std::numeric_limits<uint64_t>::quiet_NaN();
0154   unsigned int m_trackid = std::numeric_limits<unsigned int>::quiet_NaN();
0155   int m_crossing = std::numeric_limits<int>::quiet_NaN();
0156   int m_crossing_estimate = std::numeric_limits<int>::quiet_NaN();
0157   unsigned int m_tpcid = std::numeric_limits<unsigned int>::quiet_NaN();
0158   unsigned int m_silid = std::numeric_limits<unsigned int>::quiet_NaN();
0159   float m_px = std::numeric_limits<float>::quiet_NaN();
0160   float m_py = std::numeric_limits<float>::quiet_NaN();
0161   float m_pz = std::numeric_limits<float>::quiet_NaN();
0162   float m_pt = std::numeric_limits<float>::quiet_NaN();
0163   float m_eta = std::numeric_limits<float>::quiet_NaN();
0164   float m_phi = std::numeric_limits<float>::quiet_NaN();
0165   float m_deltapt = std::numeric_limits<float>::quiet_NaN();
0166   int m_charge = std::numeric_limits<int>::quiet_NaN();
0167   float m_quality = std::numeric_limits<float>::quiet_NaN();
0168   float m_chisq = std::numeric_limits<float>::quiet_NaN();
0169   float m_ndf = std::numeric_limits<float>::quiet_NaN();
0170   int m_nhits = std::numeric_limits<int>::quiet_NaN();
0171   int m_nmaps = std::numeric_limits<int>::quiet_NaN();
0172   int m_nmapsstate = std::numeric_limits<int>::quiet_NaN();
0173   int m_nintt = std::numeric_limits<int>::quiet_NaN();
0174   int m_ninttstate = std::numeric_limits<int>::quiet_NaN();
0175   int m_ntpc = std::numeric_limits<int>::quiet_NaN();
0176   int m_ntpcstate = std::numeric_limits<int>::quiet_NaN();
0177   int m_nmms = std::numeric_limits<int>::quiet_NaN();
0178   int m_nmmsstate = std::numeric_limits<int>::quiet_NaN();
0179   unsigned int m_vertexid = std::numeric_limits<unsigned int>::quiet_NaN();
0180   int m_vertex_crossing = std::numeric_limits<int>::quiet_NaN();
0181   float m_vx = std::numeric_limits<float>::quiet_NaN();
0182   float m_vy = std::numeric_limits<float>::quiet_NaN();
0183   float m_vz = std::numeric_limits<float>::quiet_NaN();
0184   int m_vertex_ntracks = std::numeric_limits<int>::quiet_NaN();
0185   float m_pcax = std::numeric_limits<float>::quiet_NaN();
0186   float m_pcay = std::numeric_limits<float>::quiet_NaN();
0187   float m_pcaz = std::numeric_limits<float>::quiet_NaN();
0188   float m_rzslope = std::numeric_limits<float>::quiet_NaN();
0189   float m_rzint = std::numeric_limits<float>::quiet_NaN();
0190   float m_xyslope = std::numeric_limits<float>::quiet_NaN();
0191   float m_xyint = std::numeric_limits<float>::quiet_NaN();
0192   float m_yzslope = std::numeric_limits<float>::quiet_NaN();
0193   float m_yzint = std::numeric_limits<float>::quiet_NaN();
0194   float m_R = std::numeric_limits<float>::quiet_NaN();
0195   float m_X0 = std::numeric_limits<float>::quiet_NaN();
0196   float m_Y0 = std::numeric_limits<float>::quiet_NaN();
0197   float m_dcaxy = std::numeric_limits<float>::quiet_NaN();
0198   float m_dcaz = std::numeric_limits<float>::quiet_NaN();
0199   float m_tracklength = std::numeric_limits<float>::quiet_NaN();
0200 
0201   float m_silseedx = std::numeric_limits<float>::quiet_NaN();
0202   float m_silseedy = std::numeric_limits<float>::quiet_NaN();
0203   float m_silseedz = std::numeric_limits<float>::quiet_NaN();
0204   float m_silseedpx = std::numeric_limits<float>::quiet_NaN();
0205   float m_silseedpy = std::numeric_limits<float>::quiet_NaN();
0206   float m_silseedpz = std::numeric_limits<float>::quiet_NaN();
0207   int m_silseedcharge = std::numeric_limits<int>::quiet_NaN();
0208   float m_silseedphi = std::numeric_limits<float>::quiet_NaN();
0209   float m_silseedeta = std::numeric_limits<float>::quiet_NaN();
0210   float m_tpcseedx = std::numeric_limits<float>::quiet_NaN();
0211   float m_tpcseedy = std::numeric_limits<float>::quiet_NaN();
0212   float m_tpcseedz = std::numeric_limits<float>::quiet_NaN();
0213   float m_tpcseedpx = std::numeric_limits<float>::quiet_NaN();
0214   float m_tpcseedpy = std::numeric_limits<float>::quiet_NaN();
0215   float m_tpcseedpz = std::numeric_limits<float>::quiet_NaN();
0216   int m_tpcseedcharge = std::numeric_limits<int>::quiet_NaN();
0217   float m_tpcseedphi = std::numeric_limits<float>::quiet_NaN();
0218   float m_tpcseedeta = std::numeric_limits<float>::quiet_NaN();
0219 
0220   float m_dedx = std::numeric_limits<float>::quiet_NaN();
0221 
0222   //! hit tree info
0223   uint32_t m_hitsetkey = std::numeric_limits<uint32_t>::quiet_NaN();
0224   float m_hitgx = std::numeric_limits<float>::quiet_NaN();
0225   float m_hitgy = std::numeric_limits<float>::quiet_NaN();
0226   float m_hitgz = std::numeric_limits<float>::quiet_NaN();
0227   int m_hitlayer = std::numeric_limits<int>::quiet_NaN();
0228   int m_sector = std::numeric_limits<int>::quiet_NaN();
0229   int m_hitpad = std::numeric_limits<int>::quiet_NaN();
0230   int m_hittbin = std::numeric_limits<int>::quiet_NaN();
0231   int m_col = std::numeric_limits<int>::quiet_NaN();
0232   int m_row = std::numeric_limits<int>::quiet_NaN();
0233   int m_strip = std::numeric_limits<int>::quiet_NaN();
0234   float m_zdriftlength = std::numeric_limits<float>::quiet_NaN();
0235   
0236   float m_mbdvtxz = std::numeric_limits<float>::quiet_NaN();
0237 
0238   int m_ntracks = std::numeric_limits<int>::quiet_NaN();
0239   int m_nvertices = std::numeric_limits<int>::quiet_NaN();
0240 
0241   //! cluster tree info
0242   float m_sclusgr = std::numeric_limits<float>::quiet_NaN();
0243   float m_sclusphi = std::numeric_limits<float>::quiet_NaN();
0244   float m_scluseta = std::numeric_limits<float>::quiet_NaN();
0245   float m_adc = std::numeric_limits<float>::quiet_NaN();
0246   float m_clusmaxadc = std::numeric_limits<float>::quiet_NaN();
0247   int m_phisize = std::numeric_limits<int>::quiet_NaN();
0248   int m_zsize = std::numeric_limits<int>::quiet_NaN();
0249   float m_scluslx = std::numeric_limits<float>::quiet_NaN();
0250   float m_scluslz = std::numeric_limits<float>::quiet_NaN();
0251   float m_sclusgx = std::numeric_limits<float>::quiet_NaN();
0252   float m_sclusgy = std::numeric_limits<float>::quiet_NaN();
0253   float m_sclusgz = std::numeric_limits<float>::quiet_NaN();
0254   int m_scluslayer = std::numeric_limits<int>::quiet_NaN();
0255   float m_scluselx = std::numeric_limits<float>::quiet_NaN();
0256   float m_scluselz = std::numeric_limits<float>::quiet_NaN();
0257   int m_clussector = std::numeric_limits<int>::quiet_NaN();
0258   int m_side = std::numeric_limits<int>::quiet_NaN();
0259   int m_staveid = std::numeric_limits<int>::quiet_NaN();
0260   int m_chipid = std::numeric_limits<int>::quiet_NaN();
0261   int m_strobeid = std::numeric_limits<int>::quiet_NaN();
0262   int m_ladderzid = std::numeric_limits<int>::quiet_NaN();
0263   int m_ladderphiid = std::numeric_limits<int>::quiet_NaN();
0264   int m_timebucket = std::numeric_limits<int>::quiet_NaN();
0265   int m_segtype = std::numeric_limits<int>::quiet_NaN();
0266   int m_tileid = std::numeric_limits<int>::quiet_NaN();
0267 
0268   //! clusters on track information
0269   std::vector<float> m_clusAdc;
0270   std::vector<float> m_clusMaxAdc;
0271   std::vector<float> m_cluslx;
0272   std::vector<float> m_cluslz;
0273   std::vector<float> m_cluselx;
0274   std::vector<float> m_cluselz;
0275   std::vector<float> m_clusgx;
0276   std::vector<float> m_clusgy;
0277   std::vector<float> m_clusgz;
0278   std::vector<float> m_clusgxunmoved;
0279   std::vector<float> m_clusgyunmoved;
0280   std::vector<float> m_clusgzunmoved;
0281   std::vector<float> m_clusgr;
0282   std::vector<int> m_clsector;
0283   std::vector<int> m_clside;
0284   std::vector<int> m_cluslayer;
0285   std::vector<int> m_clusphisize;
0286   std::vector<int> m_cluszsize;
0287   std::vector<int> m_clusedge;
0288   std::vector<int> m_clusoverlap;
0289   std::vector<uint64_t> m_cluskeys;
0290   std::vector<float> m_idealsurfcenterx;
0291   std::vector<float> m_idealsurfcentery;
0292   std::vector<float> m_idealsurfcenterz;
0293   std::vector<float> m_idealsurfnormx;
0294   std::vector<float> m_idealsurfnormy;
0295   std::vector<float> m_idealsurfnormz;
0296   std::vector<float> m_missurfcenterx;
0297   std::vector<float> m_missurfcentery;
0298   std::vector<float> m_missurfcenterz;
0299   std::vector<float> m_missurfnormx;
0300   std::vector<float> m_missurfnormy;
0301   std::vector<float> m_missurfnormz;
0302   std::vector<float> m_clusgxideal;
0303   std::vector<float> m_clusgyideal;
0304   std::vector<float> m_clusgzideal;
0305   std::vector<float> m_idealsurfalpha;
0306   std::vector<float> m_idealsurfbeta;
0307   std::vector<float> m_idealsurfgamma;
0308   std::vector<float> m_missurfalpha;
0309   std::vector<float> m_missurfbeta;
0310   std::vector<float> m_missurfgamma;
0311 
0312   //! states on track information
0313   std::vector<float> m_statelx;
0314   std::vector<float> m_statelz;
0315   std::vector<float> m_stateelx;
0316   std::vector<float> m_stateelz;
0317   std::vector<float> m_stategx;
0318   std::vector<float> m_stategy;
0319   std::vector<float> m_stategz;
0320   std::vector<float> m_statepx;
0321   std::vector<float> m_statepy;
0322   std::vector<float> m_statepz;
0323   std::vector<float> m_statepl;
0324 
0325   std::vector<float> m_statelxglobderivdx;
0326   std::vector<float> m_statelxglobderivdy;
0327   std::vector<float> m_statelxglobderivdz;
0328   std::vector<float> m_statelxglobderivdalpha;
0329   std::vector<float> m_statelxglobderivdbeta;
0330   std::vector<float> m_statelxglobderivdgamma;
0331 
0332   std::vector<float> m_statelxlocderivd0;
0333   std::vector<float> m_statelxlocderivz0;
0334   std::vector<float> m_statelxlocderivphi;
0335   std::vector<float> m_statelxlocderivtheta;
0336   std::vector<float> m_statelxlocderivqop;
0337 
0338   std::vector<float> m_statelzglobderivdx;
0339   std::vector<float> m_statelzglobderivdy;
0340   std::vector<float> m_statelzglobderivdz;
0341   std::vector<float> m_statelzglobderivdalpha;
0342   std::vector<float> m_statelzglobderivdbeta;
0343   std::vector<float> m_statelzglobderivdgamma;
0344 
0345   std::vector<float> m_statelzlocderivd0;
0346   std::vector<float> m_statelzlocderivz0;
0347   std::vector<float> m_statelzlocderivphi;
0348   std::vector<float> m_statelzlocderivtheta;
0349   std::vector<float> m_statelzlocderivqop;
0350 };
0351 
0352 #endif  // TRACKRESIDUALS_H