PWGGA/PHOSTasks/PHOS_PbPb/macros/Flow/Method1.C

   1 #include "Methods.h"
   2
   3 const Double_t kMean=0.135 ; //Approximate peak position to facilitate error estimate
   4 Int_t effcor=1; //correct on efficiency
   5
   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 ;
   9 }
  10 Double_t PeakWidth(Double_t pt){
  11   //fit to the measured peak width
  12   Double_t a=0.0068 ;
  13   Double_t b=0.0025 ;
  14   Double_t c=0.000319 ;
  15   return TMath::Sqrt(a*a+b*b/pt/pt+c*c*pt*pt) ;
  16 }
  17
  18 Double_t CB(Double_t * x, Double_t * par){
  19   //Parameterization of Real/Mixed ratio
  20   Double_t m=par[1] ;
  21   Double_t s=par[2] ;
  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) */;
  24 }
  25 Double_t CB2(Double_t * x, Double_t * par){
  26   //Another parameterization of Real/Mixed ratio
  27   Double_t m=par[1] ;
  28   Double_t s=par[2] ;
  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)*/;
  31 }
  32 Double_t CBs(Double_t * x, Double_t * par){
  33   //Parameterizatin of signal
  34   Double_t m=par[1] ;
  35   Double_t s=par[2] ;
  36   Double_t dx=(x[0]-m)/s ;
  37   return par[0]*exp(-dx*dx/2.)/TMath::Sqrt(TMath::TwoPi())/s+par[3] ;
  38 }
  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)*/;
  42 }
  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)*/;
  46 }
  47
  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
  51
  52   //Fit Real/Mixed ratio, normalize Mixed and subtract it from Real
  53
  54 //  TH1F * hCen = GetCen() ;
  55 //  TH2F * hCenPHOS = GetCenPHOS() ;
  56
  57
  58   char kind[15], w[25];
  59   sprintf(w,""); //Weight: no weight - "NW", with weight - ""
  60
  61   char PID[25] ;
  62   sprintf(PID,"Both2core") ;
  63   char key[55];
  64   char kind[15],kind2[15] ; //Reaction plane
  65   char d[15]; //detector
  66   sprintf(kind,"") ; //"" for real
  67   sprintf(kind2,"") ; //"" for mixed
  68
  69   if(side==1)sprintf(d,"V0C");
  70   else if(side==0)sprintf(d,"V0A");
  71   else sprintf(d,"TPC");
  72
  73   TH3D *h3DR, *h3DM;
  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);
  77
  78   TH1F * hCen = GetCen() ;
  79   TH2F * hCenPHOS = GetCenPHOS() ;
  80 cout<<"Finished"<<endl;
  81
  82   //  sprintf(kind,"") ; //"" for all
  83   const Double_t nphi=5;
  84   Double_t xphi[6];
  85   for(Int_t i=0;i<=nphi;i++){
  86         xphi[i]=TMath::PiOver2()*i/nphi ;
  87 //      cout<<xphi[i]<<" ";
  88   }
  89 //cout<<endl;
  90
  91   if(cen==21)
  92   sprintf(inname,"../flow11h_Apr14.root");
  93   //Read resolution
  94   TFile * f = new TFile(inname) ;
  95
  96   TH2F * hcos2AC = (TH2F*)f->Get("cos2V0AC");
  97   TH2F * hcos2FAC= (TH2F*)f->Get("cos2FAC");
  98   Double_t resAC = GetCos(hcos2AC,cen);
  99   resAC=resAC ;
 100
 101   hcos2AC = (TH2F*)f->Get("cos2V0ATPC");
 102   Double_t resAT = GetCos(hcos2AC,cen);
 103   resAT=resAT ;
 104
 105   hcos2AC = (TH2F*)f->Get("cos2V0CTPC");
 106   Double_t resCT = GetCos(hcos2AC,cen);
 107   resCT=resCT ;
 108
 109   hcos2AC = (TH2F*)f->Get("cos2AC");
 110   Double_t resT = GetRes(hcos2AC,cen);
 111
 112   //side = 1 - VOC, side = 0 - V0A
 113   //side = 2 - TPC 3 sub method, 3 - TPC 2 sub method
 114
 115   Double_t res=1;
 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);
 119   else res = resT;
 120
 121 //  res=1./res ;
 122
 123   cout<<"RP resolution = "<< res <<endl;
 124
 125   PPRstyle();
 126 cout<<"Init... ";
 127   //Fit real only
 128   //Linear Bg
 129   char key2[155];
 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) ;
 143
 144   //Quadratic Bg
 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) ;
 158
 159   TH2D * nr2int = new TH2D(key2,"Raw yield, integrated",nbin,xa,nphi,xphi) ;
 160
 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) ;
 172
 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}") ;
 180
 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) ;
 186
 187   //calculate average mass and width for all pt bins
 188   TCanvas * cav = new TCanvas("Average","Average",10,10,800,750) ;
 189   cav->Divide(3,4) ;
 190
 191 cout<<"Finished"<<endl;
 192
 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") ;
 200     hpav->Sumw2() ;
 201     hpmav->Sumw2() ;
 202     if(i<7){
 203       hpav->Rebin(2) ;
 204       hpmav->Rebin(2) ;
 205     }
 206     else{
 207       hpav->Rebin(2) ;
 208       hpmav->Rebin(2) ;
 209     }
 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") ;
 217
 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) ;
 223
 224     fit1->SetParameters(TMath::Min(0.3,0.001*i),0.14,0.01,0.008,-0.0002,0) ;
 225
 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
 229
 230     cav->cd(i) ;
 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) ;
 236
 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)) ;
 241
 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) ;
 247
 248     cav->Update() ;
 249
 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)) ;
 254
 255
 256     delete hpav ;
 257 //    delete hpcopyav ;
 258     delete hp2av ;
 259     delete hpmav;
 260     delete hpm2av;
 261   }
 262 cout<<"Finished"<<endl;
 263
 264   h3DR->GetYaxis()->SetRangeUser(0.,19.) ;
 265   h3DM->GetYaxis()->SetRangeUser(0.,19.) ;
 266
 267 cout<<"Do dNdPhi bins... ";
 268 //!!!!!!!!!!!!!!!!!DNDPHI bins:
 269
 270   for(Int_t iphi=1;iphi<=nphi;iphi++){
 271
 272     h3DR->GetZaxis()->SetRange(iphi,iphi) ;
 273     h3DM->GetZaxis()->SetRange(iphi,iphi) ;
 274
 275     TH2D * h = (TH2D*)h3DR->Project3D("yx");
 276     TH2D * hm= (TH2D*)h3DM->Project3D("yx") ;
 277     h->SetName("sliceRe") ;
 278     hm->SetName("sliceMi") ;
 279
 280     h3DR->GetZaxis()->SetRange(11-iphi,11-iphi) ;
 281     h3DM->GetZaxis()->SetRange(11-iphi,11-iphi) ;
 282
 283     TH2D * ha = (TH2D*)h3DR->Project3D("yx");
 284     TH2D * hma= (TH2D*)h3DM->Project3D("yx") ;
 285
 286     h->Add(ha) ; delete ha ;
 287     hm->Add(hma); delete hma ;
 288
 289     sprintf(key2,"mggFit%d_Signal",iphi) ;
 290     TCanvas * c3 = new TCanvas(key2,key2,10,10,800,750) ;
 291     c3->Divide(3,4) ;
 292
 293     sprintf(key2,"mggFit%d",iphi) ;
 294     TCanvas * c1 = new TCanvas(key2,key2,10,10,800,750) ;
 295     c1->Divide(3,4) ;
 296     c1->cd(0) ;
 297
 298 //    TCanvas * c2=0,*c4=0,*c5=0,*c6=0 ;
 299     TAxis * pta=h->GetYaxis() ;
 300     TAxis * ma=h->GetXaxis() ;
 301
 302 //loop over pT bins
 303     for(Int_t i=1;i<=nbin;i++){
 304       c1->cd(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) ;
 309       hp->Sumw2() ;
 310       TH1D * hpm= hm->ProjectionX("mi",imin,imax) ;
 311       hpm->Sumw2() ;
 312
 313       if(i<7){
 314         hp->Rebin(2) ;
 315         hpm->Rebin(2) ;
 316       }
 317       else{
 318         hp->Rebin(2) ;
 319         hpm->Rebin(2) ;
 320       }
 321
 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) ;
 332
 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)) ;
 336
 337 //      fit1->SetParLimits(2,0.005,0.015) ;
 338 //      fit1->SetParLimits(0,0.,1.) ;
 339
 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 ;
 343
 344       hpcopy->Fit("fit1","NQ","",rangeMin,rangeMax) ;
 345       hpcopy->Fit("fit1","MQ","",rangeMin,rangeMax) ;
 346
 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)) ;
 351
 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) ;
 358
 359       hpcopy->Fit("fit2","+NQ","",rangeMin,rangeMax) ;
 360       hpcopy->Fit("fit2","+MQ","",rangeMin,rangeMax) ;
 361
 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)) ;
 368
 369       hpcopy->GetXaxis()->SetRangeUser(0.05,0.3) ;
 370 //      hpcopy->Draw() ;
 371 /*      if(c5)
 372         c5->Update() ;
 373       else
 374         if(c2)
 375           c2->Update() ;
 376         else
 377 */
 378           c1->Update() ;
 379       //    if(getchar()=='q')return ;
 380
 381
 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));
 384
 385       hpm->Multiply(fbg1) ;
 386       hpm2->Multiply(fbg2) ;
 387       hp->Add(hpm,-1.) ;
 388       hp2->Add(hpm2,-1.) ;
 389
 390       c3->cd(i) ;
 391
 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)) ;
 396
 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) ;
 402
 403       //0.96 - due to 2 sigma
 404       nr1->SetBinContent(i,iphi,fgs->GetParameter(0)) ;
 405       nr1->SetBinError(i,iphi,fgs->GetParError(0)) ;
 406
 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) ;
 413       if(npiInt>0.){
 414         nr1int->SetBinContent(i,iphi,npiInt) ;
 415         nr1int->SetBinError(i,iphi,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
 416       }
 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) ;
 422
 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) ;
 431       if(npiInt>0.){
 432         nr2int->SetBinContent(i,iphi,npiInt) ;
 433         nr2int->SetBinError(i,iphi,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
 434       }
 435       hp2->SetTitle(key) ;
 436       hp2->Draw() ;
 437
 438       c3->Update() ;
 439
 440       delete hp ;
 441 //      delete hp2 ;
 442 //      delete hpcopy ;
 443       delete hpm ;
 444       delete hpm2 ;
 445     }
 446
 447     delete h ;
 448     delete hm ;
 449
 450   }
 451
 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) ;
 457
 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) ;
 462
 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) ;
 467
 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) ;
 472
 473
 474   TH2D * hA1 = nr1->Clone("hA1");
 475   TH2D * hA2 = nr2->Clone("hA2");
 476   TH2D * hA1int = nr1int->Clone("hA1int");
 477   TH2D * hA2int = nr2int->Clone("hA2int");
 478
 479 cout<<"Finished"<<endl;
 480
 481   Int_t col[4]={kRed,kBlue,kGreen+3,kMagenta} ;
 482   Int_t sym[4]={20,21,24,25} ;
 483
 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") ;
 489   c->Divide(3,4);
 490
 491 cout<<"Fit dNdPhi... ";
 492   for(Int_t i=1; i<=v2A1->GetNbinsX() ;i++){ //over pt
 493     c->cd(i);
 494     char hname[100];
 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]) ;
 502
 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)) ;
 509
 510
 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)) ;
 522
 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)) ;
 533
 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)) ;
 544
 545     tmp0->Draw() ;
 546     tmp1->Draw("same") ;
 547 //    tmp2->Scale(6.28) ;
 548 //    tmp3->Scale(6.28) ;
 549     tmp2->Draw("same") ;
 550     tmp3->Draw("same") ;
 551   }
 552 cout<<"Finished"<<endl;
 553
 554   TH1D *effcorrect = v2A1s->Clone("effcorrect");
 555   TH1D *effcorrect2 = v2A1s->Clone("effcorrect2");
 556
 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;}
 569
 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.;
 572
 573 cout<<"res: T="<<resT<<", resolution="<<res<<endl;
 574
 575         eff=eff*res/resT;
 576         err=err*res/resT;
 577         cout<<"Pi0 efficiency correction before resolution correction: "<<eff<<"+-"<<err<<endl;
 578   }
 579
 580   Double_t x[100],y[100],ex[100],ey[100] ;
 581
 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);
 587   }
 588
 589 /*
 590   if(effcor){
 591         v2A1s->Add(effcorrect);
 592         v2A2s->Add(effcorrect);
 593         v2A1ints->Add(effcorrect);
 594         v2A2ints->Add(effcorrect);
 595   }
 596
 597   v2A1s->Scale(res) ;
 598   v2A2s->Scale(res) ;
 599   v2A1ints->Scale(res) ;
 600   v2A2ints->Scale(res) ;
 601 */
 602
 603 //return 1;
 604
 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]) ;
 613
 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}") ;
 620   box->Draw() ;
 621   v2A1s->Draw("same") ;
 622   v2A1ints->Draw("same") ;
 623   v2A2s->Draw("same") ;
 624   v2A2ints->Draw("same") ;
 625
 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") ;
 631
 632   l->Draw() ;
 633
 634 //  return ;
 635
 636   TCanvas * cmass = new TCanvas("Mass","Pi0 mass vs pT") ;
 637
 638   mr1->SetMarkerStyle(20);
 639   mr2->SetMarkerStyle(21);
 640   mr1->Draw();
 641   mr2->SetLineColor(2);
 642   mr2->SetMarkerColor(2);
 643   mr2->Draw("same");
 644
 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") ;
 648
 649   l->Draw() ;
 650
 651   TCanvas * cwidth = new TCanvas("Width","Pi0 width vs pT") ;
 652
 653   sr1->SetMarkerStyle(20);
 654   sr2->SetMarkerStyle(21);
 655   sr1->Draw();
 656   sr2->SetLineColor(2);
 657   sr2->SetMarkerColor(2);
 658   sr2->Draw("same");
 659
 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") ;
 663
 664   l->Draw() ;
 665
 666 //!!!!!!!!!!!!!!!!!!!!!!final v2 plot
 667   TCanvas * c5 = new TCanvas("v2_res","v2 res") ;
 668
 669
 670   TH1D * v2Astat = new TH1D("v2allStat","v_{2}",nbin,xa) ;
 671   TH1D * v2Asys = new TH1D("v2allSys","v_{2}",nbin,xa) ;
 672
 673   for(Int_t i=1; i<=v2A1->GetNbinsX() ;i++){
 674     Double_t mean=
 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) ;
 679     Double_t weight=
 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) ;
 684
 685     mean/=weight ;
 686
 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))) ;
 690 cout<<rms<<endl;
 691
 692     rms/=weight;
 693
 694     Double_t ptbin = v2A1->GetBinCenter(i);
 695     Double_t sys=TMath::Sqrt(rms);
 696
 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;
 699
 700     Double_t mean_corr = mean ;// - v2ch * Ntrack_cen7080 / Ntrack_cen;
 701
 702     v2Astat->SetBinContent(i,mean_corr) ;
 703     v2Astat->SetBinError(i,statErr) ;
 704     v2Asys->SetBinContent(i,mean_corr) ;
 705     v2Asys->SetBinError(i,sys) ;
 706
 707   }
 708
 709   if(effcor){
 710         v2Astat->Add(effcorrect2);
 711         v2Asys->Add(effcorrect);
 712   }
 713
 714   v2Asys->Scale(1/res) ;
 715   v2Astat->Scale(1/res) ;
 716
 717   v2Asys->SetFillStyle(1001) ;
 718   v2Asys->SetFillColor(kYellow) ;
 719
 720   v2Astat->SetMarkerColor(kOrange+7) ;
 721   v2Astat->SetMarkerStyle(20) ;
 722   v2Astat->SetLineColor(kOrange+7) ;
 723
 724   /*
 725   v2Fsys->SetFillStyle(0) ;
 726   v2Fsys->SetFillColor(0) ;
 727   v2Fsys->SetLineColor(kBlue+3) ;
 728   */
 729
 730   v2Astat->SetMarkerColor(kOrange+7) ;
 731   v2Astat->SetMarkerStyle(20) ;
 732   v2Astat->SetLineColor(kOrange+7) ;
 733
 734   /*
 735   v2Fstat->SetMarkerColor(kBlue+3) ;
 736   v2Fstat->SetMarkerStyle(21) ;
 737   v2Fstat->SetLineColor(kBlue+3) ;
 738   */
 739
 740   v2Asys->SetXTitle("p_{t} (GeV/c)") ;
 741   v2Asys->SetYTitle("v_{2}") ;
 742
 743   box->Draw() ;
 744   v2Asys->Draw("sameE1") ;
 745   //  v2Fsys->Draw("E2same") ;
 746   v2Astat->Draw("same") ;
 747   //  v2Fstat->Draw("same") ;
 748
 749   char nname[255];
 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) ;
 755
 756   v2Asys->Write(0,TObject::kOverwrite) ;
 757   v2Astat->Write(0,TObject::kOverwrite) ;
 758   fout.Close() ;
 759 }
 760
 761
 762 //-----------------------------------------------------------------------------
 763 PPRstyle()
 764 {
 765
 766   //////////////////////////////////////////////////////////////////////
 767   //
 768   // ROOT style macro for the TRD TDR
 769   //
 770   //////////////////////////////////////////////////////////////////////
 771
 772   gStyle->SetPalette(1);
 773   gStyle->SetCanvasBorderMode(-1);
 774   gStyle->SetCanvasBorderSize(1);
 775   gStyle->SetCanvasColor(10);
 776
 777   gStyle->SetFrameFillColor(10);
 778   gStyle->SetFrameBorderSize(1);
 779   gStyle->SetFrameBorderMode(-1);
 780   gStyle->SetFrameLineWidth(1.2);
 781   gStyle->SetFrameLineColor(1);
 782
 783   gStyle->SetHistFillColor(0);
 784   gStyle->SetHistLineWidth(1);
 785   gStyle->SetHistLineColor(1);
 786
 787   gStyle->SetPadColor(10);
 788   gStyle->SetPadBorderSize(1);
 789   gStyle->SetPadBorderMode(-1);
 790
 791   gStyle->SetStatColor(10);
 792   gStyle->SetTitleColor(kBlack,"X");
 793   gStyle->SetTitleColor(kBlack,"Y");
 794
 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);
 808
 809   gStyle->SetTitleOffset(1.0,"X");
 810   gStyle->SetTitleOffset(1.4,"Y");
 811
 812   gStyle->SetFillColor(kWhite);
 813   gStyle->SetTitleFillColor(kWhite);
 814
 815   gStyle->SetOptDate(0);
 816   gStyle->SetOptTitle(1);
 817   gStyle->SetOptStat(0);
 818   gStyle->SetOptFit(0);
 819
 820 }
 821