File indexing completed on 2026-07-16 08:11:23
0001
0002 """
0003 ratio_of_ratios.py
0004 ------------------
0005 Reads the smearing systematic tables for K0S and Lambda
0006 (filtered_ks_syst.txt, filtered_lam_syst.txt) and computes the
0007 ratio of their yield ratios per pT bin:
0008
0009 R = (N_Lam_smear / N_Lam_raw) / (N_KS_smear / N_KS_raw)
0010
0011 R = 1 means smearing affects both species identically → full cancellation.
0012 R != 1 shows the residual systematic on the Lambda/K0S ratio.
0013
0014 Stat error on R propagated as:
0015 σ_R / R = sqrt( (σ_ratio_lam / ratio_lam)^2 + (σ_ratio_ks / ratio_ks)^2 )
0016
0017 Outputs:
0018 ratio_of_ratios.txt
0019 plot_ratio_of_ratios.C
0020 """
0021
0022 import numpy as np
0023
0024 KS_TABLE = "filtered_ks_syst.txt"
0025 LAM_TABLE = "filtered_lam_syst.txt"
0026
0027
0028 def read_table(path):
0029 rows = []
0030 with open(path) as fh:
0031 for line in fh:
0032 line = line.strip()
0033 if not line or line.startswith("#"):
0034 continue
0035 if line.startswith("bin_lo"):
0036 cols = line.split(",")
0037 continue
0038 parts = line.split(",")
0039 d = {c: parts[i] for i, c in enumerate(cols)}
0040 if d["ratio"] == "nan":
0041 continue
0042 rows.append(dict(
0043 lo = float(d["bin_lo"]),
0044 hi = float(d["bin_hi"]),
0045 avg_pt = float(d["avg_pt"]),
0046 nr = int(d["N_raw"]),
0047 ns = int(d["N_smeared"]),
0048 ratio = float(d["ratio"]),
0049 ratio_stat = float(d["ratio_stat_err"]),
0050 syst = float(d["syst_pct"]),
0051 ))
0052 return rows
0053
0054
0055 def main():
0056 ks_rows = read_table(KS_TABLE)
0057 lam_rows = read_table(LAM_TABLE)
0058
0059
0060 ks_map = {(r["lo"], r["hi"]): r for r in ks_rows}
0061 lam_map = {(r["lo"], r["hi"]): r for r in lam_rows}
0062
0063 common_bins = sorted(set(ks_map) & set(lam_map))
0064 if not common_bins:
0065 print("ERROR: no matching pT bins between the two tables.")
0066 return
0067
0068 results = []
0069 hdr = f"{'pT bin':>16} {'avg_pT':>7} {'R_lam':>7} {'R_ks':>7} {'R_lam/R_ks':>11} {'stat_err':>9} {'dev%':>7}"
0070 print(hdr)
0071 print("─" * len(hdr))
0072
0073 for (lo, hi) in common_bins:
0074 ks = ks_map[(lo, hi)]
0075 lam = lam_map[(lo, hi)]
0076
0077 ror = lam["ratio"] / ks["ratio"]
0078
0079
0080 rel_err = np.sqrt((lam["ratio_stat"] / lam["ratio"])**2 +
0081 (ks["ratio_stat"] / ks["ratio"])**2)
0082 ror_err = ror * rel_err
0083
0084 dev_pct = abs(ror - 1.0) * 100.0
0085
0086 avg_pt = (ks["avg_pt"] + lam["avg_pt"]) / 2.0
0087
0088 results.append(dict(
0089 lo=lo, hi=hi,
0090 avg_pt=avg_pt,
0091 ks_avg_pt=ks["avg_pt"],
0092 lam_avg_pt=lam["avg_pt"],
0093 ratio_ks=ks["ratio"],
0094 ratio_lam=lam["ratio"],
0095 ror=ror,
0096 ror_err=ror_err,
0097 dev_pct=dev_pct,
0098 ))
0099 print(f" [{lo:.1f},{hi:.1f}) {avg_pt:>7.3f} "
0100 f"{lam['ratio']:>7.4f} {ks['ratio']:>7.4f} "
0101 f"{ror:>11.4f} {ror_err:>9.4f} {dev_pct:>7.2f}%")
0102
0103
0104 out_txt = "ratio_of_ratios.txt"
0105 with open(out_txt, "w") as fh:
0106 fh.write("# ratio_of_ratios.txt\n")
0107 fh.write("# R = (N_Lam_smear/N_Lam_raw) / (N_KS_smear/N_KS_raw)\n")
0108 fh.write("# R=1 → smearing cancels identically; |R-1| = residual syst on Lam/KS ratio\n\n")
0109 fh.write("bin_lo,bin_hi,avg_pt,ratio_ks,ratio_lam,ror,ror_stat_err,dev_pct\n")
0110 for r in results:
0111 fh.write(f"{r['lo']:.2f},{r['hi']:.2f},{r['avg_pt']:.4f},"
0112 f"{r['ratio_ks']:.6f},{r['ratio_lam']:.6f},"
0113 f"{r['ror']:.6f},{r['ror_err']:.6f},{r['dev_pct']:.4f}\n")
0114
0115 valid = results
0116 dev_vals = [r["dev_pct"] for r in valid]
0117 print(f"\n R range: {min(r['ror'] for r in valid):.4f} – {max(r['ror'] for r in valid):.4f}")
0118 print(f" |R-1| range: {min(dev_vals):.2f}% – {max(dev_vals):.2f}%")
0119 print(f" Mean |R-1|: {np.mean(dev_vals):.2f}%")
0120 print(f"\nTable → {out_txt}")
0121
0122
0123 N = len(results)
0124 pts = ", ".join(f"{r['avg_pt']:.4f}" for r in results)
0125 rors = ", ".join(f"{r['ror']:.6f}" for r in results)
0126 errs = ", ".join(f"{r['ror_err']:.6f}" for r in results)
0127 devs = ", ".join(f"{r['dev_pct']:.4f}" for r in results)
0128 xerrs = ", ".join(f"{(r['hi']-r['lo'])/2:.3f}" for r in results)
0129 r_ks = ", ".join(f"{r['ratio_ks']:.6f}" for r in results)
0130 r_lam = ", ".join(f"{r['ratio_lam']:.6f}"for r in results)
0131 pt_lo = results[0]["lo"]
0132 pt_hi = results[-1]["hi"]
0133
0134 macro = f"""// plot_ratio_of_ratios.C
0135 // R = (N_Lam_smear/N_Lam_raw) / (N_KS_smear/N_KS_raw) per pT bin.
0136 // R=1: smearing cancels in the Lambda/K0S ratio.
0137 // |R-1|: residual systematic on the Lambda/K0S yield ratio from smearing.
0138 // Run with: root -l -b -q plot_ratio_of_ratios.C
0139
0140 #include <ctime>
0141 #include <sstream>
0142 #include <algorithm>
0143 #include <cmath>
0144
0145 std::string getDate2(){{
0146 std::time_t t=std::time(0); std::tm* now=std::localtime(&t);
0147 std::stringstream d;
0148 d<<(now->tm_mon+1)<<'/'<<now->tm_mday<<'/'<<(now->tm_year+1900);
0149 return d.str();
0150 }}
0151 auto addLabel=[](double x1,double y1,double x2,double y2){{
0152 TPaveText *pt=new TPaveText(x1,y1,x2,y2,"NDC");
0153 pt->SetFillStyle(0); pt->SetBorderSize(0); pt->SetTextFont(42);
0154 pt->AddText("#it{{#bf{{sPHENIX}}}} Internal, #it{{p}}+#it{{p}} #sqrt{{s}} = 200 GeV");
0155 pt->Draw();
0156 }};
0157 auto addDate=[](double x1,double y1,double x2,double y2){{
0158 TPaveText *pt=new TPaveText(x1,y1,x2,y2,"NDC");
0159 pt->SetFillStyle(0); pt->SetBorderSize(0); pt->SetTextFont(42);
0160 pt->AddText(getDate2().c_str()); pt->Draw();
0161 }};
0162
0163 void plot_ratio_of_ratios() {{
0164 gStyle->SetOptStat(0); gStyle->SetOptTitle(0);
0165
0166 const int N = {N};
0167 double avg_pt[] = {{{pts}}};
0168 double ror[] = {{{rors}}};
0169 double ror_err[] = {{{errs}}};
0170 double xerr[] = {{{xerrs}}};
0171 double dev_pct[] = {{{devs}}};
0172 double ratio_ks[] = {{{r_ks}}};
0173 double ratio_lam[] = {{{r_lam}}};
0174
0175 // ── canvas 1: ratio of ratios ─────────────────────────────────────────────
0176 TCanvas *c1 = new TCanvas("c1","Lambda/KS ratio of ratios",900,650);
0177 c1->SetLeftMargin(0.13); c1->SetBottomMargin(0.13); c1->SetRightMargin(0.06);
0178
0179 TGraphErrors *gRoR = new TGraphErrors(N,avg_pt,ror,xerr,ror_err);
0180 gRoR->SetMarkerStyle(20); gRoR->SetMarkerSize(1.4);
0181 gRoR->SetMarkerColor(kViolet+1); gRoR->SetLineColor(kViolet+1); gRoR->SetLineWidth(2);
0182 gRoR->Draw("APE");
0183 gRoR->GetXaxis()->SetTitle("#it{{p}}_{{T}} (GeV/#it{{c}})");
0184 gRoR->GetYaxis()->SetTitle("(#Lambda smear ratio) / (K_{{S}}^{{0}} smear ratio)");
0185 gRoR->GetYaxis()->SetRangeUser(0.7, 1.3);
0186 gRoR->GetXaxis()->SetTitleSize(0.05);
0187 gRoR->GetYaxis()->SetTitleSize(0.045);
0188 gRoR->GetYaxis()->SetTitleOffset(1.3);
0189
0190 TLine *unity = new TLine({pt_lo},1.0,{pt_hi},1.0);
0191 unity->SetLineStyle(2); unity->SetLineColor(kGray+1); unity->SetLineWidth(2);
0192 unity->Draw();
0193
0194 TBox *band = new TBox({pt_lo},0.95,{pt_hi},1.05);
0195 band->SetFillColorAlpha(kGray,0.30); band->SetLineWidth(0); band->Draw();
0196 gRoR->Draw("PE SAME");
0197
0198 addLabel(0.14,0.88,0.65,0.96);
0199 addDate(0.70,0.88,0.98,0.96);
0200
0201 TLegend *leg1 = new TLegend(0.14,0.72,0.80,0.87);
0202 leg1->SetBorderSize(0); leg1->SetFillStyle(0); leg1->SetTextSize(0.030);
0203 leg1->AddEntry(gRoR, "(N_{{#Lambda}}^{{smear}}/N_{{#Lambda}}^{{raw}}) / (N_{{K_{{S}}^{{0}}}}^{{smear}}/N_{{K_{{S}}^{{0}}}}^{{raw}}) (stat. err)", "lpe");
0204 leg1->AddEntry(unity, "Unity (full cancellation)", "l");
0205 leg1->AddEntry(band, "#pm5% band", "f");
0206 leg1->Draw();
0207
0208 c1->SaveAs("plot_ratio_of_ratios.pdf");
0209 c1->SaveAs("plot_ratio_of_ratios.png");
0210 std::cout << "Saved plot_ratio_of_ratios.pdf/.png\\n";
0211
0212 // ── canvas 2: individual ratios overlaid ──────────────────────────────────
0213 TCanvas *c2 = new TCanvas("c2","Individual smear ratios",900,650);
0214 c2->SetLeftMargin(0.13); c2->SetBottomMargin(0.13); c2->SetRightMargin(0.06);
0215
0216 TGraph *gKS = new TGraph(N, avg_pt, ratio_ks);
0217 TGraph *gLam = new TGraph(N, avg_pt, ratio_lam);
0218 gKS->SetMarkerStyle(20); gKS->SetMarkerSize(1.4);
0219 gKS->SetMarkerColor(kBlue+1); gKS->SetLineColor(kBlue+1); gKS->SetLineWidth(2);
0220 gLam->SetMarkerStyle(21); gLam->SetMarkerSize(1.4);
0221 gLam->SetMarkerColor(kRed+1); gLam->SetLineColor(kRed+1); gLam->SetLineWidth(2);
0222
0223 TMultiGraph *mg = new TMultiGraph();
0224 mg->Add(gKS,"PL"); mg->Add(gLam,"PL"); mg->Draw("A");
0225 mg->GetXaxis()->SetTitle("#it{{p}}_{{T}} (GeV/#it{{c}})");
0226 mg->GetYaxis()->SetTitle("N_{{smeared}} / N_{{raw}}");
0227 mg->GetYaxis()->SetRangeUser(0.7, 1.3);
0228 mg->GetXaxis()->SetTitleSize(0.05); mg->GetYaxis()->SetTitleSize(0.05);
0229
0230 TLine *unity2 = new TLine({pt_lo},1.0,{pt_hi},1.0);
0231 unity2->SetLineStyle(2); unity2->SetLineColor(kGray+1); unity2->SetLineWidth(2);
0232 unity2->Draw();
0233 TBox *band2 = new TBox({pt_lo},0.95,{pt_hi},1.05);
0234 band2->SetFillColorAlpha(kGray,0.30); band2->SetLineWidth(0); band2->Draw();
0235 gKS->Draw("PL SAME"); gLam->Draw("PL SAME");
0236
0237 TLegend *leg2 = new TLegend(0.14,0.74,0.60,0.87);
0238 leg2->SetBorderSize(0); leg2->SetFillStyle(0); leg2->SetTextSize(0.032);
0239 leg2->AddEntry(gKS, "K_{{S}}^{{0}} smear ratio","lp");
0240 leg2->AddEntry(gLam, "#Lambda^{{0}} smear ratio","lp");
0241 leg2->Draw();
0242 addLabel(0.14,0.88,0.65,0.96);
0243 addDate(0.70,0.88,0.98,0.96);
0244
0245 c2->SaveAs("plot_smear_ratios_overlay.pdf");
0246 c2->SaveAs("plot_smear_ratios_overlay.png");
0247 std::cout << "Saved plot_smear_ratios_overlay.pdf/.png\\n";
0248
0249 // ── canvas 3: residual systematic on Lambda/K0S ratio ────────────────────
0250 TCanvas *c3 = new TCanvas("c3","Residual syst on Lam/KS ratio",900,650);
0251 c3->SetLeftMargin(0.13); c3->SetBottomMargin(0.13); c3->SetRightMargin(0.06);
0252
0253 TGraph *gDev = new TGraph(N, avg_pt, dev_pct);
0254 gDev->SetMarkerStyle(20); gDev->SetMarkerSize(1.4);
0255 gDev->SetMarkerColor(kOrange+1); gDev->SetLineColor(kOrange+1); gDev->SetLineWidth(2);
0256 gDev->Draw("APL");
0257 gDev->GetXaxis()->SetTitle("#it{{p}}_{{T}} (GeV/#it{{c}})");
0258 gDev->GetYaxis()->SetTitle("Residual systematic on #Lambda/K_{{S}}^{{0}} ratio |R-1| (%)");
0259 double ymax = *std::max_element(dev_pct,dev_pct+N)*1.4 + 0.3;
0260 gDev->GetYaxis()->SetRangeUser(0, ymax);
0261 gDev->GetXaxis()->SetTitleSize(0.05); gDev->GetYaxis()->SetTitleSize(0.042);
0262 gDev->GetYaxis()->SetTitleOffset(1.4);
0263 addLabel(0.14,0.88,0.65,0.96);
0264 addDate(0.70,0.88,0.98,0.96);
0265
0266 c3->SaveAs("plot_ratio_of_ratios_dev.pdf");
0267 c3->SaveAs("plot_ratio_of_ratios_dev.png");
0268 std::cout << "Saved plot_ratio_of_ratios_dev.pdf/.png\\n";
0269 }}
0270 """
0271 out_mac = "plot_ratio_of_ratios.C"
0272 with open(out_mac, "w") as fh:
0273 fh.write(macro)
0274 print(f"Macro → {out_mac}")
0275 print(f"\nRun: root -l -b -q {out_mac}")
0276
0277
0278 if __name__ == "__main__":
0279 main()
