]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGDQ/dielectron/TakuAlberica/pair/get_spectrum_dir2.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / TakuAlberica / pair / get_spectrum_dir2.C
CommitLineData
06f630bb 1void get_spectrum_dir2(int icent=1){
2 gStyle->SetPalette(1);
3 gStyle->SetOptStat(0);
4 gStyle->SetOptFit(111111);
5 //gStyle->SetOptFit(0);
6 gStyle->SetOptTitle(0);
7 //gStyle->SetFillColor(10);
8 gStyle->SetCanvasColor(10);
9 gStyle->SetFrameBorderMode(0);
10 gStyle->SetFrameFillColor(0);
11 gStyle->SetCanvasColor(0);
12 gStyle->SetPadBorderSize(0);
13 gStyle->SetCanvasBorderSize(0);
14 //gStyle->SetPadLeftMargin(0.15);
15 gStyle->SetPadLeftMargin(0.125);
16 gStyle->SetPadBottomMargin(0.125);
17 gStyle->SetPadTopMargin(0.1);
18 gStyle->SetTitleYOffset(1.3);
19 //gStyle->SetPadLeftMargin(0.1);
20 gStyle->SetTitleW(0.7);
21 gStyle->SetTitleH(0.1);
22 cout<<"physics is always fun! "<<endl;
23
24 const int ncentbin=9;
25 // int cent_low[ncentbin] ={ 0, 0, 20, 40, 40, 0, 0, 0, 30};
26 // int cent_high[ncentbin]={100, 10, 40, 100, 80, 80, 40, 30, 100};
27
28 // int cent_low_bin[ncentbin]={0, 0, 2, 4, 4, 0, 0, 0, 3 };
29 // int cent_high_bin[ncentbin]={10, 1, 4, 10, 8, 8, 4, 3, 10};
30
31 int cent_low[ncentbin] ={ 0, 10, 20, 30};
32 int cent_high[ncentbin]={10, 20, 30, 50};
33
34 int cent_low_bin[ncentbin] ={0, 1, 2, 3};
35 int cent_high_bin[ncentbin]={1, 2, 3, 5};
36
37
38 TH2D *hmasspt_inc[11][7];
39 TH1D *hmasspt[4][5];
40 TH1D *hmasspt_ls[4][5];
41 TH1D *hdiv[5];
42
43 char outname[100];
44 //sprintf(outname,"spec_histo_loose_eid_pt04_2_cent%d_no_gst.root",icent);
45 sprintf(outname,"spec_histo_pt04_2_cent%d_no_gst.root",icent);
46 //sprintf(outname,"output/spec_histo_loose_eid_2_pt04_cent%d_6.root",icent);
47 //sprintf(outname,"spec_histo_loose_eid_pt04_cent%d_conv2.root",icent);
48 char name[100];
49 // TFile *fin = new TFile("result_histo_nmix40_pt04.root");
50 // TFile *fin = new TFile("result_histo_loose_eid_nmix10_conv2_pt04.root");
51 //TFile *fin = new TFile("output/result_histo_loose_eid_nmix40_rev1_6.root");
52 //TFile *fin = new TFile("result_histo_loose_eid_nmix40_pt04_no_gst.root");
53 TFile *fin = new TFile("result_trig1_cut1.root");
54 //TFile *fin = new TFile("result_histo_nmix20_pt04.root");
55 for(int i=0;i<11;i++){
56 for(int j=0;j<7;j++){
57 sprintf(name,"hmasspt_weight_cent%d_pair%d",i,j);
58 //sprintf(name,"hmasspt_cent%d_pair%d",i,j);
59 hmasspt_inc[i][j] = (TH2D*)fin->Get(name);
60 //cout<<i<<" "<<j<<" "<<hmasspt_inc[i][j]->GetEntries()<<endl;
61 }
62 }
63
64 TH1D *hpt = (TH1D*)hmasspt_inc[0][0]->ProjectionY("hpt");
65
66 TH1F *hstat = (TH1F*)fin->Get("hCentrality");
67 float nevents = hstat->Integral(cent_low[icent], cent_high[icent]-1);
68
69 cout<<nevents<<endl;
70
71 const int nptbin=4;
72 float pt_low[nptbin] = {0.5, 1.0, 1.5, 2.0};
73 float pt_high[nptbin] = {1.0, 1.5, 2.0, 5.0};
74 double mass_r_ptbin[nptbin+2]={0, 0.5, 1.0, 1.5, 2.0, 4.0};
75
76 for(int i=0;i<4;i++){
77 for(int j=0;j<nptbin;j++){
78 sprintf(name,"hmass2_w_pair%d_pt%d", i, j);
79 //sprintf(name,"hmass2_pair%d_pt%d", i, j);
80 hmasspt[i][j] = new TH1D(name, name, 500, 0, 5);
81 hmasspt[i][j]->Sumw2();
82
83 sprintf(name,"hmass2_ls_w_pair%d_pt%d", i, j);
84 //sprintf(name,"hmass2_pair%d_pt%d", i, j);
85 hmasspt_ls[i][j] = new TH1D(name, name, 500, 0, 5);
86 hmasspt_ls[i][j]->Sumw2();
87 }
88 }
89
90 /// first []
91 /// 0-->unlike
92 /// 1-->like
93 /// 2-->mixed unlike
94 /// 3-->mixed like
95
96 /// second [] --> pt
97 /// third [] --> ee or pp
98 float npairs[4][nptbin][2];
99
100 for(int j=0;j<nptbin;j++){
101 int binl = hpt->FindBin(pt_low[j]);
102 int binh = hpt->FindBin(pt_high[j]);
103
104 npairs[0][j][0]=0; npairs[0][j][1]=0;
105 npairs[1][j][0]=0; npairs[1][j][1]=0;
106 npairs[2][j][0]=0; npairs[2][j][1]=0;
107 npairs[3][j][0]=0; npairs[3][j][1]=0;
108
109
110
111 for(int i=cent_low_bin[icent]; i<cent_high_bin[icent]; i++){
112
113 ///////////////////////////////////////////
114 ///// pair0
115 ////////////////////////////////////////////
116 sprintf(name,"hmass_%d_%d_0", i, j);
117 hmasspt[0][j]->Add((TH1D*)hmasspt_inc[i][0]->ProjectionX(name,binl, binh-1));
118 // npairs[0][j][0] += hmasspt_inc[i][0]->ProjectionX(name,binl, binh-1)->Integral();
119 npairs[0][j][0] += hmasspt_inc[i][0]->Integral();
120
121 /////////////////////////////////////////
122 //// pair2 (like sign)
123 ////////////////////////////////////////
124 sprintf(name,"hmass_%d_%d_1", i, j);
125 hmasspt[2][j]->Add((TH1D*)hmasspt_inc[i][1]->ProjectionX(name,binl, binh-1));
126 hmasspt_ls[0][j]->Add((TH1D*)hmasspt_inc[i][1]->ProjectionX(name,binl, binh-1));
127 //npairs[2][j][0] += hmasspt_inc[i][1]->ProjectionX(name,binl, binh-1)->Integral();
128 npairs[2][j][0] += hmasspt_inc[i][1]->Integral();
129
130 sprintf(name,"hmass_%d_%d_2", i, j);
131 hmasspt[2][j]->Add((TH1D*)hmasspt_inc[i][2]->ProjectionX(name,binl, binh-1));
132 //npairs[2][j][1] += hmasspt_inc[i][2]->ProjectionX(name,binl, binh-1)->Integral();
133 hmasspt_ls[1][j]->Add((TH1D*)hmasspt_inc[i][2]->ProjectionX(name,binl, binh-1));
134 npairs[2][j][1] += hmasspt_inc[i][2]->Integral();
135
136 /////////////////////////////////////////
137 //// pair1 (mixed unlike sign)
138 /////////////////////////////////////////
139 sprintf(name,"hmass_%d_%d_1", i, j);
140 hmasspt[1][j]->Add((TH1D*)hmasspt_inc[i][3]->ProjectionX(name,binl, binh-1));
141 //npairs[1][j][0] += hmasspt_inc[i][3]->ProjectionX(name,binl, binh-1)->Integral();
142 npairs[1][j][0] += hmasspt_inc[i][3]->Integral();
143
144 sprintf(name,"hmass_%d_%d_2", i, j);
145 hmasspt[1][j]->Add((TH1D*)hmasspt_inc[i][4]->ProjectionX(name,binl, binh-1));
146 // npairs[1][j][0] += hmasspt_inc[i][4]->ProjectionX(name,binl, binh-1)->Integral();
147 npairs[1][j][0] += hmasspt_inc[i][4]->Integral();
148
149 /////////////////////////////////////////
150 //// pair3 (mixed like sign)
151 ////////////////////////////////////////
152 sprintf(name,"hmass_%d_%d_1", i, j);
153 // hmasspt[3][j]->Add((TH1D*)hmasspt_inc[i][5]->ProjectionX(name,binl, binh-1));
154 hmasspt_ls[2][j]->Add((TH1D*)hmasspt_inc[i][5]->ProjectionX(name,binl, binh-1));
155 //npairs[3][j][0] += hmasspt_inc[i][5]->ProjectionX(name,binl, binh-1)->Integral();
156 npairs[3][j][0] += hmasspt_inc[i][5]->Integral();
157
158 sprintf(name,"hmass_%d_%d_2", i, j);
159 //hmasspt[3][j]->Add((TH1D*)hmasspt_inc[i][6]->ProjectionX(name,binl, binh-1));
160 hmasspt_ls[3][j]->Add((TH1D*)hmasspt_inc[i][6]->ProjectionX(name,binl, binh-1));
161 //npairs[3][j][1] += hmasspt_inc[i][6]->ProjectionX(name,binl, binh-1)->Integral();
162 npairs[3][j][1] += hmasspt_inc[i][6]->Integral();
163
164 }
165 }
166
167 /// for mixed like sign pairs
168 /// entry should be same as like sign pairs
169 /// re-add the hist for mixed like sign pairs using weights
170 for(int j=0;j<nptbin;j++){
171 float frac1 = (npairs[2][j][0]*(npairs[3][j][0]+npairs[3][j][1]))/(npairs[3][j][0]*(npairs[2][j][0]+npairs[2][j][1]));
172 float frac2 = (npairs[2][j][1]*(npairs[3][j][0]+npairs[3][j][1]))/(npairs[3][j][1]*(npairs[2][j][0]+npairs[2][j][1]));
173 hmasspt[3][j]->Add( hmasspt_ls[2][j], frac1); ///mixed like sign
174 hmasspt[3][j]->Add( hmasspt_ls[3][j], frac2); // mixed like sign
175 //hmasspt[3][j]->Add( hmasspt_ls[2][j], 1);
176 //hmasspt[3][j]->Add( hmasspt_ls[3][j], 1);
177 }
178
179
180 for(int i=0;i<4;i++){
181 for(int j=0;j<nptbin;j++){
182 hmasspt[i][j]->Scale(1.0/nevents);
183 hmasspt_ls[i][j]->Scale(1.0/nevents);
184 cout<<hmasspt[i][j]->Integral()<<endl;
185 cout<<npairs[i][j][0]<<" "<<npairs[i][j][1]<<endl;
186 }
187 }
188
189 TCanvas *c0 = new TCanvas("c0","c0",1500,500);
190 c0->Divide(3,1);
191 for(int iptbin=1;iptbin<nptbin;iptbin++){
192 c0->cd(iptbin);
193 gPad->SetGrid(1);
194 gPad->SetLogy();
195 hmasspt[0][iptbin]->SetXTitle("M_{ee}[GeV/c^{2}]");
196 hmasspt[0][iptbin]->SetYTitle("N_{ee}/N_{evt}");
197 hmasspt[0][iptbin]->Draw();
198 }
199
200
201 TCanvas *c1 = new TCanvas("c1","c1",1500,1000);
202 c1->Divide(3,3);
203
204 TF1 *f1 = new TF1("f1","pol0",0.05,2);
205
206 double r[5];
207
208 for(int iptbin=1;iptbin<nptbin;iptbin++){
209 sprintf(name,"hdiv_pt%d",iptbin);
210 hdiv[iptbin] = new TH1D(name, name, 500,0,5);
211 hdiv[iptbin]->SetXTitle("M_{ee}[GeV/c^{2}]");
212 hdiv[iptbin]->SetYTitle("mixed unlike/mixed like");
213 hdiv[iptbin]->Sumw2();
214 //hmasspt[3][iptbin]->Scale(npairs[1][iptbin][0]/(npairs[3][iptbin][0]+npairs[3][iptbin][1]));
215 hdiv[iptbin]->Divide(hmasspt[1][iptbin], hmasspt[3][iptbin]); ///mixed unlike/mixed like
216 /*
217 for(int ii=0;ii<hdiv[iptbin]->GetNbinsX(); ii++){
218 double vent1 = hmasspt[1][iptbin]->GetBinContent(ii+1);
219 double e_vent1 = hmasspt[1][iptbin]->GetBinError(ii+1);
220
221
222 double vent2 = hmasspt_ls[2][iptbin]->GetBinContent(ii+1);
223 double e_vent2 = hmasspt_ls[2][iptbin]->GetBinError(ii+1);
224 double vent3 = hmasspt_ls[3][iptbin]->GetBinContent(ii+1);
225 double e_vent3 = hmasspt_ls[3][iptbin]->GetBinError(ii+1);
226
227
228
229 double val1 = vent1;
230 double val2 = 2*sqrt(vent2*vent3);
231
232 if(val2>0){
233 hdiv[iptbin]->SetBinContent(ii+1, val1/val2);
234 double e_val1 = e_vent1;
235 double e_val2 = val2/2*sqrt(pow(e_vent2/vent2,2)+pow(e_vent3/vent3,2));
236 hdiv[iptbin]->SetBinError(ii+1, val1/val2*sqrt(pow(e_val1/val1,2)+pow(e_val2/val2,2)));
237 }
238 }
239 */
240
241 c1->cd(iptbin);
242 gPad->SetGrid(1);
243 hmasspt[1][iptbin]->SetXTitle("M_{ee}[GeV/c^{2}]");
244 hmasspt[0][iptbin]->SetYTitle("N_{ee}/N_{evt}");
245 hmasspt[1][iptbin]->Draw();
246 hmasspt[3][iptbin]->SetLineColor(2);
247 hmasspt[3][iptbin]->Draw("same");
248
249 c1->cd(iptbin+3);
250 gPad->SetGrid(1);
251 //hdiv[iptbin]->Rebin(4); hdiv[iptbin]->Scale(1.0/4);
252
253 hdiv[iptbin]->SetAxisRange(0,3);
254 hdiv[iptbin]->SetMinimum(0);
255 hdiv[iptbin]->SetMaximum(2);
256 hdiv[iptbin]->Draw();
257 ///hdiv[iptbin]->Fit(f1,"R","");
258
259 //r[iptbin] = f1->GetParameter(0);
260 //int bin1 = hmasspt[2][iptbin]->FindBin(0.1);
261 //int bin2 = hmasspt[2][iptbin]->FindBin(1.0);
262 //r[iptbin] = hmasspt[2][iptbin]->Integral(bin1, bin2)/hmasspt[3][iptbin]->Integral(bin1, bin2);
263 //cout<<" normalization factor "<<r[iptbin]<<endl;
264
265 /*
266 for(int k=0;k<hmasspt[1][iptbin]->GetNbinsX();k++){
267 float w = f1->Eval(hmasspt[1][iptbin]->GetBinCenter(k+1));
268 float rr = hmasspt[1][iptbin]->GetBinContent(k+1);
269 float err = hmasspt[1][iptbin]->GetBinError(k+1);
270 hmasspt[1][iptbin]->SetBinContent(k+1, rr*w);
271 hmasspt[1][iptbin]->SetBinError(k+1, err*w);
272 }
273 */
274 //hmasspt[1][iptbin]->Scale(r[iptbin]);
275 hmasspt[2][iptbin]->Multiply(hdiv[iptbin]);
276 //c/out<<npairs[1][iptbin][0]/(2*sqrt(npairs[3][iptbin][0]*npairs[3][iptbin][1]))<<endl;
277 //cout<<" aaa "<<(2*sqrt(npairs[2][iptbin][0]*npairs[2][iptbin][1]))<<" "<<npairs[2][iptbin][0]+npairs[2][iptbin][1]<<" "<<hmasspt[2][iptbin]->Integral()*nevents<<endl;
278 //hmasspt[2][iptbin]->Scale((2*sqrt(npairs[2][iptbin][0]*npairs[2][iptbin][1])/(npairs[2][iptbin][0]+npairs[2][iptbin][1])));
279 //hmasspt[2][iptbin]->Scale(npairs[0][iptbin][0]/(2*sqrt(npairs[2][iptbin][0]*npairs[2][iptbin][1])));
280
281
282 for(int ii=0;ii<hdiv[iptbin]->GetNbinsX(); ii++){
283 double d_val1 = 2*sqrt(hmasspt_ls[0][iptbin]->GetBinContent(ii+1)*hmasspt_ls[1][iptbin]->GetBinContent(ii+1));
284 if(d_val1>0){
285 double d_eval1 = d_val1/2*sqrt(pow(hmasspt_ls[0][iptbin]->GetBinError(ii+1)/hmasspt_ls[0][iptbin]->GetBinContent(ii+1),2)+
286 pow(hmasspt_ls[1][iptbin]->GetBinError(ii+1)/hmasspt_ls[1][iptbin]->GetBinContent(ii+1),2));
287 // double accep = hdiv[iptbin]->GetBinContent(ii+1)-hdiv[iptbin]->GetBinError(ii+1);
288 double accep = hdiv[iptbin]->GetBinContent(ii+1);
289 hmasspt[2][iptbin]->SetBinContent(ii+1, d_val1*accep);
290 hmasspt[2][iptbin]->SetBinError(ii+1, d_eval1*accep);
291 }
292 }
293
294 /*
295 double norm1=0;
296 double norm2=0;
297 for(int ii=0;ii<hmasspt[2][iptbin]->GetNbinsX(); ii++){
298 double d_bg1 = hmasspt[1][iptbin]->GetBinContent(ii+1);
299 double d_bg_pp = hmasspt_ls[2][iptbin]->GetBinContent(ii+1);
300 double d_bg_mm = hmasspt_ls[3][iptbin]->GetBinContent(ii+1);
301 double d_fg_pp = hmasspt_ls[0][iptbin]->GetBinContent(ii+1);
302 double d_fg_mm = hmasspt_ls[1][iptbin]->GetBinContent(ii+1);
303 double d_e_fg_pp = hmasspt_ls[0][iptbin]->GetBinError(ii+1);
304 double d_e_fg_mm = hmasspt_ls[1][iptbin]->GetBinError(ii+1);
305
306 if(d_fg_pp>0 && d_fg_mm>0 && d_bg_pp>0 && d_bg_mm>0){
307 //double d_val = 2*sqrt(d_fg_pp*d_fg_mm)*d_bg1/(2*sqrt(d_bg_pp*d_bg_mm));
308 //double d_e_val = sqrt(d_fg_pp*d_fg_mm)*sqrt(pow(d_e_fg_pp/d_fg_pp,2)+pow(d_e_fg_mm/d_fg_mm,2))*d_bg1/(2*sqrt(d_bg_pp*d_bg_mm));
309 double d_val = 2*(d_fg_pp*(0.5*d_bg1/d_bg_pp)*sqrt(npairs[2][iptbin][0])+d_fg_mm*(0.5*d_bg1/d_bg_mm)*sqrt(npairs[2][iptbin][1]))/(sqrt(npairs[2][iptbin][0])+sqrt(npairs[2][iptbin][1]));
310 double d_e_val = d_e_fg_pp*(0.5*d_bg1/d_bg_pp)+d_e_fg_mm*(0.5*d_bg1/d_bg_mm);
311
312 hmasspt[2][iptbin]->SetBinContent(ii+1, d_val);
313 hmasspt[2][iptbin]->SetBinError(ii+1, d_e_val);
314
315 }
316 }
317 */
318 //hmasspt[2][iptbin]->Scale(hmasspt[0][iptbin]->Integral()/(2*sqrt(norm1*norm2)));
319 cout<<" entry ="<<hmasspt[0][iptbin]->Integral()<<" "<<hmasspt[2][iptbin]->Integral()<<" "<<hmasspt_ls[0][iptbin]->Integral()<<" "<<hmasspt_ls[1][iptbin]->Integral()<<endl;
320 c1->cd(iptbin+6);
321 hmasspt[0][iptbin]->Draw();
322 hmasspt[2][iptbin]->SetLineColor(2);
323 hmasspt[2][iptbin]->Draw("same");
324 }
325
326
327 const int nmassbin=8;
328 float mass_low[nmassbin] ={0, 0.05, 0.09, 0.14, 0.20, 0.30, 0.1, 0.09};
329 float mass_high[nmassbin]={0.03, 0.09, 0.14, 0.20, 0.30, 0.40, 0.3, 0.30};
330 float nsig[5][nmassbin];
331 float nbg[5][nmassbin];
332 float ensig[5][nmassbin];
333 float enbg[5][nmassbin];
334
335 for(int ii=0;ii<nptbin; ii++){
336 for(int j=0;j<nmassbin; j++){
337 int binlow = hmasspt[0][ii]->FindBin(mass_low[j]);
338 int binhigh = hmasspt[0][ii]->FindBin(mass_high[j]);
339 nsig[ii][j] = hmasspt[0][ii]->Integral(binlow, binhigh-1);
340 nbg[ii][j] = hmasspt[2][ii]->Integral(binlow, binhigh-1);
341 ensig[ii][j]=0;
342 enbg[ii][j]=0;
343 for(int ib=binlow; ib<binhigh;ib++){
344 ensig[ii][j] += pow(hmasspt[0][ii]->GetBinError(ib),2);
345 enbg[ii][j] += pow(hmasspt[2][ii]->GetBinError(ib),2);
346 }
347
348 ensig[ii][j] = sqrt(ensig[ii][j]);
349 enbg[ii][j] = sqrt(enbg[ii][j]);
350
351
352 float r1 = nsig[ii][j]-nbg[ii][j];
353 float r0 = nsig[ii][0]-nbg[ii][0];
354 float er1 = sqrt(pow(ensig[ii][j],2)+pow(enbg[ii][j],2));
355 float er0 = sqrt(pow(ensig[ii][0],2)+pow(enbg[ii][0],2));
356
357
358 if(r0>0 && r1>0){
359 cout<<" pt : "<<ii<<" mass "<<mass_low[j]<<"-"<<mass_high[j]<<" --> Nsig = "<<nsig[ii][j]<<", Nbg = "<<nbg[ii][j]
360 <<" : Nsig-Nbg "<<nsig[ii][j]-nbg[ii][j]<<" +- "<<sqrt(pow(ensig[ii][j],2)+pow(enbg[ii][j],2))<<" ("<<ensig[ii][j]<<" "<<enbg[ii][j]<<")"
361 <<" : Ratio = "<<r1/r0<<" +- "<<r1/r0*sqrt(pow(er1/r1,2)+pow(er0/r0,2))
362 <<endl;
363 }else{
364 cout<<" pt : "<<ii<<" mass "<<mass_low[j]<<"-"<<mass_high[j]<<" --> Nsig = "<<nsig[ii][j]<<", Nbg = "<<nbg[ii][j]
365 <<" : Nsig-Nbg "<<nsig[ii][j]-nbg[ii][j]<<" +- "<<sqrt(pow(ensig[ii][j],2)+pow(enbg[ii][j],2))<<" ("<<ensig[ii][j]<<" "<<enbg[ii][j]<<")"
366 <<endl;
367 }
368 }
369 }
370
371
372 /// for histogram
373 TH1D *hsub[5];
374 const int nbins = 55;
375 double x[nbins+1]={0, 0.01, 0.02, 0.03, 0.05, 0.09, 0.14, 0.20, 0.30, 0.40,
376 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4,
377 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4,
378 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4,
379 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4,
380 4.5, 4.6, 4.7, 4.8, 4.9, 5.0};
381
382// const int nbins = 53;
383// double x[nbins+1]={0, 0.01, 0.02, 0.03, 0.05, 0.1, 0.20, 0.40,
384// 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4,
385// 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4,
386// 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4,
387// 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4,
388// 4.5, 4.6, 4.7, 4.8, 4.9, 5.0};
389
390 int narrows[nptbin];
391 TArrow *Arrow[nptbin][10];
392
393 for(int i=0; i<nptbin; i++){
394 narrows[i]=0;
395
396
397 sprintf(name,"hsub_mass_pt%d", i);
398 hsub[i] = new TH1D(name, name, nbins, x);
399
400 for(int j=0; j<nbins; j++){
401 int binl = hmasspt[0][i]->FindBin(x[j]);
402 int binh = hmasspt[0][i]->FindBin(x[j+1]);
403 float dw = x[j+1]-x[j];
404
405 float ent = 0;
406 float e_ent = 0;
407 float sig = 0;
408 float bg = 0;
409 for(int k=binl; k<binh; k++){
410 sig += hmasspt[0][i]->GetBinContent(k);
411 bg += hmasspt[2][i]->GetBinContent(k);
412 ent += (hmasspt[0][i]->GetBinContent(k)-hmasspt[2][i]->GetBinContent(k));
413 e_ent += (pow(hmasspt[0][i]->GetBinError(k),2)+
414 pow(hmasspt[2][i]->GetBinError(k),2));
415 }
416 e_ent = sqrt(e_ent);
417 hsub[i]->SetBinContent(j+1, ent/dw);
418 hsub[i]->SetBinError(j+1, e_ent/dw);
419
420 if(sig<bg && sig>0 && x[j]<0.5){ // calculate 90% confidence level
421 double conf90=0;
422 sig = sig*nevents;
423 bg = bg*nevents;
424 TH1F *htmp = new TH1F("htmp","htmp",1000,-2*bg, 2*bg);
425 for(int kk=0;kk<100000;kk++){
426 double a1 = gRandom->PoissonD(sig);
427 double a2 = gRandom->PoissonD(bg);
428 htmp->Fill(a1-a2);
429 }
430 int bin0 = htmp->FindBin(0);
431 float ent0 = htmp->Integral(bin0, 1000);
432 htmp->Scale(1.0/ent0);
433 //htmp->Scale(1.0/100000);
434 int bincent=0;
435 for(int ib=bin0;ib<1000;ib++){
436 if(htmp->Integral(bin0, ib)>0.9){
437 bincent = ib;
438 break;
439 }
440 }
441 conf90 = htmp->GetBinCenter(bincent);
442 cout<<i<<" "<<hsub[i]->GetBinCenter(j+1)<<" : 90% conf: "<<conf90<<" from sig :"<<sig<<" bg:"<<bg<<endl;
443 Arrow[i][narrows[i]] = new TArrow(hsub[i]->GetBinCenter(j+1), 1.0e-07, hsub[i]->GetBinCenter(j+1), conf90/(dw*nevents),0.01,"<|");
444 narrows[i]++;
445 htmp->Clear();
446 }// calculate 90% confidence level
447
448 }
449 }
450
451 TH1D *hsub_clone[5];
452
453 TCanvas *c3 = new TCanvas("c3","c3",1200,600);
454 c3->Divide(4,1);
455
456 c3->cd(1);
457 gPad->SetGrid(1); gPad->SetLogy();
458 TH2F *hw = new TH2F("hw","hw",10,0,0.5,10,1.0e-07, 10);
459 hw->Draw();
460 hw->SetXTitle("M_{ee} [GeV]");
461 hw->SetYTitle("dN_{ee}/dM_{ee}");
462 c3->cd(2);
463 gPad->SetGrid(1); gPad->SetLogy();
464 hw->Draw();
465 c3->cd(3);
466 gPad->SetGrid(1); gPad->SetLogy();
467 hw->Draw();
468 c3->cd(4);
469 gPad->SetGrid(1); gPad->SetLogy();
470 hw->Draw();
471
472 for(int i=0; i<nptbin; i++){
473 hsub_clone[i] = (TH1D*)hsub[i]->Clone();
474 sprintf(name,"%s_cln",hsub[i]->GetName());
475 hsub_clone[i]->SetName(name);
476 hsub_clone[i]->SetLineColor(i+1);
477 hsub_clone[i]->SetMarkerColor(i+1);
478 hsub_clone[i]->SetMarkerStyle(20);
479 }
480
481
482 // int bin = hsub_clone[1]->FindBin(0.03);
483 // float ent = hsub_clone[1]->Integral(0,bin);
484 for(int i=0; i<nptbin; i++){
485 //float ent1 = hsub_clone[i]->Integral(0,bin);
486 //hsub_clone[i]->Scale(ent/ent1);
487 c3->cd(i+1);
488 hsub_clone[i]->Draw("same");
489 for(int k=0;k<narrows[i];k++){
490 Arrow[i][k]->SetLineColor(i+1);
491 Arrow[i][k]->Draw();
492 }
493 }
494
495 TCanvas *c5 = new TCanvas("c5","c5",1200,600);
496 c5->Divide(3,1);
497
498 c5->cd(1);
499 gPad->SetGrid(1); gPad->SetLogy();
500 TH2F *hw2 = new TH2F("hw2","hw2",10,0,0.5,10,1.0e-07, 10);
501 hw2->Draw();
502 hw2->SetXTitle("M_{ee} [GeV]");
503 hw2->SetYTitle("1/N_{evt}dN_{ee}/dM_{ee}");
504 c5->cd(2);
505 gPad->SetGrid(1); gPad->SetLogy();
506 hw2->Draw();
507 c5->cd(3);
508 gPad->SetGrid(1); gPad->SetLogy();
509 hw2->Draw();
510
511 TH1D *hmass_bin[2][nptbin];
512 for(int i=1; i<nptbin; i++){
513 c5->cd(i);
514 hmass_bin[0][i] = (TH1D*)hmasspt[0][i]->Clone();
515 hmass_bin[1][i] = (TH1D*)hmasspt[2][i]->Clone();
516
517
518 float binw = hmass_bin[0][i]->GetBinWidth(10);
519 hmass_bin[0][i]->Scale(1/binw);
520 hmass_bin[1][i]->Scale(1/binw);
521
522 hmass_bin[1][i]->SetLineColor(2);
523 hmass_bin[0][i]->Draw("same");
524 hmass_bin[1][i]->Draw("same");
525 hsub_clone[i]->Draw("same");
526 for(int k=0;k<narrows[i];k++){
527 Arrow[i][k]->SetLineColor(i+1);
528 Arrow[i][k]->Draw();
529 }
530 }
531
532
533
534
535
536
537 // TLegend *leg = new TLegend(0.5,0.7,0.875,0.875);
538 // leg->SetFillColor(10);
539 // if(nptbin==4){
540 // leg->AddEntry(hsub_clone[0],"0.5<p_{T}<1.0 GeV","pl");
541 // leg->AddEntry(hsub_clone[1],"1.0<p_{T}<1.5 GeV","pl");
542 // leg->AddEntry(hsub_clone[2],"1.5<p_{T}<2.0 GeV","pl");
543 // leg->AddEntry(hsub_clone[3],"2.0<p_{T}","pl");
544 // }else{
545 // }
546 // leg->Draw("same");
547
548 /////draw the results
549 TH1D *hmass_cocktail[4];
550 TH1D *hmass_cocktail_comp[4][6];
551 TFile *fin1 = new TFile("../../ana2/cocktail/cocktail_allpt_pt04_no_smearing.root");
552 fin1->cd();
553 hmass_cocktail[0] = (TH1D*)fin1->Get("hmass_all");
554 hmass_cocktail[0]->SetName("hmass_all_pt0");
555 hmass_cocktail_comp[0][0] = (TH1D*)fin1->Get("hmass_3");
556 hmass_cocktail_comp[0][0]->SetName("hmass_pt0_eta_prime");
557 hmass_cocktail_comp[0][1] = (TH1D*)fin1->Get("hmass_4");
558 hmass_cocktail_comp[0][1]->SetName("hmass_pt0_rho");
559 hmass_cocktail_comp[0][2] = (TH1D*)fin1->Get("hmass_7");
560 hmass_cocktail_comp[0][2]->SetName("hmass_pt0_pi0");
561 hmass_cocktail_comp[0][3] = (TH1D*)fin1->Get("hmass_8");
562 hmass_cocktail_comp[0][3]->SetName("hmass_pt0_eta");
563 hmass_cocktail_comp[0][4] = (TH1D*)fin1->Get("hmass_9");
564 hmass_cocktail_comp[0][4]->SetName("hmass_pt0_omega");
565 hmass_cocktail_comp[0][5] = (TH1D*)fin1->Get("hmass_10");
566 hmass_cocktail_comp[0][5]->SetName("hmass_pt0_phi");
567
568
569 //TFile *fin2 = new TFile("./cocktail//cocktail_pt10_15.root");
570 TFile *fin2 = new TFile("../../ana2/cocktail//cocktail_pt04_pairpt_10_15_no_smearing.root");
571 fin2->cd();
572 hmass_cocktail[1] = (TH1D*)fin2->Get("hmass_all");
573 hmass_cocktail[1]->SetName("hmass_all_pt1");
574 hmass_cocktail_comp[1][0] = (TH1D*)fin1->Get("hmass_3");
575 hmass_cocktail_comp[1][0]->SetName("hmass_pt1_eta_prime");
576 hmass_cocktail_comp[1][1] = (TH1D*)fin1->Get("hmass_4");
577 hmass_cocktail_comp[1][1]->SetName("hmass_pt1_rho");
578 hmass_cocktail_comp[1][2] = (TH1D*)fin1->Get("hmass_7");
579 hmass_cocktail_comp[1][2]->SetName("hmass_pt1_pi0");
580 hmass_cocktail_comp[1][3] = (TH1D*)fin1->Get("hmass_8");
581 hmass_cocktail_comp[1][3]->SetName("hmass_pt1_eta");
582 hmass_cocktail_comp[1][4] = (TH1D*)fin1->Get("hmass_9");
583 hmass_cocktail_comp[1][4]->SetName("hmass_pt1_omega");
584 hmass_cocktail_comp[1][5] = (TH1D*)fin1->Get("hmass_10");
585 hmass_cocktail_comp[1][5]->SetName("hmass_pt1_phi");
586
587
588
589 //TFile *fin3 = new TFile("./cocktail//cocktail_pt15_20.root");
590 //TFile *fin3 = new TFile("./cocktail//cocktail_pt15_20_no_smearing.root");
591 TFile *fin3 = new TFile("../../ana2/cocktail//cocktail_pt04_pairpt_15_20_no_smearing.root");
592 fin3->cd();
593 hmass_cocktail[2] = (TH1D*)fin3->Get("hmass_all");
594 hmass_cocktail[2]->SetName("hmass_all_pt2");
595 hmass_cocktail_comp[2][0] = (TH1D*)fin1->Get("hmass_3");
596 hmass_cocktail_comp[2][0]->SetName("hmass_pt2_eta_prime");
597 hmass_cocktail_comp[2][1] = (TH1D*)fin1->Get("hmass_4");
598 hmass_cocktail_comp[2][1]->SetName("hmass_pt2_rho");
599 hmass_cocktail_comp[2][2] = (TH1D*)fin1->Get("hmass_7");
600 hmass_cocktail_comp[2][2]->SetName("hmass_pt2_pi0");
601 hmass_cocktail_comp[2][3] = (TH1D*)fin1->Get("hmass_8");
602 hmass_cocktail_comp[2][3]->SetName("hmass_pt2_eta");
603 hmass_cocktail_comp[2][4] = (TH1D*)fin1->Get("hmass_9");
604 hmass_cocktail_comp[2][4]->SetName("hmass_pt2_omega");
605 hmass_cocktail_comp[2][5] = (TH1D*)fin1->Get("hmass_10");
606 hmass_cocktail_comp[2][5]->SetName("hmass_pt2_phi");
607
608
609
610 //TFile *fin4 = new TFile("./cocktail//cocktail_pt20_no_smearing.root");
611 //TFile *fin4 = new TFile("./cocktail//cocktail_pt20.root");
612TFile *fin4 = new TFile("../../ana2/cocktail//cocktail_pt04_pairpt_20_no_smearing.root");
613 fin4->cd();
614 hmass_cocktail[3] = (TH1D*)fin4->Get("hmass_all");
615 hmass_cocktail[3]->SetName("hmass_all_pt3");
616 hmass_cocktail_comp[3][0] = (TH1D*)fin1->Get("hmass_3");
617 hmass_cocktail_comp[3][0]->SetName("hmass_pt3_eta_prime");
618 hmass_cocktail_comp[3][1] = (TH1D*)fin1->Get("hmass_4");
619 hmass_cocktail_comp[3][1]->SetName("hmass_pt3_rho");
620 hmass_cocktail_comp[3][2] = (TH1D*)fin1->Get("hmass_7");
621 hmass_cocktail_comp[3][2]->SetName("hmass_pt3_pi0");
622 hmass_cocktail_comp[3][3] = (TH1D*)fin1->Get("hmass_8");
623 hmass_cocktail_comp[3][3]->SetName("hmass_pt3_eta");
624 hmass_cocktail_comp[3][4] = (TH1D*)fin1->Get("hmass_9");
625 hmass_cocktail_comp[3][4]->SetName("hmass_pt3_omega");
626 hmass_cocktail_comp[3][5] = (TH1D*)fin1->Get("hmass_10");
627 hmass_cocktail_comp[3][5]->SetName("hmass_pt3_phi");
628
629
630
631 for(int i=0;i<4;i++){
632 hmass_cocktail[i]->Rebin(4);
633 float ent = 0;
634 for(int j=0;j<hsub_clone[i]->GetNbinsX();j++){
635 if(hsub_clone[i]->GetBinCenter(j+1)<0.03){
636 ent += hsub_clone[i]->GetBinWidth(j+1)*hsub_clone[i]->GetBinContent(j+1);
637 }
638 }
639 float ent2=0;
640 for(int j=0;j<hmass_cocktail[i]->GetNbinsX();j++){
641 if(hmass_cocktail[i]->GetBinCenter(j+1)<0.03){
642 ent2 += hmass_cocktail[i]->GetBinWidth(j+1)*hmass_cocktail[i]->GetBinContent(j+1);
643 }
644 }
645 hmass_cocktail[i]->Scale(ent/ent2);
646 hmass_cocktail[i]->SetLineColor(6);
647 hmass_cocktail[i]->SetLineWidth(2);
648 c3->cd(i+1);
649 hmass_cocktail[i]->Draw("same,l");
650 if(i>=1){
651 c5->cd(i);
652 hmass_cocktail[i]->Draw("same,l");
653 }
654
655 for(int j=0;j<6;j++){
656 hmass_cocktail_comp[i][j]->Rebin(4);
657 hmass_cocktail_comp[i][j]->Scale(ent/ent2*0.2);
658 hmass_cocktail_comp[i][j]->SetLineColor(1);
659 //hmass_cocktail_comp[i][j]->Draw("same");
660 }
661 }
662
663 float r_cocktail[nptbin];
664 cout<<" ***************** "<<endl;
665 cout<<" mass ratio ("<<mass_low[nmassbin-1]<<"-"<<mass_high[nmassbin-1]<<"/0-30) in cocktail"<<endl;
666 TH1F *hmass_r_ck = new TH1F("hmass_r_cocktail","mass ratio cocktail", nptbin+1, mass_r_ptbin);
667 for(int i=1;i<nptbin;i++){
668 int bin1=hmass_cocktail[i]->FindBin(0);
669 int bin2=hmass_cocktail[i]->FindBin(0.03);
670 float ent1 = hmass_cocktail[i]->Integral(bin1, bin2-1);
671 bin1=hmass_cocktail[i]->FindBin(mass_low[nmassbin-1]);
672 bin2=hmass_cocktail[i]->FindBin(mass_high[nmassbin-1]);
673 float ent2 = hmass_cocktail[i]->Integral(bin1, bin2-1);
674 cout<<" pt : "<<i<<" "
675 <<ent2/ent1
676 <<endl;
677 int bin = hmass_r_ck->FindBin((0.5*(pt_low[i]+pt_high[i])));
678 hmass_r_ck->SetBinContent(bin, ent2/ent1);
679 hmass_r_ck->SetBinError(bin, ent2/ent1*0.05);
680 r_cocktail[i] = ent2/ent1;
681 }
682
683 float r_dir = 0.37; // this is from 90-300MeV
684
685 TH1F *hmass_r = new TH1F("hmass_r","mass ratio", nptbin+1, mass_r_ptbin);
686 TH1F *hrgamma = new TH1F("hrgamma","r gmamma", nptbin+1, mass_r_ptbin);
687 cout<<" ***************** "<<endl;
688 cout<<" mass ratio (100-300/0-30)"<<endl;
689 for(int i=1; i<nptbin; i++){
690 cout<<" pt : "<<i<<" "
691 <<nsig[i][nmassbin-1]-nbg[i][nmassbin-1]<<" +- "<<sqrt(pow(ensig[i][nmassbin-1],2)+pow(enbg[i][nmassbin-1],2))<<" / "
692 <<(nsig[i][0]-nbg[i][0])<<" ratio = "
693 <<(nsig[i][nmassbin-1]-nbg[i][nmassbin-1])/(nsig[i][0]-nbg[i][0])
694 <<" +- "<<sqrt(pow(ensig[i][nmassbin-1],2)+pow(enbg[i][nmassbin-1],2))/(nsig[i][0]-nbg[i][0])
695 <<endl;
696
697 int bin = hmass_r->FindBin((0.5*(pt_low[i]+pt_high[i])));
698 float r_data = (nsig[i][nmassbin-1]-nbg[i][nmassbin-1])/(nsig[i][0]-nbg[i][0]);
699 float er_data = sqrt(pow(ensig[i][nmassbin-1],2)+pow(enbg[i][nmassbin-1],2))/(nsig[i][0]-nbg[i][0]);
700
701 hmass_r->SetBinContent(bin, r_data);
702 hmass_r->SetBinError(bin, er_data);
703
704 float r_gamma = (r_data-r_cocktail[i])/(r_dir - r_cocktail[i]);
705 float e_r_gamma = (er_data)/(r_dir - r_cocktail[i]);
706 hrgamma->SetBinContent(bin, r_gamma);
707 hrgamma->SetBinError(bin, e_r_gamma);
708
709 }
710
711 hrgamma->SetXTitle("p_{T} [GeV]");
712 hrgamma->SetYTitle("R_{#gamma}^{0-30}");
713 hrgamma->SetMarkerColor(icent+1);
714 hrgamma->SetLineColor(icent+1);
715 hrgamma->SetMarkerStyle(20);
716 hrgamma->SetMinimum(0);
717 hrgamma->SetMaximum(1);
718
719 TCanvas *c4 = new TCanvas("c4","c4");
720 gPad->SetGrid(1);
721 hmass_r->SetXTitle("p_{T} [GeV]");
722 hmass_r->SetYTitle("R_{data}^{90-300}");
723 hmass_r->SetMarkerColor(icent+1);
724 hmass_r->SetLineColor(icent+1);
725 hmass_r->SetMarkerStyle(20);
726 hmass_r->SetMinimum(0);
727 hmass_r->SetMaximum(0.5);
728
729 hmass_r_ck->SetMarkerStyle(24);
730 hmass_r->Draw(); hmass_r_ck->Draw("same");
731
732
733
734
735
736
737 TFile *fout = new TFile(outname,"recreate");
738 hw->Write();
739 for(int i=0; i<nptbin; i++){
740 hsub[i]->Write();
741 hsub_clone[i]->Write();
742 }
743 for(int j=0;j<4;j++){
744 hmass_cocktail[j]->Write();
745 for(int k=0;k<6;k++){
746 hmass_cocktail_comp[j][k]->Write();
747 }
748 }
749 hmass_r->Write();
750 hmass_r_ck->Write();
751 hrgamma->Write();
752
753}
754
755
756
757
758