Scripts/TrackingPerformance/defineReconstructionPerformance.C

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2021 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008
0009 #include <array>
0010 #include <iostream>
0011 #include <string>
0012 #include <tuple>
0013 #include <vector>
0014
0015 #include <TCanvas.h>
0016 #include <TEfficiency.h>
0017 #include <TFile.h>
0018
0019 #include "CommonUtils.h"
0020 #include "TreeReader.h"
0021
0022 /// This script/function reads all the reconstructed tracks from e.g.
0023 /// 'tracksummary_ckf.root' and the (possibly selected) truth particles from
0024 /// e.g. 'track_finder_particles.root' (which contains the info of 'nHits'), and
0025 /// defines the efficiency, fake rate and duplicaiton rate. It aims to make
0026 /// custom definition and tuning of the reconstruction performance easier.
0027 /// Multiple files for the reconstructed tracks are allowed.
0028 ///
0029 /// NB: It's very likely that fiducal cuts are already imposed on the truth
0030 /// particles. Please check the selection criteria in the truth fitting example
0031 /// which writes out the 'track_finder_particles.root'. For instance, if the
0032 /// truth particles are already required to have pT > 1 GeV, it does not make
0033 /// sense to have ptMin = 0.5 GeV here.
0034 ///
0035 void defineReconstructionPerformance(
0036     const std::string& inputSimParticleFileName =
0037         "performance_track_finder.root",
0038     const std::vector<std::string>& inputTrackSummaryFileNames =
0039         {"tracksummary_ckf.root"},
0040     const std::vector<std::string>& trackSummaryFileLegends =
0041         {"CKF with reco seeds"},
0042     const std::string& simParticleTreeName = "track_finder_particles",
0043     const std::string& trackSummaryTreeName = "tracksummary_ckf",
0044     unsigned int nHitsMin = 9, unsigned int nMeasurementsMin = 6,
0045     double ptMin = 0.5, double truthMatchProbMin = 0.5) {
0046   gStyle->SetOptFit(0011);
0047   gStyle->SetOptStat(0000);
0048   gStyle->SetPadLeftMargin(0.18);
0049   gStyle->SetPadRightMargin(0.05);
0050   gStyle->SetPadTopMargin(0.1);
0051   gStyle->SetPadBottomMargin(0.15);
0052   gStyle->SetTitleSize(0.05, "xy");
0053   gStyle->SetLabelSize(0.05, "xy");
0054   gStyle->SetTitleOffset(1., "x");
0055   gStyle->SetTitleOffset(1.5, "y");
0056   gStyle->SetNdivisions(505, "y");
0057
0058   // Check the inputs are valid
0059   if (inputTrackSummaryFileNames.size() != trackSummaryFileLegends.size()) {
0060     throw std::invalid_argument(
0061         "Please specify the legends you want to show for all the track files");
0062   }
0063
0064   // Open the files for reading
0065   TFile* particleFile = TFile::Open(inputSimParticleFileName.c_str(), "read");
0066   // The number of track files to read
0067   unsigned int nTrackFiles = inputTrackSummaryFileNames.size();
0068   std::vector<TFile*> trackFiles;
0069   trackFiles.reserve(nTrackFiles);
0070   for (const auto& fileName : inputTrackSummaryFileNames) {
0071     trackFiles.push_back(TFile::Open(fileName.c_str(), "read"));
0072   }
0073
0074   // Define variables for tree reading (turn on the events sorting since we have more than one root files to read)
0075   ParticleReader pReader(
0076       (TTree*)particleFile->Get(simParticleTreeName.c_str()), true);
0077   std::vector<TrackSummaryReader> tReaders;
0078   tReaders.reserve(nTrackFiles);
0079   for (const auto& trackFile : trackFiles) {
0080     tReaders.emplace_back((TTree*)trackFile->Get(trackSummaryTreeName.c_str()), true);
0081   }
0082
0083   std::vector<std::size_t> nEvents;
0084   nEvents.reserve(nTrackFiles);
0085   for (const auto& tReader : tReaders) {
0086     std::size_t entries = tReader.tree->GetEntries();
0087     nEvents.push_back(entries);
0088   }
0089
0090   // Define the efficiency plots
0091   std::vector<TEfficiency*> trackEff_vs_eta;
0092   std::vector<TEfficiency*> fakeRate_vs_eta;
0093   std::vector<TEfficiency*> duplicateRate_vs_eta;
0094   std::vector<TEfficiency*> trackEff_vs_pt;
0095   std::vector<TEfficiency*> fakeRate_vs_pt;
0096   std::vector<TEfficiency*> duplicateRate_vs_pt;
0097
0098   for (int i = 0; i < nTrackFiles; ++i) {
0099     trackEff_vs_eta.push_back(new TEfficiency(
0100         Form("trackeff_vs_eta_%i", i), ";Truth #eta [GeV/c];Efficiency", 40, -4, 4));
0101     fakeRate_vs_eta.push_back(new TEfficiency(
0102         Form("fakerate_vs_eta_%i", i), ";#eta [GeV/c];fake rate", 40, -4, 4));
0103     duplicateRate_vs_eta.push_back(new TEfficiency(
0104         Form("duplicaterate_vs_eta_%i", i), ";#eta [GeV/c];Duplicate rate", 40, -4, 4));
0105     trackEff_vs_pt.push_back(new TEfficiency(
0106         Form("trackeff_vs_pt_%i", i), ";Truth pt [GeV/c];Efficiency", 40, 0, 100));
0107     fakeRate_vs_pt.push_back(new TEfficiency(
0108         Form("fakerate_vs_pt_%i", i), ";pt [GeV/c];fake rate", 40, 0, 100));
0109     duplicateRate_vs_pt.push_back(new TEfficiency(
0110         Form("duplicaterate_vs_pt_%i", i), ";pt [GeV/c];Duplicate rate", 40, 0, 100));
0111   }
0112
0113   // Set styles
0114   for (int i = 0; i < nTrackFiles; ++i) {
0115     auto color = i + 1;
0116     setEffStyle(trackEff_vs_eta[i], color);
0117     setEffStyle(fakeRate_vs_eta[i], color);
0118     setEffStyle(duplicateRate_vs_eta[i], color);
0119     setEffStyle(trackEff_vs_pt[i], color);
0120     setEffStyle(fakeRate_vs_pt[i], color);
0121     setEffStyle(duplicateRate_vs_pt[i], color);
0122   }
0123
0124   // The particles in each event
0125   std::map<unsigned int, std::vector<ParticleInfo>> particles;
0126
0127   // Loop over the track files
0128   for (unsigned int ifile = 0; ifile < nTrackFiles; ++ifile) {
0129     std::cout << "Processing track file: " << inputTrackSummaryFileNames[ifile]
0130               << std::endl;
0131
0132     // The container from track-particle matching info (Flushed per event)
0133     std::map<uint64_t, std::vector<RecoTrackInfo>> matchedParticles;
0134
0135     // Loop over the events to fill plots
0136     for (std::size_t i = 0; i < nEvents[ifile]; ++i) {
0137       if (i % 10 == 0) {
0138         std::cout << "Processed events: " << i << std::endl;
0139       }
0140
0141       // Get the tracks
0142       tReaders[ifile].getEntry(i);
0143
0144       // Get the particles (do nothing if the particles for this event already
0145       // read)
0146       auto it = particles.find(i);
0147       if (it == particles.end()) {
0148         particles.emplace(i, pReader.getParticles(i));
0149       }
0150
0151       // Loop over the tracks
0152       // The fake rate is defined as the ratio of selected truth-matched tracks
0153       // over all selected tracks
0154       for (std::size_t j = 0; j < tReaders[ifile].nStates->size(); ++j) {
0155         bool hasFittedParameters = tReaders[ifile].hasFittedParams->at(j);
0156         auto nMeasurements = tReaders[ifile].nMeasurements->at(j);
0157         auto nOutliers = tReaders[ifile].nOutliers->at(j);
0158         auto nHoles = tReaders[ifile].nHoles->at(j);
0159         auto theta = tReaders[ifile].eTHETA_fit->at(j);
0160         auto qop = tReaders[ifile].eQOP_fit->at(j);
0161         auto pt = std::abs(1 / qop) * std::sin(theta);
0162         auto eta = std::atanh(std::cos(theta));
0163         auto nMajorityHits = tReaders[ifile].nMajorityHits->at(j);
0164         auto majorityParticleId = tReaders[ifile].majorityParticleId->at(j);
0165
0166         // Select the track, e.g. you might also want to add cuts on the
0167         // nOutliers, nHoles
0168         if ((!hasFittedParameters) or nMeasurements < nMeasurementsMin or
0169             pt < ptMin) {
0170           continue;
0171         }
0172
0173         // Fill the fake rate plots
0174         if (nMajorityHits * 1. / nMeasurements >= truthMatchProbMin) {
0175           matchedParticles[majorityParticleId].push_back(
0176               {eta, pt, nMajorityHits, nMeasurements});
0177           fakeRate_vs_eta[ifile]->Fill(false, eta);
0178           fakeRate_vs_pt[ifile]->Fill(false, pt);
0179         } else {
0180           fakeRate_vs_eta[ifile]->Fill(true, eta);
0181           fakeRate_vs_pt[ifile]->Fill(true, pt);
0182         }
0183       }  // end of all tracks
0184
0185       // Loop over all selected and truth-matched tracks
0186       // The duplicate rate is defined as the ratio of duplicate tracks among
0187       // all the selected truth-matched tracks (only one track is 'real'; others
0188       // are 'duplicated')
0189       for (auto& [id, matchedTracks] : matchedParticles) {
0190         // Sort all tracks matched to this particle according to majority prob
0191         // and track quality
0192         std::sort(matchedTracks.begin(), matchedTracks.end(),
0193                   [](const RecoTrackInfo& lhs, const RecoTrackInfo& rhs) {
0194                     if (lhs.nMajorityHits > rhs.nMajorityHits) {
0195                       return true;
0196                     }
0197                     if (lhs.nMajorityHits < rhs.nMajorityHits) {
0198                       return false;
0199                     }
0200                     if (lhs.nMeasurements > rhs.nMeasurements) {
0201                       return true;
0202                     }
0203                     return false;
0204                   });
0205         // Fill the duplication rate plots
0206         for (std::size_t k = 0; k < matchedTracks.size(); ++k) {
0207           auto eta = matchedTracks[k].eta;
0208           auto pt = matchedTracks[k].pt;
0209           if (k == 0) {
0210             duplicateRate_vs_eta[ifile]->Fill(false, eta);
0211             duplicateRate_vs_pt[ifile]->Fill(false, pt);
0212           } else {
0213             duplicateRate_vs_eta[ifile]->Fill(true, eta);
0214             duplicateRate_vs_pt[ifile]->Fill(true, pt);
0215           }
0216         }
0217       }  // end of all selected truth-matched tracks
0218
0219       // Loop over all truth particles in this event
0220       // The efficiency is defined as the ratio of selected particles that have
0221       // been matched with reco
0222       for (const auto& particle : particles[i]) {
0223         auto nHits = particle.nHits;
0224         auto eta = particle.eta;
0225         auto pt = particle.pt;
0226         if (nHits < nHitsMin or pt < ptMin) {
0227           continue;
0228         }
0229         uint64_t id = particle.particleId;
0230
0231         // Fill the efficiency plots
0232         auto ip = matchedParticles.find(id);
0233         if (ip != matchedParticles.end()) {
0234           trackEff_vs_eta[ifile]->Fill(true, eta);
0235           trackEff_vs_pt[ifile]->Fill(true, pt);
0236         } else {
0237           trackEff_vs_eta[ifile]->Fill(false, eta);
0238           trackEff_vs_pt[ifile]->Fill(false, pt);
0239         }
0240       }  // end of all particles
0241
0242       matchedParticles.clear();
0243     }  // end of all events
0244   }    // end of all track files
0245
0246   std::cout << "All good. Now plotting..." << std::endl;
0247
0248   // The legends
0249   std::vector<TLegend*> legs;
0250   for (int i = 0; i < 6; ++i) {
0251     TLegend* legend = new TLegend(0.4, 0.8, 0.9, 0.9);
0252     legend->SetBorderSize(0);
0253     legend->SetFillStyle(0);
0254     legend->SetTextFont(42);
0255     legs.push_back(legend);
0256   }
0257   // Add entry for the legends
0258   for (int i = 0; i < nTrackFiles; ++i) {
0259     legs[0]->AddEntry(trackEff_vs_eta[i],
0260                       Form("%s", trackSummaryFileLegends[i].c_str()), "lp");
0261     legs[1]->AddEntry(fakeRate_vs_eta[i],
0262                       Form("%s", trackSummaryFileLegends[i].c_str()), "lp");
0263     legs[2]->AddEntry(duplicateRate_vs_eta[i],
0264                       Form("%s", trackSummaryFileLegends[i].c_str()), "lp");
0265     legs[3]->AddEntry(trackEff_vs_pt[i],
0266                       Form("%s", trackSummaryFileLegends[i].c_str()), "lp");
0267     legs[4]->AddEntry(fakeRate_vs_pt[i],
0268                       Form("%s", trackSummaryFileLegends[i].c_str()), "lp");
0269     legs[5]->AddEntry(duplicateRate_vs_pt[i],
0270                       Form("%s", trackSummaryFileLegends[i].c_str()), "lp");
0271   }
0272
0273   // Now draw the plots
0274   TCanvas* c1 = new TCanvas("recoPerf", " ", 1500, 800);
0275   c1->Divide(3, 2);
0276
0277   float scaleRangeMax = 1.15;
0278   for (int i = 0; i < nTrackFiles; ++i) {
0279     std::string mode = (i == 0) ? "" : "same";
0280     c1->cd(1);
0281     trackEff_vs_eta[i]->Draw(mode.c_str());
0282     if (i == nTrackFiles - 1) {
0283       legs[0]->Draw();
0284     }
0285     adaptEffRange(trackEff_vs_eta[i], 1, scaleRangeMax);
0286
0287     c1->cd(2);
0288     fakeRate_vs_eta[i]->Draw(mode.c_str());
0289     if (i == nTrackFiles - 1) {
0290       legs[1]->Draw();
0291     }
0292     adaptEffRange(fakeRate_vs_eta[i], 1, scaleRangeMax);
0293
0294     c1->cd(3);
0295     duplicateRate_vs_eta[i]->Draw(mode.c_str());
0296     if (i == nTrackFiles - 1) {
0297       legs[2]->Draw();
0298     }
0299     adaptEffRange(duplicateRate_vs_eta[i], 1, scaleRangeMax);
0300
0301     c1->cd(4);
0302     trackEff_vs_pt[i]->Draw(mode.c_str());
0303     if (i == nTrackFiles - 1) {
0304       legs[3]->Draw();
0305     }
0306     adaptEffRange(trackEff_vs_pt[i], 1, scaleRangeMax);
0307
0308     c1->cd(5);
0309     fakeRate_vs_pt[i]->Draw(mode.c_str());
0310     if (i == nTrackFiles - 1) {
0311       legs[4]->Draw();
0312     }
0313     adaptEffRange(fakeRate_vs_pt[i], 1, scaleRangeMax);
0314
0315     c1->cd(6);
0316     duplicateRate_vs_pt[i]->Draw(mode.c_str());
0317     if (i == nTrackFiles - 1) {
0318       legs[5]->Draw();
0319     }
0320     adaptEffRange(duplicateRate_vs_pt[i], 1, scaleRangeMax);
0321   }
0322
0323   c1->Update();
0324 }