Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:12:09 

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 <map>
0012 #include <string>
0013 #include <vector>
0014 
0015 #include <TCanvas.h>
0016 #include <TColor.h>
0017 #include <TDirectory.h>
0018 #include <TError.h>
0019 #include <TF1.h>
0020 #include <TFile.h>
0021 #include <TH1F.h>
0022 #include <TH2F.h>
0023 #include <TLegend.h>
0024 #include <TMath.h>
0025 #include <TProfile2D.h>
0026 #include <TStyle.h>
0027 #include <TTree.h>
0028 
0029 #include "CommonUtils.h"
0030 #include "TreeReader.h"
0031 
0032 using namespace ROOT;
0033 
0034 /// Plot the bound parameter resolutions
0035 ///
0036 /// (loc1, phi, theta, q/p, t) at all track states from root file produced by
0037 /// the RootTrackStatesWriter
0038 ///
0039 /// @param inFile the input root file
0040 /// @param treeName the input tree name (default: 'trackstates)
0041 /// @param outFile the output root file
0042 /// @param pTypes the track parameter types (prt, flt, smt)
0043 /// @param saveAs the plot saving type
0044 int boundParamResolution(const std::string& inFile, const std::string& treeName,
0045                          const std::string& outFile, bool predicted = true,
0046                          bool filtered = true, bool smoothed = true,
0047                          bool fitFiltered = true, bool fitPredicted = true,
0048                          bool fitSmoothed = true,
0049                          const std::string& saveAs = "") {
0050   // Some style options
0051   gStyle->SetOptFit(0000);
0052   gStyle->SetOptStat(0000);
0053   gStyle->SetPadLeftMargin(0.20);
0054   gStyle->SetPadRightMargin(0.1);
0055   gStyle->SetPadTopMargin(0.1);
0056   gStyle->SetPadBottomMargin(0.15);
0057 
0058   // Pull range, residual ranges are automatically determined
0059   double pullRange = 5;
0060 
0061   //
0062   // Residual ranges should be largest in predicted and smallest in smoothed,
0063   // hence reverse the order.
0064   //
0065   std::string range_tag = "smt";
0066   if (predicted) {
0067     range_tag = "prt";
0068   } else if (filtered) {
0069     range_tag = "flt";
0070   }
0071 
0072   // Section 0: file handling ---------------------------------------------
0073   //
0074   // Open root file written by RootTrackWriter
0075   // Create output root file
0076   std::cout << "Opening file: " << inFile << std::endl;
0077   TFile* file = TFile::Open(inFile.c_str(), "read");
0078 
0079   // Bail out if no tree was found
0080   if (file == nullptr) {
0081     return -1;
0082   }
0083 
0084   std::cout << "Reading tree: " << treeName << std::endl;
0085   TTree* tree = static_cast<TTree*>(file->Get(treeName.c_str()));
0086 
0087   // Bail out if no tree was found
0088   if (tree == nullptr) {
0089     return -2;
0090   }
0091 
0092   TFile* out = TFile::Open(outFile.c_str(), "recreate");
0093   out->cd();
0094 
0095   // Section 1: geometry parsing ------------------------------------------
0096   //
0097   // Gathers accessible volume_id and layer_id values
0098   // Draw the volume_id : layer_id correlation matrix
0099   std::map<int, std::vector<int>> volLayIds;
0100   TCanvas* geometryCanvas =
0101       new TCanvas("volLayCanvas", "Volume Layer Matrix", 10, 10, 450, 600);
0102   geometryCanvas->Divide(1, 2);
0103   // Volume/Layer Id
0104   geometryCanvas->cd(1);
0105   tree->Draw("layer_id:volume_id>>volID_layID(100,0.5,100.5,100,0.5,100.5)", "",
0106              "colz");
0107   auto h2_volID_layID = dynamic_cast<TH2F*>(out->Get("volID_layID"));
0108   setHistStyle(h2_volID_layID, 2);
0109   h2_volID_layID->Write();
0110   // Geometry size
0111   geometryCanvas->cd(2);
0112   std::string rz_draw_string = "(sqrt(g_x_";
0113   rz_draw_string += range_tag;
0114   rz_draw_string += "*g_x_";
0115   rz_draw_string += range_tag;
0116   rz_draw_string += "+g_y_";
0117   rz_draw_string += range_tag;
0118   rz_draw_string += "*g_y_";
0119   rz_draw_string += range_tag;
0120   rz_draw_string += ")):g_z_";
0121   rz_draw_string += range_tag;
0122   rz_draw_string += ">>geo_dim";
0123   tree->Draw(rz_draw_string.c_str());
0124   auto h2_geo_dim = dynamic_cast<TH2F*>(out->Get("geo_dim"));
0125   setHistStyle(h2_geo_dim, 2);
0126   float detectorZ = 1.15 * h2_geo_dim->GetXaxis()->GetXmax();
0127   float detectorR = 1.15 * h2_geo_dim->GetYaxis()->GetXmax();
0128   h2_geo_dim->Write();
0129 
0130   // Save the plots on screen
0131   if (!saveAs.empty()) {
0132     geometryCanvas->SaveAs((std::string("all_vol_lay_ids.") + saveAs).c_str());
0133   }
0134 
0135   // Now determine the valid volume / layer bins
0136   int volBins = h2_volID_layID->GetXaxis()->GetNbins();
0137   int layBins = h2_volID_layID->GetYaxis()->GetNbins();
0138   // Add the overall plots
0139   volLayIds[-1] = {-1};
0140   // Go through volumes and add plots per volume
0141   for (int volID = 1; volID <= volBins; ++volID) {
0142     for (int layID = 1; layID <= layBins; ++layID) {
0143       if (h2_volID_layID->GetBinContent(volID, layID) != 0.) {
0144         if (!volLayIds.contains(volID)) {
0145           // First occurrence of this layer, add -1 for volume plots
0146           volLayIds[volID] = {-1, layID};
0147         } else {
0148           volLayIds[volID].push_back(layID);
0149         }
0150       }
0151     }
0152   }
0153 
0154   // Section 2: Branch assignment -----------------------------------------
0155   //
0156   // Helper for assigning branches to the input tree
0157   TrackStatesReader tsReader(tree, false);
0158 
0159   // Section 3: Histogram booking -----------------------------------------
0160 
0161   TCanvas* rangeCanvas =
0162       new TCanvas("rangeCanvas", "Range Estimation", 10, 10, 900, 600);
0163   rangeCanvas->Divide(3, 2);
0164 
0165   std::vector<std::string> res_ranges = {std::string("res_eLOC0_") + range_tag,
0166                                          std::string("res_eLOC1_") + range_tag,
0167                                          std::string("res_ePHI_") + range_tag,
0168                                          std::string("res_eTHETA_") + range_tag,
0169                                          std::string("res_eQOP_") + range_tag,
0170                                          std::string("res_eT_") + range_tag};
0171 
0172   /// Helper method to make a volLayId string for identification & cut
0173   ///
0174   /// ensures a unique identification
0175   auto volLayIdCut = [](int vol, int lay) -> std::array<std::string, 2> {
0176     if (vol < 0 && lay < 0) {
0177       return {std::string("all_"), ""};
0178     }
0179 
0180     std::string vlstr = "vol" + std::to_string(vol);
0181     std::string vlcut = "volume_id == " + std::to_string(vol);
0182     if (lay > 0) {
0183       vlstr += (std::string("_lay") + std::to_string(lay));
0184       vlcut += (std::string(" && layer_id == ") + std::to_string(lay));
0185     }
0186     vlstr += std::string("_");
0187     return {vlstr, vlcut};
0188   };
0189 
0190   /// Helper method to estimate the ranges
0191   ///
0192   /// Takes identification & cut from the first one
0193   auto histRanges =
0194       [&](const std::array<std::string, 2>& vlIdCut) -> std::array<float, 6> {
0195     std::array<float, 6> ranges = {0., 0., 0., 0., 0., 0.};
0196     for (unsigned int ir = 0; ir < 6; ++ir) {
0197       rangeCanvas->cd(ir + 1);
0198       std::string drawString = res_ranges[ir] + ">>";
0199       std::string histString =
0200           vlIdCut[0] + std::string("range_") + res_ranges[ir];
0201       tree->Draw((drawString + histString).c_str(), vlIdCut[1].c_str());
0202       auto h1_range = dynamic_cast<TH1F*>(out->Get(histString.c_str()));
0203       h1_range->Write();
0204       float range = pullRange * h1_range->GetRMS();
0205       ranges[ir] = range;
0206     }
0207     if (!saveAs.empty()) {
0208       rangeCanvas->SaveAs(
0209           (vlIdCut[0] + std::string("res_ranges.") + saveAs).c_str());
0210     }
0211     return ranges;
0212   };
0213 
0214   // Track parameter name
0215   std::vector<std::string> paramNames = {"loc0",   "loc1", "#phi",
0216                                          "#theta", "q/p",  "t"};
0217 
0218   // Book histograms (with adapted range):
0219   //
0220   // Global profile histograms : residuals/pulls
0221   std::array<TProfile2D*, 6> p2d_res_zr_prt{};
0222   std::array<TProfile2D*, 6> p2d_res_zr_flt{};
0223   std::array<TProfile2D*, 6> p2d_res_zr_smt{};
0224   std::array<TProfile2D*, 6> p2d_pull_zr_prt{};
0225   std::array<TProfile2D*, 6> p2d_pull_zr_flt{};
0226   std::array<TProfile2D*, 6> p2d_pull_zr_smt{};
0227 
0228   for (unsigned int ipar = 0; ipar < paramNames.size(); ++ipar) {
0229     const auto& par = paramNames[ipar];
0230 
0231     if (predicted) {
0232       p2d_res_zr_prt[ipar] =
0233           new TProfile2D(Form("all_prof_res_prt_%s", par.c_str()),
0234                          Form("residual profile of %s", par.c_str()), 100,
0235                          -detectorZ, detectorZ, 50, 0., detectorR);
0236 
0237       p2d_pull_zr_prt[ipar] =
0238           new TProfile2D(Form("all_prof_pull_prt_%s", par.c_str()),
0239                          Form("pull profile of %s", par.c_str()), 100,
0240                          -detectorZ, detectorZ, 50, 0., detectorR);
0241 
0242       p2d_res_zr_prt[ipar]->SetErrorOption("s");
0243       p2d_res_zr_prt[ipar]->GetXaxis()->SetTitle("z [mm]");
0244       p2d_res_zr_prt[ipar]->GetYaxis()->SetTitle("R [mm]");
0245       p2d_res_zr_prt[ipar]->GetZaxis()->SetTitle(Form("r_{%s}", par.c_str()));
0246       setHistStyle(p2d_res_zr_prt[ipar], 1);
0247 
0248       p2d_pull_zr_prt[ipar]->SetErrorOption("s");
0249       p2d_pull_zr_prt[ipar]->GetXaxis()->SetTitle("z [mm]");
0250       p2d_pull_zr_prt[ipar]->GetYaxis()->SetTitle("R [mm]");
0251       p2d_pull_zr_prt[ipar]->GetZaxis()->SetTitle(
0252           Form("pull_{%s}", par.c_str()));
0253       setHistStyle(p2d_pull_zr_prt[ipar], 1);
0254     }
0255 
0256     if (filtered) {
0257       p2d_res_zr_flt[ipar] =
0258           new TProfile2D(Form("all_prof_res_flt_%s", par.c_str()),
0259                          Form("residual profile of %s", par.c_str()), 100,
0260                          -detectorZ, detectorZ, 50, 0., detectorR);
0261       p2d_pull_zr_flt[ipar] =
0262           new TProfile2D(Form("all_prof_pull_flt_%s", par.c_str()),
0263                          Form("pull profile of %s", par.c_str()), 100,
0264                          -detectorZ, detectorZ, 50, 0., detectorR);
0265 
0266       p2d_res_zr_flt[ipar]->SetErrorOption("s");
0267       p2d_res_zr_flt[ipar]->GetXaxis()->SetTitle("z [mm]");
0268       p2d_res_zr_flt[ipar]->GetYaxis()->SetTitle("R [mm]");
0269       p2d_res_zr_flt[ipar]->GetZaxis()->SetTitle(Form("r_{%s}", par.c_str()));
0270       setHistStyle(p2d_res_zr_flt[ipar], 2);
0271 
0272       p2d_pull_zr_flt[ipar]->SetErrorOption("s");
0273       p2d_pull_zr_flt[ipar]->GetXaxis()->SetTitle("z [mm]");
0274       p2d_pull_zr_flt[ipar]->GetYaxis()->SetTitle("R [mm]");
0275       p2d_pull_zr_flt[ipar]->GetZaxis()->SetTitle(
0276           Form("pull_{%s}", par.c_str()));
0277       setHistStyle(p2d_pull_zr_flt[ipar], 2);
0278     }
0279 
0280     if (smoothed) {
0281       p2d_res_zr_smt[ipar] =
0282           new TProfile2D(Form("all_prof_res_smt_%s", par.c_str()),
0283                          Form("residual profile of %s", par.c_str()), 100,
0284                          -detectorZ, detectorZ, 50, 0., detectorR);
0285 
0286       p2d_pull_zr_smt[ipar] =
0287           new TProfile2D(Form("all_prof_pull_smt_%s", par.c_str()),
0288                          Form("pull profile of %s", par.c_str()), 100,
0289                          -detectorZ, detectorZ, 50, 0., detectorR);
0290 
0291       p2d_res_zr_smt[ipar]->SetErrorOption("s");
0292       p2d_pull_zr_smt[ipar]->SetErrorOption("s");
0293 
0294       p2d_res_zr_smt[ipar]->GetXaxis()->SetTitle("z [mm]");
0295       p2d_res_zr_smt[ipar]->GetYaxis()->SetTitle("R [mm]");
0296       p2d_res_zr_smt[ipar]->GetZaxis()->SetTitle(Form("r_{%s}", par.c_str()));
0297       setHistStyle(p2d_res_zr_smt[ipar], 4);
0298 
0299       p2d_pull_zr_smt[ipar]->GetXaxis()->SetTitle("z [mm]");
0300       p2d_pull_zr_smt[ipar]->GetYaxis()->SetTitle("R [mm]");
0301       p2d_pull_zr_smt[ipar]->GetZaxis()->SetTitle(
0302           Form("pull_{%s}", par.c_str()));
0303       setHistStyle(p2d_pull_zr_smt[ipar], 4);
0304     }
0305   }
0306 
0307   // Residual / Pull histograms
0308   std::map<std::string, TH1F*> res_prt;
0309   std::map<std::string, TH1F*> res_flt;
0310   std::map<std::string, TH1F*> res_smt;
0311   std::map<std::string, TH1F*> pull_prt;
0312   std::map<std::string, TH1F*> pull_flt;
0313   std::map<std::string, TH1F*> pull_smt;
0314 
0315   // - per layer (identified by vol, lay)
0316   // - per volume (identified by vol, -1)
0317   // - overall (identified by lay)
0318   for (const auto& [vol, layers] : volLayIds) {
0319     for (const auto& lay : layers) {
0320       // Estimate the ranges from smoothed
0321       auto vlIdCut = volLayIdCut(vol, lay);
0322       auto ranges = histRanges(vlIdCut);
0323 
0324       // Residual range
0325       std::map<std::pair<std::string, std::string>, double> paramResidualRange =
0326           {{{"loc0", "l_{0}"}, ranges[0]}, {{"loc1", "l_{1}"}, ranges[1]},
0327            {{"#phi", "#phi"}, ranges[2]},  {{"#theta", "#theta"}, ranges[3]},
0328            {{"q/p", "q/p"}, ranges[4]},    {{"t", "t"}, ranges[5]}};
0329 
0330       // Create the hists and set up for them
0331       for (const auto& [partwin, resRange] : paramResidualRange) {
0332         // histogram creation
0333         std::string par = partwin.first;
0334         std::string id_par = vlIdCut[0] + par;
0335 
0336         TString par_string(partwin.second.c_str());
0337         TString res_string =
0338             par_string + TString("^{rec} - ") + par_string + TString("^{true}");
0339 
0340         TString pull_string = TString("(") + res_string +
0341                               TString(")/#sigma_{") + par_string + TString("}");
0342 
0343         if (predicted) {
0344           res_prt[id_par] =
0345               new TH1F(Form((vlIdCut[0] + std::string("res_prt_%s")).c_str(),
0346                             par.c_str()),
0347                        Form("residual of %s", par.c_str()), 100, -1 * resRange,
0348                        resRange);
0349           res_prt[id_par]->GetXaxis()->SetTitle(res_string.Data());
0350           res_prt[id_par]->GetYaxis()->SetTitle("Entries");
0351           setHistStyle(res_prt[id_par], kRed);
0352 
0353           pull_prt[id_par] = new TH1F(
0354               Form((vlIdCut[0] + std::string("pull_prt_%s")).c_str(),
0355                    par.c_str()),
0356               Form("pull of %s", par.c_str()), 100, -1 * pullRange, pullRange);
0357           pull_prt[id_par]->GetXaxis()->SetTitle(pull_string.Data());
0358           pull_prt[id_par]->GetYaxis()->SetTitle("Arb. Units");
0359           setHistStyle(pull_prt[id_par], kRed);
0360         }
0361 
0362         if (filtered) {
0363           res_flt[id_par] =
0364               new TH1F(Form((vlIdCut[0] + std::string("res_flt_%s")).c_str(),
0365                             par.c_str()),
0366                        Form("residual of %s", par.c_str()), 100, -1 * resRange,
0367                        resRange);
0368           res_flt[id_par]->GetXaxis()->SetTitle(res_string.Data());
0369           res_flt[id_par]->GetYaxis()->SetTitle("Entries");
0370           setHistStyle(res_flt[id_par], kBlue);
0371 
0372           pull_flt[id_par] = new TH1F(
0373               Form((vlIdCut[0] + std::string("pull_flt_%s")).c_str(),
0374                    par.c_str()),
0375               Form("pull of %s", par.c_str()), 100, -1 * pullRange, pullRange);
0376           pull_flt[id_par]->GetXaxis()->SetTitle(pull_string.Data());
0377           pull_flt[id_par]->GetYaxis()->SetTitle("Arb. Units");
0378           setHistStyle(pull_flt[id_par], kBlue);
0379         }
0380 
0381         if (smoothed) {
0382           res_smt[id_par] =
0383               new TH1F(Form((vlIdCut[0] + std::string("res_smt_%s")).c_str(),
0384                             par.c_str()),
0385                        Form("residual of %s", par.c_str()), 100, -1 * resRange,
0386                        resRange);
0387 
0388           res_smt[id_par]->GetXaxis()->SetTitle(res_string.Data());
0389           res_smt[id_par]->GetYaxis()->SetTitle("Entries");
0390           setHistStyle(res_smt[id_par], kBlack);
0391 
0392           pull_smt[id_par] = new TH1F(
0393               Form((vlIdCut[0] + std::string("pull_smt_%s")).c_str(),
0394                    par.c_str()),
0395               Form("pull of %s", par.c_str()), 100, -1 * pullRange, pullRange);
0396 
0397           pull_smt[id_par]->GetXaxis()->SetTitle(pull_string.Data());
0398           pull_smt[id_par]->GetYaxis()->SetTitle("Arb. Units");
0399           setHistStyle(pull_smt[id_par], kBlack);
0400         }
0401       }
0402     }
0403   }
0404 
0405   // Section 4: Histogram filling -----------------------------------------
0406   //
0407   // - Running through the entries and filling the histograms
0408   int entries = tree->GetEntries();
0409   for (int j = 0; j < entries; j++) {
0410     tsReader.getEntry(j);
0411 
0412     for (unsigned int i = 0; i < tsReader.nMeasurements; i++) {
0413       // global profile filling
0414       if (predicted && tsReader.predicted->at(i)) {
0415         auto x_prt = tsReader.g_x_prt->at(i);
0416         auto y_prt = tsReader.g_y_prt->at(i);
0417         auto r_prt = std::hypot(x_prt, y_prt);
0418         auto z_prt = tsReader.g_z_prt->at(i);
0419         p2d_res_zr_prt[0]->Fill(z_prt, r_prt, tsReader.res_LOC0_prt->at(i));
0420         p2d_res_zr_prt[1]->Fill(z_prt, r_prt, tsReader.res_LOC1_prt->at(i));
0421         p2d_res_zr_prt[2]->Fill(z_prt, r_prt, tsReader.res_PHI_prt->at(i));
0422         p2d_res_zr_prt[3]->Fill(z_prt, r_prt, tsReader.res_THETA_prt->at(i));
0423         p2d_res_zr_prt[4]->Fill(z_prt, r_prt, tsReader.res_QOP_prt->at(i));
0424         p2d_res_zr_prt[5]->Fill(z_prt, r_prt, tsReader.res_T_prt->at(i));
0425         p2d_pull_zr_prt[0]->Fill(z_prt, r_prt, tsReader.pull_LOC0_prt->at(i));
0426         p2d_pull_zr_prt[1]->Fill(z_prt, r_prt, tsReader.pull_LOC1_prt->at(i));
0427         p2d_pull_zr_prt[2]->Fill(z_prt, r_prt, tsReader.pull_PHI_prt->at(i));
0428         p2d_pull_zr_prt[3]->Fill(z_prt, r_prt, tsReader.pull_THETA_prt->at(i));
0429         p2d_pull_zr_prt[4]->Fill(z_prt, r_prt, tsReader.pull_QOP_prt->at(i));
0430         p2d_pull_zr_prt[5]->Fill(z_prt, r_prt, tsReader.pull_T_prt->at(i));
0431       }
0432       if (filtered && tsReader.filtered->at(i)) {
0433         auto x_flt = tsReader.g_x_flt->at(i);
0434         auto y_flt = tsReader.g_y_flt->at(i);
0435         auto r_flt = std::hypot(x_flt, y_flt);
0436         auto z_flt = tsReader.g_z_flt->at(i);
0437         p2d_res_zr_flt[0]->Fill(z_flt, r_flt, tsReader.res_LOC0_flt->at(i));
0438         p2d_res_zr_flt[1]->Fill(z_flt, r_flt, tsReader.res_LOC1_flt->at(i));
0439         p2d_res_zr_flt[2]->Fill(z_flt, r_flt, tsReader.res_PHI_flt->at(i));
0440         p2d_res_zr_flt[3]->Fill(z_flt, r_flt, tsReader.res_THETA_flt->at(i));
0441         p2d_res_zr_flt[4]->Fill(z_flt, r_flt, tsReader.res_QOP_flt->at(i));
0442         p2d_res_zr_flt[5]->Fill(z_flt, r_flt, tsReader.res_T_flt->at(i));
0443         p2d_pull_zr_flt[0]->Fill(z_flt, r_flt, tsReader.pull_LOC0_flt->at(i));
0444         p2d_pull_zr_flt[1]->Fill(z_flt, r_flt, tsReader.pull_LOC1_flt->at(i));
0445         p2d_pull_zr_flt[2]->Fill(z_flt, r_flt, tsReader.pull_PHI_flt->at(i));
0446         p2d_pull_zr_flt[3]->Fill(z_flt, r_flt, tsReader.pull_THETA_flt->at(i));
0447         p2d_pull_zr_flt[4]->Fill(z_flt, r_flt, tsReader.pull_QOP_flt->at(i));
0448         p2d_pull_zr_flt[5]->Fill(z_flt, r_flt, tsReader.pull_T_flt->at(i));
0449       }
0450       if (smoothed && tsReader.smoothed->at(i)) {
0451         auto x_smt = tsReader.g_x_smt->at(i);
0452         auto y_smt = tsReader.g_y_smt->at(i);
0453         auto r_smt = std::hypot(x_smt, y_smt);
0454         auto z_smt = tsReader.g_z_smt->at(i);
0455         p2d_res_zr_smt[0]->Fill(z_smt, r_smt, tsReader.res_LOC0_smt->at(i));
0456         p2d_res_zr_smt[1]->Fill(z_smt, r_smt, tsReader.res_LOC1_smt->at(i));
0457         p2d_res_zr_smt[2]->Fill(z_smt, r_smt, tsReader.res_PHI_smt->at(i));
0458         p2d_res_zr_smt[3]->Fill(z_smt, r_smt, tsReader.res_THETA_smt->at(i));
0459         p2d_res_zr_smt[4]->Fill(z_smt, r_smt, tsReader.res_QOP_smt->at(i));
0460         p2d_res_zr_smt[5]->Fill(z_smt, r_smt, tsReader.res_T_smt->at(i));
0461         p2d_pull_zr_smt[0]->Fill(z_smt, r_smt, tsReader.pull_LOC0_smt->at(i));
0462         p2d_pull_zr_smt[1]->Fill(z_smt, r_smt, tsReader.pull_LOC1_smt->at(i));
0463         p2d_pull_zr_smt[2]->Fill(z_smt, r_smt, tsReader.pull_PHI_smt->at(i));
0464         p2d_pull_zr_smt[3]->Fill(z_smt, r_smt, tsReader.pull_THETA_smt->at(i));
0465         p2d_pull_zr_smt[4]->Fill(z_smt, r_smt, tsReader.pull_QOP_smt->at(i));
0466         p2d_pull_zr_smt[5]->Fill(z_smt, r_smt, tsReader.pull_T_smt->at(i));
0467       }
0468 
0469       int vol = tsReader.volume_id->at(i);
0470       int lay = tsReader.layer_id->at(i);
0471 
0472       /// Always fill (-1,-1), (vol, -1), (vol, lay)
0473       std::vector<std::array<int, 2>> fillIds = {
0474           {-1, -1}, {vol, -1}, {vol, lay}};
0475 
0476       for (const auto& fid : fillIds) {
0477         auto vlID = volLayIdCut(fid[0], fid[1])[0];
0478         // Fill predicated parameters
0479         if (predicted && tsReader.predicted->at(i)) {
0480           res_prt[vlID + paramNames[0]]->Fill(tsReader.res_LOC0_prt->at(i), 1);
0481           res_prt[vlID + paramNames[1]]->Fill(tsReader.res_LOC1_prt->at(i), 1);
0482           res_prt[vlID + paramNames[2]]->Fill(tsReader.res_PHI_prt->at(i), 1);
0483           res_prt[vlID + paramNames[3]]->Fill(tsReader.res_THETA_prt->at(i), 1);
0484           res_prt[vlID + paramNames[4]]->Fill(tsReader.res_QOP_prt->at(i), 1);
0485           res_prt[vlID + paramNames[5]]->Fill(tsReader.res_T_prt->at(i), 1);
0486           pull_prt[vlID + paramNames[0]]->Fill(tsReader.pull_LOC0_prt->at(i),
0487                                                1);
0488           pull_prt[vlID + paramNames[1]]->Fill(tsReader.pull_LOC1_prt->at(i),
0489                                                1);
0490           pull_prt[vlID + paramNames[2]]->Fill(tsReader.pull_PHI_prt->at(i), 1);
0491           pull_prt[vlID + paramNames[3]]->Fill(tsReader.pull_THETA_prt->at(i),
0492                                                1);
0493           pull_prt[vlID + paramNames[4]]->Fill(tsReader.pull_QOP_prt->at(i), 1);
0494           pull_prt[vlID + paramNames[5]]->Fill(tsReader.pull_T_prt->at(i), 1);
0495         }
0496         // Fill filtered parameters
0497         if (filtered && tsReader.filtered->at(i)) {
0498           res_flt[vlID + paramNames[0]]->Fill(tsReader.res_LOC0_flt->at(i), 1);
0499           res_flt[vlID + paramNames[1]]->Fill(tsReader.res_LOC1_flt->at(i), 1);
0500           res_flt[vlID + paramNames[2]]->Fill(tsReader.res_PHI_flt->at(i), 1);
0501           res_flt[vlID + paramNames[3]]->Fill(tsReader.res_THETA_flt->at(i), 1);
0502           res_flt[vlID + paramNames[4]]->Fill(tsReader.res_QOP_flt->at(i), 1);
0503           res_flt[vlID + paramNames[5]]->Fill(tsReader.res_T_flt->at(i), 1);
0504           pull_flt[vlID + paramNames[0]]->Fill(tsReader.pull_LOC0_flt->at(i),
0505                                                1);
0506           pull_flt[vlID + paramNames[1]]->Fill(tsReader.pull_LOC1_flt->at(i),
0507                                                1);
0508           pull_flt[vlID + paramNames[2]]->Fill(tsReader.pull_PHI_flt->at(i), 1);
0509           pull_flt[vlID + paramNames[3]]->Fill(tsReader.pull_THETA_flt->at(i),
0510                                                1);
0511           pull_flt[vlID + paramNames[4]]->Fill(tsReader.pull_QOP_flt->at(i), 1);
0512           pull_flt[vlID + paramNames[5]]->Fill(tsReader.pull_T_flt->at(i), 1);
0513         }
0514         // Fill smoothed parameters
0515         if (smoothed && tsReader.smoothed->at(i)) {
0516           res_smt[vlID + paramNames[0]]->Fill(tsReader.res_LOC0_smt->at(i), 1);
0517           res_smt[vlID + paramNames[1]]->Fill(tsReader.res_LOC1_smt->at(i), 1);
0518           res_smt[vlID + paramNames[2]]->Fill(tsReader.res_PHI_smt->at(i), 1);
0519           res_smt[vlID + paramNames[3]]->Fill(tsReader.res_THETA_smt->at(i), 1);
0520           res_smt[vlID + paramNames[4]]->Fill(tsReader.res_QOP_smt->at(i), 1);
0521           res_smt[vlID + paramNames[5]]->Fill(tsReader.res_T_smt->at(i), 1);
0522           pull_smt[vlID + paramNames[0]]->Fill(tsReader.pull_LOC0_smt->at(i),
0523                                                1);
0524           pull_smt[vlID + paramNames[1]]->Fill(tsReader.pull_LOC1_smt->at(i),
0525                                                1);
0526           pull_smt[vlID + paramNames[2]]->Fill(tsReader.pull_PHI_smt->at(i), 1);
0527           pull_smt[vlID + paramNames[3]]->Fill(tsReader.pull_THETA_smt->at(i),
0528                                                1);
0529           pull_smt[vlID + paramNames[4]]->Fill(tsReader.pull_QOP_smt->at(i), 1);
0530           pull_smt[vlID + paramNames[5]]->Fill(tsReader.pull_T_smt->at(i), 1);
0531         }
0532       }
0533     }
0534   }
0535 
0536   // Section 5: Histogram plotting
0537 
0538   // Plotting global profiles
0539   auto respull_mean_prf =
0540       new TCanvas("respull_mean_prf",
0541                   "Residual/Pull Distributions - mean profiles", 1800, 800);
0542   respull_mean_prf->Divide(3, 2);
0543 
0544   auto respull_var_prf =
0545       new TCanvas("respull_var_prf",
0546                   "Residual/Pull Distributions - variance profiles", 1800, 800);
0547   respull_var_prf->Divide(3, 2);
0548 
0549   auto plotProfiles = [&](std::array<TProfile2D*, 6>& profiles,
0550                           const std::string& res_pull,
0551                           const std::string& type) -> void {
0552     // Mean
0553     for (std::size_t ipar = 0; ipar < paramNames.size(); ++ipar) {
0554       respull_mean_prf->cd(ipar + 1);
0555 
0556       if (res_pull == "pull") {
0557         profiles[ipar]->GetZaxis()->SetRangeUser(-1. * pullRange, pullRange);
0558       }
0559       adaptColorPalette(profiles[ipar], -1. * pullRange, pullRange, 0., 0.25,
0560                         104);
0561       profiles[ipar]->Draw("colz");
0562       profiles[ipar]->Write();
0563     }
0564     // Save the canvas: mean
0565     if (!saveAs.empty()) {
0566       respull_mean_prf->SaveAs((std::string("all_") + res_pull +
0567                                 std::string("_mean_prf_") + type +
0568                                 std::string(".") + saveAs)
0569                                    .c_str());
0570     }
0571 
0572     // Variance
0573     for (std::size_t ipar = 0; ipar < paramNames.size(); ++ipar) {
0574       respull_var_prf->cd(ipar + 1);
0575       auto zAxis = profiles[ipar]->GetXaxis();
0576       auto rAxis = profiles[ipar]->GetYaxis();
0577       auto binsZ = zAxis->GetNbins();
0578       auto binsR = rAxis->GetNbins();
0579       std::string hist_name = "all_";
0580       hist_name += res_pull;
0581       hist_name += "_";
0582       hist_name += paramNames[ipar];
0583       hist_name += "_";
0584       hist_name += type;
0585       TH2F* var =
0586           new TH2F(hist_name.c_str(),
0587                    (res_pull + std::string(" ") + paramNames[ipar]).c_str(),
0588                    binsZ, zAxis->GetXmin(), zAxis->GetXmax(), binsR,
0589                    rAxis->GetXmin(), rAxis->GetXmax());
0590       for (int iz = 1; iz <= binsZ; ++iz) {
0591         for (int ir = 1; ir <= binsR; ++ir) {
0592           if (profiles[ipar]->GetBinContent(iz, ir) != 0.) {
0593             var->SetBinContent(iz, ir, profiles[ipar]->GetBinError(iz, ir));
0594           }
0595         }
0596       }
0597       if (res_pull == "pull") {
0598         var->GetZaxis()->SetRangeUser(0., pullRange);
0599       }
0600       adaptColorPalette(var, 0., pullRange, 1., 0.5, 104);
0601       var->Draw("colz");
0602       var->Write();
0603     }
0604     // Save the canvas: pulls
0605     if (!saveAs.empty()) {
0606       respull_var_prf->SaveAs((std::string("all_") + res_pull +
0607                                std::string("_var_prf_") + type +
0608                                std::string(".") + saveAs)
0609                                   .c_str());
0610     }
0611   };
0612 
0613   // Plotting profiles: res/pulls
0614   if (predicted) {
0615     plotProfiles(p2d_res_zr_prt, "res", "prt");
0616     plotProfiles(p2d_pull_zr_prt, "pull", "prt");
0617   }
0618   if (filtered) {
0619     plotProfiles(p2d_res_zr_flt, "res", "flt");
0620     plotProfiles(p2d_pull_zr_flt, "pull", "flt");
0621   }
0622   if (smoothed) {
0623     plotProfiles(p2d_res_zr_smt, "res", "smt");
0624     plotProfiles(p2d_pull_zr_smt, "pull", "smt");
0625   }
0626 
0627   // Plotting residual
0628   auto residuals =
0629       new TCanvas("residuals", "Residual Distributions", 1200, 800);
0630   residuals->Divide(3, 2);
0631 
0632   auto pulls = new TCanvas("pulls", "Pull distributions", 1200, 800);
0633   pulls->Divide(3, 2);
0634 
0635   for (const auto& [vol, layers] : volLayIds) {
0636     for (const auto& lay : layers) {
0637       auto vlID = volLayIdCut(vol, lay)[0];
0638 
0639       // Residual plotting
0640       for (std::size_t ipar = 0; ipar < paramNames.size(); ipar++) {
0641         residuals->cd(ipar + 1);
0642 
0643         auto legend = new TLegend(0.7, 0.7, 0.9, 0.9);
0644 
0645         const std::string name = vlID + paramNames.at(ipar);
0646 
0647         // Draw them
0648         if (smoothed) {
0649           res_smt[name]->DrawNormalized("");
0650           res_smt[name]->Write();
0651           legend->AddEntry(res_smt[name], "smoothed", "l");
0652         }
0653         if (filtered) {
0654           std::string drawOptions = smoothed ? "same" : "";
0655           res_flt[name]->DrawNormalized(drawOptions.c_str());
0656           res_flt[name]->Write();
0657           legend->AddEntry(res_flt[name], "filtered", "l");
0658         }
0659         if (predicted) {
0660           std::string drawOptions = (smoothed || filtered) ? "same" : "";
0661           res_prt[name]->DrawNormalized(drawOptions.c_str());
0662           res_prt[name]->Write();
0663           legend->AddEntry(res_prt[name], "predicted", "l");
0664         }
0665 
0666         legend->SetBorderSize(0);
0667         legend->SetFillStyle(0);
0668         legend->SetTextFont(42);
0669         legend->Draw();
0670       }
0671       if (!saveAs.empty()) {
0672         residuals->SaveAs((vlID + std::string("residuals.") + saveAs).c_str());
0673       }
0674 
0675       // Pull plotting & writing
0676       for (std::size_t ipar = 0; ipar < paramNames.size(); ipar++) {
0677         const std::string name = vlID + paramNames.at(ipar);
0678 
0679         pulls->cd(ipar + 1);
0680 
0681         auto legend = new TLegend(0.7, 0.5, 0.9, 0.9);
0682 
0683         // Fit the smoothed distribution
0684         if (smoothed) {
0685           auto drawOptions = "pe";
0686 
0687           auto scale = 1. / pull_smt[name]->Integral("width");
0688           pull_smt[name]->Scale(scale);
0689           pull_smt[name]->Draw(drawOptions);
0690           pull_smt[name]->Write();
0691 
0692           legend->AddEntry(pull_smt[name], "smoothed", "pe");
0693 
0694           if (fitSmoothed) {
0695             pull_smt[name]->Fit("gaus", "q");
0696             TF1* gauss = pull_smt[name]->GetFunction("gaus");
0697             gauss->SetLineColorAlpha(kBlack, 0.5);
0698 
0699             auto mu = gauss->GetParameter(1);
0700             auto sigma = gauss->GetParameter(2);
0701             auto mu_fit_info = TString::Format("#mu = %.3f", mu);
0702             auto su_fit_info = TString::Format("#sigma = %.3f", sigma);
0703 
0704             legend->AddEntry(gauss, mu_fit_info.Data(), "l");
0705             legend->AddEntry(pull_prt[name], su_fit_info.Data(), "");
0706           }
0707         }
0708 
0709         if (filtered) {
0710           auto drawOptions = smoothed ? "pe same" : "pe";
0711 
0712           auto scale = 1. / pull_flt[name]->Integral("width");
0713           pull_flt[name]->Scale(scale);
0714           pull_flt[name]->Draw(drawOptions);
0715           pull_flt[name]->Write();
0716 
0717           legend->AddEntry(pull_flt[name], "filtered", "pe");
0718 
0719           if (fitFiltered) {
0720             pull_flt[name]->Fit("gaus", "q", "same");
0721             TF1* gauss = pull_flt[name]->GetFunction("gaus");
0722             gauss->SetLineColorAlpha(kBlue, 0.5);
0723 
0724             auto mu = gauss->GetParameter(1);
0725             auto sigma = gauss->GetParameter(2);
0726             auto mu_fit_info = TString::Format("#mu = %.3f", mu);
0727             auto su_fit_info = TString::Format("#sigma = %.3f", sigma);
0728 
0729             legend->AddEntry(gauss, mu_fit_info.Data(), "l");
0730             legend->AddEntry(pull_prt[name], su_fit_info.Data(), "");
0731           }
0732         }
0733 
0734         if (predicted) {
0735           auto drawOptions = (smoothed || filtered) ? "pe same" : "pe";
0736 
0737           auto scale = 1. / pull_prt[name]->Integral("width");
0738           pull_prt[name]->Scale(scale);
0739           pull_prt[name]->Draw(drawOptions);
0740           pull_prt[name]->Write();
0741 
0742           legend->AddEntry(pull_prt[name], "predicted", "pe");
0743 
0744           if (fitPredicted) {
0745             pull_prt[name]->Fit("gaus", "q", "same");
0746             TF1* gauss = pull_prt[name]->GetFunction("gaus");
0747             gauss->SetLineColorAlpha(kRed, 0.5);
0748 
0749             auto mu = gauss->GetParameter(1);
0750             auto sigma = gauss->GetParameter(2);
0751             auto mu_fit_info = TString::Format("#mu = %.3f", mu);
0752             auto su_fit_info = TString::Format("#sigma = %.3f", sigma);
0753 
0754             legend->AddEntry(gauss, mu_fit_info.Data(), "l");
0755             legend->AddEntry(pull_prt[name], su_fit_info.Data(), "");
0756           }
0757         }
0758 
0759         // Reference standard normal pdf
0760         auto ref = new TF1("ref", "ROOT::Math::normal_pdf(x,1,0)", -5, 5);
0761         ref->SetLineColor(kGreen);
0762         ref->Draw("same");
0763         legend->AddEntry(ref, "#mu = 0 #sigma = 1", "l");
0764 
0765         legend->SetBorderSize(0);
0766         legend->SetFillStyle(0);
0767         legend->SetTextFont(42);
0768         legend->Draw();
0769       }
0770 
0771       // Save the Canvases as pictures
0772       if (!saveAs.empty()) {
0773         pulls->SaveAs((vlID + std::string("pulls.") + saveAs).c_str());
0774       }
0775     }
0776   }
0777 
0778   if (out != nullptr) {
0779     out->Close();
0780   }
0781   return 0;
0782 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!