Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-05-23 08:12:15 

0001 #include "analysisHelper.h"
0002 #include "RooUnfoldResponse.h"
0003 
0004 // ─────────────────────────────────────────────────────────────────────────────
0005 //  drawClosure.C
0006 //
0007 //  Draws diagnostic plots from the outputs of:
0008 //    jet_pT_doUnfolding.C  → jetPTFile
0009 //    wEEC_doUnfolding.C    → wEECFile
0010 //    merged response file  → respFile  (truth reference)
0011 //
0012 //  Truth reference for wEEC per-bin plots:
0013 //    kFull — hWEEC3D_truth_{k} from respFile (fullClosure response)
0014 //    kHalf — hWEEC3D_meas_truth_{k} from respFile (odd-half truth = correct closure target)
0015 //    kData — hWEEC3D_truth_{k} from respFile (fullClosure response, full sim)
0016 //
0017 //  Dijet pT plots:
0018 //    jetPT-{label}.png              — 2D truth/unfolded/ratio  (kFull/kHalf)
0019 //    jetPT_projections-{label}.png  — 1D lead/subl pT          (kFull/kHalf)
0020 //    jetPT_recoVsUnfolded-{label}   — reco vs unfolded          (kData)
0021 //
0022 //  wEEC plots:
0023 //    wEEC_inclusive-{label}.png     — inclusive unfolded vs truth
0024 //    wEEC_perBin-{label}.png        — grid: unfolded vs truth per dijet pT bin
0025 //    wEEC_ratio-{label}.png         — grid: unfolded/truth ratio per dijet pT bin
0026 //    wEEC_measVsUnf-{label}.png     — grid: measured vs unfolded (vs truth)
0027 //                                     per dijet pT bin
0028 //
0029 //  All wEEC histograms normalized by NormalizeWEEC (integral then bin width).
0030 // ─────────────────────────────────────────────────────────────────────────────
0031 
0032 static const int kColTruth    = kRed   + 1;
0033 static const int kColUnfolded = kBlue  + 1;
0034 static const int kColMeas     = kGreen + 2;
0035 static const int kMarTruth    = 21;
0036 static const int kMarUnfolded = 20;
0037 static const int kMarMeas     = 33;
0038 
0039 std::pair<TPad*,TPad*> MakeRatioPads(TVirtualPad* parent)
0040 {
0041     parent->cd();
0042     TPad* pTop = new TPad("pTop","",0,0.30,1,1.0);
0043     TPad* pBot = new TPad("pBot","",0,0.00,1,0.30);
0044     pTop->SetBottomMargin(0.03); pBot->SetTopMargin(0.03);
0045     pBot->SetBottomMargin(0.35);
0046     pTop->Draw(); pBot->Draw();
0047     return {pTop, pBot};
0048 }
0049 
0050 void StyleRatio(TH1D* h, const char* xTitle,
0051                 const char* yTitle = "Unf/Truth",
0052                 double lo = 0.8, double hi = 1.2)
0053 {
0054     h->SetLineColor(kBlack); h->SetMarkerColor(kBlack); h->SetMarkerStyle(20);
0055     h->GetYaxis()->SetRangeUser(lo,hi);
0056     h->GetYaxis()->SetTitle(yTitle);
0057     h->GetYaxis()->SetTitleSize(0.12); h->GetYaxis()->SetTitleOffset(0.38);
0058     h->GetYaxis()->SetLabelSize(0.10); h->GetYaxis()->SetNdivisions(504);
0059     h->GetXaxis()->SetTitle(xTitle);
0060     h->GetXaxis()->SetTitleSize(0.12); h->GetXaxis()->SetLabelSize(0.10);
0061     h->SetTitle("");
0062 }
0063 
0064 void DrawRefLine(TH1D* h, double y = 1.0)
0065 {
0066     TLine* l = new TLine(h->GetXaxis()->GetXmin(),y,h->GetXaxis()->GetXmax(),y);
0067     l->SetLineStyle(2); l->SetLineColor(kGray+1); l->Draw();
0068 }
0069 
0070 void StyleLine(TH1D* h, int col, int mar)
0071 { h->SetLineColor(col); h->SetMarkerColor(col); h->SetMarkerStyle(mar); }
0072 
0073 TH1D* MakeRatioHist(TH1D* hNum, TH1D* hDen, const char* name)
0074 {
0075     TH1D* hR = (TH1D*)hNum->Clone(name);
0076     hR->SetDirectory(0);
0077     for (int k=1; k<=hR->GetNbinsX(); ++k) {
0078         double n=hNum->GetBinContent(k), ne=hNum->GetBinError(k);
0079         double d=hDen->GetBinContent(k), de=hDen->GetBinError(k);
0080         if (d>0) {
0081             double r=n/d;
0082             hR->SetBinContent(k,r);
0083             hR->SetBinError(k, r*std::sqrt((n>0?(ne/n)*(ne/n):0)+(de/d)*(de/d)));
0084         } else { hR->SetBinContent(k,0); hR->SetBinError(k,0); }
0085     }
0086     return hR;
0087 }
0088 
0089 // ─────────────────────────────────────────────────────────────────────────────
0090 void drawClosure(const char* outDir      = ".",
0091                  const char* respFile    = nullptr,
0092                  Mode        mode        = Mode::kFull)
0093 {
0094     gROOT->SetBatch(true);
0095     gStyle->SetOptStat(0); gStyle->SetOptTitle(0);
0096     gStyle->SetPaintTextFormat(".2f");
0097 
0098     const bool doTrim    = true;
0099     const bool usePWMean = true;   // true  → luminosity-weighted mean pair weight per bin
0100                                    // false → arithmetic bin-center approximation
0101     const bool projAccessibleOnly = false;  // true  → skip truth bins below reco threshold
0102     const bool       useDirectCov = true; // true  → covDirectWEEC3D_{k} (RooUnfold Errunfold, preferred)
0103                                         // false → covWEEC3D_{k}       (hand-built M*Vmeas*M^T)
0104 
0105     // ── Analysis-quality bin mask ─────────────────────────────────────────
0106     // Bins listed here are grayed out in the *_masked variants of the grid
0107     // plots.  Indices are 0-based (iL, iS).  Edit this list to change the
0108     // mask without touching anything else.
0109     const std::vector<std::pair<int,int>> maskedBins = {
0110         {3, 1},  // 40-50 GeV lead / 15-20 GeV subl
0111         {4, 1},  // 50-60 GeV lead / 15-20 GeV subl
0112         {5, 1},  // 60-80 GeV lead / 15-20 GeV subl
0113         {4, 2},  // 50-60 GeV lead / 20-30 GeV subl
0114         {5, 2},  // 60-80 GeV lead / 20-30 GeV subl
0115         {5, 3},  // (reserved — update if bin structure changes)
0116     };
0117                                            //          in wEEC_projections (cleaner closure)
0118                                            // false → include all truth bins (shows full
0119                                            //          unfolding scope incl. miss-only bins)
0120     const bool isData    = (mode == Mode::kData);
0121     const bool isHalf    = (mode == Mode::kHalf);
0122     const bool isClosure = !isData;
0123 
0124     // Name of the per-ΔΦ truth histogram to use as the closure reference:
0125     //   kFull/kData → hWEEC3D_truth_{k}       (training-half or full-sim truth)
0126     //   kHalf       → hWEEC3D_meas_truth_{k}  (odd-half truth = correct target)
0127     auto truth3DName = [&](int k) -> std::string {
0128         return isHalf
0129             ? std::format("hWEEC3D_meas_truth_{}", k)
0130             : std::format("hWEEC3D_truth_{}", k);
0131     };
0132     const std::string label   = ModeLabel(mode);
0133     const std::string plotDir = std::format("{}/Plots", outDir);
0134 
0135     std::string wEECFile  = std::format("{}/wEEC-{}.root",   outDir, label);
0136 
0137     // ════════════════════════════════════════════════════════════════════════
0138     //  Dijet pT plots
0139     //  Jet pT histograms are written by wEEC_doUnfolding.C into the same
0140     //  wEEC output file, so only one input file is needed.
0141     // ════════════════════════════════════════════════════════════════════════
0142     TFile* fWEEC = TFile::Open(wEECFile.c_str(),"READ");
0143     if (!fWEEC || fWEEC->IsZombie()) { std::cerr<<"Cannot open "<<wEECFile<<"\n"; return; }
0144 
0145     // Open response file immediately — needed for truth reference throughout
0146     TFile* fResp = nullptr;
0147     if (respFile) {
0148         fResp = TFile::Open(respFile,"READ");
0149         if (!fResp || fResp->IsZombie()) {
0150             std::cerr<<"Cannot open resp file "<<respFile<<"\n"; fResp=nullptr;
0151         }
0152     }
0153 
0154     // ── Load pair weight mean histograms (used when usePWMean == true) ────
0155     // hWEEC3D_pwNum_{k} = sum(evtWeight * pairWeight) per truth flat3D bin
0156     // hWEEC3D_pwDen_{k} = sum(evtWeight)              per truth flat3D bin
0157     // Both are nullptr if usePWMean is false or the histograms are absent.
0158     std::vector<TH1D*> hPWNum(nDphi, nullptr), hPWDen(nDphi, nullptr);
0159     if (usePWMean && fResp) {
0160         for (int k = 0; k < nDphi; ++k) {
0161             hPWNum[k] = (TH1D*)fResp->Get(std::format("hWEEC3D_pwNum_{}", k).c_str());
0162             hPWDen[k] = (TH1D*)fResp->Get(std::format("hWEEC3D_pwDen_{}", k).c_str());
0163             if (hPWNum[k]) hPWNum[k]->SetDirectory(0);
0164             if (hPWDen[k]) hPWDen[k]->SetDirectory(0);
0165         }
0166     }
0167 
0168     // ── Load unfolding covariance matrices ────────────────────────────────
0169     // Two covariance options written by wEEC_doUnfolding.C:
0170     //   covDirectWEEC3D_{k} — directly from RooUnfoldBayes::Errunfold()
0171     //                         (kCovariance mode); full off-diagonal propagation
0172     //                         through the Bayesian iterations.  Preferred.
0173     //   covWEEC3D_{k}       — hand-built analytic M*Vmeas*M^T; correct only
0174     //                         in the linear (low-iteration) limit and ignores
0175     //                         correlations introduced by the Bayesian updates.
0176     // Selected via useDirectCov (default: true → covDirectWEEC3D).
0177     // nullptr entries mean the covariance is unavailable — projections fall
0178     // back to naive diagonal (GetBinError) propagation automatically.
0179     std::vector<TMatrixD*> hCov(nDphi, nullptr);
0180     {
0181         const std::string covKeyFmt = useDirectCov
0182             ? "covDirectWEEC3D_{}"
0183             : "covWEEC3D_{}";
0184         const std::string covLabel  = useDirectCov
0185             ? "covDirectWEEC3D (RooUnfold Errunfold)"
0186             : "covWEEC3D (M*Vmeas*M^T)";
0187         int nLoaded = 0;
0188         for (int k = 0; k < nDphi; ++k) {
0189             hCov[k] = (TMatrixD*)fWEEC->Get(
0190                 std::format(useDirectCov
0191             ? "covDirectWEEC3D_{}"
0192             : "covWEEC3D_{}", k).c_str());
0193             if (hCov[k]) ++nLoaded;
0194         }
0195         std::cout << std::format("Loaded {}/{} covariance matrices [{}] from {}\n",
0196             nLoaded, nDphi, covLabel, wEECFile);
0197         if (nLoaded > 0 && hCov[0])
0198             std::cout << std::format("  cov[0] dimensions: {}×{}\n",
0199                 hCov[0]->GetNrows(), hCov[0]->GetNcols());
0200     }
0201 
0202     // ════════════════════════════════════════════════════════════════════════
0203     //  Project hWEEC3D_unfolded_{k} → per-dijet-pT-bin and inclusive wEEC
0204     //  This replaces loading pre-projected hWEEC_bin_{ft} / hWEEC_inclusive
0205     //  which are no longer written by wEEC_doUnfolding.C.
0206     // ════════════════════════════════════════════════════════════════════════
0207     std::vector<TH1D*> hUnfBin(nTrueFlat(), nullptr);
0208     for (int ft=0; ft<nTrueFlat(); ++ft) {
0209         int iL,iS; TrueFlatToIJ(ft,iL,iS);
0210         hUnfBin[ft] = new TH1D(
0211             std::format("hUnfBin_{}",ft).c_str(),
0212             std::format("Unfolded iL{} iS{};#Delta#phi;wEEC",iL,iS).c_str(),
0213             nDphi, dPhiBins.data());
0214         hUnfBin[ft]->Sumw2();
0215         hUnfBin[ft]->SetDirectory(0);
0216     }
0217     TH1D* hInclUnfolded = new TH1D("hWEEC_inclusive",";#Delta#phi;wEEC",nDphi,dPhiBins.data());
0218     hInclUnfolded->Sumw2();
0219     hInclUnfolded->SetDirectory(0);
0220 
0221     for (int k=0; k<nDphi; ++k) {
0222         TH1D* hU3D = (TH1D*)fWEEC->Get(std::format("hWEEC3D_unfolded_{}",k).c_str());
0223         if (!hU3D) continue;
0224         // Inclusive projection — IsRecoAccessible gate is inside ProjectWEEC3DSlice
0225         double ival, ierr;
0226         if (usePWMean && hPWNum[k] && hPWDen[k])
0227             ProjectWEEC3DSlice(hU3D, ival, ierr, hPWNum[k], hPWDen[k], hCov[k]);
0228         else
0229             ProjectWEEC3DSlice(hU3D, ival, ierr, hCov[k]);
0230         hInclUnfolded->SetBinContent(k+1, ival);
0231         hInclUnfolded->SetBinError  (k+1, ierr);
0232         // Per-bin projection — skip bins with no reco support
0233         for (int ft=0; ft<nTrueFlat(); ++ft) {
0234             int iL,iS; TrueFlatToIJ(ft,iL,iS);
0235             if (!IsRecoAccessible(iL,iS)) continue;
0236             double val,err;
0237             if (usePWMean && hPWNum[k] && hPWDen[k])
0238                 ProjectWEEC3DSliceForBin(hU3D, iL, iS, val, err, hPWNum[k], hPWDen[k], hCov[k]);
0239             else
0240                 ProjectWEEC3DSliceForBin(hU3D, iL, iS, val, err, hCov[k]);
0241             hUnfBin[ft]->SetBinContent(k+1, val);
0242             hUnfBin[ft]->SetBinError  (k+1, err);
0243         }
0244     }
0245     NormalizeWEEC(hInclUnfolded);
0246     for (int ft=0; ft<nTrueFlat(); ++ft) NormalizeWEEC(hUnfBin[ft]);
0247 
0248     TH1D* hJetPt_unfolded = (TH1D*)fWEEC->Get("hJetPt_unfolded");
0249     TH1D* hJetPt_meas     = (TH1D*)fWEEC->Get("hJetPt_meas");
0250     TH1D* hJetPt_truth    = isClosure ? (TH1D*)fWEEC->Get("hJetPt_truth") : nullptr;
0251     if (!hJetPt_unfolded || !hJetPt_meas) {
0252         std::cerr<<"Missing jet pT histograms in "<<wEECFile<<"\n"; return;
0253     }
0254     hJetPt_unfolded->SetDirectory(0); hJetPt_meas->SetDirectory(0);
0255     if (hJetPt_truth) hJetPt_truth->SetDirectory(0);
0256 
0257     // Unfolded and truth live in truth-index space; meas lives in reco-index space
0258     TH2D* h2Unf  = new TH2D("h2Unf", ";Lead p_{T} (GeV);Subl p_{T} (GeV)",
0259                               nTrueLead,trueLeadPtBins.data(),nTrueSubl,trueSublPtBins.data());
0260     TH2D* h2Meas = new TH2D("h2Meas",";Lead p_{T} (GeV);Subl p_{T} (GeV)",
0261                               nRecoLead,recoLeadPtBins.data(),nRecoSubl,recoSublPtBins.data());
0262     TH2D* h2Truth = isClosure
0263         ? new TH2D("h2Truth",";Lead p_{T} (GeV);Subl p_{T} (GeV)",
0264                    nTrueLead,trueLeadPtBins.data(),nTrueSubl,trueSublPtBins.data()) : nullptr;
0265     TH2D* h2Ratio = isClosure
0266         ? new TH2D("h2Ratio",";Lead p_{T} (GeV);Subl p_{T} (GeV)",
0267                    nTrueLead,trueLeadPtBins.data(),nTrueSubl,trueSublPtBins.data()) : nullptr;
0268     h2Unf ->SetDirectory(0);
0269     h2Meas->SetDirectory(0);
0270     if (h2Truth) h2Truth->SetDirectory(0);
0271     if (h2Ratio) h2Ratio->SetDirectory(0);
0272 
0273     // Fill unfolded and truth (truth-index space)
0274     for (int f=0; f<nTrueFlat(); ++f) {
0275         int iL,iS; TrueFlatToIJ(f,iL,iS);
0276         double u=hJetPt_unfolded->GetBinContent(f+1), ue=hJetPt_unfolded->GetBinError(f+1);
0277         h2Unf->SetBinContent(iL+1,iS+1,u); h2Unf->SetBinError(iL+1,iS+1,ue);
0278         if (isClosure && hJetPt_truth) {
0279             double t=hJetPt_truth->GetBinContent(f+1), te=hJetPt_truth->GetBinError(f+1);
0280             h2Truth->SetBinContent(iL+1,iS+1,t); h2Truth->SetBinError(iL+1,iS+1,te);
0281             if (t>0) {
0282                 double r=u/t, re=r*std::sqrt((u>0?(ue/u)*(ue/u):0)+(te/t)*(te/t));
0283                 h2Ratio->SetBinContent(iL+1,iS+1,r); h2Ratio->SetBinError(iL+1,iS+1,re);
0284             }
0285         }
0286     }
0287     // Fill meas (reco-index space)
0288     for (int f=0; f<nRecoFlat(); ++f) {
0289         int iL,iS; RecoFlatToIJ(f,iL,iS);
0290         double m=hJetPt_meas->GetBinContent(f+1), me=hJetPt_meas->GetBinError(f+1);
0291         h2Meas->SetBinContent(iL+1,iS+1,m); h2Meas->SetBinError(iL+1,iS+1,me);
0292     }
0293 
0294     if (isClosure) {
0295         TCanvas* c1=new TCanvas("c1","Jet pT 2D",1600,550); c1->Divide(3,1);
0296         c1->cd(1); gPad->SetLogz(); h2Truth->SetTitle("Truth"); h2Truth->Draw("COLZ");
0297         c1->cd(2); gPad->SetLogz(); h2Unf->SetTitle("Unfolded"); h2Unf->Draw("COLZ");
0298         c1->cd(3); h2Ratio->SetTitle("Unfolded/Truth");
0299         h2Ratio->GetZaxis()->SetRangeUser(0.8,1.2); h2Ratio->Draw("COLZ TEXT");
0300         c1->SaveAs(std::format("{}/jetPT-{}.png",plotDir,label).c_str()); c1->Close(); delete c1;
0301 
0302         TH1D* hLT=(TH1D*)h2Truth->ProjectionX("hLT"), *hLU=(TH1D*)h2Unf->ProjectionX("hLU");
0303         TH1D* hST=(TH1D*)h2Truth->ProjectionY("hST"), *hSU=(TH1D*)h2Unf->ProjectionY("hSU");
0304         StyleLine(hLT,kColTruth,kMarTruth); StyleLine(hLU,kColUnfolded,kMarUnfolded);
0305         StyleLine(hST,kColTruth,kMarTruth); StyleLine(hSU,kColUnfolded,kMarUnfolded);
0306         TH1D* hLR=MakeRatioHist(hLU,hLT,"hLR"), *hSR=MakeRatioHist(hSU,hST,"hSR");
0307         TCanvas* c2=new TCanvas("c2","Jet pT Proj",1200,600); c2->Divide(2,1);
0308         for (int side=0; side<2; ++side) {
0309             auto [pTop,pBot]=MakeRatioPads(c2->cd(side+1));
0310             TH1D* hT=(side==0)?hLT:hST, *hU=(side==0)?hLU:hSU, *hR=(side==0)?hLR:hSR;
0311             const char* xL=(side==0)?"Lead p_{T} (GeV)":"Subl p_{T} (GeV)";
0312             pTop->cd(); pTop->SetLogy();
0313             hT->GetYaxis()->SetRangeUser(0.1,std::max(hT->GetMaximum(),hU->GetMaximum())*1.3);
0314             hT->GetXaxis()->SetLabelSize(0); hT->Draw("E"); hU->Draw("E SAME");
0315             TLegend* lg=new TLegend(0.55,0.72,0.88,0.88);
0316             lg->AddEntry(hT,"Truth","lp"); lg->AddEntry(hU,"Unfolded","lp"); lg->Draw();
0317             pBot->cd(); StyleRatio(hR,xL); hR->Draw("E"); DrawRefLine(hR);
0318         }
0319         c2->SaveAs(std::format("{}/jetPT_projections-{}.png",plotDir,label).c_str()); c2->Close(); delete c2;
0320     } else {
0321         TH1D* hLM=(TH1D*)h2Meas->ProjectionX("hLM"), *hLU=(TH1D*)h2Unf->ProjectionX("hLU");
0322         TH1D* hSM=(TH1D*)h2Meas->ProjectionY("hSM"), *hSU=(TH1D*)h2Unf->ProjectionY("hSU");
0323         StyleLine(hLM,kColMeas,kMarMeas); StyleLine(hLU,kColUnfolded,kMarUnfolded);
0324         StyleLine(hSM,kColMeas,kMarMeas); StyleLine(hSU,kColUnfolded,kMarUnfolded);
0325         TH1D* hLR=MakeRatioHist(hLU,hLM,"hLR"), *hSR=MakeRatioHist(hSU,hSM,"hSR");
0326         TCanvas* c2=new TCanvas("c2","Jet pT Reco vs Unf",1200,600); c2->Divide(2,1);
0327         for (int side=0; side<2; ++side) {
0328             auto [pTop,pBot]=MakeRatioPads(c2->cd(side+1));
0329             TH1D* hM=(side==0)?hLM:hSM, *hU=(side==0)?hLU:hSU, *hR=(side==0)?hLR:hSR;
0330             const char* xL=(side==0)?"Lead p_{T} (GeV)":"Subl p_{T} (GeV)";
0331             pTop->cd(); pTop->SetLogy();
0332             hM->GetYaxis()->SetRangeUser(0.1,std::max(hM->GetMaximum(),hU->GetMaximum())*1.3);
0333             hM->GetXaxis()->SetLabelSize(0); hM->Draw("E"); hU->Draw("E SAME");
0334             TLegend* lg=new TLegend(0.55,0.72,0.88,0.88);
0335             lg->AddEntry(hM,"Measured","lp"); lg->AddEntry(hU,"Unfolded","lp"); lg->Draw();
0336             pBot->cd(); StyleRatio(hR,xL,"Unf/Meas",0.5,1.5); hR->Draw("E"); DrawRefLine(hR);
0337         }
0338         c2->SaveAs(std::format("{}/jetPT_recoVsUnfolded-{}.png",plotDir,label).c_str());
0339         c2->Close(); delete c2;
0340     }
0341 
0342     // ════════════════════════════════════════════════════════════════════════
0343     //  wEEC plots — fWEEC is already open from the jet pT section above
0344     // ════════════════════════════════════════════════════════════════════════
0345 
0346     // ── Build per-dijet-pT-bin truth reference ────────────────────────────
0347     // For each ft, project hWEEC3D_truth_{k} for every k using
0348     // ProjectWEEC3DSliceForBin, assemble into a TH1D vs ΔΦ, then normalize.
0349     // For kData respFile points to the fullClosure response (full sim truth).
0350     std::vector<TH1D*> hTruthBin(nTrueFlat(), nullptr);
0351     if (fResp) {
0352         for (int ft=0; ft<nTrueFlat(); ++ft) {
0353             int iL,iS; TrueFlatToIJ(ft,iL,iS);
0354             if (!IsRecoAccessible(iL,iS)) { hTruthBin[ft] = nullptr; continue; }
0355             TH1D* hT = new TH1D(
0356                 std::format("hTruthBin_{}",ft).c_str(),
0357                 std::format("Truth iL{} iS{};#Delta#phi;wEEC",iL,iS).c_str(),
0358                 nDphi, dPhiBins.data());
0359             hT->Sumw2();
0360             hT->SetDirectory(0);
0361             for (int k=0; k<nDphi; ++k) {
0362                 TH1D* hT3D=(TH1D*)fResp->Get(truth3DName(k).c_str());
0363                 if (!hT3D) continue;
0364                 double val,err;
0365                 if (usePWMean && hPWNum[k] && hPWDen[k])
0366                     ProjectWEEC3DSliceForBin(hT3D, iL, iS, val, err, hPWNum[k], hPWDen[k]);
0367                 else
0368                     ProjectWEEC3DSliceForBin(hT3D, iL, iS, val, err);
0369                 hT->SetBinContent(k+1,val);
0370                 hT->SetBinError  (k+1,err);
0371             }
0372             NormalizeWEEC(hT);
0373             hTruthBin[ft] = hT;
0374         }
0375     }
0376 
0377     // ── Build inclusive truth reference ───────────────────────────────────
0378     TH1D* hInclTruth = nullptr;
0379     if (fResp) {
0380         hInclTruth = new TH1D("hInclTruth",";#Delta#phi;wEEC",nDphi,dPhiBins.data());
0381         hInclTruth->Sumw2();
0382         hInclTruth->SetDirectory(0);
0383         for (int k=0; k<nDphi; ++k) {
0384                 TH1D* hT3D=(TH1D*)fResp->Get(truth3DName(k).c_str());
0385                 if (!hT3D) continue;
0386                 double val,err;
0387                 // ProjectWEEC3DSlice already gates on IsRecoAccessible
0388                 if (usePWMean && hPWNum[k] && hPWDen[k])
0389                     ProjectWEEC3DSlice(hT3D, val, err, hPWNum[k], hPWDen[k]);
0390                 else
0391                     ProjectWEEC3DSlice(hT3D, val, err);
0392                 hInclTruth->SetBinContent(k+1,val);
0393                 hInclTruth->SetBinError  (k+1,err);
0394             }
0395         NormalizeWEEC(hInclTruth);
0396     }
0397 
0398     // ── Load measured per-bin histograms and normalize ────────────────────
0399     // For kFull/kHalf: hWEEC3D_meas_{k} from response file, project per bin
0400     // For kData: hWEEC3D_data_{k} from measFile (not available here, skip)
0401     std::vector<TH1D*> hMeasBin(nRecoFlat(), nullptr);
0402     {
0403         // Use fResp for kFull/kHalf (meas is in response file),
0404         // skip for kData (measured data file not passed to drawClosure)
0405         TFile* fMeasSrc = (!isData && fResp) ? fResp : nullptr;
0406         if (fMeasSrc) {
0407             for (int ft=0; ft<nRecoFlat(); ++ft) {
0408                 int iL,iS; RecoFlatToIJ(ft,iL,iS);
0409                 TH1D* hM = new TH1D(
0410                     std::format("hMeasBin_{}",ft).c_str(),
0411                     std::format("Meas iL{} iS{};#Delta#phi;wEEC",iL,iS).c_str(),
0412                     nDphi, dPhiBins.data());
0413                 hM->Sumw2();
0414                 hM->SetDirectory(0);
0415                 for (int k=0; k<nDphi; ++k) {
0416                     TH1D* hM3D=(TH1D*)fMeasSrc->Get(
0417                         std::format("hWEEC3D_meas_{}",k).c_str());
0418                     if (!hM3D) continue;
0419                     double val,err;
0420                     ProjectWEEC3DSliceForBin(hM3D,iL,iS,val,err);
0421                     hM->SetBinContent(k+1,val);
0422                     hM->SetBinError  (k+1,err);
0423                 }
0424                 NormalizeWEEC(hM);
0425                 hMeasBin[ft] = hM;
0426             }
0427         }
0428     }
0429 
0430     // ════════════════════════════════════════════════════════════════════════
0431     //  Inclusive wEEC
0432     // ════════════════════════════════════════════════════════════════════════
0433     if (hInclUnfolded) {
0434         TCanvas* cI=new TCanvas("cI","Inclusive wEEC",800,600);
0435         if (hInclTruth) {
0436             TPad* pTop=new TPad("pTop","",0,0.30,1,1.0);
0437             TPad* pBot=new TPad("pBot","",0,0.00,1,0.30);
0438             pTop->SetBottomMargin(0.03); pBot->SetTopMargin(0.03);
0439             pBot->SetBottomMargin(0.35); pTop->Draw(); pBot->Draw();
0440             pTop->cd(); pTop->SetLogy();
0441             StyleLine(hInclTruth,   kColTruth,   kMarTruth);
0442             StyleLine(hInclUnfolded,kColUnfolded,kMarUnfolded);
0443             double ymin=1e30,ymax=-1e30;
0444             for (TH1D* h:{hInclTruth,hInclUnfolded})
0445                 for (int i=1;i<=h->GetNbinsX();++i) {
0446                     double v=h->GetBinContent(i);
0447                     if (v>0) ymin=std::min(ymin,v);
0448                     ymax=std::max(ymax,v);
0449                 }
0450             hInclTruth->GetYaxis()->SetRangeUser(ymin*0.3,ymax*5.0);
0451             hInclTruth->GetYaxis()->SetTitle("1/N dN/d#Delta#phi");
0452             hInclTruth->GetXaxis()->SetLabelSize(0);
0453             hInclTruth->Draw("E"); hInclUnfolded->Draw("E SAME");
0454             TLegend* lg=new TLegend(0.55,0.72,0.90,0.88);
0455             lg->AddEntry(hInclTruth,"Truth (truth towers)","lp");
0456             lg->AddEntry(hInclUnfolded,"Unfolded","lp"); lg->Draw();
0457             pBot->cd();
0458             TH1D* hR=MakeRatioHist(hInclUnfolded,hInclTruth,"hInclRatio");
0459             StyleRatio(hR,"#Delta#phi"); hR->Draw("E"); DrawRefLine(hR);
0460         } else {
0461             cI->cd(); cI->SetLogy();
0462             StyleLine(hInclUnfolded,kColUnfolded,kMarUnfolded);
0463             hInclUnfolded->GetYaxis()->SetTitle("1/N dN/d#Delta#phi");
0464             hInclUnfolded->GetXaxis()->SetTitle("#Delta#phi");
0465             hInclUnfolded->Draw("E");
0466         }
0467         cI->SaveAs(std::format("{}/wEEC_inclusive-{}.png",plotDir,label).c_str());
0468         cI->Close(); delete cI;
0469     }
0470 
0471     // ════════════════════════════════════════════════════════════════════════
0472     //  Per-dijet-pT-bin grid helper
0473     // ════════════════════════════════════════════════════════════════════════
0474     const int padW=280, padH=220;
0475 
0476     // Compute shared y range across all populated bins (log scale)
0477     auto yRange1 = [](const std::vector<TH1D*>& hA,
0478                       const std::vector<TH1D*>& hB,
0479                       const std::vector<TH1D*>& hC)
0480         -> std::pair<double,double>
0481     {
0482         double yLo=std::numeric_limits<double>::max();
0483         double yHi=-std::numeric_limits<double>::max();
0484         for (int ft=0; ft<nTrueFlat(); ++ft)
0485             for (TH1D* h : {hA[ft], hB[ft], hC[ft]}) {
0486                 if (!h) continue;
0487                 for (int i=1;i<=h->GetNbinsX();++i) {
0488                     double v=h->GetBinContent(i);
0489                     if (v>0) yLo=std::min(yLo,v);
0490                     yHi=std::max(yHi,v);
0491                 }
0492             }
0493         if (yLo>=yHi) return {1e-6,1.0};
0494         return {yLo*0.3, yHi*5.0};
0495     };
0496 
0497     // Grid drawing — use a plain nested scope instead of a lambda to avoid
0498     // Cling interpreter crashes on lambda destruction at function exit.
0499     // Called three times below via a simple struct-less inline pattern.
0500     auto DrawGrid = [=](const char* canName, const char* outFile,
0501                         const std::vector<TH1D*>& hA,
0502                         const std::vector<TH1D*>& hB,
0503                         const std::vector<TH1D*>& hC,
0504                         bool isRatio,
0505                         double yLo, double yHi,
0506                         const char* legA, const char* legB, const char* legC,
0507                         const std::vector<std::pair<int,int>>& binMask = {})
0508     {
0509         TCanvas* cG=new TCanvas(canName,canName,nTrueLead*padW,nTrueSubl*padH);
0510         cG->Divide(nTrueLead,nTrueSubl,0,0);
0511         for (int iS=0; iS<nTrueSubl; ++iS)
0512         for (int iL=0; iL<nTrueLead; ++iL)
0513         {
0514             cG->cd(iS*nTrueLead+iL+1);
0515             gPad->SetLeftMargin(0.18); gPad->SetBottomMargin(0.18);
0516             gPad->SetRightMargin(0.04); gPad->SetTopMargin(0.10);
0517             TLatex tex; tex.SetNDC(); tex.SetTextSize(0.09); tex.SetTextAlign(13);
0518             std::string bl=std::format("{:.0f}-{:.0f}/{:.0f}-{:.0f} GeV",
0519                 trueLeadPtBins[iL],trueLeadPtBins[iL+1],trueSublPtBins[iS],trueSublPtBins[iS+1]);
0520 
0521             int ft = TrueFlatIndex(iL,iS);
0522 
0523             // Gray out kinematically forbidden bins
0524             if (ft < 0 || !hA[ft]) {
0525                 gPad->SetFillColor(kGray);
0526                 tex.DrawLatex(0.22,0.88,bl.c_str()); continue;
0527             }
0528 
0529             // Gray out analysis-quality masked bins (darker than forbidden)
0530             bool isMasked = std::any_of(binMask.begin(), binMask.end(),
0531                 [&](const std::pair<int,int>& p){ return p.first==iL && p.second==iS; });
0532             if (isMasked) {
0533                 gPad->SetFillColor(kGray+2);
0534                 tex.DrawLatex(0.22,0.88,bl.c_str()); continue;
0535             }
0536 
0537             if (isRatio) {
0538                 TH1D* hH=hA[ft];
0539                 hH->SetLineColor(kBlack); hH->SetMarkerColor(kBlack);
0540                 hH->SetMarkerStyle(20);
0541                 hH->GetYaxis()->SetRangeUser(yLo,yHi);
0542                 hH->GetYaxis()->SetTitle("Unf/Truth");
0543                 hH->GetYaxis()->SetTitleSize(0.09); hH->GetYaxis()->SetLabelSize(0.08);
0544                 hH->GetYaxis()->SetNdivisions(504);
0545                 hH->GetXaxis()->SetTitle("#Delta#phi");
0546                 hH->GetXaxis()->SetTitleSize(0.09); hH->GetXaxis()->SetLabelSize(0.08);
0547                 hH->SetTitle("");
0548                 hH->Draw("E"); DrawRefLine(hH);
0549             } else {
0550                 gPad->SetLogy();
0551                 TH1D* hBase = hB[ft] ? hB[ft] : (hC[ft] ? hC[ft] : hA[ft]);
0552                 hBase->GetYaxis()->SetRangeUser(yLo,yHi);
0553                 hBase->GetYaxis()->SetTitle("dwEEC/d#Delta#phi");
0554                 hBase->GetYaxis()->SetTitleSize(0.09); hBase->GetYaxis()->SetTitleOffset(0.9);
0555                 hBase->GetYaxis()->SetLabelSize(0.08);
0556                 hBase->GetXaxis()->SetTitle("#Delta#phi");
0557                 hBase->GetXaxis()->SetTitleSize(0.09); hBase->GetXaxis()->SetLabelSize(0.08);
0558                 hBase->SetTitle("");
0559                 StyleLine(hA[ft],kColUnfolded,kMarUnfolded);
0560                 if (hB[ft]) StyleLine(hB[ft],kColTruth,kMarTruth);
0561                 if (hC[ft]) StyleLine(hC[ft],kColMeas,  kMarMeas);
0562                 hBase->Draw("E");
0563                 hA[ft]->Draw("E SAME");
0564                 if (hB[ft] && hB[ft]!=hBase) hB[ft]->Draw("E SAME");
0565                 if (hC[ft] && hC[ft]!=hBase) hC[ft]->Draw("E SAME");
0566                 TLegend* lg=new TLegend(0.40,0.72,0.96,0.88);
0567                 lg->SetTextSize(0.07);
0568                 lg->AddEntry(hA[ft],legA,"lp");
0569                 if (hB[ft]) lg->AddEntry(hB[ft],legB,"lp");
0570                 if (hC[ft]) lg->AddEntry(hC[ft],legC,"lp");
0571                 lg->Draw();
0572             }
0573             tex.DrawLatex(0.22,0.88,bl.c_str());
0574         }
0575         cG->cd();
0576         TLatex ct; ct.SetNDC(); ct.SetTextSize(0.022); ct.SetTextAlign(22);
0577         ct.DrawLatex(0.5,0.003,"Lead jet p_{T} bin #rightarrow");
0578         ct.SetTextAngle(90); ct.DrawLatex(0.006,0.5,"Subl jet p_{T} bin #rightarrow");
0579         cG->SaveAs(outFile); cG->Close(); delete cG;
0580     };
0581 
0582     // ════════════════════════════════════════════════════════════════════════
0583     //  Plot 1: unfolded vs truth  (wEEC_perBin)
0584     // ════════════════════════════════════════════════════════════════════════
0585     {
0586         std::vector<TH1D*> hNull(nTrueFlat(), nullptr);
0587         auto [yLo,yHi] = yRange1(hUnfBin, hTruthBin, hNull);
0588         DrawGrid("cPerBin",
0589                  std::format("{}/wEEC_perBin-{}.png",plotDir,label).c_str(),
0590                  hUnfBin, hTruthBin, hNull,
0591                  false, yLo, yHi,
0592                  "Unfolded", "Truth (truth towers)", "");
0593         DrawGrid("cPerBinMasked",
0594                  std::format("{}/wEEC_perBin_masked-{}.png",plotDir,label).c_str(),
0595                  hUnfBin, hTruthBin, hNull,
0596                  false, yLo, yHi,
0597                  "Unfolded", "Truth (truth towers)", "", maskedBins);
0598 
0599         for(int ft=0; ft<nTrueFlat(); ++ft) {
0600             if(!hUnfBin[ft] || !hTruthBin[ft]) continue;
0601 
0602             TCanvas *cBinByBin = new TCanvas("cBinByBin","",padW,padH);
0603             cBinByBin->SetLogy();
0604 
0605             int iL, iS;
0606             TrueFlatToIJ(ft, iL, iS);
0607             gPad->SetLeftMargin(0.18); gPad->SetBottomMargin(0.18);
0608             gPad->SetRightMargin(0.04); gPad->SetTopMargin(0.10);
0609             TLatex tex; tex.SetNDC(); tex.SetTextSize(0.09); tex.SetTextAlign(13);
0610             std::string bl=std::format("{:.0f}-{:.0f}/{:.0f}-{:.0f} GeV",
0611                 trueLeadPtBins[iL],trueLeadPtBins[iL+1],trueSublPtBins[iS],trueSublPtBins[iS+1]);
0612 
0613             hUnfBin[ft]->GetYaxis()->SetTitle("dwEEC/d#Delta#phi");
0614             hUnfBin[ft]->GetYaxis()->SetTitleSize(0.09); hUnfBin[ft]->GetYaxis()->SetTitleOffset(0.9);
0615             hUnfBin[ft]->GetYaxis()->SetLabelSize(0.08);
0616             hUnfBin[ft]->GetXaxis()->SetTitle("#Delta#phi");
0617             hUnfBin[ft]->GetXaxis()->SetTitleSize(0.09); hUnfBin[ft]->GetXaxis()->SetLabelSize(0.08);
0618             hUnfBin[ft]->SetTitle("");
0619             
0620             StyleLine(hUnfBin[ft],kColUnfolded,kMarUnfolded);
0621             StyleLine(hTruthBin[ft],kColTruth,kMarTruth);
0622 
0623             hUnfBin[ft]->Draw("E");
0624             hTruthBin[ft]->Draw("E SAME");
0625 
0626             TLegend* lg=new TLegend(0.40,0.72,0.96,0.88);
0627             lg->SetTextSize(0.07);
0628             lg->AddEntry(hUnfBin[ft],"Unfolded","lp");
0629             lg->AddEntry(hTruthBin[ft],"Truth (truth towers)","lp");
0630             lg->Draw();
0631 
0632             tex.DrawLatex(0.22,0.88,bl.c_str());
0633 
0634             cBinByBin->SaveAs(std::format("{}/wEEC_singleBinBin-{}-{}-{}.png",plotDir,iL,iS,label).c_str());
0635 
0636             delete cBinByBin;
0637         }
0638     }
0639 
0640     // ════════════════════════════════════════════════════════════════════════
0641     //  Plot 2: ratio unfolded/truth  (wEEC_ratio)  — only if truth available
0642     // ════════════════════════════════════════════════════════════════════════
0643     {
0644         std::vector<TH1D*> hRatioBin(nTrueFlat(), nullptr);
0645         for (int ft=0; ft<nTrueFlat(); ++ft) {
0646             if (hUnfBin[ft] && hTruthBin[ft])
0647                 hRatioBin[ft] = MakeRatioHist(hUnfBin[ft], hTruthBin[ft],
0648                                     std::format("hRatioBin_{}",ft).c_str());
0649         }
0650         std::vector<TH1D*> hNull(nTrueFlat(), nullptr);
0651         DrawGrid("cRatio",
0652                  std::format("{}/wEEC_ratio-{}.png",plotDir,label).c_str(),
0653                  hRatioBin, hNull, hNull,
0654                  true, 0.5, 1.5,
0655                  "Unf/Truth", "", "");
0656         DrawGrid("cRatioMasked",
0657                  std::format("{}/wEEC_ratio_masked-{}.png",plotDir,label).c_str(),
0658                  hRatioBin, hNull, hNull,
0659                  true, 0.5, 1.5,
0660                  "Unf/Truth", "", "", maskedBins);
0661 
0662         for(int ft=0; ft<nTrueFlat(); ++ft) {
0663             if(!hUnfBin[ft] || !hTruthBin[ft]) continue;
0664 
0665             TCanvas *cRBinByBin = new TCanvas("cRBinByBin","",padW,padH);
0666             
0667             int iL, iS;
0668             TrueFlatToIJ(ft, iL, iS);
0669             gPad->SetLeftMargin(0.18); gPad->SetBottomMargin(0.18);
0670             gPad->SetRightMargin(0.04); gPad->SetTopMargin(0.10);
0671             TLatex tex; tex.SetNDC(); tex.SetTextSize(0.09); tex.SetTextAlign(13);
0672             std::string bl=std::format("{:.0f}-{:.0f}/{:.0f}-{:.0f} GeV",
0673                 trueLeadPtBins[iL],trueLeadPtBins[iL+1],trueSublPtBins[iS],trueSublPtBins[iS+1]);
0674 
0675             hRatioBin[ft]->GetYaxis()->SetRangeUser(0.5,1.5);
0676             hRatioBin[ft]->GetYaxis()->SetTitle("Unf/Truth");
0677             hRatioBin[ft]->GetYaxis()->SetTitleSize(0.09); hRatioBin[ft]->GetYaxis()->SetTitleOffset(0.9);
0678             hRatioBin[ft]->GetYaxis()->SetLabelSize(0.08);
0679             hRatioBin[ft]->GetXaxis()->SetTitle("#Delta#phi");
0680             hRatioBin[ft]->GetXaxis()->SetTitleSize(0.09); hRatioBin[ft]->GetXaxis()->SetLabelSize(0.08);
0681             hRatioBin[ft]->SetTitle("");
0682             
0683             StyleLine(hRatioBin[ft],kColUnfolded,kMarUnfolded);
0684 
0685             hRatioBin[ft]->Draw("E");
0686             DrawRefLine(hRatioBin[ft]);
0687 
0688             tex.DrawLatex(0.22,0.88,bl.c_str());
0689 
0690             cRBinByBin->SaveAs(std::format("{}/wEEC_ratio_singleBinBin-{}-{}-{}.png",plotDir,iL,iS,label).c_str());
0691 
0692             delete cRBinByBin;
0693         }
0694         for (int ft=0; ft<nTrueFlat(); ++ft) delete hRatioBin[ft];
0695     }
0696 
0697     // ════════════════════════════════════════════════════════════════════════
0698     //  Plot 3: measured vs unfolded vs truth  (wEEC_measVsUnf)
0699     // ════════════════════════════════════════════════════════════════════════
0700     {
0701         // Expand reco-indexed meas into truth-indexed space for DrawGrid
0702         std::vector<TH1D*> hMeasBinTrue(nTrueFlat(), nullptr);
0703         for (int ftR=0; ftR<nRecoFlat(); ++ftR) {
0704             if (!hMeasBin[ftR]) continue;
0705             int iLr,iSr; RecoFlatToIJ(ftR,iLr,iSr);
0706             int iLt=FindBin(recoLeadPtBins[iLr],trueLeadPtBins);
0707             int iSt=FindBin(recoSublPtBins[iSr], trueSublPtBins);
0708             if (iLt>=0 && iSt>=0) hMeasBinTrue[TrueFlatIndex(iLt,iSt)]=hMeasBin[ftR];
0709         }
0710         std::vector<TH1D*> hNull(nTrueFlat(), nullptr);
0711         auto [yLo,yHi] = yRange1(hUnfBin, hTruthBin, !isData ? hMeasBinTrue : hNull);
0712         DrawGrid("cMeasVsUnf",
0713                  std::format("{}/wEEC_measVsUnf-{}.png",plotDir,label).c_str(),
0714                  hUnfBin,
0715                  hTruthBin,
0716                  !isData ? hMeasBinTrue : hNull,
0717                  false, yLo, yHi,
0718                  "Unfolded", "Truth (truth towers)",
0719                  !isData ? "Measured" : "");
0720         DrawGrid("cMeasVsUnfMasked",
0721                  std::format("{}/wEEC_measVsUnf_masked-{}.png",plotDir,label).c_str(),
0722                  hUnfBin,
0723                  hTruthBin,
0724                  !isData ? hMeasBinTrue : hNull,
0725                  false, yLo, yHi,
0726                  "Unfolded", "Truth (truth towers)",
0727                  !isData ? "Measured" : "", maskedBins);
0728     }
0729 
0730     // ════════════════════════════════════════════════════════════════════════
0731     //  1D closure projections: lead pT, subl pT, pair weight
0732     //  wEEC_projections-{label}.png
0733     // ════════════════════════════════════════════════════════════════════════
0734     //  wEEC projections: one image per ΔΦ bin
0735     //  For each ΔΦ bin k, project hWEEC3D_unfolded_{k} and hWEEC3D_truth_{k}
0736     //  onto lead pT, subl pT, and pair weight by summing the flat3D bins.
0737     //  All (iL,iS) truth bins are included — no IsRecoAccessible gate — so
0738     //  the full closure of each independently unfolded ΔΦ bin is visible.
0739     //  Saved to plotDir/wEEC_projections_{k}-{label}.png
0740     // ════════════════════════════════════════════════════════════════════════
0741     if (isClosure && fResp) {
0742         for (int k=0; k<nDphi; ++k) {
0743             TH1D* hU3D = (TH1D*)fWEEC->Get(std::format("hWEEC3D_unfolded_{}",k).c_str());
0744             TH1D* hT3D = (TH1D*)fResp->Get(truth3DName(k).c_str());
0745             if (!hU3D || !hT3D) continue;
0746 
0747             TH1D* hUnfLead = new TH1D(std::format("hUnfLead_{}",k).c_str(),"",
0748                 nTrueLead,trueLeadPtBins.data());
0749             TH1D* hUnfSubl = new TH1D(std::format("hUnfSubl_{}",k).c_str(),"",
0750                 nTrueSubl,trueSublPtBins.data());
0751             TH1D* hUnfPW   = new TH1D(std::format("hUnfPW_{}",k).c_str(),"",
0752                 nPairWeight,pairWeightBins.data());
0753             TH1D* hTruLead = new TH1D(std::format("hTruLead_{}",k).c_str(),"",
0754                 nTrueLead,trueLeadPtBins.data());
0755             TH1D* hTruSubl = new TH1D(std::format("hTruSubl_{}",k).c_str(),"",
0756                 nTrueSubl,trueSublPtBins.data());
0757             TH1D* hTruPW   = new TH1D(std::format("hTruPW_{}",k).c_str(),"",
0758                 nPairWeight,pairWeightBins.data());
0759             for (TH1D* h:{hUnfLead,hUnfSubl,hUnfPW,hTruLead,hTruSubl,hTruPW})
0760                 { h->SetDirectory(0); h->Sumw2(); }
0761 
0762             // Project flat3D bins onto lead pT, subl pT, and pair weight axes.
0763             // For the unfolded histograms, use the full covariance matrix so
0764             // off-diagonal correlations are included when bins are summed.
0765             // Truth bins use independent Sumw2 errors — quadrature is correct there.
0766 
0767             // First pass: accumulate values and covariance sums for unfolded.
0768             // For each projected bin B, Var(sum_B x_i) = sum_{i,j in B} C_{ij}.
0769             // We accumulate the variance incrementally as we visit each (i,j) pair.
0770             // Map from projected histogram bin index → variance accumulator.
0771             std::map<int,double> varL, varS, varP;
0772 
0773             for (int ft=0; ft<nTrueFlat(); ++ft)
0774             for (int iPw=0; iPw<nPairWeight; ++iPw) {
0775                 int iL, iS; TrueFlatToIJ(ft, iL, iS);
0776                 if (projAccessibleOnly && !IsRecoAccessible(iL, iS)) continue;
0777                 int iF = ft * nPairWeight + iPw;
0778                 double lCen = 0.5*(trueLeadPtBins[iL]+trueLeadPtBins[iL+1]);
0779                 double sCen = 0.5*(trueSublPtBins[iS]+trueSublPtBins[iS+1]);
0780                 double den   = (usePWMean && hPWDen[k]) ? hPWDen[k]->GetBinContent(iF+1) : 0.0;
0781                 double wMean = (den > 0)
0782                     ? hPWNum[k]->GetBinContent(iF+1) / den
0783                     : 0.5*(pairWeightBins[iPw]+pairWeightBins[iPw+1]);
0784 
0785                 int bL = hUnfLead->FindBin(lCen);
0786                 int bS = hUnfSubl->FindBin(sCen);
0787                 int bP = hUnfPW  ->FindBin(wMean);
0788 
0789                 // Accumulate values (weights=1 for lead/subl, wMean for pair weight)
0790                 double uV = hU3D->GetBinContent(iF+1);
0791                 hUnfLead->AddBinContent(bL, uV);
0792                 hUnfSubl->AddBinContent(bS, uV);
0793                 hUnfPW  ->AddBinContent(bP, uV);
0794 
0795                 // Truth: independent bins, quadrature sum is correct
0796                 double tV = hT3D->GetBinContent(iF+1), tE = hT3D->GetBinError(iF+1);
0797                 double tBLE = hTruLead->GetBinError(bL);
0798                 double tBSE = hTruSubl->GetBinError(bS);
0799                 double tBPE = hTruPW  ->GetBinError(bP);
0800                 hTruLead->AddBinContent(bL, tV);
0801                 hTruSubl->AddBinContent(bS, tV);
0802                 hTruPW  ->AddBinContent(bP, tV);
0803                 hTruLead->SetBinError(bL, std::hypot(tBLE, tE));
0804                 hTruSubl->SetBinError(bS, std::hypot(tBSE, tE));
0805                 hTruPW  ->SetBinError(bP, std::hypot(tBPE, tE));
0806 
0807                 // Unfolded variance: sum C_{ij} over all j in the same projected bin.
0808                 // We do this as a second inner loop over j.
0809                 for (int ft2=0; ft2<nTrueFlat(); ++ft2)
0810                 for (int jPw=0; jPw<nPairWeight; ++jPw) {
0811                     int jL, jS; TrueFlatToIJ(ft2, jL, jS);
0812                     if (projAccessibleOnly && !IsRecoAccessible(jL, jS)) continue;
0813                     int jF = ft2 * nPairWeight + jPw;
0814                     double jlCen = 0.5*(trueLeadPtBins[jL]+trueLeadPtBins[jL+1]);
0815                     double jsCen = 0.5*(trueSublPtBins[jS]+trueSublPtBins[jS+1]);
0816                     double jden   = (usePWMean && hPWDen[k]) ? hPWDen[k]->GetBinContent(jF+1) : 0.0;
0817                     double jwMean = (jden > 0)
0818                         ? hPWNum[k]->GetBinContent(jF+1) / jden
0819                         : 0.5*(pairWeightBins[jPw]+pairWeightBins[jPw+1]);
0820                     int jbL = hUnfLead->FindBin(jlCen);
0821                     int jbS = hUnfSubl->FindBin(jsCen);
0822                     int jbP = hUnfPW  ->FindBin(jwMean);
0823 
0824                     double cij = hCov[k] ? (*hCov[k])(iF, jF)
0825                                          : (iF==jF ? std::pow(hU3D->GetBinError(iF+1),2) : 0.0);
0826 
0827                     // Only accumulate if both bins project to the same output bin
0828                     if (bL == jbL) varL[bL] += cij;
0829                     if (bS == jbS) varS[bS] += cij;
0830                     if (bP == jbP) varP[bP] += cij;
0831                 }
0832             }
0833             // Set unfolded errors from accumulated variances
0834             for (auto& [b, v] : varL) hUnfLead->SetBinError(b, std::sqrt(std::max(v,0.0)));
0835             for (auto& [b, v] : varS) hUnfSubl->SetBinError(b, std::sqrt(std::max(v,0.0)));
0836             for (auto& [b, v] : varP) hUnfPW  ->SetBinError(b, std::sqrt(std::max(v,0.0)));
0837 
0838             struct ProjSpec { TH1D* hU; TH1D* hT; const char* xL; bool logx; };
0839             std::vector<ProjSpec> projs = {
0840                 {hUnfLead,hTruLead,"Lead p_{T} (GeV)",false},
0841                 {hUnfSubl,hTruSubl,"Subl p_{T} (GeV)",false},
0842                 {hUnfPW,  hTruPW,  "Pair weight",      true }
0843             };
0844             std::string dphiStr = std::format("#Delta#phi bin {}: [{:.2f},{:.2f})",
0845                 k, dPhiBins[k], dPhiBins[k+1]);
0846             TCanvas* cProj = new TCanvas(
0847                 std::format("cProj_{}",k).c_str(), "", 1800, 600);
0848             cProj->Divide(3,1);
0849             for (int p=0; p<3; ++p) {
0850                 auto [hU,hT,xL,logx] = projs[p];
0851                 StyleLine(hT,kColTruth,kMarTruth);
0852                 StyleLine(hU,kColUnfolded,kMarUnfolded);
0853                 auto [pTop,pBot] = MakeRatioPads(cProj->cd(p+1));
0854                 pTop->cd(); pTop->SetLogy(); if (logx) pTop->SetLogx();
0855                 double ymax = std::max(hT->GetMaximum(), hU->GetMaximum());
0856                 double ymin = 1e30;
0857                 for (int b=1; b<=hT->GetNbinsX(); ++b)
0858                     if (hT->GetBinContent(b)>0) ymin=std::min(ymin,hT->GetBinContent(b));
0859                 if (ymin>ymax) ymin=ymax*1e-4;
0860                 hT->GetYaxis()->SetRangeUser(ymin*0.3, ymax*3.0);
0861                 hT->GetXaxis()->SetLabelSize(0);
0862                 hT->GetYaxis()->SetTitle("Pairs (arb.)");
0863                 hT->SetTitle(dphiStr.c_str());
0864                 hT->Draw("E"); hU->Draw("E SAME");
0865                 TLegend* lg = new TLegend(0.55,0.72,0.88,0.88);
0866                 lg->AddEntry(hT,"Truth","lp"); lg->AddEntry(hU,"Unfolded","lp");
0867                 lg->Draw();
0868                 pBot->cd(); if (logx) pBot->SetLogx();
0869                 TH1D* hR = MakeRatioHist(hU,hT,std::format("hProjR_{}_{}",k,p).c_str());
0870                 StyleRatio(hR,xL); hR->Draw("E"); DrawRefLine(hR);
0871             }
0872             cProj->SaveAs(std::format("{}/wEEC_projections_{}-{}.png",
0873                 plotDir, k, label).c_str());
0874             cProj->Close(); delete cProj;
0875             for (TH1D* h:{hUnfLead,hUnfSubl,hUnfPW,hTruLead,hTruSubl,hTruPW})
0876                 delete h;
0877         }
0878     }
0879 
0880     // ════════════════════════════════════════════════════════════════════════
0881     //  Pair weight closure: per ΔΦ bin grid (pairWeight_closure_dPhi)
0882     //  and per (ΔΦ, dijet pT) bin PDF (pairWeight_closure_perBin)
0883     // ════════════════════════════════════════════════════════════════════════
0884     if (isClosure && fResp) {
0885         // Build pair weight projections: hPWUnf[k] and hPWTru[k] (summed over dijet pT)
0886         // and hPWUnfBin[k*nTrueFlat()+ft] / hPWTruBin[k*nTrueFlat()+ft] (per dijet pT bin)
0887         std::vector<TH1D*> hPWUnf(nDphi,nullptr), hPWTru(nDphi,nullptr);
0888         std::vector<TH1D*> hPWUnfBin(nDphi*nTrueFlat(),nullptr);
0889         std::vector<TH1D*> hPWTruBin(nDphi*nTrueFlat(),nullptr);
0890         std::cout << "nPairWeight: " << nPairWeight << std::endl;
0891         for(int k=0; k<=nPairWeight; k++)
0892         {
0893             std::cout << "   pw " << k << ": " << pairWeightBins[k] << std::endl;
0894         }
0895         for (int k=0; k<nDphi; ++k) {
0896             hPWUnf[k]=new TH1D(std::format("hPWUnf_{}",k).c_str(),"",nPairWeight,pairWeightBins.data());
0897             hPWTru[k]=new TH1D(std::format("hPWTru_{}",k).c_str(),"",nPairWeight,pairWeightBins.data());
0898             hPWUnf[k]->SetDirectory(0); hPWTru[k]->SetDirectory(0);
0899             hPWUnf[k]->Sumw2();        hPWTru[k]->Sumw2();
0900             for (int ft=0; ft<nTrueFlat(); ++ft) {
0901                 hPWUnfBin[k*nTrueFlat()+ft]=new TH1D(
0902                     std::format("hPWUnfBin_{}_{}",k,ft).c_str(),"",nPairWeight,pairWeightBins.data());
0903                 hPWTruBin[k*nTrueFlat()+ft]=new TH1D(
0904                     std::format("hPWTruBin_{}_{}",k,ft).c_str(),"",nPairWeight,pairWeightBins.data());
0905                 hPWUnfBin[k*nTrueFlat()+ft]->SetDirectory(0);
0906                 hPWTruBin[k*nTrueFlat()+ft]->SetDirectory(0);
0907                 hPWUnfBin[k*nTrueFlat()+ft]->Sumw2();
0908                 hPWTruBin[k*nTrueFlat()+ft]->Sumw2();
0909             }
0910         }
0911 
0912         // Read unfolded 3D histograms directly from already-open fWEEC
0913         for (int k=0; k<nDphi; ++k) {
0914                 TH1D* hU=(TH1D*)fWEEC->Get(std::format("hWEEC3D_unfolded_{}",k).c_str());
0915                 TH1D* hT=(TH1D*)fResp->Get(std::format("hWEEC3D_truth_{}",k).c_str());
0916                 for (int ft=0; ft<nTrueFlat(); ++ft)
0917                 for (int iPw=0; iPw<nPairWeight; ++iPw) {
0918                     int iL, iS; TrueFlatToIJ(ft, iL, iS);
0919                     if (!IsRecoAccessible(iL, iS)) continue;
0920                     int iF = ft * nPairWeight + iPw;
0921                     double den   = (usePWMean && hPWDen[k]) ? hPWDen[k]->GetBinContent(iF+1) : 0.0;
0922                     double wMean = (den > 0)
0923                         ? hPWNum[k]->GetBinContent(iF+1) / den
0924                         : 0.5*(pairWeightBins[iPw]+pairWeightBins[iPw+1]);
0925                     if (hU) {
0926                             double v=hU->GetBinContent(iF+1), e=hU->GetBinError(iF+1);
0927                             int bP=hPWUnf[k]->FindBin(wMean);
0928                             int bPb=hPWUnfBin[k*nTrueFlat()+ft]->FindBin(wMean);
0929                             double bPe = hPWUnf[k]->GetBinError(bP);
0930                             double bPbe = hPWUnfBin[k*nTrueFlat()+ft]->GetBinError(bPbe);
0931                             hPWUnf[k]->Fill(wMean,v);
0932                             hPWUnfBin[k*nTrueFlat()+ft]->Fill(wMean,v);
0933                             hPWUnf[k]->SetBinError(bP,std::hypot(bPe,e));
0934                             hPWUnfBin[k*nTrueFlat()+ft]->SetBinError(bPb,std::hypot(bPbe,e));
0935                     }
0936                     if (hT) {
0937                             double v=hT->GetBinContent(iF+1);
0938                             double e=hT->GetBinError(iF+1);
0939                             int bP = hPWTru[k]->FindBin(wMean);
0940                             int bPb = hPWTruBin[k*nTrueFlat()+ft]->FindBin(wMean);
0941                             double bPe = hPWTru[k]->GetBinError(bP);
0942                             double bPbe = hPWTruBin[k*nTrueFlat()+ft]->GetBinError(bPb);
0943                             hPWTru[k]->Fill(wMean,v);
0944                             hPWTruBin[k*nTrueFlat()+ft]->Fill(wMean,v);
0945                             hPWTru[k]->SetBinError(bP,std::hypot(bPe,e));
0946                             hPWTruBin[k*nTrueFlat()+ft]->SetBinError(bPb,std::hypot(bPbe,e));
0947                     }
0948                 }
0949         }
0950 
0951         // ── Plot 1: 4×8 ΔΦ grid ──────────────────────────────────────────
0952         const int nColsPW=4, nRowsPW=(nDphi+nColsPW-1)/nColsPW;
0953         TCanvas* cPWGrid=new TCanvas("cPWGrid","PW closure dPhi",nColsPW*250,nRowsPW*500);
0954         cPWGrid->Divide(nColsPW,nRowsPW,0,0);
0955         for (int k=0; k<nDphi; ++k) {
0956             StyleLine(hPWTru[k],kColTruth,kMarTruth);
0957             StyleLine(hPWUnf[k],kColUnfolded,kMarUnfolded);
0958             auto [pTop,pBot]=MakeRatioPads(cPWGrid->cd(k+1));
0959             pTop->cd(); pTop->SetLogy(); pTop->SetLogx();
0960             double ymax=std::max(hPWTru[k]->GetMaximum(),hPWUnf[k]->GetMaximum());
0961             double ymin=1e30;
0962             for (int b=1;b<=hPWTru[k]->GetNbinsX();++b)
0963                 if (hPWTru[k]->GetBinContent(b)>0) ymin=std::min(ymin,hPWTru[k]->GetBinContent(b));
0964             if (ymin>ymax) ymin=ymax*1e-4;
0965             hPWTru[k]->GetYaxis()->SetRangeUser(ymin*0.3,ymax*3.0);
0966             hPWTru[k]->GetXaxis()->SetLabelSize(0);
0967             hPWTru[k]->GetYaxis()->SetTitle("Pairs");
0968             hPWTru[k]->GetYaxis()->SetTitleSize(0.07);
0969             hPWTru[k]->SetTitle(std::format("{:.2f}-{:.2f}",dPhiBins[k],dPhiBins[k+1]).c_str());
0970             hPWTru[k]->Draw("E"); hPWUnf[k]->Draw("E SAME");
0971             TLegend* lg=new TLegend(0.55,0.72,0.92,0.88);
0972             lg->SetTextSize(0.07);
0973             lg->AddEntry(hPWTru[k],"Truth","lp"); lg->AddEntry(hPWUnf[k],"Unfolded","lp"); lg->Draw();
0974             pBot->cd(); pBot->SetLogx();
0975             TH1D* hR=MakeRatioHist(hPWUnf[k],hPWTru[k],std::format("hPWRatio_{}",k).c_str());
0976             StyleRatio(hR,"Pair weight"); hR->Draw("E"); DrawRefLine(hR);
0977         }
0978         cPWGrid->SaveAs(std::format("{}/pairWeight_closure_dPhi-{}.png",plotDir,label).c_str());
0979         cPWGrid->Close(); delete cPWGrid;
0980 
0981         // ── Plot 2: PDF — one page per ΔΦ bin, nTrueLead×nTrueSubl grid ──
0982             
0983         TH2D *hFrameA = new TH2D("hFrameA","",1,4e-8,2.0,1,1e-3,1e12);
0984         TH2D *hFrameB = new TH2D("hFrameB","",1,4e-8,2.0,1,0.5,1.5);        
0985         std::string pdfName=std::format("{}/pairWeight_closure_perBin-{}.pdf",plotDir,label);
0986         for (int k=0; k<nDphi; ++k) {
0987             TCanvas* cPage=new TCanvas(
0988                 std::format("cPWPage_{}",k).c_str(),"",nTrueLead*220,nTrueSubl*220);
0989 
0990             cPage->Divide(nTrueLead,nTrueSubl,0,0);
0991             for (int iS=0; iS<nTrueSubl; ++iS)
0992             for (int iL=0; iL<nTrueLead; ++iL) {
0993                 int ft=TrueFlatIndex(iL,iS);
0994                 int cell=iS*nTrueLead+iL;
0995                 TVirtualPad* parent=cPage->cd(cell+1);
0996                 std::string bl=std::format("{:.0f}-{:.0f}/{:.0f}-{:.0f} GeV",
0997                     trueLeadPtBins[iL],trueLeadPtBins[iL+1],
0998                     trueSublPtBins[iS],trueSublPtBins[iS+1]);
0999 
1000                 // Gray out kinematically forbidden bins — TRatioPlot needs
1001                 // valid histograms so handle these before attempting to build it.
1002                 if (ft < 0) {
1003                     parent->SetFillColor(kGray);
1004                     parent->cd();
1005                     TLatex tx; tx.SetNDC(); tx.SetTextSize(0.08);
1006                     tx.DrawLatex(0.15,0.85,bl.c_str()); continue;
1007                 }
1008 
1009                 TH1D* hU=hPWUnfBin[k*nTrueFlat()+ft];
1010                 TH1D* hT=hPWTruBin[k*nTrueFlat()+ft];
1011 
1012                 if (hT->Integral()==0) {
1013                     parent->SetFillColor(kGray+1);
1014                     parent->cd();
1015                     TLatex tx; tx.SetNDC(); tx.SetTextSize(0.08);
1016                     tx.DrawLatex(0.15,0.85,bl.c_str()); continue;
1017                 }
1018 
1019                 // Clone so each pad owns its histograms independently
1020                 TH1D* hUc=(TH1D*)hU->Clone(std::format("hUc_{}_{}",k,ft).c_str());
1021                 TH1D* hTc=(TH1D*)hT->Clone(std::format("hTc_{}_{}",k,ft).c_str());
1022                 StyleLine(hTc,kColTruth,kMarTruth);
1023                 StyleLine(hUc,kColUnfolded,kMarUnfolded);
1024 
1025                 double ymax=std::max(hTc->GetMaximum(),hUc->GetMaximum());
1026                 double ymin=1e30;
1027                 for (int b=1;b<=hTc->GetNbinsX();++b)
1028                     if (hTc->GetBinContent(b)>0) ymin=std::min(ymin,hTc->GetBinContent(b));
1029                 if (ymin>ymax) ymin=ymax*1e-4;
1030 
1031                 // Manual split pads with unique names — TRatioPlot is not used
1032                 // because the first pair weight bin starts at 0, which makes
1033                 // log-x impossible for TRatioPlot internally.
1034                 // Both pads use the same histogram bin edges so axes align exactly.
1035                 parent->cd();
1036                 TPad* pTop2 = new TPad(
1037                     std::format("pTop2_{}_{}_{}", k, iL, iS).c_str(), "",
1038                     0, 0.30, 1, 1.0);
1039                 TPad* pBot2 = new TPad(
1040                     std::format("pBot2_{}_{}_{}", k, iL, iS).c_str(), "",
1041                     0, 0.00, 1, 0.30);
1042                 pTop2->SetBottomMargin(0.02);
1043                 pTop2->SetTopMargin(0.05);
1044                 pBot2->SetTopMargin(0.02);
1045                 pBot2->SetBottomMargin(0.38);
1046                 pTop2->Draw(); pBot2->Draw();
1047 
1048                 pTop2->cd();
1049                 pTop2->SetLogy();
1050                 pTop2->SetLogx();
1051                 hFrameA->GetYaxis()->SetRangeUser(ymin*0.3, ymax*3.0);
1052                 hFrameA->GetYaxis()->SetTitleSize(0.08);
1053                 hFrameA->GetYaxis()->SetLabelSize(0.07);
1054                 hFrameA->GetYaxis()->SetTitle("Pairs");
1055                 hFrameA->GetXaxis()->SetLabelSize(0);
1056                 hFrameA->SetTitle("");
1057                 hFrameA->Draw();
1058                 hTc->Draw("pSAME"); hUc->Draw("p SAME");
1059                 TLatex tx; tx.SetNDC(); tx.SetTextSize(0.08); tx.SetTextAlign(13);
1060                 tx.DrawLatex(0.15, 0.88, bl.c_str());
1061 
1062                 pBot2->cd();
1063                 pBot2->SetLogx();
1064                 TH1D* hR = MakeRatioHist(hUc, hTc,
1065                     std::format("hPWBinR_{}_{}",k,ft).c_str());
1066                 hR->SetDirectory(0);
1067                 StyleLine(hR, kColUnfolded, kMarUnfolded);
1068                 hFrameB->GetYaxis()->SetRangeUser(0.8, 1.2);
1069                 hFrameB->GetYaxis()->SetTitle("Unf/Truth");
1070                 hFrameB->GetYaxis()->SetTitleSize(0.12);
1071                 hFrameB->GetYaxis()->SetTitleOffset(0.38);
1072                 hFrameB->GetYaxis()->SetLabelSize(0.10);
1073                 hFrameB->GetYaxis()->SetNdivisions(504);
1074                 hFrameB->GetXaxis()->SetTitle("Pair weight bin");
1075                 hFrameB->GetXaxis()->SetTitleSize(0.12);
1076                 hFrameB->GetXaxis()->SetLabelSize(0.10);
1077                 hFrameB->SetTitle("");
1078                 hFrameB->Draw();
1079                 hR->Draw("pSAME");
1080                 DrawRefLine(hR);
1081                 // Do not delete hR — ROOT needs it alive until SaveAs repaints
1082             }
1083             cPage->cd();
1084             TLatex ct; ct.SetNDC(); ct.SetTextSize(0.015); ct.SetTextAlign(22);
1085             ct.DrawLatex(0.5,0.002,std::format("#Delta#phi: {:.2f}-{:.2f} (bin {})",
1086                 dPhiBins[k],dPhiBins[k+1],k).c_str());
1087             cPage->Update();  // flush all TRatioPlot pads before saving
1088             std::string pageOpt=(k==0)?(pdfName+"("):((k==nDphi-1)?(pdfName+")"):pdfName);
1089             cPage->SaveAs(pageOpt.c_str());
1090             cPage->Close(); delete cPage;
1091         }
1092 
1093         // Cleanup pair weight histograms
1094         for (int k=0; k<nDphi; ++k) {
1095             delete hPWUnf[k]; delete hPWTru[k];
1096             for (int ft=0; ft<nTrueFlat(); ++ft) {
1097                 delete hPWUnfBin[k*nTrueFlat()+ft];
1098                 delete hPWTruBin[k*nTrueFlat()+ft];
1099             }
1100         }
1101     }
1102 
1103 
1104     // ════════════════════════════════════════════════════════════════════════
1105     //  Response matrix visualization
1106     //  Two plots, both read from fResp:
1107     //
1108     //  respMatrix_pairWeight-{label}.png
1109     //    Grid of 8 representative ΔΦ bins showing the pair weight response:
1110     //    the nPairWeight × nPairWeight matrix obtained by marginalizing the
1111     //    3D response over all (iL, iS) dijet pT bins.  Reco pair weight on X,
1112     //    truth pair weight on Y.  Column-normalized so each reco bin sums to 1,
1113     //    making the migration probability visible.
1114     //
1115     //  respMatrix_dijetPt-{label}.png
1116     //    The nFlat × nFlat = 49 × 49 dijet pT response matrix (lead pT × subl pT
1117     //    in flat index) summed over all ΔΦ bins and pair weight bins.
1118     //    Column-normalized.  Shows the jet energy scale migration.
1119     // ════════════════════════════════════════════════════════════════════════
1120     if (fResp) {
1121         // ── Full 3D response matrix grid: all ΔΦ bins ────────────────────
1122         // Each pad shows the full nFlat3D×nFlat3D response matrix for one
1123         // ΔΦ bin, with reco flat3D index on X and truth flat3D index on Y.
1124         // Saved as PDF for lossless high-resolution rendering.
1125         const int nCols3D = 4, nRows3D = (nDphi + nCols3D - 1) / nCols3D;
1126         const int pad3D = 400;
1127         TCanvas* cResp = new TCanvas("cResp", "Full 3D response matrix",
1128                                      nCols3D*pad3D, nRows3D*pad3D);
1129         cResp->Divide(nCols3D, nRows3D, 0, 0);
1130 
1131         //vector<TH2D*> hDraw(nDphi, nullptr);
1132 
1133         for (int k = 0; k < nDphi; ++k) {
1134             cResp->cd(k+1);
1135             gPad->SetLeftMargin(0.12); gPad->SetBottomMargin(0.12);
1136             gPad->SetRightMargin(0.15); gPad->SetTopMargin(0.10);
1137             gPad->SetLogz();
1138 
1139             RooUnfoldResponse* resp =
1140                 (RooUnfoldResponse*)fResp->Get(
1141                     std::format("response_wEEC3D_{}", k).c_str());
1142             if (!resp) continue;
1143 
1144             TH2D* hFull = (TH2D*)resp->Hresponse();
1145             if (!hFull || hFull->Integral() == 0) continue;
1146 
1147             // Clone so we can set title/style without modifying the stored object
1148             TH2D *hDraw = (TH2D*)hFull->Clone(
1149                 std::format("hResp3D_{}",k).c_str());
1150             hDraw->SetDirectory(0);
1151             hDraw->SetTitle(
1152                 std::format("{:.2f}<#Delta#phi<{:.2f};Reco flat3D;Truth flat3D",
1153                     dPhiBins[k], dPhiBins[k+1]).c_str());
1154             hDraw->GetXaxis()->SetTitleSize(0.05);
1155             hDraw->GetYaxis()->SetTitleSize(0.05);
1156             hDraw->GetXaxis()->SetLabelSize(0.04);
1157             hDraw->GetYaxis()->SetLabelSize(0.04);
1158             hDraw->GetZaxis()->SetLabelSize(0.04);
1159             //double zMax = hDraw->GetMaximum();
1160             //if (zMax > 0) hDraw->GetZaxis()->SetRangeUser(zMax*1e-4, zMax);
1161             hDraw->Draw("COLZ");
1162 
1163             TLatex tex; tex.SetNDC(); tex.SetTextSize(0.06); tex.SetTextAlign(13);
1164             tex.DrawLatex(0.14, 0.92,
1165                 std::format("{:.2f}-{:.2f}",
1166                     dPhiBins[k], dPhiBins[k+1]).c_str());
1167             //delete hDraw;
1168         }
1169 
1170         cResp->SaveAs(std::format("{}/respMatrix_full3D-{}.pdf",
1171                                   plotDir, label).c_str());
1172 
1173         cResp->Close(); delete cResp;
1174 
1175         if(doTrim)
1176         {
1177             TCanvas* cRespTrimmed = new TCanvas("cRespTrimmed", "Full 3D response matrix Trimmed",
1178                                         nCols3D*pad3D, nRows3D*pad3D);
1179             cRespTrimmed->Divide(nCols3D, nRows3D, 0, 0);
1180 
1181             //vector<TH2D*> hDraw(nDphi, nullptr);
1182 
1183             for (int k = 0; k < nDphi; ++k) {
1184                 cRespTrimmed->cd(k+1);
1185                 gPad->SetLeftMargin(0.12); gPad->SetBottomMargin(0.12);
1186                 gPad->SetRightMargin(0.15); gPad->SetTopMargin(0.10);
1187                 gPad->SetLogz();
1188 
1189                 RooUnfoldResponse* resp =
1190                     (RooUnfoldResponse*)fWEEC->Get(
1191                         std::format("response_wEEC3D_trimmed_{}", k).c_str());
1192                 if (!resp) continue;
1193 
1194                 TH2D* hFull = (TH2D*)resp->Hresponse();
1195                 if (!hFull || hFull->Integral() == 0) continue;
1196 
1197                 // Clone so we can set title/style without modifying the stored object
1198                 TH2D *hDraw = (TH2D*)hFull->Clone(
1199                     std::format("hResp3D_trimmed_{}",k).c_str());
1200                 hDraw->SetDirectory(0);
1201                 hDraw->SetTitle(
1202                     std::format("{:.2f}<#Delta#phi<{:.2f};Reco flat3D;Truth flat3D",
1203                         dPhiBins[k], dPhiBins[k+1]).c_str());
1204                 hDraw->GetXaxis()->SetTitleSize(0.05);
1205                 hDraw->GetYaxis()->SetTitleSize(0.05);
1206                 hDraw->GetXaxis()->SetLabelSize(0.04);
1207                 hDraw->GetYaxis()->SetLabelSize(0.04);
1208                 hDraw->GetZaxis()->SetLabelSize(0.04);
1209                 //double zMax = hDraw->GetMaximum();
1210                 //if (zMax > 0) hDraw->GetZaxis()->SetRangeUser(zMax*1e-4, zMax);
1211                 hDraw->Draw("COLZ");
1212 
1213                 TLatex tex; tex.SetNDC(); tex.SetTextSize(0.06); tex.SetTextAlign(13);
1214                 tex.DrawLatex(0.14, 0.92,
1215                     std::format("{:.2f}-{:.2f}",
1216                         dPhiBins[k], dPhiBins[k+1]).c_str());
1217                 //delete hDraw;
1218             }
1219 
1220             cRespTrimmed->SaveAs(std::format("{}/respMatrixTrimmed_full3D-{}.pdf",
1221                                     plotDir, label).c_str());
1222         }
1223 
1224         // ── Cleanup ───────────────────────────────────────────────────────────
1225         fWEEC->Close(); delete fWEEC;
1226 
1227         /*
1228         // ── dijet pT response: sum over all ΔΦ and pair weight bins ──────
1229         TH2D* hJetResp = new TH2D("hJetResp",
1230             "Dijet pT response (from wEEC resp, marginalized);"
1231             "Reco flat dijet pT bin;Truth flat dijet pT bin",
1232             nRecoFlat(), 0, nRecoFlat(), nTrueFlat(), 0, nTrueFlat());
1233         hJetResp->SetDirectory(0);
1234 
1235         for (int k=0; k<nDphi; ++k) {
1236             RooUnfoldResponse* resp =
1237                 (RooUnfoldResponse*)fResp->Get(
1238                     std::format("response_wEEC3D_{}", k).c_str());
1239             if (!resp) continue;
1240             TH2D* hFull = (TH2D*)resp->Hresponse();
1241             if (!hFull || hFull->Integral() == 0) continue;
1242 
1243             for (int ftR=0; ftR<nRecoFlat(); ++ftR)
1244             for (int ftT=0; ftT<nTrueFlat(); ++ftT) {
1245                 double sum = 0;
1246                 for (int rIPw=0; rIPw<nPairWeight; ++rIPw)
1247                 for (int tIPw=0; tIPw<nPairWeight; ++tIPw) {
1248                     int rFlat = ftR * nPairWeight + rIPw;
1249                     int tFlat = ftT * nPairWeight + tIPw;
1250                     sum += hFull->GetBinContent(rFlat+1, tFlat+1);
1251                 }
1252                 if (sum > 0)
1253                     hJetResp->AddBinContent(
1254                         hJetResp->GetBin(ftR+1, ftT+1), sum);
1255             }
1256         }
1257 
1258         // Column-normalize
1259         for (int rBin=1; rBin<=nRecoFlat(); ++rBin) {
1260             double colSum=0;
1261             for (int tBin=1; tBin<=nTrueFlat(); ++tBin)
1262                 colSum += hJetResp->GetBinContent(rBin,tBin);
1263             if (colSum>0)
1264                 for (int tBin=1; tBin<=nTrueFlat(); ++tBin)
1265                     hJetResp->SetBinContent(rBin,tBin,
1266                         hJetResp->GetBinContent(rBin,tBin)/colSum);
1267         }
1268 
1269         TCanvas* cJR = new TCanvas("cJR","Dijet pT response (marginalized)",800,700);
1270         cJR->cd(); gPad->SetLogz();
1271         gPad->SetLeftMargin(0.12); gPad->SetBottomMargin(0.12);
1272         gPad->SetRightMargin(0.15); gPad->SetTopMargin(0.08);
1273         double jrMax = hJetResp->GetMaximum();
1274         if (jrMax > 0) hJetResp->GetZaxis()->SetRangeUser(jrMax*1e-3, jrMax);
1275         hJetResp->GetXaxis()->SetTitle("Reco flat dijet pT bin");
1276         hJetResp->GetYaxis()->SetTitle("Truth flat dijet pT bin");
1277         hJetResp->Draw("COLZ");
1278         cJR->SaveAs(std::format("{}/respMatrix_dijetPt-{}.png",
1279                                 plotDir,label).c_str());
1280         cJR->Close(); delete cJR;
1281         delete hJetResp;
1282         */
1283 
1284         RooUnfoldResponse *jetResp = (RooUnfoldResponse*)fResp->Get("response_jet_pT");
1285         TH2D *hJetResp = (TH2D*)jetResp->Hresponse();
1286         
1287         TCanvas* cJR = new TCanvas("cJR","Dijet pT response",800,700);
1288         cJR->cd(); gPad->SetLogz();
1289         gPad->SetLeftMargin(0.12); gPad->SetBottomMargin(0.12);
1290         gPad->SetRightMargin(0.15); gPad->SetTopMargin(0.08);
1291         double jrMax = hJetResp->GetMaximum();
1292         if (jrMax > 0) hJetResp->GetZaxis()->SetRangeUser(jrMax*1e-5, jrMax);
1293         hJetResp->GetXaxis()->SetTitle("Reco flat dijet pT bin");
1294         hJetResp->GetYaxis()->SetTitle("Truth flat dijet pT bin");
1295         hJetResp->Draw("COLZ");
1296         cJR->SaveAs(std::format("{}/respMatrix_dijetPt-{}.png",
1297                                 plotDir,label).c_str());
1298         cJR->Close(); delete cJR;
1299         delete hJetResp;
1300         delete jetResp;
1301 
1302         fResp->Close(); delete fResp;
1303     }
1304 
1305     for (int ft=0; ft<nTrueFlat(); ++ft) {
1306         delete hTruthBin[ft];
1307         delete hUnfBin[ft];
1308     }
1309     for (int ft=0; ft<nRecoFlat(); ++ft)
1310         delete hMeasBin[ft];
1311 
1312     std::cout << "Done.\n";
1313     gSystem->Exit(0);
1314 }
[ 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!