SCorrelatorUtilities/src/TrkTools.cc

0001 /// ---------------------------------------------------------------------------
0002 /*! \file   TrkTools.cc
0003  *  \author Derek Anderson
0004  *  \date   02.29.2024
0005  *
0006  *  Collection of frequent track-related methods utilized
0007  *  in the sPHENIX Cold QCD Energy-Energy Correlator analysis.
0008  */
0009 /// ---------------------------------------------------------------------------
0010
0011 #define SCORRELATORUTILITIES_TRKTOOLS_CC
0012
0013 // namespace definition
0014 #include "TrkTools.h"
0015
0016 // make common namespaces implicit
0017 using namespace std;
0018
0019
0020
0021 // track methods ==============================================================
0022
0023 namespace SColdQcdCorrelatorAnalysis {
0024
0025   // --------------------------------------------------------------------------
0026   //! Get no. of layers hit from a specific subsystem
0027   // --------------------------------------------------------------------------
0028   int Tools::GetNumLayer(SvtxTrack* track, const int16_t sys) {
0029
0030     // grab track seed
0031     TrackSeed* seed = GetTrackSeed(track, sys);
0032
0033     // set min no. of layers
0034     const int minInttLayer = Const::NMvtxLayer();
0035     const int minTpcLayer  = Const::NMvtxLayer() + Const::NInttLayer();
0036
0037     // initialize hit flags
0038     vector<bool> isMvtxLayerHit( Const::NMvtxLayer(), false );
0039     vector<bool> isInttLayerHit( Const::NInttLayer(), false );
0040     vector<bool> isTpcLayerHit(  Const::NTpcLayer(),  false  );
0041
0042     // determine which layers were hit
0043     int layer     = 0;
0044     int mvtxLayer = 0;
0045     int inttLayer = 0;
0046     int tpcLayer  = 0;
0047     for (
0048       auto itClustKey = (seed -> begin_cluster_keys());
0049       itClustKey != (seed -> end_cluster_keys());
0050       ++itClustKey
0051     ) {
0052
0053       // grab layer number
0054       layer = TrkrDefs::getLayer(*itClustKey);
0055
0056       // increment accordingly
0057       switch (sys) {
0058         case Const::Subsys::Mvtx:
0059           if (layer < Const::NMvtxLayer()) {
0060             mvtxLayer                 = layer;
0061             isMvtxLayerHit[mvtxLayer] = true;
0062           }
0063           break;
0064         case Const::Subsys::Intt:
0065           if ((layer >= minInttLayer) && (layer < minTpcLayer)) {
0066             inttLayer                 = layer - minInttLayer;
0067             isInttLayerHit[inttLayer] = true;
0068           }
0069           break;
0070         case Const::Subsys::Tpc:
0071           if (layer >= minTpcLayer) {
0072             tpcLayer                = layer - minTpcLayer;
0073             isTpcLayerHit[tpcLayer] = true;
0074           }
0075           break;
0076         default:
0077           break;
0078       }
0079     }  // end cluster loop
0080
0081     // get the relevant no. of layers
0082     int nLayer = 0;
0083     switch (sys) {
0084       case Const::Subsys::Mvtx:
0085         for (int iMvtxLayer = 0; iMvtxLayer < Const::NMvtxLayer(); iMvtxLayer++) {
0086           if (isMvtxLayerHit[iMvtxLayer]) ++nLayer;
0087         }
0088         break;
0089       case Const::Subsys::Intt:
0090         for (int iInttLayer = 0; iInttLayer < Const::NInttLayer(); iInttLayer++) {
0091           if (isInttLayerHit[iInttLayer]) ++nLayer;
0092         }
0093         break;
0094       case Const::Subsys::Tpc:
0095         for (int iTpcLayer = 0; iTpcLayer < Const::NTpcLayer(); iTpcLayer++) {
0096           if (isTpcLayerHit[iTpcLayer]) ++nLayer;
0097         }
0098         break;
0099       default:
0100         break;
0101     }
0102     return nLayer;
0103
0104   }  // end 'GetNumLayer(SvtxTrack*, int16_t)'
0105
0106
0107
0108   // --------------------------------------------------------------------------
0109   //! Get no. of clusters from a specific subsystem
0110   // --------------------------------------------------------------------------
0111   int Tools::GetNumClust(SvtxTrack* track, const int16_t sys) {
0112
0113     // grab track seed
0114     TrackSeed* seed = GetTrackSeed(track, sys);
0115
0116     // set min no. of layers
0117     const int minInttLayer = Const::NMvtxLayer();
0118     const int minTpcLayer  = Const::NMvtxLayer() + Const::NInttLayer();
0119
0120     // determine no. of clusters for a given layer
0121     int layer    = 0;
0122     int nCluster = 0;
0123     for (auto itClustKey = (seed -> begin_cluster_keys()); itClustKey != (seed -> end_cluster_keys()); ++itClustKey) {
0124
0125       // grab layer number
0126       layer = TrkrDefs::getLayer(*itClustKey);
0127
0128       // increment accordingly
0129       switch (sys) {
0130         case Const::Subsys::Mvtx:
0131           if (layer < Const::NMvtxLayer()) {
0132             ++nCluster;
0133           }
0134           break;
0135         case Const::Subsys::Intt:
0136           if ((layer >= minInttLayer) && (layer < minTpcLayer)) {
0137             ++nCluster;
0138           }
0139           break;
0140         case Const::Subsys::Tpc:
0141           if (layer >= minTpcLayer) {
0142             ++nCluster;
0143           }
0144           break;
0145         default:
0146           break;
0147       }
0148     }  // end cluster loop
0149     return nCluster;
0150
0151   }  // end 'GetNumClust(SvtxTrack*, int16_t)'
0152
0153
0154
0155   // --------------------------------------------------------------------------
0156   //! Get barcode of matching truth particle
0157   // --------------------------------------------------------------------------
0158   int Tools::GetMatchID(SvtxTrack* track, SvtxTrackEval* trackEval) {
0159
0160     // get best match from truth particles
0161     PHG4Particle* bestMatch = trackEval -> max_truth_particle_by_nclusters(track);
0162
0163     // grab barcode of best match
0164     int matchID;
0165     if (bestMatch) {
0166       matchID = bestMatch -> get_barcode();
0167     } else {
0168       matchID = numeric_limits<int>::max();
0169     }
0170     return matchID;
0171
0172   }  // end 'GetMatchID(SvtxTrack*)'
0173
0174
0175
0176   // --------------------------------------------------------------------------
0177   //! Does track have the relevant seed?
0178   // --------------------------------------------------------------------------
0179   bool Tools::IsGoodTrackSeed(SvtxTrack* track, const bool requireSiSeeds) {
0180
0181     // get track seeds
0182     TrackSeed* trkSiSeed  = track -> get_silicon_seed();
0183     TrackSeed* trkTpcSeed = track -> get_tpc_seed();
0184
0185     // check if one or both seeds are present as needed
0186     bool isSeedGood = (trkSiSeed && trkTpcSeed);
0187     if (!requireSiSeeds) {
0188       isSeedGood = (trkSiSeed || trkTpcSeed);
0189     }
0190     return isSeedGood;
0191
0192   }  // end 'IsGoodSeed(SvtxTrack*)'
0193
0194
0195
0196   // --------------------------------------------------------------------------
0197   //! Is track connected to the primary vertex?
0198   // --------------------------------------------------------------------------
0199   bool Tools::IsFromPrimaryVtx(SvtxTrack* track, PHCompositeNode* topNode) {
0200
0201     // get id of vertex associated with track
0202     const int vtxID = (int) track -> get_vertex_id();
0203
0204     // get id of primary vertex
0205     GlobalVertex* primVtx   = Interfaces::GetGlobalVertex(topNode);
0206     const int     primVtxID = primVtx -> get_id();
0207
0208     // check if from vertex and return
0209     const bool isFromPrimVtx = (vtxID == primVtxID);
0210     return isFromPrimVtx;
0211
0212   }  // end 'IsFromPrimaryVtx(SvtTrack*, PHCompsiteNode*)'
0213
0214
0215
0216   // --------------------------------------------------------------------------
0217   //! Get uncertainty on track pt from fit
0218   // --------------------------------------------------------------------------
0219   double Tools::GetTrackDeltaPt(SvtxTrack* track) {
0220
0221     // grab covariances
0222     const float trkCovXX = track -> get_error(3, 3);
0223     const float trkCovXY = track -> get_error(3, 4);
0224     const float trkCovYY = track -> get_error(4, 4);
0225
0226     // grab momentum
0227     const float trkPx = track -> get_px();
0228     const float trkPy = track -> get_py();
0229     const float trkPt = track -> get_pt();
0230
0231     // calculate delta-pt
0232     const float numer    = (trkCovXX * trkPx * trkPx) + 2 * (trkCovXY * trkPx * trkPy) + (trkCovYY * trkPy * trkPy);
0233     const float denom    = (trkPx * trkPx) + (trkPy * trkPy);
0234     const float ptDelta2 = numer / denom;
0235     const float ptDelta  = sqrt(ptDelta2) / trkPt;
0236     return ptDelta;
0237
0238   }  // end 'GetTrackDeltaPt(SvtxTrack*)'
0239
0240
0241
0242   // --------------------------------------------------------------------------
0243   //! Get track XY and Z DCA
0244   // --------------------------------------------------------------------------
0245   pair<double, double> Tools::GetTrackDcaPair(SvtxTrack* track, PHCompositeNode* topNode) {
0246
0247     // get global vertex and convert to acts vector
0248     GlobalVertex* sphxVtx = Interfaces::GetGlobalVertex(topNode);
0249     Acts::Vector3 actsVtx = Acts::Vector3(sphxVtx -> get_x(), sphxVtx -> get_y(), sphxVtx -> get_z());
0250
0251     // return dca
0252     const auto dcaAndErr = TrackAnalysisUtils::get_dca(track, actsVtx);
0253     return make_pair(dcaAndErr.first.first, dcaAndErr.second.first);
0254
0255   }  // end 'GetTrackDcaPair(SvtxTrack*, PHCompositeNode*)'
0256
0257
0258
0259   // --------------------------------------------------------------------------
0260   //! Get track seed
0261   // --------------------------------------------------------------------------
0262   TrackSeed* Tools::GetTrackSeed(SvtxTrack* track, const int16_t sys) {
0263
0264     // get both track seeds
0265     TrackSeed* trkSiSeed  = track -> get_silicon_seed();
0266     TrackSeed* trkTpcSeed = track -> get_tpc_seed();
0267
0268     // select relevant seed
0269     const bool hasBothSeeds   = (trkSiSeed  && trkTpcSeed);
0270     const bool hasOnlyTpcSeed = (!trkSiSeed && trkTpcSeed);
0271
0272     TrackSeed* seed = NULL;
0273     switch (sys) {
0274       case Const::Subsys::Mvtx:
0275         if (hasBothSeeds)   seed = trkSiSeed;
0276         if (hasOnlyTpcSeed) seed = trkTpcSeed;
0277         break;
0278       case Const::Subsys::Intt:
0279         if (hasBothSeeds)   seed = trkSiSeed;
0280         if (hasOnlyTpcSeed) seed = trkTpcSeed;
0281         break;
0282       case Const::Subsys::Tpc:
0283         seed = trkTpcSeed;
0284         break;
0285       default:
0286         seed = NULL;
0287         break;
0288     }
0289     return seed;
0290
0291   }  // end 'GetTrackSeed(SvtxTrack*, int16_t)'
0292
0293
0294
0295   // --------------------------------------------------------------------------
0296   //! Get track's associated vertex
0297   // --------------------------------------------------------------------------
0298   ROOT::Math::XYZVector Tools::GetTrackVertex(SvtxTrack* track, PHCompositeNode* topNode) {
0299
0300     // get vertex associated with track
0301     const int     vtxID = (int) track -> get_vertex_id();
0302     GlobalVertex* vtx   = Interfaces::GetGlobalVertex(topNode, vtxID);
0303
0304     // return vertex 3-vector
0305     ROOT::Math::XYZVector xyzVtx = ROOT::Math::XYZVector(vtx -> get_x(), vtx -> get_y(), vtx -> get_z());
0306     return xyzVtx;
0307
0308   }  // end 'GetTrackVertex(SvtxTrack*, PHCompositeNode*)'
0309
0310 }  // end SColdQcdCorrealtorAnalysis namespace
0311
0312 // end ------------------------------------------------------------------------