analysis_code/final_asymmetries/CreateSamePhenixPlot.C

0001 #include "/sphenix/u/virgilemahaut/style/sPhenixStyle_Greg.C"
0002
0003 #include <TCanvas.h>
0004 #include <TGraphErrors.h>
0005 #include <TLegend.h>
0006 #include <TAxis.h>
0007 #include <TStyle.h>
0008
0009 #include <fstream>
0010 #include <sstream>
0011 #include <string>
0012 #include <vector>
0013 #include <algorithm>
0014 #include <iostream>
0015
0016 double GetAbsMaxInRange(TGraphErrors *graph, double xMin, double xMax)
0017 {
0018   if (!graph ) return -1; // Error handling
0019
0020   double maxVal = -1;
0021   int n = graph->GetN();
0022
0023   for (int i = 0; i < n; i++)
0024   {
0025     double x, y, ey, yVal;
0026     graph->GetPoint(i, x, y);
0027     ey = graph->GetErrorY(i);
0028     yVal = std::max(std::abs(y-ey), std::abs(y+ey));
0029     if (x >= xMin && x <= xMax)
0030     {
0031       if (yVal > maxVal)
0032       {
0033         maxVal = yVal;
0034       }
0035     }
0036   }
0037   return maxVal;
0038 }
0039
0040 void parse_theory_file(const std::string& filename,
0041                        std::vector<double>& x,
0042                        std::vector<double>& y,
0043                        std::vector<double>& ey)
0044 {
0045   std::ifstream fin(filename);
0046   if (!fin.is_open()) {
0047     std::cerr << "Cannot open file: " << filename << std::endl;
0048     return;
0049   }
0050
0051   std::string line;
0052   // Skip first line  (labels)
0053   std::getline(fin, line);
0054   // Read individual entries
0055   while (std::getline(fin, line)) {
0056     // Skip empty lines
0057     if (line.empty()) continue;
0058
0059     // Replace '&' with spaces so stringstream can parse numbers
0060     std::replace(line.begin(), line.end(), '&', ' ');
0061
0062     std::stringstream ss(line);
0063     double xv, yv, eyv;
0064     if (!(ss >> xv >> yv >> eyv)) {
0065       std::cerr << "Skipping malformed line: " << line << std::endl;
0066       continue;
0067     }
0068
0069     x.push_back(xv);
0070     y.push_back(yv);
0071     ey.push_back(eyv);
0072   }
0073 }
0074
0075 void parse_file(const std::string& filename,
0076                 std::vector<double>& x,
0077                 std::vector<double>& y,
0078                 std::vector<double>& estat,
0079                 std::vector<double>& esyst,
0080                 std::vector<double>& exsyst,
0081                 const double width = 0.15
0082                 )
0083 {
0084   std::ifstream fin(filename);
0085   if (!fin.is_open()) {
0086     std::cerr << "Cannot open file: " << filename << std::endl;
0087     return;
0088   }
0089
0090   std::string line;
0091   while (std::getline(fin, line)) {
0092     // Skip empty lines
0093     if (line.empty()) continue;
0094
0095     // Replace '&' with spaces so stringstream can parse numbers
0096     std::replace(line.begin(), line.end(), '&', ' ');
0097
0098     std::stringstream ss(line);
0099     double xv, yv, statv, systv;
0100     if (!(ss >> xv >> yv >> statv >> systv)) {
0101       std::cerr << "Skipping malformed line: " << line << std::endl;
0102       continue;
0103     }
0104
0105     x.push_back(xv);
0106     y.push_back(yv);
0107     estat.push_back(statv);
0108     esyst.push_back(systv);
0109
0110     // Constant half-width for the systematic boxes in x.
0111     // Choose something visually reasonable for your spacing.
0112     exsyst.push_back(width);
0113   }
0114 }
0115
0116 void parse_phenix_file(const std::string& filename,
0117                        std::vector<double>& x,
0118                        std::vector<double>& y,
0119                        std::vector<double>& estat,
0120                        std::vector<double>& esyst,
0121                        std::vector<double>& exsyst,
0122                        const double width = 0.15
0123                 )
0124 {
0125   std::ifstream fin(filename);
0126   if (!fin.is_open()) {
0127     std::cerr << "Cannot open file: " << filename << std::endl;
0128     return;
0129   }
0130
0131   std::string line;
0132   // Skip first lines
0133   for (int i = 0; i < 8; i++) std::getline(fin, line);
0134   while (std::getline(fin, line)) {
0135     // Skip empty lines
0136     if (line.empty()) continue;
0137
0138     // Replace ',' with spaces so stringstream can parse numbers
0139     std::replace(line.begin(), line.end(), ',', ' ');
0140
0141     std::stringstream ss(line);
0142     double xv, yv, statv, systv, dummy;
0143     if (!(ss >> xv >> dummy >> dummy >> yv >> statv >> dummy >> dummy >> dummy >> dummy >> dummy >> systv)) {
0144       std::cerr << "Skipping malformed line: " << line << std::endl;
0145       continue;
0146     }
0147
0148     x.push_back(xv);
0149     y.push_back(yv);
0150     estat.push_back(statv);
0151     esyst.push_back(systv);
0152
0153     // Constant half-width for the systematic boxes in x.
0154     // Choose something visually reasonable for your spacing.
0155     exsyst.push_back(width);
0156   }
0157 }
0158
0159 void CreateSamePhenixPlot()
0160 {
0161   std::string plots_folder_pt = "pt_asymmetries/PHENIX_ACCEPTANCE/";
0162
0163   gSystem->Exec(("mkdir -p " + plots_folder_pt).c_str());
0164
0165   const std::string inputfolder = "asym_values/";
0166   const std::string inputfolder_phenix = "PHENIX21_AN/";
0167   const int nParticles = 2;
0168   const std::string particles[nParticles] = {"pi0", "eta"};
0169   for (int iP = 0; iP < nParticles; iP++) {
0170     std::stringstream inputfilename;
0171     inputfilename << inputfolder << "asym_summary_" << particles[iP] << "_pT_phenix" << ".txt";
0172     std::stringstream canvas_name;
0173     canvas_name << "canvas_asym_" << particles[iP] << "_pT";
0174     std::vector<double> x, y, estat, esyst, exsyst;
0175     parse_file(inputfilename.str(), x, y, estat, esyst, exsyst, 0.12);
0176
0177     // Parse phenix plots
0178     std::stringstream th_filename;
0179     th_filename << inputfolder_phenix << "/phenix_asym_" << particles[iP] << ".csv";
0180     std::vector<double> x_ph, y_ph, estat_ph, esyst_ph, exsyst_ph;
0181     parse_phenix_file(th_filename.str(), x_ph, y_ph, estat_ph, esyst_ph, exsyst_ph, 0.12);
0182
0183     const int n = x.size();
0184     auto gSyst = new TGraphErrors(n);
0185     auto gStat = new TGraphErrors(n);
0186     auto gSyst_ph = new TGraphErrors();
0187     auto gStat_ph = new TGraphErrors();
0188     int i_th = 0;
0189     for (int i = 0; i < n; ++i) {
0190       gSyst->SetPoint(i, x[i], y[i]);
0191       gSyst->SetPointError(i, exsyst[i], esyst[i]);
0192
0193       gStat->SetPoint(i, x[i], y[i]);
0194       gStat->SetPointError(i, 0.0, estat[i]);
0195     }
0196     const int n_ph = x_ph.size();
0197     for (int i = 0; i < n_ph; ++i) {
0198       gSyst_ph->SetPoint(i, x_ph[i], y_ph[i]);
0199       gSyst_ph->SetPointError(i, exsyst[0], esyst_ph[i]);
0200
0201       gStat_ph->SetPoint(i, x_ph[i], y_ph[i]);
0202       gStat_ph->SetPointError(i, 0.0, estat_ph[i]);
0203     }
0204
0205     SetsPhenixStyle();
0206
0207     TCanvas *canvas = new TCanvas(canvas_name.str().c_str(), "", 1600, 900);
0208     gStat->SetTitle(";p_{T} GeV;A_{N}");
0209
0210     double y_bound = 0;
0211     if (iP == 0){
0212       y_bound = GetAbsMaxInRange(gStat, 2.0, 10);
0213     } else {
0214       y_bound = GetAbsMaxInRange(gStat, 2.0, 12);
0215     }
0216
0217     int color = (iP == 0 ? kRed : kBlue);
0218
0219     gStat->SetLineWidth(2);
0220     gStat->SetLineColor(color);
0221     gStat->SetMarkerColor(color);
0222     gStat->SetMarkerStyle(kFullCircle);
0223     gStat->SetMarkerSize(2.2);
0224     gStat->SetMinimum(-2 * y_bound);
0225     gStat->SetMaximum(2 * y_bound);
0226     gStat->Draw("AP E1");
0227
0228     // Style systematic boxes (transparent fill)
0229     if (iP == 0) {
0230       gSyst->SetFillColorAlpha(kPink-9, 0.45); // alpha controls transparency
0231       gSyst->SetLineColor(kPink-9);
0232     } else {
0233       gSyst->SetFillColorAlpha(kAzure-4, 0.45); // alpha controls transparency
0234       gSyst->SetLineColor(kAzure-4);
0235     }
0236     gSyst->SetMarkerStyle(0);
0237     gSyst->SetMarkerColor(color);
0238     gSyst->Draw("P E2 SAME");
0239
0240     std::cout << "PHENIX stats:" << std::endl;
0241     gStat_ph->Print();
0242     std::cout << "sPHENIX stats:" << std::endl;
0243     gStat->Print();
0244
0245
0246
0247     gSyst_ph->Print();
0248
0249     gStat_ph->SetLineWidth(2);
0250     gStat_ph->SetLineColor(kGray+2);
0251     gStat_ph->SetMarkerColor(kGray+2);
0252     gStat_ph->SetMarkerStyle(59);
0253     gStat_ph->SetMarkerSize(2.4);
0254     gStat_ph->SetMinimum(-2 * y_bound);
0255     gStat_ph->SetMaximum(2 * y_bound);
0256     gStat_ph->Draw("P E1 SAME");
0257
0258     gSyst_ph->SetMarkerStyle(0);
0259     gSyst_ph->SetMarkerColor(kGray+2);
0260     gSyst_ph->SetFillColorAlpha(kGray, 0.45); // alpha controls transparency
0261     gSyst_ph->SetLineColor(kGray);
0262     gSyst_ph->Draw("P E2 SAME");
0263
0264     if (iP == 0) gStat->GetXaxis()->SetLimits(2.0, 20.0);
0265     else gStat->GetXaxis()->SetLimits(2.0, 20.0);
0266
0267     TLegend *legend = new TLegend(0.2, 0.25, 0.4, 0.35);
0268     legend->SetTextSize(0.04);
0269     legend->SetBorderSize(0);
0270     legend->SetFillStyle(0);
0271     legend->AddEntry(gStat, "sPHENIX 2024, |#eta| < 0.35");
0272     legend->AddEntry(gStat_ph, "PHENIX PRD 103 (2021) 052009, |#eta| < 0.35");
0273     legend->Draw();
0274
0275     gPad->Modified();
0276     gPad->Update();
0277     double min_x = gPad->GetUxmin();
0278     double max_x = gPad->GetUxmax();
0279
0280     // Add "sPHENIX internal"
0281     std::stringstream stream;
0282     stream.str("");
0283     TLatex latex;
0284     latex.SetNDC();
0285     latex.SetTextColor(kBlack);
0286     latex.DrawLatex(0.67, 0.87, "#font[72]{sPHENIX} Internal");
0287     latex.DrawLatex(0.67, 0.77, "p^{#uparrow}+p #sqrt{s} = 200 GeV");
0288     latex.SetTextSize(0.03);
0289     latex.DrawLatex(0.64, 0.7, "7% polarization scale uncertainty not shown");
0290     latex.SetTextSize(0.05);
0291
0292     if (iP == 0) {
0293       latex.DrawLatex(0.45, 0.85, "p^{#uparrow}+p #rightarrow #pi^{0} X");
0294     } else if (iP == 1) {
0295       latex.DrawLatex(0.45, 0.85, "p^{#uparrow}+p #rightarrow #eta X");
0296     }
0297
0298     TLine *tline = new TLine();
0299     tline->SetLineWidth(2);
0300     tline->SetLineColor(kBlack);
0301     tline->SetLineStyle(kDashed);
0302     tline->DrawLine(min_x, 0, max_x, 0);
0303
0304     canvas->Update();
0305     canvas->Draw();
0306     canvas->SaveAs((plots_folder_pt +  "/" + canvas_name.str() + ".png").c_str());
0307     canvas->SaveAs((plots_folder_pt +  "/" + canvas_name.str() + ".pdf").c_str());
0308   }
0309   gSystem->Exit(0);
0310 }