Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-11-19 23:20:34

0001 #include "performance.h"
0002 
0003 enum { pT5 = 7, pT3 = 5, T5 = 4, pLS = 8 };
0004 
0005 //__________________________________________________________________________________________________________________________________________________________________________
0006 int main(int argc, char** argv) {
0007   // Parse arguments
0008   parseArguments(argc, argv);
0009 
0010   // Initialize input and output root files
0011   initializeInputsAndOutputs();
0012 
0013   // Set of pdgids
0014   std::vector<int> pdgids = {0, 11, 211, 13, 321};
0015 
0016   // Set of charges
0017   std::vector<int> charges = {0, 1, -1};
0018 
0019   // Set of extra selections for efficiency plots
0020   std::vector<TString> selnames = {
0021       "base",    // default baseline that is more inline with MTV
0022       "loweta",  // When the eta cut is restricted to 2.4
0023       "xtr",     // When the eta cut is restricted to transition regions
0024       "vtr"      // When the eta cut is vetoing transition regions
0025   };
0026   std::vector<std::function<bool(unsigned int)>> sels = {
0027       [&](unsigned int isim) { return 1.; },
0028       [&](unsigned int isim) { return abs(lstEff.sim_eta().at(isim)) < 2.4; },
0029       [&](unsigned int isim) { return abs(lstEff.sim_eta().at(isim)) > 1.1 and abs(lstEff.sim_eta().at(isim)) < 1.7; },
0030       [&](unsigned int isim) {
0031         return (abs(lstEff.sim_eta().at(isim)) < 1.1 or abs(lstEff.sim_eta().at(isim)) > 1.7) and
0032                abs(lstEff.sim_eta().at(isim)) < 2.4;
0033       }};
0034   pdgids.insert(pdgids.end(), ana.pdgids.begin(), ana.pdgids.end());
0035 
0036   std::vector<SimTrackSetDefinition> list_effSetDef;
0037 
0038   // creating a set of efficiency plots for each pdgids being considered
0039   for (auto& pdgid : pdgids) {
0040     for (auto& charge : charges) {
0041       for (unsigned int isel = 0; isel < sels.size(); ++isel) {
0042         list_effSetDef.push_back(
0043             SimTrackSetDefinition(/* name  */
0044                                   TString("TC_") + selnames[isel],
0045                                   /* pdgid */ pdgid,
0046                                   /* q     */ charge,
0047                                   /* pass  */ [&](unsigned int isim) { return lstEff.sim_TC_matched().at(isim) > 0; },
0048                                   /* sel   */ sels[isel]));
0049         list_effSetDef.push_back(SimTrackSetDefinition(
0050             /* name  */
0051             TString("pT5_") + selnames[isel],
0052             /* pdgid */ pdgid,
0053             /* q     */ charge,
0054             /* pass  */ [&](unsigned int isim) { return lstEff.sim_TC_matched_mask().at(isim) & (1 << pT5); },
0055             /* sel   */ sels[isel]));
0056         list_effSetDef.push_back(SimTrackSetDefinition(
0057             /* name  */
0058             TString("pT3_") + selnames[isel],
0059             /* pdgid */ pdgid,
0060             /* q     */ charge,
0061             /* pass  */ [&](unsigned int isim) { return lstEff.sim_TC_matched_mask().at(isim) & (1 << pT3); },
0062             /* sel   */ sels[isel]));
0063         list_effSetDef.push_back(SimTrackSetDefinition(
0064             /* name  */
0065             TString("T5_") + selnames[isel],
0066             /* pdgid */ pdgid,
0067             /* q     */ charge,
0068             /* pass  */ [&](unsigned int isim) { return lstEff.sim_TC_matched_mask().at(isim) & (1 << T5); },
0069             /* sel   */ sels[isel]));
0070         list_effSetDef.push_back(SimTrackSetDefinition(
0071             /* name  */
0072             TString("pLS_") + selnames[isel],
0073             /* pdgid */ pdgid,
0074             /* q     */ charge,
0075             /* pass  */ [&](unsigned int isim) { return lstEff.sim_TC_matched_mask().at(isim) & (1 << pLS); },
0076             /* sel   */ sels[isel]));
0077 
0078         if (ana.do_lower_level) {
0079           //lower objects - name will have pT5_lower_, T5_lower_, pT3_lower_
0080           list_effSetDef.push_back(SimTrackSetDefinition(
0081               /* name  */
0082               TString("pT5_lower_") + selnames[isel],
0083               /* pdgid */ pdgid,
0084               /* q     */ charge,
0085               /* pass  */ [&](unsigned int isim) { return lstEff.sim_pT5_matched().at(isim) > 0; },
0086               /* sel   */ sels[isel]));
0087           list_effSetDef.push_back(
0088               SimTrackSetDefinition(/* name  */
0089                                     TString("T5_lower_") + selnames[isel],
0090                                     /* pdgid */ pdgid,
0091                                     /* q     */ charge,
0092                                     /* pass  */ [&](unsigned int isim) { return lstEff.sim_T5_matched().at(isim) > 0; },
0093                                     /* sel   */ sels[isel]));
0094           list_effSetDef.push_back(SimTrackSetDefinition(
0095               /* name  */
0096               TString("pT3_lower_") + selnames[isel],
0097               /* pdgid */ pdgid,
0098               /* q     */ charge,
0099               /* pass  */ [&](unsigned int isim) { return lstEff.sim_pT3_matched().at(isim) > 0; },
0100               /* sel   */ sels[isel]));
0101         }
0102       }
0103     }
0104   }
0105 
0106   bookEfficiencySets(list_effSetDef);
0107 
0108   // creating a set of fake rate plots
0109   std::vector<RecoTrackSetDefinition> list_FRSetDef;
0110   list_FRSetDef.push_back(
0111       RecoTrackSetDefinition(/* name  */
0112                              "TC",
0113                              /* pass  */ [&](unsigned int itc) { return lstEff.tc_isFake().at(itc) > 0; },
0114                              /* sel   */ [&](unsigned int itc) { return 1; },
0115                              /* pt    */ tas::tc_pt,
0116                              /* eta   */ tas::tc_eta,
0117                              /* phi   */ tas::tc_phi,
0118                              /* type  */ tas::tc_type));
0119   list_FRSetDef.push_back(
0120       RecoTrackSetDefinition(/* name  */
0121                              "pT5",
0122                              /* pass  */ [&](unsigned int itc) { return lstEff.tc_isFake().at(itc) > 0; },
0123                              /* sel   */ [&](unsigned int itc) { return lstEff.tc_type().at(itc) == pT5; },
0124                              /* pt    */ tas::tc_pt,
0125                              /* eta   */ tas::tc_eta,
0126                              /* phi   */ tas::tc_phi,
0127                              /* type  */ tas::tc_type));
0128   list_FRSetDef.push_back(
0129       RecoTrackSetDefinition(/* name  */
0130                              "pT3",
0131                              /* pass  */ [&](unsigned int itc) { return lstEff.tc_isFake().at(itc) > 0; },
0132                              /* sel   */ [&](unsigned int itc) { return lstEff.tc_type().at(itc) == pT3; },
0133                              /* pt    */ tas::tc_pt,
0134                              /* eta   */ tas::tc_eta,
0135                              /* phi   */ tas::tc_phi,
0136                              /* type  */ tas::tc_type));
0137   list_FRSetDef.push_back(
0138       RecoTrackSetDefinition(/* name  */
0139                              "T5",
0140                              /* pass  */ [&](unsigned int itc) { return lstEff.tc_isFake().at(itc) > 0; },
0141                              /* sel   */ [&](unsigned int itc) { return lstEff.tc_type().at(itc) == T5; },
0142                              /* pt    */ tas::tc_pt,
0143                              /* eta   */ tas::tc_eta,
0144                              /* phi   */ tas::tc_phi,
0145                              /* type  */ tas::tc_type));
0146   list_FRSetDef.push_back(
0147       RecoTrackSetDefinition(/* name  */
0148                              "pLS",
0149                              /* pass  */ [&](unsigned int itc) { return lstEff.tc_isFake().at(itc) > 0; },
0150                              /* sel   */ [&](unsigned int itc) { return lstEff.tc_type().at(itc) == pLS; },
0151                              /* pt    */ tas::tc_pt,
0152                              /* eta   */ tas::tc_eta,
0153                              /* phi   */ tas::tc_phi,
0154                              /* type  */ tas::tc_type));
0155 
0156   if (ana.do_lower_level) {
0157     list_FRSetDef.push_back(RecoTrackSetDefinition(
0158         /* name  */
0159         "pT5_lower",
0160         /* pass  */ [&](unsigned int ipT5) { return lstEff.pT5_isFake().at(ipT5) > 0; },
0161         /* sel   */ [&](unsigned int ipT5) { return 1; },
0162         /* pt    */ tas::pT5_pt,
0163         /* eta   */ tas::pT5_eta,
0164         /* phi   */ tas::pT5_phi,
0165         /* type  */ [&]() { return static_cast<const std::vector<int>>(std::vector<int>(tas::pT5_pt().size(), 1)); }));
0166     list_FRSetDef.push_back(RecoTrackSetDefinition(
0167         /* name  */
0168         "T5_lower",
0169         /* pass  */ [&](unsigned int it5) { return lstEff.t5_isFake().at(it5) > 0; },
0170         /* sel   */ [&](unsigned int it5) { return 1; },
0171         /* pt    */ tas::t5_pt,
0172         /* eta   */ tas::t5_eta,
0173         /* phi   */ tas::t5_phi,
0174         /* type  */ [&]() { return static_cast<const std::vector<int>>(std::vector<int>(tas::t5_pt().size(), 1)); }));
0175     list_FRSetDef.push_back(RecoTrackSetDefinition(
0176         /* name  */
0177         "pT3_lower",
0178         /* pass  */ [&](unsigned int ipT3) { return lstEff.pT3_isFake().at(ipT3) > 0; },
0179         /* sel   */ [&](unsigned int ipT3) { return 1; },
0180         /* pt    */ tas::pT3_pt,
0181         /* eta   */ tas::pT3_eta,
0182         /* phi   */ tas::pT3_phi,
0183         /* type  */ [&]() { return static_cast<const std::vector<int>>(std::vector<int>(tas::pT3_pt().size(), 1)); }));
0184   }
0185 
0186   bookFakeRateSets(list_FRSetDef);
0187 
0188   // creating a set of duplicate rate plots
0189   std::vector<RecoTrackSetDefinition> list_DRSetDef;
0190   list_DRSetDef.push_back(
0191       RecoTrackSetDefinition(/* name  */
0192                              "TC",
0193                              /* pass  */ [&](unsigned int itc) { return lstEff.tc_isDuplicate().at(itc) > 0; },
0194                              /* sel   */ [&](unsigned int itc) { return 1; },
0195                              /* pt    */ tas::tc_pt,
0196                              /* eta   */ tas::tc_eta,
0197                              /* phi   */ tas::tc_phi,
0198                              /* type  */ tas::tc_type));
0199   list_DRSetDef.push_back(
0200       RecoTrackSetDefinition(/* name  */
0201                              "pT5",
0202                              /* pass  */ [&](unsigned int itc) { return lstEff.tc_isDuplicate().at(itc) > 0; },
0203                              /* sel   */ [&](unsigned int itc) { return lstEff.tc_type().at(itc) == pT5; },
0204                              /* pt    */ tas::tc_pt,
0205                              /* eta   */ tas::tc_eta,
0206                              /* phi   */ tas::tc_phi,
0207                              /* type  */ tas::tc_type));
0208   list_DRSetDef.push_back(
0209       RecoTrackSetDefinition(/* name  */
0210                              "pT3",
0211                              /* pass  */ [&](unsigned int itc) { return lstEff.tc_isDuplicate().at(itc) > 0; },
0212                              /* sel   */ [&](unsigned int itc) { return lstEff.tc_type().at(itc) == pT3; },
0213                              /* pt    */ tas::tc_pt,
0214                              /* eta   */ tas::tc_eta,
0215                              /* phi   */ tas::tc_phi,
0216                              /* type  */ tas::tc_type));
0217   list_DRSetDef.push_back(
0218       RecoTrackSetDefinition(/* name  */
0219                              "T5",
0220                              /* pass  */ [&](unsigned int itc) { return lstEff.tc_isDuplicate().at(itc) > 0; },
0221                              /* sel   */ [&](unsigned int itc) { return lstEff.tc_type().at(itc) == T5; },
0222                              /* pt    */ tas::tc_pt,
0223                              /* eta   */ tas::tc_eta,
0224                              /* phi   */ tas::tc_phi,
0225                              /* type  */ tas::tc_type));
0226   list_DRSetDef.push_back(
0227       RecoTrackSetDefinition(/* name  */
0228                              "pLS",
0229                              /* pass  */ [&](unsigned int itc) { return lstEff.tc_isDuplicate().at(itc) > 0; },
0230                              /* sel   */ [&](unsigned int itc) { return lstEff.tc_type().at(itc) == pLS; },
0231                              /* pt    */ tas::tc_pt,
0232                              /* eta   */ tas::tc_eta,
0233                              /* phi   */ tas::tc_phi,
0234                              /* type  */ tas::tc_type));
0235 
0236   if (ana.do_lower_level) {
0237     list_DRSetDef.push_back(RecoTrackSetDefinition(
0238         /* name  */
0239         "pT5_lower",
0240         /* pass  */ [&](unsigned int ipT5) { return lstEff.pT5_isDuplicate().at(ipT5) > 0; },
0241         /* sel   */ [&](unsigned int ipT5) { return 1; },
0242         /* pt    */ tas::pT5_pt,
0243         /* eta   */ tas::pT5_eta,
0244         /* phi   */ tas::pT5_phi,
0245         /* type  */ [&]() { return static_cast<const std::vector<int>>(std::vector<int>(tas::pT5_pt().size(), 1)); }));
0246     list_DRSetDef.push_back(RecoTrackSetDefinition(
0247         /* name  */
0248         "T5_lower",
0249         /* pass  */ [&](unsigned int it5) { return lstEff.t5_isDuplicate().at(it5) > 0; },
0250         /* sel   */ [&](unsigned int it5) { return 1; },
0251         /* pt    */ tas::t5_pt,
0252         /* eta   */ tas::t5_eta,
0253         /* phi   */ tas::t5_phi,
0254         /* type  */ [&]() { return static_cast<const std::vector<int>>(std::vector<int>(tas::t5_pt().size(), 1)); }));
0255     list_DRSetDef.push_back(RecoTrackSetDefinition(
0256         /* name  */
0257         "pT3_lower",
0258         /* pass  */ [&](unsigned int ipT3) { return lstEff.pT3_isDuplicate().at(ipT3) > 0; },
0259         /* sel   */ [&](unsigned int ipT3) { return 1; },
0260         /* pt    */ tas::pT3_pt,
0261         /* eta   */ tas::pT3_eta,
0262         /* phi   */ tas::pT3_phi,
0263         /* type  */ [&]() { return static_cast<const std::vector<int>>(std::vector<int>(tas::pT3_pt().size(), 1)); }));
0264   }
0265 
0266   bookDuplicateRateSets(list_DRSetDef);
0267 
0268   // Book Histograms
0269   ana.cutflow.bookHistograms(ana.histograms);  // if just want to book everywhere
0270 
0271   int nevts = 0;
0272 
0273   // Looping input file
0274   while (ana.looper.nextEvent()) {
0275     // If splitting jobs are requested then determine whether to process the event or not based on remainder
0276     if (ana.job_index != -1 and ana.nsplit_jobs != -1) {
0277       if (ana.looper.getNEventsProcessed() % ana.nsplit_jobs != (unsigned int)ana.job_index)
0278         continue;
0279     }
0280 
0281     // Reset all temporary variables necessary for histogramming
0282     ana.tx.clear();
0283 
0284     // Compute all temporary variables and pack the vector quantities that will get filled to the histograms
0285     fillEfficiencySets(list_effSetDef);
0286     fillFakeRateSets(list_FRSetDef);
0287     fillDuplicateRateSets(list_DRSetDef);
0288 
0289     // Reset all temporary variables necessary for histogramming
0290     ana.cutflow.fill();
0291 
0292     // Counting number of events processed
0293     nevts++;
0294   }
0295 
0296   // Write number of events processed
0297   ana.output_tfile->cd();
0298   TH1F* h_nevts = new TH1F("nevts", "nevts", 1, 0, 1);
0299   h_nevts->SetBinContent(1, nevts);
0300   h_nevts->Write();
0301 
0302   // Writing output file
0303   ana.cutflow.saveOutput();
0304 
0305   // The below can be sometimes crucial
0306   delete ana.output_tfile;
0307 
0308   // Done
0309   return 0;
0310 }
0311 
0312 // ---------------------------------------------------------=============================================-------------------------------------------------------------------
0313 // ---------------------------------------------------------=============================================-------------------------------------------------------------------
0314 // ---------------------------------------------------------=============================================-------------------------------------------------------------------
0315 // ---------------------------------------------------------=============================================-------------------------------------------------------------------
0316 // ---------------------------------------------------------=============================================-------------------------------------------------------------------
0317 
0318 //__________________________________________________________________________________________________________________________________________________________________________
0319 void bookEfficiencySets(std::vector<SimTrackSetDefinition>& effsets) {
0320   for (auto& effset : effsets)
0321     bookEfficiencySet(effset);
0322 }
0323 
0324 //__________________________________________________________________________________________________________________________________________________________________________
0325 void bookEfficiencySet(SimTrackSetDefinition& effset) {
0326   TString category_name = TString::Format("%s_%d_%d", effset.set_name.Data(), effset.pdgid, effset.q);
0327 
0328   // Denominator tracks' quantities
0329   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_denom_pt");
0330   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_denom_eta");
0331   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_denom_dxy");
0332   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_denom_vxy");
0333   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_denom_dz");
0334   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_denom_phi");
0335 
0336   // Numerator tracks' quantities
0337   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_numer_pt");
0338   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_numer_eta");
0339   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_numer_dxy");
0340   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_numer_vxy");
0341   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_numer_dz");
0342   ana.tx.createBranch<std::vector<float>>(category_name + "_ef_numer_phi");
0343 
0344   // Inefficiencies
0345   ana.tx.createBranch<std::vector<float>>(category_name + "_ie_numer_pt");
0346   ana.tx.createBranch<std::vector<float>>(category_name + "_ie_numer_eta");
0347   ana.tx.createBranch<std::vector<float>>(category_name + "_ie_numer_dxy");
0348   ana.tx.createBranch<std::vector<float>>(category_name + "_ie_numer_vxy");
0349   ana.tx.createBranch<std::vector<float>>(category_name + "_ie_numer_dz");
0350   ana.tx.createBranch<std::vector<float>>(category_name + "_ie_numer_phi");
0351 
0352   ana.histograms.addVecHistogram(category_name + "_ef_denom_pt", getPtBounds(0), [&, category_name]() {
0353     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_denom_pt");
0354   });
0355   ana.histograms.addVecHistogram(category_name + "_ef_denom_ptlow", getPtBounds(4), [&, category_name]() {
0356     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_denom_pt");
0357   });
0358   ana.histograms.addVecHistogram(category_name + "_ef_denom_ptmtv", getPtBounds(9), [&, category_name]() {
0359     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_denom_pt");
0360   });
0361   ana.histograms.addVecHistogram(category_name + "_ef_denom_ptflatbin", 180, 0., 100, [&, category_name]() {
0362     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_denom_pt");
0363   });
0364   ana.histograms.addVecHistogram(category_name + "_ef_denom_eta", 180, -4.5, 4.5, [&, category_name]() {
0365     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_denom_eta");
0366   });
0367   ana.histograms.addVecHistogram(category_name + "_ef_denom_dxy", 180, -30., 30., [&, category_name]() {
0368     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_denom_dxy");
0369   });
0370   ana.histograms.addVecHistogram(category_name + "_ef_denom_vxy", 180, -30., 30., [&, category_name]() {
0371     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_denom_vxy");
0372   });
0373   ana.histograms.addVecHistogram(category_name + "_ef_denom_dz", 180, -30., 30., [&, category_name]() {
0374     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_denom_dz");
0375   });
0376   ana.histograms.addVecHistogram(category_name + "_ef_denom_phi", 180, -M_PI, M_PI, [&, category_name]() {
0377     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_denom_phi");
0378   });
0379 
0380   ana.histograms.addVecHistogram(category_name + "_ef_numer_pt", getPtBounds(0), [&, category_name]() {
0381     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_numer_pt");
0382   });
0383   ana.histograms.addVecHistogram(category_name + "_ef_numer_ptlow", getPtBounds(4), [&, category_name]() {
0384     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_numer_pt");
0385   });
0386   ana.histograms.addVecHistogram(category_name + "_ef_numer_ptmtv", getPtBounds(9), [&, category_name]() {
0387     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_numer_pt");
0388   });
0389   ana.histograms.addVecHistogram(category_name + "_ef_numer_ptflatbin", 180, 0., 100, [&, category_name]() {
0390     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_numer_pt");
0391   });
0392   ana.histograms.addVecHistogram(category_name + "_ef_numer_eta", 180, -4.5, 4.5, [&, category_name]() {
0393     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_numer_eta");
0394   });
0395   ana.histograms.addVecHistogram(category_name + "_ef_numer_dxy", 180, -30., 30., [&, category_name]() {
0396     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_numer_dxy");
0397   });
0398   ana.histograms.addVecHistogram(category_name + "_ef_numer_vxy", 180, -30., 30., [&, category_name]() {
0399     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_numer_vxy");
0400   });
0401   ana.histograms.addVecHistogram(category_name + "_ef_numer_dz", 180, -30., 30., [&, category_name]() {
0402     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_numer_dz");
0403   });
0404   ana.histograms.addVecHistogram(category_name + "_ef_numer_phi", 180, -M_PI, M_PI, [&, category_name]() {
0405     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ef_numer_phi");
0406   });
0407 
0408   ana.histograms.addVecHistogram(category_name + "_ie_numer_pt", getPtBounds(0), [&, category_name]() {
0409     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ie_numer_pt");
0410   });
0411   ana.histograms.addVecHistogram(category_name + "_ie_numer_ptlow", getPtBounds(4), [&, category_name]() {
0412     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ie_numer_pt");
0413   });
0414   ana.histograms.addVecHistogram(category_name + "_ie_numer_ptmtv", getPtBounds(9), [&, category_name]() {
0415     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ie_numer_pt");
0416   });
0417   ana.histograms.addVecHistogram(category_name + "_ie_numer_ptflatbin", 180, 0., 100, [&, category_name]() {
0418     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ie_numer_pt");
0419   });
0420   ana.histograms.addVecHistogram(category_name + "_ie_numer_eta", 180, -4.5, 4.5, [&, category_name]() {
0421     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ie_numer_eta");
0422   });
0423   ana.histograms.addVecHistogram(category_name + "_ie_numer_dxy", 180, -30., 30., [&, category_name]() {
0424     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ie_numer_dxy");
0425   });
0426   ana.histograms.addVecHistogram(category_name + "_ie_numer_vxy", 180, -30., 30., [&, category_name]() {
0427     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ie_numer_vxy");
0428   });
0429   ana.histograms.addVecHistogram(category_name + "_ie_numer_dz", 180, -30., 30., [&, category_name]() {
0430     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ie_numer_dz");
0431   });
0432   ana.histograms.addVecHistogram(category_name + "_ie_numer_phi", 180, -M_PI, M_PI, [&, category_name]() {
0433     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_ie_numer_phi");
0434   });
0435 }
0436 
0437 //__________________________________________________________________________________________________________________________________________________________________________
0438 void bookDuplicateRateSets(std::vector<RecoTrackSetDefinition>& DRsets) {
0439   for (auto& DRset : DRsets) {
0440     bookDuplicateRateSet(DRset);
0441   }
0442 }
0443 
0444 //__________________________________________________________________________________________________________________________________________________________________________
0445 void bookDuplicateRateSet(RecoTrackSetDefinition& DRset) {
0446   TString category_name = DRset.set_name;
0447 
0448   // Denominator tracks' quantities
0449   ana.tx.createBranch<std::vector<float>>(category_name + "_dr_denom_pt");
0450   ana.tx.createBranch<std::vector<float>>(category_name + "_dr_denom_eta");
0451   ana.tx.createBranch<std::vector<float>>(category_name + "_dr_denom_phi");
0452 
0453   // Numerator tracks' quantities
0454   ana.tx.createBranch<std::vector<float>>(category_name + "_dr_numer_pt");
0455   ana.tx.createBranch<std::vector<float>>(category_name + "_dr_numer_eta");
0456   ana.tx.createBranch<std::vector<float>>(category_name + "_dr_numer_phi");
0457 
0458   // Histogram utility object that is used to define the histograms
0459   ana.histograms.addVecHistogram(category_name + "_dr_denom_pt", getPtBounds(0), [&, category_name]() {
0460     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_dr_denom_pt");
0461   });
0462   ana.histograms.addVecHistogram(category_name + "_dr_denom_ptlow", getPtBounds(4), [&, category_name]() {
0463     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_dr_denom_pt");
0464   });
0465   ana.histograms.addVecHistogram(category_name + "_dr_denom_ptmtv", getPtBounds(9), [&, category_name]() {
0466     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_dr_denom_pt");
0467   });
0468   ana.histograms.addVecHistogram(category_name + "_dr_denom_eta", 180, -4.5, 4.5, [&, category_name]() {
0469     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_dr_denom_eta");
0470   });
0471   ana.histograms.addVecHistogram(category_name + "_dr_denom_phi", 180, -M_PI, M_PI, [&, category_name]() {
0472     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_dr_denom_phi");
0473   });
0474   ana.histograms.addVecHistogram(category_name + "_dr_numer_pt", getPtBounds(0), [&, category_name]() {
0475     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_dr_numer_pt");
0476   });
0477   ana.histograms.addVecHistogram(category_name + "_dr_numer_ptlow", getPtBounds(4), [&, category_name]() {
0478     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_dr_numer_pt");
0479   });
0480   ana.histograms.addVecHistogram(category_name + "_dr_numer_ptmtv", getPtBounds(9), [&, category_name]() {
0481     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_dr_numer_pt");
0482   });
0483   ana.histograms.addVecHistogram(category_name + "_dr_numer_eta", 180, -4.5, 4.5, [&, category_name]() {
0484     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_dr_numer_eta");
0485   });
0486   ana.histograms.addVecHistogram(category_name + "_dr_numer_phi", 180, -M_PI, M_PI, [&, category_name]() {
0487     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_dr_numer_phi");
0488   });
0489 }
0490 
0491 //__________________________________________________________________________________________________________________________________________________________________________
0492 void bookFakeRateSets(std::vector<RecoTrackSetDefinition>& FRsets) {
0493   for (auto& FRset : FRsets) {
0494     bookFakeRateSet(FRset);
0495   }
0496 }
0497 
0498 //__________________________________________________________________________________________________________________________________________________________________________
0499 void bookFakeRateSet(RecoTrackSetDefinition& FRset) {
0500   TString category_name = FRset.set_name;
0501 
0502   // Denominator tracks' quantities
0503   ana.tx.createBranch<std::vector<float>>(category_name + "_fr_denom_pt");
0504   ana.tx.createBranch<std::vector<float>>(category_name + "_fr_denom_eta");
0505   ana.tx.createBranch<std::vector<float>>(category_name + "_fr_denom_phi");
0506 
0507   // Numerator tracks' quantities
0508   ana.tx.createBranch<std::vector<float>>(category_name + "_fr_numer_pt");
0509   ana.tx.createBranch<std::vector<float>>(category_name + "_fr_numer_eta");
0510   ana.tx.createBranch<std::vector<float>>(category_name + "_fr_numer_phi");
0511 
0512   // Histogram utility object that is used to define the histograms
0513   ana.histograms.addVecHistogram(category_name + "_fr_denom_pt", getPtBounds(0), [&, category_name]() {
0514     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_fr_denom_pt");
0515   });
0516   ana.histograms.addVecHistogram(category_name + "_fr_denom_ptlow", getPtBounds(4), [&, category_name]() {
0517     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_fr_denom_pt");
0518   });
0519   ana.histograms.addVecHistogram(category_name + "_fr_denom_ptmtv", getPtBounds(9), [&, category_name]() {
0520     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_fr_denom_pt");
0521   });
0522   ana.histograms.addVecHistogram(category_name + "_fr_denom_eta", 180, -4.5, 4.5, [&, category_name]() {
0523     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_fr_denom_eta");
0524   });
0525   ana.histograms.addVecHistogram(category_name + "_fr_numer_phi", 180, -M_PI, M_PI, [&, category_name]() {
0526     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_fr_numer_phi");
0527   });
0528   ana.histograms.addVecHistogram(category_name + "_fr_numer_pt", getPtBounds(0), [&, category_name]() {
0529     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_fr_numer_pt");
0530   });
0531   ana.histograms.addVecHistogram(category_name + "_fr_numer_ptlow", getPtBounds(4), [&, category_name]() {
0532     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_fr_numer_pt");
0533   });
0534   ana.histograms.addVecHistogram(category_name + "_fr_numer_ptmtv", getPtBounds(9), [&, category_name]() {
0535     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_fr_numer_pt");
0536   });
0537   ana.histograms.addVecHistogram(category_name + "_fr_numer_eta", 180, -4.5, 4.5, [&, category_name]() {
0538     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_fr_numer_eta");
0539   });
0540   ana.histograms.addVecHistogram(category_name + "_fr_denom_phi", 180, -M_PI, M_PI, [&, category_name]() {
0541     return ana.tx.getBranchLazy<std::vector<float>>(category_name + "_fr_denom_phi");
0542   });
0543 }
0544 
0545 // ---------------------------------------------------------=============================================-------------------------------------------------------------------
0546 // ---------------------------------------------------------=============================================-------------------------------------------------------------------
0547 // ---------------------------------------------------------=============================================-------------------------------------------------------------------
0548 // ---------------------------------------------------------=============================================-------------------------------------------------------------------
0549 // ---------------------------------------------------------=============================================-------------------------------------------------------------------
0550 
0551 //__________________________________________________________________________________________________________________________________________________________________________
0552 void fillEfficiencySets(std::vector<SimTrackSetDefinition>& effsets) {
0553   for (auto& effset : effsets) {
0554     for (unsigned int isimtrk = 0; isimtrk < lstEff.sim_pt().size(); ++isimtrk) {
0555       fillEfficiencySet(isimtrk, effset);
0556     }
0557   }
0558 }
0559 
0560 //__________________________________________________________________________________________________________________________________________________________________________
0561 void fillEfficiencySet(int isimtrk, SimTrackSetDefinition& effset) {
0562   //=========================================================
0563   // NOTE: The following is not applied as the LSTNtuple no longer writes this.
0564   // const int &bunch = lstEff.sim_bunchCrossing()[isimtrk];
0565   // const int &event = lstEff.sim_event()[isimtrk];
0566   // if (bunch != 0)
0567   //     return;
0568   // if (event != 0)
0569   //     return;
0570   //=========================================================
0571 
0572   const float& pt = lstEff.sim_pt()[isimtrk];
0573   const float& eta = lstEff.sim_eta()[isimtrk];
0574   const float& dz = lstEff.sim_pca_dz()[isimtrk];
0575   const float& dxy = lstEff.sim_pca_dxy()[isimtrk];
0576   const float& phi = lstEff.sim_phi()[isimtrk];
0577   const int& pdgidtrk = lstEff.sim_pdgId()[isimtrk];
0578   const int& q = lstEff.sim_q()[isimtrk];
0579   const float& vtx_x = lstEff.sim_vx()[isimtrk];
0580   const float& vtx_y = lstEff.sim_vy()[isimtrk];
0581   const float& vtx_z = lstEff.sim_vz()[isimtrk];
0582   const float& vtx_perp = sqrt(vtx_x * vtx_x + vtx_y * vtx_y);
0583   bool pass = effset.pass(isimtrk);
0584   bool sel = effset.sel(isimtrk);
0585 
0586   if (effset.pdgid != 0) {
0587     if (abs(pdgidtrk) != abs(effset.pdgid))
0588       return;
0589   }
0590 
0591   if (effset.q != 0) {
0592     if (q != effset.q)
0593       return;
0594   }
0595 
0596   if (effset.pdgid == 0 and q == 0)
0597     return;
0598 
0599   if (not sel)
0600     return;
0601 
0602   TString category_name = TString::Format("%s_%d_%d", effset.set_name.Data(), effset.pdgid, effset.q);
0603 
0604   // https://github.com/cms-sw/cmssw/blob/7cbdb18ec6d11d5fd17ca66c1153f0f4e869b6b0/SimTracker/Common/python/trackingParticleSelector_cfi.py
0605   // https://github.com/cms-sw/cmssw/blob/7cbdb18ec6d11d5fd17ca66c1153f0f4e869b6b0/SimTracker/Common/interface/TrackingParticleSelector.h#L122-L124
0606   const float vtx_z_thresh = 30;
0607   const float vtx_perp_thresh = 2.5;
0608 
0609   // N minus eta cut
0610   if (pt > ana.pt_cut and abs(vtx_z) < vtx_z_thresh and abs(vtx_perp) < vtx_perp_thresh) {
0611     // vs. eta plot
0612     ana.tx.pushbackToBranch<float>(category_name + "_ef_denom_eta", eta);
0613     if (pass)
0614       ana.tx.pushbackToBranch<float>(category_name + "_ef_numer_eta", eta);
0615     else
0616       ana.tx.pushbackToBranch<float>(category_name + "_ie_numer_eta", eta);
0617   }
0618 
0619   // N minus pt cut
0620   if (abs(eta) < ana.eta_cut and abs(vtx_z) < vtx_z_thresh and abs(vtx_perp) < vtx_perp_thresh) {
0621     // vs. pt plot
0622     ana.tx.pushbackToBranch<float>(category_name + "_ef_denom_pt", pt);
0623     if (pass)
0624       ana.tx.pushbackToBranch<float>(category_name + "_ef_numer_pt", pt);
0625     else
0626       ana.tx.pushbackToBranch<float>(category_name + "_ie_numer_pt", pt);
0627   }
0628 
0629   // N minus dxy cut
0630   if (abs(eta) < ana.eta_cut and pt > ana.pt_cut and abs(vtx_z) < vtx_z_thresh) {
0631     // vs. dxy plot
0632     ana.tx.pushbackToBranch<float>(category_name + "_ef_denom_dxy", dxy);
0633     ana.tx.pushbackToBranch<float>(category_name + "_ef_denom_vxy", vtx_perp);
0634     if (pass) {
0635       ana.tx.pushbackToBranch<float>(category_name + "_ef_numer_dxy", dxy);
0636       ana.tx.pushbackToBranch<float>(category_name + "_ef_numer_vxy", vtx_perp);
0637     } else {
0638       ana.tx.pushbackToBranch<float>(category_name + "_ie_numer_dxy", dxy);
0639       ana.tx.pushbackToBranch<float>(category_name + "_ie_numer_vxy", vtx_perp);
0640     }
0641   }
0642 
0643   // N minus dz cut
0644   if (abs(eta) < ana.eta_cut and pt > ana.pt_cut and abs(vtx_perp) < vtx_perp_thresh) {
0645     // vs. dz plot
0646     ana.tx.pushbackToBranch<float>(category_name + "_ef_denom_dz", dz);
0647     if (pass)
0648       ana.tx.pushbackToBranch<float>(category_name + "_ef_numer_dz", dz);
0649     else
0650       ana.tx.pushbackToBranch<float>(category_name + "_ie_numer_dz", dz);
0651   }
0652 
0653   // All phase-space cuts
0654   if (abs(eta) < ana.eta_cut and pt > ana.pt_cut and abs(vtx_z) < vtx_z_thresh and abs(vtx_perp) < vtx_perp_thresh) {
0655     // vs. Phi plot
0656     ana.tx.pushbackToBranch<float>(category_name + "_ef_denom_phi", phi);
0657     if (pass)
0658       ana.tx.pushbackToBranch<float>(category_name + "_ef_numer_phi", phi);
0659     else
0660       ana.tx.pushbackToBranch<float>(category_name + "_ie_numer_phi", phi);
0661   }
0662 }
0663 
0664 //__________________________________________________________________________________________________________________________________________________________________________
0665 void fillFakeRateSets(std::vector<RecoTrackSetDefinition>& FRsets) {
0666   for (auto& FRset : FRsets) {
0667     for (unsigned int itc = 0; itc < FRset.pt().size(); ++itc) {
0668       fillFakeRateSet(itc, FRset);
0669     }
0670   }
0671 }
0672 
0673 //__________________________________________________________________________________________________________________________________________________________________________
0674 void fillFakeRateSet(int itc, RecoTrackSetDefinition& FRset) {
0675   float pt = FRset.pt().at(itc);
0676   float eta = FRset.eta().at(itc);
0677   float phi = FRset.phi().at(itc);
0678   TString category_name = FRset.set_name;
0679   bool pass = FRset.pass(itc);
0680   bool sel = FRset.sel(itc);
0681 
0682   if (not sel)
0683     return;
0684 
0685   if (pt > ana.pt_cut) {
0686     ana.tx.pushbackToBranch<float>(category_name + "_fr_denom_eta", eta);
0687     if (pass)
0688       ana.tx.pushbackToBranch<float>(category_name + "_fr_numer_eta", eta);
0689   }
0690   if (abs(eta) < ana.eta_cut) {
0691     ana.tx.pushbackToBranch<float>(category_name + "_fr_denom_pt", pt);
0692     if (pass)
0693       ana.tx.pushbackToBranch<float>(category_name + "_fr_numer_pt", pt);
0694   }
0695   if (abs(eta) < ana.eta_cut and pt > ana.pt_cut) {
0696     ana.tx.pushbackToBranch<float>(category_name + "_fr_denom_phi", phi);
0697     if (pass)
0698       ana.tx.pushbackToBranch<float>(category_name + "_fr_numer_phi", phi);
0699   }
0700 }
0701 
0702 //__________________________________________________________________________________________________________________________________________________________________________
0703 void fillDuplicateRateSets(std::vector<RecoTrackSetDefinition>& DRsets) {
0704   for (auto& DRset : DRsets) {
0705     for (unsigned int itc = 0; itc < DRset.pt().size(); ++itc) {
0706       fillDuplicateRateSet(itc, DRset);
0707     }
0708   }
0709 }
0710 
0711 //__________________________________________________________________________________________________________________________________________________________________________
0712 void fillDuplicateRateSet(int itc, RecoTrackSetDefinition& DRset) {
0713   float pt = DRset.pt().at(itc);
0714   float eta = DRset.eta().at(itc);
0715   float phi = DRset.phi().at(itc);
0716   TString category_name = DRset.set_name;
0717   bool pass = DRset.pass(itc);
0718   bool sel = DRset.sel(itc);
0719 
0720   if (not sel)
0721     return;
0722 
0723   if (pt > ana.pt_cut) {
0724     ana.tx.pushbackToBranch<float>(category_name + "_dr_denom_eta", eta);
0725     if (pass)
0726       ana.tx.pushbackToBranch<float>(category_name + "_dr_numer_eta", eta);
0727   }
0728   if (abs(eta) < ana.eta_cut) {
0729     ana.tx.pushbackToBranch<float>(category_name + "_dr_denom_pt", pt);
0730     if (pass)
0731       ana.tx.pushbackToBranch<float>(category_name + "_dr_numer_pt", pt);
0732   }
0733   if (abs(eta) < ana.eta_cut and pt > ana.pt_cut) {
0734     ana.tx.pushbackToBranch<float>(category_name + "_dr_denom_phi", phi);
0735     if (pass)
0736       ana.tx.pushbackToBranch<float>(category_name + "_dr_numer_phi", phi);
0737   }
0738 }