3 const Double_t kMean=0.135 ; //Approximate peak position to facilitate error estimate
4 Int_t effcor=1; //correct on efficiency
6 Double_t PeakPosition(Double_t pt){
7 //Fit to the measured peak position
8 return 4.99292e-003*exp(-pt*9.32300e-001)+1.34944e-001 ;
10 Double_t PeakWidth(Double_t pt){
11 //fit to the measured peak width
15 return TMath::Sqrt(a*a+b*b/pt/pt+c*c*pt*pt) ;
18 Double_t CB(Double_t * x, Double_t * par){
19 //Parameterization of Real/Mixed ratio
22 Double_t dx=(x[0]-m)/s ;
23 return par[0]*exp(-dx*dx/2.)+par[3] + par[4]*(x[0]-kMean) /*+ par[5]*(x[0]-kMean)*(x[0]-kMean) */;
25 Double_t CB2(Double_t * x, Double_t * par){
26 //Another parameterization of Real/Mixed ratio
29 Double_t dx=(x[0]-m)/s ;
30 return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean) + par[5]*(x[0]-kMean)*(x[0]-kMean) /*+ par[6]*(x[0]-kMean)*(x[0]-kMean)*(x[0]-kMean)*/;
32 Double_t CBs(Double_t * x, Double_t * par){
33 //Parameterizatin of signal
36 Double_t dx=(x[0]-m)/s ;
37 return par[0]*exp(-dx*dx/2.)/TMath::Sqrt(TMath::TwoPi())/s+par[3] ;
39 Double_t BG1(Double_t * x, Double_t * par){
40 //Normalizatino of Mixed
41 return par[0] + par[1]*(x[0]-kMean) /*+ par[2]*(x[0]-kMean)*(x[0]-kMean)*/;
43 Double_t BG2(Double_t * x, Double_t * par){
44 //Another normalization of Mixed
45 return par[0]+par[1]*(x[0]-kMean) + par[2]*(x[0]-kMean)*(x[0]-kMean) /*+ par[3]*(x[0]-kMean)*(x[0]-kMean)*(x[0]-kMean)*/;
48 void Method1(Int_t cen=1,Int_t side=1){
49 //side = 1 - VOC, side = 0 - V0A
50 //side = 2 - TPC 3 sub method, 3 - TPC 2 sub method
52 //Fit Real/Mixed ratio, normalize Mixed and subtract it from Real
54 // TH1F * hCen = GetCen() ;
55 // TH2F * hCenPHOS = GetCenPHOS() ;
59 sprintf(w,""); //Weight: no weight - "NW", with weight - ""
62 sprintf(PID,"Both2core") ;
64 char kind[15],kind2[15] ; //Reaction plane
65 char d[15]; //detector
66 sprintf(kind,"") ; //"" for real
67 sprintf(kind2,"") ; //"" for mixed
69 if(side==1)sprintf(d,"V0C");
70 else if(side==0)sprintf(d,"V0A");
71 else sprintf(d,"TPC");
74 cout<<"Getting histos from data..."<<endl;
75 h3DR = (TH3D*)GetRealMixed(cen,w,d,kind,kind2,PID,1);
76 h3DM = (TH3D*)GetRealMixed(cen,w,d,kind,kind2,PID,0);
78 TH1F * hCen = GetCen() ;
79 TH2F * hCenPHOS = GetCenPHOS() ;
80 cout<<"Finished"<<endl;
82 // sprintf(kind,"") ; //"" for all
83 const Double_t nphi=5;
85 for(Int_t i=0;i<=nphi;i++){
86 xphi[i]=TMath::PiOver2()*i/nphi ;
87 // cout<<xphi[i]<<" ";
92 sprintf(inname,"../flow11h_Apr14.root");
94 TFile * f = new TFile(inname) ;
96 TH2F * hcos2AC = (TH2F*)f->Get("cos2V0AC");
97 TH2F * hcos2FAC= (TH2F*)f->Get("cos2FAC");
98 Double_t resAC = GetCos(hcos2AC,cen);
101 hcos2AC = (TH2F*)f->Get("cos2V0ATPC");
102 Double_t resAT = GetCos(hcos2AC,cen);
105 hcos2AC = (TH2F*)f->Get("cos2V0CTPC");
106 Double_t resCT = GetCos(hcos2AC,cen);
109 hcos2AC = (TH2F*)f->Get("cos2AC");
110 Double_t resT = GetRes(hcos2AC,cen);
112 //side = 1 - VOC, side = 0 - V0A
113 //side = 2 - TPC 3 sub method, 3 - TPC 2 sub method
116 if(side==1) res = TMath::Sqrt(resAC*resCT/resAT);
117 else if(side==0)res = TMath::Sqrt(resAC*resAT/resCT);
118 else if(side==2)res = TMath::Sqrt(resAT*resCT/resAC);
123 cout<<"RP resolution = "<< res <<endl;
130 sprintf(key,"PID%s_cen%d",kind,cen) ;
131 sprintf(key2,"%s_mr1",key) ;
132 TH1D * mr1 = new TH1D(key2,"Mass",nbin,xa) ;
133 sprintf(key2,"%s_sr1",key) ;
134 TH1D * sr1 = new TH1D(key2,"Width",nbin,xa) ;
135 sprintf(key2,"%s_ar1",key) ;
136 TH2D * ar1 = new TH2D(key2,"a",nbin,xa,nphi,xphi) ;
137 sprintf(key2,"%s_br1",key) ;
138 TH2D * br1 = new TH2D(key2,"b",nbin,xa,nphi,xphi) ;
139 sprintf(key2,"%s_yr1",key) ;
140 TH2D * nr1 = new TH2D(key2,"Raw yield",nbin,xa,nphi,xphi) ;
141 sprintf(key2,"%s_yr1int",key) ;
142 TH2D * nr1int = new TH2D(key2,"Raw yield, integrated",nbin,xa,nphi,xphi) ;
145 sprintf(key2,"%s_mr2",key) ;
146 TH1D * mr2 = new TH1D(key2,"Mass",nbin,xa) ;
147 sprintf(key2,"%s_sr2",key) ;
148 TH1D * sr2 = new TH1D(key2,"Width",nbin,xa) ;
149 sprintf(key2,"%s_ar2",key) ;
150 TH2D * ar2 = new TH2D(key2,"a",nbin,xa,nphi,xphi) ;
151 sprintf(key2,"%s_br2",key) ;
152 TH2D * br2 = new TH2D(key2,"b",nbin,xa,nphi,xphi) ;
153 sprintf(key2,"%s_cr2",key) ;
154 TH2D * cr2 = new TH2D(key2,"c",nbin,xa,nphi,xphi) ;
155 sprintf(key2,"%s_yr2",key) ;
156 TH2D * nr2 = new TH2D(key2,"Raw yield",nbin,xa,nphi,xphi) ;
157 sprintf(key2,"%s_yr2int",key) ;
159 TH2D * nr2int = new TH2D(key2,"Raw yield, integrated",nbin,xa,nphi,xphi) ;
161 TF1 * fit1 = new TF1("fit1",CB,0.,1.,5) ;
162 fit1->SetParName(0,"A") ;
163 fit1->SetParName(1,"m_{0}") ;
164 fit1->SetParName(2,"#sigma") ;
165 fit1->SetParName(3,"a_{0}") ;
166 fit1->SetParName(4,"a_{1}") ;
167 fit1->SetLineWidth(2) ;
168 fit1->SetLineColor(2) ;
169 TF1 * fgs = new TF1("gs",CBs,0.,1.,4) ;
170 fgs->SetLineColor(2) ;
171 fgs->SetLineWidth(1) ;
173 TF1 * fit2 = new TF1("fit2",CB2,0.,1.,6) ;
174 fit2->SetParName(0,"A") ;
175 fit2->SetParName(1,"m_{0}") ;
176 fit2->SetParName(2,"#sigma") ;
177 fit2->SetParName(3,"a_{0}") ;
178 fit2->SetParName(4,"a_{1}") ;
179 fit2->SetParName(5,"a_{2}") ;
181 fit2->SetLineWidth(2) ;
182 fit2->SetLineColor(4) ;
183 fit2->SetLineStyle(2) ;
184 TF1 * fbg1 = new TF1("bg1",BG1,0.,1.,2) ;
185 TF1 * fbg2 = new TF1("bg2",BG2,0.,1.,3) ;
187 //calculate average mass and width for all pt bins
188 TCanvas * cav = new TCanvas("Average","Average",10,10,800,750) ;
191 cout<<"Finished"<<endl;
193 cout<<"Do avarage minv... ";
194 for(Int_t i=1;i<=nbin;i++){
195 h3DR->GetYaxis()->SetRangeUser(xa[i-1],xa[i]) ;
196 h3DM->GetYaxis()->SetRangeUser(xa[i-1],xa[i]) ;
197 Double_t pt=(xa[i]+xa[i-1])/2. ;
198 TH2D * hpav = (TH2D*)h3DR->Project3D("x");
199 TH2D * hpmav= (TH2D*)h3DM->Project3D("x") ;
210 for(Int_t ib=1; ib<=hpav->GetNbinsX();ib++){if(hpav->GetBinContent(ib)==0)hpav->SetBinError(ib,1.);}
211 for(Int_t ib=1; ib<=hpav->GetNbinsX();ib++){if(hpmav->GetBinContent(ib)==0)hpmav->SetBinError(ib,1.);}
212 TH1D * hpm2av = (TH1D*)hpmav->Clone("Bg1av") ;
213 TH1D * hpcopyav = (TH1D*)hpav->Clone("hpcopyav") ;
214 TH1D * hp2av = (TH1D*)hpav->Clone("hp2av") ;
215 TH1D * hpm3av = (TH1D*)hpmav->Clone("Bg2av") ;
216 TH1D * hp3av = (TH1D*)hpav->Clone("hp3av") ;
218 hpcopyav->Divide(hpmav) ;
219 sprintf(key,"%3.1f<p_{t}<%3.1f GeV/c",xa[i-1],xa[i]) ;
220 hpcopyav->SetTitle(key) ;
221 hpcopyav->SetMarkerStyle(20) ;
222 hpcopyav->SetMarkerSize(0.7) ;
224 fit1->SetParameters(TMath::Min(0.3,0.001*i),0.14,0.01,0.008,-0.0002,0) ;
226 fit1->SetParLimits(2,0.003,0.03) ; //peak width
227 fit1->SetParLimits(1,0.1,0.2); //peak mean
228 fit1->SetParLimits(0,0.,10.) ; //peak height
231 Double_t rangeMin=0.1 ;
232 Double_t rangeMax=0.25 ; //TMath::Min(0.15+i*0.05,0.3) ;
233 // if(cen==0)rangeMax=0.30; //TMath::Min(0.15+i*0.05,0.3) ;
234 hpcopyav->Fit("fit1","Q","",rangeMin,rangeMax) ;
235 hpcopyav->Fit("fit1","MQ","",rangeMin,rangeMax) ;
237 mr1->SetBinContent(i,fit1->GetParameter(1)) ;
238 mr1->SetBinError(i,fit1->GetParError(1)) ;
239 sr1->SetBinContent(i,TMath::Abs(fit1->GetParameter(2))) ;
240 sr1->SetBinError(i,fit1->GetParError(2)) ;
242 fit2->SetParameters(fit1->GetParameters()) ;
243 fit2->SetParameter(5,0) ;
244 hpcopyav->Fit("fit2","+Q","",rangeMin,rangeMax) ;
245 hpcopyav->Fit("fit2","+MQ","",rangeMin,rangeMax) ;
246 hpcopyav->GetXaxis()->SetRangeUser(0.05,0.3) ;
250 mr2->SetBinContent(i,fit2->GetParameter(1)) ;
251 mr2->SetBinError(i,fit2->GetParError(1)) ;
252 sr2->SetBinContent(i,TMath::Abs(fit2->GetParameter(2))) ;
253 sr2->SetBinError(i,fit2->GetParError(2)) ;
262 cout<<"Finished"<<endl;
264 h3DR->GetYaxis()->SetRangeUser(0.,19.) ;
265 h3DM->GetYaxis()->SetRangeUser(0.,19.) ;
267 cout<<"Do dNdPhi bins... ";
268 //!!!!!!!!!!!!!!!!!DNDPHI bins:
270 for(Int_t iphi=1;iphi<=nphi;iphi++){
272 h3DR->GetZaxis()->SetRange(iphi,iphi) ;
273 h3DM->GetZaxis()->SetRange(iphi,iphi) ;
275 TH2D * h = (TH2D*)h3DR->Project3D("yx");
276 TH2D * hm= (TH2D*)h3DM->Project3D("yx") ;
277 h->SetName("sliceRe") ;
278 hm->SetName("sliceMi") ;
280 h3DR->GetZaxis()->SetRange(11-iphi,11-iphi) ;
281 h3DM->GetZaxis()->SetRange(11-iphi,11-iphi) ;
283 TH2D * ha = (TH2D*)h3DR->Project3D("yx");
284 TH2D * hma= (TH2D*)h3DM->Project3D("yx") ;
286 h->Add(ha) ; delete ha ;
287 hm->Add(hma); delete hma ;
289 sprintf(key2,"mggFit%d_Signal",iphi) ;
290 TCanvas * c3 = new TCanvas(key2,key2,10,10,800,750) ;
293 sprintf(key2,"mggFit%d",iphi) ;
294 TCanvas * c1 = new TCanvas(key2,key2,10,10,800,750) ;
298 // TCanvas * c2=0,*c4=0,*c5=0,*c6=0 ;
299 TAxis * pta=h->GetYaxis() ;
300 TAxis * ma=h->GetXaxis() ;
303 for(Int_t i=1;i<=nbin;i++){
305 Int_t imin=pta->FindBin(xa[i-1]+0.0001);
306 Int_t imax=pta->FindBin(xa[i]-0.0001) ;
307 Double_t pt=(xa[i]+xa[i-1])/2. ;
308 TH1D * hp = h->ProjectionX("re",imin,imax) ;
310 TH1D * hpm= hm->ProjectionX("mi",imin,imax) ;
322 for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hp->GetBinContent(ib)==0)hp->SetBinError(ib,1.);}
323 for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hpm->GetBinContent(ib)==0)hpm->SetBinError(ib,1.);}
324 TH1D * hpm2 = (TH1D*)hpm->Clone("Bg1") ;
325 TH1D * hpcopy = (TH1D*)hp->Clone("hpcopy") ;
326 TH1D * hp2 = (TH1D*)hp->Clone("hp2") ;
327 hpcopy->Divide(hpm) ;
328 sprintf(key,"%3.1f<p_{t}<%3.1f GeV/c",xa[i-1],xa[i]) ;
329 hpcopy->SetTitle(key) ;
330 hpcopy->SetMarkerStyle(20) ;
331 hpcopy->SetMarkerSize(0.7) ;
333 fit1->SetParameters(TMath::Min(0.3,0.0001*i),0.135,0.008,0.008,-0.0002,0) ;
334 fit1->FixParameter(1,mr1->GetBinContent(i)) ;
335 fit1->FixParameter(2,sr1->GetBinContent(i)) ;
337 // fit1->SetParLimits(2,0.005,0.015) ;
338 // fit1->SetParLimits(0,0.,1.) ;
340 Double_t rangeMin=0.08 ;
341 Double_t rangeMax=0.25 ; //TMath::Min(0.15+i*0.05,0.3) ;
342 // if(cen==0)rangeMax=0.25 ;
344 hpcopy->Fit("fit1","NQ","",rangeMin,rangeMax) ;
345 hpcopy->Fit("fit1","MQ","",rangeMin,rangeMax) ;
347 ar1->SetBinContent(i,iphi,fit1->GetParameter(3)) ;
348 ar1->SetBinError(i,iphi,fit1->GetParError(3)) ;
349 br1->SetBinContent(i,iphi,fit1->GetParameter(4)) ;
350 br1->SetBinError(i,iphi,fit1->GetParError(4)) ;
352 fit2->SetParameters(fit1->GetParameters()) ;
353 fit2->FixParameter(1,mr2->GetBinContent(i)) ;
354 fit2->FixParameter(2,sr2->GetBinContent(i)) ;
355 // fit2->SetParLimits(2,0.005,0.015) ;
356 // fit2->SetParLimits(0,0.,1.) ;
357 fit2->SetParameter(5,0) ;
359 hpcopy->Fit("fit2","+NQ","",rangeMin,rangeMax) ;
360 hpcopy->Fit("fit2","+MQ","",rangeMin,rangeMax) ;
362 ar2->SetBinContent(i,iphi,fit2->GetParameter(3)) ;
363 ar2->SetBinError(i,iphi,fit2->GetParError(3)) ;
364 br2->SetBinContent(i,iphi,fit2->GetParameter(4)) ;
365 br2->SetBinError(i,iphi,fit2->GetParError(4)) ;
366 cr2->SetBinContent(i,iphi,fit2->GetParameter(5)) ;
367 cr2->SetBinError(i,iphi,fit2->GetParError(5)) ;
369 hpcopy->GetXaxis()->SetRangeUser(0.05,0.3) ;
379 // if(getchar()=='q')return ;
382 fbg1->SetParameters(fit1->GetParameter(3),fit1->GetParameter(4),fit1->GetParameter(5));
383 fbg2->SetParameters(fit2->GetParameter(3),fit2->GetParameter(4),fit2->GetParameter(5),fit2->GetParameter(6));
385 hpm->Multiply(fbg1) ;
386 hpm2->Multiply(fbg2) ;
392 fgs->SetParameters(hp->Integral(13,15)/3.,fit1->GetParameter(1),fit1->GetParameter(2),0.) ;
393 fgs->SetParLimits(1,0.05,0.15);
394 fgs->FixParameter(1,mr1->GetBinContent(i)) ;
395 fgs->FixParameter(2,sr1->GetBinContent(i)) ;
397 hp->Fit(fgs,"Q","",rangeMin,rangeMax) ;
398 hp->SetMaximum(hp2->GetMaximum()*1.1) ;
399 hp->SetMinimum(hp2->GetMinimum()*1.1) ;
400 hp->SetMarkerStyle(20) ;
401 hp->SetMarkerSize(0.7) ;
403 //0.96 - due to 2 sigma
404 nr1->SetBinContent(i,iphi,fgs->GetParameter(0)) ;
405 nr1->SetBinError(i,iphi,fgs->GetParError(0)) ;
407 Int_t intBinMin=hp->GetXaxis()->FindBin(fgs->GetParameter(1)-4.*TMath::Abs(fgs->GetParameter(2))) ;
408 Int_t intBinMax=hp->GetXaxis()->FindBin(fgs->GetParameter(1)+4.*TMath::Abs(fgs->GetParameter(2))) ;
409 Double_t errStat=hpm->Integral(intBinMin,intBinMax);
410 Double_t npiInt=hp->Integral(intBinMin,intBinMax) ;
411 Double_t norm=fbg1->GetParameter(0) ;
412 Double_t normErr=fbg1->GetParError(0) ;
414 nr1int->SetBinContent(i,iphi,npiInt) ;
415 nr1int->SetBinError(i,iphi,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
417 hp2->GetXaxis()->SetRangeUser(0.05,0.3) ;
418 hp2->SetMaximum(hp2->GetMaximum()*1.1) ;
419 hp2->SetMinimum(hp2->GetMinimum()*1.1) ;
420 hp2->SetMarkerStyle(20) ;
421 hp2->SetMarkerSize(0.7) ;
423 fgs->FixParameter(1,mr2->GetBinContent(i)) ;
424 fgs->FixParameter(2,sr2->GetBinContent(i)) ;
425 hp2->Fit(fgs,"Q","",rangeMin,rangeMax) ;
426 nr2->SetBinContent(i,iphi,fgs->GetParameter(0)) ;
427 nr2->SetBinError(i,iphi,fgs->GetParError(0)) ;
428 npiInt=hp2->Integral(intBinMin,intBinMax) ;
429 norm=fbg2->GetParameter(0) ;
430 normErr=fbg2->GetParError(0) ;
432 nr2int->SetBinContent(i,iphi,npiInt) ;
433 nr2int->SetBinError(i,iphi,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
452 cout<<"Fit histos init... ";
453 TH1D * v2A1 = new TH1D("v2all1","v_{2} two harmonics",nbin,xa) ;
454 TH1D * v2A2 = new TH1D("v2all2","v_{2} two harmonics",nbin,xa) ;
455 TH1D * v2A1int = new TH1D("v2all1int","v_{2} two harmonics",nbin,xa) ;
456 TH1D * v2A2int = new TH1D("v2all2int","v_{2} two harmonics",nbin,xa) ;
458 TH1D * v2F1 = new TH1D("v2gap1","v_{2} two harmonics",nbin,xa) ;
459 TH1D * v2F2 = new TH1D("v2gap2","v_{2} two harmonics",nbin,xa) ;
460 TH1D * v2F1int = new TH1D("v2gap1int","v_{2} two harmonics",nbin,xa) ;
461 TH1D * v2F2int = new TH1D("v2gap2int","v_{2} two harmonics",nbin,xa) ;
463 TH1D * v2A1s= new TH1D("v2all1S","v_{2} single harmonic",nbin,xa) ;
464 TH1D * v2A2s= new TH1D("v2all2S","v_{2} single harmonic",nbin,xa) ;
465 TH1D * v2A1ints= new TH1D("v2all1intS","v_{2} single harmonic",nbin,xa) ;
466 TH1D * v2A2ints= new TH1D("v2all2intS","v_{2} single harmonic",nbin,xa) ;
468 TH1D * v2F1s= new TH1D("v2gap1S","v_{2} single harmonic",nbin,xa) ;
469 TH1D * v2F2s= new TH1D("v2gap2S","v_{2} single harmonic",nbin,xa) ;
470 TH1D * v2F1ints= new TH1D("v2gap1intS","v_{2} single harmonic",nbin,xa) ;
471 TH1D * v2F2ints= new TH1D("v2gap2intS","v_{2} single harmonic",nbin,xa) ;
474 TH2D * hA1 = nr1->Clone("hA1");
475 TH2D * hA2 = nr2->Clone("hA2");
476 TH2D * hA1int = nr1int->Clone("hA1int");
477 TH2D * hA2int = nr2int->Clone("hA2int");
479 cout<<"Finished"<<endl;
481 Int_t col[4]={kRed,kBlue,kGreen+3,kMagenta} ;
482 Int_t sym[4]={20,21,24,25} ;
484 // TF1 * fit = new TF1("fit","[0]*(1+2.*[1]*cos(2*(x-[2])))",0.,5.) ;
485 TF1 * v2fit = new TF1("v2fit","[0]*(1+2.*[1]*cos(2.*x)+2.*[2]*cos(4.*x) )",0.,5.) ;
486 v2fit->SetLineColor(2) ;
487 TF1 * v2fit2 = new TF1("v2fit","[0]*(1+2.*[1]*cos(2.*x) )",0.,5.) ;
488 TCanvas * c = new TCanvas("fit","fit") ;
491 cout<<"Fit dNdPhi... ";
492 for(Int_t i=1; i<=v2A1->GetNbinsX() ;i++){ //over pt
495 sprintf(hname,"A1%d",i);
496 TH1D * tmp0 = hA1->ProjectionY(hname,i,i) ;
497 tmp0->SetTitle(Form("%3.1f<p_{t}<%3.1f GeV",xa[i-1],xa[i])) ;
498 tmp0->SetXTitle("#Delta#phi (rad)") ;
499 tmp0->SetYTitle("dN/d#phi (a.u.)") ;
500 tmp0->SetMarkerStyle(sym[0]) ;
501 tmp0->SetMarkerColor(col[0]) ;
503 v2fit2->SetLineStyle(1) ;
504 v2fit2->SetLineColor(col[0]) ;
505 v2fit2->SetParameters(tmp0->GetBinContent(5),0.1) ;
506 tmp0->Fit(v2fit2,"Qi") ;
507 v2A1s->SetBinContent(i,v2fit2->GetParameter(1)) ;
508 v2A1s->SetBinError(i,v2fit2->GetParError(1)) ;
511 sprintf(hname,"A2%d",i);
512 //tmp2d = (TH2D*)hA2->Clone(hname);
513 TH1D * tmp1 = hA2->ProjectionY(hname,i,i) ;
514 tmp1->SetMarkerStyle(sym[1]) ;
515 tmp1->SetMarkerColor(col[1]) ;
516 v2fit2->SetLineColor(col[1]) ;
517 v2fit2->SetParameters(tmp1->GetBinContent(5),0.1) ;
518 v2fit2->SetLineColor(col[1]) ;
519 tmp1->Fit(v2fit2,"Qi") ;
520 v2A2s->SetBinContent(i,v2fit2->GetParameter(1)) ;
521 v2A2s->SetBinError(i,v2fit2->GetParError(1)) ;
523 sprintf(hname,"A1int%d",i);
524 TH1D * tmp2 = hA1int->ProjectionY(hname,i,i) ;
525 tmp2->SetMarkerStyle(sym[2]) ;
526 tmp2->SetMarkerColor(col[2]) ;
527 v2fit2->SetParameters(tmp2->GetBinContent(5),0.1) ;
528 v2fit2->SetLineStyle(5) ;
529 v2fit2->SetLineColor(col[2]) ;
530 tmp2->Fit(v2fit2,"i") ;
531 v2A1ints->SetBinContent(i,v2fit2->GetParameter(1)) ;
532 v2A1ints->SetBinError(i,v2fit2->GetParError(1)) ;
534 sprintf(hname,"A2int%d",i);
535 //tmp2d = (TH2D*)hA2int->Clone(hname);
536 TH1D * tmp3 = hA2int->ProjectionY(hname,i,i) ;
537 tmp3->SetMarkerStyle(sym[3]) ;
538 tmp3->SetMarkerColor(col[3]) ;
539 v2fit2->SetParameters(tmp3->GetBinContent(5),0.1) ;
540 v2fit2->SetLineColor(col[3]) ;
541 tmp3->Fit(v2fit2,"Qi") ;
542 v2A2ints->SetBinContent(i,v2fit2->GetParameter(1)) ;
543 v2A2ints->SetBinError(i,v2fit2->GetParError(1)) ;
547 // tmp2->Scale(6.28) ;
548 // tmp3->Scale(6.28) ;
552 cout<<"Finished"<<endl;
554 TH1D *effcorrect = v2A1s->Clone("effcorrect");
555 TH1D *effcorrect2 = v2A1s->Clone("effcorrect2");
557 Double_t eff=0, err=0;
558 if(strcmp(PID,"Both2core")==0){
559 if(cen==0) {eff=.03; err=.04;}
560 if(cen==1) {eff=.08; err=.04;}
561 if(cen==2) {eff=.03; err=.03;}
562 if(cen==3) {eff=.02; err=.02;}
563 if(cen==4) {eff=.02; err=.02;}
564 if(cen==5) {eff=.01; err=.02;}
565 if(cen==10) {eff=.05; err=.02;}
566 if(cen==11) {eff=.02; err=.02;}
567 if(cen==20) {eff=.02; err=.02;}
568 if(cen==21) {eff=.05; err=.02;}
570 eff=eff*2.*TMath::Pi()/5./sin(2.*TMath::Pi()/5.)/4.;
571 err=err*2.*TMath::Pi()/5./sin(2.*TMath::Pi()/5.)/4.;
573 cout<<"res: T="<<resT<<", resolution="<<res<<endl;
577 cout<<"Pi0 efficiency correction before resolution correction: "<<eff<<"+-"<<err<<endl;
580 Double_t x[100],y[100],ex[100],ey[100] ;
582 for(Int_t i=0;i<effcorrect->GetNbinsX();i++){
583 effcorrect->SetBinContent(i,eff);
584 effcorrect->SetBinError(i,err);
585 effcorrect2->SetBinContent(i,eff);
586 effcorrect2->SetBinError(i,0);
591 v2A1s->Add(effcorrect);
592 v2A2s->Add(effcorrect);
593 v2A1ints->Add(effcorrect);
594 v2A2ints->Add(effcorrect);
599 v2A1ints->Scale(res) ;
600 v2A2ints->Scale(res) ;
605 v2A1s->SetMarkerStyle(sym[0]) ;
606 v2A1s->SetMarkerColor(col[0]) ;
607 v2A2s->SetMarkerStyle(sym[1]) ;
608 v2A2s->SetMarkerColor(col[1]) ;
609 v2A1ints->SetMarkerStyle(sym[2]) ;
610 v2A1ints->SetMarkerColor(col[2]) ;
611 v2A2ints->SetMarkerStyle(sym[3]) ;
612 v2A2ints->SetMarkerColor(col[3]) ;
614 TCanvas * c2 = new TCanvas("v2","v2") ;
615 TH1D * box = new TH1D("box","V_{2}{RP} with different fit methods",150,0.,20.) ;
616 box->SetMinimum(0.) ;
617 box->SetMaximum(0.5) ;
618 box->SetXTitle("p_{t}^{#pi} (GeV/c)") ;
619 box->SetYTitle("v{2}") ;
621 v2A1s->Draw("same") ;
622 v2A1ints->Draw("same") ;
623 v2A2s->Draw("same") ;
624 v2A2ints->Draw("same") ;
626 TLegend * l = new TLegend(0.6,0.7,0.9,0.9) ;
627 l->AddEntry(v2A1s,"Fit, pol1","p") ;
628 l->AddEntry(v2A2s,"Fit, pol2","p") ;
629 l->AddEntry(v2A1ints,"Int, pol1","p") ;
630 l->AddEntry(v2A2ints,"Int, pol2","p") ;
636 TCanvas * cmass = new TCanvas("Mass","Pi0 mass vs pT") ;
638 mr1->SetMarkerStyle(20);
639 mr2->SetMarkerStyle(21);
641 mr2->SetLineColor(2);
642 mr2->SetMarkerColor(2);
645 TLegend * l = new TLegend(0.6,0.7,0.9,0.9) ;
646 l->AddEntry(mr1,"pol1","p") ;
647 l->AddEntry(mr2,"pol2","p") ;
651 TCanvas * cwidth = new TCanvas("Width","Pi0 width vs pT") ;
653 sr1->SetMarkerStyle(20);
654 sr2->SetMarkerStyle(21);
656 sr2->SetLineColor(2);
657 sr2->SetMarkerColor(2);
660 TLegend * l = new TLegend(0.6,0.7,0.9,0.9) ;
661 l->AddEntry(mr1,"pol1","p") ;
662 l->AddEntry(mr2,"pol2","p") ;
666 //!!!!!!!!!!!!!!!!!!!!!!final v2 plot
667 TCanvas * c5 = new TCanvas("v2_res","v2 res") ;
670 TH1D * v2Astat = new TH1D("v2allStat","v_{2}",nbin,xa) ;
671 TH1D * v2Asys = new TH1D("v2allSys","v_{2}",nbin,xa) ;
673 for(Int_t i=1; i<=v2A1->GetNbinsX() ;i++){
675 v2A1s->GetBinContent(i)/v2A1s->GetBinError(i)/v2A1s->GetBinError(i)
676 +v2A2s->GetBinContent(i)/v2A2s->GetBinError(i)/v2A2s->GetBinError(i)
677 +v2A1ints->GetBinContent(i)/v2A1ints->GetBinError(i)/v2A1ints->GetBinError(i)
678 +v2A2ints->GetBinContent(i)/v2A2ints->GetBinError(i)/v2A2ints->GetBinError(i) ;
680 1./v2A1s->GetBinError(i)/v2A1s->GetBinError(i)
681 +1./v2A2s->GetBinError(i)/v2A2s->GetBinError(i)
682 +1./v2A1ints->GetBinError(i)/v2A1ints->GetBinError(i)
683 +1./v2A2ints->GetBinError(i)/v2A2ints->GetBinError(i) ;
687 Double_t rms = (v2A1s->GetBinContent(i) - mean)*(v2A1s->GetBinContent(i) - mean)/v2A1s->GetBinError(i)/v2A1s->GetBinError(i) + (v2A2s->GetBinContent(i) - mean)*(v2A2s->GetBinContent(i) - mean)/v2A2s->GetBinError(i)/v2A2s->GetBinError(i) + (v2A1ints->GetBinContent(i) - mean)*(v2A1ints->GetBinContent(i) - mean)/v2A1ints->GetBinError(i)/v2A1ints->GetBinError(i) + (v2A2ints->GetBinContent(i) - mean)*(v2A2ints->GetBinContent(i) - mean)/v2A2ints->GetBinError(i)/v2A2ints->GetBinError(i);
688 Double_t statErr = TMath::Min(v2A1s->GetBinError(i),v2A2s->GetBinError(i)) ;
689 // TMath::Min(v2A1ints->GetBinError(i),v2A2ints->GetBinError(i))) ;
694 Double_t ptbin = v2A1->GetBinCenter(i);
695 Double_t sys=TMath::Sqrt(rms);
697 cout<<"Sys error for pT="<<ptbin<<" is "<<sys<<" in percent: "<<sys/mean<<endl;
698 cout<<"Sys error due to eff corr: "<<eff/res<<" in percent: "<<eff/res/mean<<endl;
700 Double_t mean_corr = mean ;// - v2ch * Ntrack_cen7080 / Ntrack_cen;
702 v2Astat->SetBinContent(i,mean_corr) ;
703 v2Astat->SetBinError(i,statErr) ;
704 v2Asys->SetBinContent(i,mean_corr) ;
705 v2Asys->SetBinError(i,sys) ;
710 v2Astat->Add(effcorrect2);
711 v2Asys->Add(effcorrect);
714 v2Asys->Scale(1/res) ;
715 v2Astat->Scale(1/res) ;
717 v2Asys->SetFillStyle(1001) ;
718 v2Asys->SetFillColor(kYellow) ;
720 v2Astat->SetMarkerColor(kOrange+7) ;
721 v2Astat->SetMarkerStyle(20) ;
722 v2Astat->SetLineColor(kOrange+7) ;
725 v2Fsys->SetFillStyle(0) ;
726 v2Fsys->SetFillColor(0) ;
727 v2Fsys->SetLineColor(kBlue+3) ;
730 v2Astat->SetMarkerColor(kOrange+7) ;
731 v2Astat->SetMarkerStyle(20) ;
732 v2Astat->SetLineColor(kOrange+7) ;
735 v2Fstat->SetMarkerColor(kBlue+3) ;
736 v2Fstat->SetMarkerStyle(21) ;
737 v2Fstat->SetLineColor(kBlue+3) ;
740 v2Asys->SetXTitle("p_{t} (GeV/c)") ;
741 v2Asys->SetYTitle("v_{2}") ;
744 v2Asys->Draw("sameE1") ;
745 // v2Fsys->Draw("E2same") ;
746 v2Astat->Draw("same") ;
747 // v2Fstat->Draw("same") ;
750 TFile fout("v2_method1_QA.root","update");
751 sprintf(nname,"v2sys_m1_%s_%s_%d_%d",w,PID,cen,side) ;
752 v2Asys->SetName(nname) ;
753 sprintf(nname,"v2stat_m1_%s_%s_%d_%d",w,PID,cen,side) ;
754 v2Astat->SetName(nname) ;
756 v2Asys->Write(0,TObject::kOverwrite) ;
757 v2Astat->Write(0,TObject::kOverwrite) ;
762 //-----------------------------------------------------------------------------
766 //////////////////////////////////////////////////////////////////////
768 // ROOT style macro for the TRD TDR
770 //////////////////////////////////////////////////////////////////////
772 gStyle->SetPalette(1);
773 gStyle->SetCanvasBorderMode(-1);
774 gStyle->SetCanvasBorderSize(1);
775 gStyle->SetCanvasColor(10);
777 gStyle->SetFrameFillColor(10);
778 gStyle->SetFrameBorderSize(1);
779 gStyle->SetFrameBorderMode(-1);
780 gStyle->SetFrameLineWidth(1.2);
781 gStyle->SetFrameLineColor(1);
783 gStyle->SetHistFillColor(0);
784 gStyle->SetHistLineWidth(1);
785 gStyle->SetHistLineColor(1);
787 gStyle->SetPadColor(10);
788 gStyle->SetPadBorderSize(1);
789 gStyle->SetPadBorderMode(-1);
791 gStyle->SetStatColor(10);
792 gStyle->SetTitleColor(kBlack,"X");
793 gStyle->SetTitleColor(kBlack,"Y");
795 gStyle->SetLabelSize(0.04,"X");
796 gStyle->SetLabelSize(0.04,"Y");
797 gStyle->SetLabelSize(0.04,"Z");
798 gStyle->SetTitleSize(0.04,"X");
799 gStyle->SetTitleSize(0.04,"Y");
800 gStyle->SetTitleSize(0.04,"Z");
801 gStyle->SetTitleFont(42,"X");
802 gStyle->SetTitleFont(42,"Y");
803 gStyle->SetTitleFont(42,"X");
804 gStyle->SetLabelFont(42,"X");
805 gStyle->SetLabelFont(42,"Y");
806 gStyle->SetLabelFont(42,"Z");
807 gStyle->SetStatFont(42);
809 gStyle->SetTitleOffset(1.0,"X");
810 gStyle->SetTitleOffset(1.4,"Y");
812 gStyle->SetFillColor(kWhite);
813 gStyle->SetTitleFillColor(kWhite);
815 gStyle->SetOptDate(0);
816 gStyle->SetOptTitle(1);
817 gStyle->SetOptStat(0);
818 gStyle->SetOptFit(0);