Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2026-04-04 08:08:14

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 CreateTheoryPlot()
0117 {
0118   std::string plots_folder_pt = "pt_asymmetries/THEORY/";
0119   std::string plots_folder_eta = "eta_asymmetries/THEORY/";
0120   std::string plots_folder_xf = "xf_asymmetries/THEORY/";
0121 
0122   gSystem->Exec(("mkdir -p " + plots_folder_pt).c_str());
0123   gSystem->Exec(("mkdir -p " + plots_folder_eta).c_str());
0124   gSystem->Exec(("mkdir -p " + plots_folder_xf).c_str());
0125 
0126   
0127   const std::string inputfolder = "asym_values/";
0128   const std::string inputfolder_theory_qgq = "sPHENIX26_AN/";
0129   const int nParticles = 2;
0130   const std::string particles[nParticles] = {"pi0", "eta"};
0131   const int nDirections = 2;
0132   const std::string directions[nDirections] = {"forward", "backward"};
0133   for (int iP = 0; iP < nParticles; iP++) {
0134     for (int iDir = 0; iDir < nDirections; iDir++) {
0135       std::stringstream inputfilename;
0136       inputfilename << inputfolder << "asym_summary_" << particles[iP] << "_" << directions[iDir] << "_pT" << ".txt";
0137       std::stringstream canvas_name;
0138       canvas_name << "canvas_asym_" << particles[iP] << "_" << directions[iDir] << "_pT";
0139       std::vector<double> x, y, estat, esyst, exsyst;
0140       parse_file(inputfilename.str(), x, y, estat, esyst, exsyst, 0.12);
0141 
0142       // Parse theory plots
0143       std::stringstream th_filename;
0144       th_filename << inputfolder_theory_qgq << "/" << particles[iP] << "_AN_vs_pT-" << directions[iDir] << "-200.dat";
0145       std::vector<double> x_th, y_th, ey_th;
0146       parse_theory_file(th_filename.str(), x_th, y_th, ey_th);
0147 
0148       const int n = x.size();
0149       auto gSyst = new TGraphErrors(n);
0150       auto gStat = new TGraphErrors(n);
0151       auto gTh_qgq = new TGraphErrors(n-1);
0152       int i_th = 0;
0153       for (int i = 0; i < n; ++i) {
0154         gSyst->SetPoint(i, x[i], y[i]);
0155         gSyst->SetPointError(i, exsyst[i], esyst[i]);
0156         
0157         gStat->SetPoint(i, x[i], y[i]);
0158         gStat->SetPointError(i, 0.0, estat[i]);
0159 
0160         if (iP == 0 && x[i] < 10 ||
0161             iP == 1 && x[i] > 2) {
0162           gTh_qgq->SetPoint(i_th, x_th[i_th], y_th[i_th]);
0163           gTh_qgq->SetPointError(i_th, 0, ey_th[i_th]);
0164           i_th++;
0165         }
0166       }
0167       gStat->Print();
0168       gSyst->Print();
0169       gTh_qgq->Print();
0170       SetsPhenixStyle();
0171 
0172       TCanvas *canvas = new TCanvas(canvas_name.str().c_str(), "", 1600, 900);
0173       gStat->SetTitle(";p_{T} GeV;A_{N}");
0174 
0175       double y_bound = 0;
0176       if (iP == 0){
0177         y_bound = GetAbsMaxInRange(gStat, 0.5, 10);
0178       } else {
0179         y_bound = GetAbsMaxInRange(gStat, 1.5, 12);
0180       }
0181 
0182       int color = (iP == 0 ? kRed : kBlue);
0183 
0184       gStat->SetLineWidth(2);
0185       gStat->SetLineColor(color);
0186       gStat->SetMarkerColor(color);
0187       gStat->SetMarkerStyle(kFullCircle);
0188       gStat->SetMarkerSize(2.2);
0189       gStat->SetMinimum(-2 * y_bound);
0190       gStat->SetMaximum(2 * y_bound);
0191       gStat->Draw("AP E1");
0192       
0193       // Style systematic boxes (transparent fill)
0194       if (iP == 0) {
0195         gSyst->SetFillColorAlpha(kPink-9, 0.35); // alpha controls transparency
0196         gSyst->SetLineColor(kPink-9);
0197       } else {
0198         gSyst->SetFillColorAlpha(kAzure-4, 0.35); // alpha controls transparency
0199         gSyst->SetLineColor(kAzure-4);
0200       }
0201       gSyst->SetMarkerStyle(0);
0202       gSyst->Draw("P E2 SAME");
0203 
0204       // Theory plots
0205       gTh_qgq->SetLineWidth(3);
0206       gTh_qgq->SetLineStyle(5);
0207       gTh_qgq->SetLineColor(kGreen);
0208       gTh_qgq->SetMarkerColor(kGreen);
0209       gTh_qgq->SetMarkerSize(0.0);
0210       gTh_qgq->Draw("SAME");
0211       
0212       if (iP == 0) gStat->GetXaxis()->SetLimits(1.0, 10.0);
0213       else gStat->GetXaxis()->SetLimits(2.0, 13.0);
0214 
0215       TLegend *legend = new TLegend(0.2, 0.25, 0.4, 0.35);
0216       legend->SetTextSize(0.04);
0217       legend->SetBorderSize(0);
0218       legend->SetFillStyle(0);
0219       legend->AddEntry(gTh_qgq, "qgq contribution");
0220       legend->Draw();
0221 
0222       gPad->Modified();
0223       gPad->Update();
0224       double min_x = gPad->GetUxmin();
0225       double max_x = gPad->GetUxmax();
0226 
0227       // Add "sPHENIX internal"
0228       TPad *p = new TPad("p","p",0.,0.,1.,1.); p->SetFillStyle(0); p->Draw(); p->cd();
0229       TBox *whiteBox = new TBox(0.17, 0.72, 0.46, 0.90);
0230       whiteBox->Draw();
0231       canvas->cd();
0232       whiteBox->SetFillColorAlpha(kWhite, 1);
0233       std::stringstream stream;
0234       stream.str("");
0235       TLatex latex;
0236       latex.SetNDC();
0237       latex.SetTextColor(kBlack);
0238       latex.DrawLatex(0.22, 0.85, "#font[72]{sPHENIX} Internal");
0239       latex.DrawLatex(0.22, 0.75, "p^{#uparrow}+p #sqrt{s} = 200 GeV");
0240       latex.SetTextSize(0.03);
0241       latex.DrawLatex(0.19, 0.68, "7% polarization scale uncertainty not shown");
0242       latex.SetTextSize(0.05);
0243 
0244       if (iP == 0 && iDir == 0) {
0245         latex.DrawLatex(0.5, 0.85, "p^{#uparrow}+p #rightarrow #pi^{0} X,   x_{F}>0");
0246       } else if (iP == 1 && iDir == 0) {
0247         latex.DrawLatex(0.5, 0.85, "p^{#uparrow}+p #rightarrow #eta X,   x_{F}>0");
0248       } else if (iP == 0 && iDir == 1) {
0249         latex.DrawLatex(0.5, 0.85, "p^{#uparrow}+p #rightarrow #pi^{0} X,   x_{F}<0");
0250       } else if (iP == 1 && iDir == 1) {
0251         latex.DrawLatex(0.5, 0.85, "p^{#uparrow}+p #rightarrow #eta X,   x_{F}<0");
0252       }
0253 
0254       TLine *tline = new TLine();
0255       tline->SetLineWidth(2);
0256       tline->SetLineColor(kBlack);
0257       tline->SetLineStyle(kDashed);
0258       tline->DrawLine(min_x, 0, max_x, 0);
0259 
0260       canvas->Update();
0261       canvas->Draw();
0262       canvas->SaveAs((plots_folder_pt +  "/" + canvas_name.str() + ".png").c_str());
0263       canvas->SaveAs((plots_folder_pt +  "/" + canvas_name.str() + ".pdf").c_str());
0264     }
0265     {
0266       std::stringstream inputfilename;
0267       inputfilename << inputfolder << "asym_summary_" << particles[iP] << "_eta" << ".txt";
0268       std::stringstream canvas_name;
0269       canvas_name << "canvas_asym_" << particles[iP] << "_eta";
0270       std::vector<double> x, y, estat, esyst, exsyst;
0271 
0272       parse_file(inputfilename.str(), x, y, estat, esyst, exsyst, 0.05);
0273 
0274       // Parse theory plots
0275       std::stringstream th_filename;
0276       th_filename << inputfolder_theory_qgq << "/" << particles[iP] << "_AN_vs_eta-200.dat";
0277       std::vector<double> x_th, y_th, ey_th;
0278       parse_theory_file(th_filename.str(), x_th, y_th, ey_th);
0279 
0280       const int n = x.size();
0281       auto gSyst = new TGraphErrors(n);
0282       auto gStat = new TGraphErrors(n);
0283       auto gTh_qgq = new TGraphErrors(n);
0284 
0285       for (int i = 0; i < n; ++i) {
0286         gSyst->SetPoint(i, x[i], y[i]);
0287         gSyst->SetPointError(i, exsyst[i], esyst[i]);
0288         
0289         gStat->SetPoint(i, x[i], y[i]);
0290         gStat->SetPointError(i, 0.0, estat[i]);
0291 
0292         gTh_qgq->SetPoint(i, x_th[i], y_th[i]);
0293         gTh_qgq->SetPointError(i, 0, ey_th[i]);
0294       }
0295       gStat->Print();
0296       gSyst->Print();
0297       SetsPhenixStyle();
0298 
0299       TCanvas *canvas = new TCanvas(canvas_name.str().c_str(), "", 1600, 900);
0300       gStat->SetTitle(";#eta;A_{N}");
0301 
0302       double y_bound = 0;
0303       y_bound = GetAbsMaxInRange(gStat, -3, 3);
0304 
0305       int color = (iP == 0 ? kRed : kBlue);
0306 
0307       gStat->SetLineWidth(2);
0308       gStat->SetLineColor(color);
0309       gStat->SetMarkerColor(color);
0310       gStat->SetMarkerStyle(kFullCircle);
0311       gStat->SetMarkerSize(2.2);
0312       gStat->SetMinimum(-2 * y_bound);
0313       gStat->SetMaximum(2 * y_bound);
0314       gStat->Draw("AP E1");
0315 
0316       // Theory plots
0317       gTh_qgq->SetLineWidth(3);
0318       gTh_qgq->SetLineStyle(5);
0319       gTh_qgq->SetLineColor(kGreen);
0320       gTh_qgq->SetMarkerColor(kGreen);
0321       gTh_qgq->SetMarkerSize(0.0);
0322       gTh_qgq->Draw("SAME");
0323 
0324        // Style systematic boxes (transparent fill)
0325       if (iP == 0) {
0326         gSyst->SetFillColorAlpha(kPink-9, 0.35); // alpha controls transparency
0327         gSyst->SetLineColor(kPink-9);
0328       } else {
0329         gSyst->SetFillColorAlpha(kAzure-4, 0.35); // alpha controls transparency
0330         gSyst->SetLineColor(kAzure-4);
0331       }
0332       gSyst->SetMarkerStyle(0);
0333       gSyst->Draw("P E2 SAME");
0334 
0335       TLegend *legend = new TLegend(0.2, 0.25, 0.4, 0.35);
0336       legend->SetTextSize(0.04);
0337       legend->SetBorderSize(0);
0338       legend->SetFillStyle(0);
0339       legend->AddEntry(gTh_qgq, "qgq contribution");
0340       legend->Draw();
0341 
0342       gPad->Modified();
0343       gPad->Update();
0344       double min_x = gPad->GetUxmin();
0345       double max_x = gPad->GetUxmax();
0346 
0347       // Add "sPHENIX internal"
0348       TPad *p = new TPad("p","p",0.,0.,1.,1.); p->SetFillStyle(0); p->Draw(); p->cd();
0349       TBox *whiteBox = new TBox(0.17, 0.72, 0.46, 0.90);
0350       whiteBox->Draw();
0351       canvas->cd();
0352       whiteBox->SetFillColorAlpha(kWhite, 1);
0353       std::stringstream stream;
0354       stream.str("");
0355       TLatex latex;
0356       latex.SetNDC();
0357       latex.SetTextColor(kBlack);
0358       latex.DrawLatex(0.22, 0.85, "#font[72]{sPHENIX} Internal");
0359       latex.DrawLatex(0.22, 0.75, "p^{#uparrow}+p #sqrt{s} = 200 GeV");
0360       latex.SetTextSize(0.03);
0361       latex.DrawLatex(0.19, 0.68, "7% polarization scale uncertainty not shown");
0362       latex.SetTextSize(0.05);
0363 
0364       if (iP == 0) {
0365         latex.DrawLatex(0.5, 0.85, "p^{#uparrow}+p #rightarrow #pi^{0} X");
0366       } else {
0367         latex.DrawLatex(0.5, 0.85, "p^{#uparrow}+p #rightarrow #eta X");
0368       }
0369 
0370       TLine *tline = new TLine();
0371       tline->SetLineWidth(2);
0372       tline->SetLineColor(kBlack);
0373       tline->SetLineStyle(kDashed);
0374       tline->DrawLine(min_x, 0, max_x, 0);
0375 
0376       canvas->Update();
0377       canvas->Draw();
0378       canvas->SaveAs((plots_folder_eta +  "/" + canvas_name.str() + ".png").c_str());
0379       canvas->SaveAs((plots_folder_eta +  "/" + canvas_name.str() + ".pdf").c_str());
0380     }
0381     {
0382       std::stringstream inputfilename;
0383       inputfilename << inputfolder << "asym_summary_" << particles[iP] << "_xf" << ".txt";
0384       std::stringstream canvas_name;
0385       canvas_name << "canvas_asym_" << particles[iP] << "_xf";
0386       std::vector<double> x, y, estat, esyst, exsyst;
0387 
0388       parse_file(inputfilename.str(), x, y, estat, esyst, exsyst, 0.002);
0389 
0390       // Parse theory plots
0391       std::stringstream th_filename;
0392       th_filename << inputfolder_theory_qgq << "/" << particles[iP] << "_AN_vs_xF-200.dat";
0393       std::vector<double> x_th, y_th, ey_th;
0394       parse_theory_file(th_filename.str(), x_th, y_th, ey_th);
0395 
0396       const int n = x.size();
0397       auto gSyst = new TGraphErrors(n);
0398       auto gStat = new TGraphErrors(n);
0399       auto gTh_qgq = new TGraphErrors(n);
0400 
0401       int ind_th = 0;
0402       for (int i = 0; i < n; ++i) {
0403         gSyst->SetPoint(i, x[i], y[i]);
0404         gSyst->SetPointError(i, exsyst[i], esyst[i]);
0405         
0406         gStat->SetPoint(i, x[i], y[i]);
0407         gStat->SetPointError(i, 0.0, estat[i]);
0408 
0409         gTh_qgq->SetPoint(i, x_th[i], y_th[i]);
0410         gTh_qgq->SetPointError(i, 0, ey_th[i]);
0411       }
0412       gStat->Print();
0413       gSyst->Print();
0414       SetsPhenixStyle();
0415 
0416       TCanvas *canvas = new TCanvas(canvas_name.str().c_str(), "", 1600, 900);
0417       gStat->SetTitle(";x_{F};A_{N}");
0418 
0419       double y_bound = 0;
0420       y_bound = GetAbsMaxInRange(gStat, -0.3, 0.3);
0421 
0422       int color = (iP == 0 ? kRed : kBlue);
0423 
0424       gStat->SetLineWidth(2);
0425       gStat->SetLineColor(color);
0426       gStat->SetMarkerColor(color);
0427       gStat->SetMarkerStyle(kFullCircle);
0428       gStat->SetMarkerSize(2.2);
0429       gStat->SetMinimum(-2 * y_bound);
0430       gStat->SetMaximum(2 * y_bound);
0431       gStat->Draw("AP E1");
0432 
0433        // Style systematic boxes (transparent fill)
0434       if (iP == 0) {
0435         gSyst->SetFillColorAlpha(kPink-9, 0.35); // alpha controls transparency
0436         gSyst->SetLineColor(kPink-9);
0437       } else {
0438         gSyst->SetFillColorAlpha(kAzure-4, 0.35); // alpha controls transparency
0439         gSyst->SetLineColor(kAzure-4);
0440       }
0441       gSyst->SetMarkerStyle(0);
0442       gSyst->Draw("P E2 SAME");
0443 
0444       // Theory plots
0445       gTh_qgq->SetLineWidth(3);
0446       gTh_qgq->SetLineStyle(5);
0447       gTh_qgq->SetLineColor(kGreen);
0448       gTh_qgq->SetMarkerColor(kGreen);
0449       gTh_qgq->SetMarkerSize(0.0);
0450       gTh_qgq->Draw("SAME");
0451 
0452       TLegend *legend = new TLegend(0.2, 0.25, 0.4, 0.35);
0453       legend->SetTextSize(0.04);
0454       legend->SetBorderSize(0);
0455       legend->SetFillStyle(0);
0456       legend->AddEntry(gTh_qgq, "qgq contribution");
0457       legend->Draw();
0458 
0459       gPad->Modified();
0460       gPad->Update();
0461       double min_x = gPad->GetUxmin();
0462       double max_x = gPad->GetUxmax();
0463 
0464       // Add "sPHENIX internal"
0465       TPad *p = new TPad("p","p",0.,0.,1.,1.); p->SetFillStyle(0); p->Draw(); p->cd();
0466       TBox *whiteBox = new TBox(0.17, 0.72, 0.46, 0.90);
0467       whiteBox->Draw();
0468       canvas->cd();
0469       whiteBox->SetFillColorAlpha(kWhite, 1);
0470       std::stringstream stream;
0471       stream.str("");
0472       TLatex latex;
0473       latex.SetNDC();
0474       latex.SetTextColor(kBlack);
0475       latex.DrawLatex(0.22, 0.85, "#font[72]{sPHENIX} Internal");
0476       latex.DrawLatex(0.22, 0.75, "p^{#uparrow}+p #sqrt{s} = 200 GeV");
0477       latex.SetTextSize(0.03);
0478       latex.DrawLatex(0.19, 0.68, "7% polarization scale uncertainty not shown");
0479       latex.SetTextSize(0.05);
0480 
0481       if (iP == 0) {
0482         latex.DrawLatex(0.5, 0.85, "p^{#uparrow}+p #rightarrow #pi^{0} X");
0483       } else {
0484         latex.DrawLatex(0.5, 0.85, "p^{#uparrow}+p #rightarrow #eta X");
0485       }
0486 
0487       TLine *tline = new TLine();
0488       tline->SetLineWidth(2);
0489       tline->SetLineColor(kBlack);
0490       tline->SetLineStyle(kDashed);
0491       tline->DrawLine(min_x, 0, max_x, 0);
0492 
0493       canvas->Update();
0494       canvas->Draw();
0495       canvas->SaveAs((plots_folder_xf +  "/" + canvas_name.str() + ".png").c_str());
0496       canvas->SaveAs((plots_folder_xf +  "/" + canvas_name.str() + ".pdf").c_str());
0497     }
0498   }
0499   gSystem->Exit(0);
0500 }