systematics/trackMomScale/compare_mass_ptbins.C

0001 // compare_mass_ptbins.C
0002 // Shows invariant mass distributions in each smearing pT bin for:
0003 //   1. Data                       (cuts pre-applied)
0004 //   2. SV sim filtered            (cuts pre-applied)
0005 //   3. SV sim filtered + smeared  (cuts pre-applied)
0006 // All histograms normalised to unit area within each bin.
0007 // Produces one canvas per particle (K0S and Lambda), 4x2 layout:
0008 //   7 pads = 7 pT bins,  8th pad = shared legend.
0009 //
0010 // Run with:  root -l -b -q compare_mass_ptbins.C
0011
0012 #include "TFile.h"
0013 #include "TTree.h"
0014 #include "TH1F.h"
0015 #include "TCanvas.h"
0016 #include "TLegend.h"
0017 #include "TLatex.h"
0018 #include "TPaveText.h"
0019 #include "TLine.h"
0020 #include "TStyle.h"
0021 #include <cmath>
0022 #include <ctime>
0023 #include <sstream>
0024 #include <iostream>
0025
0026 // ── pT bins (same as used in smearing systematic) ────────────────────────────
0027 const int    NBINS   = 7;
0028 const double PT_LO[] = {0.5, 0.8, 1.1, 1.4, 1.8, 2.2, 3.0};
0029 const double PT_HI[] = {0.8, 1.1, 1.4, 1.8, 2.2, 3.0, 4.0};
0030
0031 // ── date ─────────────────────────────────────────────────────────────────────
0032 std::string getDateStr()
0033 {
0034     std::time_t t = std::time(0); std::tm* n = std::localtime(&t);
0035     std::stringstream ss;
0036     ss << (n->tm_mon+1) << '/' << n->tm_mday << '/' << (n->tm_year+1900);
0037     return ss.str();
0038 }
0039
0040 // ── style ─────────────────────────────────────────────────────────────────────
0041 void hStyle(TH1F* h, Color_t col, Style_t ls, float lw=2)
0042 {
0043     h->SetLineColor(col); h->SetLineWidth(lw); h->SetLineStyle(ls);
0044     h->SetFillStyle(0);
0045 }
0046
0047 // ── fill mass-in-pT-bins from any file (no cuts applied) ─────────────────────
0048 void fillMassFiltered(TH1F* h[], const char* fname,
0049                       const char* pt_br, const char* mass_br)
0050 {
0051     TFile* f = TFile::Open(fname,"READ");
0052     if (!f || f->IsZombie()) { std::cerr << "Cannot open " << fname << "\n"; return; }
0053     TTree* t = (TTree*)f->Get("DecayTree");
0054     if (!t) { f->Close(); return; }
0055
0056     float pt, mass;
0057     t->SetBranchAddress(pt_br,   &pt);
0058     t->SetBranchAddress(mass_br, &mass);
0059
0060     Long64_t N = t->GetEntries();
0061     for (Long64_t i = 0; i < N; ++i) {
0062         t->GetEntry(i);
0063         float m_mev = mass * 1000.f;
0064         for (int b = 0; b < NBINS; ++b)
0065             if (pt >= PT_LO[b] && pt < PT_HI[b]) { h[b]->Fill(m_mev); break; }
0066     }
0067     f->Close();
0068     std::cout << "  " << fname << " : " << N << " events\n";
0069 }
0070
0071 // ── draw one pad ─────────────────────────────────────────────────────────────
0072 void drawBinPad(TVirtualPad* pad, TH1F* hD, TH1F* hF, TH1F* hS,
0073                 int bin_idx, const char* mass_title, bool fit_gauss = false)
0074 {
0075     pad->cd();
0076     pad->SetLeftMargin(0.14); pad->SetRightMargin(0.04);
0077     pad->SetBottomMargin(0.16); pad->SetTopMargin(0.08);
0078
0079     // normalise to unit area
0080     for (TH1F* h : {hD, hF, hS})
0081         if (h->Integral() > 0) h->Scale(1.0 / h->Integral());
0082
0083     hStyle(hD, kBlack,   1, 2);
0084     hStyle(hF, kAzure+7, 2, 2);
0085     hStyle(hS, kRed+1,   1, 2);
0086
0087     double ymax = 1.25 * std::max({hD->GetMaximum(), hF->GetMaximum(), hS->GetMaximum()});
0088     if (ymax == 0) ymax = 1.0;
0089
0090     hD->GetXaxis()->SetTitle(mass_title);
0091     hD->GetYaxis()->SetTitle("Norm. counts");
0092     hD->GetYaxis()->SetRangeUser(0, ymax);
0093     hD->GetXaxis()->SetTitleSize(0.07);
0094     hD->GetYaxis()->SetTitleSize(0.065);
0095     hD->GetXaxis()->SetLabelSize(0.065);
0096     hD->GetYaxis()->SetLabelSize(0.065);
0097     hD->GetXaxis()->SetTitleOffset(1.0);
0098     hD->GetYaxis()->SetTitleOffset(0.95);
0099     hD->GetYaxis()->SetNdivisions(505);
0100
0101     hD->Draw("HIST");
0102     hF->Draw("HIST SAME");
0103     hS->Draw("HIST SAME");
0104
0105     // pT bin label
0106     TLatex tex; tex.SetNDC(); tex.SetTextSize(0.065); tex.SetTextFont(42);
0107     char label[64];
0108     if (PT_HI[bin_idx] < 9)
0109         snprintf(label, sizeof(label), "%.1f < p_{T} < %.1f GeV", PT_LO[bin_idx], PT_HI[bin_idx]);
0110     else
0111         snprintf(label, sizeof(label), "p_{T} > %.1f GeV", PT_LO[bin_idx]);
0112     tex.DrawLatex(0.18, 0.82, label);
0113
0114     if (fit_gauss) {
0115         const double FIT_LO = 475, FIT_HI = 520;
0116         struct Entry { TH1F* h; Color_t col; const char* tag; };
0117         Entry entries[] = {
0118             {hD, kBlack,   "Data"},
0119             {hF, kAzure+7, "Sim "},
0120             {hS, kRed+1,   "Smr "}
0121         };
0122         TLatex stex; stex.SetNDC(); stex.SetTextFont(42); stex.SetTextSize(0.058);
0123         double y = 0.71;
0124         for (auto& e : entries) {
0125             if (e.h->Integral() <= 0) { y -= 0.095; continue; }
0126             TF1* g = new TF1(Form("gfit_%d_%s", bin_idx, e.tag),
0127                              "gaus", FIT_LO, FIT_HI);
0128             e.h->Fit(g, "RQ0");
0129             stex.SetTextColor(e.col);
0130             char buf[64];
0131             snprintf(buf, sizeof(buf), "%s #sigma = %.1f MeV", e.tag, g->GetParameter(2));
0132             stex.DrawLatex(0.18, y, buf);
0133             y -= 0.095;
0134         }
0135     }
0136 }
0137
0138 // ── draw legend pad ───────────────────────────────────────────────────────────
0139 void drawLegendPad(TVirtualPad* pad)
0140 {
0141     pad->cd();
0142     pad->SetLeftMargin(0.05); pad->SetRightMargin(0.05);
0143
0144     TLegend* leg = new TLegend(0.08, 0.28, 0.92, 0.72);
0145     leg->SetBorderSize(0); leg->SetFillStyle(0); leg->SetTextSize(0.075);
0146     // dummy hists for legend entries
0147     TH1F* dD = new TH1F("dD","",1,0,1); hStyle(dD, kBlack,   1, 2);
0148     TH1F* dF = new TH1F("dF","",1,0,1); hStyle(dF, kAzure+7, 2, 2);
0149     TH1F* dS = new TH1F("dS","",1,0,1); hStyle(dS, kRed+1,   1, 2);
0150     leg->AddEntry(dD, "Data",                    "l");
0151     leg->AddEntry(dF, "SV sim (filtered)",       "l");
0152     leg->AddEntry(dS, "SV sim (filt. + smeared)","l");
0153     leg->Draw();
0154
0155     // sPHENIX label
0156     TLatex tex; tex.SetNDC(); tex.SetTextFont(42); tex.SetTextSize(0.085);
0157     tex.DrawLatex(0.08, 0.85, "#it{#bf{sPHENIX}} Internal");
0158     tex.SetTextSize(0.075);
0159     tex.DrawLatex(0.08, 0.76, "#it{p}+#it{p}  #sqrt{s} = 200 GeV");
0160     tex.SetTextColor(kGray+2); tex.SetTextSize(0.070);
0161     tex.DrawLatex(0.08, 0.18, getDateStr().c_str());
0162 }
0163
0164 // ── build and save one full canvas ───────────────────────────────────────────
0165 void makeParticleCanvas(
0166     TH1F* hD[], TH1F* hF[], TH1F* hS[],
0167     const char* cname, const char* save_base, const char* mass_title,
0168     bool fit_gauss = false)
0169 {
0170     TCanvas* c = new TCanvas(cname, cname, 1600, 800);
0171     c->Divide(4, 2, 0.003, 0.003);
0172
0173     for (int b = 0; b < NBINS; ++b)
0174         drawBinPad(c->GetPad(b+1), hD[b], hF[b], hS[b], b, mass_title, fit_gauss);
0175
0176     drawLegendPad(c->GetPad(8));
0177
0178     TString pdf = TString(save_base) + ".pdf";
0179     TString png = TString(save_base) + ".png";
0180     c->SaveAs(pdf); c->SaveAs(png);
0181     std::cout << "Saved " << save_base << ".pdf/.png\n";
0182 }
0183
0184 // ── main ──────────────────────────────────────────────────────────────────────
0185 void compare_mass_ptbins()
0186 {
0187     gStyle->SetOptStat(0);
0188     gStyle->SetOptTitle(0);
0189
0190     const char* KS_DATA   = "KShort6RunCombined.root";
0191     const char* KS_FILT   = "outputKFParticleKShortRecoSV_filtered.root";
0192     //const char* KS_SMEAR  = "outputKFParticleKShortRecoSV_filtered_smeared.root";
0193     const char* KS_SMEAR = "outputKFParticleKShortRecoSV_filtered_masssmeared.root";
0194     const char* LAM_DATA  = "Lambda6RunCombined.root";
0195     const char* LAM_FILT  = "outputKFParticleLambda0SV_filtered.root";
0196     const char* LAM_SMEAR = "outputKFParticleLambda0SV_filtered_smeared.root";
0197
0198     // ── K0S ─────────────────────────────────────────────────────────────────
0199     std::cout << "\n=== K0S mass in pT bins ===\n";
0200
0201     const double KS_LO = 450, KS_HI = 550;
0202     const int    KS_NB = 50;   // 2 MeV/bin
0203
0204     TH1F* ksD[NBINS], *ksF[NBINS], *ksS[NBINS];
0205     for (int b = 0; b < NBINS; ++b) {
0206         ksD[b] = new TH1F(Form("ksD%d",b), "", KS_NB, KS_LO, KS_HI);
0207         ksF[b] = new TH1F(Form("ksF%d",b), "", KS_NB, KS_LO, KS_HI);
0208         ksS[b] = new TH1F(Form("ksS%d",b), "", KS_NB, KS_LO, KS_HI);
0209     }
0210
0211     fillMassFiltered(ksD, KS_DATA,  "K_S0_pT", "K_S0_mass");
0212     fillMassFiltered(ksF, KS_FILT,  "K_S0_pT", "K_S0_mass");
0213     fillMassFiltered(ksS, KS_SMEAR, "K_S0_pT", "K_S0_mass");
0214
0215     makeParticleCanvas(ksD, ksF, ksS,
0216                        "cKS", "compare_ks_mass_ptbins",
0217                        "K_{S}^{0} mass (MeV/#it{c}^{2})", true);
0218
0219     // ── Lambda ───────────────────────────────────────────────────────────────
0220     std::cout << "\n=== Lambda mass in pT bins ===\n";
0221
0222     const double LAM_LO = 1090, LAM_HI = 1150;
0223     const int    LAM_NB = 60;   // 1 MeV/bin
0224
0225     TH1F* lmD[NBINS], *lmF[NBINS], *lmS[NBINS];
0226     for (int b = 0; b < NBINS; ++b) {
0227         lmD[b] = new TH1F(Form("lmD%d",b), "", LAM_NB, LAM_LO, LAM_HI);
0228         lmF[b] = new TH1F(Form("lmF%d",b), "", LAM_NB, LAM_LO, LAM_HI);
0229         lmS[b] = new TH1F(Form("lmS%d",b), "", LAM_NB, LAM_LO, LAM_HI);
0230     }
0231
0232     fillMassFiltered(lmD, LAM_DATA,  "Lambda0_pT", "Lambda0_mass");
0233     fillMassFiltered(lmF, LAM_FILT,  "Lambda0_pT", "Lambda0_mass");
0234     fillMassFiltered(lmS, LAM_SMEAR, "Lambda0_pT", "Lambda0_mass");
0235
0236     makeParticleCanvas(lmD, lmF, lmS,
0237                        "cLam", "compare_lam_mass_ptbins",
0238                        "#Lambda^{0} mass (MeV/#it{c}^{2})");
0239
0240     std::cout << "\nAll done.\n";
0241 }