File indexing completed on 2026-07-16 08:11:24
0001 float lower_bin_bounds[] = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.8, 2.1, 2.4, 2.7, 3, 4};
0002 const unsigned int n_variable_bins = sizeof(lower_bin_bounds)/sizeof(lower_bin_bounds[0]) - 1;
0003
0004 string plotPath = "plots_looser_cuts/";
0005
0006 const unsigned int Kshort_polynomial_order = 5;
0007 const unsigned int Lambda_polynomial_order = 8;
0008
0009 struct fit_ranges
0010 {
0011 float min;
0012 float ignore_min;
0013 float ignore_max;
0014 float max;
0015 };
0016
0017 fit_ranges Kshort_ranges{0.35, 0.47, 0.535, 0.65};
0018 fit_ranges Lambda_ranges{1.09, 1.105, 1.125, 1.15};
0019
0020 float Kshort_xVals[n_variable_bins], Kshort_xErrs[n_variable_bins], Kshort_yVals[n_variable_bins], Kshort_yErrs[n_variable_bins];
0021 float Lambda_xVals[n_variable_bins], Lambda_xErrs[n_variable_bins], Lambda_yVals[n_variable_bins], Lambda_yErrs[n_variable_bins];
0022
0023 TF1 *fitFunc, *fitFunc_fullBkg;
0024
0025 template <typename T>
0026 string to_string_with_precision(const T a_value, const int n = 0)
0027 {
0028 ostringstream out;
0029 out.precision(n);
0030 out << fixed << a_value;
0031 return out.str();
0032 }
0033
0034 class fitVals
0035 {
0036 public:
0037 string pTmin;
0038 string pTmax;
0039 string yield;
0040 string error;
0041
0042 fitVals(float min, float max, float fitYield, float fitErr)
0043 {
0044 pTmin = to_string_with_precision(min, 1);
0045 pTmax = to_string_with_precision(max, 1);
0046 yield = to_string_with_precision(fitYield, 0);
0047 error = to_string_with_precision(fitErr, 0);
0048 };
0049 ~fitVals();
0050
0051 };
0052
0053
0054 TH1F makeHisto(int nBins, float min, float max, string type, string xAxisTitle, string unit, int precision)
0055 {
0056 string histo_name = type;
0057
0058 TH1F myHisto(histo_name.c_str(), histo_name.c_str(), nBins, min, max);
0059
0060 if (unit != "") xAxisTitle += " [" + unit + "]";
0061 myHisto.GetXaxis()->SetTitle(xAxisTitle.c_str());
0062
0063 float binWidth = (float) (max - min)/nBins;
0064 string yAxisTitle = "Candidates";
0065 if (unit != "") yAxisTitle += " / " + to_string_with_precision(binWidth, precision) + " " + unit;
0066 myHisto.GetYaxis()->SetTitle(yAxisTitle.c_str());
0067
0068 return myHisto;
0069 }
0070
0071 Double_t fitfunc_Kshort(Double_t *x, Double_t *par)
0072 {
0073 if (x[0] > Kshort_ranges.ignore_min && x[0] < Kshort_ranges.ignore_max)
0074 {
0075 TF1::RejectPoint();
0076 return 0;
0077 }
0078
0079 double value = 0;
0080 for (unsigned int i = 0; i <= Kshort_polynomial_order; ++i) value += par[i]*pow(x[0], i);
0081
0082 return value;
0083 }
0084
0085 Double_t fitfunc_Lambda(Double_t *x, Double_t *par)
0086 {
0087 if (x[0] > Lambda_ranges.ignore_min && x[0] < Lambda_ranges.ignore_max)
0088 {
0089 TF1::RejectPoint();
0090 return 0;
0091 }
0092
0093 double value = 0;
0094 for (unsigned int i = 0; i <= Lambda_polynomial_order; ++i) value += par[i]*pow(x[0], i);
0095
0096 return value;
0097 }
0098
0099 Double_t fitfunc_noGap(Double_t *x, Double_t *par)
0100 {
0101 double value = 0;
0102 for (unsigned int i = 0; i <= Kshort_polynomial_order; ++i) value += par[i]*pow(x[0], i);
0103
0104 return value;
0105 }
0106
0107 Double_t fitfunc_Lambda_noGap(Double_t *x, Double_t *par)
0108 {
0109 double value = 0;
0110 for (unsigned int i = 0; i <= Lambda_polynomial_order; ++i) value += par[i]*pow(x[0], i);
0111
0112 return value;
0113 }
0114
0115 template <typename T>
0116 void savePlots(T myPlot, string plotName, float xMin, float xMax, float yield = 0, float error = 0, float bkg = 0, float bkgErr = 0, float total = 0)
0117 {
0118 TGaxis::SetMaxDigits(3);
0119 string makeDirectory = "mkdir -p " + plotPath;
0120 system(makeDirectory.c_str());
0121
0122 TCanvas *c1 = new TCanvas("myCanvas", "myCanvas",800,800);
0123
0124 myPlot.GetXaxis()->SetNdivisions(505);
0125 myPlot.GetYaxis()->SetRangeUser(0, myPlot.GetMaximum()*1.6);
0126 myPlot.Sumw2();
0127 myPlot.Draw("PE1");
0128 fitFunc->Draw("SAME");
0129 fitFunc_fullBkg->Draw("SAME");
0130
0131 TPaveText *pt;
0132 pt = new TPaveText(0.15,0.74,0.95,0.89, "NDC");
0133 pt->SetFillColor(0);
0134 pt->SetFillStyle(0);
0135 pt->SetTextFont(42);
0136 TText *pt_LaTex;
0137 pt->AddText("#it{#bf{sPHENIX}} Internal, #sqrt{s} = 200 GeV, pp");
0138 string range = to_string_with_precision(xMin, 1) + " #leq p_{T} [GeV] < " + to_string_with_precision(xMax, 1) + ", N. Cand. = " + to_string_with_precision(total, 0);
0139 string result = "Signal yield = " + to_string_with_precision(yield, 0) + " #pm " + to_string_with_precision(error, 0);
0140 string result_bkg = "Background yield = " + to_string_with_precision(bkg, 0) + " #pm " + to_string_with_precision(bkgErr, 0);
0141 pt->AddText(range.c_str());
0142 pt->AddText(result.c_str());
0143 pt->AddText(result_bkg.c_str());
0144 pt->SetBorderSize(0);
0145 pt->Draw();
0146 gPad->Modified();
0147
0148 string extensions[] = {".C", ".pdf", ".png", ".root"};
0149 for (auto extension : extensions)
0150 {
0151 string output = plotPath + plotName + extension;
0152 c1->SaveAs(output.c_str());
0153 }
0154 }
0155
0156 void processData(string type = "Kshort2pipi")
0157 {
0158 bool processingKshort = type == "Kshort2pipi" ? true : false;
0159
0160 float xVals[n_variable_bins], xErrs[n_variable_bins], yVals[n_variable_bins], yErrs[n_variable_bins];
0161
0162 string dir = "/sphenix/tg/tg01/hf/cdean/LF_analysis/data_nTuples/";
0163 string fileName = processingKshort ? dir + "output_Kshort_run3pp_looseCuts_20260608.root"
0164 : dir + "output_Lambda0_run3pp_looseCuts_20260608.root";
0165
0166 TFile *file = new TFile(fileName.c_str());
0167 TTree* data = (TTree*)file->Get("DecayTree");
0168 string massBranch = processingKshort ? "K_S0_mass" : "Lambda0_mass";
0169 string pTBranch = processingKshort ? "K_S0_pT" : "Lambda0_pT";
0170 string yBranch = processingKshort ? "K_S0_rapidity" : "Lambda0_rapidity";
0171
0172 string mass_string = processingKshort ? "m_{#pi#pi}" : "m_{p#pi}";
0173
0174 float mass_min = processingKshort ? Kshort_ranges.min : Lambda_ranges.min;
0175 float mass_max = processingKshort ? Kshort_ranges.max : Lambda_ranges.max;
0176
0177 for (int i = 0; i < n_variable_bins; ++i)
0178 {
0179 float min = lower_bin_bounds[i];
0180 float max = lower_bin_bounds[i+1];
0181 float signal_yield = 0;
0182 float signal_error = 0;
0183
0184 string title = type + "_pT_range_" + to_string_with_precision(min,1) + "_to_" + to_string_with_precision(max,1);
0185 TH1F binnedMass = makeHisto(50, mass_min, mass_max, title.c_str(), mass_string.c_str(), "GeV", 3);
0186 string fillString = massBranch + " >> " + title;
0187 string cutString = to_string_with_precision(min,1) + " <= " + pTBranch + " && " + pTBranch + " < " + to_string_with_precision(max,1)
0188 + " && " + to_string(mass_min) + " <= " + massBranch + " && " + massBranch + " <= " + to_string(mass_max)
0189 + " && min(track_1_pT, track_2_pT) >= 0.2 && abs(" + yBranch + ") <= 0.5";
0190 data->Draw(fillString.c_str(), cutString.c_str());
0191
0192 if (processingKshort) fitFunc = new TF1("fit", fitfunc_Kshort, mass_min, mass_max, Kshort_polynomial_order);
0193 else fitFunc = new TF1("fit", fitfunc_Lambda, mass_min, mass_max, Lambda_polynomial_order);
0194 fitFunc->SetLineColor(kRed);
0195 TFitResultPtr r = binnedMass.Fit(fitFunc, "RS");
0196
0197
0198 if (processingKshort) fitFunc_fullBkg = new TF1("fit", fitfunc_noGap, mass_min, mass_max, Kshort_polynomial_order);
0199 else fitFunc_fullBkg = new TF1("fit", fitfunc_Lambda_noGap, mass_min, mass_max, Lambda_polynomial_order);
0200 fitFunc_fullBkg->SetLineColor(kBlue);
0201 for (int j = 0; j < fitFunc_fullBkg->GetNpar(); ++j) fitFunc_fullBkg->SetParameter(j, fitFunc->GetParameter(j));
0202
0203 float bkg_area = fitFunc_fullBkg->Integral(mass_min, mass_max);
0204 float bkg_areaErr = fitFunc->IntegralError(mass_min, mass_max, r->GetParams(), r->GetCovarianceMatrix().GetMatrixArray());
0205
0206 float binWidth = binnedMass.GetBinWidth(1);
0207 float bkg_yield = bkg_area / binWidth;
0208 float bkg_yieldErr = bkg_areaErr / binWidth;
0209
0210 signal_yield = (float) binnedMass.GetEntries() - bkg_yield;
0211 signal_error = signal_yield*(bkg_yieldErr/bkg_yield);
0212
0213 float nEntries = (float) binnedMass.GetEntries();
0214 cout << "Number of entries in histogram = " << nEntries << endl;
0215 cout << "Background yield from total integral = " << bkg_yield << " +/- " << bkg_yieldErr << endl;
0216 cout << "Signal yield = " << signal_yield << " +/- " << signal_error << endl;
0217
0218 savePlots(binnedMass, title.c_str(), min, max, signal_yield, signal_error, bkg_yield, bkg_yieldErr, nEntries);
0219
0220 if (processingKshort)
0221 {
0222 Kshort_xVals[i] = xVals[i] = (max + min)/2;
0223 Kshort_xErrs[i] = xErrs[i] = (max - min)/2;
0224 Kshort_yVals[i] = yVals[i] = signal_yield;
0225 Kshort_yErrs[i] = yErrs[i] = signal_error;
0226 }
0227 else
0228 {
0229 Lambda_xVals[i] = xVals[i] = (max + min)/2;
0230 Lambda_xErrs[i] = xErrs[i] = (max - min)/2;
0231 Lambda_yVals[i] = yVals[i] = signal_yield;
0232 Lambda_yErrs[i] = yErrs[i] = signal_error;
0233 }
0234 }
0235
0236 TCanvas *c1 = new TCanvas("c1","A Simple Graph",800,800);
0237
0238 TGraphErrors *gr = new TGraphErrors(n_variable_bins,xVals,yVals,xErrs,yErrs);
0239 string graph_x_axis_title = processingKshort ? "K^{0}_{S} p_{T} [GeV]" : "#Lambda^{0} p_{T} [GeV]";
0240 gr->GetXaxis()->SetTitle(graph_x_axis_title.c_str());
0241 gr->GetYaxis()->SetTitle("Yield");
0242 gr->SetMarkerSize(2);
0243 gr->Draw("A*");
0244
0245 string extensions[] = {".C", ".pdf", ".png", ".root"};
0246 for (auto extension : extensions)
0247 {
0248 string output = plotPath + type + "_yield" + extension;
0249 c1->SaveAs(output.c_str());
0250 }
0251 }
0252
0253 void extractYields()
0254 {
0255 processData();
0256 processData("Lambda02ppi");
0257
0258 TCanvas *c1 = new TCanvas("c1","A Simple Graph",800,800);
0259
0260 float xVals[n_variable_bins], xErrs[n_variable_bins], yVals[n_variable_bins], yErrs[n_variable_bins];
0261
0262 for (int i = 0; i < n_variable_bins; ++i)
0263 {
0264 xVals[i] = Kshort_xVals[i];
0265 xErrs[i] = Kshort_xErrs[i];
0266
0267 float rawRatio = Lambda_yVals[i]/(2*Kshort_yVals[i]);
0268 float rawError = rawRatio*(TMath::Sqrt(TMath::Power(Lambda_yErrs[i]/Lambda_yVals[i], 2) + TMath::Power(Kshort_yErrs[i]/Kshort_yVals[i], 2)));
0269 yVals[i] = rawRatio;
0270 yErrs[i] = rawError;
0271
0272 cout << "pT bin mean: " << xVals[i]
0273 << ", Lambda yield = " << to_string_with_precision(Lambda_yVals[i], 0) << " +/- " << to_string_with_precision(Lambda_yErrs[i], 0)
0274 << ", Kshort yield = " << to_string_with_precision(Kshort_yVals[i], 0) << " +/- " << to_string_with_precision(Kshort_yErrs[i], 0)
0275 << ", raw ratio = " << yVals[i] << " +/- " << yErrs[i] << endl;
0276 }
0277
0278 TGraphErrors *gr = new TGraphErrors(n_variable_bins,xVals,yVals,xErrs,yErrs);
0279 gr->GetXaxis()->SetTitle("p_{T} [GeV]");
0280 gr->GetYaxis()->SetTitle("#Lambda^{0}/2K^{0}_{S} Raw Ratio");
0281 gr->GetXaxis()->SetRangeUser(0.5, 4);
0282 gr->GetYaxis()->SetRangeUser(0, 0.415);
0283 gr->SetMarkerSize(2);
0284 gr->Draw("A*");
0285
0286 string extensions[] = {".C", ".pdf", ".png", ".root"};
0287 for (auto extension : extensions)
0288 {
0289 string output = plotPath + "RawRatio" + extension;
0290 c1->SaveAs(output.c_str());
0291 }
0292 }
