Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
upsilon_raa_sphenix_estimator.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file upsilon_raa_sphenix_estimator.C
1 
2 
3 void upsilon_raa_sphenix_estimator() {
4 
5  //gROOT->LoadMacro("sPHENIXStyle/sPhenixStyle.C");
6  //SetsPhenixStyle();
7 
8 gStyle->SetOptStat(0);
9 gStyle->SetOptTitle(0);
10 gStyle->SetOptFit(0);
11 
12  bool yrs13 = true;
13  bool yrs15 = true;
14  bool ups3s = true;
15 
16  double offset = 0.0;
17  if(yrs13 && yrs15)
18  offset = 3;
19 
20  // Input assumptions about level of suppression are taken from the theory paper:
21  // Michael Strickland, Dennis Bazow, Nucl.Phys. A879 (2012) 25-58
22  //==============================================
23  // Read in Strickland curves for different eta/s values
24  // These are exclusive - i.e. they are the
25  // suppression of that state with no feed-down
26  // effects included
27  //==============================================
28  // The curves vs eta/s read in below were provided by Strickland privately, since
29  // the paper contains only the eta/s = 1 cases for RHIC energy
30 
31  // Get the exclusive raa values for the Y1S, Y2S, Y3S, chib1, chib2
32  // Checked that this is being read in correctly
33  double str_npart[101];
34  double str_raa[5][3][101];
35  for(int istate=0;istate<5;istate++)
36  {
37  for(int ietas=0;ietas<3;ietas++)
38  {
39  if(istate < 3)
40  {
41  char fname[500];
42  sprintf(fname,"./strickland_calculations/Y%is-potb-eta%i-npart.dat",istate+1,ietas+1);
43  ifstream fin(fname);
44 
45  for(int inpart=0;inpart<101;inpart++)
46  {
47  fin >> str_npart[inpart] >> str_raa[istate][ietas][inpart];
48  //cout << istate << " " << ietas << " " << inpart << " " << str_npart[inpart] << " " << str_raa[istate][ietas][inpart] << endl;
49  }
50  fin.close();
51  }
52 
53  if(istate > 2)
54  {
55  char fname[500];
56  sprintf(fname,"./strickland_calculations/chib%i-potb-eta%i-npart.dat",istate-2,ietas+1);
57  ifstream fin(fname);
58 
59  for(int inpart=0;inpart<101;inpart++)
60  {
61  fin >> str_npart[inpart] >> str_raa[istate][ietas][inpart];
62  //cout << istate << " " << ietas << " " << inpart << " " << str_npart[inpart] << " " << str_raa[istate][ietas][inpart] << endl;
63  }
64  fin.close();
65  }
66  }
67 
68  }
69 
70  // Now construct the inclusive RAA values from Strickland's exclusive modifications
71  // ff's from table II of the paper for the 1S state are as follows
72  // They add up to 1.0
73  double ff1S[5] = {0.51, 0.107, 0.008, 0.27, 0.105}; // Y1S, Y2S, Y3S, chib1, chib2 respectively
74  // These are from arXiv:1210.7512
75  double ff2S[2] = {0.5, 0.5}; // Y2S and Y3S respectively
76 
77  double str_raa_inclusive[3][3][101];
78 
79  // checked the Y1S inclusive result against figure 10 of arXiv:1112.2761
80  // There is no plot available for the 2S and 3S
81  for(int ietas=0;ietas<3;ietas++)
82  for(int inpart=0;inpart<101;inpart++)
83  {
84  str_raa_inclusive[0][ietas][inpart] =
85  str_raa[0][ietas][inpart] * ff1S[0]
86  + str_raa[1][ietas][inpart] * ff1S[1]
87  + str_raa[2][ietas][inpart] * ff1S[2]
88  + str_raa[3][ietas][inpart] * ff1S[3]
89  + str_raa[4][ietas][inpart] * ff1S[4];
90 
91  str_raa_inclusive[1][ietas][inpart] =
92  str_raa[1][ietas][inpart] * ff2S[0]
93  + str_raa[2][ietas][inpart] * ff2S[1];
94 
95  str_raa_inclusive[2][ietas][inpart] = str_raa[2][ietas][inpart];
96  }
97 
98  TGraph *grth[3][3];
99  for(int is = 0;is<3;is++)
100  {
101  for(int ieta =0; ieta < 3; ++ieta)
102  {
103  grth[is][ieta]=new TGraph(101,str_npart,str_raa_inclusive[is][ieta]);
104  }
105  }
106 
107 
108  static const int NCENT = 7;
109  // static const int NCENT3S = 4;
110  static const int NCENT3S = 3;
111 
112  double Ncoll[NCENT] = {1067, 858, 610, 378, 224, 94.2, 14.5};
113  double Npart[NCENT] = {351, 302, 236, 168, 116, 61.6, 14.4};
114  double Npart_plot[2][NCENT];
115  for(int i=0;i<NCENT;++i)
116  {
117  Npart_plot[0][i] = Npart[i] + offset;
118  Npart_plot[1][i] = Npart[i] - offset;
119  }
120 
121  // separate binning for Y(3S)
122  //double Npart3S[NCENT3S] = {281, 142, 62, 14.4};
123  double Npart3S[NCENT3S] = {243, 80.0, 14.4};
124  double Npart3S_plot[2][NCENT3S];
125  for(int i=0;i<NCENT3S;++i)
126  {
127  Npart3S_plot[0][i] = Npart3S[i] + offset;
128  Npart3S_plot[1][i] = Npart3S[i] - offset;
129  }
130 
131 
132  double raa1[2][NCENT];
133  double raa2[2][NCENT];
134  for(int iset=0;iset<2;++iset)
135  for(int i=0; i<NCENT; ++i)
136  {
137  raa1[iset][i] = grth[0][1]->Eval(Npart_plot[iset][i]);
138  raa2[iset][i] = grth[1][1]->Eval(Npart_plot[iset][i]);
139  cout << " icent " << i << " Npart " << Npart_plot[iset][i] << " raa1 " << raa1[iset][i] << " raa2 " << raa2[iset][i] << endl;
140  }
141 
142  double raa3[2][NCENT3S];
143  for(int iset=0;iset<2;++iset)
144  for(int i=0; i<NCENT3S; ++i)
145  {
146  raa3[iset][i] = grth[2][1]->Eval(Npart3S_plot[iset][i]);
147  cout << " 3S: icent " << i << " Npart " << Npart3S_plot[iset][i] << " raa1 " << " raa3 " << raa3[iset][i] << endl;
148  }
149 
150  // output the RAA values to feed back into the previous macro
151  cout << " double raa1[" << NCENT << "] = {";
152  for(int i=0; i< NCENT; ++i)
153  {
154  if(i != NCENT-1)
155  cout << raa1[0][i] << ", ";
156  else
157  cout << raa1[0][i] << "}; ";
158  }
159  cout << endl;
160 
161  // output the RAA values to feed back into the previous macro
162  cout << " double raa2[" << NCENT << "] = {";
163  for(int i=0; i< NCENT; ++i)
164  {
165  if(i != NCENT-1)
166  cout << raa2[0][i] << ", ";
167  else
168  cout << raa2[0][i] << "}; ";
169  }
170  cout << endl;
171 
172  // output the RAA values to feed back into the previous macro
173  cout << " double raa3[" << NCENT3S << "] = {";
174  for(int i=0; i< NCENT3S; ++i)
175  {
176  if(i != NCENT3S-1)
177  cout << raa3[0][i] << ", ";
178  else
179  cout << raa3[0][i] << "}; ";
180  }
181  cout << endl;
182 
183  // These errors are the only thing that change with different luminosity
184  // updated 8/26/20
185  double yrs13_erraa1[NCENT] = {0.021, 0.0227, 0.0211, 0.0251, 0.0315, 0.0365, 0.0665};
186  double yrs13_erraa2[NCENT] = {0.0327, 0.0424, 0.0335, 0.0401, 0.0537, 0.0662, 0.157};
187  double yrs13_erraa3[NCENT3S] = {0.042 ,0.094, 0.346};
188 
189  //double yrs15_erraa1[NCENT] = {0.012, 0.013, 0.012, 0.015, 0.019, 0.022, 0.042};
190  double yrs15_erraa1[NCENT] = {0.0136, 0.0147, 0.0137, 0.0163, 0.0203, 0.0235, 0.0435};
191  //double yrs15_erraa2[NCENT] = {0.022, 0.024, 0.021, 0.025, 0.034, 0.041, 0.097};
192  double yrs15_erraa2[NCENT] = {0.0235, 0.0233, 0.0207, 0.0277, 0.0364, 0.0424, 0.103};
193  double yrs15_erraa3[NCENT3S] = {0.034, 0.0547, 0.208};
194 
195  double erraa1[2][NCENT];
196  double erraa2[2][NCENT];
197  for(int i=0; i<NCENT;++i)
198  {
199  erraa1[0][i] = yrs13_erraa1[i];
200  erraa2[0][i] = yrs13_erraa2[i];
201  erraa1[1][i] = yrs15_erraa1[i];
202  erraa2[1][i] = yrs15_erraa2[i];
203  }
204  double erraa3[2][NCENT3S];
205  for(int i=0; i<NCENT3S;++i)
206  {
207  erraa3[0][i] = yrs13_erraa3[i];
208  erraa3[1][i] = yrs15_erraa3[i];
209  }
210 
211 
212  TCanvas *craa = new TCanvas("craa","craa",0,0,1000,800);
213  TH1D *h = new TH1D("h","h",100,-10,400);
214  h->GetXaxis()->SetTitle("N_{part}");
215  h->GetYaxis()->SetTitle("R_{AA}");
216  h->SetMaximum(1.2);
217  h->SetMinimum(0);
218 
219  h->Draw();
220 
221  TGraphErrors *gr1[2];
222  TGraphErrors *gr2[2];
223  TGraphErrors *gr3[2];
224  for(int i=0; i< 2;i++)
225  {
226  gr1[i] = new TGraphErrors(NCENT,Npart_plot[i], raa1[i], 0, erraa1[i]);
227  gr1[i]->SetMarkerStyle(20);
228  gr1[i]->SetMarkerSize(1.5);
229 
230  gr2[i] = new TGraphErrors(NCENT,Npart_plot[i], raa2[i], 0, erraa2[i]);
231  gr2[i]->SetMarkerStyle(21);
232  gr2[i]->SetMarkerSize(1.5);
233  gr2[i]->SetMarkerColor(kRed);
234  gr2[i]->SetLineColor(kRed);
235 
236  gr3[i] = new TGraphErrors(NCENT3S,Npart3S_plot[i], raa3[i], 0, erraa3[i]);
237  gr3[i]->SetMarkerStyle(20);
238  gr3[i]->SetMarkerSize(1.5);
239  gr3[i]->SetMarkerColor(kBlue);
240  gr3[i]->SetLineColor(kBlue);
241  }
242 
243  if(yrs13)
244  {
245  gr1[0]->Draw("p");
246  gr2[0]->Draw("p");
247  if(ups3s) gr3[0]->Draw("p");
248  }
249 
250  if(yrs15)
251  {
252  if(yrs13 && yrs15)
253  {
254  gr1[0]->SetMarkerStyle(24);
255  gr2[0]->SetMarkerStyle(25);
256  gr3[0]->SetMarkerStyle(25);
257  }
258  gr1[1]->Draw("p");
259  gr2[1]->Draw("p");
260  if(ups3s) gr3[1]->Draw("p");
261  }
262 
263 
264  TLegend *lups = new TLegend(0.38, 0.72, 0.72, 0.87,"#bf{#it{sPHENIX}} Projection");
265  lups->SetBorderSize(0);
266  lups->SetFillColor(0);
267  lups->SetFillStyle(0);
268  if(yrs13 && yrs15)
269  {
270  lups->AddEntry(gr1[1],"Y(1S)","p");
271  lups->AddEntry(gr2[1],"Y(2S)","p");
272  if(ups3s) lups->AddEntry(gr3[1],"Y(3S)","p");
273  }
274  else
275  {
276  lups->AddEntry(gr1[0],"Y(1S)","p");
277  lups->AddEntry(gr2[0],"Y(2S)","p");
278  if(ups3s) lups->AddEntry(gr3[0],"Y(3S)","p");
279  }
280  lups->Draw();
281 
282  TLatex *lat1[2];
283  lat1[0] = new TLatex(0.66, 0.833,"#font[42]{Au+Au, Years 1-3}");
284  lat1[0]->SetTextSize(0.038);
285  lat1[0]->SetNDC(1);
286  lat1[1] = new TLatex(0.66, 0.833,"#font[42]{Au+Au, Years 1-5}");
287  lat1[1]->SetTextSize(0.038);
288  lat1[1]->SetNDC(1);
289 
290  if(yrs13 && !yrs15)
291  lat1[0]->Draw();
292  if(yrs15 || (yrs13 && yrs15))
293  lat1[1]->Draw();
294 
295  TLatex *lat2[2];
296  //lat2[0] = new TLatex(0.66, 0.763,"#font[42]{#splitline{62 pb^{-1} trig. p+p}{142B rec. Au+Au}}");
297  lat2[0] = new TLatex(0.65, 0.763,"#font[42]{#splitline{62 pb^{-1} trig. p+p}{21 nb^{-1} rec. Au+Au}}");
298  lat2[0]->SetTextSize(0.038);
299  lat2[0]->SetNDC(1);
300  //lat2[1] = new TLatex(0.66, 0.763,"#font[42]{#splitline{142 pb^{-1} trig. p+p}{348B rec. Au+Au}}");
301  lat2[1] = new TLatex(0.65, 0.763,"#font[42]{#splitline{142 pb^{-1} trig. p+p}{51 nb^{-1} rec. Au+Au}}");
302  lat2[1]->SetTextSize(0.038);
303  lat2[1]->SetNDC(1);
304 
305  if(yrs13 && !yrs15)
306  lat2[0]->Draw();
307 
308  if(yrs15 || (yrs13 && yrs15))
309  lat2[1]->Draw();
310 
311  TLatex *lat4 = new TLatex(0.57,0.625,"#font[42]{#it{open markers: Years 1-3}}");
312  lat4->SetTextSize(0.038);
313  lat4->SetNDC(1);
314  if(yrs15 || (yrs13 && yrs15))
315  lat4->Draw();
316 
317  // Plot the theory curves also
318  double col[3] = {kBlack, kRed, kBlue};
319  int lstyle[3] = {kDotted, kSolid, kDashDotted};
320  int nstates = 2;
321  if(ups3s) nstates = 3;
322  for(int is = 0;is<nstates;is++)
323  {
324  for(int ieta =0; ieta < 3; ++ieta)
325  {
326  grth[is][ieta]->SetLineColor(col[is]);
327  grth[is][ieta]->SetMarkerSize(0.5);
328  grth[is][ieta]->SetMarkerColor(col[is]);
329  grth[is][ieta]->SetMarkerStyle(20);
330  grth[is][ieta]->SetLineStyle(lstyle[ieta]);
331  grth[is][ieta]->SetLineWidth(2);
332  grth[is][ieta]->Draw("l");
333  }
334  }
335 
336  TLatex *lat3 = new TLatex(0.4, 0.675,"#font[42]{Strickland,Bazow, N.P. A879 (2012) 25}");
337  lat3->SetTextSize(0.038);
338  lat3->SetNDC(1);
339  lat3->Draw();
340 
341  TLegend *letas = new TLegend(0.14, 0.14, 0.31, 0.28,"");
342  letas->SetBorderSize(0);
343  letas->SetFillColor(0);
344  letas->SetFillStyle(0);
345  letas->AddEntry(grth[0][0], "4#pi#eta/S=1","l");
346  letas->AddEntry(grth[0][1], "4#pi#eta/S=2","l");
347  letas->AddEntry(grth[0][2], "4#pi#eta/S=3","l");
348  letas->Draw();
349 }
350 
351 
352