Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
newplotbg_newsupp_2022_60pct.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file newplotbg_newsupp_2022_60pct.C
1 double CBFunction(double *x, double *p)
2 {
3  double norm = p[0];
4  double alpha = p[1]; // tail start
5  double n = p[2]; // tail shape
6  double mu = p[3]; // upsilon mass
7  double sigma = p[4]; // upsilon width
8 
9  double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
10  double B = n/fabs(alpha) - fabs(alpha);
11  double k = (x[0]-mu)/sigma;
12 
13  double val;
14  if( k > -alpha )
15  val = norm*TMath::Exp(-0.5*pow(k,2));
16  else
17  val = norm*A*pow(B-k,-n);
18 
19  if( TMath::IsNaN(val) ) val = 0.0;
20 
21  return val;
22 }
23 
24 double TripleCBFunction(double *x, double *p)
25 {
26  double norm1 = p[0]; // amplitude of Y(1S) peak
27  double alpha = p[1]; // tail start
28  double n = p[2]; // tail shape
29  double mu1 = p[3]; // upsilon (1S) mass
30  double sigma = p[4]; // upsilon width
31  double norm2 = p[5];
32  double norm3 = p[6];
33  double mu2 = mu1*1.0595;
34  double mu3 = mu1*1.0946;
35 
36  double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
37  double B = n/fabs(alpha) - fabs(alpha);
38  double k1 = (x[0]-mu1)/sigma;
39  double k2 = (x[0]-mu2)/sigma;
40  double k3 = (x[0]-mu3)/sigma;
41 
42  double val,val1,val2,val3;
43 
44  if( k1 > -alpha ) { val1 = norm1*TMath::Exp(-0.5*pow(k1,2)); }
45  else { val1 = norm1*A*pow(B-k1,-n); }
46  if( k2 > -alpha ) { val2 = norm2*TMath::Exp(-0.5*pow(k2,2)); }
47  else { val2 = norm2*A*pow(B-k2,-n); }
48  if( k3 > -alpha ) { val3 = norm3*TMath::Exp(-0.5*pow(k3,2)); }
49  else { val3 = norm3*A*pow(B-k3,-n); }
50 
51  val = val1 + val2 + val3;
52 
53  if( TMath::IsNaN(val) ) val = 0.0;
54 
55  return val;
56 }
57 
58 double SandB_CBFunction(double *x, double *p)
59 {
60  double norm1 = p[0]; // amplitude of Y(1S) peak
61  double alpha = p[1]; // tail start
62  double n = p[2]; // tail shape
63  double mu1 = p[3]; // upsilon (1S) mass
64  double sigma = p[4]; // upsilon width
65  double norm2 = p[5];
66  double norm3 = p[6];
67  double mu2 = mu1*1.0595;
68  double mu3 = mu1*1.0946;
69 
70  double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
71  double B = n/fabs(alpha) - fabs(alpha);
72  double k1 = (x[0]-mu1)/sigma;
73  double k2 = (x[0]-mu2)/sigma;
74  double k3 = (x[0]-mu3)/sigma;
75 
76  double val,val1,val2,val3;
77 
78  if( k1 > -alpha ) { val1 = norm1*TMath::Exp(-0.5*pow(k1,2)); }
79  else { val1 = norm1*A*pow(B-k1,-n); }
80  if( k2 > -alpha ) { val2 = norm2*TMath::Exp(-0.5*pow(k2,2)); }
81  else { val2 = norm2*A*pow(B-k2,-n); }
82  if( k3 > -alpha ) { val3 = norm3*TMath::Exp(-0.5*pow(k3,2)); }
83  else { val3 = norm3*A*pow(B-k3,-n); }
84 
85  double bgnorm1 = p[7];
86  double bgslope1 = p[8];
87 
88  double bg = exp(bgnorm1+x[0]*bgslope1);
89 
90  val = val1 + val2 + val3 + bg;
91  if( TMath::IsNaN(val) ) val = 0.0;
92 
93  return val;
94 }
95 
96 //----------------------------------------------------------------
97 //----------------------------------------------------------------
98 //----------------------------------------------------------------
99 
100 void newplotbg_newsupp_2022_60pct() {
101 
102 //gROOT->LoadMacro("sPHENIXStyle/sPhenixStyle.C");
103 //SetsPhenixStyle();
104 
105 gStyle->SetOptStat(0);
106 gStyle->SetOptFit(0);
107 
108 TRandom* myrandom = new TRandom3();
109 const int nbins = 15;
110 double eideff = 0.9;
111 string times = "*";
112 TLatex* tl[15];
113 char tlchar[999];
114 
115 double statscale = 1.4; // from 10B to 14B
116 double statscale_lowlim = 7.0;
117 double statscale_uplim = 14.0;
118 
119 double ppcorr = (2400./14.)/962.; // 2400B pp and 140B MB AuAu
120 
121  TF1* fCBpp = new TF1("fCBpp",CBFunction,5.,14.,5);
122  TF1* fCBauau = new TF1("fCBauau",CBFunction,5.,14.,5);
123  TF1* fCB1s = new TF1("fCB1s",CBFunction,5.,14.,5);
124  TF1* fCB2s = new TF1("fCB2s",CBFunction,5.,14.,5);
125  TF1* fCB3s = new TF1("fCB3s",CBFunction,5.,14.,5);
126  TF1* fTCB = new TF1("fTCB",TripleCBFunction,5.,14.,7);
127  TF1* fTCBpp = new TF1("fTCBpp",TripleCBFunction,5.,14.,7);
128  TF1* fTCBauau = new TF1("fTCBauau",TripleCBFunction,5.,14.,7);
129  TF1* fSandB = new TF1("fSandB",SandB_CBFunction,5.,14.,9);
130  TF1* fSandBfordave = new TF1("fSandBfordave",SandB_CBFunction,5.,14.,9);
131  TF1* fSandBpp = new TF1("fSandBpp",SandB_CBFunction,5.,14.,9);
132  TF1* fSandBauau = new TF1("fSandBauau",SandB_CBFunction,5.,14.,9);
133 
134 //---------------------------------------------------------
135 // Upsilons
136 //---------------------------------------------------------
137 
138 int whichsuppression = 1; // 0=CMS suppression, 1=Strickland suppression
139 
140 int npts1 = 10;
141 int npts2 = 8;
142 int npts3 = 7;
143 int npts3_rebin = 1;
144 if(whichsuppression==0) {npts3_rebin=4;}
145 
146 string str_UpsilonPt = "(2.0*3.14159*x*[0]*pow((1 + x*x/(4*[1]) ),-[2]))"; // dN/dpT
147 string str_UpsilonXPt = "(2.0*3.14159*x*x*[0]*pow((1 + x*x/(4*[1]) ),-[2]))"; // need this for mean pT calculation
148 TF1* fUpsilonPt = new TF1("fUpsilonPt",str_UpsilonPt.c_str(),0.,20.);
149 TF1* fUpsilonXPt = new TF1("fUpsilonXPt",str_UpsilonXPt.c_str(),0.,20.);
150 fUpsilonPt->SetParameters(72.1, 26.516, 10.6834);
151 fUpsilonXPt->SetParameters(72.1, 26.516, 10.6834);
152 double upsnorm = fUpsilonPt->Integral(0.,20.);
153 
154 double xx1[nbins+1]; for(int i=0; i<nbins+1; i++) {xx1[i] = 0.5 + double(i);}
155 double xx2[nbins+1]; for(int i=0; i<nbins+1; i++) {xx2[i] = 0.5 + double(i);}
156 double xx3[nbins+1]; for(int i=0; i<nbins+1; i++) {xx3[i] = 0.5 + double(i);}
157 double xx3_rebin[nbins+1];
158 if(npts3_rebin==2) {
159  xx3_rebin[0] = fUpsilonXPt->Integral(0.,4.)/fUpsilonPt->Integral(0.,4.);
160  xx3_rebin[1] = fUpsilonXPt->Integral(4.,10.)/fUpsilonPt->Integral(4.,10.);
161 } else if(npts3_rebin==4) {
162  xx3_rebin[0] = fUpsilonXPt->Integral(0.,2.)/fUpsilonPt->Integral(0.,2.);
163  xx3_rebin[1] = fUpsilonXPt->Integral(2.,4.)/fUpsilonPt->Integral(2.,4.);
164  xx3_rebin[2] = fUpsilonXPt->Integral(4.,6.)/fUpsilonPt->Integral(4.,6.);
165  xx3_rebin[3] = fUpsilonXPt->Integral(6.,10.)/fUpsilonPt->Integral(6.,10.);
166 } else if(npts3_rebin==1) {
167  xx3_rebin[0] = fUpsilonXPt->Integral(0.,10.)/fUpsilonPt->Integral(0.,10.);
168 } else {cerr << "ERROR: Wrong Y(3S) binning." << endl; return;}
169 
170 
171 // Raa for 0-60% centrality bin
172 double raa1s[nbins+1];
173 double raa2s[nbins+1];
174 double raa3s[nbins+1];
175 for(int i=0; i<nbins+1; i++) {raa1s[i]=0.654; raa2s[i]=0.277; raa3s[i]=0.075;}
176 if(whichsuppression==0) {for(int i=0; i<nbins+1; i++) {raa3s[i]=0.138;}}
177 
178 TGraph* grRAA1S = new TGraph(nbins,xx1,raa1s);
179 TGraph* grRAA2S = new TGraph(nbins,xx2,raa2s);
180 TGraph* grRAA3S = new TGraph(nbins,xx3,raa3s);
181 
182 int nchan=400;
183 double start=0.0;
184 double stop=20.0;
185 
186 double frac[3]; // upsilons fractions
187  frac[0] = 0.7117;
188  frac[1] = 0.1851;
189  frac[2] = 0.1032;
190 double scale[3]; // mass scaling
191  scale[0] = 1.0;
192  scale[1] = 1.0595;
193  scale[2] = 1.0946;
194 
195 //Npart for 0-60% centrality = 160.2
196 //Ncoll for 0-60% centrality = 388.8
197 //Ncoll for 0-10% centrality = 955.
198 //Npionpairs for 0-60% centrality = 0.406667
199 //Number of Upsilons in 0-60% centrality = 19405 1920 549
200 double combbg_corr_60pct = 0.406667 * 6.;
201 double corrbg_corr_60pct = 388.8/955. * 6.;
202 
203 //---------- Use fixed numbers for 2022 estimate --------------------------------------
204 // Numbers from Marzia for 140 B MB Au+Au and 2400 B p+p events
205 // RAA for Y(3S) is 8.5% Number of Upsilons for 0-60% centrality bin
206  int Nups1 = 19405;
207  int Nups2 = 1920;
208  int Nups3 = 549; // with CMS suppression
209  if(whichsuppression==1) {Nups3 = 549*(3.5/8.5);} // with Strickland's suppression
210  int Nups23 = Nups2+Nups3;
211  int Nups1pp = 2.86e+03;
212  int Nups2pp = 7.16e+02;
213  int Nups3pp = 3.98e+02;
214  int Nups23pp = Nups2pp+Nups3pp;
215 
216 //-------------------------------------------------------------------------------------
217 
218 //====================================================================================
219 //====================================================================================
220 //====================================================================================
221 
222 //
223 // these histograms (hhups*) are randomly generated using the fit above to single
224 // peak (fCB function) and used for the RAA uncertainty calculation
225 //
226 // FORCE specific RESOLUTION and tail shape
227  double tonypar1 = 0.98; // Tony's 100 pion simulation April 2019
228  double tonypar2 = 0.93; // Tony's 100 pion simulation April 2019
229  //double tonypar3 = 9.437; // Tony's 100 pion simulation April 2019
230  double tonypar3 = 9.448; // benchmark
231  double tonypar4 = 0.100; // benchmark
232  double tonypar4pp = 0.089; // Tony's 100 pion simulation April 2019
233  fCBpp->SetParameter(0,1000.);
234  fCBpp->SetParameter(1,tonypar1);
235  fCBpp->SetParameter(2,tonypar2);
236  fCBpp->SetParameter(3,tonypar3);
237  fCBpp->SetParameter(4,tonypar4pp);
238  fCBauau->SetParameter(0,1000.);
239  fCBauau->SetParameter(1,tonypar1);
240  fCBauau->SetParameter(2,tonypar2);
241  fCBauau->SetParameter(3,tonypar3);
242  fCBauau->SetParameter(4,tonypar4);
243 
244 
245 char hhname[999];
246 TH1D* hhups[nbins+1];
247 TH1D* hhups1[nbins+1];
248 TH1D* hhups2[nbins+1];
249 TH1D* hhups3[nbins+1];
250 TH1D* hhupspp[nbins+1];
251 TH1D* hhups1pp[nbins+1];
252 TH1D* hhups2pp[nbins+1];
253 TH1D* hhups3pp[nbins+1];
254 for(int i=0; i<nbins+1; i++) {
255  sprintf(hhname,"hhups_%d",i);
256  hhups[i] = new TH1D(hhname,"",nchan,start,stop);
257  hhups[i]->Sumw2();
258  sprintf(hhname,"hhups1_%d",i);
259  hhups1[i] = new TH1D(hhname,"",nchan,start,stop);
260  hhups1[i]->Sumw2();
261  sprintf(hhname,"hhups2_%d",i);
262  hhups2[i] = new TH1D(hhname,"",nchan,start,stop);
263  hhups2[i]->Sumw2();
264  sprintf(hhname,"hhups3_%d",i);
265  hhups3[i] = new TH1D(hhname,"",nchan,start,stop);
266  hhups3[i]->Sumw2();
267  hhups[i]->SetLineWidth(2);
268  hhups1[i]->SetLineWidth(2);
269  hhups2[i]->SetLineWidth(2);
270  hhups3[i]->SetLineWidth(2);
271  sprintf(hhname,"hhupspp_%d",i);
272  hhupspp[i] = new TH1D(hhname,"",nchan,start,stop);
273  hhupspp[i]->Sumw2();
274  sprintf(hhname,"hhups1pp_%d",i);
275  hhups1pp[i] = new TH1D(hhname,"",nchan,start,stop);
276  hhups1pp[i]->Sumw2();
277  sprintf(hhname,"hhups2pp_%d",i);
278  hhups2pp[i] = new TH1D(hhname,"",nchan,start,stop);
279  hhups2pp[i]->Sumw2();
280  sprintf(hhname,"hhups3pp_%d",i);
281  hhups3pp[i] = new TH1D(hhname,"",nchan,start,stop);
282  hhups3pp[i]->Sumw2();
283  hhupspp[i]->SetLineWidth(2);
284  hhups1pp[i]->SetLineWidth(2);
285  hhups2pp[i]->SetLineWidth(2);
286  hhups3pp[i]->SetLineWidth(2);
287 }
288 
289 for(int j=0; j<nbins; j++) {
290 
291  double s1 = j*1.0;
292  double s2 = s1 + 1.0;
293  double tmpnups1 = Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm;
294  double tmpnups2 = Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm;
295  double tmpnups3 = Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm;
296  fCBauau->SetParameter(3,tonypar3);
297  for(int i=0; i<int(tmpnups1+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups1[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
298  fCBauau->SetParameter(3,tonypar3*scale[1]);
299  for(int i=0; i<int(tmpnups2+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups2[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
300  fCBauau->SetParameter(3,tonypar3*scale[2]);
301  for(int i=0; i<int(tmpnups3+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups3[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
302 
303  double tmpnups1pp = Nups1pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
304  double tmpnups2pp = Nups2pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
305  double tmpnups3pp = Nups3pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
306  fCBpp->SetParameter(3,tonypar3);
307  for(int i=0; i<int(tmpnups1pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups1pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
308  fCBpp->SetParameter(3,tonypar3*scale[1]);
309  for(int i=0; i<int(tmpnups2pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups2pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
310  fCBpp->SetParameter(3,tonypar3*scale[2]);
311  for(int i=0; i<int(tmpnups3pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups3pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
312 
313 } // end loop over pT bins
314 
315 // all pT
316 
317  fCBpp->SetParameter(3,tonypar3);
318  fCBauau->SetParameter(3,tonypar3);
319  for(int i=0; i<int(Nups1+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups1[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
320  for(int i=0; i<int(Nups1pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups1pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
321  fCBpp->SetParameter(3,tonypar3*scale[1]);
322  fCBauau->SetParameter(3,tonypar3*scale[1]);
323  for(int i=0; i<int(Nups2+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups2[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
324  for(int i=0; i<int(Nups2pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups2pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
325  fCBpp->SetParameter(3,tonypar3*scale[2]);
326  fCBauau->SetParameter(3,tonypar3*scale[2]);
327  for(int i=0; i<int(Nups3+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups3[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
328  for(int i=0; i<int(Nups3pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups3pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
329 
330 // fit and draw pp no bg
331 
332 TCanvas* cupspp = new TCanvas("cupspp","Upsilons in p+p",100,100,600,600);
333  fTCBpp->SetParameter(0,2000.);
334  fTCBpp->FixParameter(1,tonypar1);
335  fTCBpp->FixParameter(2,tonypar2);
336  fTCBpp->SetParameter(3,tonypar3);
337  fTCBpp->FixParameter(4,tonypar4); // fix width from the single peak fit
338  fTCBpp->SetParameter(5,500.);
339  fTCBpp->SetParameter(6,100.);
340  hhupspp[nbins]->Fit(fTCBpp,"rl","",7.,11.);
341  hhupspp[nbins]->SetAxisRange(7.,11.);
342  hhupspp[nbins]->SetMarkerSize(1.0);
343  hhupspp[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
344  hhupspp[nbins]->GetXaxis()->SetTitleOffset(1.0);
345  double tmpamp1 = hhupspp[nbins]->GetFunction("fTCBpp")->GetParameter(0);
346  double tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
347  double tmpamp6 = tmpamp1*frac[2]/frac[0];
348  hhupspp[nbins]->Draw();
349 
350  fCB1s->SetLineColor(kBlue);
351  fCB1s->SetLineWidth(1);
352  fCB1s->SetParameter(0,fTCBpp->GetParameter(0));
353  fCB1s->SetParameter(1,fTCBpp->GetParameter(1));
354  fCB1s->SetParameter(2,fTCBpp->GetParameter(2));
355  fCB1s->SetParameter(3,fTCBpp->GetParameter(3)*scale[0]);
356  fCB1s->SetParameter(4,fTCBpp->GetParameter(4));
357  fCB2s->SetLineColor(kRed);
358  fCB2s->SetLineWidth(1);
359  fCB2s->SetParameter(0,tmpamp5);
360  fCB2s->SetParameter(1,fTCBpp->GetParameter(1));
361  fCB2s->SetParameter(2,fTCBpp->GetParameter(2));
362  fCB2s->SetParameter(3,fTCBpp->GetParameter(3)*scale[1]);
363  fCB2s->SetParameter(4,fTCBpp->GetParameter(4));
364  fCB3s->SetLineColor(kGreen+2);
365  fCB3s->SetLineWidth(1);
366  fCB3s->SetParameter(0,tmpamp6);
367  fCB3s->SetParameter(1,fTCBpp->GetParameter(1));
368  fCB3s->SetParameter(2,fTCBpp->GetParameter(2));
369  fCB3s->SetParameter(3,fTCBpp->GetParameter(3)*scale[2]);
370  fCB3s->SetParameter(4,fTCBpp->GetParameter(4));
371  fCB1s->Draw("same");
372  fCB2s->Draw("same");
373  fCB3s->Draw("same");
374 
375 // fit and draw AuAu no bg
376 
377 TCanvas* cupsauau = new TCanvas("cupsauau","Upsilons in Au+Au",100,100,600,600);
378  fTCBauau->SetParameter(0,2000.);
379  fTCBauau->FixParameter(1,tonypar1);
380  fTCBauau->FixParameter(2,tonypar2);
381  fTCBauau->SetParameter(3,tonypar3);
382  fTCBauau->FixParameter(4,tonypar4); // fix width from the single peak fit
383  fTCBauau->SetParameter(5,500.);
384  fTCBauau->SetParameter(6,100.);
385  hhups[nbins]->Fit(fTCBauau,"rl","",7.,11.);
386  hhups[nbins]->SetAxisRange(7.,11.);
387  hhups[nbins]->SetMarkerSize(1.0);
388  hhups[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
389  hhups[nbins]->GetXaxis()->SetTitleOffset(1.0);
390  tmpamp1 = hhups[nbins]->GetFunction("fTCBauau")->GetParameter(0);
391  tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
392  tmpamp6 = tmpamp1*frac[2]/frac[0];
393  hhups[nbins]->Draw();
394 
395 TCanvas* cupsvspt = new TCanvas("cupsvspt","Upsilons vs. p_{T}",100,100,1200,900);
396 cupsvspt->Divide(4,3);
397 for(int i=0; i<12; i++) {
398  cupsvspt->cd(i+1);
399  hhups[i]->SetAxisRange(7.0,11.0); hhups[i]->Draw();
400  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhups[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
401 }
402 
403 //----------------------------------------------------
404 // Backgrounds
405 //----------------------------------------------------
406 
407 TH1D* hhall[nbins+1];
408 TH1D* hhall_scaled[nbins+1];
409 
410 TH1D* hhtotbg[nbins+1];
411 TH1D* hhtotbg_scaled[nbins+1];
412 TH1D* hhcombbg[nbins+1];
413 TH1D* hhcombbg_scaled[nbins+1];
414 TH1D* hhfakefake[nbins+1];
415 TH1D* hhfakehf[nbins+1];
416 TH1D* hhbottom[nbins+1];
417 TH1D* hhcharm[nbins+1];
418 TH1D* hhdy[nbins+1];
419 TH1D* hhcorrbg[nbins+1];
420 TH1D* hhcorrbg_scaled[nbins+1];
421 TH1D* hhfit[nbins+1];
422 char tmpname[999];
423 
424 //----------------------------------------------------------------------------------------
425 // Correlated BG calculated for 10B central AuAu events
426 //----------------------------------------------------------------------------------------
427 
428 double corrbgfitpar0;
429 double corrbgfitpar1;
430 
431 TFile* f=new TFile("ccbb_eideff09.root");
432 for(int i=0; i<nbins+1; i++) {
433  sprintf(tmpname,"hhbottom_%d",i);
434  hhbottom[i] = (TH1D*)f->Get(tmpname);
435  hhbottom[i]->SetDirectory(gROOT);
436  sprintf(tmpname,"hhcharm_%d",i);
437  hhcharm[i] = (TH1D*)f->Get(tmpname);
438  hhcharm[i]->SetDirectory(gROOT);
439  sprintf(tmpname,"hhdy_%d",i);
440  hhdy[i] = (TH1D*)f->Get(tmpname);
441  hhdy[i]->SetDirectory(gROOT);
442  sprintf(tmpname,"hhcorrbg_%d",i);
443  hhcorrbg[i] = (TH1D*)hhbottom[i]->Clone(tmpname);
444  hhcorrbg[i]->Add(hhcharm[i]);
445  hhcorrbg[i]->Add(hhdy[i]);
446  sprintf(tmpname,"hhcorrbg_scaled_%d",i);
447  hhcorrbg_scaled[i] = (TH1D*)hhcorrbg[i]->Clone(tmpname);
448  hhcorrbg[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
449  hhbottom[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
450  hhdy[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
451  if(i==nbins) {
452  corrbgfitpar0 = hhcorrbg[i]->GetFunction("expo")->GetParameter(0);
453  corrbgfitpar1 = hhcorrbg[i]->GetFunction("expo")->GetParameter(1);
454  }
455  for(int k=1; k<=hhcorrbg[i]->GetNbinsX(); k++) {
456  if(hhcorrbg[i]->GetBinLowEdge(k)<statscale_lowlim || (hhcorrbg[i]->GetBinLowEdge(k)+hhcorrbg[i]->GetBinWidth(k))>statscale_uplim) {
457  hhcorrbg_scaled[i]->SetBinContent(k,0.);
458  hhcorrbg_scaled[i]->SetBinError(k,0.);
459  }
460  else {
461  double tmp = corrbg_corr_60pct * statscale * hhcorrbg[i]->GetFunction("expo")->Eval(hhcorrbg[i]->GetBinCenter(k));
462  double tmprnd = myrandom->Poisson(tmp);
463  if(tmprnd<0.) { tmprnd=0.; }
464  hhcorrbg_scaled[i]->SetBinContent(k,tmprnd);
465  hhcorrbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
466  }
467  }
468  hhcorrbg_scaled[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
469  hhcorrbg[i]->SetDirectory(gROOT);
470  hhcorrbg_scaled[i]->SetDirectory(gROOT);
471 }
472 f->Close();
473 
474 TCanvas* c2 = new TCanvas("c2","Correlated BG (all p_{T})",100,100,600,600);
475 
476 hhbottom[nbins]->SetAxisRange(7.0,14.0);
477 hhbottom[nbins]->SetLineColor(kBlue);
478 hhbottom[nbins]->SetLineWidth(2);
479 
480 hhbottom[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
481 hhbottom[nbins]->GetXaxis()->SetTitleOffset(1.0);
482 hhbottom[nbins]->GetXaxis()->SetTitleColor(1);
483 hhbottom[nbins]->GetXaxis()->SetTitleSize(0.040);
484 hhbottom[nbins]->GetXaxis()->SetLabelSize(0.040);
485 //hhbottom[nbins]->GetYaxis()->SetTitle("Correlated background");
486 hhbottom[nbins]->GetYaxis()->SetTitleOffset(1.3);
487 hhbottom[nbins]->GetYaxis()->SetTitleSize(0.040);
488 hhbottom[nbins]->GetYaxis()->SetLabelSize(0.040);
489 
490 hhcorrbg[nbins]->SetAxisRange(7.0,14.0);
491 hhcorrbg[nbins]->SetMinimum(0.1);
492 hhcorrbg[nbins]->SetLineColor(kBlack);
493 hhcorrbg[nbins]->SetLineWidth(2);
494 hhcorrbg[nbins]->Draw("hist");
495 
496 hhbottom[nbins]->SetMarkerColor(kBlue);
497 hhbottom[nbins]->SetLineColor(kBlue);
498 hhbottom[nbins]->Draw("same");
499 
500 hhdy[nbins]->SetMarkerColor(kGreen+2);
501 hhdy[nbins]->SetLineColor(kGreen+2);
502 hhdy[nbins]->SetLineWidth(2);
503 hhdy[nbins]->Draw("same");
504 
505 hhcharm[nbins]->SetMarkerColor(kRed);
506 hhcharm[nbins]->SetLineColor(kRed);
507 hhcharm[nbins]->SetLineWidth(2);
508 hhcharm[nbins]->Draw("same");
509 
510 TCanvas* c0 = new TCanvas("c0","Correlated BG vs. p_{T} 10B events",100,100,1200,900);
511 c0->Divide(4,3);
512 for(int i=0; i<nbins; i++) {
513  c0->cd(i+1);
514  hhcorrbg[i]->SetAxisRange(7.,14.);
515  hhcorrbg[i]->GetFunction("expo")->SetLineColor(kBlack);
516  hhbottom[i]->GetFunction("expo")->SetLineColor(kBlue);
517  hhdy[i]->GetFunction("expo")->SetLineColor(kGreen+2);
518  hhcorrbg[i]->SetLineColor(kBlack);
519  hhbottom[i]->SetLineColor(kBlue);
520  hhdy[i]->SetLineColor(kGreen+2);
521  hhcorrbg[i]->Draw();
522  hhbottom[i]->Draw("same");
523  hhdy[i]->Draw("same");
524  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhcorrbg[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
525 }
526 
527 TCanvas* c0scaled = new TCanvas("c0scaled","SCALED Correlated BG vs. p_{T}",100,100,1200,900);
528 c0scaled->Divide(4,3);
529 for(int i=0; i<nbins; i++) {
530  c0scaled->cd(i+1);
531  hhcorrbg_scaled[i]->SetAxisRange(7.0,14.0); hhcorrbg_scaled[i]->Draw();
532  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhcorrbg_scaled[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
533 }
534 
535 //-----------------------------------------------------------------------------------------
536 // Correlated bg in p+p
537 //-----------------------------------------------------------------------------------------
538 
539 TH1D* hhbottom_pp[nbins+1];
540 TH1D* hhdy_pp[nbins+1];
541 TH1D* hhcorrbg_pp[nbins+1];
542 TH1D* hhall_pp[nbins+1];
543 
544 TF1* fbottom_nosup_corr = new TF1("fbottom_nosup_corr","[0]+[1]*x",5.,14.);
545 fbottom_nosup_corr->SetParameters(-2.13861, 0.683323);
546 
547 for(int i=0; i<nbins+1; i++) {
548 
549  sprintf(tmpname,"hhbottom_pp_%d",i);
550  hhbottom_pp[i] = (TH1D*)hhbottom[i]->Clone(tmpname);
551  for(int k=1; k<=hhbottom_pp[i]->GetNbinsX(); k++) {
552  if(hhbottom_pp[i]->GetBinLowEdge(k)<statscale_lowlim || (hhbottom_pp[i]->GetBinLowEdge(k)+hhbottom_pp[i]->GetBinWidth(k))>statscale_uplim) {
553  hhbottom_pp[i]->SetBinContent(k,0.);
554  hhbottom_pp[i]->SetBinError(k,0.);
555  }
556  else {
557  double tmp = corrbg_corr_60pct * ppcorr * fbottom_nosup_corr->Eval(hhbottom[i]->GetBinCenter(k)) * hhbottom[i]->GetFunction("expo")->Eval(hhbottom[i]->GetBinCenter(k));
558  double tmprnd = myrandom->Poisson(tmp);
559  if(tmprnd<0.) { tmprnd=0.; }
560  hhbottom_pp[i]->SetBinContent(k,tmprnd);
561  hhbottom_pp[i]->SetBinError(k,sqrt(tmprnd));
562  }
563  }
564 
565  sprintf(tmpname,"hhdy_pp_%d",i);
566  hhdy_pp[i] = (TH1D*)hhdy[i]->Clone(tmpname);
567  for(int k=1; k<=hhdy_pp[i]->GetNbinsX(); k++) {
568  if(hhdy_pp[i]->GetBinLowEdge(k)<statscale_lowlim || (hhdy_pp[i]->GetBinLowEdge(k)+hhdy_pp[i]->GetBinWidth(k))>statscale_uplim) {
569  hhdy_pp[i]->SetBinContent(k,0.);
570  hhdy_pp[i]->SetBinError(k,0.);
571  }
572  else {
573  double tmp = ppcorr * hhdy[i]->GetFunction("expo")->Eval(hhdy[i]->GetBinCenter(k));
574  double tmprnd = myrandom->Poisson(tmp);
575  if(tmprnd<0.) { tmprnd=0.; }
576  hhdy_pp[i]->SetBinContent(k,tmprnd);
577  hhdy_pp[i]->SetBinError(k,sqrt(tmprnd));
578  }
579  }
580 
581  sprintf(tmpname,"hhcorrbg_pp_%d",i);
582  hhcorrbg_pp[i] = (TH1D*)hhbottom_pp[i]->Clone(tmpname);
583  hhcorrbg_pp[i]->Add(hhdy_pp[i]);
584  hhcorrbg_pp[i]->SetMarkerColor(kBlack);
585  hhcorrbg_pp[i]->SetLineColor(kBlack);
586  hhbottom_pp[i]->SetLineColor(kBlue);
587  hhdy_pp[i]->SetLineColor(kGreen+2);
588  sprintf(tmpname,"hhall_pp_%d",i);
589  hhall_pp[i] = (TH1D*)hhcorrbg_pp[i]->Clone(tmpname);
590  hhall_pp[i]->Add(hhupspp[i]);
591  hhall_pp[i]->SetLineColor(kMagenta);
592  hhall_pp[i]->SetMarkerColor(kMagenta);
593 
594  hhcorrbg_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
595  hhcorrbg_pp[i]->GetFunction("expo")->SetLineColor(kBlack);
596  hhbottom_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
597  hhbottom_pp[i]->GetFunction("expo")->SetLineColor(kBlue);
598  hhdy_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
599  hhdy_pp[i]->GetFunction("expo")->SetLineColor(kGreen+2);
600 
601 } // end loop over pT bins
602 
603 TCanvas* cbginpp = new TCanvas("cbginpp","corr bg in pp",10,10,700,700);
604  hhall_pp[nbins]->SetAxisRange(7.,12.);
605  hhcorrbg_pp[nbins]->Draw("pehist");
606  hhbottom_pp[nbins]->Draw("histsame");
607  hhdy_pp[nbins]->Draw("histsame");
608 
609 
610 TCanvas* cpp = new TCanvas("cpp","corr bg + sig in pp",100,100,700,700);
611  hhall_pp[nbins]->SetAxisRange(7.,12.);
612  hhall_pp[nbins]->Draw("pehist");
613  hhcorrbg_pp[nbins]->Draw("pesame");
614  hhbottom_pp[nbins]->Draw("same");
615  hhdy_pp[nbins]->Draw("same");
616 
617 TCanvas* cpp_vspt = new TCanvas("cpp_vspt","corr bg + sig vs pt in pp",50,50,1200,900);
618 cpp_vspt->Divide(4,3);
619 for(int i=0; i<nbins; i++) {
620  cpp_vspt->cd(i+1);
621  hhall_pp[i]->SetAxisRange(7.0,12.0); hhall_pp[i]->Draw("hist"); hhcorrbg_pp[i]->Draw("histsame");
622  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhcorrbg_pp[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
623 }
624 
625 //-----------------------------------------------------------------------------------------
626 // Combinatorial BG calculated for 10B central AuAu events
627 //-----------------------------------------------------------------------------------------
628 TCanvas* cdummy = new TCanvas("cdummy","cdummy",0,0,500,500);
629 
630 f = new TFile("fakee_eideff09.root");
631 for(int i=0; i<nbins+1; i++) {
632  sprintf(tmpname,"hhfakefake_%d",i);
633  hhfakefake[i] = (TH1D*)f->Get(tmpname);
634  hhfakefake[i]->SetDirectory(gROOT);
635 }
636 f->Close();
637 
638 f = new TFile("crossterms_eideff09.root");
639 for(int i=0; i<nbins+1; i++) {
640  sprintf(tmpname,"hhfakehf_%d",i);
641  hhfakehf[i] = (TH1D*)f->Get(tmpname);
642  hhfakehf[i]->SetDirectory(gROOT);
643 }
644 f->Close();
645 
646 TF1* fbg = new TF1("fbg","exp([0]+[1]*x)+exp([2]+[3]*x)",8.,11.);
647 fbg->SetParameters(10., -1.0, 4., -0.1);
648 fbg->SetParLimits(1.,-999.,0.);
649 fbg->SetParLimits(3.,-999.,0.);
650 
651 for(int i=0; i<nbins+1; i++) {
652  sprintf(tmpname,"hhcombbg_%d",i);
653  hhcombbg[i] = (TH1D*)hhfakefake[i]->Clone(tmpname);
654  hhcombbg[i]->Add(hhfakehf[i]);
655  sprintf(tmpname,"hhcombbg_scaled_%d",i);
656  hhcombbg_scaled[i] = (TH1D*)hhcombbg[i]->Clone(tmpname);
657  if(i==nbins) { fbg->SetParameters(10., -1.0, 4., -0.1); }
658  hhcombbg[i]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
659 
660  for(int k=1; k<=hhcombbg[i]->GetNbinsX(); k++) {
661  if(hhcombbg[i]->GetBinLowEdge(k)<statscale_lowlim || (hhcombbg[i]->GetBinLowEdge(k)+hhcombbg[i]->GetBinWidth(k))>statscale_uplim) {
662  hhcombbg_scaled[i]->SetBinContent(k,0.);
663  hhcombbg_scaled[i]->SetBinError(k,0.);
664  }
665  else {
666  double tmp = combbg_corr_60pct * statscale * hhcombbg[i]->GetFunction("fbg")->Eval(hhcombbg[i]->GetBinCenter(k));
667  double tmprnd = myrandom->Poisson(tmp);
668  if(tmprnd<0.) { tmprnd=0.; }
669  hhcombbg_scaled[i]->SetBinContent(k,tmprnd);
670  hhcombbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
671  }
672  }
673  hhcombbg_scaled[i]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
674 }
675 
676 delete cdummy;
677 
678 TCanvas* C1 = new TCanvas("C1","Combinatorial BG (ALL p_{T})",100,100,600,600);
679 C1->SetLogy();
680  hhfakefake[nbins]->SetAxisRange(7.0,14.0);
681  hhfakefake[nbins]->SetMinimum(0.1);
682  hhfakefake[nbins]->SetMaximum(5000.);
683  hhfakefake[nbins]->SetLineColor(kGreen+2);
684  hhfakefake[nbins]->SetLineWidth(2);
685  hhfakefake[nbins]->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
686  hhfakefake[nbins]->GetXaxis()->SetTitleOffset(1.0);
687  hhfakefake[nbins]->GetXaxis()->SetTitleColor(1);
688  hhfakefake[nbins]->GetXaxis()->SetTitleSize(0.040);
689  hhfakefake[nbins]->GetXaxis()->SetLabelSize(0.040);
690  hhfakefake[nbins]->GetYaxis()->SetTitle("Combinatorial background");
691  hhfakefake[nbins]->GetYaxis()->SetTitleOffset(1.3);
692  hhfakefake[nbins]->GetYaxis()->SetTitleSize(0.040);
693  hhfakefake[nbins]->GetYaxis()->SetLabelSize(0.040);
694  hhfakefake[nbins]->Draw("e");
695 
696  hhfakehf[nbins]->SetLineColor(kOrange+4);
697  hhfakehf[nbins]->SetLineWidth(2);
698  hhfakehf[nbins]->Draw("esame");
699 
700  hhcombbg[nbins]->SetLineColor(kBlack);
701  hhcombbg[nbins]->SetLineWidth(2);
702  hhcombbg[nbins]->Draw("esame");
703 
704 TCanvas* C1sc = new TCanvas("C1sc","SCALED Combinatorial BG (ALL p_{T})",100,100,600,600);
705 C1sc->SetLogy();
706  hhcombbg_scaled[nbins]->SetAxisRange(7.,14.);
707  hhcombbg_scaled[nbins]->Draw("esame");
708 
709 TCanvas* c00 = new TCanvas("c00","Combinatorial BG vs. p_{T}",150,150,1200,900);
710 c00->Divide(4,3);
711 
712 for(int i=0; i<nbins; i++) {
713  c00->cd(i+1);
714  hhcombbg[i]->SetAxisRange(7.0,14.0); hhcombbg[i]->Draw();
715  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhcombbg[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
716 }
717 
718 TCanvas* c00scaled = new TCanvas("c00scaled","SCALED Combinatorial BG vs. p_{T}",150,150,1200,900);
719 c00scaled->Divide(4,3);
720 
721 for(int i=0; i<nbins; i++) {
722  c00scaled->cd(i+1);
723  hhcombbg_scaled[i]->SetAxisRange(statscale_lowlim,statscale_uplim); hhcombbg_scaled[i]->Draw();
724  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhcombbg_scaled[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
725 }
726 
727 //---------------------------------------------
728 // Signal + BG
729 //---------------------------------------------
730 
731 for(int i=0; i<nbins+1; i++) {
732  sprintf(tmpname,"hhtotbg_scaled_%d",i);
733  hhtotbg_scaled[i] = (TH1D*)hhcombbg_scaled[i]->Clone(tmpname);
734  hhtotbg_scaled[i]->Add(hhcorrbg_scaled[i]);
735 }
736 
737 for(int i=0; i<nbins+1; i++) {
738 sprintf(tmpname,"hhall_scaled_%d",i);
739 hhall_scaled[i] = (TH1D*)hhtotbg_scaled[i]->Clone(tmpname);
740 hhall_scaled[i]->Add(hhups[i]);
741 }
742 
743 TCanvas* c000 = new TCanvas("c000","Signal + Background vs. p_{T}",200,200,1200,900);
744 c000->Divide(4,3);
745 for(int i=0; i<nbins; i++) {
746  c000->cd(i+1);
747  hhall_scaled[i]->SetAxisRange(7.0,14.0); hhall_scaled[i]->SetMarkerStyle(1); hhall_scaled[i]->Draw("pehist");
748  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhall_scaled[i]->GetMaximum()*0.9,tlchar); tl[i]->Draw();
749 }
750 
751 // Fit Signal + Correlated BG
752 
753 // hhfit[] is signal + correlated bg
754 for(int i=0; i<nbins+1; i++) {
755  sprintf(tmpname,"hhfit_%d",i);
756  hhfit[i] = (TH1D*)hhall_scaled[i]->Clone(tmpname);
757  for(int j=1; j<=hhall_scaled[i]->GetNbinsX(); j++) {
758  hhfit[i]->SetBinContent(j,hhfit[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j)));
759  hhfit[i]->SetBinError(j,sqrt(pow(hhfit[i]->GetBinError(j),2)+hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j))));
760  }
761 }
762 
763 //---------------------------------------------------------------------
764 // plot signal + correlated bg for all pT
765 //---------------------------------------------------------------------
766 
767  double tmppar0 = corrbgfitpar0+TMath::Log(statscale);
768  double tmppar1 = corrbgfitpar1;
769 //cout << "###### " << tmppar0 << " " << tmppar1 << endl;
770 
771 double ppauauscale = fTCBauau->GetParameter(0)/fTCBpp->GetParameter(0)*0.9783;
772 fSandBpp->SetParameter(0,fTCBpp->GetParameter(0)*ppauauscale);
773 fSandBpp->SetParameter(1,fTCBpp->GetParameter(1));
774 fSandBpp->SetParameter(2,fTCBpp->GetParameter(2));
775 fSandBpp->SetParameter(3,fTCBpp->GetParameter(3));
776 fSandBpp->SetParameter(4,fTCBpp->GetParameter(4));
777 fSandBpp->SetParameter(5,fTCBpp->GetParameter(5)*ppauauscale);
778 fSandBpp->SetParameter(6,fTCBpp->GetParameter(6)*ppauauscale);
779 fSandBpp->SetParameter(7,tmppar0);
780 fSandBpp->SetParameter(8,tmppar1);
781 fSandBpp->SetLineColor(kBlue);
782 fSandBpp->SetLineStyle(2);
783 
784 //cout << "######### " << fTCBauau->GetParameter(0) << " " << fTCBauau->GetParameter(5) << " " << fTCBauau->GetParameter(6) << endl;
785 fSandBauau->SetParameter(0,fTCBauau->GetParameter(0));
786 fSandBauau->SetParameter(1,fTCBauau->GetParameter(1));
787 fSandBauau->SetParameter(2,fTCBauau->GetParameter(2));
788 fSandBauau->SetParameter(3,fTCBauau->GetParameter(3));
789 fSandBauau->SetParameter(4,fTCBauau->GetParameter(4));
790 fSandBauau->SetParameter(5,fTCBauau->GetParameter(5));
791 fSandBauau->SetParameter(6,fTCBauau->GetParameter(6));
792 fSandBauau->SetParameter(7,tmppar0);
793 fSandBauau->SetParameter(8,tmppar1);
794 fSandBauau->SetLineColor(kRed);
795 
796 
797 TCanvas* cfitall = new TCanvas("cfitall","FIT all pT",270,270,600,600);
798  hhfit[nbins]->SetAxisRange(7.0,14.);
799  hhfit[nbins]->GetXaxis()->CenterTitle();
800  hhfit[nbins]->GetXaxis()->SetTitle("Mass(e^{+}e^{-}) [GeV/c^2]");
801  hhfit[nbins]->GetXaxis()->SetTitleOffset(1.1);
802  hhfit[nbins]->GetXaxis()->SetLabelSize(0.045);
803  hhfit[nbins]->GetYaxis()->CenterTitle();
804  hhfit[nbins]->GetYaxis()->SetLabelSize(0.045);
805  hhfit[nbins]->GetYaxis()->SetTitle("Events / (50 MeV/c^{2})");
806  hhfit[nbins]->GetYaxis()->SetTitleOffset(1.5);
807  hhfit[nbins]->Draw("pehist");
808  fSandBpp->Draw("same");
809  fSandBauau->Draw("same");
810  //hhcorrbg_scaled[nbins]->SetLineColor(kRed);
811  //hhcorrbg_scaled[nbins]->Scale(0.82);
812  //hhcorrbg_scaled[nbins]->Scale(0.70);
813  //hhcorrbg_scaled[nbins]->Draw("histsame");
814 
815  TF1* fmycorrbg = new TF1("fmycorrbg","exp([0]+[1]*x)",7.,14.);
816  fmycorrbg->SetParameters(tmppar0,tmppar1);
817  fmycorrbg->SetLineStyle(2);
818  fmycorrbg->SetLineColor(kRed);
819  fmycorrbg->Draw("same");
820 
821 
822 double myheight = fTCBauau->GetParameter(0);
823 TLatex* ld1 = new TLatex(10.1,myheight,"sPHENIX Simulation");
824 ld1->SetTextSize(0.035);
825 ld1->Draw();
826 TLatex* ld2 = new TLatex(10.1,myheight-100.,"0-10% Au+Au #sqrt{s} = 200 GeV");
827 ld2->SetTextSize(0.035);
828 ld2->Draw();
829 
830 TCanvas* cfitall2 = new TCanvas("cfitall2","FIT all pT",270,270,600,600);
831 TH1D* hhfit_tmp = (TH1D*)hhfit[nbins]->Clone("hhfit_tmp");
832 hhfit_tmp->SetAxisRange(8.0,11.);
833 hhfit_tmp->Draw("pehist");
834 fSandBauau->Draw("same");
835 fmycorrbg->Draw("same");
836 
837 
838 //----------------------------------------------------------------------
839 
840 // plot signal + both backgrounds for all pT
841 
842 TCanvas* callpt = new TCanvas("callpt","Signal+BG (all p_{T})",300,300,600,600);
843 
844  hhall_scaled[nbins]->GetXaxis()->SetTitle("Invariant mass GeV/c");
845  hhall_scaled[nbins]->SetLineColor(kBlack);
846  hhall_scaled[nbins]->SetMarkerColor(kBlack);
847  hhall_scaled[nbins]->SetMarkerStyle(20);
848  hhall_scaled[nbins]->SetAxisRange(8.0,10.8);
849  hhall_scaled[nbins]->Draw("pehist");
850  hhcombbg_scaled[nbins]->SetLineColor(kBlue);
851  hhcombbg_scaled[nbins]->Draw("histsame");
852  hhcorrbg_scaled[nbins]->SetLineColor(kRed);
853  hhcorrbg_scaled[nbins]->Draw("histsame");
854 
855 //----------------------------------------------------------------
856 // Calculate RAA
857 //----------------------------------------------------------------
858 
859 double u1start = 9.25;
860 double u1stop = 9.65;
861 double u2start = 9.80;
862 double u2stop = 10.20;
863 double u3start = 10.20;
864 double u3stop = 10.55;
865 
866  double raa1[nbins+1],raa2[nbins+1],raa3[nbins+1],erraa1[nbins+1],erraa2[nbins+1],erraa3[nbins+1];
867  for(int i=0; i<nbins; i++) {
868  raa1[i] = grRAA1S->Eval(0.5+i*1.0);
869  raa2[i] = grRAA2S->Eval(0.5+i*1.0);
870  raa3[i] = grRAA3S->Eval(0.5+i*1.0);
871  }
872  int fbin1 = hhall_scaled[nbins]->FindBin(u1start + 0.001);
873  int lbin1 = hhall_scaled[nbins]->FindBin(u1stop - 0.001);
874  int fbin2 = hhall_scaled[nbins]->FindBin(u2start + 0.001);
875  int lbin2 = hhall_scaled[nbins]->FindBin(u2stop - 0.001);
876  int fbin3 = hhall_scaled[nbins]->FindBin(u3start + 0.001);
877  int lbin3 = hhall_scaled[nbins]->FindBin(u3stop - 0.001);
878  cout << "Y(1S) bin range: " << fbin1 << " - " << lbin1 << endl;
879  cout << "Y(1S) inv. mass range: " << u1start << " - " << u1stop << endl;
880  cout << "Y(2S) bin range: " << fbin2 << " - " << lbin2 << endl;
881  cout << "Y(2S) inv. mass range: " << u2start << " - " << u2stop << endl;
882  cout << "Y(3S) bin range: " << fbin3 << " - " << lbin3 << endl;
883  cout << "Y(3S) inv. mass range: " << u3start << " - " << u3stop << endl;
884 
885  double sum1[99] = {0.};
886  double truesum1[99] = {0.};
887  double ersum1[99] = {0.};
888  double sumpp1[99] = {0.};
889  double ersumpp1[99] = {0.};
890  double sum2[99] = {0.};
891  double truesum2[99] = {0.};
892  double ersum2[99] = {0.};
893  double sumpp2[99] = {0.};
894  double ersumpp2[99] = {0.};
895  double sum3[99] = {0.};
896  double truesum3[99] = {0.};
897  double ersum3[99] = {0.};
898  double sumpp3[99] = {0.};
899  double ersumpp3[99] = {0.};
900 
901  double sumsum1[99] = {0.};
902  double sumsum2[99] = {0.};
903  double sumsum3[99] = {0.};
904  double sumsum1pp[99] = {0.};
905  double sumsum2pp[99] = {0.};
906  double sumsum3pp[99] = {0.};
907 
908  for(int i=0; i<nbins+1; i++) {
909 
910  double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
911 
912  for(int j=fbin1; j<=lbin1; j++) {
913  sum1[i] += (hhall_scaled[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhall_scaled[i]->GetBinCenter(j)) - hhcorrbg_scaled[i]->GetFunction("expo")->Eval(hhall_scaled[i]->GetBinCenter(j)));
914  truesum1[i] += hhups1[i]->GetBinContent(j);
915  ersum1[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
916  sumpp1[i] += hhups1pp[i]->GetBinContent(j);
917  ersumpp1[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
918  }
919  sumsum1[i] = truesum1[i]; // direct count in mass range
920  sumsum1pp[i] = sumpp1[i];
921  //sumsum1[i] = hhups1[i]->GetEntries(); // total number of upsilons in pT bin (rounded up)
922  //sumsum1pp[i] = hhups1pp[i]->GetEntries();
923  //sumsum1[i] = Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm; // total number of upsilons in pT bin
924  //sumsum1pp[i] = Nups1pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
925 
926  if(sumsum1[i]>0. && sumsum1pp[i]>0.) {
927  erraa1[i] = raa1[i]*sqrt(ersum1[i]/sumsum1[i]/sumsum1[i] + ersumpp1[i]/sumsum1pp[i]/sumsum1pp[i]);
928  } else {raa1[i]=-1.0; erraa1[i] = 999.; }
929 
930  for(int j=fbin2; j<=lbin2; j++) {
931  sum2[i] += (hhall_scaled[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhall_scaled[i]->GetBinCenter(j)) - hhcorrbg_scaled[i]->GetFunction("expo")->Eval(hhall_scaled[i]->GetBinCenter(j)));
932  truesum2[i] += hhups2[i]->GetBinContent(j);
933  ersum2[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
934  sumpp2[i] += hhups2pp[i]->GetBinContent(j);
935  ersumpp2[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
936  }
937  sumsum2[i] = truesum2[i]; // direct count in mass range
938  sumsum2pp[i] = sumpp2[i];
939  //sumsum2[i] = hhups2[i]->GetEntries(); // total number of upsilons in pT bin (rounded up)
940  //sumsum2pp[i] = hhups2pp[i]->GetEntries();
941  //sumsum2[i] = Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm; // total number of upsilons in pT bin
942  //sumsum2pp[i] = Nups2pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
943 
944  if(sumsum2[i]>0. && sumsum2pp[i]>0.) {
945  erraa2[i] = raa2[i]*sqrt(ersum2[i]/sumsum2[i]/sumsum2[i] + ersumpp2[i]/sumsum2pp[i]/sumsum2pp[i]);
946  } else {raa2[i]=-1.0; erraa2[i] = 999.; }
947 
948  for(int j=fbin3; j<=lbin3; j++) {
949  sum3[i] += (hhall_scaled[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhall_scaled[i]->GetBinCenter(j)) - hhcorrbg_scaled[i]->GetFunction("expo")->Eval(hhall_scaled[i]->GetBinCenter(j)));
950  truesum3[i] += hhups3[i]->GetBinContent(j);
951  ersum3[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
952  sumpp3[i] += hhups3pp[i]->GetBinContent(j);
953  ersumpp3[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
954  }
955  sumsum3[i] = truesum3[i]; // direct count in mass range
956  sumsum3pp[i] = sumpp3[i];
957  //sumsum3[i] = hhups3[i]->GetEntries(); // total number of upsilons in pT bin (rounded up)
958  //sumsum3pp[i] = hhups3pp[i]->GetEntries();
959  //sumsum3[i] = Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm; // total number of upsilons in pT bin
960  //sumsum3pp[i] = Nups3pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
961 
962  if(truesum3[i]>0. && sumpp3[i]>0.) {
963  erraa3[i] = raa3[i]*sqrt(ersum3[i]/sumsum3[i]/sumsum3[i] + ersumpp3[i]/sumsum3pp[i]/sumsum3pp[i]);
964  } else {raa3[i]=-1.0; erraa3[i] = 999.; }
965 
966  }
967 
968 double raa2_rebin[9],raapt_rebin2[9],erraa2_rebin[9];
969 double raa3_rebin[9],raapt_rebin3[9],erraa3_rebin[9];
970 double sum2_rebin[9],ersum2_rebin[9],sum2pp_rebin[9],ersumpp2_rebin[9];
971 double sum3_rebin[9],ersum3_rebin[9],sum3pp_rebin[9],ersumpp3_rebin[9];
972 
973 // Rebin Y(3S) by 2
974 if(npts3_rebin==4) {
975 raapt_rebin3[0] = 1.;
976 raapt_rebin3[1] = 3.;
977 raapt_rebin3[2] = 5.;
978 raapt_rebin3[3] = 7.;
979 raa3_rebin[0] = grRAA3S->Eval(raapt_rebin2[0]);
980 raa3_rebin[1] = grRAA3S->Eval(raapt_rebin2[1]);
981 raa3_rebin[2] = grRAA3S->Eval(raapt_rebin2[2]);
982 raa3_rebin[3] = grRAA3S->Eval(raapt_rebin2[3]);
983 sum3_rebin[0] = truesum3[0]+truesum3[1];
984 sum3_rebin[1] = truesum3[2]+truesum3[3];
985 sum3_rebin[2] = truesum3[4]+truesum3[5];
986 sum3_rebin[3] = truesum3[6]+truesum3[7]+truesum3[8]+truesum3[9];
987 ersum3_rebin[0] = ersum3[0]+ersum3[1];
988 ersum3_rebin[1] = ersum3[2]+ersum3[3];
989 ersum3_rebin[2] = ersum3[4]+ersum3[5];
990 ersum3_rebin[3] = ersum3[6]+ersum3[7]+ersum3[8]+ersum3[9];
991 sum3pp_rebin[0] = sumpp3[0]+sumpp3[1];
992 sum3pp_rebin[1] = sumpp3[2]+sumpp3[3];
993 sum3pp_rebin[2] = sumpp3[4]+sumpp3[5];
994 sum3pp_rebin[3] = sumpp3[6]+sumpp3[7]+sumpp3[8]+sumpp3[9];
995 ersumpp3_rebin[0] = ersumpp3[0]+ersumpp3[1];
996 ersumpp3_rebin[1] = ersumpp3[2]+ersumpp3[3];
997 ersumpp3_rebin[2] = ersumpp3[4]+ersumpp3[5];
998 ersumpp3_rebin[3] = ersumpp3[6]+ersumpp3[7]+ersumpp3[8]+ersumpp3[9];
999  erraa3_rebin[0] = raa3[0]*sqrt(ersum3_rebin[0]/sum3_rebin[0]/sum3_rebin[0] + ersumpp3_rebin[0]/sum3pp_rebin[0]/sum3pp_rebin[0]);
1000  erraa3_rebin[1] = raa3[1]*sqrt(ersum3_rebin[1]/sum3_rebin[1]/sum3_rebin[1] + ersumpp3_rebin[1]/sum3pp_rebin[1]/sum3pp_rebin[1]);
1001  erraa3_rebin[2] = raa3[2]*sqrt(ersum3_rebin[2]/sum3_rebin[2]/sum3_rebin[2] + ersumpp3_rebin[2]/sum3pp_rebin[2]/sum3pp_rebin[2]);
1002  erraa3_rebin[3] = raa3[3]*sqrt(ersum3_rebin[3]/sum3_rebin[3]/sum3_rebin[3] + ersumpp3_rebin[3]/sum3pp_rebin[3]/sum3pp_rebin[3]);
1003 }
1004 else if(npts3_rebin==2) {
1005 // Rebin Y(3S) by 4
1006 raapt_rebin3[0] = 2.;
1007 raapt_rebin3[1] = 6.;
1008 raa3_rebin[0] = grRAA3S->Eval(raapt_rebin3[0]);
1009 raa3_rebin[1] = grRAA3S->Eval(raapt_rebin3[1]);
1010 sum3_rebin[0] = truesum3[0]+truesum3[1]+truesum3[2]+truesum3[3];
1011 sum3_rebin[1] = truesum3[4]+truesum3[5]+truesum3[6]+truesum3[7]+truesum3[8]+truesum3[9];
1012 ersum3_rebin[0] = ersum3[0]+ersum3[1]+ersum3[2]+ersum3[3];
1013 ersum3_rebin[1] = ersum3[4]+ersum3[5]+ersum3[6]+ersum3[7]+ersum3[8]+ersum3[9];
1014 sum3pp_rebin[0] = sumpp3[0]+sumpp3[1]+sumpp3[3]+sumpp3[3];
1015 sum3pp_rebin[1] = sumpp3[4]+sumpp3[5]+sumpp3[6]+sumpp3[7]+sumpp3[8]+sumpp3[9];
1016 ersumpp3_rebin[0] = ersumpp3[0]+ersumpp3[1]+ersumpp3[2]+ersumpp3[3];
1017 ersumpp3_rebin[1] = ersumpp3[4]+ersumpp3[5]+ersumpp3[6]+ersumpp3[7]+ersumpp3[8]+ersumpp3[9];
1018  erraa3_rebin[0] = raa3[0]*sqrt(ersum3_rebin[0]/sum3_rebin[0]/sum3_rebin[0] + ersumpp3_rebin[0]/sum3pp_rebin[0]/sum3pp_rebin[0]);
1019  erraa3_rebin[1] = raa3[1]*sqrt(ersum3_rebin[1]/sum3_rebin[1]/sum3_rebin[1] + ersumpp3_rebin[1]/sum3pp_rebin[1]/sum3pp_rebin[1]);
1020 }
1021 else if(npts3_rebin==1) {
1022 // Rebin Y(3S) by 8
1023 raa3_rebin[0] = grRAA3S->Eval(raapt_rebin3[0]);
1024 sum3_rebin[0] = truesum3[0]+truesum3[1]+truesum3[2]+truesum3[3]+
1025  truesum3[4]+truesum3[5]+truesum3[6]+truesum3[7]+truesum3[8]+truesum3[9];
1026 ersum3_rebin[0] = ersum3[0]+ersum3[1]+ersum3[2]+ersum3[3]+
1027  ersum3[4]+ersum3[5]+ersum3[6]+ersum3[7]+ersum3[8]+ersum3[9];
1028 sum3pp_rebin[0] = sumpp3[0]+sumpp3[1]+sumpp3[3]+sumpp3[3]+
1029  sumpp3[4]+sumpp3[5]+sumpp3[6]+sumpp3[7]+sumpp3[8]+sumpp3[9];
1030 ersumpp3_rebin[0] = ersumpp3[0]+ersumpp3[1]+ersumpp3[2]+ersumpp3[3]+
1031  ersumpp3[4]+ersumpp3[5]+ersumpp3[6]+ersumpp3[7]+ersumpp3[8]+ersumpp3[9];
1032 erraa3_rebin[0] = raa3[0]*sqrt(ersum3_rebin[0]/sum3_rebin[0]/sum3_rebin[0] + ersumpp3_rebin[0]/sum3pp_rebin[0]/sum3pp_rebin[0]);
1033 } else {cerr << "ERROR: Wring Y(3S) binning." << endl; return;}
1034 
1035 cout << "====== Y(1S):" << endl;
1036  for(int i=0; i<nbins+1; i++) {
1037  double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1038  cout << " " << i << " " << truesum1[i] << "(" << Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- "
1039  << sqrt(ersum1[i]) << " \t\t pp: " << sumpp1[i] << " +- " << sqrt(ersumpp1[i]) << endl;
1040  }
1041 cout << "====== Y(2S):" << endl;
1042  for(int i=0; i<nbins+1; i++) {
1043  double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1044  cout << " " << i << " " << truesum2[i] << "(" << Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- "
1045  << sqrt(ersum2[i]) << " \t\t pp: " << sumpp2[i] << " +- " << sqrt(ersumpp2[i]) << endl;
1046  }
1047 cout << "====== Y(3S):" << endl;
1048  for(int i=0; i<nbins+1; i++) {
1049 
1050  double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1051  cout << " " << i << " " << truesum3[i] << "(" << Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- "
1052  << sqrt(ersum3[i]) << " \t\t pp: " << sumpp3[i] << " +- " << sqrt(ersumpp3[i]) << endl;
1053  }
1054 
1055 //-------------------------------------------------
1056 // Plot RAA
1057 //-------------------------------------------------
1058 
1059 TCanvas* craa = new TCanvas("craa","R_{AA}",120,120,850,600);
1060 TH2F* hh2 = new TH2F("hh2"," ",10,0.,float(npts1),10,0.,1.1);
1061 hh2->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
1062 hh2->GetXaxis()->SetTitleOffset(1.0);
1063 hh2->GetXaxis()->SetTitleColor(1);
1064 hh2->GetXaxis()->SetTitleSize(0.040);
1065 hh2->GetXaxis()->SetLabelSize(0.040);
1066 hh2->GetYaxis()->SetTitle("R_{AA}");
1067 hh2->GetYaxis()->SetTitleOffset(0.7);
1068 hh2->GetYaxis()->SetTitleSize(0.050);
1069 hh2->GetYaxis()->SetLabelSize(0.040);
1070 hh2->Draw();
1071 
1072 TGraphErrors* gr1 = new TGraphErrors(npts1-1,xx1,raa1,0,erraa1);
1073 gr1->SetMarkerStyle(20);
1074 gr1->SetMarkerColor(kBlack);
1075 gr1->SetLineColor(kBlack);
1076 gr1->SetLineWidth(2);
1077 gr1->SetMarkerSize(1.5);
1078 gr1->SetName("gr1");
1079 gr1->Draw("p");
1080 
1081 TGraphErrors* gr2 = new TGraphErrors(npts2-1,xx2,raa2,0,erraa2);
1082 gr2->SetMarkerStyle(20);
1083 gr2->SetMarkerColor(kRed);
1084 gr2->SetLineColor(kRed);
1085 gr2->SetLineWidth(2);
1086 gr2->SetMarkerSize(1.5);
1087 gr2->SetName("gr2");
1088 gr2->Draw("p");
1089 
1090 //--- 3S state -------------------
1091 TGraphErrors* gr3_rebin;
1092 if(npts3_rebin==4) {
1093  gr3_rebin = new TGraphErrors(npts3_rebin-1,xx3_rebin,raa3_rebin,0,erraa3_rebin);
1094 } else {
1095  gr3_rebin = new TGraphErrors(npts3_rebin,xx3_rebin,raa3_rebin,0,erraa3_rebin);
1096 }
1097  gr3_rebin->SetMarkerStyle(20);
1098  gr3_rebin->SetMarkerColor(kBlue);
1099  gr3_rebin->SetLineColor(kBlue);
1100  gr3_rebin->SetLineWidth(2);
1101  gr3_rebin->SetMarkerSize(1.5);
1102  gr3_rebin->SetName("gr3");
1103  gr3_rebin->Draw("p");
1104 
1105 //TArrow* aa[9];
1106 //TLine* ll[9];
1107 //erraa3[3] = erraa3[3]*1.5;
1108 //erraa3[4] = erraa3[4]*0.6;
1109 //erraa3[4] = erraa3[4]*1.5;
1110 //erraa3[6] = erraa3[6]*1.2;
1111 /*
1112 for(int i=0; i<npts3; i++) {
1113  aa[i] = new TArrow(xx3[i],1.64*erraa3[i],xx3[i],-0.02,0.02);
1114  aa[i]->SetLineColor(kBlue);
1115  aa[i]->SetLineWidth(2);
1116  aa[i]->Draw();
1117  ll[i] = new TLine(xx3[i]-0.15,1.64*erraa3[i],xx3[i]+0.15,1.64*erraa3[i]);
1118  ll[i]->SetLineColor(kBlue);
1119  ll[i]->SetLineWidth(2);
1120  ll[i]->Draw();
1121 }
1122 */
1123 //--- end 3S state --------------
1124 
1125 
1126 TLegend *leg = new TLegend(0.72,0.70,0.88,0.88);
1127 leg->SetBorderSize(0);
1128 leg->SetFillColor(10);
1129 leg->SetFillStyle(1001);
1130 TLegendEntry *entry1=leg->AddEntry("gr1","Y(1S)","p");
1131 TLegendEntry *entry2=leg->AddEntry("gr2","Y(2S)","p");
1132 //TLegendEntry *entry3=leg->AddEntry("gr3","Y(3S)","l");
1133 TLegendEntry *entry3=leg->AddEntry("gr3","Y(3S)","p");
1134 leg->Draw();
1135 
1136 TLatex* l1 = new TLatex(0.5,1.02,"#font[72]{sPHENIX} Projection"); l1->SetTextSize(0.045); l1->SetTextFont(42); l1->Draw();
1137 TLatex* l11 = new TLatex(0.5,0.93,"0-60% cent. Au+Au, Years 1-3"); l11->SetTextSize(0.045); l11->SetTextFont(42); l11->Draw();
1138 TLatex* l2 = new TLatex(0.5,0.84,"21 nb^{-1} rec. Au+Au, 62 pb^{-1} samp. #it{p+p}"); l2->SetTextSize(0.045); l2->SetTextFont(42); l2->Draw();
1139 TLatex* l3 = new TLatex(1.8,0.48,"STAR Prelim. 0-60% Y(2S+3S)"); l3->SetTextColor(kGreen+2); l3->SetTextSize(0.045); l3->SetTextFont(42); l3->Draw();
1140 TMarker* m1 = new TMarker(1.55,0.50,29); m1->SetMarkerSize(2); m1->SetMarkerColor(kGreen+2); m1->Draw();
1141 
1142 TLine* lll = new TLine(0.6,0.64,1.3,0.64);
1143 lll->SetLineColor(kBlue);
1144 lll->SetLineWidth(2);
1145 //lll->Draw();
1146 
1147 //==================================================================================
1148 // STAR
1149 
1150 double star_x[3],star_erx[3],star[3],star_stater[3],star_syser[3];
1151 star_x[0] = 1.;
1152 star_x[1] = 3.;
1153 star_x[2] = 7.;
1154 star_erx[0] = 0.16;
1155 star_erx[1] = 0.16;
1156 star_erx[2] = 0.16;
1157 star[0] = 0.384;
1158 star[1] = 0.249;
1159 star[2] = 0.449;
1160 star_stater[0] = 0.147;
1161 star_stater[1] = 0.113;
1162 star_stater[2] = 0.164;
1163 star_syser[0] = 0.130;
1164 star_syser[1] = 0.086;
1165 star_syser[2] = 0.128;
1166 
1167 TGraphErrors* gr_star = new TGraphErrors(3,star_x,star,0,star_stater);
1168 gr_star->SetMarkerStyle(kFullStar);
1169 gr_star->SetMarkerColor(kGreen+2);
1170 gr_star->SetLineColor(kGreen+2);
1171 gr_star->SetLineWidth(2);
1172 gr_star->SetMarkerSize(2);
1173 gr_star->SetName("gr_star");
1174 
1175 TGraphErrors* grsys_star = new TGraphErrors(3,star_x,star,star_erx,star_syser);
1176 grsys_star->SetMarkerStyle(kFullStar);
1177 grsys_star->SetMarkerColor(kGreen+2);
1178 grsys_star->SetLineColor(kGreen+2);
1179 grsys_star->SetLineWidth(2);
1180 grsys_star->SetMarkerSize(2);
1181 grsys_star->SetName("gr_star");
1182 grsys_star->Draw("5");
1183 gr_star->Draw("p");
1184  gr3_rebin->Draw("p");
1185 
1186 /*
1187 fout = new TFile("test.root","recreate");
1188 for(int i=0; i<nbins+1; i++) {
1189  sprintf(tmpname,"hhfit_%d",i);
1190  hhfit[i]->Write();
1191  hhups1pp[i]->Write();
1192  hhups2pp[i]->Write();
1193  hhups3pp[i]->Write();
1194  hhupspp[i]->Write();
1195  hhups1[i]->Write();
1196  hhups2[i]->Write();
1197  hhups3[i]->Write();
1198  hhups[i]->Write();
1199  hhall_pp[i]->Write();
1200  hhcorrbg_pp[i]->Write();
1201  hhcorrbg_scaled[i]->Write();
1202 }
1203 fout->Write();
1204 fout->Close();
1205 */
1206 
1207 }
1208