systematics/trackMomScale/compare_daughter_pt_ptbins.C

0001 // compare_daughter_pt_ptbins.C
0002 // Shows daughter track pT distributions in each smearing pT bin for:
0003 //   1. Data                       (quality cuts 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 4 canvases (4x2 layout, 7 pT bins + legend pad):
0008 //   compare_ks_t1pt_ptbins   — K0S  track_1_pT (π⁺) in each K0S pT bin
0009 //   compare_ks_t2pt_ptbins   — K0S  track_2_pT (π⁻) in each K0S pT bin
0010 //   compare_lam_t1pt_ptbins  — Λ⁰   track_1_pT (π⁻) in each Λ⁰ pT bin
0011 //   compare_lam_t2pt_ptbins  — Λ⁰   track_2_pT (p)  in each Λ⁰ pT bin
0012 //
0013 // Run with:  root -l -b -q compare_daughter_pt_ptbins.C
0014
0015 #include "TFile.h"
0016 #include "TTree.h"
0017 #include "TH1F.h"
0018 #include "TCanvas.h"
0019 #include "TLegend.h"
0020 #include "TLatex.h"
0021 #include "TPaveText.h"
0022 #include "TStyle.h"
0023 #include <cmath>
0024 #include <ctime>
0025 #include <sstream>
0026 #include <iostream>
0027
0028 // ── pT bins ───────────────────────────────────────────────────────────────────
0029 const int    NBINS   = 7;
0030 const double PT_LO[] = {0.5, 0.8, 1.1, 1.4, 1.8, 2.2, 3.0};
0031 const double PT_HI[] = {0.8, 1.1, 1.4, 1.8, 2.2, 3.0, 4.0};
0032
0033 // ── daughter pT histogram range ───────────────────────────────────────────────
0034 const int    DPT_NB = 60;
0035 const double DPT_LO = 0.0, DPT_HI = 3.0;   // pion
0036 const double PRO_LO = 0.0, PRO_HI = 4.0;   // proton
0037
0038 // ── date ─────────────────────────────────────────────────────────────────────
0039 std::string getDateStr3()
0040 {
0041     std::time_t t = std::time(0); std::tm* n = std::localtime(&t);
0042     std::stringstream ss;
0043     ss << (n->tm_mon+1) << '/' << n->tm_mday << '/' << (n->tm_year+1900);
0044     return ss.str();
0045 }
0046
0047 // ── style ─────────────────────────────────────────────────────────────────────
0048 void hSty3(TH1F* h, Color_t col, Style_t ls, float lw=2)
0049 {
0050     h->SetLineColor(col); h->SetLineWidth(lw); h->SetLineStyle(ls);
0051     h->SetFillStyle(0);
0052 }
0053
0054 // ── fill daughter pT from any file (no cuts applied) ─────────────────────────
0055 // h1[b] = track_1_pT in parent-pT bin b
0056 // h2[b] = track_2_pT in parent-pT bin b
0057 void fillDptFiltered(TH1F* h1[], TH1F* h2[], const char* fname,
0058                      const char* parent_pt_br)
0059 {
0060     TFile* f = TFile::Open(fname,"READ");
0061     if (!f||f->IsZombie()) { std::cerr<<"Cannot open "<<fname<<"\n"; return; }
0062     TTree* t = (TTree*)f->Get("DecayTree");
0063     if (!t) { f->Close(); return; }
0064
0065     float ppt, dpt1, dpt2;
0066     t->SetBranchAddress(parent_pt_br,  &ppt);
0067     t->SetBranchAddress("track_1_pT",  &dpt1);
0068     t->SetBranchAddress("track_2_pT",  &dpt2);
0069
0070     Long64_t N = t->GetEntries();
0071     for (Long64_t i = 0; i < N; ++i) {
0072         t->GetEntry(i);
0073         for (int b = 0; b < NBINS; ++b) {
0074             if (ppt >= PT_LO[b] && ppt < PT_HI[b]) {
0075                 h1[b]->Fill(dpt1);
0076                 h2[b]->Fill(dpt2);
0077                 break;
0078             }
0079         }
0080     }
0081     f->Close();
0082     std::cout << "  " << fname << " : " << N << " events\n";
0083 }
0084
0085 // ── draw one pad ─────────────────────────────────────────────────────────────
0086 void drawDptPad(TVirtualPad* pad, TH1F* hD, TH1F* hF, TH1F* hS,
0087                 int bin_idx, const char* xtitle)
0088 {
0089     pad->cd();
0090     pad->SetLeftMargin(0.14); pad->SetRightMargin(0.04);
0091     pad->SetBottomMargin(0.16); pad->SetTopMargin(0.08);
0092
0093     for (TH1F* h : {hD, hF, hS})
0094         if (h->Integral() > 0) h->Scale(1.0 / h->Integral());
0095
0096     hSty3(hD, kBlack,   1, 2);
0097     hSty3(hF, kAzure+7, 2, 2);
0098     hSty3(hS, kRed+1,   1, 2);
0099
0100     double ymax = 1.25 * std::max({hD->GetMaximum(), hF->GetMaximum(), hS->GetMaximum()});
0101     if (ymax == 0) ymax = 1.0;
0102
0103     hD->GetXaxis()->SetTitle(xtitle);
0104     hD->GetYaxis()->SetTitle("Norm. counts");
0105     hD->GetYaxis()->SetRangeUser(0, ymax);
0106     hD->GetXaxis()->SetTitleSize(0.07);
0107     hD->GetYaxis()->SetTitleSize(0.065);
0108     hD->GetXaxis()->SetLabelSize(0.065);
0109     hD->GetYaxis()->SetLabelSize(0.065);
0110     hD->GetXaxis()->SetTitleOffset(1.0);
0111     hD->GetYaxis()->SetTitleOffset(0.95);
0112     hD->GetYaxis()->SetNdivisions(505);
0113
0114     hD->Draw("HIST");
0115     hF->Draw("HIST SAME");
0116     hS->Draw("HIST SAME");
0117
0118     TLatex tex; tex.SetNDC(); tex.SetTextSize(0.065); tex.SetTextFont(42);
0119     char label[64];
0120     snprintf(label, sizeof(label), "%.1f < p_{T} < %.1f GeV", PT_LO[bin_idx], PT_HI[bin_idx]);
0121     tex.DrawLatex(0.18, 0.82, label);
0122 }
0123
0124 // ── legend pad ────────────────────────────────────────────────────────────────
0125 void drawLegendPad3(TVirtualPad* pad, const char* track_label)
0126 {
0127     pad->cd();
0128     pad->SetLeftMargin(0.05); pad->SetRightMargin(0.05);
0129
0130     TLegend* leg = new TLegend(0.08, 0.28, 0.92, 0.72);
0131     leg->SetBorderSize(0); leg->SetFillStyle(0); leg->SetTextSize(0.075);
0132     TH1F* dD = new TH1F(Form("dD_%s",track_label),"",1,0,1); hSty3(dD,kBlack,  1,2);
0133     TH1F* dF = new TH1F(Form("dF_%s",track_label),"",1,0,1); hSty3(dF,kAzure+7,2,2);
0134     TH1F* dS = new TH1F(Form("dS_%s",track_label),"",1,0,1); hSty3(dS,kRed+1,  1,2);
0135     leg->AddEntry(dD,"Data",                    "l");
0136     leg->AddEntry(dF,"SV sim (filtered)",       "l");
0137     leg->AddEntry(dS,"SV sim (filt. + smeared)","l");
0138     leg->Draw();
0139
0140     TLatex tex; tex.SetNDC(); tex.SetTextFont(42);
0141     tex.SetTextSize(0.085);
0142     tex.DrawLatex(0.08, 0.86, "#it{#bf{sPHENIX}} Internal");
0143     tex.SetTextSize(0.075);
0144     tex.DrawLatex(0.08, 0.78, "#it{p}+#it{p}  #sqrt{s} = 200 GeV");
0145     tex.SetTextColor(kGray+2); tex.SetTextSize(0.070);
0146     tex.DrawLatex(0.08, 0.18, getDateStr3().c_str());
0147 }
0148
0149 // ── build canvas ──────────────────────────────────────────────────────────────
0150 void makeCanvas3(TH1F* hD[], TH1F* hF[], TH1F* hS[],
0151                  const char* cname, const char* save_base,
0152                  const char* xtitle, const char* track_label)
0153 {
0154     TCanvas* c = new TCanvas(cname, cname, 1600, 800);
0155     c->Divide(4, 2, 0.003, 0.003);
0156
0157     for (int b = 0; b < NBINS; ++b)
0158         drawDptPad(c->GetPad(b+1), hD[b], hF[b], hS[b], b, xtitle);
0159
0160     drawLegendPad3(c->GetPad(8), track_label);
0161
0162     c->SaveAs(TString(save_base)+".pdf");
0163     c->SaveAs(TString(save_base)+".png");
0164     std::cout << "Saved " << save_base << ".pdf/.png\n";
0165 }
0166
0167 // ── main ──────────────────────────────────────────────────────────────────────
0168 void compare_daughter_pt_ptbins()
0169 {
0170     gStyle->SetOptStat(0);
0171     gStyle->SetOptTitle(0);
0172
0173     const char* KS_DATA   = "KShort6RunCombined.root";
0174     const char* KS_FILT   = "outputKFParticleKShortRecoSV_filtered.root";
0175     const char* KS_SMEAR  = "outputKFParticleKShortRecoSV_filtered_smeared.root";
0176     const char* LAM_DATA  = "Lambda6RunCombined.root";
0177     const char* LAM_FILT  = "outputKFParticleLambda0SV_filtered.root";
0178     const char* LAM_SMEAR = "outputKFParticleLambda0SV_filtered_smeared.root";
0179
0180     // ── K0S ──────────────────────────────────────────────────────────────────
0181     std::cout << "\n=== K0S daughter pT in K0S pT bins ===\n";
0182
0183     TH1F *ksD1[NBINS], *ksD2[NBINS], *ksF1[NBINS], *ksF2[NBINS], *ksS1[NBINS], *ksS2[NBINS];
0184     for (int b = 0; b < NBINS; ++b) {
0185         ksD1[b] = new TH1F(Form("ksD1_%d",b),"", DPT_NB, DPT_LO, DPT_HI);
0186         ksD2[b] = new TH1F(Form("ksD2_%d",b),"", DPT_NB, DPT_LO, DPT_HI);
0187         ksF1[b] = new TH1F(Form("ksF1_%d",b),"", DPT_NB, DPT_LO, DPT_HI);
0188         ksF2[b] = new TH1F(Form("ksF2_%d",b),"", DPT_NB, DPT_LO, DPT_HI);
0189         ksS1[b] = new TH1F(Form("ksS1_%d",b),"", DPT_NB, DPT_LO, DPT_HI);
0190         ksS2[b] = new TH1F(Form("ksS2_%d",b),"", DPT_NB, DPT_LO, DPT_HI);
0191     }
0192
0193     fillDptFiltered(ksD1, ksD2, KS_DATA,  "K_S0_pT");
0194     fillDptFiltered(ksF1, ksF2, KS_FILT,  "K_S0_pT");
0195     fillDptFiltered(ksS1, ksS2, KS_SMEAR, "K_S0_pT");
0196
0197     makeCanvas3(ksD1, ksF1, ksS1, "cKS1",
0198                 "compare_ks_t1pt_ptbins",
0199                 "#pi^{+} #it{p}_{T} (GeV/#it{c})", "ks1");
0200
0201     makeCanvas3(ksD2, ksF2, ksS2, "cKS2",
0202                 "compare_ks_t2pt_ptbins",
0203                 "#pi^{-} #it{p}_{T} (GeV/#it{c})", "ks2");
0204
0205     // ── Lambda ───────────────────────────────────────────────────────────────
0206     std::cout << "\n=== Lambda daughter pT in Lambda pT bins ===\n";
0207
0208     TH1F *lmD1[NBINS], *lmD2[NBINS], *lmF1[NBINS], *lmF2[NBINS], *lmS1[NBINS], *lmS2[NBINS];
0209     for (int b = 0; b < NBINS; ++b) {
0210         lmD1[b] = new TH1F(Form("lmD1_%d",b),"", DPT_NB, DPT_LO, DPT_HI);
0211         lmD2[b] = new TH1F(Form("lmD2_%d",b),"", DPT_NB, PRO_LO, PRO_HI);
0212         lmF1[b] = new TH1F(Form("lmF1_%d",b),"", DPT_NB, DPT_LO, DPT_HI);
0213         lmF2[b] = new TH1F(Form("lmF2_%d",b),"", DPT_NB, PRO_LO, PRO_HI);
0214         lmS1[b] = new TH1F(Form("lmS1_%d",b),"", DPT_NB, DPT_LO, DPT_HI);
0215         lmS2[b] = new TH1F(Form("lmS2_%d",b),"", DPT_NB, PRO_LO, PRO_HI);
0216     }
0217
0218     fillDptFiltered(lmD1, lmD2, LAM_DATA, "Lambda0_pT");
0219     fillDptFiltered(lmF1, lmF2, LAM_FILT,  "Lambda0_pT");
0220     fillDptFiltered(lmS1, lmS2, LAM_SMEAR, "Lambda0_pT");
0221
0222     makeCanvas3(lmD1, lmF1, lmS1, "cLam1",
0223                 "compare_lam_t1pt_ptbins",
0224                 "#pi^{-} #it{p}_{T} (GeV/#it{c})", "lam1");
0225
0226     makeCanvas3(lmD2, lmF2, lmS2, "cLam2",
0227                 "compare_lam_t2pt_ptbins",
0228                 "p #it{p}_{T} (GeV/#it{c})", "lam2");
0229
0230     std::cout << "\nAll done.\n";
0231 }