Scripts/TrackingPerformance/trackSummaryAnalysis.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 <bitset>
0010 #include <fstream>
0011 #include <iostream>
0012 #include <limits>
0013 #include <map>
0014 #include <string>
0015 #include <numbers>
0016 #include <vector>
0017
0018 #include <TCanvas.h>
0019 #include <TChain.h>
0020 #include <TF1.h>
0021 #include <TFile.h>
0022 #include <TH1F.h>
0023 #include <TH2F.h>
0024 #include <TH3F.h>
0025 #include <TLegend.h>
0026 #include <TMath.h>
0027 #include <TStyle.h>
0028 #include <TTree.h>
0029 #include <TVectorF.h>
0030
0031 #include "CommonUtils.h"
0032 #include "TreeReader.h"
0033
0034 using namespace ROOT;
0035
0036 /// This ROOT script will plot the residual and pull of perigee track parameters
0037 /// (d0, z0, phi, theta, q/p, pT, t) from root file produced by the
0038 /// RootTrackSummaryWriter
0039 ///
0040 /// @param inFiles the list of input files
0041 /// @param inTree the name of the input tree
0042 /// @param outFile the name of the output file
0043 /// @param inConfig the (optional) input configuration JSON file
0044 /// @param outConfig the (optional) output configuration JSON file
0045 /// @param residualPulls the bitset of the parameters set
0046 ///
0047 ///
0048 /// @note A lot of this could be done with creating appropriate TH3F histograms
0049 /// and chose the relevant projections, but this would need one loop through the
0050 /// entries per parameter and would create a very long analysis turn-around.
0051 /// That's why the data/histograms are prepared in handles and then filled in a
0052 /// single event analysis loop.
0053 ///
0054 int trackSummaryAnalysis(
0055     const std::vector<std::string>& inFiles, const std::string& inTree,
0056     const std::string& outFile, const std::string& inConfig = "",
0057     const std::string& outConfig = "", unsigned long nEntries = 0,
0058     unsigned int nPeakEntries = 0, float pullRange = 6.,
0059     unsigned int nHistBins = 61, unsigned int nPhiBins = 10,
0060     const std::array<float, 2>& phiRange = {-std::numbers::pi_v<float>, std::numbers::pi_v<float>},
0061     unsigned int nEtaBins = 10, const std::array<float, 2>& etaRange = {-3, 3},
0062     const std::vector<double>& ptBorders =
0063         {0., std::numeric_limits<double>::infinity()},
0064     const std::bitset<7>& residualPulls = std::bitset<7>{"1111111"},
0065     const std::bitset<5>& auxiliary = std::bitset<5>{"11111"}) {
0066   // Load the tree chain
0067   TChain* treeChain = new TChain(inTree.c_str());
0068   for (const auto& inFile : inFiles) {
0069     treeChain->Add(inFile.c_str());
0070     // Open root file written by RootTrackWriter
0071     std::cout << "*** Adding file: " << inFile << std::endl;
0072   }
0073
0074   if (treeChain->GetEntries() == 0) {
0075     std::cout << "xxx No entries found ... bailing out." << std::endl;
0076     return -1;
0077   }
0078
0079   TCanvas* rangeCanvas =
0080       new TCanvas("rangeCanvas", "Range Estimation", 100, 100, 620, 400);
0081
0082   TrackSummaryReader tracks(treeChain, true);
0083
0084   // Loop over the entries of the keys
0085   unsigned long entries = estimateEntries(*tracks.tree, nEntries);
0086   unsigned long peakEntries = estimateEntries(*tracks.tree, nPeakEntries);
0087
0088   // Temporary counter
0089   unsigned int histBarcode = 0;
0090
0091   // Deduction
0092   unsigned int nPtBins = ptBorders.size() - 1;
0093
0094   // One time initialization of the residual and pull handles
0095   using ResidualPulls = std::vector<ResidualPullHandle>;
0096   ResidualPulls baseResidualPulls = {};
0097
0098   if (residualPulls.test(0)) {
0099     // A standard d0Handle
0100     ResidualPullHandle d0Handle;
0101     d0Handle.tag = "d0";
0102     d0Handle.residualStr = "d_{0}^{rec} - d_{0}^{true}";
0103     d0Handle.residualUnit = "[mm]";
0104     d0Handle.errorStr = "#sigma(d_{0})";
0105     d0Handle.rangeDrawStr = "eLOC0_fit-t_d0";
0106     d0Handle.rangeMaxStr = "(1000,-10,10)";
0107     d0Handle.value = ResidualAccessor{tracks.eLOC0_fit, tracks.t_d0};
0108     d0Handle.error = DirectAccessor<float>{tracks.err_eLOC0_fit};
0109     d0Handle.accept = AcceptAll{};
0110     // Push it
0111     baseResidualPulls.push_back(d0Handle);
0112   }
0113
0114   if (residualPulls.test(1)) {
0115     // A standard z0Handle
0116     ResidualPullHandle z0Handle;
0117     z0Handle.tag = "z0";
0118     z0Handle.residualStr = "z_{0}^{rec} - z_{0}^{true}";
0119     z0Handle.residualUnit = "[mm]";
0120     z0Handle.errorStr = "#sigma(z_{0})";
0121     z0Handle.rangeDrawStr = "eLOC1_fit-t_z0";
0122     z0Handle.rangeMaxStr = "(1000,-10,10)";
0123     z0Handle.value = ResidualAccessor{tracks.eLOC1_fit, tracks.t_z0};
0124     z0Handle.error = DirectAccessor<float>{tracks.err_eLOC1_fit};
0125     z0Handle.accept = AcceptAll{};
0126     // Push it
0127     baseResidualPulls.push_back(z0Handle);
0128   }
0129
0130   if (residualPulls.test(2)) {
0131     // A standard phi0Handle
0132     ResidualPullHandle phi0Handle;
0133     phi0Handle.tag = "phi0";
0134     phi0Handle.residualStr = "#phi_{0}^{rec} - #phi_{0}^{true}";
0135     phi0Handle.errorStr = "#sigma(phi_{0})";
0136     phi0Handle.rangeDrawStr = "ePHI_fit-t_phi";
0137     phi0Handle.rangeMaxStr = "(1000,-0.01,0.01)";
0138     phi0Handle.value = ResidualAccessor{tracks.ePHI_fit, tracks.t_phi};
0139     phi0Handle.error = DirectAccessor<float>{tracks.err_ePHI_fit};
0140     phi0Handle.accept = AcceptAll{};
0141     // Push it
0142     baseResidualPulls.push_back(phi0Handle);
0143   }
0144
0145   if (residualPulls.test(3)) {
0146     // A standard theta0Handle
0147     ResidualPullHandle theta0Handle;
0148     theta0Handle.tag = "theta0";
0149     theta0Handle.residualStr = "#theta_{0}^{rec} - #theta_{0}^{true}";
0150     theta0Handle.errorStr = "#sigma(theta_{0})";
0151     theta0Handle.rangeDrawStr = "eTHETA_fit-t_theta";
0152     theta0Handle.rangeMaxStr = "(1000,-0.01,0.01)";
0153     theta0Handle.value = ResidualAccessor{tracks.eTHETA_fit, tracks.t_theta};
0154     theta0Handle.error = DirectAccessor<float>{tracks.err_eTHETA_fit};
0155     theta0Handle.accept = AcceptAll{};
0156     // Push it
0157     baseResidualPulls.push_back(theta0Handle);
0158   }
0159
0160   if (residualPulls.test(4)) {
0161     // The standard qop Handle
0162     ResidualPullHandle qopHandle;
0163     qopHandle.tag = "qop";
0164     qopHandle.residualStr = "q/p^{rec} - q/p^{true}";
0165     qopHandle.residualUnit = "[GeV^{-1}]";
0166     qopHandle.errorStr = "#sigma(q/p)";
0167     qopHandle.rangeDrawStr = "eQOP_fit-t_charge/t_p";
0168     qopHandle.rangeMaxStr = "(1000,-0.1,0.1)";
0169     qopHandle.value =
0170         QopResidualAccessor{tracks.eQOP_fit, tracks.t_charge, tracks.t_p};
0171     qopHandle.error = DirectAccessor<float>{tracks.err_eQOP_fit};
0172     qopHandle.accept = AcceptAll{};
0173     // Push it
0174     baseResidualPulls.push_back(qopHandle);
0175   }
0176
0177   if (residualPulls.test(5)) {
0178     // The pt measurement
0179     ResidualPullHandle tHandle;
0180     tHandle.tag = "t";
0181     tHandle.residualStr = "t^{rec} - t^{true}";
0182     tHandle.residualUnit = "[ns]";
0183     tHandle.errorStr = "#sigma(t})";
0184     tHandle.rangeDrawStr = "eT_fit - t_time";
0185     tHandle.rangeMaxStr = "(1000,-10.,10.)";
0186     tHandle.value = ResidualAccessor{tracks.eT_fit, tracks.t_time};
0187     tHandle.error = DirectAccessor<float>{tracks.err_eT_fit};
0188     tHandle.accept = AcceptAll{};
0189     // Push it
0190     baseResidualPulls.push_back(tHandle);
0191   }
0192
0193   if (residualPulls.test(6)) {
0194     // The pt measurement
0195     ResidualPullHandle ptHandle;
0196     ptHandle.tag = "pt";
0197     ptHandle.residualStr = "p_{T}^{rec} - p_{T}^{true}";
0198     ptHandle.residualUnit = "[GeV]";
0199     ptHandle.errorStr = "#sigma(p_{T})";
0200     ptHandle.rangeDrawStr = "1./abs(eQOP_fit) * sin(eTHETA_fit) - (1./t_pT)";
0201     ptHandle.rangeMaxStr = "(1000,-10.,10.)";
0202     ptHandle.value =
0203         PtResidualAccessor{tracks.eQOP_fit, tracks.eTHETA_fit, tracks.t_pT};
0204     ptHandle.error = PtErrorAccessor{tracks.eQOP_fit, tracks.err_eQOP_fit,
0205                                      tracks.eTHETA_fit, tracks.err_eTHETA_fit};
0206     ptHandle.accept = AcceptAll{};
0207     // Push it
0208     baseResidualPulls.push_back(ptHandle);
0209   }
0210
0211   using Auxiliaries = std::vector<SingleHandle>;
0212   Auxiliaries baseAuxilaries;
0213
0214   if (auxiliary.test(0)) {
0215     // Chi2/ndf
0216     SingleHandle chi2ndf;
0217     chi2ndf.tag = "chi2ndf";
0218     chi2ndf.label = "#Chi^{2}/ndf";
0219     chi2ndf.bins = nHistBins;
0220     chi2ndf.range = {0., 5.};
0221     chi2ndf.value =
0222         DivisionAccessor<float, unsigned int>{tracks.chi2Sum, tracks.NDF};
0223     chi2ndf.accept = AcceptAll{};
0224     // Push it
0225     baseAuxilaries.push_back(chi2ndf);
0226   }
0227
0228   if (auxiliary.test(1)) {
0229     // Measurements
0230     SingleHandle measurements;
0231     measurements.tag = "measurements";
0232     measurements.label = "#(measurements)";
0233     measurements.rangeDrawStr = "nMeasurements";
0234     measurements.value = DirectAccessor<unsigned int>{tracks.nMeasurements};
0235     measurements.accept = AcceptAll{};
0236     estimateIntegerRange(measurements, *rangeCanvas, *tracks.tree, peakEntries,
0237                          250, 5, ++histBarcode);
0238     // Push it
0239     baseAuxilaries.push_back(measurements);
0240   }
0241
0242   if (auxiliary.test(2)) {
0243     // Holes
0244     SingleHandle holes;
0245     holes.tag = "holes";
0246     holes.label = "#(holes)";
0247     holes.rangeDrawStr = "nHoles";
0248     holes.value = DirectAccessor<unsigned int>{tracks.nHoles};
0249     holes.accept = AcceptAll{};
0250     estimateIntegerRange(holes, *rangeCanvas, *tracks.tree, peakEntries, 10, 2,
0251                          ++histBarcode);
0252     // Push it
0253     baseAuxilaries.push_back(holes);
0254   }
0255
0256   if (auxiliary.test(3)) {
0257     // Holes
0258     SingleHandle outliers;
0259     outliers.tag = "outliers";
0260     outliers.label = "#(outliers)";
0261     outliers.rangeDrawStr = "nOutliers";
0262     outliers.value = DirectAccessor<unsigned int>{tracks.nOutliers};
0263     outliers.accept = AcceptAll{};
0264     estimateIntegerRange(outliers, *rangeCanvas, *tracks.tree, peakEntries, 10,
0265                          2, ++histBarcode);
0266     // Push it
0267     baseAuxilaries.push_back(outliers);
0268   }
0269
0270   if (auxiliary.test(4)) {
0271     // Holes
0272     SingleHandle shared;
0273     shared.tag = "shared";
0274     shared.label = "#(shared)";
0275     shared.rangeDrawStr = "nSharedHits";
0276     shared.value = DirectAccessor<unsigned int>{tracks.nSharedHits};
0277     shared.accept = AcceptAll{};
0278     estimateIntegerRange(shared, *rangeCanvas, *tracks.tree, peakEntries, 10, 2,
0279                          ++histBarcode);
0280     // Push it
0281     baseAuxilaries.push_back(shared);
0282   }
0283
0284   // Preparation phase for handles
0285 #ifdef NLOHMANN_AVAILABLE
0286   nlohmann::json handle_configs;
0287   if (! inConfig.empty()) {
0288     std::ifstream ifs(inConfig.c_str());
0289     handle_configs = nlohmann::json::parse(ifs);
0290   }
0291
0292   /// Helper method to handle the range, it is either read in or estimated
0293   ///
0294   /// It attempts to read in from a json input, if that does not work,
0295   /// it will estimate it (time consuming)
0296   ///
0297   /// @param handle the parameter handle in question
0298   /// @param handleTag the unique tangle tag
0299   /// @param peakE the number of entries used for range peaking
0300   auto handleRange = [&](ResidualPullHandle& handle, const TString& handleTag, unsigned int peakE) -> void {
0301     bool rangeDetermined = false;
0302     if (! inConfig.empty()) {
0303       if (handle_configs.contains((handleTag).Data())) {
0304         auto handle_config = handle_configs[(handleTag).Data()];
0305         handle.range = handle_config["range"].get<std::array<float, 2>>();
0306         rangeDetermined = true;
0307       }
0308     }
0309     if (! rangeDetermined) {
0310       estimateResiudalRange(handle, *rangeCanvas, *tracks.tree, peakE,
0311                             ++histBarcode);
0312     }
0313
0314     if (! outConfig.empty()) {
0315       nlohmann::json range_config;
0316       range_config["range"] = handle.range;
0317       handle_configs[(handleTag).Data()] = range_config;
0318     }
0319   };
0320 #else
0321   /// Helper method to handle range without possibility to read in
0322   ///
0323   /// @param handle the parameter handle in question
0324   /// @param handleTag the unique tangle tag
0325   /// @param peakEntries the number of entries used for range peaking
0326   auto handleRange = [&](ResidualPullHandle& handle, const TString& handleTag,
0327                          unsigned long peakEntries) -> void {
0328     estimateResiudalRange(handle, *rangeCanvas, *tracks.tree, peakEntries,
0329                           ++histBarcode);
0330   };
0331 #endif
0332
0333   // Full Range handles - they accept all tracks
0334   ResidualPulls fullResidualPulls = baseResidualPulls;
0335   // Range Estimation and booking histogram
0336   for (auto& fHandle : fullResidualPulls) {
0337     // The full tag
0338     TString handleTag(fHandle.tag + std::string("_all"));
0339     handleRange(fHandle, handleTag, peakEntries);
0340     // Book histograms
0341     bookHistograms(fHandle, pullRange, nHistBins, ++histBarcode);
0342     // The Histogram names
0343     TString residualN = TString("res_") + handleTag;
0344     TString pullN = TString("pull_") + handleTag;
0345     // Style and name
0346     setHistStyle(fHandle.residualHist);
0347     fHandle.residualHist->SetName(residualN);
0348     setHistStyle(fHandle.pullHist);
0349     fHandle.pullHist->SetName(pullN);
0350   }
0351
0352   // Regional/Binned  handles
0353   using ResidualPullsVector = std::vector<ResidualPulls>;
0354   using ResidualPullsMatrix = std::vector<ResidualPullsVector>;
0355
0356   // Eta-Pt residual/pull handles
0357   ResidualPullsVector ptResidualPulls =
0358       ResidualPullsVector(nPtBins, baseResidualPulls);
0359   ResidualPullsMatrix etaPtResidualPulls =
0360       ResidualPullsMatrix(nEtaBins, ptResidualPulls);
0361
0362   // Eta-Phi residual/pull handles
0363   ResidualPullsVector phiResidualPulls =
0364       ResidualPullsVector(nPhiBins, baseResidualPulls);
0365   ResidualPullsMatrix etaPhiResidualPulls =
0366       ResidualPullsMatrix(nEtaBins, phiResidualPulls);
0367
0368   Auxiliaries fullAuxiliaries = baseAuxilaries;
0369
0370   // Histogram booking for full range auxiliaries
0371   for (auto& fAuxiliary : fullAuxiliaries) {
0372     fAuxiliary.hist =
0373         new TH1F(fAuxiliary.tag.c_str(), fAuxiliary.tag.c_str(),
0374                  fAuxiliary.bins, fAuxiliary.range[0], fAuxiliary.range[1]);
0375     fAuxiliary.hist->GetXaxis()->SetTitle(fAuxiliary.label.c_str());
0376     fAuxiliary.hist->GetYaxis()->SetTitle("Entries");
0377     setHistStyle(fAuxiliary.hist);
0378   }
0379
0380   using AuxiliariesVector = std::vector<Auxiliaries>;
0381   using AuxiliariesMatrix = std::vector<AuxiliariesVector>;
0382
0383   // Eta-Pt auxiliaries
0384   AuxiliariesVector ptAuxiliaries = AuxiliariesVector(nPtBins, baseAuxilaries);
0385   AuxiliariesMatrix etaPtAuxiliaries =
0386       AuxiliariesMatrix(nEtaBins, ptAuxiliaries);
0387
0388   // Eta-Phi auxiliaries
0389   AuxiliariesVector phiAuxiliaries =
0390       AuxiliariesVector(nPhiBins, baseAuxilaries);
0391   AuxiliariesMatrix etaPhiAuxiliaries =
0392       AuxiliariesMatrix(nEtaBins, phiAuxiliaries);
0393
0394   // Loop over the binned handles & fill the acceptors
0395   float phiStep = (phiRange[1] - phiRange[0]) / nPhiBins;
0396   float etaStep = (etaRange[1] - etaRange[0]) / nEtaBins;
0397
0398   TVectorF phiVals(nPhiBins + 1);
0399   TVectorF etaVals(nEtaBins + 1);
0400   TVectorF ptVals(nPtBins + 1);
0401
0402   phiVals[0] = phiRange[0];
0403   etaVals[0] = etaRange[0];
0404   ptVals[0] = ptBorders[0];
0405
0406 #ifdef BOOST_AVAILABLE
0407   std::cout << "*** Handle Preparation: " << std::endl;
0408   progress_display handle_preparation_progress(
0409       nPhiBins * nEtaBins * nPtBins * baseResidualPulls.size());
0410 #endif
0411
0412   std::string land = " && ";
0413
0414   /// Preparation of handles / acceptance range
0415   for (unsigned int iphi = 0; iphi < nPhiBins; ++iphi) {
0416     // Prepare the phi range for this batch
0417     float phiMin = phiRange[0] + iphi * phiStep;
0418     float phiMax = phiRange[0] + (iphi + 1) * phiStep;
0419     phiVals[iphi + 1] = phiMax;
0420
0421     // Acceptance range
0422     AcceptRange phiBin{tracks.t_phi, {phiMin, phiMax}};
0423     // Name tag
0424     TString phiTag = "_phi";
0425     phiTag += iphi;
0426     // Range cut string
0427     std::string phiCut = "t_phi >= ";
0428     phiCut += std::to_string(phiMin);
0429     phiCut += land;
0430     phiCut += std::string("t_phi < ");
0431     phiCut += std::to_string(phiMax);
0432
0433     for (unsigned int ieta = 0; ieta < nEtaBins; ++ieta) {
0434       // Prepare the eta ragne for this batch
0435       float etaMin = etaRange[0] + ieta * etaStep;
0436       float etaMax = etaRange[0] + (ieta + 1) * etaStep;
0437       etaVals[ieta + 1] = etaMax;
0438       // Acceptance range
0439       AcceptRange etaBin{tracks.t_eta, {etaMin, etaMax}};
0440       // Name tag
0441       TString etaTag = "_eta";
0442       etaTag += ieta;
0443       // Range cut string
0444       std::string etaCut = "t_eta >= ";
0445       etaCut += std::to_string(etaMin);
0446       etaCut += land;
0447       etaCut += std::string("t_eta < ");
0448       etaCut += std::to_string(etaMax);
0449
0450       // Combined eta/phi acceptance
0451       AcceptCombination etaPhiBin{etaBin, phiBin};
0452       TString etaPhiTag = etaTag + phiTag;
0453
0454       for (unsigned int ipt = 0; ipt < nPtBins; ++ipt) {
0455         // Acceptance range for this pT bin
0456         float ptMin = static_cast<float>(ptBorders[ipt]);
0457         float ptMax = static_cast<float>(ptBorders[ipt + 1]);
0458         AcceptRange ptBin{tracks.t_pT, {ptMin, ptMax}};
0459
0460         float upperPtBorder =
0461             ptBorders[ipt + 1] == std::numeric_limits<double>::infinity()
0462                 ? 100.
0463                 : ptBorders[ipt + 1];
0464         ptVals[ipt + 1] = upperPtBorder;
0465         // Name tag
0466         TString ptTag = "_pt";
0467         ptTag += ipt;
0468
0469         // Range cut string
0470         std::string ptCut = "t_pT >= ";
0471         ptCut += std::to_string(ptMin);
0472         ptCut += land;
0473         ptCut += std::string("t_pT < ");
0474         ptCut += std::to_string(ptMax);
0475
0476         // Combined eta/pt acceptance
0477         AcceptCombination etaPtBin{etaBin, ptBin};
0478
0479         for (unsigned int iresp = 0; iresp < baseResidualPulls.size();
0480              ++iresp) {
0481           // Eta-Pt handles -- restrict for iphi == 0
0482           if (iphi == 0) {
0483             auto& etaPtHandle = etaPtResidualPulls[ieta][ipt][iresp];
0484             // Create the handle tag
0485             TString handleTag(etaPtHandle.tag + etaTag + ptTag);
0486             // Accept range and range cut
0487             etaPtHandle.accept = etaPtBin;
0488             etaPtHandle.rangeCutStr = ptCut + land + etaCut;
0489             handleRange(etaPtHandle, handleTag, peakEntries);
0490             bookHistograms(etaPtHandle, pullRange, nHistBins, ++histBarcode);
0491
0492             // Set name and style/
0493             TString residualN = TString("res_") + handleTag;
0494             TString pullN = TString("pull_") + handleTag;
0495             // Range histogram does not exist when read in from configuration
0496             if (etaPtHandle.rangeHist != nullptr) {
0497               TString rangeN = TString("range_") + handleTag;
0498               setHistStyle(etaPtHandle.rangeHist);
0499               etaPtHandle.rangeHist->SetName(rangeN);
0500             }
0501             setHistStyle(etaPtHandle.residualHist);
0502             etaPtHandle.residualHist->SetName(residualN);
0503             setHistStyle(etaPtHandle.pullHist);
0504             etaPtHandle.pullHist->SetName(pullN);
0505           }
0506           // Eta-Phi handles --- restrice for ipt == 0
0507           if (ipt == 0) {
0508             auto& etaPhiHandle = etaPhiResidualPulls[ieta][iphi][iresp];
0509             // Create the handle tag
0510             TString handleTag(etaPhiHandle.tag + etaTag + phiTag);
0511             etaPhiHandle.accept = etaPhiBin;
0512             handleRange(etaPhiHandle, handleTag, peakEntries);
0513             bookHistograms(etaPhiHandle, pullRange, nHistBins, ++histBarcode);
0514
0515             // Set name and style
0516             TString residualN = TString("res_") + handleTag;
0517             TString pullN = TString("pull_") + handleTag;
0518
0519             // Range histogram does not exist when read in from configuration
0520             if (etaPhiHandle.rangeHist != nullptr) {
0521               TString rangeN = TString("range_") + handleTag;
0522               setHistStyle(etaPhiHandle.rangeHist);
0523               etaPhiHandle.rangeHist->SetName(rangeN);
0524             }
0525             setHistStyle(etaPhiHandle.residualHist);
0526             etaPhiHandle.residualHist->SetName(residualN);
0527             setHistStyle(etaPhiHandle.pullHist);
0528             etaPhiHandle.pullHist->SetName(pullN);
0529           }
0530
0531 #ifdef BOOST_AVAILABLE
0532           ++handle_preparation_progress;
0533 #endif
0534         }
0535
0536         // Auxiliary plots
0537         for (unsigned int iaux = 0; iaux < baseAuxilaries.size(); ++iaux) {
0538           // Eta-Pt handles - restrict to iphi == 0
0539           if (iphi == 0) {
0540             auto& etaPtAux = etaPtAuxiliaries[ieta][ipt][iaux];
0541             etaPtAux.accept = etaPtBin;
0542             TString handleTag(etaPtAux.tag + etaTag + ptTag);
0543             etaPtAux.hist =
0544                 new TH1F(handleTag.Data(), etaPtAux.tag.c_str(), etaPtAux.bins,
0545                          etaPtAux.range[0], etaPtAux.range[1]);
0546             etaPtAux.hist->GetXaxis()->SetTitle(etaPtAux.label.c_str());
0547             etaPtAux.hist->GetYaxis()->SetTitle("Entries");
0548             setHistStyle(etaPtAux.hist);
0549           }
0550
0551           // Eta-Phi handles - restrict to ipt == 0
0552           if (ipt == 0) {
0553             auto& etaPhiAux = etaPhiAuxiliaries[ieta][iphi][iaux];
0554             etaPhiAux.accept = etaPhiBin;
0555             TString handleTag(etaPhiAux.tag + etaTag + phiTag);
0556             etaPhiAux.hist = new TH1F(handleTag.Data(), etaPhiAux.tag.c_str(),
0557                                       etaPhiAux.bins, etaPhiAux.range[0],
0558                                       etaPhiAux.range[1]);
0559             etaPhiAux.hist->GetXaxis()->SetTitle(etaPhiAux.label.c_str());
0560             etaPhiAux.hist->GetYaxis()->SetTitle("Entries");
0561             setHistStyle(etaPhiAux.hist);
0562           }
0563         }
0564       }
0565     }
0566   }
0567
0568 #ifdef NLOHMANN_AVAILABLE
0569   if (! outConfig.empty()) {
0570     std::ofstream config_out;
0571     config_out.open(outConfig.c_str());
0572     config_out << handle_configs.dump(4);
0573   }
0574 #endif
0575
0576 #ifdef BOOST_AVAILABLE
0577   std::cout << "*** Event Loop: " << std::endl;
0578   progress_display event_loop_progress(entries);
0579 #endif
0580
0581   for (unsigned long ie = 0; ie < entries; ++ie) {
0582 #ifdef BOOST_AVAILABLE
0583     ++event_loop_progress;
0584 #endif
0585
0586     // Make sure you have the entry
0587     tracks.tree->GetEntry(ie);
0588     std::size_t nTracks = tracks.hasFittedParams->size();
0589     for (std::size_t it = 0; it < nTracks; ++it) {
0590       if (tracks.hasFittedParams->at(it)) {
0591         // Residual handlesL
0592         // Full range handles
0593         for (auto& fHandle : fullResidualPulls) {
0594           fHandle.fill(it);
0595         }
0596
0597         // Eta-Pt handles
0598         for (auto& etaBatch : etaPtResidualPulls) {
0599           for (auto& ptBatch : etaBatch) {
0600             for (auto& bHandle : ptBatch) {
0601               bHandle.fill(it);
0602             }
0603           }
0604         }
0605
0606         // Eta-Phi handles
0607         for (auto& etaBatch : etaPhiResidualPulls) {
0608           for (auto& phiBatch : etaBatch) {
0609             for (auto& bHandle : phiBatch) {
0610               bHandle.fill(it);
0611             }
0612           }
0613         }
0614       }
0615
0616       // Auxiliary handles:
0617       // Full range handles
0618       for (auto& fAuxiliary : fullAuxiliaries) {
0619         fAuxiliary.fill(it);
0620       }
0621
0622       // Eta-Pt handles
0623       for (auto& etaBatch : etaPtAuxiliaries) {
0624         for (auto& ptBatch : etaBatch) {
0625           for (auto& bHandle : ptBatch) {
0626             bHandle.fill(it);
0627           }
0628         }
0629       }
0630
0631       // Eta-Phi handles
0632       for (auto& etaBatch : etaPhiAuxiliaries) {
0633         for (auto& phiBatch : etaBatch) {
0634           for (auto& bHandle : phiBatch) {
0635             bHandle.fill(it);
0636           }
0637         }
0638       }
0639     }
0640   }
0641
0642   // The output file section
0643   auto output = TFile::Open(outFile.c_str(), "recreate");
0644   output->cd();
0645
0646   // Full range handles : residual and pulls
0647   for (auto& fHandle : fullResidualPulls) {
0648     if (fHandle.rangeHist != nullptr) {
0649       fHandle.rangeHist->Write();
0650     }
0651     fHandle.residualHist->Write();
0652     fHandle.pullHist->Write();
0653   }
0654
0655   // Full range handles : auxiliaries
0656   for (auto& fAuxiliary : fullAuxiliaries) {
0657     fAuxiliary.hist->SetName(fAuxiliary.hist->GetName() + TString("_all"));
0658     fAuxiliary.hist->Write();
0659   }
0660
0661   struct SummaryHistograms {
0662     TH2F* fillStats = nullptr;
0663
0664     std::vector<TH2F*> residualRMS;
0665     std::vector<TH2F*> residualMean;
0666     std::vector<TH2F*> pullSigma;
0667     std::vector<TH2F*> pullMean;
0668
0669     std::vector<TH2F*> auxiliaries;
0670   };
0671
0672   /// Helper method to analyse bins
0673   /// @param residualPullsMatrix the 2D matrix of handles
0674   /// @param auxiliaryMatrix the 2D matrix of the auxiliary handles
0675   /// @param matrixTag the identification tag for the matrix
0676   /// @param outputBorders the border vector for the outer bins
0677   /// @param innerBorders the border vector for the inner bins
0678   /// @param fXTitle the title of the x axis of the first projection
0679   /// @param sXTitle the title of the x axis of the second projection
0680   ///
0681   /// @note this is a void function
0682   auto analyseBins = [&](ResidualPullsMatrix& residualPullsMatrix,
0683                          AuxiliariesMatrix& auxiliaryMatrix,
0684                          const TString& matrixTag, const TVectorF& outerBorders,
0685                          const TVectorF& innerBorders,
0686                          const TString& fXTitle = "#eta",
0687                          const TString& sXTitle = "#phi") -> void {
0688     // The summary histogram set
0689     SummaryHistograms summary;
0690
0691     // 2D handles ---------------------------
0692     unsigned int nOuterBins = outerBorders.GetNrows() - 1;
0693     auto outerValues = outerBorders.GetMatrixArray();
0694     unsigned int nInnerBins = innerBorders.GetNrows() - 1;
0695     auto innerValues = innerBorders.GetMatrixArray();
0696
0697     TString statN = TString("entries") + matrixTag;
0698     summary.fillStats =
0699         new TH2F(statN, "", nOuterBins, outerValues, nInnerBins, innerValues);
0700
0701 #ifdef BOOST_AVAILABLE
0702     progress_display analysis_progress(nOuterBins * nInnerBins);
0703 #endif
0704
0705     // Prepare by looping over the base bHandles - residuals
0706     for (auto& bHandle : baseResidualPulls) {
0707       // Create a unique handle tag
0708       TString handleTag = TString(bHandle.tag) + matrixTag;
0709       // ... and names
0710       TString residualRMSN = TString("res_rms_") + handleTag;
0711       TString residualMeanN = TString("res_mean_") + handleTag;
0712       TString pullSigmaN = TString("pull_sigma_") + handleTag;
0713       TString pullMeanN = TString("pull_mean_") + handleTag;
0714
0715       TH2F* residualRMS = new TH2F(residualRMSN, "", nOuterBins, outerValues,
0716                                    nInnerBins, innerValues);
0717       TH2F* residualMean = new TH2F(residualMeanN, "", nOuterBins, outerValues,
0718                                     nInnerBins, innerValues);
0719       TH2F* pullSigma = new TH2F(pullSigmaN, "", nOuterBins, outerValues,
0720                                  nInnerBins, innerValues);
0721       TH2F* pullMean = new TH2F(pullMeanN, "", nOuterBins, outerValues,
0722                                 nInnerBins, innerValues);
0723       // Booked & ready
0724       summary.residualRMS.push_back(residualRMS);
0725       summary.residualMean.push_back(residualMean);
0726       summary.pullSigma.push_back(pullSigma);
0727       summary.pullMean.push_back(pullMean);
0728     }
0729
0730     // Prepare by looping over the base handles - auxiliaries
0731     for (auto& aHandle : baseAuxilaries) {
0732       // Create a unique handle tag
0733       TString auxiliaryTag = TString(aHandle.tag) + matrixTag;
0734       TH2F* auxHist = new TH2F(auxiliaryTag, auxiliaryTag, nOuterBins,
0735                                  outerValues, nInnerBins, innerValues);
0736       summary.auxiliaries.push_back(auxHist);
0737     }
0738
0739     unsigned int io = 0;
0740     for (auto& outerBatch : residualPullsMatrix) {
0741       unsigned int ii = 0;
0742       for (auto& innerBatch : outerBatch) {
0743         // residual/pull loop
0744         unsigned int iresp = 0;
0745         for (auto& bHandle : innerBatch) {
0746           // Range estimates / could be empty if read from configuration
0747           if (bHandle.rangeHist != nullptr) {
0748             bHandle.rangeHist->Write();
0749           }
0750           // Fill the stats
0751           if (iresp == 0) {
0752             summary.fillStats->SetBinContent(io + 1, ii + 1, bHandle.accepted);
0753           }
0754
0755           // Residuals
0756           // Get RMS/ RMSError
0757           float rrms = bHandle.residualHist->GetRMS();
0758           float rrerr = bHandle.residualHist->GetRMSError();
0759           float rmean = bHandle.residualHist->GetMean();
0760           float rmerr = bHandle.residualHist->GetMeanError();
0761           summary.residualRMS[iresp]->SetBinContent(io + 1, ii + 1, rrms);
0762           summary.residualRMS[iresp]->SetBinError(io + 1, ii + 1, rrerr);
0763           summary.residualMean[iresp]->SetBinContent(io + 1, ii + 1, rmean);
0764           summary.residualMean[iresp]->SetBinError(io + 1, ii + 1, rmerr);
0765           bHandle.residualHist->Write();
0766           // Pulls
0767           bHandle.pullHist->Fit("gaus", "q");
0768           TF1* gauss = bHandle.pullHist->GetFunction("gaus");
0769           if (gauss != nullptr) {
0770             float pmu = gauss->GetParameter(1);
0771             float pmerr = gauss->GetParError(1);
0772             float psigma = gauss->GetParameter(2);
0773             float pserr = gauss->GetParError(2);
0774             summary.pullSigma[iresp]->SetBinContent(io + 1, ii + 1, psigma);
0775             summary.pullSigma[iresp]->SetBinError(io + 1, ii + 1, pserr);
0776             summary.pullMean[iresp]->SetBinContent(io + 1, ii + 1, pmu);
0777             summary.pullMean[iresp]->SetBinError(io + 1, ii + 1, pmerr);
0778           }
0779           bHandle.pullHist->Write();
0780           ++iresp;
0781         }
0782
0783         // auxilaiary loop
0784         auto auxiliaryBatch = auxiliaryMatrix[io][ii];
0785         unsigned int iaux = 0;
0786         for (auto& aHandle : auxiliaryBatch) {
0787           float value = aHandle.hist->GetMean();
0788           float error = aHandle.hist->GetMeanError();
0789           summary.auxiliaries[iaux]->SetBinContent(io + 1, ii + 1, value);
0790           summary.auxiliaries[iaux]->SetBinError(io + 1, ii + 1, error);
0791           aHandle.hist->Write();
0792
0793           ++iaux;
0794         }
0795 #ifdef BOOST_AVAILABLE
0796         ++analysis_progress;
0797 #endif
0798         ++ii;
0799       }
0800       ++io;
0801     }
0802
0803     /// Write out the projection histogram
0804     ///
0805     /// @param h2 the 2D histogram as source for the projects
0806     /// @param fXTitle the title of the x axis of the first projection
0807     /// @param fYTitle the title of the y axis of the first projection
0808     /// @param sXTitle the title of the x axis of the second projection
0809     /// @param sYTitle the title of the y axis of the second projection
0810     auto writeProjections = [](const TH2F& h2, const TString& fXTitleP = "#eta",
0811                                const TString& fYTitleP = "sigma",
0812                                const TString& sXTitleP = "#phi",
0813                                const TString& sYTitleP = "sigma") -> void {
0814       const TString& fTag = "_pX";
0815       const TString& sTag = "_pY";
0816
0817       unsigned int nBinsX = h2.GetXaxis()->GetNbins();
0818       unsigned int nBinsY = h2.GetYaxis()->GetNbins();
0819       if (fTag != "") {
0820         TH1D* pX =
0821             dynamic_cast<TH1D*>(h2.ProjectionX((h2.GetName() + fTag).Data()));
0822         setHistStyle(pX);
0823         if (pX != nullptr) {
0824           pX->GetXaxis()->SetTitle(fXTitleP.Data());
0825           pX->GetYaxis()->SetTitle(fYTitleP.Data());
0826           pX->Write();
0827         }
0828         // Bin-wise projections
0829         for (unsigned int iy = 1; iy <= nBinsY; ++iy) {
0830           pX = dynamic_cast<TH1D*>(h2.ProjectionX(
0831               (h2.GetName() + fTag + sTag + (iy - 1)).Data(), iy, iy));
0832           setHistStyle(pX);
0833           if (pX != nullptr) {
0834             pX->GetXaxis()->SetTitle(fXTitleP.Data());
0835             pX->GetYaxis()->SetTitle(fYTitleP.Data());
0836             pX->Write();
0837           }
0838         }
0839       }
0840       if (sTag != "") {
0841         TH1D* pY =
0842             dynamic_cast<TH1D*>(h2.ProjectionY((h2.GetName() + sTag).Data()));
0843         setHistStyle(pY);
0844         if (pY != nullptr) {
0845           pY->GetXaxis()->SetTitle(sXTitleP.Data());
0846           pY->GetYaxis()->SetTitle(sYTitleP.Data());
0847           pY->Write();
0848         }
0849         // Bin-wise projections
0850         for (unsigned int ix = 1; ix <= nBinsX; ++ix) {
0851           pY = dynamic_cast<TH1D*>(h2.ProjectionY(
0852               (h2.GetName() + sTag + fTag + (ix - 1)).Data(), ix, ix));
0853           setHistStyle(pY);
0854           if (pY != nullptr) {
0855             pY->GetXaxis()->SetTitle(sXTitleP.Data());
0856             pY->GetYaxis()->SetTitle(sYTitleP.Data());
0857             pY->Write();
0858           }
0859         }
0860       }
0861     };
0862
0863     setHistStyle(summary.fillStats);
0864     summary.fillStats->Write();
0865
0866     // Write mapped residual/pull histograms and their projections
0867     for (unsigned int iresp = 0; iresp < baseResidualPulls.size(); ++iresp) {
0868       // Get the handle for writing out
0869       auto bHandle = baseResidualPulls[iresp];
0870
0871       TString rrms = TString("RMS[") + bHandle.residualStr + TString("] ") +
0872                      bHandle.residualUnit;
0873       TString rmu = TString("#mu[") + bHandle.residualStr + TString("] ") +
0874                     bHandle.residualUnit;
0875
0876       TString psigma = TString("#sigma[(") + bHandle.residualStr +
0877                        TString(")/") + bHandle.errorStr + TString("]");
0878       TString pmu = TString("#mu[(") + bHandle.residualStr + TString(")/") +
0879                     bHandle.errorStr + TString("]");
0880
0881       // 2D map histograms
0882       setHistStyle(summary.residualRMS[iresp]);
0883       summary.residualRMS[iresp]->GetXaxis()->SetTitle(fXTitle);
0884       summary.residualRMS[iresp]->GetYaxis()->SetTitle(sXTitle);
0885       summary.residualRMS[iresp]->GetZaxis()->SetTitle(rrms);
0886       summary.residualRMS[iresp]->Write();
0887
0888       setHistStyle(summary.residualMean[iresp]);
0889       summary.residualMean[iresp]->GetXaxis()->SetTitle(fXTitle);
0890       summary.residualMean[iresp]->GetYaxis()->SetTitle(sXTitle);
0891       summary.residualMean[iresp]->GetZaxis()->SetTitle(rmu);
0892       summary.residualMean[iresp]->Write();
0893
0894       setHistStyle(summary.pullSigma[iresp]);
0895       adaptColorPalette(summary.pullSigma[iresp], 0., 4., 1., 0.1, 104);
0896       summary.pullSigma[iresp]->GetXaxis()->SetTitle(fXTitle);
0897       summary.pullSigma[iresp]->GetYaxis()->SetTitle(sXTitle);
0898       summary.pullSigma[iresp]->GetZaxis()->SetRangeUser(0., 4.);
0899       summary.pullSigma[iresp]->GetZaxis()->SetTitle(psigma);
0900       summary.pullSigma[iresp]->Write();
0901
0902       setHistStyle(summary.pullMean[iresp]);
0903       adaptColorPalette(summary.pullMean[iresp], -1., 1., 0., 0.1, 104);
0904       summary.pullMean[iresp]->GetXaxis()->SetTitle(fXTitle);
0905       summary.pullMean[iresp]->GetYaxis()->SetTitle(sXTitle);
0906       summary.pullMean[iresp]->GetZaxis()->SetRangeUser(-1., 1.);
0907       summary.pullMean[iresp]->GetZaxis()->SetTitle(pmu);
0908       summary.pullMean[iresp]->Write();
0909
0910       // Write the projection histograms
0911       writeProjections(*summary.residualRMS[iresp], fXTitle, rrms, sXTitle,
0912                        rrms);
0913       writeProjections(*summary.residualMean[iresp], fXTitle, rmu, sXTitle,
0914                        rmu);
0915       writeProjections(*summary.pullSigma[iresp], fXTitle, psigma, sXTitle,
0916                        psigma);
0917       writeProjections(*summary.pullMean[iresp], fXTitle, pmu, sXTitle, pmu);
0918     }
0919
0920     // Write mapped auxiliary histograms and their projections
0921     for (unsigned int iaux = 0; iaux < baseAuxilaries.size(); ++iaux) {
0922       setHistStyle(summary.auxiliaries[iaux]);
0923       summary.auxiliaries[iaux]->GetXaxis()->SetTitle(fXTitle);
0924       summary.auxiliaries[iaux]->GetYaxis()->SetTitle(sXTitle);
0925       summary.auxiliaries[iaux]->GetZaxis()->SetTitle(
0926           baseAuxilaries[iaux].label.c_str());
0927       summary.auxiliaries[iaux]->Write();
0928
0929       writeProjections(*summary.auxiliaries[iaux], fXTitle,
0930                        baseAuxilaries[iaux].label, sXTitle,
0931                        baseAuxilaries[iaux].label);
0932     }
0933
0934     return;
0935   };
0936
0937 // The handle matrices
0938 #ifdef BOOST_AVAILABLE
0939   std::cout << "*** Bin/Projection Analysis: " << std::endl;
0940 #endif
0941   analyseBins(etaPtResidualPulls, etaPtAuxiliaries, TString("_eta_pt"), etaVals,
0942               ptVals, "#eta", "p_{T} [GeV]");
0943   analyseBins(etaPhiResidualPulls, etaPhiAuxiliaries, TString("_eta_phi"),
0944               etaVals, phiVals, "#eta", "#phi");
0945
0946   output->Close();
0947
0948   return 1;
0949 }