macros/wEEC_3D_unfolding_kinematics/draw_3DClosure.C

0001 #include "analysisHelper.h"
0002 #include "RooUnfoldResponse.h"
0003
0004 void draw_3DClosure(const char* outDir   = "/sphenix/tg/tg01/jets/bkimelman/VandyDSTs_wEEC_3D_unfolding_kinematics_Apr26_2026/")
0005 {
0006     gROOT->SetBatch(true);
0007     gStyle->SetOptStat(0); gStyle->SetOptTitle(0);
0008     gStyle->SetPaintTextFormat(".2f");
0009     const std::string label   = ModeLabel(Mode::kFull);
0010
0011     const std::string plotDir = std::format("{}/Plots", outDir);
0012
0013     std::string wEECFile  = std::format("{}/wEEC-{}.root",   outDir, label);
0014
0015     TFile* fWEEC = new TFile(wEECFile.c_str(),"READ");
0016     if (!fWEEC || fWEEC->IsZombie()) { std::cerr<<"Cannot open "<<wEECFile<<"\n"; return; }
0017
0018     // Open response file immediately — needed for truth reference throughout
0019     TFile* fResp = new TFile(std::format("{}/response-all-fullClosure.root",outDir).c_str(),"READ");
0020
0021     std::vector<TH1D*> unfHist(nDphi, nullptr);
0022     std::vector<TH1D*> truthHist(nDphi, nullptr);
0023     std::vector<TH1D*> measHist(nDphi, nullptr);
0024     std::vector<TH1D*> respTruthHist(nDphi, nullptr);
0025     std::vector<TH1D*> respMeasHist(nDphi, nullptr);
0026
0027     for(int i=0; i<nDphi; i++)
0028     {
0029         unfHist[i] = (TH1D*)fWEEC->Get(std::format("hWEEC3D_unfolded_{}",i).c_str());
0030         truthHist[i] = (TH1D*)fResp->Get(std::format("hWEEC3D_truth_{}",i).c_str());
0031         measHist[i] = (TH1D*)fResp->Get(std::format("hWEEC3D_meas_{}",i).c_str());
0032         RooUnfoldResponse *resp = (RooUnfoldResponse*)fResp->Get(std::format("response_wEEC3D_{}",i).c_str());
0033         respTruthHist[i] = (TH1D*)resp->Htruth();
0034         respMeasHist[i] = (TH1D*)resp->Hmeasured();
0035         delete resp;
0036     }
0037
0038     TCanvas *c1 = new TCanvas("c1","",8*592,4*592);
0039
0040     //std::pair<vector<TH1D*>, vector<TH1D*>> combos[3] = {std::make_pair(unfHist, truthHist), std::make_pair(truthHist, respTruthHist), std::make_pair(measHist, respMeasHist)};
0041     std::string names[3] = {"unfolded_truth", "truth_respTruth", "meas_respMeas"};
0042
0043     for(int i=0; i<3; i++)
0044     {
0045         c1->Clear();
0046         c1->Divide(8,4);
0047
0048         vector<TH1D*> rat(nDphi, nullptr);
0049         vector<TGraph*> gr(nDphi, nullptr);
0050         vector<TGraph*> grG(nDphi, nullptr);
0051
0052         for(int j=0; j<nDphi; j++)
0053         {
0054             c1->cd(j+1);
0055             c1->cd(j+1)->SetLogy(0);
0056             if(i == 0)
0057             {
0058                 rat[j] = (TH1D*)unfHist[j]->Clone();
0059                 rat[j]->Divide(truthHist[j]);
0060             }
0061             else if(i == 1)
0062             {
0063                 rat[j] = (TH1D*)truthHist[j]->Clone();
0064                 rat[j]->Divide(respTruthHist[j]);
0065             }
0066              else if(i == 2)
0067             {
0068                 rat[j] = (TH1D*)measHist[j]->Clone();
0069                 rat[j]->Divide(respMeasHist[j]);
0070             }
0071             rat[j]->SetMarkerStyle(20);
0072
0073
0074             gr[j] = new TGraph();
0075             grG[j] = new TGraph();
0076             double minV = 1.0;
0077             double maxV = 1.0;
0078
0079             for(int k=1; k<=rat[j]->GetNbinsX(); k++)
0080             {
0081                 double bc = rat[j]->GetBinContent(k);
0082                 if(bc == 0.0) continue;
0083
0084                 if(bc < minV) minV = bc;
0085                 if(bc > maxV) maxV = bc;
0086
0087                 if(bc == 1.0)
0088                 {
0089                     grG[j]->AddPoint(rat[j]->GetXaxis()->GetBinCenter(k), bc);
0090                 }
0091                 else
0092                 {
0093                     gr[j]->AddPoint(rat[j]->GetXaxis()->GetBinCenter(k), bc);
0094                 }
0095             }
0096             gr[j]->SetMarkerStyle(21);
0097             gr[j]->SetMarkerColor(kRed);
0098
0099             grG[j]->SetMarkerStyle(22);
0100             grG[j]->SetMarkerColor(kGreen+2);
0101
0102             rat[j]->GetYaxis()->SetRangeUser(0.75*minV, 1.25*maxV);
0103             rat[j]->Draw("P");
0104             if(gr[j]->GetN()>0) gr[j]->Draw("PSAME");
0105             if(grG[j]->GetN()>0) grG[j]->Draw("PSAME");
0106         }
0107
0108         c1->SaveAs(std::format("{}/ratio_{}.pdf",plotDir,names[i]).c_str());
0109     }
0110
0111
0112
0113 }