analysis_code/final_asymmetries/CreateSysPlot.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
0041 void parse_file(const std::string& filename,
0042                 std::vector<double>& x,
0043                 std::vector<double>& y,
0044                 std::vector<double>& estat,
0045                 std::vector<double>& esyst,
0046                 std::vector<double>& exsyst,
0047                 const double width = 0.15
0048                 )
0049 {
0050   std::ifstream fin(filename);
0051   if (!fin.is_open()) {
0052     std::cerr << "Cannot open file: " << filename << std::endl;
0053     return;
0054   }
0055
0056   std::string line;
0057   while (std::getline(fin, line)) {
0058     // Skip empty lines
0059     if (line.empty()) continue;
0060
0061     // Replace '&' with spaces so stringstream can parse numbers
0062     std::replace(line.begin(), line.end(), '&', ' ');
0063
0064     std::stringstream ss(line);
0065     double xv, yv, statv, systv;
0066     if (!(ss >> xv >> yv >> statv >> systv)) {
0067       std::cerr << "Skipping malformed line: " << line << std::endl;
0068       continue;
0069     }
0070
0071     x.push_back(xv);
0072     y.push_back(yv);
0073     estat.push_back(statv);
0074     esyst.push_back(systv);
0075
0076     // Constant half-width for the systematic boxes in x.
0077     // Choose something visually reasonable for your spacing.
0078     exsyst.push_back(width);
0079   }
0080 }
0081
0082 void CreateSysPlot()
0083 {
0084   std::string plots_folder_pt = "pt_asymmetries/SYSTEMATIC/";
0085   std::string plots_folder_eta = "eta_asymmetries/SYSTEMATIC/";
0086   std::string plots_folder_xf = "xf_asymmetries/SYSTEMATIC/";
0087
0088   gSystem->Exec(("mkdir -p " + plots_folder_pt).c_str());
0089   gSystem->Exec(("mkdir -p " + plots_folder_eta).c_str());
0090   gSystem->Exec(("mkdir -p " + plots_folder_xf).c_str());
0091
0092   const std::string inputfolder = "asym_values/";
0093   const int nParticles = 2;
0094   const std::string particles[nParticles] = {"pi0", "eta"};
0095   const int nDirections = 2;
0096   const std::string directions[nDirections] = {"forward", "backward"};
0097
0098   float pTBounds[nParticles][2] = {{-0.010, 0.020}, {-0.09, 0.07}};
0099
0100   float etaBounds[nParticles][2] = {{-0.002, 0.004}, {-0.03, 0.03}};
0101
0102   float xFBounds[nParticles][2] = {{-0.002, 0.004}, {-0.03, 0.03}};
0103
0104   for (int iP = 0; iP < nParticles; iP++) {
0105     for (int iDir = 0; iDir < nDirections; iDir++) {
0106       std::stringstream inputfilename;
0107       inputfilename << inputfolder << "asym_summary_" << particles[iP] << "_" << directions[iDir] << "_pT" << ".txt";
0108       std::stringstream canvas_name;
0109       canvas_name << "canvas_asym_" << particles[iP] << "_" << directions[iDir] << "_pT";
0110       std::vector<double> x, y, estat, esyst, exsyst;
0111
0112       parse_file(inputfilename.str(), x, y, estat, esyst, exsyst, 0.12);
0113
0114       const int n = x.size();
0115       auto gSyst = new TGraphErrors(n);
0116       auto gStat = new TGraphErrors(n);
0117
0118       for (int i = 0; i < n; ++i) {
0119         gSyst->SetPoint(i, x[i], y[i]);
0120         gSyst->SetPointError(i, exsyst[i], esyst[i]);
0121
0122         gStat->SetPoint(i, x[i], y[i]);
0123         gStat->SetPointError(i, 0.0, estat[i]);
0124       }
0125       gStat->Print();
0126       gSyst->Print();
0127       SetsPhenixStyle();
0128
0129       TCanvas *canvas = new TCanvas(canvas_name.str().c_str(), "", 1600, 900);
0130       gStat->SetTitle(";p_{T} [GeV];A_{N}");
0131
0132       double y_bound = 0;
0133       if (iP == 0){
0134         y_bound = GetAbsMaxInRange(gStat, 0.5, 10);
0135       } else {
0136         y_bound = GetAbsMaxInRange(gStat, 1.5, 12);
0137       }
0138
0139       int color = (iP == 0 ? kRed : kBlue);
0140
0141       gStat->SetLineWidth(2);
0142       gStat->SetLineColor(color);
0143       gStat->SetMarkerColor(color);
0144       gStat->SetMarkerStyle(kFullSquare);
0145       gStat->SetMarkerSize(2.4);
0146       // gStat->SetMinimum(-2 * y_bound);
0147       // gStat->SetMaximum(2 * y_bound);
0148       gStat->SetMinimum(pTBounds[iP][0]);
0149       gStat->SetMaximum(pTBounds[iP][1]);
0150       if (iP == 0) gStat->GetXaxis()->SetNdivisions(9, 5, 0);
0151       else gStat->GetXaxis()->SetNdivisions(11, 5, 0);
0152       gStat->Draw("AP E1");
0153
0154       // Style systematic boxes (transparent fill)
0155       if (iP == 0) {
0156         gSyst->SetFillColorAlpha(kPink-9, 0.35); // alpha controls transparency
0157         gSyst->SetLineColor(kPink-9);
0158       } else {
0159         gSyst->SetFillColorAlpha(kAzure-4, 0.35); // alpha controls transparency
0160         gSyst->SetLineColor(kAzure-4);
0161       }
0162       gSyst->SetMarkerStyle(0);
0163       gSyst->Draw("P E2 SAME");
0164
0165       if (iP == 0) gStat->GetXaxis()->SetLimits(1.0, 10.0);
0166       else gStat->GetXaxis()->SetLimits(2.0, 13.0);
0167
0168       // TLegend *legend = new TLegend(0.62, 0.75, 0.88, 0.88);
0169       // legend->SetBorderSize(0);
0170       // legend->SetFillStyle(0);
0171
0172       gPad->Modified();
0173       gPad->Update();
0174       double min_x = gPad->GetUxmin();
0175       double max_x = gPad->GetUxmax();
0176
0177       // Add "sPHENIX internal"
0178       TPad *p = new TPad("p","p",0.,0.,1.,1.); p->SetFillStyle(0); p->Draw(); p->cd();
0179       TBox *whiteBox = new TBox(0.185, 0.72, 0.49, 0.90);
0180       whiteBox->Draw();
0181       canvas->cd();
0182       whiteBox->SetFillColorAlpha(kWhite, 1);
0183       std::stringstream stream;
0184       stream.str("");
0185       TLatex latex;
0186       latex.SetNDC();
0187       latex.SetTextColor(kBlack);
0188       latex.DrawLatex(0.22, 0.85, "#font[72]{sPHENIX} Internal");
0189       latex.DrawLatex(0.22, 0.79, "p^{#uparrow}+p #sqrt{s} = 200 GeV");
0190       latex.SetTextSize(0.03);
0191       latex.DrawLatex(0.19, 0.73, "7% polarization scale uncertainty not shown");
0192       latex.SetTextSize(0.07);
0193       if (iP == 0 && iDir == 0) {
0194         latex.DrawLatex(0.57, 0.84, "p^{#uparrow}+p #rightarrow #pi^{0} X,   x_{F} > 0");
0195       } else if (iP == 1 && iDir == 0) {
0196         latex.DrawLatex(0.57, 0.84, "p^{#uparrow}+p #rightarrow #eta X,   x_{F} > 0");
0197       } else if (iP == 0 && iDir == 1) {
0198         latex.DrawLatex(0.57, 0.84, "p^{#uparrow}+p #rightarrow #pi^{0} X,   x_{F} < 0");
0199       } else if (iP == 1 && iDir == 1) {
0200         latex.DrawLatex(0.57, 0.84, "p^{#uparrow}+p #rightarrow #eta X,   x_{F} < 0");
0201       }
0202
0203       TLine *tline = new TLine();
0204       tline->SetLineWidth(2);
0205       tline->SetLineColor(kBlack);
0206       tline->SetLineStyle(kDashed);
0207       tline->DrawLine(min_x, 0, max_x, 0);
0208
0209       canvas->Update();
0210       canvas->Draw();
0211       canvas->SaveAs((plots_folder_pt +  "/" + canvas_name.str() + ".png").c_str());
0212       canvas->SaveAs((plots_folder_pt +  "/" + canvas_name.str() + ".pdf").c_str());
0213     }
0214     {
0215       std::stringstream inputfilename;
0216       inputfilename << inputfolder << "asym_summary_" << particles[iP] << "_eta" << ".txt";
0217       std::stringstream canvas_name;
0218       canvas_name << "canvas_asym_" << particles[iP] << "_eta";
0219       std::vector<double> x, y, estat, esyst, exsyst;
0220
0221       parse_file(inputfilename.str(), x, y, estat, esyst, exsyst, 0.05);
0222
0223       const int n = x.size();
0224       auto gSyst = new TGraphErrors(n);
0225       auto gStat = new TGraphErrors(n);
0226
0227       for (int i = 0; i < n; ++i) {
0228         gSyst->SetPoint(i, x[i], y[i]);
0229         gSyst->SetPointError(i, exsyst[i], esyst[i]);
0230
0231         gStat->SetPoint(i, x[i], y[i]);
0232         gStat->SetPointError(i, 0.0, estat[i]);
0233       }
0234       gStat->Print();
0235       gSyst->Print();
0236       SetsPhenixStyle();
0237
0238       TCanvas *canvas = new TCanvas(canvas_name.str().c_str(), "", 1600, 900);
0239       gStat->SetTitle(";#eta;A_{N}");
0240
0241       double y_bound = 0;
0242       y_bound = GetAbsMaxInRange(gStat, -3, 3);
0243
0244       int color = (iP == 0 ? kRed : kBlue);
0245
0246       gStat->SetLineWidth(2);
0247       gStat->SetLineColor(color);
0248       gStat->SetMarkerColor(color);
0249       gStat->SetMarkerStyle(kFullSquare);
0250       gStat->SetMarkerSize(2.4);
0251       // gStat->SetMinimum(-2 * y_bound);
0252       // gStat->SetMaximum(2 * y_bound);
0253       gStat->SetMinimum(etaBounds[iP][0]);
0254       gStat->SetMaximum(etaBounds[iP][1]);
0255       gStat->Draw("AP E1");
0256
0257        // Style systematic boxes (transparent fill)
0258       if (iP == 0) {
0259         gSyst->SetFillColorAlpha(kPink-9, 0.35); // alpha controls transparency
0260         gSyst->SetLineColor(kPink-9);
0261       } else {
0262         gSyst->SetFillColorAlpha(kAzure-4, 0.35); // alpha controls transparency
0263         gSyst->SetLineColor(kAzure-4);
0264       }
0265       gSyst->SetMarkerStyle(0);
0266       gSyst->Draw("P E2 SAME");
0267
0268       gPad->Modified();
0269       gPad->Update();
0270       double min_x = gPad->GetUxmin();
0271       double max_x = gPad->GetUxmax();
0272
0273       // Add "sPHENIX internal"
0274       TPad *p = new TPad("p","p",0.,0.,1.,1.); p->SetFillStyle(0); p->Draw(); p->cd();
0275       TBox *whiteBox = new TBox(0.185, 0.72, 0.49, 0.90);
0276       whiteBox->Draw();
0277       canvas->cd();
0278       whiteBox->SetFillColorAlpha(kWhite, 1);
0279       std::stringstream stream;
0280       stream.str("");
0281       TLatex latex;
0282       latex.SetNDC();
0283       latex.SetTextColor(kBlack);
0284       latex.DrawLatex(0.22, 0.85, "#font[72]{sPHENIX} Internal");
0285       latex.DrawLatex(0.22, 0.79, "p^{#uparrow}+p #sqrt{s} = 200 GeV");
0286       latex.SetTextSize(0.03);
0287       latex.DrawLatex(0.19, 0.73, "7% polarization scale uncertainty not shown");
0288       latex.SetTextSize(0.07);
0289
0290       if (iP == 0) {
0291         latex.DrawLatex(0.57, 0.84, "p^{#uparrow}+p #rightarrow #pi^{0} X");
0292       } else {
0293         latex.DrawLatex(0.57, 0.84, "p^{#uparrow}+p #rightarrow #eta X");
0294       }
0295
0296       TLine *tline = new TLine();
0297       tline->SetLineWidth(2);
0298       tline->SetLineColor(kBlack);
0299       tline->SetLineStyle(kDashed);
0300       tline->DrawLine(min_x, 0, max_x, 0);
0301
0302       canvas->Update();
0303       canvas->Draw();
0304       canvas->SaveAs((plots_folder_eta +  "/" + canvas_name.str() + ".png").c_str());
0305       canvas->SaveAs((plots_folder_eta +  "/" + canvas_name.str() + ".pdf").c_str());
0306     }
0307     {
0308       std::stringstream inputfilename;
0309       inputfilename << inputfolder << "asym_summary_" << particles[iP] << "_xf" << ".txt";
0310       std::stringstream canvas_name;
0311       canvas_name << "canvas_asym_" << particles[iP] << "_xf";
0312       std::vector<double> x, y, estat, esyst, exsyst;
0313
0314       parse_file(inputfilename.str(), x, y, estat, esyst, exsyst, 0.002);
0315
0316       const int n = x.size();
0317       auto gSyst = new TGraphErrors(n);
0318       auto gStat = new TGraphErrors(n);
0319
0320       for (int i = 0; i < n; ++i) {
0321         gSyst->SetPoint(i, x[i], y[i]);
0322         gSyst->SetPointError(i, exsyst[i], esyst[i]);
0323
0324         gStat->SetPoint(i, x[i], y[i]);
0325         gStat->SetPointError(i, 0.0, estat[i]);
0326       }
0327       gStat->Print();
0328       gSyst->Print();
0329       SetsPhenixStyle();
0330
0331       TCanvas *canvas = new TCanvas(canvas_name.str().c_str(), "", 1600, 900);
0332       gStat->SetTitle(";x_{F};A_{N}");
0333
0334       double y_bound = 0;
0335       y_bound = GetAbsMaxInRange(gStat, -0.3, 0.3);
0336
0337       int color = (iP == 0 ? kRed : kBlue);
0338
0339       gStat->SetLineWidth(2);
0340       gStat->SetLineColor(color);
0341       gStat->SetMarkerColor(color);
0342       gStat->SetMarkerStyle(kFullSquare);
0343       gStat->SetMarkerSize(2.4);
0344       // gStat->SetMinimum(-2 * y_bound);
0345       // gStat->SetMaximum(2 * y_bound);
0346       gStat->SetMinimum(xFBounds[iP][0]);
0347       gStat->SetMaximum(xFBounds[iP][1]);
0348       gStat->Draw("AP E1");
0349
0350        // Style systematic boxes (transparent fill)
0351       if (iP == 0) {
0352         gSyst->SetFillColorAlpha(kPink-9, 0.35); // alpha controls transparency
0353         gSyst->SetLineColor(kPink-9);
0354       } else {
0355         gSyst->SetFillColorAlpha(kAzure-4, 0.35); // alpha controls transparency
0356         gSyst->SetLineColor(kAzure-4);
0357       }
0358       gSyst->SetMarkerStyle(0);
0359       gSyst->Draw("P E2 SAME");
0360
0361       gPad->Modified();
0362       gPad->Update();
0363       double min_x = gPad->GetUxmin();
0364       double max_x = gPad->GetUxmax();
0365
0366       // Add "sPHENIX internal"
0367       TPad *p = new TPad("p","p",0.,0.,1.,1.); p->SetFillStyle(0); p->Draw(); p->cd();
0368       TBox *whiteBox = new TBox(0.185, 0.72, 0.49, 0.90);
0369       whiteBox->Draw();
0370       canvas->cd();
0371       whiteBox->SetFillColorAlpha(kWhite, 1);
0372       std::stringstream stream;
0373       stream.str("");
0374       TLatex latex;
0375       latex.SetNDC();
0376       latex.SetTextColor(kBlack);
0377       latex.DrawLatex(0.22, 0.85, "#font[72]{sPHENIX} Internal");
0378       latex.DrawLatex(0.22, 0.79, "p^{#uparrow}+p #sqrt{s} = 200 GeV");
0379       latex.SetTextSize(0.03);
0380       latex.DrawLatex(0.19, 0.73, "7% polarization scale uncertainty not shown");
0381       latex.SetTextSize(0.07);
0382
0383       if (iP == 0) {
0384         latex.DrawLatex(0.57, 0.84, "p^{#uparrow}+p #rightarrow #pi^{0} X");
0385       } else {
0386         latex.DrawLatex(0.57, 0.84, "p^{#uparrow}+p #rightarrow #eta X");
0387       }
0388
0389       TLine *tline = new TLine();
0390       tline->SetLineWidth(2);
0391       tline->SetLineColor(kBlack);
0392       tline->SetLineStyle(kDashed);
0393       tline->DrawLine(min_x, 0, max_x, 0);
0394
0395       canvas->Update();
0396       canvas->Draw();
0397       canvas->SaveAs((plots_folder_xf +  "/" + canvas_name.str() + ".png").c_str());
0398       canvas->SaveAs((plots_folder_xf +  "/" + canvas_name.str() + ".pdf").c_str());
0399     }
0400   }
0401   gSystem->Exit(0);
0402 }