macros/wEEC_3D_unfolding_kinematics/makeVtxWeights.C

0001 #include "analysisHelper.h"
0002
0003 // ─────────────────────────────────────────────────────────────────────────────
0004 //  makeVtxWeights.C
0005 //
0006 //  Forms the vtx_z reweighting factors needed to correct the MC vtx_z
0007 //  distribution to match data.
0008 //
0009 //  Inputs:
0010 //    mcFile   — merged MC file (output of hadd over all kVtx response files)
0011 //               must contain hVtxMC
0012 //    dataFile — merged data file (output of hadd over all fillData.C outputs)
0013 //               must contain hVtxData
0014 //
0015 //  Procedure:
0016 //    1. Load hVtxMC and hVtxData from the two merged files.
0017 //    2. Normalize both to unit area within [-10, 10) cm.
0018 //    3. Divide: hVtxWeight = hVtxData_norm / hVtxMC_norm.
0019 //    4. Print the bin contents as a C++ array to paste into analysisHelper.h.
0020 //    5. Write all histograms and a comparison canvas to outFile.
0021 //
0022 //  Usage:
0023 //    root -l 'makeVtxWeights.C("vtxWeights.root",
0024 //                               "mc_vtx_merged.root",
0025 //                               "data_merged.root")'
0026 //
0027 //  After running, replace the vtxWeights[] array in analysisHelper.h with
0028 //  the printed output, then rerun fillResponse.C in kData mode.
0029 // ─────────────────────────────────────────────────────────────────────────────
0030
0031 void makeVtxWeights(const char* outDir  = "vtxWeights.root")
0032 {
0033
0034     std::string MCFile = std::format("{}/response-all-vtx.root",outDir);
0035     std::string dataFile = std::format("{}/data_measured-all.root",outDir);
0036     std::string outFile = std::format("{}/vtxWeights.root",outDir);
0037
0038     // ── load histograms ───────────────────────────────────────────────────
0039     TFile* fMC = TFile::Open(MCFile.c_str(), "READ");
0040     if (!fMC || fMC->IsZombie()) {
0041         std::cerr << "Cannot open MC file: " << MCFile << "\n"; return; }
0042
0043     TH1D* hMC = (TH1D*)fMC->Get("hVtxMC");
0044     if (!hMC) {
0045         std::cerr << "hVtxMC not found in " << MCFile << "\n"; return; }
0046     hMC = (TH1D*)hMC->Clone("hVtxMC_sum"); hMC->SetDirectory(0);
0047     fMC->Close(); delete fMC;
0048
0049     TFile* fData = TFile::Open(dataFile.c_str(), "READ");
0050     if (!fData || fData->IsZombie()) {
0051         std::cerr << "Cannot open data file: " << dataFile << "\n"; return; }
0052
0053     TH1D* hData = (TH1D*)fData->Get("hVtxData");
0054     if (!hData) {
0055         std::cerr << "hVtxData not found in " << dataFile << "\n"; return; }
0056     hData = (TH1D*)hData->Clone("hVtxData_sum"); hData->SetDirectory(0);
0057     fData->Close(); delete fData;
0058
0059     if (hMC  ->Integral() <= 0) { std::cerr << "MC vtx histogram is empty\n";   return; }
0060     if (hData->Integral() <= 0) { std::cerr << "Data vtx histogram is empty\n"; return; }
0061
0062     // ── normalize both to unit area within [-10, 10) cm ──────────────────
0063     // Use Integral() which sums all bins — since both histograms are already
0064     // restricted to this range, this is the full integral.
0065     TH1D* hMC_norm   = (TH1D*)hMC  ->Clone("hVtxMC_norm");
0066     TH1D* hData_norm = (TH1D*)hData->Clone("hVtxData_norm");
0067     hMC_norm  ->Scale(1.0 / hMC_norm  ->Integral());
0068     hData_norm->Scale(1.0 / hData_norm->Integral());
0069
0070     // ── form ratio ────────────────────────────────────────────────────────
0071     TH1D* hWeight = (TH1D*)hData_norm->Clone("hVtxWeight");
0072     hWeight->SetTitle("vtx_z reweight: data/MC;v_{tz} (cm);data/MC");
0073
0074     for (int i = 1; i <= nVtxBins; ++i) {
0075         double mc  = hMC_norm  ->GetBinContent(i);
0076         double dat = hData_norm->GetBinContent(i);
0077         double mc_err  = hMC_norm  ->GetBinError(i);
0078         double dat_err = hData_norm->GetBinError(i);
0079         if (mc > 0) {
0080             double w   = dat / mc;
0081             double err = w * std::sqrt(std::pow(dat_err/dat, 2) +
0082                                        std::pow(mc_err /mc,  2));
0083             hWeight->SetBinContent(i, w);
0084             hWeight->SetBinError  (i, err);
0085         } else {
0086             // MC has no entries in this bin — weight undefined, set to 1
0087             // (should not occur within [-10,10) cm for a large MC sample).
0088             std::cerr << std::format(
0089                 "Warning: MC vtx bin {} ([{:.0f},{:.0f}) cm) is empty — weight set to 1\n",
0090                 i, vtxZMin + (i-1)*vtxBinWidth, vtxZMin + i*vtxBinWidth);
0091             hWeight->SetBinContent(i, 1.0);
0092             hWeight->SetBinError  (i, 0.0);
0093         }
0094     }
0095
0096     // ── print array for analysisHelper.h ─────────────────────────────────
0097     std::cout << "\n";
0098     std::cout << "// ── paste into analysisHelper.h vtxWeights[] ──────────\n";
0099     std::cout << "constexpr double vtxWeights[nVtxBins] = {\n";
0100     for (int i = 1; i <= nVtxBins; ++i) {
0101         double w   = hWeight->GetBinContent(i);
0102         double lo  = vtxZMin + (i-1) * vtxBinWidth;
0103         double hi  = lo + vtxBinWidth;
0104         std::cout << std::format("    {:.6f},  // bin {:2d}: [{:.0f}, {:.0f}) cm\n",
0105             w, i-1, lo, hi);
0106     }
0107     std::cout << "};\n\n";
0108
0109     // ── draw comparison ───────────────────────────────────────────────────
0110     gStyle->SetOptStat(0);
0111     TCanvas* c = new TCanvas("cVtx", "vtx_z distributions", 1200, 500);
0112     c->Divide(3, 1);
0113
0114     c->cd(1);
0115     hMC_norm  ->SetLineColor(kBlue+1);   hMC_norm  ->SetLineWidth(2);
0116     hData_norm->SetLineColor(kRed+1);    hData_norm->SetLineWidth(2);
0117     double ymax = std::max(hMC_norm->GetMaximum(), hData_norm->GetMaximum()) * 1.15;
0118     hMC_norm->GetYaxis()->SetRangeUser(0, ymax);
0119     hMC_norm->SetTitle("vtx_z: MC vs data (normalized);v_{tz} (cm);a.u.");
0120     hMC_norm->Draw("HIST");
0121     hData_norm->Draw("E SAME");
0122     TLegend* leg = new TLegend(0.15, 0.75, 0.55, 0.90);
0123     leg->AddEntry(hMC_norm,   "MC (xsec weighted)", "l");
0124     leg->AddEntry(hData_norm, "Data",               "lep");
0125     leg->Draw();
0126
0127     c->cd(2);
0128     hWeight->SetLineColor(kBlack); hWeight->SetLineWidth(2);
0129     hWeight->SetMarkerStyle(20);
0130     hWeight->GetYaxis()->SetRangeUser(0, 2);
0131     hWeight->Draw("E");
0132     TLine* unity = new TLine(vtxZMin, 1, vtxZMax, 1);
0133     unity->SetLineStyle(2); unity->SetLineColor(kGray+1); unity->Draw();
0134
0135     c->cd(3);
0136     // Show MC before and after reweighting as a closure check
0137     TH1D* hMC_rw = (TH1D*)hMC_norm->Clone("hVtxMC_reweighted");
0138     for (int i = 1; i <= nVtxBins; ++i)
0139         hMC_rw->SetBinContent(i,
0140             hMC_norm->GetBinContent(i) * hWeight->GetBinContent(i));
0141     hMC_rw->Scale(1.0 / hMC_rw->Integral());
0142     hMC_rw->SetLineColor(kGreen+2); hMC_rw->SetLineWidth(2);
0143     hMC_rw->SetTitle("Closure: reweighted MC vs data;v_{tz} (cm);a.u.");
0144     hData_norm->GetYaxis()->SetRangeUser(0, ymax);
0145     hData_norm->Draw("E");
0146     hMC_rw->Draw("HIST SAME");
0147     TLegend* leg2 = new TLegend(0.15, 0.75, 0.65, 0.90);
0148     leg2->AddEntry(hData_norm, "Data",              "lep");
0149     leg2->AddEntry(hMC_rw,    "MC (reweighted)",   "l");
0150     leg2->Draw();
0151
0152     // ── write output ──────────────────────────────────────────────────────
0153     TFile* fOut = new TFile(outFile.c_str(), "RECREATE");
0154     hMC      ->Write("hVtxMC");
0155     hData    ->Write("hVtxData");
0156     hMC_norm ->Write("hVtxMC_norm");
0157     hData_norm->Write("hVtxData_norm");
0158     hWeight  ->Write("hVtxWeight");
0159     hMC_rw   ->Write("hVtxMC_reweighted");
0160     c->Write("cVtxComparison");
0161     fOut->Close();
0162
0163     std::cout << "Written: " << outFile << "\n";
0164 }