Scripts/TrackingPerformance/defineReconstructionPerformance.C

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