SCorrelatorQAMaker/src/SCheckTrackPairs.cc

0001 // ----------------------------------------------------------------------------
0002 // 'SCheckTrackPairs.cc'
0003 // Derek Anderson
0004 // 102.21.2024
0005 //
0006 // SCorrelatorQAMaker plugin to iterate through
0007 // all pairs of tracks in an event and fill
0008 // tuples/histograms comparing them.
0009 // ----------------------------------------------------------------------------
0010
0011 #define SCORRELATORQAMAKER_SCHECKTRACKPAIRS_CC
0012
0013 // plugin definition
0014 #include "SCheckTrackPairs.h"
0015
0016 // make common namespaces implicit
0017 using namespace std;
0018 using namespace findNode;
0019
0020
0021
0022 namespace SColdQcdCorrelatorAnalysis {
0023
0024   // SCheckTrackPairs public methods ------------------------------------------
0025
0026   int SCheckTrackPairs::Init(PHCompositeNode* topNode) {
0027
0028     InitOutput();
0029     InitTuples();
0030     return Fun4AllReturnCodes::EVENT_OK;
0031
0032   }  // end 'Init(PHCompositeNode*)'
0033
0034
0035
0036   int SCheckTrackPairs::process_event(PHCompositeNode* topNode) {
0037
0038     ResetVectors();
0039     DoDoubleTrackLoop(topNode);
0040     return Fun4AllReturnCodes::EVENT_OK;
0041
0042   }  // end 'process_event(PHCompositeNode* topNode)'
0043
0044
0045
0046   int SCheckTrackPairs::End(PHCompositeNode* topNode) {
0047
0048     SaveOutput();
0049     CloseOutput();
0050     return Fun4AllReturnCodes::EVENT_OK;
0051
0052   }  // end 'End(PHCompositeNode*)'
0053
0054
0055
0056   // SCheckTrackPairs internal methods ----------------------------------------
0057
0058   void SCheckTrackPairs::InitTuples() {
0059
0060     if (m_isDebugOn && (m_verbosity > 2)) {
0061       cout << "SColdQcdCorrelatorAnalysis::SCheckTrackPairs::InitTuples(): initializing output tuple." << endl;
0062     }
0063
0064     // create leaf lists and add leaves
0065     vector<string> vecTrkPairLeavesA = Types::TrkInfo::GetListOfMembers();
0066     vector<string> vecTrkPairLeavesB = Types::TrkInfo::GetListOfMembers();
0067     Interfaces::AddTagToLeaves("_a", vecTrkPairLeavesA);
0068     Interfaces::AddTagToLeaves("_b", vecTrkPairLeavesB);
0069     vecTrkPairLeavesA.push_back("nClustKey_a");
0070     vecTrkPairLeavesB.push_back("nClustKey_b");
0071
0072     // combine lists
0073     vector<string> vecTrkPairLeaves;
0074     Interfaces::CombineLeafLists(vecTrkPairLeavesA, vecTrkPairLeaves);
0075     Interfaces::CombineLeafLists(vecTrkPairLeavesB, vecTrkPairLeaves);
0076
0077     // add leaves for no. of same cluster keys b/n pair and
0078     // distance between pair
0079     vecTrkPairLeaves.push_back("nSameClustKey");
0080     vecTrkPairLeaves.push_back("trackDeltaR");
0081
0082     // compress leaves into a colon-separated list
0083     string argTrkPairLeaves = Interfaces::FlattenLeafList(vecTrkPairLeaves);
0084
0085     // create tuple and return
0086     m_ntTrackPairs = new TNtuple("ntTrackPairs", "Pairs of tracks", argTrkPairLeaves.data());
0087     m_vecTrackPairLeaves.reserve(vecTrkPairLeaves.size());
0088     return;
0089
0090   }  // end 'InitTuples()'
0091
0092
0093
0094   void SCheckTrackPairs::SaveOutput() {
0095
0096     if (m_isDebugOn && (m_verbosity > 2)) {
0097       cout << "SColdQcdCorrelatorAnalysis::SCheckTrackPairs::SaveOutput(): saving output." << endl;
0098     }
0099
0100     m_outDir       -> cd();
0101     m_ntTrackPairs -> Write();
0102     return;
0103
0104   }  // end 'SaveOutput()'
0105
0106
0107
0108   void SCheckTrackPairs::ResetVectors() {
0109
0110     if (m_isDebugOn && (m_verbosity > 4)) {
0111       cout << "SColdQcdCorrelatorAnalysis::SCheckTrackPairs::ResetVectors(): resetting vectors." << endl;
0112     }
0113
0114     // set leaf values to a default
0115     const size_t nLeaves = m_vecTrackPairLeaves.size();
0116     for (size_t iLeaf = 0; iLeaf < nLeaves; iLeaf++) {
0117       m_vecTrackPairLeaves[iLeaf] = -999.;
0118     }
0119
0120     // clear cluster keys
0121     m_vecClustKeysA.clear();
0122     m_vecClustKeysB.clear();
0123     return;
0124
0125   }  // end 'ResetVectors()'
0126
0127
0128
0129   void SCheckTrackPairs::DoDoubleTrackLoop(PHCompositeNode* topNode) {
0130
0131     if (m_isDebugOn && (m_verbosity > 2)) {
0132       cout << "SColdQcdCorrelatorAnalysis::SCheckTrackPairs::DoDoubleTrackLoop(): looping over all pairs of tracks." << endl;
0133     }
0134
0135     // loop over tracks
0136     SvtxTrack*    trackA  = NULL;
0137     SvtxTrack*    trackB  = NULL;
0138     SvtxTrackMap* mapTrks = Interfaces::GetTrackMap(topNode);
0139     for (
0140       SvtxTrackMap::Iter itTrkA = mapTrks -> begin();
0141       itTrkA != mapTrks -> end();
0142       ++itTrkA
0143     ) {
0144
0145       // grab track A and skip if bad
0146       trackA = itTrkA -> second;
0147       if (!trackA) continue;
0148
0149       const bool isGoodTrackA = IsGoodTrack(trackA, topNode);
0150       if (!isGoodTrackA) continue;
0151
0152       // grab track info
0153       Types::TrkInfo trkInfoA(trackA, topNode);
0154
0155       // loop over tracks again
0156       for (SvtxTrackMap::Iter itTrkB = mapTrks -> begin(); itTrkB != mapTrks -> end(); ++itTrkB) {
0157
0158         // grab track B and skip if bad
0159         trackB = itTrkB -> second;
0160         if (!trackB) continue;
0161
0162         const bool isGoodTrackB = IsGoodTrack(trackB, topNode);
0163         if (!isGoodTrackB) continue;
0164
0165         // skip if same as track A
0166         if ((trackA -> get_id()) == (trackB -> get_id())) continue;
0167
0168         // grab track info
0169         Types::TrkInfo trkInfoB(trackB, topNode);
0170
0171         // calculate delta-R
0172         const double dfTrkAB = trkInfoA.GetPhi() - trkInfoB.GetPhi();
0173         const double dhTrkAB = trkInfoA.GetEta() - trkInfoB.GetEta();
0174         const double drTrkAB = sqrt((dfTrkAB * dfTrkAB) + (dhTrkAB * dhTrkAB));
0175
0176         // clear vectors for checking cluster keys
0177         m_vecClustKeysA.clear();
0178         m_vecClustKeysB.clear();
0179
0180         // loop over clusters from track A
0181         auto seedTpcA = trackA -> get_tpc_seed();
0182         if (seedTpcA) {
0183           for (auto local_iterA = seedTpcA -> begin_cluster_keys(); local_iterA != seedTpcA -> end_cluster_keys(); ++local_iterA) {
0184             TrkrDefs::cluskey cluster_keyA = *local_iterA;
0185             m_vecClustKeysA.push_back(cluster_keyA);
0186           }
0187         }
0188
0189         // loop over clusters from track B
0190         auto seedTpcB = trackB -> get_tpc_seed();
0191         if (seedTpcB) {
0192           for (auto local_iterB = seedTpcB -> begin_cluster_keys(); local_iterB != seedTpcB -> end_cluster_keys(); ++local_iterB) {
0193             TrkrDefs::cluskey cluster_keyB = *local_iterB;
0194             m_vecClustKeysB.push_back(cluster_keyB);
0195           }
0196         }
0197
0198         // calculate no. of same cluster keys
0199         uint64_t nSameKey = 0;
0200         for (auto keyA : m_vecClustKeysA) {
0201           for (auto keyB : m_vecClustKeysB) {
0202             if (keyA == keyB) {
0203               ++nSameKey;
0204               break;
0205             }
0206           }  // end cluster key B loop
0207         }  // end cluster key A loop
0208
0209         // set tuple leaves
0210         m_vecTrackPairLeaves[0]  = (float) trkInfoA.GetID();
0211         m_vecTrackPairLeaves[1]  = (float) trkInfoA.GetNMvtxLayer();
0212         m_vecTrackPairLeaves[2]  = (float) trkInfoA.GetNInttLayer();
0213         m_vecTrackPairLeaves[3]  = (float) trkInfoA.GetNTpcLayer();
0214         m_vecTrackPairLeaves[4]  = (float) trkInfoA.GetNMvtxClust();
0215         m_vecTrackPairLeaves[5]  = (float) trkInfoA.GetNInttClust();
0216         m_vecTrackPairLeaves[6]  = (float) trkInfoA.GetNTpcClust();
0217         m_vecTrackPairLeaves[7]  = (float) trkInfoA.GetEta();
0218         m_vecTrackPairLeaves[8]  = (float) trkInfoA.GetPhi();
0219         m_vecTrackPairLeaves[9]  = (float) trkInfoA.GetPX();
0220         m_vecTrackPairLeaves[10] = (float) trkInfoA.GetPY();
0221         m_vecTrackPairLeaves[11] = (float) trkInfoA.GetPZ();
0222         m_vecTrackPairLeaves[12] = (float) trkInfoA.GetPT();
0223         m_vecTrackPairLeaves[13] = (float) trkInfoA.GetEne();
0224         m_vecTrackPairLeaves[14] = (float) trkInfoA.GetDcaXY();
0225         m_vecTrackPairLeaves[15] = (float) trkInfoA.GetDcaZ();
0226         m_vecTrackPairLeaves[16] = (float) trkInfoA.GetPtErr();
0227         m_vecTrackPairLeaves[17] = (float) trkInfoA.GetQuality();
0228         m_vecTrackPairLeaves[18] = (float) trkInfoA.GetVX();
0229         m_vecTrackPairLeaves[19] = (float) trkInfoA.GetVY();
0230         m_vecTrackPairLeaves[20] = (float) trkInfoA.GetVZ();
0231         m_vecTrackPairLeaves[21] = (float) m_vecClustKeysA.size();
0232         m_vecTrackPairLeaves[22] = (float) trkInfoB.GetID();
0233         m_vecTrackPairLeaves[23] = (float) trkInfoB.GetNMvtxLayer();
0234         m_vecTrackPairLeaves[24] = (float) trkInfoB.GetNInttLayer();
0235         m_vecTrackPairLeaves[25] = (float) trkInfoB.GetNTpcLayer();
0236         m_vecTrackPairLeaves[26] = (float) trkInfoB.GetNMvtxClust();
0237         m_vecTrackPairLeaves[27] = (float) trkInfoB.GetNInttClust();
0238         m_vecTrackPairLeaves[28] = (float) trkInfoB.GetNTpcClust();
0239         m_vecTrackPairLeaves[29] = (float) trkInfoB.GetEta();
0240         m_vecTrackPairLeaves[30] = (float) trkInfoB.GetPhi();
0241         m_vecTrackPairLeaves[31] = (float) trkInfoB.GetPX();
0242         m_vecTrackPairLeaves[32] = (float) trkInfoB.GetPY();
0243         m_vecTrackPairLeaves[33] = (float) trkInfoB.GetPZ();
0244         m_vecTrackPairLeaves[34] = (float) trkInfoB.GetPT();
0245         m_vecTrackPairLeaves[35] = (float) trkInfoB.GetEne();
0246         m_vecTrackPairLeaves[36] = (float) trkInfoB.GetDcaXY();
0247         m_vecTrackPairLeaves[37] = (float) trkInfoB.GetDcaZ();
0248         m_vecTrackPairLeaves[38] = (float) trkInfoB.GetPtErr();
0249         m_vecTrackPairLeaves[39] = (float) trkInfoB.GetQuality();
0250         m_vecTrackPairLeaves[40] = (float) trkInfoB.GetVX();
0251         m_vecTrackPairLeaves[41] = (float) trkInfoB.GetVY();
0252         m_vecTrackPairLeaves[42] = (float) trkInfoB.GetVZ();
0253         m_vecTrackPairLeaves[43] = (float) m_vecClustKeysB.size();
0254         m_vecTrackPairLeaves[44] = (float) nSameKey;
0255         m_vecTrackPairLeaves[45] = (float) drTrkAB;
0256
0257         // fill track pair tuple
0258         m_ntTrackPairs -> Fill(m_vecTrackPairLeaves.data());
0259
0260       }  // end 2nd track loop
0261     }  // end track loop
0262     return;
0263
0264   }  // end 'DoDoubleTrackLoop(PHCompositeNode*)'
0265
0266
0267
0268   bool SCheckTrackPairs::IsGoodTrack(SvtxTrack* track, PHCompositeNode* topNode) {
0269
0270     // print debug statement
0271     if (m_isDebugOn && (m_verbosity > 4)) {
0272       cout << "SCheckTrackPairs::IsGoodTrack(SvtxTrack* track, PHCompositeNode*) Checking if track is good..." << endl;
0273     }
0274
0275     // grab info
0276     Types::TrkInfo info(track, topNode);
0277
0278     // if needed, check if dca is in pt-dependent range
0279     bool isInDcaSigma = true;
0280     if (m_config.doDcaSigCut) {
0281       isInDcaSigma = info.IsInSigmaDcaCut(m_config.nSigCut, m_config.ptFitMax, m_config.fSigDca);
0282     }
0283
0284     // if needed, check if track is from primary vertex
0285     const bool isFromPrimVtx = m_config.useOnlyPrimVtx ? Tools::IsFromPrimaryVtx(track, topNode) : true;
0286
0287     // check if seed is good & track is in acceptance
0288     const bool isSeedGood = Tools::IsGoodTrackSeed(track, m_config.requireSiSeed);
0289     const bool isInAccept = info.IsInAcceptance(m_config.trkAccept);
0290
0291     // return overall goodness of track
0292     return (isFromPrimVtx && isInDcaSigma && isSeedGood && isInAccept);
0293
0294   }  // end 'IsGoodTrack(SvtxTrack*, PHCompositeNode* topNode)'
0295
0296 }  // end SColdQcdCorrelatorAnalysis namespace
0297
0298 // end ------------------------------------------------------------------------
0299