SCorrelatorQAMaker/src/SCheckCstPairs.cc

0001 // ----------------------------------------------------------------------------
0002 // 'SCheckCstPairs.cc'
0003 // Derek Anderson
0004 // 03.10.2024
0005 //
0006 // SCorrelatorQAMaker plugin to iterate through all pairs of constituents in
0007 // an event and fill tuples/histograms comparing them.
0008 //
0009 // This is similar to the `SCheckTrackPairs` plugin, which specifically looks
0010 // at pairs of tracks off the node tree. This plugin compares constituents
0011 // of any type off the correlator jet tree.
0012 // ----------------------------------------------------------------------------
0013
0014 #define SCORRELATORQAMAKER_SCHECKCSTPAIRS_CC
0015
0016 // plugin definition
0017 #include "SCheckCstPairs.h"
0018
0019 // make common namespaces implicit
0020 using namespace std;
0021
0022
0023
0024 namespace SColdQcdCorrelatorAnalysis {
0025
0026   // SCheckCstPairs public methods --------------------------------------------
0027
0028   void SCheckCstPairs::Init() {
0029
0030     // announce start
0031     cout << "\n Starting constituent pair check!" << endl;
0032
0033     // initialize output
0034     InitOutput();
0035
0036     // run internal routines
0037     InitInput();
0038     InitTree();
0039     InitHists();
0040     return;
0041
0042   }  // end 'Init()'
0043
0044
0045
0046   void SCheckCstPairs::Analyze() {
0047
0048     DoDoubleCstLoop();
0049     return;
0050
0051   }  // end 'Analyze()'
0052
0053
0054
0055   void SCheckCstPairs::End() {
0056
0057     // run internal routines
0058     SaveOutput();
0059     CloseInput();
0060
0061      // close output
0062     CloseOutput();
0063
0064     // announce end and exit
0065     cout << "  Finished constituent-pair check!\n" << endl;
0066     return;
0067
0068   }  // end 'End()'
0069
0070
0071
0072   // SCheckCstPairs internal methods ------------------------------------------
0073
0074   void SCheckCstPairs::InitInput() {
0075
0076
0077     // open input files
0078     m_fInput = new TFile(m_config.inFileName.data(), "read");
0079     if (!m_fInput) {
0080       cerr << "PANIC: couldn't open input file:\n"
0081            << "       \"" << m_config.inFileName << "\"\n"
0082            << endl;
0083       assert(m_fInput);
0084     }
0085     cout << "    Opened input file:\n"
0086          << "      \"" << m_config.inFileName << "\""
0087          << endl;
0088
0089     // grab input tree
0090     m_cInput = (TChain*) m_fInput -> Get(m_config.inChainName.data());
0091     if (!m_tInput) {
0092       cerr << "PANIC: couldn't grab input chain \"" << m_config.inChainName << "\"!\n" << endl;
0093       assert(m_cInput);
0094     }
0095     cout << "    Grabbed input chain \"" << m_config.inChainName << "\"." << endl;
0096
0097     // exit routine
0098     return;
0099
0100   }  // end 'InitInput()'
0101
0102
0103
0104   void SCheckCstPairs::InitTree() {
0105
0106     // set up input chain
0107     m_fCurrent = -1;
0108     m_cInput -> SetMakeClass(1);
0109
0110     // set truth vs. reco branch addresses
0111     if (m_config.isInChainTruth) {
0112       m_cInput -> SetBranchAddress("Parton3_ID",   &m_partonID.first,    &m_brPartonID.first);
0113       m_cInput -> SetBranchAddress("Parton4_ID",   &m_partonID.second,   &m_brPartonID.second);
0114       m_cInput -> SetBranchAddress("Parton3_MomX", &m_partonMomX.first,  &m_brPartonMomX.first);
0115       m_cInput -> SetBranchAddress("Parton3_MomY", &m_partonMomY.first,  &m_brPartonMomY.first);
0116       m_cInput -> SetBranchAddress("Parton3_MomZ", &m_partonMomZ.first,  &m_brPartonMomZ.first);
0117       m_cInput -> SetBranchAddress("Parton4_MomX", &m_partonMomX.second, &m_brPartonMomX.second);
0118       m_cInput -> SetBranchAddress("Parton4_MomY", &m_partonMomY.second, &m_brPartonMomY.second);
0119       m_cInput -> SetBranchAddress("Parton4_MomZ", &m_partonMomZ.second, &m_brPartonMomZ.second);
0120       m_cInput -> SetBranchAddress("EvtSumParEne", &m_evtSumPar,         &m_brEvtSumPar);
0121       m_cInput -> SetBranchAddress("CstID",        &m_cstID,             &m_brCstID);
0122       m_cInput -> SetBranchAddress("CstEmbedID",   &m_cstEmbedID,        &m_brCstEmbedID);
0123     } else {
0124       m_cInput -> SetBranchAddress("EvtNumTrks",    &m_evtNumTrks, &m_brEvtNumTrks);
0125       m_cInput -> SetBranchAddress("EvtSumECalEne", &m_evtSumECal, &m_brEvtSumECal);
0126       m_cInput -> SetBranchAddress("EvtSumHCalEne", &m_evtSumHCal, &m_brEvtSumHCal);
0127       m_cInput -> SetBranchAddress("CstMatchID",    &m_cstMatchID, &m_brCstMatchID);
0128     }
0129
0130     // set generic branch addresses
0131     m_cInput -> SetBranchAddress("EvtVtxX",    &m_evtVtxX,    &m_brEvtVtxX);
0132     m_cInput -> SetBranchAddress("EvtVtxY",    &m_evtVtxY,    &m_brEvtVtxY);
0133     m_cInput -> SetBranchAddress("EvtVtxZ",    &m_evtVtxZ,    &m_brEvtVtxZ);
0134     m_cInput -> SetBranchAddress("EvtNumJets", &m_evtNumJets, &m_brEvtNumJets);
0135     m_cInput -> SetBranchAddress("JetNumCst",  &m_jetNumCst,  &m_brJetNumCst);
0136     m_cInput -> SetBranchAddress("JetID",      &m_jetID,      &m_brJetID);
0137     m_cInput -> SetBranchAddress("JetEnergy",  &m_jetEnergy,  &m_brJetEnergy);
0138     m_cInput -> SetBranchAddress("JetPt",      &m_jetPt,      &m_brJetPt);
0139     m_cInput -> SetBranchAddress("JetEta",     &m_jetEta,     &m_brJetEta);
0140     m_cInput -> SetBranchAddress("JetPhi",     &m_jetPhi,     &m_brJetPhi);
0141     m_cInput -> SetBranchAddress("JetArea",    &m_jetArea,    &m_brJetArea);
0142     m_cInput -> SetBranchAddress("CstZ",       &m_cstZ,       &m_brCstZ);
0143     m_cInput -> SetBranchAddress("CstDr",      &m_cstDr,      &m_brCstDr);
0144     m_cInput -> SetBranchAddress("CstEnergy",  &m_cstEnergy,  &m_brCstEnergy);
0145     m_cInput -> SetBranchAddress("CstJt",      &m_cstPt,      &m_brCstPt);
0146     m_cInput -> SetBranchAddress("CstEta",     &m_cstEta,     &m_brCstEta);
0147     m_cInput -> SetBranchAddress("CstPhi",     &m_cstPhi,     &m_brCstPhi);
0148     cout << "    Initialized input chain." << endl;
0149
0150     // exit routine
0151     return;
0152
0153   }  // end 'InitInput()'
0154
0155
0156
0157   void SCheckCstPairs::InitHists() {
0158
0159     // FIXME remove dependence on EECLongestSide
0160     vector<double> drBinEdges  = m_eecLongSide[0] -> bin_edges();
0161     size_t         nDrBinEdges = drBinEdges.size();
0162
0163     double drBinEdgeArray[nDrBinEdges];
0164     for (size_t iDrEdge = 0; iDrEdge < nDrBinEdges; iDrEdge++) {
0165       drBinEdgeArray[iDrEdge] = drBinEdges.at(iDrEdge);
0166     }
0167     hCstPtOneVsDr      = new TH2D("hCstPtOneVsDr",      "", m_nBinsDr, drBinEdgeArray, 200,  0.,   100.);
0168     hCstPtTwoVsDr      = new TH2D("hCstPtTwoVsDr",      "", m_nBinsDr, drBinEdgeArray, 200,  0.,   100.);
0169     hCstPtFracVsDr     = new TH2D("hCstPtFracVsDr",     "", m_nBinsDr, drBinEdgeArray, 500,  0.,   5.);
0170     hCstPhiOneVsDr     = new TH2D("hCstPhiOneVsDr",     "", m_nBinsDr, drBinEdgeArray, 360, -3.15, 3.15);;
0171     hCstPhiTwoVsDr     = new TH2D("hCstPhiTwoVsDr",     "", m_nBinsDr, drBinEdgeArray, 360, -3.15, 3.15);
0172     hCstEtaOneVsDr     = new TH2D("hCstEtaOneVsDr",     "", m_nBinsDr, drBinEdgeArray, 400, -2.,   2.);
0173     hCstEtaTwoVsDr     = new TH2D("hCstEtaTwoVsDr",     "", m_nBinsDr, drBinEdgeArray, 400, -2.,   2.);
0174     hDeltaPhiOneVsDr   = new TH2D("hDeltaPhiOneVsDr",   "", m_nBinsDr, drBinEdgeArray, 720, -6.30, 6.30);
0175     hDeltaPhiTwoVsDr   = new TH2D("hDeltaPhiTwoVsDr",   "", m_nBinsDr, drBinEdgeArray, 720, -6.30, 6.30);
0176     hDeltaEtaOneVsDr   = new TH2D("hDeltaEtaOneVsDr",   "", m_nBinsDr, drBinEdgeArray, 800, -4.,   4.);
0177     hDeltaEtaTwoVsDr   = new TH2D("hDeltaEtaTwoVsDr",   "", m_nBinsDr, drBinEdgeArray, 800, -4.,   4.);
0178     hJetPtFracOneVsDr  = new TH2D("hJetPtFracOneVsDr",  "", m_nBinsDr, drBinEdgeArray, 500, 0.,    5.);
0179     hJetPtFracTwoVsDr  = new TH2D("hJetPtFracTwoVsDr",  "", m_nBinsDr, drBinEdgeArray, 500, 0.,    5.);
0180     hCstPairWeightVsDr = new TH2D("hCstPairWeightVsDr", "", m_nBinsDr, drBinEdgeArray, 100,  0.,   1.);
0181     cout << "    Initialized histograms." << endl;
0182
0183     // exit routine
0184     return;
0185
0186   }  // end 'InitHists()'
0187
0188
0189
0190   void SCheckCstPairs::SaveOutput() {
0191
0192     m_outDir           -> cd();
0193     hCstPtOneVsDr      -> Write();
0194     hCstPtTwoVsDr      -> Write();
0195     hCstPtFracVsDr     -> Write();
0196     hCstPhiOneVsDr     -> Write();
0197     hCstPhiTwoVsDr     -> Write();
0198     hCstEtaOneVsDr     -> Write();
0199     hCstEtaTwoVsDr     -> Write();
0200     hDeltaPhiOneVsDr   -> Write();
0201     hDeltaPhiTwoVsDr   -> Write();
0202     hDeltaEtaOneVsDr   -> Write();
0203     hDeltaEtaTwoVsDr   -> Write();
0204     hJetPtFracOneVsDr  -> Write();
0205     hJetPtFracTwoVsDr  -> Write();
0206     hCstPairWeightVsDr -> Write();
0207     return;
0208
0209   }  // end 'SaveOutput()'
0210
0211
0212
0213
0214   void SCheckCstPairs::CloseInput() {
0215
0216     m_fInput -> cd();
0217     m_fInput -> Close();
0218     return;
0219
0220   }  // end 'CloseInput()'
0221
0222
0223
0224   void SCheckCstPairs::DoDoubleCstLoop() {
0225
0226     // announce start of event loop
0227     const uint64_t nEvts = m_cInput -> GetEntriesFast();
0228
0229     // event loop
0230     uint64_t nBytes = 0;
0231     for (uint64_t iEvt = 0; iEvt < nEvts; iEvt++) {
0232
0233       const uint64_t entry = Interfaces::LoadTree(m_cInput, iEvt, m_fCurrent);
0234       if (entry < 0) break;
0235
0236       const uint64_t bytes = Interfaces::GetEntry(m_cInput, iEvt);
0237       if (bytes < 0) {
0238         break;
0239       } else {
0240         nBytes += bytes;
0241       }
0242
0243       // announce progress
0244       const int64_t iProg = iEvt + 1;
0245       if (iProg == nEvents) {
0246         cout << "      Processing entry " << iEvt << "/" << nEvts << "..." << endl;
0247       } else {
0248         cout << "      Processing entry " << iEvt << "/" << nEvts << "...\r" << flush;
0249       }
0250
0251       // jet loop
0252       uint64_t nJets = (int) m_evtNumJets;
0253       for (uint64_t iJet = 0; iJet < nJets; iJet++) {
0254
0255         // clear vector for correlator
0256         m_jetCstVector.clear();
0257
0258         // get jet info
0259         const uint64_t nCsts  = m_jetNumCst -> at(iJet);
0260         const double   ptJet  = m_jetPt     -> at(iJet);
0261         const double   etaJet = m_jetEta    -> at(iJet);
0262
0263         // select jet pt bin & apply jet cuts
0264         const uint32_t iPtJetBin = GetJetPtBin(ptJet);
0265         const bool     isGoodJet = ApplyJetCuts(ptJet, etaJet);
0266         if (!isGoodJet) continue;
0267
0268         // constituent loop
0269         for (uint64_t iCst = 0; iCst < nCsts; iCst++) {
0270
0271           // get cst info
0272           const double drCst  = (m_cstDr  -> at(iJet)).at(iCst);
0273           const double etaCst = (m_cstEta -> at(iJet)).at(iCst);
0274           const double phiCst = (m_cstPhi -> at(iJet)).at(iCst);
0275           const double ptCst  = (m_cstPt  -> at(iJet)).at(iCst);
0276
0277           // for weird cst check
0278           for (uint64_t jCst = 0; jCst < nCsts; jCst++) {
0279
0280             // skip over the same cst
0281             if (jCst == iCst) continue;
0282
0283             // get cst info
0284             const double etaCstB = (m_cstEta -> at(iJet)).at(jCst);
0285             const double phiCstB = (m_cstPhi -> at(iJet)).at(jCst);
0286             const double ptCstB  = (m_cstPt  -> at(iJet)).at(jCst);
0287
0288             // calculate separation and pt-weight
0289             const double dhCstAB  = (etaCst - etaCstB);
0290             const double dfCstAB  = (phiCst - phiCstB);
0291             const double drCstAB  = sqrt((dhCstAB * dhCstAB) + (dfCstAB * dfCstAB));
0292             const double ptFrac   = ptCst / ptCstB;
0293             const double ztJetA   = ptCst / ptJet;
0294             const double ztJetB   = ptCstB / ptJet;
0295             const double ptWeight = (ptCst * ptCstB) / (ptJet * ptJet);
0296             hCstPtOneVsDr      -> Fill(drCstAB, ptCst);
0297             hCstPtTwoVsDr      -> Fill(drCstAB, ptCstB);
0298             hCstPtFracVsDr     -> Fill(drCstAB, ptFrac);
0299             hCstPhiOneVsDr     -> Fill(drCstAB, phiCst);
0300             hCstPhiTwoVsDr     -> Fill(drCstAB, phiCstB);
0301             hCstEtaOneVsDr     -> Fill(drCstAB, etaCst);
0302             hCstEtaTwoVsDr     -> Fill(drCstAB, etaCstB);
0303             hDeltaPhiOneVsDr   -> Fill(drCstAB, dfCstAB);
0304             hDeltaPhiTwoVsDr   -> Fill(drCstAB, dfCstAB);
0305             hDeltaEtaOneVsDr   -> Fill(drCstAB, dhCstAB);
0306             hDeltaEtaTwoVsDr   -> Fill(drCstAB, dhCstAB);
0307             hJetPtFracOneVsDr  -> Fill(drCstAB, ztJetA);
0308             hJetPtFracTwoVsDr  -> Fill(drCstAB, ztJetB);
0309             hCstPairWeightVsDr -> Fill(drCstAB, ptWeight);
0310           }  // end 2nd cst loop
0311         }  // end 1st cst loop
0312       }  // end jet loop
0313     }  // end event loop
0314     cout << "    Event loop finished!" << endl;
0315
0316     // exit routine
0317     return;
0318
0319   }  // end 'DoDoubleCstLoop()'
0320
0321
0322
0323   bool SCheckCstPairs::IsGoodJet() {
0324
0325     /* TODO fill in */
0326     return true;
0327
0328   }  // end 'IsGoodJet'
0329
0330
0331
0332   bool SCheckCstPairs::IsGoodCst() {
0333
0334     /* TODO fill in */
0335     return true;
0336
0337   }  // end 'IsGoodCst()'
0338
0339 }  // end SColdQcdCorrelatorAnalysis namespace
0340
0341 #endif
0342
0343 // end ------------------------------------------------------------------------