File indexing completed on 2026-07-16 08:11:23
0001 #include <algorithm>
0002
0003 float lower_bin_bounds_pion[] = {0, 0.1, 0.2, 0.3, 0.4, 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};
0004 const unsigned int n_variable_bins_pion = sizeof(lower_bin_bounds_pion)/sizeof(lower_bin_bounds_pion[0]) - 1;
0005 float lower_bin_bounds[] = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
0006 const unsigned int n_variable_bins = sizeof(lower_bin_bounds)/sizeof(lower_bin_bounds[0]) - 1;
0007
0008 float min_track_pT_proton = 1.0;
0009 float min_track_pT_kaon = 0.75;
0010 int nBins_eta = 50;
0011 int nBins_phi = 20;
0012 int nBins_pT = 20;
0013 float eta_max = 1.2;
0014 float phi_max = M_PI;
0015
0016 string plotPath = "plots_negative_track_to_positive_track_with_DCA/";
0017
0018 string makeDirectory = "mkdir -p " + plotPath;
0019 system(makeDirectory.c_str());
0020
0021 template <typename T>
0022 string to_string_with_precision(const T a_value, const int n = 0)
0023 {
0024 ostringstream out;
0025 out.precision(n);
0026 out << fixed << a_value;
0027 return out.str();
0028 }
0029
0030
0031 TH1F makeHisto(int nBins, float min, float max, string xAxisTitle, string unit, string yAxisTitle)
0032 {
0033 TH1F myHisto(yAxisTitle.c_str(), yAxisTitle.c_str(), nBins, min, max);
0034 if (unit != "") xAxisTitle += " [" + unit + "]";
0035 myHisto.GetXaxis()->SetTitle(xAxisTitle.c_str());
0036 myHisto.GetYaxis()->SetTitle(yAxisTitle.c_str());
0037
0038 return myHisto;
0039 }
0040
0041 TH1F makeHisto(int nBins, float* min, string xAxisTitle, string unit, string yAxisTitle = "Raw Ratio")
0042 {
0043 TH1F myHisto(xAxisTitle.c_str(), xAxisTitle.c_str(), nBins, min);
0044
0045 if (unit != "") xAxisTitle += " [" + unit + "]";
0046 myHisto.GetXaxis()->SetTitle(xAxisTitle.c_str());
0047 myHisto.GetYaxis()->SetTitle(yAxisTitle.c_str());
0048
0049 return myHisto;
0050 }
0051
0052 TH1F pT_proton_histo = makeHisto(nBins_pT, 0, min_track_pT_proton, "p_{T}", "GeV", "p Yield");
0053 TH1F pT_antiproton_histo = makeHisto(nBins_pT, 0, min_track_pT_proton, "p_{T}", "GeV", "#bar{p} Yield");
0054 TH1F pT_proton_ratio = makeHisto(nBins_pT, 0, min_track_pT_proton, "p_{T}", "GeV", "#bar{p}/p Raw Ratio");
0055 TH1F pT_piplus_histo = makeHisto(nBins_pT, 0, 4, "p_{T}", "GeV", "#pi^{+} Yield");
0056 TH1F pT_piminus_histo = makeHisto(nBins_pT, 0, 4, "p_{T}", "GeV", "#pi^{-} Yield");
0057 TH1F pT_pion_ratio = makeHisto(nBins_pT, 0, 4, "p_{T}", "GeV", "#pi^{-}/#pi^{+} Raw Ratio");
0058 TH1F pT_Kplus_histo = makeHisto(nBins_pT, 0, min_track_pT_kaon, "p_{T}", "GeV", "K^{+} Yield");
0059 TH1F pT_Kminus_histo = makeHisto(nBins_pT, 0, min_track_pT_kaon, "p_{T}", "GeV", "K^{-} Yield");
0060 TH1F pT_Kaon_ratio = makeHisto(nBins_pT, 0, min_track_pT_kaon, "p_{T}", "GeV", "K^{-}/K^{+} Raw Ratio");
0061 TH1F eta_proton_histo = makeHisto(nBins_eta, -1*eta_max, eta_max, "#eta", "", "p Yield");
0062 TH1F eta_antiproton_histo = makeHisto(nBins_eta, -1*eta_max, eta_max, "#eta", "", "#bar{p} Yield");
0063 TH1F eta_proton_ratio = makeHisto(nBins_eta, -1*eta_max, eta_max, "#eta", "", "#bar{p}/p Raw Ratio");
0064 TH1F eta_piplus_histo = makeHisto(nBins_eta, -1*eta_max, eta_max, "#eta", "", "#pi^{+} Yield");
0065 TH1F eta_piminus_histo = makeHisto(nBins_eta, -1*eta_max, eta_max, "#eta", "", "#pi^{-} Yield");
0066 TH1F eta_pion_ratio = makeHisto(nBins_eta, -1*eta_max, eta_max, "#eta", "", "#pi^{-}/#pi^{+} Raw Ratio");
0067 TH1F eta_Kplus_histo = makeHisto(nBins_eta, -1*eta_max, eta_max, "#eta", "", "K^{+} Yield");
0068 TH1F eta_Kminus_histo = makeHisto(nBins_eta, -1*eta_max, eta_max, "#eta", "", "K^{-} Yield");
0069 TH1F eta_Kaon_ratio = makeHisto(nBins_eta, -1*eta_max, eta_max, "#eta", "", "K^{-}/K^{+} Raw Ratio");
0070 TH1F phi_proton_histo = makeHisto(nBins_phi, -1*phi_max, phi_max, "#phi", "", "p Yield");
0071 TH1F phi_antiproton_histo = makeHisto(nBins_phi, -1*phi_max, phi_max, "#phi", "", "#bar{p} Yield");
0072 TH1F phi_proton_ratio = makeHisto(nBins_phi, -1*phi_max, phi_max, "#phi", "", "#bar{p}/p Raw Ratio");
0073 TH1F phi_piplus_histo = makeHisto(nBins_phi, -1*phi_max, phi_max, "#phi", "", "#pi^{+} Yield");
0074 TH1F phi_piminus_histo = makeHisto(nBins_phi, -1*phi_max, phi_max, "#phi", "", "#pi^{-} Yield");
0075 TH1F phi_pion_ratio = makeHisto(nBins_phi, -1*phi_max, phi_max, "#phi", "", "#pi^{-}/#pi^{+} Raw Ratio");
0076 TH1F phi_Kplus_histo = makeHisto(nBins_phi, -1*phi_max, phi_max, "#phi", "", "K^{+} Yield");
0077 TH1F phi_Kminus_histo = makeHisto(nBins_phi, -1*phi_max, phi_max, "#phi", "", "K^{-} Yield");
0078 TH1F phi_Kaon_ratio = makeHisto(nBins_phi, -1*phi_max, phi_max, "#phi", "", "K^{-}/K^{+} Raw Ratio");
0079
0080 template <typename T>
0081 void savePlots(T myPlot, string plotName)
0082 {
0083 TGaxis::SetMaxDigits(3);
0084
0085 TCanvas *c1 = new TCanvas("myCanvas", "myCanvas",800,800);
0086
0087 myPlot.GetXaxis()->SetNdivisions(505);
0088 myPlot.GetYaxis()->SetRangeUser(0, myPlot.GetMaximum()*1.2);
0089
0090 myPlot.Draw("PE1");
0091
0092 TPaveText *pt;
0093 pt = new TPaveText(0.15,0.84,0.95,0.89, "NDC");
0094 pt->SetFillColor(0);
0095 pt->SetFillStyle(0);
0096 pt->SetTextFont(42);
0097 TText *pt_LaTex;
0098 pt->AddText("#it{#bf{sPHENIX}} Internal, #sqrt{s} = 200 GeV, pp");
0099 pt->SetBorderSize(0);
0100 pt->Draw();
0101 gPad->Modified();
0102
0103 string extensions[] = {".C", ".pdf", ".png"};
0104 for (auto extension : extensions)
0105 {
0106 string output = plotPath + plotName + extension;
0107 c1->SaveAs(output.c_str());
0108 }
0109 }
0110
0111 template <typename T>
0112 void saveManyPlots(vector<T> myPlot, string plotName, bool logY = false, vector<string> legendText = {""})
0113 {
0114 TCanvas *c1 = new TCanvas("myCanvas", "myCanvas",800,800);
0115
0116 if (logY) gPad->SetLogy();
0117
0118 int colors[] = {4, 2, 6, 1, 3, 6, 7, 8, 9};
0119 Size_t markerStyle[] = {8, 21, 33, 34};
0120 Size_t markerSize = 1;
0121
0122 float legendYpos[2] = {0.65, 0.85};
0123 float legendTextSize = ((legendYpos[1] - legendYpos[0])/(1.6*myPlot.size()));
0124 TLegend *legend = new TLegend(0.65,legendYpos[0],0.80, legendYpos[1]);
0125 legend->SetTextFont(42);
0126
0127 for (unsigned int i = 0; i < myPlot.size(); ++i)
0128 {
0129 float transparency = 1;
0130 myPlot[i].SetFillColorAlpha(colors[i], transparency);
0131 myPlot[i].SetLineColor(colors[i]);
0132 myPlot[i].SetMarkerColorAlpha(colors[i], transparency);
0133 myPlot[i].SetMarkerSize(markerSize);
0134 myPlot[i].SetMarkerStyle(markerStyle[i]);
0135 if (legendText.size() < 2) legend->AddEntry(&myPlot[i], myPlot[i].GetXaxis()->GetTitle(), "PL");
0136 else legend->AddEntry(&myPlot[i], legendText[i].c_str(), "PL");
0137 }
0138
0139 myPlot[0].GetXaxis()->SetNdivisions(505);
0140 float min = std::max(myPlot[0].GetMinimum(), 1.);
0141 myPlot[0].SetMinimum(min*0.8);
0142 float max = logY ? myPlot[0].GetMaximum()*1e3 : myPlot[0].GetMaximum()*1.2;
0143 myPlot[0].SetMaximum(max);
0144 if (strncmp(typeid(myPlot[0]).name(), "4TH2F", 5) == 0)
0145 {
0146 myPlot[0].Draw();
0147 for (unsigned int i = 0; i < myPlot.size(); ++i) myPlot[i].Draw("SAME");
0148 }
0149 else
0150 {
0151 myPlot[0].Draw("PE1");
0152 for (unsigned int i = 0; i < myPlot.size(); ++i) myPlot[i].Draw("PE1SAME");
0153 }
0154
0155 TPaveText *pt;
0156 pt = new TPaveText(0.15,0.84,0.95,0.89, "NDC");
0157 pt->SetFillColor(0);
0158 pt->SetFillStyle(0);
0159 pt->SetTextFont(42);
0160 TText *pt_LaTex;
0161 pt->AddText("#it{#bf{sPHENIX}} Internal, #sqrt{s} = 200 GeV, pp");
0162 pt->SetBorderSize(0);
0163 pt->Draw();
0164
0165 legend->SetFillColor(0);
0166 legend->SetFillStyle(0);
0167 legend->SetBorderSize(0);
0168 legend->Draw();
0169
0170 gPad->Modified();
0171
0172 string extensions[] = {".C", ".pdf", ".png"};
0173 for (auto extension : extensions)
0174 {
0175 string output = plotPath + "/" + plotName + extension;
0176 c1->SaveAs(output.c_str());
0177 }
0178 }
0179 bool isInRange(float min, float value, float max)
0180 {
0181 return min <= value && value <= max;
0182 }
0183
0184 void processData()
0185 {
0186 string fileName = "/gpfs/mnt/gpfs02/sphenix/user/cdean/software/TrackingAnalysis/Physics_Val_TF/data/macro/evalFile_00079516_withDCA.root";
0187
0188 TFile *file = new TFile(fileName.c_str());
0189 TNtuple* data = (TNtuple*)file->Get("ntp_track");
0190
0191 float track_eta; data->SetBranchAddress("eta", &track_eta);
0192 float track_phi; data->SetBranchAddress("phi", &track_phi);
0193 float track_pT; data->SetBranchAddress("pt", &track_pT);
0194 float track_charge; data->SetBranchAddress("charge", &track_charge);
0195 float track_nmaps; data->SetBranchAddress("nmaps", &track_nmaps);
0196 float track_dedx; data->SetBranchAddress("dedx", &track_dedx);
0197 float track_expected_proton_dedx; data->SetBranchAddress("prdedx", &track_expected_proton_dedx);
0198 float track_expected_pion_dedx; data->SetBranchAddress("pidedx", &track_expected_pion_dedx);
0199 float track_expected_kaon_dedx; data->SetBranchAddress("kdedx", &track_expected_kaon_dedx);
0200 float track_dca3dxy; data->SetBranchAddress("dca3dxy", &track_dca3dxy);
0201
0202 int tmp = 0;
0203 int barWidth = 50;
0204 uint64_t nEntries = data->GetEntries();
0205
0206 for (uint64_t i = 0; i < nEntries; i++)
0207 {
0208 if (tmp != (int)100*i/nEntries)
0209 {
0210 tmp = (int)100*i/nEntries;
0211 if ((tmp%1) == 0)
0212 {
0213 cout << "[";
0214 int pos = barWidth * tmp/100;
0215 for (int k = 0; k < barWidth; ++k)
0216 {
0217 if (k < pos) cout << "=";
0218 else if (k == pos) cout << ">";
0219 else cout << " ";
0220 }
0221 cout << "] " << tmp << " %\r";
0222 cout.flush();
0223 }
0224 }
0225
0226 data->GetEntry(i);
0227
0228 if (track_nmaps < 2) continue;
0229 if (abs(track_dca3dxy) > 0.1) continue;
0230
0231 bool accept_proton = isInRange(0.4, track_dedx/track_expected_proton_dedx, 1.8);
0232 bool accept_pion = isInRange(0.01, track_dedx/track_expected_pion_dedx, 1.2);
0233 bool accept_kaon = isInRange(0.2, track_dedx/track_expected_kaon_dedx, 1.8);
0234
0235 if (track_charge > 0)
0236 {
0237 if (accept_proton && track_pT <= min_track_pT_proton)
0238 {
0239 pT_proton_histo.Fill(track_pT);
0240 eta_proton_histo.Fill(track_eta);
0241 phi_proton_histo.Fill(track_phi);
0242 }
0243 if (accept_pion)
0244 {
0245 pT_piplus_histo.Fill(track_pT);
0246 eta_piplus_histo.Fill(track_eta);
0247 phi_piplus_histo.Fill(track_phi);
0248 }
0249 if (accept_kaon && track_pT <= min_track_pT_kaon)
0250 {
0251 pT_Kplus_histo.Fill(track_pT);
0252 eta_Kplus_histo.Fill(track_eta);
0253 phi_Kplus_histo.Fill(track_phi);
0254 }
0255 }
0256 else
0257 {
0258 if (accept_proton && track_pT <= min_track_pT_proton)
0259 {
0260 pT_antiproton_histo.Fill(track_pT);
0261 eta_antiproton_histo.Fill(track_eta);
0262 phi_antiproton_histo.Fill(track_phi);
0263 pT_proton_ratio.Fill(track_pT);
0264 eta_proton_ratio.Fill(track_eta);
0265 phi_proton_ratio.Fill(track_phi);
0266 }
0267 if (accept_pion)
0268 {
0269 pT_piminus_histo.Fill(track_pT);
0270 eta_piminus_histo.Fill(track_eta);
0271 phi_piminus_histo.Fill(track_phi);
0272 pT_pion_ratio.Fill(track_pT);
0273 eta_pion_ratio.Fill(track_eta);
0274 phi_pion_ratio.Fill(track_phi);
0275 }
0276 if (accept_kaon && track_pT <= min_track_pT_kaon)
0277 {
0278 pT_Kminus_histo.Fill(track_pT);
0279 eta_Kminus_histo.Fill(track_eta);
0280 phi_Kminus_histo.Fill(track_phi);
0281 pT_Kaon_ratio.Fill(track_pT);
0282 eta_Kaon_ratio.Fill(track_eta);
0283 phi_Kaon_ratio.Fill(track_phi);
0284 }
0285 }
0286 }
0287
0288 cout << "[";
0289 int pos = barWidth * tmp;
0290 for (int i = 0; i < barWidth; ++i)
0291 {
0292 if (i < pos) cout << "=";
0293 else if (i == pos) cout << ">";
0294 else cout << " ";
0295 }
0296 cout << "] 100 %\r";
0297 cout.flush();
0298 cout<<endl;
0299 }
0300
0301 void charged_hadron_ratios()
0302 {
0303 processData();
0304
0305 pT_proton_ratio.Sumw2();
0306 eta_proton_ratio.Sumw2();
0307 phi_proton_ratio.Sumw2();
0308
0309 pT_pion_ratio.Sumw2();
0310 eta_pion_ratio.Sumw2();
0311 phi_pion_ratio.Sumw2();
0312
0313 pT_Kaon_ratio.Sumw2();
0314 eta_Kaon_ratio.Sumw2();
0315 phi_Kaon_ratio.Sumw2();
0316
0317 pT_proton_histo.Sumw2();
0318 eta_proton_histo.Sumw2();
0319 phi_proton_histo.Sumw2();
0320
0321 pT_piplus_histo.Sumw2();
0322 eta_piplus_histo.Sumw2();
0323 phi_piplus_histo.Sumw2();
0324
0325 pT_Kplus_histo.Sumw2();
0326 eta_Kplus_histo.Sumw2();
0327 phi_Kplus_histo.Sumw2();
0328
0329 pT_antiproton_histo.Sumw2();
0330 eta_antiproton_histo.Sumw2();
0331 phi_antiproton_histo.Sumw2();
0332
0333 pT_piminus_histo.Sumw2();
0334 eta_piminus_histo.Sumw2();
0335 phi_piminus_histo.Sumw2();
0336
0337 pT_Kminus_histo.Sumw2();
0338 eta_Kminus_histo.Sumw2();
0339 phi_Kminus_histo.Sumw2();
0340
0341 vector<TH1F> pT_proton_histos = {pT_antiproton_histo, pT_proton_histo};
0342 saveManyPlots(pT_proton_histos, "pT_proton_histos", true, {"#bar{p}", "p"});
0343
0344 vector<TH1F> pT_pion_histos = {pT_piminus_histo, pT_piplus_histo};
0345 saveManyPlots(pT_pion_histos, "pT_pion_histos", true, {"#pi^{-}", "#pi^{+}"});
0346
0347 vector<TH1F> pT_Kaon_histos = {pT_Kminus_histo, pT_Kplus_histo};
0348 saveManyPlots(pT_Kaon_histos, "pT_Kaon_histos", true, {"K^{-}", "K^{+}"});
0349
0350 vector<TH1F> eta_proton_histos = {eta_antiproton_histo, eta_proton_histo};
0351 saveManyPlots(eta_proton_histos, "eta_proton_histos", false, {"#bar{p}", "p"});
0352
0353 vector<TH1F> eta_pion_histos = {eta_piminus_histo, eta_piplus_histo};
0354 saveManyPlots(eta_pion_histos, "eta_pion_histos", false, {"#pi^{-}", "#pi^{+}"});
0355
0356 vector<TH1F> eta_Kaon_histos = {eta_Kminus_histo, eta_Kplus_histo};
0357 saveManyPlots(eta_Kaon_histos, "eta_Kaon_histos", false, {"K^{-}", "K^{+}"});
0358
0359 vector<TH1F> phi_proton_histos = {phi_antiproton_histo, phi_proton_histo};
0360 saveManyPlots(phi_proton_histos, "phi_proton_histos", false, {"#bar{p}", "p"});
0361
0362 vector<TH1F> phi_pion_histos = {phi_piminus_histo, phi_piplus_histo};
0363 saveManyPlots(phi_pion_histos, "phi_pion_histos", false, {"#pi^{-}", "#pi^{+}"});
0364
0365 vector<TH1F> phi_Kaon_histos = {phi_Kminus_histo, phi_Kplus_histo};
0366 saveManyPlots(phi_Kaon_histos, "phi_Kaon_histos", false, {"K^{-}", "K^{+}"});
0367
0368 pT_proton_ratio.Divide(&pT_proton_histo);
0369 eta_proton_ratio.Divide(&eta_proton_histo);
0370 phi_proton_ratio.Divide(&phi_proton_histo);
0371
0372 pT_pion_ratio.Divide(&pT_piplus_histo);
0373 eta_pion_ratio.Divide(&eta_piplus_histo);
0374 phi_pion_ratio.Divide(&phi_piplus_histo);
0375
0376 pT_Kaon_ratio.Divide(&pT_Kplus_histo);
0377 eta_Kaon_ratio.Divide(&eta_Kplus_histo);
0378 phi_Kaon_ratio.Divide(&phi_Kplus_histo);
0379
0380 savePlots(pT_proton_ratio, "pT_proton_ratio");
0381 savePlots(eta_proton_ratio, "eta_proton_ratio");
0382 savePlots(phi_proton_ratio, "phi_proton_ratio");
0383
0384 savePlots(pT_pion_ratio, "pT_pion_ratio");
0385 savePlots(eta_pion_ratio, "eta_pion_ratio");
0386 savePlots(phi_pion_ratio, "phi_pion_ratio");
0387
0388 savePlots(pT_Kaon_ratio, "pT_Kaon_ratio");
0389 savePlots(eta_Kaon_ratio, "eta_Kaon_ratio");
0390 savePlots(phi_Kaon_ratio, "phi_Kaon_ratio");
0391 }
