analysis/NNScaleMap/plotInjectionScaleResiduals.C

0001 #include <TCanvas.h>
0002 #include <TFile.h>
0003 #include <TGaxis.h>
0004 #include <TH2.h>
0005 #include <TLatex.h>
0006 #include <TObject.h>
0007 #include <TPad.h>
0008 #include <TStyle.h>
0009 #include <TSystem.h>
0010
0011 #include <algorithm>
0012 #include <cmath>
0013 #include <iostream>
0014 #include <memory>
0015 #include <string>
0016 #include <vector>
0017
0018 void drawResidualPadLabel(const double pt, const int charge, const bool grid)
0019 {
0020     TLatex label;
0021     label.SetNDC();
0022     label.SetTextFont(42);
0023     label.SetTextSize(grid ? 0.045 : 0.035);
0024
0025     const double x = gPad->GetLeftMargin();
0026     const double y = 1 - gPad->GetTopMargin() + 0.025;
0027     label.DrawLatex(x, y, Form("p_{T} = %.2f GeV, 1/p_{T} = %.2f GeV^{-1}, q = %+d", pt, 1.0 / pt, charge));
0028 }
0029
0030 void plotInjectionScaleResiduals(std::string injectionClosureFile = "/home/hjheng/Documents/sphenix/NN-Momentum-Calibration/Injection-Pythia/closure_diagnostics-kshort/kshort_closure.root",
0031     std::string learnedMapFile = "/home/hjheng/Documents/sphenix/NN-Momentum-Calibration/NN-training/calib_out_pythiaInjection_20260805/kappa_map.root", std::string outdir = "")
0032 {
0033     const std::vector<double> plotPt = {0.3, 0.5, 1.0, 2.0, 3.0};
0034     const std::vector<std::string> ptTag = {"pt0p30", "pt0p50", "pt1p00", "pt2p00", "pt3p00"};
0035     const int charge[2] = {+1, -1};
0036     const int injectionChargeIndex[2] = {0, 1};
0037     const std::string learnedChargeName[2] = {"qplus", "qminus"};
0038
0039     if (gSystem->AccessPathName(learnedMapFile.c_str(), kReadPermission))
0040     {
0041         const std::string singularName = "kappa_map.root";
0042         const std::size_t pos = learnedMapFile.rfind(singularName);
0043         if (pos != std::string::npos && pos + singularName.size() == learnedMapFile.size())
0044         {
0045             std::string fallback = learnedMapFile;
0046             fallback.replace(pos, singularName.size(), "kappa_maps.root");
0047             if (!gSystem->AccessPathName(fallback.c_str(), kReadPermission))
0048             {
0049                 std::cout << "Requested learned map " << learnedMapFile << " was not found; using " << fallback << " instead." << std::endl;
0050                 learnedMapFile = fallback;
0051             }
0052         }
0053     }
0054
0055     if (outdir.empty())
0056     {
0057         std::string learnedDir = ".";
0058         const std::size_t slash = learnedMapFile.find_last_of('/');
0059         if (slash != std::string::npos)
0060             learnedDir = learnedMapFile.substr(0, slash);
0061         outdir = learnedDir + "/injection_scale_residuals";
0062     }
0063     gSystem->mkdir(outdir.c_str(), true);
0064
0065     TGaxis::SetMaxDigits(3);
0066     gStyle->SetOptStat(0);
0067     gStyle->SetPalette(kLightTemperature);
0068
0069     std::unique_ptr<TFile> injection(TFile::Open(injectionClosureFile.c_str(), "READ"));
0070     if (!injection || injection->IsZombie())
0071     {
0072         std::cerr << "Could not open injection closure file: " << injectionClosureFile << std::endl;
0073         return;
0074     }
0075
0076     std::unique_ptr<TFile> learned(TFile::Open(learnedMapFile.c_str(), "READ"));
0077     if (!learned || learned->IsZombie())
0078     {
0079         std::cerr << "Could not open learned scale-map file: " << learnedMapFile << std::endl;
0080         return;
0081     }
0082
0083     const std::string residualRootFile = outdir + "/injection_scale_residuals.root";
0084     std::unique_ptr<TFile> output(TFile::Open(residualRootFile.c_str(), "RECREATE"));
0085     if (!output || output->IsZombie())
0086     {
0087         std::cerr << "Could not create residual output file: " << residualRootFile << std::endl;
0088         return;
0089     }
0090
0091     std::vector<TH2 *> residuals;
0092     std::vector<int> residualCharge;
0093     std::vector<double> residualPt;
0094     std::vector<std::string> residualStem;
0095     double maxAbsResidual = 0.0;
0096
0097     for (std::size_t ipt = 0; ipt < plotPt.size(); ++ipt)
0098     {
0099         for (int icharge = 0; icharge < 2; ++icharge)
0100         {
0101             const std::string injectionKey = Form("diagnostics/h_pt_scale_q%d_pt%d", injectionChargeIndex[icharge], static_cast<int>(ipt));
0102             const std::string learnedKey = "kappa_" + learnedChargeName[icharge] + "_" + ptTag[ipt];
0103             TH2 *injectedPtScale = dynamic_cast<TH2 *>(injection->Get(injectionKey.c_str()));
0104             TH2 *learnedPtScale = dynamic_cast<TH2 *>(learned->Get(learnedKey.c_str()));
0105             if (!injectedPtScale || !learnedPtScale)
0106             {
0107                 std::cerr << "Missing histogram pair: " << injectionKey << " and " << learnedKey << std::endl;
0108                 return;
0109             }
0110
0111             if (injectedPtScale->GetNbinsX() != learnedPtScale->GetNbinsX() || injectedPtScale->GetNbinsY() != learnedPtScale->GetNbinsY() ||
0112                 std::abs(injectedPtScale->GetXaxis()->GetXmin() - learnedPtScale->GetXaxis()->GetXmin()) > 1.0e-9 ||
0113                 std::abs(injectedPtScale->GetXaxis()->GetXmax() - learnedPtScale->GetXaxis()->GetXmax()) > 1.0e-9 ||
0114                 std::abs(injectedPtScale->GetYaxis()->GetXmin() - learnedPtScale->GetYaxis()->GetXmin()) > 1.0e-9 ||
0115                 std::abs(injectedPtScale->GetYaxis()->GetXmax() - learnedPtScale->GetYaxis()->GetXmax()) > 1.0e-9)
0116             {
0117                 std::cerr << "Injection and learned scale maps do not have matching eta-phi binning for " << learnedKey << std::endl;
0118                 return;
0119             }
0120
0121             const std::string stem = "residual_" + learnedChargeName[icharge] + "_" + ptTag[ipt];
0122             TH2 *residual = static_cast<TH2 *>(learnedPtScale->Clone(stem.c_str()));
0123             residual->SetDirectory(nullptr);
0124             residual->Reset("ICES");
0125             residual->SetStats(false);
0126             residual->SetTitle(Form("Injected #times learned p_{T} scale residual, q = %+d, p_{T} = %.2f GeV;#eta;#phi;injected scale #times learned scale - 1", charge[icharge], plotPt[ipt]));
0127
0128             for (int ix = 1; ix <= residual->GetNbinsX(); ++ix)
0129             {
0130                 for (int iy = 1; iy <= residual->GetNbinsY(); ++iy)
0131                 {
0132                     const double injectedMultiplier = 1.0 + injectedPtScale->GetBinContent(ix, iy) / 100.0;
0133                     const double value = injectedMultiplier * learnedPtScale->GetBinContent(ix, iy) - 1.0;
0134                     residual->SetBinContent(ix, iy, value);
0135                     maxAbsResidual = std::max(maxAbsResidual, std::abs(value));
0136                 }
0137             }
0138
0139             output->cd();
0140             residual->Write("", TObject::kOverwrite);
0141             residuals.push_back(residual);
0142             residualCharge.push_back(charge[icharge]);
0143             residualPt.push_back(plotPt[ipt]);
0144             residualStem.push_back(stem);
0145         }
0146     }
0147
0148     if (maxAbsResidual <= 0.0)
0149         maxAbsResidual = 1.0e-12;
0150     // const double zmax = (maxAbsResidual > 1.0e-2) ? 1.05 * maxAbsResidual : 1.0e-2;
0151     const double zmax = maxAbsResidual;
0152
0153     for (std::size_t i = 0; i < residuals.size(); ++i)
0154     {
0155         TCanvas canvas(Form("c_%s", residualStem[i].c_str()), residuals[i]->GetTitle(), 850, 740);
0156         canvas.SetTicks(1, 1);
0157         canvas.SetLeftMargin(0.15);
0158         canvas.SetRightMargin(0.22);
0159         // canvas.SetBottomMargin(0.13);
0160         canvas.SetTopMargin(0.07);
0161
0162         residuals[i]->SetMinimum(-zmax);
0163         residuals[i]->SetMaximum(+zmax);
0164         residuals[i]->SetContour(255);
0165         residuals[i]->GetZaxis()->SetTitle("scale residual");
0166         residuals[i]->GetZaxis()->SetTitleOffset(1.6);
0167         residuals[i]->Draw("COLZ");
0168         drawResidualPadLabel(residualPt[i], residualCharge[i], false);
0169
0170         canvas.SaveAs((outdir + "/" + residualStem[i] + ".png").c_str());
0171         canvas.SaveAs((outdir + "/" + residualStem[i] + ".pdf").c_str());
0172     }
0173
0174     TCanvas grid("c_injection_scale_residuals", "Injection scale residuals", 360 * plotPt.size(), 680);
0175     grid.Divide(static_cast<int>(plotPt.size()), 2, 0.001, 0.001);
0176     for (std::size_t i = 0; i < residuals.size(); ++i)
0177     {
0178         int ptIndex = 0;
0179         for (std::size_t ipt = 0; ipt < plotPt.size(); ++ipt)
0180         {
0181             if (std::abs(plotPt[ipt] - residualPt[i]) < 1.0e-9)
0182             {
0183                 ptIndex = static_cast<int>(ipt);
0184                 break;
0185             }
0186         }
0187
0188         const int row = residualCharge[i] > 0 ? 0 : 1;
0189         grid.cd(row * static_cast<int>(plotPt.size()) + ptIndex + 1);
0190         gPad->SetTicks(1, 1);
0191         gPad->SetLeftMargin(0.13);
0192         gPad->SetRightMargin(0.18);
0193         // gPad->SetTopMargin(0.12);
0194         gPad->SetBottomMargin(0.12);
0195
0196         residuals[i]->SetMinimum(-zmax);
0197         residuals[i]->SetMaximum(+zmax);
0198         residuals[i]->GetZaxis()->SetTitleOffset(1.5);
0199         residuals[i]->Draw("COLZ");
0200         drawResidualPadLabel(residualPt[i], residualCharge[i], true);
0201     }
0202
0203     grid.SaveAs((outdir + "/all_injection_scale_residuals.png").c_str());
0204     grid.SaveAs((outdir + "/all_injection_scale_residuals.pdf").c_str());
0205     output->cd();
0206     grid.Write("", TObject::kOverwrite);
0207     output->Close();
0208
0209     std::cout << "Saved " << residuals.size() << " residual maps to " << outdir << std::endl;
0210     std::cout << "Residual ROOT file: " << residualRootFile << std::endl;
0211     std::cout << "Compared (1 + injected p_{T} scale[%]/100) * learned p_{T} scale - 1." << std::endl;
0212
0213     for (TH2 *residual : residuals)
0214         delete residual;
0215 }