File indexing completed on 2023-03-17 11:22:38
0001 #include "PlotValidation.hh"
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031 PlotValidation::PlotValidation(const TString& inName,
0032 const TString& outName,
0033 const Bool_t cmsswComp,
0034 const int algo,
0035 const Bool_t mvInput,
0036 const Bool_t rmSuffix,
0037 const Bool_t saveAs,
0038 const TString& outType)
0039 : fInName(inName),
0040 fOutName(outName),
0041 fCmsswComp(cmsswComp),
0042 fAlgo(algo),
0043 fMvInput(mvInput),
0044 fRmSuffix(rmSuffix),
0045 fSaveAs(saveAs),
0046 fOutType(outType) {
0047
0048 PlotValidation::SetupStyle();
0049 if (fAlgo > 0 && !fRmSuffix)
0050 fOutName = fOutName + "_iter" + algo;
0051 PlotValidation::MakeOutDir(fOutName);
0052 PlotValidation::SetupBins();
0053 PlotValidation::SetupCommonVars();
0054
0055
0056 fInRoot = TFile::Open(fInName.Data());
0057 if (fInRoot == (TFile*)NULL) {
0058 std::cerr << "File: " << fInName.Data() << " does not exist!!! Exiting..." << std::endl;
0059 exit(1);
0060 }
0061 gROOT->cd();
0062 efftree = (TTree*)fInRoot->Get((fCmsswComp ? "cmsswefftree" : "efftree"));
0063 frtree = (TTree*)fInRoot->Get((fCmsswComp ? "cmsswfrtree" : "frtree"));
0064 if (algo > 0)
0065 frtree = frtree->CopyTree(Form("algorithm==%i", algo));
0066
0067 fOutRoot = new TFile(fOutName + "/plots.root", "RECREATE");
0068 }
0069
0070 PlotValidation::~PlotValidation() {
0071 delete efftree;
0072 delete frtree;
0073 delete fInRoot;
0074 delete fOutRoot;
0075 }
0076
0077 void PlotValidation::Validation(int algo) {
0078 std::cout << "Computing Efficiency, Inefficiency, and Duplicate Rate ..." << std::endl;
0079 PlotValidation::PlotEffTree(algo);
0080
0081 std::cout << "Computing Fake Rate, <nHits/track>, and kinematic diffs to " << fSRefTitle.Data() << " tracks ..."
0082 << std::endl;
0083 PlotValidation::PlotFRTree(algo);
0084
0085 std::cout << "Printing Totals ..." << std::endl;
0086 PlotValidation::PrintTotals(algo);
0087
0088 if (fMvInput)
0089 PlotValidation::MoveInput();
0090 }
0091
0092
0093 void PlotValidation::PlotEffTree(int algo) {
0094
0095
0096
0097
0098 const TStrVec rates = {"eff", "dr", "ineff"};
0099 const TStrVec srates = {"Efficiency", "Duplicate Rate", "Inefficiency"};
0100 const UInt_t nrates = rates.size();
0101
0102 const TStrVec regs = {"brl", "trans", "ec"};
0103 const TStrVec sregs = {"Barrel", "Transition", "Endcap"};
0104 const FltVec etacuts = {0, 0.9, 1.7, 2.45};
0105 const UInt_t nregs = regs.size();
0106
0107
0108
0109
0110
0111 TEffRefMap plots;
0112 for (auto i = 0U; i < fNVars; i++)
0113 {
0114 const auto& var = fVars[i];
0115 const auto& svar = fSVars[i];
0116 const auto& sunit = fSUnits[i];
0117
0118
0119 const auto& varbins = fVarBins[i];
0120 const Double_t* bins = &varbins[0];
0121
0122 for (auto j = 0U; j < fNTrks; j++)
0123 {
0124 const auto& trk = fTrks[j];
0125 const auto& strk = fSTrks[j];
0126
0127 for (auto k = 0U; k < fNPtCuts; k++)
0128 {
0129 const auto& sptcut = fSPtCuts[k];
0130 const auto& hptcut = fHPtCuts[k];
0131
0132 for (auto l = 0U; l < nrates; l++)
0133 {
0134 const auto& rate = rates[l];
0135 const auto& srate = srates[l];
0136
0137
0138 const TString plotkey = Form("%i_%i_%i_%i", i, j, k, l);
0139 const TString plotname = Form("%s_", fCmsswComp ? "cmssw" : "sim") + var + "_" + trk + "_pt" + hptcut;
0140 const TString plottitle = strk + " Track " + srate + " vs " + fSRefTitle + " " + svar + " {" + fSVarPt +
0141 " > " + sptcut + " " + fSUnitPt + "};" + svar + sunit + ";" + srate;
0142
0143
0144 if (l < 2) {
0145 const TString tmpname = rate + "_" + plotname;
0146 plots[plotkey] = new TEfficiency(tmpname.Data(), plottitle.Data(), varbins.size() - 1, bins);
0147 } else
0148 {
0149 for (auto m = 0U; m < nregs; m++)
0150 {
0151 const auto& reg = regs[m];
0152 const auto& sreg = sregs[m];
0153
0154 const TString tmpkey = Form("%s_%i", plotkey.Data(), m);
0155 const TString tmpname = rate + "_" + reg + "_" + plotname;
0156 const TString tmptitle = strk + " Track " + srate + " vs " + fSRefTitle + " " + svar + "{" + fSVarPt +
0157 " > " + sptcut + " " + fSUnitPt + ", " + sreg + "};" + svar + sunit + ";" +
0158 srate;
0159
0160 plots[tmpkey] = new TEfficiency(tmpname.Data(), tmptitle.Data(), varbins.size() - 1, bins);
0161 }
0162 }
0163 }
0164 }
0165 }
0166 }
0167
0168
0169
0170
0171
0172
0173 FltVec vars_ref(fNVars);
0174 TBrRefVec vars_ref_br(fNVars);
0175 for (auto i = 0U; i < fNVars; i++)
0176 {
0177 const auto& var = fVars[i];
0178 auto& var_ref = vars_ref[i];
0179 auto& var_ref_br = vars_ref_br[i];
0180
0181
0182 var_ref = 0.;
0183 var_ref_br = 0;
0184
0185
0186 efftree->SetBranchAddress(
0187 var + "_" + ((fSRefVar == "cmssw" || var != "nLayers") ? fSRefVar : fSRefVarTrk), &var_ref, &var_ref_br);
0188 }
0189
0190
0191 IntVec refmask_trks(fNTrks);
0192 TBrRefVec refmask_trks_br(fNTrks);
0193
0194 IntVec duplmask_trks(fNTrks);
0195 TBrRefVec duplmask_trks_br(fNTrks);
0196
0197 std::vector<ULong64_t> itermask_trks(fNTrks);
0198 TBrRefVec itermask_trks_br(fNTrks);
0199
0200 std::vector<ULong64_t> iterduplmask_trks(fNTrks);
0201 TBrRefVec iterduplmask_trks_br(fNTrks);
0202
0203 ULong64_t algoseed_trk;
0204 TBranch* algoseed_trk_br;
0205
0206 for (auto j = 0U; j < fNTrks; j++)
0207 {
0208 const auto& trk = fTrks[j];
0209 auto& refmask_trk = refmask_trks[j];
0210 auto& refmask_trk_br = refmask_trks_br[j];
0211 auto& duplmask_trk = duplmask_trks[j];
0212 auto& duplmask_trk_br = duplmask_trks_br[j];
0213 auto& itermask_trk = itermask_trks[j];
0214 auto& itermask_trk_br = itermask_trks_br[j];
0215
0216 auto& iterduplmask_trk = iterduplmask_trks[j];
0217 auto& iterduplmask_trk_br = iterduplmask_trks_br[j];
0218
0219
0220 refmask_trk = 0;
0221 refmask_trk_br = 0;
0222
0223
0224 duplmask_trk = 0;
0225 duplmask_trk_br = 0;
0226
0227
0228 itermask_trk = 0;
0229 itermask_trk_br = 0;
0230
0231 iterduplmask_trk = 0;
0232 iterduplmask_trk_br = 0;
0233
0234 algoseed_trk = 0;
0235 algoseed_trk_br = 0;
0236
0237
0238 efftree->SetBranchAddress(fSRefMask + "mask_" + trk, &refmask_trk, &refmask_trk_br);
0239 efftree->SetBranchAddress("duplmask_" + trk, &duplmask_trk, &duplmask_trk_br);
0240 efftree->SetBranchAddress("itermask_" + trk, &itermask_trk, &itermask_trk_br);
0241 efftree->SetBranchAddress("iterduplmask_" + trk, &iterduplmask_trk, &iterduplmask_trk_br);
0242 }
0243 efftree->SetBranchAddress("algo_seed", &algoseed_trk, &algoseed_trk_br);
0244
0245
0246
0247
0248
0249 const auto nentries = efftree->GetEntries();
0250 for (auto e = 0U; e < nentries; e++) {
0251
0252 for (auto i = 0U; i < fNVars; i++) {
0253 auto& var_ref_br = vars_ref_br[i];
0254
0255 var_ref_br->GetEntry(e);
0256 }
0257 for (auto j = 0U; j < fNTrks; j++) {
0258 auto& refmask_trk_br = refmask_trks_br[j];
0259 auto& duplmask_trk_br = duplmask_trks_br[j];
0260 auto& itermask_trk_br = itermask_trks_br[j];
0261 auto& iterduplmask_trk_br = iterduplmask_trks_br[j];
0262
0263 refmask_trk_br->GetEntry(e);
0264 duplmask_trk_br->GetEntry(e);
0265 itermask_trk_br->GetEntry(e);
0266 iterduplmask_trk_br->GetEntry(e);
0267 }
0268 algoseed_trk_br->GetEntry(e);
0269
0270 const auto pt_ref = vars_ref[0];
0271
0272
0273 for (auto k = 0U; k < fNPtCuts; k++)
0274 {
0275 const auto ptcut = fPtCuts[k];
0276
0277 if (pt_ref < ptcut)
0278 continue;
0279
0280 for (auto i = 0U; i < fNVars; i++)
0281 {
0282 const auto var_ref = vars_ref[i];
0283
0284 for (auto j = 0U; j < fNTrks; j++)
0285 {
0286 const auto refmask_trk = refmask_trks[j];
0287 const auto duplmask_trk = duplmask_trks[j];
0288 const auto itermask_trk = itermask_trks[j];
0289 const auto iterduplmask_trk = iterduplmask_trks[j];
0290
0291 const auto effIteration = algo > 0 ? ((itermask_trk >> algo) & 1) : 1;
0292 const auto oneIteration = algo > 0 ? ((iterduplmask_trk >> algo) & 1) : 1;
0293 const auto ineffIteration = algo > 0 ? (((itermask_trk >> algo) & 1) == 0) : (refmask_trk == 0);
0294 const auto seedalgo_flag = (algoseed_trk > 0 && algo > 0) ? ((algoseed_trk >> algo) & 1) : 1;
0295
0296
0297 const TString basekey = Form("%i_%i_%i", i, j, k);
0298
0299
0300 if (refmask_trk != -1 && seedalgo_flag)
0301 plots[basekey + "_0"]->Fill((refmask_trk == 1) && effIteration,
0302 var_ref);
0303
0304
0305 if (duplmask_trk != -1 && effIteration && seedalgo_flag)
0306 plots[basekey + "_1"]->Fill((duplmask_trk == 1) && oneIteration,
0307 var_ref);
0308
0309
0310 if (refmask_trk != -1) {
0311 for (auto m = 0U; m < regs.size(); m++) {
0312 const auto eta_ref = std::abs(vars_ref[1]);
0313 const auto etalow = etacuts[m];
0314 const auto etaup = etacuts[m + 1];
0315
0316
0317 if ((eta_ref >= etalow) && (eta_ref < etaup))
0318 plots[Form("%s_2_%i", basekey.Data(), m)]->Fill(ineffIteration, var_ref);
0319 }
0320 }
0321
0322 }
0323 }
0324 }
0325 }
0326
0327
0328
0329
0330
0331
0332 TStrVec dirnames = {"efficiency", "duplicaterate", "inefficiency"};
0333 for (auto& dirname : dirnames)
0334 dirname += fSRefDir;
0335
0336 TDirRefVec subdirs(nrates);
0337 for (auto l = 0U; l < nrates; l++)
0338 subdirs[l] = PlotValidation::MakeSubDirs(dirnames[l]);
0339
0340
0341 for (auto i = 0U; i < fNVars; i++) {
0342 for (auto j = 0U; j < fNTrks; j++) {
0343 for (auto k = 0U; k < fNPtCuts; k++) {
0344 for (auto l = 0U; l < nrates; l++) {
0345 const auto& dirname = dirnames[l];
0346 auto& subdir = subdirs[l];
0347
0348 const TString plotkey = Form("%i_%i_%i_%i", i, j, k, l);
0349 if (l < 2)
0350 {
0351 PlotValidation::DrawWriteSavePlot(plots[plotkey], subdir, dirname, "AP");
0352 delete plots[plotkey];
0353 } else {
0354 for (auto m = 0U; m < nregs; m++) {
0355 const TString tmpkey = Form("%s_%i", plotkey.Data(), m);
0356 PlotValidation::DrawWriteSavePlot(plots[tmpkey], subdir, dirname, "AP");
0357 delete plots[tmpkey];
0358 }
0359 }
0360 }
0361 }
0362 }
0363 }
0364 }
0365
0366
0367 void PlotValidation::PlotFRTree(int algo) {
0368
0369
0370
0371
0372
0373 const TStrVec colls = {"allreco", "fake", "allmatch", "bestmatch"};
0374 const TStrVec scolls = {"All Reco", "Fake", "All Match", "Best Match"};
0375 const UInt_t ncolls = colls.size();
0376
0377
0378 const DblVecVec trkqualbins = {fNHitsBins, fFracHitsBins, fScoreBins};
0379
0380
0381 const TStrVec dvars = {"dnHits", "dinvpt", "deta", "dphi"};
0382 const TStrVec sdvars = {"nHits", "1/p_{T}", "#eta", "#phi"};
0383 const UInt_t ndvars = dvars.size();
0384
0385
0386 const DblVecVec dvarbins = {fDNHitsBins, fDInvPtBins, fDEtaBins, fDPhiBins};
0387
0388
0389
0390
0391
0392 TEffRefMap plots;
0393 TH1FRefMap hists;
0394 for (auto j = 0U; j < fNTrks; j++)
0395 {
0396 const auto& trk = fTrks[j];
0397 const auto& strk = fSTrks[j];
0398
0399 for (auto k = 0U; k < fNPtCuts; k++)
0400 {
0401 const auto& sptcut = fSPtCuts[k];
0402 const auto& hptcut = fHPtCuts[k];
0403
0404
0405 for (auto i = 0U; i < fNVars; i++)
0406 {
0407 const auto& var = fVars[i];
0408 const auto& svar = fSVars[i];
0409 const auto& sunit = fSUnits[i];
0410
0411
0412 const TString plotkey = Form("%i_%i_%i", i, j, k);
0413 const TString plotname = "fr_reco_" + var + "_" + trk + "_pt" + hptcut;
0414 const TString plottitle = strk + " Track Fake Rate vs Reco " + svar + " {" + fSVarPt + " > " + sptcut + " " +
0415 fSUnitPt + "};" + svar + sunit + ";Fake Rate";
0416
0417
0418 const auto& varbins = fVarBins[i];
0419 const Double_t* bins = &varbins[0];
0420
0421 plots[plotkey] = new TEfficiency(plotname.Data(), plottitle.Data(), varbins.size() - 1, bins);
0422 }
0423
0424
0425 for (auto n = 0U; n < fNTrkQual; n++)
0426 {
0427 const auto& trkqual = fTrkQual[n];
0428 const auto& strkqual = fSTrkQual[n];
0429
0430
0431 const auto& varbins = trkqualbins[n];
0432 const Double_t* bins = &varbins[0];
0433
0434 for (auto o = 0U; o < ncolls; o++)
0435 {
0436 const auto& coll = colls[o];
0437 const auto& scoll = scolls[o];
0438
0439
0440 const TString histkey = Form("%i_%i_%i_%i", j, k, n, o);
0441 const TString histname = "h_" + trkqual + "_" + coll + "_" + trk + "_pt" + hptcut;
0442 const TString histtitle = scoll + " " + strk + " Track vs " + strkqual + " {" + fSVarPt + " > " + sptcut +
0443 " " + fSUnitPt + "};" + strkqual + ";nTracks";
0444
0445
0446 hists[histkey] = new TH1F(histname.Data(), histtitle.Data(), varbins.size() - 1, bins);
0447 hists[histkey]->Sumw2();
0448 }
0449 }
0450
0451
0452 for (auto p = 0U; p < ndvars; p++)
0453 {
0454 const auto& dvar = dvars[p];
0455 const auto& sdvar = sdvars[p];
0456
0457
0458 const auto& varbins = dvarbins[p];
0459 const Double_t* bins = &varbins[0];
0460
0461
0462 for (auto o = 2U; o < ncolls; o++) {
0463 const auto& coll = colls[o];
0464 const auto& scoll = scolls[o];
0465
0466
0467 const TString histkey = Form("%i_%i_d_%i_%i", j, k, p, o);
0468 const TString histname = "h_" + dvar + "_" + coll + "_" + trk + "_pt" + hptcut;
0469 const TString histtitle = "#Delta" + sdvar + "(" + scoll + " " + strk + "," + fSRefTitle + ") {" + fSVarPt +
0470 " > " + sptcut + " " + fSUnitPt + "};" + sdvar + "^{" + scoll + " " + strk + "}-" +
0471 sdvar + "^{" + fSRefTitle + "};nTracks";
0472
0473
0474 hists[histkey] = new TH1F(histname.Data(), histtitle.Data(), varbins.size() - 1, bins);
0475 hists[histkey]->Sumw2();
0476 }
0477 }
0478
0479 }
0480 }
0481
0482
0483
0484
0485
0486
0487 FltVecVec vars_trks(fNVars);
0488 TBrRefVecVec vars_trks_br(fNVars);
0489 for (auto i = 0U; i < fNVars; i++)
0490 {
0491 const auto& var = fVars[i];
0492 auto& var_trks = vars_trks[i];
0493 auto& var_trks_br = vars_trks_br[i];
0494
0495 var_trks.resize(fNTrks);
0496 var_trks_br.resize(fNTrks);
0497
0498 for (auto j = 0U; j < fNTrks; j++)
0499 {
0500 const auto& trk = fTrks[j];
0501 auto& var_trk = var_trks[j];
0502 auto& var_trk_br = var_trks_br[j];
0503
0504
0505 var_trk = 0.;
0506 var_trk_br = 0;
0507
0508
0509 frtree->SetBranchAddress(var + "_" + trk, &var_trk, &var_trk_br);
0510 }
0511 }
0512
0513
0514 IntVec refmask_trks(fNTrks);
0515 TBrRefVec refmask_trks_br(fNTrks);
0516 IntVec iTkMatches_trks(fNTrks);
0517 TBrRefVec iTkMatches_trks_br(fNTrks);
0518
0519
0520 IntVec nHits_trks(fNTrks);
0521 TBrRefVec nHits_trks_br(fNTrks);
0522 FltVec fracHits_trks(fNTrks);
0523 TBrRefVec fracHits_trks_br(fNTrks);
0524 IntVec score_trks(fNTrks);
0525 TBrRefVec score_trks_br(fNTrks);
0526
0527
0528 FltVec nLayers_ref_trks(fNTrks);
0529 TBrRefVec nLayers_ref_trks_br(fNTrks);
0530 FltVec pt_ref_trks(fNTrks);
0531 TBrRefVec pt_ref_trks_br(fNTrks);
0532 FltVec eta_ref_trks(fNTrks);
0533 TBrRefVec eta_ref_trks_br(fNTrks);
0534 FltVec dphi_trks(fNTrks);
0535 TBrRefVec dphi_trks_br(fNTrks);
0536
0537
0538 for (auto j = 0U; j < fNTrks; j++)
0539 {
0540 const auto& trk = fTrks[j];
0541 auto& refmask_trk = refmask_trks[j];
0542 auto& refmask_trk_br = refmask_trks_br[j];
0543 auto& iTkMatches_trk = iTkMatches_trks[j];
0544 auto& iTkMatches_trk_br = iTkMatches_trks_br[j];
0545 auto& nHits_trk = nHits_trks[j];
0546 auto& nHits_trk_br = nHits_trks_br[j];
0547 auto& fracHits_trk = fracHits_trks[j];
0548 auto& fracHits_trk_br = fracHits_trks_br[j];
0549 auto& score_trk = score_trks[j];
0550 auto& score_trk_br = score_trks_br[j];
0551 auto& nLayers_ref_trk = nLayers_ref_trks[j];
0552 auto& nLayers_ref_trk_br = nLayers_ref_trks_br[j];
0553 auto& pt_ref_trk = pt_ref_trks[j];
0554 auto& pt_ref_trk_br = pt_ref_trks_br[j];
0555 auto& eta_ref_trk = eta_ref_trks[j];
0556 auto& eta_ref_trk_br = eta_ref_trks_br[j];
0557 auto& dphi_trk = dphi_trks[j];
0558 auto& dphi_trk_br = dphi_trks_br[j];
0559
0560
0561 refmask_trk = 0;
0562 refmask_trk_br = 0;
0563 iTkMatches_trk = 0;
0564 iTkMatches_trk_br = 0;
0565
0566
0567 nHits_trk = 0;
0568 nHits_trk_br = 0;
0569 fracHits_trk = 0.f;
0570 fracHits_trk_br = 0;
0571 score_trk = 0;
0572 score_trk_br = 0;
0573
0574
0575 nLayers_ref_trk = 0;
0576 nLayers_ref_trk_br = 0;
0577 pt_ref_trk = 0.f;
0578 pt_ref_trk_br = 0;
0579 eta_ref_trk = 0.f;
0580 eta_ref_trk_br = 0;
0581 dphi_trk = 0.f;
0582 dphi_trk_br = 0;
0583
0584
0585 frtree->SetBranchAddress(fSRefMask + "mask_" + trk, &refmask_trk, &refmask_trk_br);
0586 frtree->SetBranchAddress("iTkMatches_" + trk, &iTkMatches_trk, &iTkMatches_trk_br);
0587
0588 frtree->SetBranchAddress("nHits_" + trk, &nHits_trk, &nHits_trk_br);
0589 frtree->SetBranchAddress("fracHitsMatched_" + trk, &fracHits_trk, &fracHits_trk_br);
0590 frtree->SetBranchAddress("score_" + trk, &score_trk, &score_trk_br);
0591
0592 frtree->SetBranchAddress("nLayers_" + fSRefVarTrk + "_" + trk, &nLayers_ref_trk, &nLayers_ref_trk_br);
0593 frtree->SetBranchAddress("pt_" + fSRefVarTrk + "_" + trk, &pt_ref_trk, &pt_ref_trk_br);
0594 frtree->SetBranchAddress("eta_" + fSRefVarTrk + "_" + trk, &eta_ref_trk, &eta_ref_trk_br);
0595 frtree->SetBranchAddress("dphi_" + trk, &dphi_trk, &dphi_trk_br);
0596 }
0597
0598
0599
0600
0601
0602
0603 const UInt_t nentries = frtree->GetEntries();
0604 for (auto e = 0U; e < nentries; e++) {
0605
0606 for (auto i = 0U; i < fNVars; i++) {
0607 auto& var_trks_br = vars_trks_br[i];
0608 for (auto j = 0U; j < fNTrks; j++) {
0609 auto& var_trk_br = var_trks_br[j];
0610
0611 var_trk_br->GetEntry(e);
0612 }
0613 }
0614 for (auto j = 0U; j < fNTrks; j++) {
0615 auto& refmask_trk_br = refmask_trks_br[j];
0616 auto& iTkMatches_trk_br = iTkMatches_trks_br[j];
0617 auto& nHits_trk_br = nHits_trks_br[j];
0618 auto& fracHits_trk_br = fracHits_trks_br[j];
0619 auto& score_trk_br = score_trks_br[j];
0620 auto& nLayers_ref_trk_br = nLayers_ref_trks_br[j];
0621 auto& pt_ref_trk_br = pt_ref_trks_br[j];
0622 auto& eta_ref_trk_br = eta_ref_trks_br[j];
0623 auto& dphi_trk_br = dphi_trks_br[j];
0624
0625 refmask_trk_br->GetEntry(e);
0626 iTkMatches_trk_br->GetEntry(e);
0627
0628 nHits_trk_br->GetEntry(e);
0629 fracHits_trk_br->GetEntry(e);
0630 score_trk_br->GetEntry(e);
0631
0632 nLayers_ref_trk_br->GetEntry(e);
0633 pt_ref_trk_br->GetEntry(e);
0634 eta_ref_trk_br->GetEntry(e);
0635 dphi_trk_br->GetEntry(e);
0636 }
0637
0638
0639 for (auto j = 0U; j < fNTrks; j++)
0640 {
0641 const auto pt_trk = vars_trks[0][j];
0642 const auto eta_trk = vars_trks[1][j];
0643 const auto phi_trk = vars_trks[2][j];
0644
0645 const auto refmask_trk = refmask_trks[j];
0646 const auto iTkMatches_trk = iTkMatches_trks[j];
0647 const auto nHits_trk = nHits_trks[j];
0648 const auto fracHits_trk = fracHits_trks[j];
0649 const auto score_trk = score_trks[j];
0650
0651 const auto nLayers_ref_trk = nLayers_ref_trks[j];
0652 const auto pt_ref_trk = pt_ref_trks[j];
0653 const auto eta_ref_trk = eta_ref_trks[j];
0654 const auto dphi_trk = dphi_trks[j];
0655
0656 for (auto k = 0U; k < fNPtCuts; k++)
0657 {
0658 const auto ptcut = fPtCuts[k];
0659
0660 if (pt_trk < ptcut)
0661 continue;
0662
0663
0664 for (auto i = 0U; i < fNVars; i++)
0665 {
0666 const auto var_trk = vars_trks[i][j];
0667
0668
0669 const TString plotkey = Form("%i_%i_%i", i, j, k);
0670
0671
0672 if (refmask_trk >= 0)
0673 plots[plotkey]->Fill((refmask_trk == 0), var_trk);
0674 }
0675
0676
0677 const TString basekey = Form("%i_%i", j, k);
0678
0679
0680 const TString nhitkey = Form("%s_0", basekey.Data());
0681 const TString frackey = Form("%s_1", basekey.Data());
0682 const TString scorekey = Form("%s_2", basekey.Data());
0683
0684 const TString dnhitkey = Form("%s_d_0", basekey.Data());
0685 const TString dinvptkey = Form("%s_d_1", basekey.Data());
0686 const TString detakey = Form("%s_d_2", basekey.Data());
0687 const TString dphikey = Form("%s_d_3", basekey.Data());
0688
0689
0690 hists[Form("%s_0", nhitkey.Data())]->Fill(nHits_trk);
0691 hists[Form("%s_0", frackey.Data())]->Fill(fracHits_trk);
0692 hists[Form("%s_0", scorekey.Data())]->Fill(score_trk);
0693
0694 if (refmask_trk == 0)
0695 {
0696 hists[Form("%s_1", nhitkey.Data())]->Fill(nHits_trk);
0697 hists[Form("%s_1", frackey.Data())]->Fill(fracHits_trk);
0698 hists[Form("%s_1", scorekey.Data())]->Fill(score_trk);
0699 } else if (refmask_trk == 1)
0700 {
0701 hists[Form("%s_2", nhitkey.Data())]->Fill(nHits_trk);
0702 hists[Form("%s_2", frackey.Data())]->Fill(fracHits_trk);
0703 hists[Form("%s_2", scorekey.Data())]->Fill(score_trk);
0704
0705 hists[Form("%s_2", dnhitkey.Data())]->Fill(nHits_trk - (Int_t)nLayers_ref_trk);
0706 hists[Form("%s_2", dinvptkey.Data())]->Fill(1.f / pt_trk - 1.f / pt_ref_trk);
0707 hists[Form("%s_2", detakey.Data())]->Fill(eta_trk - eta_ref_trk);
0708 hists[Form("%s_2", dphikey.Data())]->Fill(dphi_trk);
0709
0710 if (iTkMatches_trk == 0)
0711 {
0712 hists[Form("%s_3", nhitkey.Data())]->Fill(nHits_trk);
0713 hists[Form("%s_3", frackey.Data())]->Fill(fracHits_trk);
0714 hists[Form("%s_3", scorekey.Data())]->Fill(score_trk);
0715
0716 hists[Form("%s_3", dnhitkey.Data())]->Fill(nHits_trk - (Int_t)nLayers_ref_trk);
0717 hists[Form("%s_3", dinvptkey.Data())]->Fill(1.f / pt_trk - 1.f / pt_ref_trk);
0718 hists[Form("%s_3", detakey.Data())]->Fill(eta_trk - eta_ref_trk);
0719 hists[Form("%s_3", dphikey.Data())]->Fill(dphi_trk);
0720 }
0721 }
0722 }
0723 }
0724 }
0725
0726
0727
0728
0729
0730
0731 TStrVec dirnames = {"fakerate", "quality", "kindiffs"};
0732 for (auto& dirname : dirnames)
0733 dirname += fSRefDir;
0734 const UInt_t ndirs = dirnames.size();
0735
0736 TDirRefVec subdirs(ndirs);
0737 for (auto q = 0U; q < ndirs; q++)
0738 subdirs[q] = PlotValidation::MakeSubDirs(dirnames[q]);
0739
0740
0741 for (auto j = 0U; j < fNTrks; j++)
0742 {
0743 for (auto k = 0U; k < fNPtCuts; k++)
0744 {
0745
0746 for (auto i = 0U; i < fNVars; i++)
0747 {
0748 const Int_t diridx = 0;
0749 const TString plotkey = Form("%i_%i_%i", i, j, k);
0750 PlotValidation::DrawWriteSavePlot(plots[plotkey], subdirs[diridx], dirnames[diridx], "AP");
0751 delete plots[plotkey];
0752 }
0753
0754
0755 for (auto n = 0U; n < fNTrkQual; n++)
0756 {
0757 for (auto o = 0U; o < ncolls; o++)
0758 {
0759 const Int_t diridx = 1;
0760 const TString histkey = Form("%i_%i_%i_%i", j, k, n, o);
0761 PlotValidation::DrawWriteSavePlot(hists[histkey], subdirs[diridx], dirnames[diridx], "");
0762 delete hists[histkey];
0763 }
0764 }
0765
0766
0767 for (auto p = 0U; p < ndvars; p++)
0768 {
0769 for (auto o = 2U; o < ncolls; o++)
0770 {
0771 const Int_t diridx = 2;
0772 const TString histkey = Form("%i_%i_d_%i_%i", j, k, p, o);
0773 PlotValidation::DrawWriteSavePlot(hists[histkey], subdirs[diridx], dirnames[diridx], "");
0774 delete hists[histkey];
0775 }
0776 }
0777
0778 }
0779 }
0780 }
0781
0782 void PlotValidation::PrintTotals(int algo) {
0783
0784
0785
0786
0787 Int_t Nevents = 0;
0788 Int_t evtID = 0;
0789 TBranch* b_evtID = 0;
0790 efftree->SetBranchAddress("evtID", &evtID, &b_evtID);
0791 const UInt_t nentries = efftree->GetEntries();
0792 for (auto e = 0U; e < nentries; e++) {
0793 b_evtID->GetEntry(e);
0794 if (evtID > Nevents)
0795 Nevents = evtID;
0796 }
0797
0798 const Int_t NtracksMC = efftree->GetEntries();
0799 const Float_t ntkspevMC = Float_t(NtracksMC) / Float_t(Nevents);
0800 const Int_t NtracksReco = frtree->GetEntries();
0801 const Float_t ntkspevReco = Float_t(NtracksReco) / Float_t(Nevents);
0802
0803
0804
0805
0806
0807
0808 const TStrVec rates = {"eff", "fr", "dr"};
0809 const TStrVec srates = {"Efficiency", "Fake Rate", "Duplicate Rate"};
0810 const TStrVec dirnames = {"efficiency", "fakerate", "duplicaterate"};
0811 const TStrVec types = (fCmsswComp ? TStrVec{"cmssw", "reco", "cmssw"}
0812 : TStrVec{"sim", "reco", "sim"});
0813 const UInt_t nrates = rates.size();
0814
0815 const TStrVec snumers = {
0816 fSRefTitle + " Tracks Matched", "Unmatched Reco Tracks", fSRefTitle + " Tracks Matched (nTimes>1)"};
0817 const TStrVec sdenoms = {
0818 "Eligible " + fSRefTitle + " Tracks", "Eligible Reco Tracks", "Eligible " + fSRefTitle + " Tracks"};
0819
0820 TEffRefMap plots;
0821 for (auto j = 0U; j < fNTrks; j++) {
0822 const auto& trk = fTrks[j];
0823
0824 for (auto k = 0U; k < fNPtCuts; k++) {
0825 const auto& hptcut = fHPtCuts[k];
0826
0827 for (auto l = 0U; l < nrates; l++) {
0828 const auto& rate = rates[l];
0829 const auto& type = types[l];
0830 const auto& dirname = dirnames[l];
0831
0832 const TString plotkey = Form("%i_%i_%i", j, k, l);
0833 const TString plotname = dirname + fSRefDir + "/" + rate + "_" + type + "_phi_" + trk + "_pt" + hptcut;
0834 plots[plotkey] = (TEfficiency*)fOutRoot->Get(plotname.Data());
0835 }
0836 }
0837 }
0838
0839
0840 const TStrVec colls = {"allreco", "fake", "bestmatch"};
0841 const TStrVec scolls = {"All Reco", "Fake", "Best Match"};
0842 const UInt_t ncolls = colls.size();
0843
0844 TH1FRefMap hists;
0845 for (auto j = 0U; j < fNTrks; j++) {
0846 const auto& trk = fTrks[j];
0847
0848 for (auto k = 0U; k < fNPtCuts; k++) {
0849 const auto& hptcut = fHPtCuts[k];
0850
0851 for (auto n = 0U; n < fNTrkQual; n++) {
0852 const auto& trkqual = fTrkQual[n];
0853
0854 for (auto o = 0U; o < ncolls; o++) {
0855 const auto& coll = colls[o];
0856
0857 const TString histkey = Form("%i_%i_%i_%i", j, k, n, o);
0858 const TString histname = "quality" + fSRefDir + "/h_" + trkqual + "_" + coll + "_" + trk + "_pt" + hptcut;
0859 hists[histkey] = (TH1F*)fOutRoot->Get(histname.Data());
0860 }
0861 }
0862 }
0863 }
0864
0865
0866 const TString outfilename = fOutName + "/totals_" + fOutName + fSRefOut + ".txt";
0867 std::ofstream totalsout(outfilename.Data());
0868
0869 std::cout << "--------Track Reconstruction Summary--------" << std::endl;
0870 std::cout << "nEvents: " << Nevents << Form(" n%sTracks/evt: ", fSRefTitle.Data()) << ntkspevMC
0871 << " nRecoTracks/evt: " << ntkspevReco << std::endl;
0872 std::cout << "++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
0873 std::cout << std::endl;
0874
0875 totalsout << "--------Track Reconstruction Summary--------" << std::endl;
0876 totalsout << "nEvents: " << Nevents << Form(" n%sTracks/evt: ", fSRefTitle.Data()) << ntkspevMC
0877 << " nRecoTracks/evt: " << ntkspevReco << std::endl;
0878 totalsout << "++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
0879 totalsout << std::endl;
0880
0881 for (auto k = 0U; k < fNPtCuts; k++) {
0882 const auto& ptcut = fPtCuts[k];
0883
0884 std::cout << Form("xxxxxxxxxx Track pT > %3.1f Cut xxxxxxxxxx", ptcut) << std::endl;
0885 std::cout << std::endl;
0886
0887 totalsout << Form("xxxxxxxxxx Track pT > %3.1f Cut xxxxxxxxxx", ptcut) << std::endl;
0888 totalsout << std::endl;
0889
0890 for (auto j = 0U; j < fNTrks; j++) {
0891 const auto& strk = fSTrks[j];
0892
0893 std::cout << strk.Data() << " Tracks" << std::endl;
0894 std::cout << "++++++++++++++++++++++++++++++++++++++++++" << std::endl << std::endl;
0895 std::cout << "Quality Info for " << strk.Data() << " Track Collections" << std::endl;
0896 std::cout << "==========================================" << std::endl;
0897
0898 totalsout << strk.Data() << " Tracks" << std::endl;
0899 totalsout << "++++++++++++++++++++++++++++++++++++++++++" << std::endl << std::endl;
0900 totalsout << "Quality Info for " << strk.Data() << " Track Collections" << std::endl;
0901 totalsout << "==========================================" << std::endl;
0902 for (auto o = 0U; o < ncolls; o++) {
0903 const auto& scoll = scolls[o];
0904
0905 const Float_t nHits_mean = hists[Form("%i_%i_0_%i", j, k, o)]->GetMean(1);
0906 const Float_t nHits_mean_unc = hists[Form("%i_%i_0_%i", j, k, o)]->GetMeanError(1);
0907 const Float_t fracHits_mean = hists[Form("%i_%i_1_%i", j, k, o)]->GetMean(1);
0908 const Float_t fracHits_mean_unc = hists[Form("%i_%i_1_%i", j, k, o)]->GetMeanError(1);
0909 const Float_t score_mean = hists[Form("%i_%i_2_%i", j, k, o)]->GetMean(1);
0910 const Float_t score_mean_unc = hists[Form("%i_%i_2_%i", j, k, o)]->GetMeanError(1);
0911
0912 std::cout << scoll.Data() << " Tracks" << std::endl;
0913 std::cout << "Mean nHits / Track = " << nHits_mean << " +/- " << nHits_mean_unc << std::endl;
0914 std::cout << "Mean Shared Hits / Track = " << fracHits_mean << " +/- " << fracHits_mean_unc << std::endl;
0915 std::cout << "Mean Track Score = " << score_mean << " +/- " << score_mean_unc << std::endl;
0916 std::cout << "------------------------------------------" << std::endl;
0917
0918 totalsout << scoll.Data() << " Tracks" << std::endl;
0919 totalsout << "Mean nHits / Track = " << nHits_mean << " +/- " << nHits_mean_unc << std::endl;
0920 totalsout << "Mean Shared Hits / Track = " << fracHits_mean << " +/- " << fracHits_mean_unc << std::endl;
0921 totalsout << "Mean Track Score = " << score_mean << " +/- " << score_mean_unc << std::endl;
0922 totalsout << "------------------------------------------" << std::endl;
0923 }
0924
0925 std::cout << std::endl << "Rates for " << strk.Data() << " Tracks" << std::endl;
0926 std::cout << "==========================================" << std::endl;
0927
0928 totalsout << std::endl << "Rates for " << strk.Data() << " Tracks" << std::endl;
0929 totalsout << "==========================================" << std::endl;
0930 for (auto l = 0U; l < nrates; l++) {
0931 const auto& snumer = snumers[l];
0932 const auto& sdenom = sdenoms[l];
0933 const auto& srate = srates[l];
0934
0935 EffStruct effs;
0936 PlotValidation::GetTotalEfficiency(plots[Form("%i_%i_%i", j, k, l)], effs);
0937
0938 std::cout << snumer.Data() << ": " << effs.passed_ << std::endl;
0939 std::cout << sdenom.Data() << ": " << effs.total_ << std::endl;
0940 std::cout << "------------------------------------------" << std::endl;
0941 std::cout << srate.Data() << ": " << effs.eff_ << ", -" << effs.elow_ << ", +" << effs.eup_ << std::endl;
0942 std::cout << "------------------------------------------" << std::endl;
0943
0944 totalsout << snumer.Data() << ": " << effs.passed_ << std::endl;
0945 totalsout << sdenom.Data() << ": " << effs.total_ << std::endl;
0946 totalsout << "------------------------------------------" << std::endl;
0947 totalsout << srate.Data() << ": " << effs.eff_ << ", -" << effs.elow_ << ", +" << effs.eup_ << std::endl;
0948 totalsout << "------------------------------------------" << std::endl;
0949 }
0950 std::cout << std::endl << std::endl;
0951 totalsout << std::endl << std::endl;
0952 }
0953 }
0954
0955
0956 for (auto& hist : hists)
0957 delete hist.second;
0958 for (auto& plot : plots)
0959 delete plot.second;
0960 }
0961
0962 template <typename T>
0963 void PlotValidation::DrawWriteSavePlot(T*& plot, TDirectory*& subdir, const TString& subdirname, const TString& option) {
0964
0965 subdir->cd();
0966 plot->SetDirectory(subdir);
0967 plot->Write(plot->GetName(), TObject::kWriteDelete);
0968
0969
0970 if (fSaveAs) {
0971 auto canv = new TCanvas();
0972 canv->cd();
0973 plot->Draw(option.Data());
0974
0975
0976 canv->SetLogy(1);
0977 canv->SaveAs(Form("%s/%s/log/%s.%s", fOutName.Data(), subdirname.Data(), plot->GetName(), fOutType.Data()));
0978
0979
0980 canv->SetLogy(0);
0981 canv->SaveAs(Form("%s/%s/lin/%s.%s", fOutName.Data(), subdirname.Data(), plot->GetName(), fOutType.Data()));
0982
0983 delete canv;
0984 }
0985 }
0986
0987 void PlotValidation::GetTotalEfficiency(const TEfficiency* eff, EffStruct& effs) {
0988 effs.passed_ = eff->GetPassedHistogram()->Integral();
0989 effs.total_ = eff->GetTotalHistogram()->Integral();
0990
0991 auto tmp_eff = new TEfficiency("tmp_eff", "tmp_eff", 1, 0, 1);
0992 tmp_eff->SetTotalEvents(1, effs.total_);
0993 tmp_eff->SetPassedEvents(1, effs.passed_);
0994
0995 effs.eff_ = tmp_eff->GetEfficiency(1);
0996 effs.elow_ = tmp_eff->GetEfficiencyErrorLow(1);
0997 effs.eup_ = tmp_eff->GetEfficiencyErrorUp(1);
0998
0999 delete tmp_eff;
1000 }
1001
1002 void PlotValidation::MakeOutDir(const TString& outdirname) {
1003
1004 FileStat_t dummyFileStat;
1005 if (gSystem->GetPathInfo(outdirname.Data(), dummyFileStat) == 1) {
1006 const TString mkDir = "mkdir -p " + outdirname;
1007 gSystem->Exec(mkDir.Data());
1008 }
1009 }
1010
1011 void PlotValidation::MoveInput() {
1012 const TString mvin = "mv " + fInName + " " + fOutName;
1013 gSystem->Exec(mvin.Data());
1014 }
1015
1016 TDirectory* PlotValidation::MakeSubDirs(const TString& subdirname) {
1017 PlotValidation::MakeOutDir(fOutName + "/" + subdirname);
1018 PlotValidation::MakeOutDir(fOutName + "/" + subdirname + "/lin");
1019 PlotValidation::MakeOutDir(fOutName + "/" + subdirname + "/log");
1020
1021 return fOutRoot->mkdir(subdirname.Data());
1022 }
1023
1024 void PlotValidation::SetupStyle() {
1025
1026 gROOT->Reset();
1027 gStyle->SetOptStat("emou");
1028 gStyle->SetTitleFontSize(0.04);
1029 gStyle->SetOptFit(1011);
1030 gStyle->SetStatX(0.9);
1031 gStyle->SetStatW(0.1);
1032 gStyle->SetStatY(1.0);
1033 gStyle->SetStatH(0.08);
1034 }
1035
1036 void PlotValidation::SetupBins() {
1037
1038 PlotValidation::SetupVariableBins(
1039 "0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.5 3 3.5 4 4.5 5 5 6 7 8 9 10 15 20 25 30 40 50 100 200 500 1000", fPtBins);
1040
1041
1042 PlotValidation::SetupFixedBins(60, -3, 3, fEtaBins);
1043
1044
1045 PlotValidation::SetupFixedBins(70, -3.5, 3.5, fPhiBins);
1046
1047
1048 PlotValidation::SetupFixedBins(26, -0.5, 25.5, fNLayersBins);
1049
1050
1051 PlotValidation::SetupFixedBins(40, 0, 40, fNHitsBins);
1052
1053
1054 PlotValidation::SetupFixedBins(110, 0, 1.1, fFracHitsBins);
1055
1056
1057 PlotValidation::SetupFixedBins(50, -500, 5000, fScoreBins);
1058
1059
1060 PlotValidation::SetupFixedBins(40, -20, 20, fDNHitsBins);
1061
1062
1063 PlotValidation::SetupFixedBins(45, -1.0, 1.0, fDInvPtBins);
1064
1065
1066 PlotValidation::SetupFixedBins(45, -0.1, 0.1, fDPhiBins);
1067
1068
1069 PlotValidation::SetupFixedBins(45, -0.1, 0.1, fDEtaBins);
1070 }
1071
1072 void PlotValidation::SetupVariableBins(const std::string& s_bins, DblVec& bins) {
1073 std::stringstream ss(s_bins);
1074 Double_t boundary;
1075 while (ss >> boundary)
1076 bins.emplace_back(boundary);
1077 }
1078
1079 void PlotValidation::SetupFixedBins(const UInt_t nBins, const Double_t low, const Double_t high, DblVec& bins) {
1080 const Double_t width = (high - low) / nBins;
1081
1082 for (auto i = 0U; i <= nBins; i++)
1083 bins.emplace_back(i * width + low);
1084 }
1085
1086 void PlotValidation::SetupCommonVars() {
1087
1088 fVars = {"pt", "eta", "phi", "nLayers"};
1089 fSVars = {"p_{T}", "#eta", "#phi", "Number of layers"};
1090 fSUnits = {"GeV/c", "", "", ""};
1091 fNVars = fVars.size();
1092
1093 fSVarPt = fSVars[0];
1094 fSUnitPt = fSUnits[0];
1095
1096
1097 for (auto& sunit : fSUnits) {
1098 if (!sunit.EqualTo("")) {
1099 sunit.Prepend(" [");
1100 sunit.Append("]");
1101 }
1102 }
1103
1104
1105 fVarBins = {fPtBins, fEtaBins, fPhiBins, fNLayersBins};
1106
1107
1108 fTrks = (fCmsswComp ? TStrVec{"build", "fit"} : TStrVec{"seed", "build", "fit"});
1109 fSTrks = (fCmsswComp ? TStrVec{"Build", "Fit"}
1110 : TStrVec{"Seed", "Build", "Fit"});
1111 fNTrks = fTrks.size();
1112
1113
1114 fPtCuts = {0.f, 0.9f, 2.f};
1115 for (const auto ptcut : fPtCuts) {
1116 fSPtCuts.emplace_back(Form("%3.1f", ptcut));
1117 }
1118 for (const auto& sptcut : fSPtCuts) {
1119 TString hptcut = sptcut;
1120 hptcut.ReplaceAll(".", "p");
1121 fHPtCuts.emplace_back(hptcut);
1122 }
1123 fNPtCuts = fPtCuts.size();
1124
1125
1126 fTrkQual = {"nHits", "fracHitsMatched", "score"};
1127 fSTrkQual = {"nHits / Track", "Highest Fraction of Matched Hits / Track", "Track Score"};
1128 fNTrkQual = fTrkQual.size();
1129
1130
1131 fSRefTitle = (fCmsswComp ? "CMSSW" : "Sim");
1132 fSRefVar = (fCmsswComp ? "cmssw" : "mc_gen");
1133 fSRefMask = (fCmsswComp ? "cmssw" : "mc");
1134 fSRefVarTrk = (fCmsswComp ? "cmssw" : "mc");
1135 fSRefDir = (fCmsswComp ? "_cmssw" : "");
1136 fSRefOut = (fCmsswComp ? "_cmssw" : "");
1137 }