SDeltaPtCutStudy/src/SDeltaPtCutStudy.ana.h

0001 // ----------------------------------------------------------------------------
0002 // 'SDeltaPtCutStudy.ana.h'
0003 // Derek Anderson
0004 // 07.07.2023
0005 //
0006 // Reads in the 'ntp_track' Ntuple
0007 // generated by the SVtxEvaluator
0008 // class and studies how deltapt/pt
0009 // varies with quality cuts.
0010 // ----------------------------------------------------------------------------
0011
0012 #pragma once
0013
0014 using namespace std;
0015
0016
0017
0018 // analysis methods -----------------------------------------------------------
0019
0020 void SDeltaPtCutStudy::ApplyFlatDeltaPtCuts() {
0021
0022   // announce start of track loop
0023   cout << "      First loop over reco. tracks:" << endl;
0024
0025   // 1st track loop
0026   uint64_t nBytesTrk = 0;
0027   for (uint64_t iTrk = 0; iTrk < nTrks; iTrk++) {
0028
0029     // grab entry
0030     const uint64_t bytesTrk = ntTrack -> GetEntry(iTrk);
0031     if (bytesTrk < 0.) {
0032       cerr << "WARNING: something wrong with track #" << iTrk << "! Aborting loop!" << endl;
0033       break;
0034     }
0035     nBytesTrk += bytesTrk;
0036
0037     // announce progress
0038     const uint64_t iProgTrk = iTrk + 1;
0039     if (iProgTrk == nTrks) {
0040       cout << "        Processing track " << iProgTrk << "/" << nTrks << "..." << endl;
0041     } else {
0042       cout << "        Processing track " << iProgTrk << "/" << nTrks << "...\r" << flush;
0043     }
0044
0045     // do calculations
0046     const double ptFrac  = trk_pt / trk_gpt;
0047     const double ptDelta = trk_deltapt / trk_pt;
0048
0049     // apply trk cuts
0050     const bool isInZVtxCut = (abs(trk_vz) <  vzTrkMax);
0051     const bool isInInttCut = (trk_nintt   >= nInttTrkMin);
0052     const bool isInMVtxCut = (trk_nlmaps  >  nMVtxTrkMin);
0053     const bool isInTpcCut  = (trk_ntpc    >  nTpcTrkMin);
0054     const bool isInPtCut   = (trk_pt      >  ptTrkMin);
0055     const bool isInQualCut = (trk_quality <  qualTrkMax);
0056     const bool isGoodTrk   = (isInZVtxCut && isInInttCut && isInMVtxCut && isInTpcCut && isInPtCut && isInQualCut);
0057     if (!isGoodTrk) continue;
0058
0059     // fill histograms
0060     hPtDelta        -> Fill(ptDelta);
0061     hPtTrack        -> Fill(trk_pt);
0062     hPtFrac         -> Fill(ptFrac);
0063     hPtTrkTru       -> Fill(trk_gpt);
0064     hPtDeltaVsFrac  -> Fill(ptFrac,  ptDelta);
0065     hPtDeltaVsTrue  -> Fill(trk_gpt, ptDelta);
0066     hPtDeltaVsTrack -> Fill(trk_pt,  ptDelta);
0067     hPtTrueVsTrack  -> Fill(trk_pt,  trk_gpt);
0068
0069     // apply delta-pt cuts
0070     const bool isNormalTrk = ((ptFrac > normRange[0]) && (ptFrac < normRange[1]));
0071     for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
0072       const bool isInDeltaPtCut = (ptDelta < ptDeltaMax[iCut]);
0073       if (isInDeltaPtCut) {
0074
0075         // fill histograms
0076         hPtDeltaCut[iCut]        -> Fill(ptDelta);
0077         hPtTrackCut[iCut]        -> Fill(trk_pt);
0078         hPtFracCut[iCut]         -> Fill(ptFrac);
0079         hPtTrkTruCut[iCut]       -> Fill(trk_gpt);
0080         hPtDeltaVsFracCut[iCut]  -> Fill(ptFrac,  ptDelta);
0081         hPtDeltaVsTrueCut[iCut]  -> Fill(trk_gpt, ptDelta);
0082         hPtDeltaVsTrackCut[iCut] -> Fill(trk_pt,  ptDelta);
0083         hPtTrueVsTrackCut[iCut]  -> Fill(trk_pt,  trk_gpt);
0084
0085         // increment counters
0086         if (isNormalTrk) {
0087           ++nNormCut[iCut];
0088         } else {
0089           ++nWeirdCut[iCut];
0090         }
0091       }
0092     }  // end delta-pt cut
0093   }  // end 1st track loop
0094
0095   cout << "      First loop over reco. tracks finished!" << endl;
0096   return;
0097
0098 }  // end 'ApplyFlatDeltaPtCuts()'
0099
0100
0101
0102 void SDeltaPtCutStudy::ApplyPtDependentDeltaPtCuts() {
0103
0104   // announce start of track loop
0105   cout << "      Second loop over reco. tracks:" << endl;
0106
0107   // 2nd track loop
0108   uint64_t nBytesTrk = 0;
0109   for (uint64_t iTrk = 0; iTrk < nTrks; iTrk++) {
0110
0111     // grab entry
0112     const uint64_t bytesTrk = ntTrack -> GetEntry(iTrk);
0113     if (bytesTrk < 0.) {
0114       cerr << "WARNING: something wrong with track #" << iTrk << "! Aborting loop!" << endl;
0115       break;
0116     }
0117     nBytesTrk += bytesTrk;
0118
0119     // announce progress
0120     const uint64_t iProgTrk = iTrk + 1;
0121     if (iProgTrk == nTrks) {
0122       cout << "        Processing track " << iProgTrk << "/" << nTrks << "..." << endl;
0123     } else {
0124       cout << "        Processing track " << iProgTrk << "/" << nTrks << "...\r" << flush;
0125     }
0126
0127     // do calculations
0128     const double ptFrac  = trk_pt / trk_gpt;
0129     const double ptDelta = trk_deltapt / trk_pt;
0130
0131     // apply trk cuts
0132     const bool isInZVtxCut = (abs(trk_vz) <  vzTrkMax);
0133     const bool isInInttCut = (trk_nintt   >= nInttTrkMin);
0134     const bool isInMVtxCut = (trk_nlmaps  >  nMVtxTrkMin);
0135     const bool isInTpcCut  = (trk_ntpc    >  nTpcTrkMin);
0136     const bool isInPtCut   = (trk_pt      >  ptTrkMin);
0137     const bool isInQualCut = (trk_quality <  qualTrkMax);
0138     const bool isGoodTrk   = (isInZVtxCut && isInInttCut && isInMVtxCut && isInTpcCut && isInPtCut && isInQualCut);
0139     if (!isGoodTrk) continue;
0140
0141     // apply delta-pt cuts
0142     const bool isNormalTrk = ((ptFrac > normRange[0]) && (ptFrac < normRange[1]));
0143     for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
0144
0145       // get bounds
0146       const float ptDeltaMin = fMuLoProj[iSig] -> Eval(trk_pt);
0147       const float ptDeltaMax = fMuHiProj[iSig] -> Eval(trk_pt);
0148
0149       const bool isInDeltaPtSigma = ((ptDelta >= ptDeltaMin) && (ptDelta <= ptDeltaMax));
0150       if (isInDeltaPtSigma) {
0151
0152         // fill histograms
0153         hPtDeltaSig[iSig]        -> Fill(ptDelta);
0154         hPtTrackSig[iSig]        -> Fill(trk_pt);
0155         hPtFracSig[iSig]         -> Fill(ptFrac);
0156         hPtTrkTruSig[iSig]       -> Fill(trk_gpt);
0157         hPtDeltaVsFracSig[iSig]  -> Fill(ptFrac,  ptDelta);
0158         hPtDeltaVsTrueSig[iSig]  -> Fill(trk_gpt, ptDelta);
0159         hPtDeltaVsTrackSig[iSig] -> Fill(trk_pt,  ptDelta);
0160         hPtTrueVsTrackSig[iSig]  -> Fill(trk_pt,  trk_gpt);
0161
0162         // increment counters
0163         if (isNormalTrk) {
0164           ++nNormSig[iSig];
0165         } else {
0166           ++nWeirdSig[iSig];
0167         }
0168       }
0169     }  // end delta-pt cut
0170   }  // end 1st track loop
0171
0172   cout << "      Second loop over reco. tracks finished!" << endl;
0173   return;
0174
0175 }  // end 'ApplyPtDependentDeltaptCuts()'
0176
0177
0178
0179 void SDeltaPtCutStudy::FillTruthHistograms() {
0180
0181   // announce start of truth loop
0182   cout << "      Loop over particles:" << endl;
0183
0184   // truth loop
0185   uint64_t nBytesTru = 0;
0186   for (uint64_t iTru = 0; iTru < nTrus; iTru++) {
0187
0188     // grab entry
0189     const uint64_t bytesTru = ntTruth -> GetEntry(iTru);
0190     if (bytesTru < 0.) {
0191       cerr << "WARNING: something wrong with particle #" << iTru << "! Aborting loop!" << endl;
0192       break;
0193     }
0194     nBytesTru += bytesTru;
0195
0196     // announce progress
0197     const uint64_t iProgTru = iTru + 1;
0198     if (iProgTru == nTrus) {
0199       cout << "        Processing particle " << iProgTru << "/" << nTrus << "..." << endl;
0200     } else {
0201       cout << "        Processing particle " << iProgTru << "/" << nTrus << "...\r" << flush;
0202     }
0203
0204     // fill truth histogram
0205     const bool isPrimary = (tru_gprimary == 1);
0206     if (isPrimary) {
0207       hPtTruth -> Fill(tru_gpt);
0208     }
0209   }  // end track loop
0210   cout << "      Loop over particles finished!" << endl;
0211
0212 }  // end 'FillTruthHistograms()'
0213
0214
0215
0216 void SDeltaPtCutStudy::CreateSigmaGraphs() {
0217
0218   // for graph names
0219   const TString sMuHiBase = "MeanPlusSigma";
0220   const TString sMuLoBase = "MeanMinusSigma";
0221   const TString sSigBase  = "ProjectionSigma";
0222   const TString sMuBase   = "ProjectionMean";
0223
0224   // projection fit names
0225   vector<TString> sFitProj(nProj);
0226   for (size_t iProj = 0; iProj < nProj; iProj++) {
0227     sFitProj[iProj] = "f";
0228     sFitProj[iProj].Append(sPtProjBase.Data());
0229     sFitProj[iProj].Append(sProjSuffix[iProj].Data());
0230   }
0231
0232   // project slices of delta-pt and get sigmas
0233   const uint32_t fWidFit = 2;
0234   const uint32_t fLinFit = 1;
0235   for (size_t iProj = 0; iProj < nProj; iProj++) {
0236
0237     // do projection
0238     const uint32_t iBinProj = hPtDeltaVsTrack -> GetXaxis() -> FindBin(ptProj[iProj]);
0239     hPtDeltaProj[iProj]   = hPtDeltaVsTrack -> ProjectionY(sPtProj[iProj], iBinProj, iBinProj, "");
0240
0241     // get initial values for fit
0242     const float ampGuess = hPtDeltaProj[iProj] -> GetMaximum();
0243     const float muGuess  = hPtDeltaProj[iProj] -> GetMean();
0244     const float sigGuess = hPtDeltaProj[iProj] -> GetRMS();
0245
0246     // fit with gaussian
0247     fPtDeltaProj[iProj] = new TF1(sFitProj[iProj].Data(), "gaus", deltaFitRange[0], deltaFitRange[1]);
0248     fPtDeltaProj[iProj] -> SetLineColor(fColFit[iProj]);
0249     fPtDeltaProj[iProj] -> SetLineStyle(fLinFit);
0250     fPtDeltaProj[iProj] -> SetLineWidth(fWidFit);
0251     fPtDeltaProj[iProj] -> SetParameter(0, ampGuess);
0252     fPtDeltaProj[iProj] -> SetParameter(1, muGuess);
0253     fPtDeltaProj[iProj] -> SetParameter(2, sigGuess);
0254     hPtDeltaProj[iProj] -> Fit(sFitProj[iProj].Data(), "R");
0255
0256     // add values to arrays
0257     muProj[iProj]  = fPtDeltaProj[iProj] -> GetParameter(1);
0258     sigProj[iProj] = fPtDeltaProj[iProj] -> GetParameter(2);
0259     for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
0260       muHiProj[iSig][iProj] = muProj[iProj] + (ptDeltaSig[iSig] * sigProj[iProj]);
0261       muLoProj[iSig][iProj] = muProj[iProj] - (ptDeltaSig[iSig] * sigProj[iProj]);
0262     }
0263   }  // end projection loop
0264   cout << "      Obtained delta-pt projections, fits, and sigmas." << endl;
0265
0266   // sigma graph names
0267   TString sMuProj("gr");
0268   TString sSigProj("gr");
0269   sMuProj.Append(sMuBase.Data());
0270   sSigProj.Append(sSigBase.Data());
0271
0272   vector<TString> sGrMuHiProj(nSigCuts);
0273   vector<TString> sGrMuLoProj(nSigCuts);
0274   vector<TString> sFnMuHiProj(nSigCuts);
0275   vector<TString> sFnMuLoProj(nSigCuts);
0276   for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
0277     sGrMuHiProj[iSig] = "gr";
0278     sGrMuLoProj[iSig] = "gr";
0279     sFnMuHiProj[iSig] = "f";
0280     sFnMuLoProj[iSig] = "f";
0281     sGrMuHiProj[iSig].Append(sMuHiBase.Data());
0282     sGrMuLoProj[iSig].Append(sMuLoBase.Data());
0283     sFnMuHiProj[iSig].Append(sMuHiBase.Data());
0284     sFnMuLoProj[iSig].Append(sMuLoBase.Data());
0285     sGrMuHiProj[iSig].Append(sSigSuffix[iSig].Data());
0286     sGrMuLoProj[iSig].Append(sSigSuffix[iSig].Data());
0287     sFnMuHiProj[iSig].Append(sSigSuffix[iSig].Data());
0288     sFnMuLoProj[iSig].Append(sSigSuffix[iSig].Data());
0289   }
0290
0291   // turn std::vectors into TVectors
0292   TVectorD tvecPtProj(ptProj.size(), ptProj.data());
0293   TVectorD tvecMuProj(muProj.size(), muProj.data());
0294   TVectorD tvecSigProj(sigProj.size(), sigProj.data());
0295
0296   vector<TVectorD> tvecMuHiProj;
0297   vector<TVectorD> tvecMuLoProj;
0298   for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
0299     tvecMuHiProj.push_back(TVectorD(muHiProj[iSig].size(), muHiProj[iSig].data()));
0300     tvecMuLoProj.push_back(TVectorD(muLoProj[iSig].size(), muLoProj[iSig].data()));
0301   }
0302
0303   // construct sigma graphs
0304   grMuProj  = new TGraph(tvecPtProj, tvecMuProj);
0305   grSigProj = new TGraph(tvecPtProj, tvecSigProj);
0306   grMuProj  -> SetName(sMuProj);
0307   grSigProj -> SetName(sSigProj);
0308
0309   // fit sigma graphs
0310   for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
0311
0312     // create graphs
0313     grMuHiProj[iSig] = new TGraph(tvecPtProj, tvecMuHiProj[iSig]);
0314     grMuLoProj[iSig] = new TGraph(tvecPtProj, tvecMuLoProj[iSig]);
0315     grMuHiProj[iSig] -> SetName(sGrMuHiProj[iSig].Data());
0316     grMuLoProj[iSig] -> SetName(sGrMuLoProj[iSig].Data());
0317
0318     // create fit functions
0319     fMuHiProj[iSig] = new TF1(sFnMuHiProj[iSig].Data(), "pol2", rPtRange[0], rPtRange[1]);
0320     fMuLoProj[iSig] = new TF1(sFnMuLoProj[iSig].Data(), "pol2", rPtRange[0], rPtRange[1]);
0321     fMuHiProj[iSig] -> SetLineColor(fColSigFit[iSig]);
0322     fMuLoProj[iSig] -> SetLineColor(fColSigFit[iSig]);
0323     fMuHiProj[iSig] -> SetLineStyle(fLinFit);
0324     fMuLoProj[iSig] -> SetLineStyle(fLinFit);
0325     fMuHiProj[iSig] -> SetLineWidth(fWidFit);
0326     fMuLoProj[iSig] -> SetLineWidth(fWidFit);
0327     fMuHiProj[iSig] -> SetParameter(0, sigHiGuess[0]);
0328     fMuLoProj[iSig] -> SetParameter(0, sigLoGuess[0]);
0329     fMuHiProj[iSig] -> SetParameter(1, sigHiGuess[1]);
0330     fMuLoProj[iSig] -> SetParameter(1, sigLoGuess[1]);
0331     fMuHiProj[iSig] -> SetParameter(2, sigHiGuess[2]);
0332     fMuLoProj[iSig] -> SetParameter(2, sigLoGuess[2]);
0333
0334     // do fitting
0335     grMuHiProj[iSig] -> Fit(sFnMuHiProj[iSig].Data(), "", "", ptFitRange[0], ptFitRange[1]);
0336     grMuLoProj[iSig] -> Fit(sFnMuLoProj[iSig].Data(), "", "", ptFitRange[0], ptFitRange[1]);
0337   }
0338
0339   cout << "      Created and fit sigma graphs."  << endl;
0340   return;
0341
0342 }  // end 'CreateSigmaGraphs()'
0343
0344
0345
0346 void SDeltaPtCutStudy::CalculateRejectionFactors() {
0347
0348   // for graph names
0349   const TString sRejCutBase = "Reject_flatDPtCut";
0350   const TString sRejSigBase = "Reject_sigmaCut";
0351
0352   // calculate flat delta-pt rejection factors
0353   for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
0354     rejCut[iCut] = (double) nNormCut[iCut] / (double) nWeirdCut[iCut];
0355   }
0356   cout << "      Calculated flat delta-pt rejection factors." << endl;
0357
0358   // calculate pt-dependent delta-pt rejection factors
0359   for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
0360     rejSig[iSig] = (double) nNormSig[iSig] / (double) nWeirdSig[iSig];
0361   }
0362   cout << "      Calculated pt-depdendent delta-pt rejection factors\n"
0363        << "      Rejection factors:\n"
0364        << "        Flat delta-pt cuts"
0365        << endl;
0366
0367   // announce flat delta-pt rejection factors
0368   for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
0369     cout << "          n(Norm, Weird) = (" << nNormCut[iCut] << ", " << nWeirdCut[iCut] << "), rejection = " << rejCut[iCut] << endl;
0370   }
0371
0372   // announce pt-dependent delta-pt rejection factors
0373   cout << "        Pt-dependent delta-pt cuts" << endl;
0374   for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
0375     cout << "          n(Norm, Weird) = (" << nNormSig[iSig] << ", " << nWeirdSig[iSig] << "), rejection = " << rejSig[iSig] << endl;
0376   }
0377
0378   // graph names
0379   TString sRejCut("gr");
0380   TString sRejSig("gr");
0381   sRejCut.Append(sRejCutBase.Data());
0382   sRejSig.Append(sRejSigBase.Data());
0383
0384   // convert vectors to TVectors
0385   TVectorD tvecPtDeltaMax(ptDeltaMax.size(), ptDeltaMax.data());
0386   TVectorD tvecPtDeltaSig(ptDeltaSig.size(), ptDeltaSig.data());
0387   TVectorD tvecRejCut(rejCut.size(), rejCut.data());
0388   TVectorD tvecRejSig(rejSig.size(), rejSig.data());
0389
0390   // make rejection graphs
0391   grRejCut = new TGraph(tvecPtDeltaMax, tvecRejCut);
0392   grRejSig = new TGraph(tvecPtDeltaSig, tvecRejSig);
0393   grRejCut -> SetName(sRejCut.Data());
0394   grRejSig -> SetName(sRejSig.Data());
0395
0396   cout << "      Made rejection factor graph." << endl;
0397   return;
0398
0399 }  // end 'CalculateRejectionFactors()'
0400
0401
0402
0403 void SDeltaPtCutStudy::CalculateEfficiencies() {
0404
0405   // for histogram names
0406   const TString sEffBase = "Efficiency";
0407
0408   // rebin histograms if needed
0409   if (doEffRebin) {
0410     hPtTruth  -> Rebin(nEffRebin);
0411     hPtTrkTru -> Rebin(nEffRebin);
0412     for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
0413       hPtTrkTruCut[iCut] -> Rebin(nEffRebin);
0414     }
0415     for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
0416       hPtTrkTruSig[iSig] -> Rebin(nEffRebin);
0417     }
0418     cout << "      Rebinned efficiency histograms." << endl;
0419   }
0420
0421   TString sEff("h");
0422   sEff.Append(sEffBase.Data());
0423
0424   // create flat delta-pt cut efficiency names
0425   vector<TString> sEffCut(nDPtCuts);
0426   for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
0427     sEffCut[iCut] = "h";
0428     sEffCut[iCut].Append(sEffBase.Data());
0429     sEffCut[iCut].Append(sDPtSuffix[iCut].Data());
0430   }
0431
0432   // create pt-dependent delta-pt cut efficiency names
0433   vector<TString> sEffSig(nSigCuts);
0434   for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
0435     sEffSig[iSig] = "h";
0436     sEffSig[iSig].Append(sEffBase.Data());
0437     sEffSig[iSig].Append(sSigSuffix[iSig].Data());
0438   }
0439
0440   hEff = (TH1D*) hPtTruth -> Clone();
0441   hEff -> SetName(sEff.Data());
0442   hEff -> Reset("ICES");
0443   hEff -> Divide(hPtTrkTru, hPtTruth, 1., 1.);
0444
0445   // calculate flat delta-pt cut efficiencies
0446   for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
0447     hEffCut[iCut] = (TH1D*) hPtTruth -> Clone();
0448     hEffCut[iCut] -> SetName(sEffCut[iCut].Data());
0449     hEffCut[iCut] -> Reset("ICES");
0450     hEffCut[iCut] -> Divide(hPtTrkTruCut[iCut], hPtTruth, 1., 1.);
0451   }
0452
0453   // calculate pt-dependent delta-pt cut efficiencies
0454   for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
0455     hEffSig[iSig] = (TH1D*) hPtTruth -> Clone();
0456     hEffSig[iSig] -> SetName(sEffSig[iSig].Data());
0457     hEffSig[iSig] -> Reset("ICES");
0458     hEffSig[iSig] -> Divide(hPtTrkTruSig[iSig], hPtTruth, 1., 1.);
0459   }
0460
0461   cout << "      Calculated efficiencies." << endl;
0462   return;
0463
0464 }  // end 'CalculateEfficiencies()'
0465
0466 // end ------------------------------------------------------------------------