Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
makeV1_BUP2020.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file makeV1_BUP2020.C
1 #include <fstream>
2 #include <iostream>
3 #include "stdio.h"
4 
5 #include <TChain.h>
6 #include <TFile.h>
7 #include <TGraphAsymmErrors.h>
8 #include <TGraphErrors.h>
9 #include <TLatex.h>
10 #include <TLegend.h>
11 #include <TLine.h>
12 #include <TROOT.h>
13 #include <TMath.h>
14 #include <TString.h>
15 #include <TTree.h>
16 #include <TVectorD.h>
17 #include <TVirtualFitter.h>
18 #include <cassert>
19 #include <cmath>
20 #include <vector>
21 
22 #include "SaveCanvas.C"
23 #include "sPhenixStyle.C"
24 
25 using namespace std;
26 
27 Double_t v1_err(Double_t sig, Double_t v1, Double_t Res)
28 {
29  const Double_t Pi = 3.1415927;
30  return 2.3 / sig * sqrt(1 - 16 * v1 * v1 / Pi / Pi) / Res; // 2.3 - coefficient from simple two rapidity bin linear fit
31 }
32 
33 TGraphErrors *GraphShiftScaling(TGraphErrors *gr_src, const double x_shift, const double err_scaling)
34 {
35  assert(gr_src);
36 
37  const int npoint = gr_src->GetN();
38 
39  TVectorD vx(npoint);
40  TVectorD vy(npoint);
41  TVectorD vex(npoint);
42  TVectorD vey(npoint);
43 
44  int nfilled = 0;
45  for (int i = 0; i < npoint; ++i)
46  {
47  const double &x = gr_src->GetX()[i];
48  // if (x<x_min or x>x_max) continue;
49 
50  vx[nfilled] = x + x_shift;
51  vy[nfilled] = gr_src->GetY()[i];
52  vex[nfilled] = gr_src->GetEX()[i];
53  vey[nfilled] = gr_src->GetEY()[i] * err_scaling;
54 
55  ++nfilled;
56  }
57 
58  TGraphErrors *gr = new TGraphErrors(nfilled, vx.GetMatrixArray(), vy.GetMatrixArray(),
59  vex.GetMatrixArray(), vey.GetMatrixArray());
60 
61  gr->SetMarkerColor(gr_src->GetMarkerColor());
62  gr->SetMarkerStyle(gr_src->GetMarkerStyle());
63  gr->SetMarkerSize(gr_src->GetMarkerSize());
64  gr->SetLineWidth(gr_src->GetLineWidth());
65  gr->SetLineColor(gr_src->GetLineColor());
66 
67  return gr;
68 }
69 
70 void makeV1_BUP2020(const Bool_t sPHENIX = 1)
71 {
72  // gROOT->LoadMacro("sPhenixStyle.C");
73  SetsPhenixStyle();
74 
75  const Double_t MassD = 1.86484;
76 
77  // STAR data points
78  const Int_t NS = 4;
79  const Double_t offset = 0.02;
80  Double_t yp[NS], v1p[NS], v1pe[NS], v1pes[NS];
81  Double_t yn[NS], v1n[NS], v1ne[NS], v1nes[NS];
82  ifstream inData;
83  inData.open("D0v1_STAR_final.txt");
84  for (int i = 0; i < NS; i++)
85  {
86  inData >> yp[i] >> v1p[i] >> v1pe[i] >> v1pes[i];
87  yp[i] -= offset;
88  }
89  inData.close();
90  inData.open("D0barv1_STAR_final.txt");
91  for (int i = 0; i < NS; i++)
92  {
93  inData >> yn[i] >> v1n[i] >> v1ne[i] >> v1nes[i];
94  yn[i] += offset;
95  }
96  inData.close();
97 
98  // TGraphErrors *gr_D_STAR = new TGraphErrors("D0v1_STAR_final.txt","%lg %lg %lg");
99  TGraphErrors *gr_D_STAR = new TGraphErrors(NS, yp, v1p, 0, v1pe);
100  gr_D_STAR->SetMarkerColor(kRed - 6);
101  gr_D_STAR->SetMarkerStyle(kFullStar);
102  gr_D_STAR->SetMarkerSize(2);
103  gr_D_STAR->SetLineWidth(4);
104  gr_D_STAR->SetLineColor(kRed - 6);
105 
106  // TGraphErrors *gr_Dbar_STAR = new TGraphErrors("D0barv1_STAR_final.txt","%lg %lg %lg");
107  TGraphErrors *gr_Dbar_STAR = new TGraphErrors(NS, yn, v1n, 0, v1ne);
108  gr_Dbar_STAR->SetMarkerColor(kRed - 6);
109  gr_Dbar_STAR->SetMarkerStyle(kOpenStar);
110  gr_Dbar_STAR->SetMarkerSize(2);
111  gr_Dbar_STAR->SetLineWidth(4);
112  gr_Dbar_STAR->SetLineColor(kRed - 6);
113 
114  TFile *fin = new TFile("D0_significance.root");
115  TGraph *gr_sig_D = (TGraph *) fin->Get("gProD0_0_80_noPid");
116  double sig2_tot = 0.;
117  for (int i = 0; i < gr_sig_D->GetN(); i++)
118  {
119  sig2_tot += TMath::Power(gr_sig_D->GetY()[i], 2.0);
120  }
121  double sig_tot = sqrt(sig2_tot);
122  cout << " Total D/Dbar significance in |y|<1 = " << sig_tot << endl;
123  double sig_tot_D = sig_tot / sqrt(2.); // one charge sign
124 
125  const Double_t EPRes = 0.37; // STAR 1st EP resolution
126 
127  // theory curves
128  // Greco - D/Dbar v1 difference
129  TGraph *gr_D_Greco = new TGraph("D0v1_RHIC_Greco.txt", "%lg %lg");
130  TGraph *gr_Dbar_Greco = new TGraph("D0barv1_RHIC_Greco.txt", "%lg %lg"); // be careful, this calculation has a sign flip
131  TGraph *gr_DDbar_Bozek_tmp = new TGraph("DDbarv1_eta_1712.01180.txt", "%lg %lg"); // v1 vs eta, in percentage
132  double x[100], y[100];
133  for (int i = 0; i < gr_DDbar_Bozek_tmp->GetN(); i++)
134  {
135  double eta = gr_DDbar_Bozek_tmp->GetX()[i];
136  x[i] = eta; // questionable!!
137  y[i] = gr_DDbar_Bozek_tmp->GetY()[i] / 100.;
138  }
139  TGraph *gr_DDbar_Bozek = new TGraph(gr_DDbar_Bozek_tmp->GetN(), x, y);
140  cout << " AAAA " << endl;
141 
142  // for projections v1
143  const Int_t N_D = 8;
144  double y_D[N_D], v1_D[N_D], v1_Dbar[N_D], v1_err_D[N_D], v1_err_Dbar[N_D];
145 
146  for (int i = 0; i < N_D; i++)
147  {
148  y_D[i] = (i + 0.5) * 2.0 / N_D - 1.0;
149  v1_D[i] = -1.0 * gr_D_Greco->Eval(y_D[i]) + gr_DDbar_Bozek->Eval(y_D[i]); // there is a sign flip in Greco's calculation
150  v1_Dbar[i] = -1.0 * gr_Dbar_Greco->Eval(y_D[i]) + gr_DDbar_Bozek->Eval(y_D[i]);
151 
152  double sig_per_bin = sig_tot_D / sqrt((double) N_D); // //not combine symmetric rapidity bins
153  v1_err_D[i] = v1_err(sig_per_bin, v1_D[i], EPRes);
154  v1_err_Dbar[i] = v1_err(sig_per_bin, v1_Dbar[i], EPRes);
155  cout << y_D[i] << " " << v1_D[i] << "+/-" << v1_err_D[i] << endl;
156  }
157 
158  TGraph *gr_D = new TGraph(N_D, y_D, v1_D);
159  gr_D->SetMarkerSize(2);
160  gr_D->SetMarkerColor(1);
161  gr_D->SetMarkerStyle(20);
162  gr_D->SetLineWidth(2.);
163  gr_D->SetLineStyle(1);
164  gr_D->SetLineColor(1);
165 
166  TGraph *gr_Dbar = new TGraph(N_D, y_D, v1_Dbar);
167  gr_Dbar->SetMarkerSize(2);
168  gr_Dbar->SetMarkerColor(4);
169  gr_Dbar->SetMarkerStyle(24);
170  gr_Dbar->SetLineWidth(2.);
171  gr_Dbar->SetLineStyle(1);
172  gr_Dbar->SetLineColor(4);
173 
174  TGraphErrors *gr_D_proj = new TGraphErrors(N_D / 2, y_D + 4, v1_D + 4, 0, v1_err_D + 4);
175  gr_D_proj->SetMarkerSize(2);
176  gr_D_proj->SetMarkerColor(1);
177  gr_D_proj->SetMarkerStyle(20);
178  gr_D_proj->SetLineWidth(4.);
179  gr_D_proj->SetLineStyle(1);
180  gr_D_proj->SetLineColor(1);
181 
182  TGraphErrors *gr_D_proj_2 = new TGraphErrors(N_D / 2, y_D, v1_D, 0, v1_err_D);
183  gr_D_proj_2->SetMarkerSize(2);
184  gr_D_proj_2->SetMarkerColor(1);
185  gr_D_proj_2->SetMarkerStyle(20);
186  gr_D_proj_2->SetLineWidth(4.);
187  gr_D_proj_2->SetLineStyle(1);
188  gr_D_proj_2->SetLineColor(1);
189 
190  TGraphErrors *gr_Dbar_proj = new TGraphErrors(N_D / 2, y_D + 4, v1_Dbar + 4, 0, v1_err_Dbar + 4);
191  gr_Dbar_proj->SetMarkerSize(2);
192  gr_Dbar_proj->SetMarkerColor(4);
193  gr_Dbar_proj->SetMarkerStyle(20);
194  gr_Dbar_proj->SetLineWidth(4.);
195  gr_Dbar_proj->SetLineStyle(1);
196  gr_Dbar_proj->SetLineColor(4);
197 
198  TGraphErrors *gr_Dbar_proj_2 = new TGraphErrors(N_D / 2, y_D, v1_Dbar, 0, v1_err_Dbar);
199  gr_Dbar_proj_2->SetMarkerSize(2);
200  gr_Dbar_proj_2->SetMarkerColor(4);
201  gr_Dbar_proj_2->SetMarkerStyle(20);
202  gr_Dbar_proj_2->SetLineWidth(4.);
203  gr_Dbar_proj_2->SetLineStyle(1);
204  gr_Dbar_proj_2->SetLineColor(4);
205 
207  // BUP2020 lumi scaling
209 
210  const double refAuAuMB = 240e9;
211  const double refAuAuXSec = 6.8252; // b
212 
213  const double AuAu_Ncoll_C0_10 = 960.2; // [DOI:?10.1103/PhysRevC.87.034911?]
214  const double AuAu_Ncoll_C0_20 = 770.6; // [DOI:?10.1103/PhysRevC.91.064904?]
215  const double AuAu_Ncoll_60_70 = 29.8; //PHYSICAL REVIEW C 87, 034911 (2013)
216  const double AuAu_Ncoll_70_80 = 12.6; //PHYSICAL REVIEW C 87, 034911 (2013)
217  const double AuAu_Ncoll_C0_100 = 238.5; // BUP2020
218  const double pAu_Ncoll_C0_100 = 4.7; // pb^-1 [sPH-TRG-000]
219 
220  const double AuAu_rec_3year = (5.7 + 15) * 1e9; // BUP2020
221  const double AuAu_rec_5year = AuAu_rec_3year + 30e9; // BUP2020
222 
223  const double error_scale_3yr = sqrt((refAuAuMB) / (AuAu_rec_3year * refAuAuXSec));
224  const double error_scale_5yr = sqrt((refAuAuMB) / (AuAu_rec_5year * refAuAuXSec));
225 
226  TGraphErrors *gr_D_proj_3yr = GraphShiftScaling(gr_D_proj, 0.05, error_scale_3yr);
227  TGraphErrors *gr_D_proj_2_3yr = GraphShiftScaling(gr_D_proj_2, 0.05, error_scale_3yr);
228  TGraphErrors *gr_Dbar_proj_3yr = GraphShiftScaling(gr_Dbar_proj, 0.05, error_scale_3yr);
229  TGraphErrors *gr_Dbar_proj_2_3yr = GraphShiftScaling(gr_Dbar_proj_2, 0.05, error_scale_3yr);
230 
231  TGraphErrors *gr_D_proj_5yr = GraphShiftScaling(gr_D_proj, 0., error_scale_5yr);
232  TGraphErrors *gr_D_proj_2_5yr = GraphShiftScaling(gr_D_proj_2, 0., error_scale_5yr);
233  TGraphErrors *gr_Dbar_proj_5yr = GraphShiftScaling(gr_Dbar_proj, 0., error_scale_5yr);
234  TGraphErrors *gr_Dbar_proj_2_5yr = GraphShiftScaling(gr_Dbar_proj_2, 0., error_scale_5yr);
235 
236  gr_D_proj_3yr->SetMarkerStyle(kOpenCircle);
237  gr_D_proj_2_3yr->SetMarkerStyle(kOpenCircle);
238  gr_Dbar_proj_3yr->SetMarkerStyle(kOpenCircle);
239  gr_Dbar_proj_2_3yr->SetMarkerStyle(kOpenCircle);
240 
241  gr_D_proj_3yr->SetMarkerColor(kGray + 2);
242  gr_D_proj_2_3yr->SetMarkerColor(kGray + 2);
243  gr_Dbar_proj_3yr->SetMarkerColor(kBlue - 6);
244  gr_Dbar_proj_2_3yr->SetMarkerColor(kBlue - 6);
245 
246  gr_D_proj_3yr->SetLineColor(kGray + 2);
247  gr_D_proj_2_3yr->SetLineColor(kGray + 2);
248  gr_Dbar_proj_3yr->SetLineColor(kBlue - 6);
249  gr_Dbar_proj_2_3yr->SetLineColor(kBlue - 6);
250 
251  // TCanvas *c1 = new TCanvas("v1_5y_compare", "v1_5y_compare", 0, 0, 800, 600);
252 
253  TCanvas *c1 = new TCanvas("D0_BUP2020_AuAu_v1_5y_compare", "D0_BUP2020_AuAu_v1_5y_compare", 1100, 800);
254  c1->Divide(1, 1);
255  int idx = 1;
256  TPad *p;
257 
258  p = (TPad *) c1->cd(idx++);
259  c1->Update();
260  // gStyle->SetOptFit(0);
261  // gStyle->SetOptStat(0);
262  // gStyle->SetEndErrorSize(0.01);
263  // c1->SetFillColor(10);
264  // c1->SetBorderMode(0);
265  // c1->SetBorderSize(2);
266  // c1->SetFrameFillColor(0);
267  // c1->SetFrameBorderMode(0);
268 
269  // // c1->SetGridx();
270  // // c1->SetGridy();
271  // c1->SetTickx();
272  // c1->SetTicky();
273 
274  // // c1->SetLogx();
275  // // c1->SetLogy();
276 
277  // c1->SetLeftMargin(0.13);
278  // c1->SetBottomMargin(0.16);
279  // c1->SetTopMargin(0.02);
280  // c1->SetRightMargin(0.06);
281 
282  double x1 = -1.2;
283  double x2 = 1.2;
284  double y1 = -0.1;
285  double y2 = 0.1;
286 
287  TH1D *d0 = new TH1D("d0", "", 1, x1, x2);
288  d0->SetMinimum(y1);
289  d0->SetMaximum(y2);
290  // d0->GetXaxis()->SetNdivisions(208);
291  // d0->GetXaxis()->CenterTitle();
292  d0->GetXaxis()->SetTitle("Rapidity");
293  // d0->GetXaxis()->SetTitleOffset(1.2);
294  // d0->GetXaxis()->SetTitleSize(0.06);
295  // d0->GetXaxis()->SetLabelOffset(0.01);
296  // d0->GetXaxis()->SetLabelSize(0.045);
297  // d0->GetXaxis()->SetLabelFont(42);
298  // d0->GetXaxis()->SetTitleFont(42);
299  // d0->GetYaxis()->SetNdivisions(505);
300  // d0->GetYaxis()->CenterTitle();
301  d0->GetYaxis()->SetTitle("v_{1}");
302  // d0->GetYaxis()->SetTitleOffset(1.0);
303  // d0->GetYaxis()->SetTitleSize(0.06);
304  // d0->GetYaxis()->SetLabelOffset(0.005);
305  // d0->GetYaxis()->SetLabelSize(0.045);
306  // d0->GetYaxis()->SetLabelFont(42);
307  // d0->GetYaxis()->SetTitleFont(42);
308  // d0->SetLineWidth(2);
309  d0->Draw("c");
310 
311  TLine *l0 = new TLine(x1, 0, x2, 0);
312  l0->SetLineWidth(1);
313  l0->SetLineColor(kBlack);
314  l0->SetLineStyle(1);
315  l0->Draw("same");
316 
317  gr_D_STAR->Draw("p");
318  gr_Dbar_STAR->Draw("p");
319 
320  if (sPHENIX)
321  {
322  gr_D->Draw("c");
323  gr_Dbar->Draw("c");
324 
325  gr_D_proj_2_3yr->Draw("p");
326  gr_Dbar_proj_2_3yr->Draw("p");
327 
328  gr_D_proj_3yr->Draw("p");
329  gr_Dbar_proj_3yr->Draw("p");
330 
331  gr_D_proj_2_5yr->Draw("p");
332  gr_Dbar_proj_2_5yr->Draw("p");
333 
334  gr_D_proj_5yr->Draw("p");
335  gr_Dbar_proj_5yr->Draw("p");
336  }
337 
338  TLegend *leg = new TLegend(.3, .77, .9, .9);
339  leg->SetFillStyle(0);
340  leg->AddEntry("", "#it{#bf{sPHENIX}} Projection, Years 1-5", "");
341  leg->AddEntry("", ("0-80% Au+Au, Res(#Psi_{1})=0.37"), "");
342  // leg->AddEntry("", Form("%.0f nb^{-1} str. O+O", OO_rec_5year/1e9), "");
343  // leg->AddEntry("", Form("%.0f nb^{-1} str. Ar+Ar", ArAr_rec_5year/1e9), "");
344  leg->Draw();
345 
346  leg = new TLegend(.2, .18, .4, .4, "D^{0}");
347  leg->SetFillStyle(0);
348  leg->AddEntry(gr_D_proj_3yr, " ", "pe");
349  leg->AddEntry(gr_D_proj_5yr, " ", "pe");
350  leg->AddEntry(gr_D_STAR, " ", "pe");
351  leg->Draw();
352 
353  leg = new TLegend(.25, .18, .5, .4, " #bar{D}^{0}");
354  leg->SetFillStyle(0);
355  leg->AddEntry(gr_Dbar_proj_3yr, Form("sPHENIX, Year 1-3, %.0fnb^{-1} rec. Au+Au", AuAu_rec_3year / 1e9), "pe");
356  leg->AddEntry(gr_Dbar_proj_5yr, Form("sPHENIX, Year 1-5, %.0fnb^{-1} rec.+str. Au+Au", AuAu_rec_5year / 1e9), "pe");
357  leg->AddEntry(gr_Dbar_STAR, "STAR, 10-80% Au+Au, PRL#bf{123}, 162301", "pe");
358  leg->Draw();
359 
360  // TLegend *leg = new TLegend(0.12, 0.2, 0.52, 0.4);
361  // leg->SetFillColor(10);
362  // leg->SetFillStyle(10);
363  // leg->SetLineStyle(4000);
364  // leg->SetLineColor(10);
365  // leg->SetLineWidth(0.);
366  // leg->SetTextFont(42);
367  // leg->SetTextSize(0.045);
368  // leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation", "");
369  // leg->AddEntry("", "Au+Au #sqrt{s_{NN}}=200 GeV", "");
370  // leg->AddEntry("", "240B MB, EPRes=0.37", "");
371  // if (sPHENIX) leg->Draw();
372  //
373  // // TLegend *leg;
374  // if (sPHENIX)
375  // leg = new TLegend(0.62, 0.66, 0.98, 0.92);
376  // else
377  // leg = new TLegend(0.62, 0.79, 0.98, 0.92);
378  // leg->SetFillColor(10);
379  // leg->SetFillStyle(10);
380  // leg->SetLineStyle(4000);
381  // leg->SetLineColor(10);
382  // leg->SetLineWidth(0.);
383  // leg->SetTextFont(42);
384  // leg->SetTextSize(0.04);
385  // leg->AddEntry(gr_D_STAR, "D^{0} - STAR", "pl");
386  // leg->AddEntry(gr_Dbar_STAR, "#bar{D}^{0} - STAR", "pl");
387  // if (sPHENIX)
388  // {
389  // leg->AddEntry(gr_D_proj, "D^{0} - sPHENIX proj.", "pl");
390  // leg->AddEntry(gr_Dbar_proj, "#bar{D}^{0} - sPHENIX proj.", "pl");
391  // }
392  // // leg->AddEntry(gr_proj_D,"Uncert. of D^{0}","p");
393  // // leg->AddEntry(gr_proj_D_B,"Uncert. of D^{0} from B","p");
394  // leg->Draw();
395 
396  // c1->SaveAs(Form("fig/v1_proj_240B_%d.pdf", sPHENIX));
397  // c1->SaveAs(Form("fig/v1_proj_240B_%d.png", sPHENIX));
398 
399  SaveCanvas(c1, "fig_BUP2020/" + TString(c1->GetName()), kTRUE);
400 }