analysis_code/final_asymmetries/CreatePhenixPlot.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 CreatePhenixPlot()
0160 {
0161   std::string plots_folder_pt = "pt_asymmetries/PHENIX/";
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
0170   float pTBounds[nParticles][2] = {{-0.026, 0.026}, {-0.09, 0.08}};
0171
0172   for (int iP = 0; iP < nParticles; iP++) {
0173     std::stringstream inputfilename;
0174     inputfilename << inputfolder << "asym_summary_" << particles[iP] << "_pT" << ".txt";
0175     std::stringstream canvas_name;
0176     canvas_name << "canvas_asym_" << particles[iP] << "_pT";
0177     std::vector<double> x, y, estat, esyst, exsyst;
0178     parse_file(inputfilename.str(), x, y, estat, esyst, exsyst, 0.20);
0179
0180     // Parse phenix plots
0181     std::stringstream th_filename;
0182     th_filename << inputfolder_phenix << "/phenix_asym_" << particles[iP] << ".csv";
0183     std::vector<double> x_ph, y_ph, estat_ph, esyst_ph, exsyst_ph;
0184     parse_phenix_file(th_filename.str(), x_ph, y_ph, estat_ph, esyst_ph, exsyst_ph, 0.20);
0185
0186     const int n = x.size();
0187     auto gSyst = new TGraphErrors(n);
0188     auto gStat = new TGraphErrors(n);
0189     auto gSyst_ph = new TGraphErrors();
0190     auto gStat_ph = new TGraphErrors();
0191     int i_th = 0;
0192     for (int i = (iP == 0 ? 0 : 1); i < (iP == 0 ? n-1 : n); ++i) {
0193       gSyst->SetPoint(i, x[i], y[i]);
0194       gSyst->SetPointError(i, exsyst[i], esyst[i]);
0195
0196       gStat->SetPoint(i, x[i], y[i]);
0197       gStat->SetPointError(i, 0.0, estat[i]);
0198     }
0199     const int n_ph = x_ph.size();
0200     for (int i = 0; i < n_ph; ++i) {
0201       gSyst_ph->SetPoint(i, x_ph[i], y_ph[i]);
0202       gSyst_ph->SetPointError(i, exsyst[0], esyst_ph[i]);
0203
0204       gStat_ph->SetPoint(i, x_ph[i], y_ph[i]);
0205       gStat_ph->SetPointError(i, 0.0, estat_ph[i]);
0206     }
0207
0208     SetsPhenixStyle();
0209
0210     TCanvas *canvas = new TCanvas(canvas_name.str().c_str(), "", 1600, 900);
0211     gStat->SetTitle(";p_{T} [GeV];A_{N}");
0212
0213     double y_bound = 0;
0214     if (iP == 0){
0215       y_bound = GetAbsMaxInRange(gStat, 0.5, 10);
0216     } else {
0217       y_bound = GetAbsMaxInRange(gStat, 1.5, 12);
0218     }
0219
0220     int color = (iP == 0 ? kRed : kBlue);
0221
0222     gStat->SetLineWidth(2);
0223     gStat->SetLineColor(color);
0224     gStat->SetMarkerColor(color);
0225     gStat->SetMarkerStyle(kFullSquare);
0226     gStat->SetMarkerSize(2.4);
0227     // gStat->SetMinimum(-2 * y_bound);
0228     // gStat->SetMaximum(2 * y_bound);
0229     gStat->SetMinimum(pTBounds[iP][0]);
0230     gStat->SetMaximum(pTBounds[iP][1]);
0231     gStat->Draw("AP E1");
0232
0233     // Style systematic boxes (transparent fill)
0234     if (iP == 0) {
0235       gSyst->SetFillColorAlpha(kPink-9, 0.45); // alpha controls transparency
0236       gSyst->SetLineColor(kPink-9);
0237     } else {
0238       gSyst->SetFillColorAlpha(kAzure-4, 0.45); // alpha controls transparency
0239       gSyst->SetLineColor(kAzure-4);
0240     }
0241     gSyst->SetMarkerStyle(0);
0242     gSyst->SetMarkerColor(color);
0243     gSyst->Draw("P E2 SAME");
0244
0245     std::cout << "PHENIX stats:" << std::endl;
0246     gStat_ph->Print();
0247     std::cout << "sPHENIX stats:" << std::endl;
0248     gStat->Print();
0249
0250
0251
0252     gSyst_ph->Print();
0253
0254     gStat_ph->SetLineWidth(2);
0255     gStat_ph->SetLineColor(kGray+2);
0256     gStat_ph->SetMarkerColor(kGray+2);
0257     gStat_ph->SetMarkerStyle(59);
0258     gStat_ph->SetMarkerSize(2.6);
0259     gStat_ph->SetMinimum(-2 * y_bound);
0260     gStat_ph->SetMaximum(2 * y_bound);
0261     gStat_ph->Draw("P E1 SAME");
0262
0263     gSyst_ph->SetMarkerStyle(0);
0264     gSyst_ph->SetMarkerColor(kGray+2);
0265     gSyst_ph->SetFillColorAlpha(kGray, 0.45); // alpha controls transparency
0266     gSyst_ph->SetLineColor(kGray);
0267     gSyst_ph->Draw("P E2 SAME");
0268
0269     if (iP == 0) gStat->GetXaxis()->SetLimits(1.0, 20.0);
0270     else gStat->GetXaxis()->SetLimits(2.0, 20.0);
0271
0272     TLegend *legend = new TLegend(0.2, 0.25, 0.5, 0.35);
0273     legend->SetFillColor(kOrange);
0274     legend->SetTextSize(0.04);
0275     legend->SetBorderSize(0);
0276     legend->SetFillStyle(0);
0277     legend->AddEntry(gStat, "sPHENIX 2024, |#eta| < 2.0");
0278     legend->AddEntry(gStat_ph, "PHENIX PRD 103 (2021) 052009, |#eta| < 0.35");
0279     legend->Draw();
0280
0281     gPad->Modified();
0282     gPad->Update();
0283     double min_x = gPad->GetUxmin();
0284     double max_x = gPad->GetUxmax();
0285
0286     // Add "sPHENIX internal"
0287     TPad *p = new TPad("p","p",0.,0.,1.,1.); p->SetFillStyle(0); p->Draw(); p->cd();
0288     TBox *whiteBox = new TBox(0.17, 0.72, 0.46, 0.90);
0289     whiteBox->Draw();
0290     canvas->cd();
0291     whiteBox->SetFillColorAlpha(kWhite, 1);
0292     std::stringstream stream;
0293     stream.str("");
0294     TLatex latex;
0295     latex.SetNDC();
0296     latex.SetTextColor(kBlack);
0297     latex.DrawLatex(0.22, 0.85, "#font[72]{sPHENIX} Internal");
0298     latex.DrawLatex(0.22, 0.79, "p^{#uparrow}+p #sqrt{s} = 200 GeV");
0299     latex.SetTextSize(0.03);
0300     latex.DrawLatex(0.19, 0.73, "7% polarization scale uncertainty not shown");
0301     latex.SetTextSize(0.07);
0302
0303     if (iP == 0) {
0304       latex.DrawLatex(0.57, 0.84, "p^{#uparrow}+p #rightarrow #pi^{0} X");
0305     } else if (iP == 1) {
0306       latex.DrawLatex(0.57, 0.84, "p^{#uparrow}+p #rightarrow #eta X");
0307     }
0308
0309     TLine *tline = new TLine();
0310     tline->SetLineWidth(2);
0311     tline->SetLineColor(kBlack);
0312     tline->SetLineStyle(kDashed);
0313     tline->DrawLine(min_x, 0, max_x, 0);
0314
0315     canvas->Update();
0316     canvas->Draw();
0317     canvas->SaveAs((plots_folder_pt +  "/" + canvas_name.str() + ".png").c_str());
0318     canvas->SaveAs((plots_folder_pt +  "/" + canvas_name.str() + ".pdf").c_str());
0319   }
0320   gSystem->Exit(0);
0321 }