PWGGA/PHOSTasks/PHOS_PbPb/macros/MakeMmixPi0.C

   1 /* $Id$ */
   2
   3 const Double_t kMean=0.136 ; //Approximate peak position to facilitate error estimate
   4
   5 //-----------------------------------------------------------------------------
   6 void MakeMmixPi0(const TString filename = "Pi0Flow_000167920.root",
   7                  const TString listPath = "PHOSPi0Flow/PHOSPi0FlowCoutput1",
   8                  const Int_t centrality=0,
   9                  const char* pid="CPV")
  10 {
  11   //Fit Real/Mixed ratio, normalize Mixed and subtract it from Real
  12
  13   TFile * f = new TFile(filename) ;
  14   //TList *histoList = (TList*)f->Get("PHOSPi0Flow");
  15   TList *histoList = (TList*)f->Get(listPath); // lego train
  16
  17   char key[125] ;
  18
  19   TH1F * hev = (TH1F*)histoList->FindObject("hTotSelEvents") ;
  20   TH2F * hCentrality  = (TH2F*)histoList->FindObject("hCenPHOSCells") ;
  21   TH1D * hCentrality1 = hCentrality->ProjectionX();
  22
  23   printf("TotSelEvents: %f \n",hev->GetBinContent(4)) ;
  24   printf("Centrality:   %f \n",hCentrality1->Integral()) ;
  25
  26   TString inputKey;
  27   TString outputKey = Form("%s_cent%d",pid,centrality);
  28
  29   TH2F *h , *hAdd;
  30   TH2F *hm, *hmAdd;
  31
  32   if (centrality >= 0 && centrality < 4) {
  33     printf("\tCentrality %d\n",centrality);
  34     inputKey = Form("hPi0%s_cen%d"  ,pid,centrality);
  35     TH2F *h = (TH2F*)histoList->FindObject(inputKey) ;
  36     inputKey = Form("hMiPi0%s_cen%d",pid,centrality);
  37     TH2F *hm= (TH2F*)histoList->FindObject(inputKey) ;
  38     if (h==0) {
  39       printf("Missing histogram %s\n",inputKey);
  40       return;
  41     }
  42   }
  43   else {
  44     printf("Wrong centrality %d. Allowed values are 0,1,2,3,10.\n",centrality);
  45     return;
  46   }
  47
  48   Int_t nPadX = 3, nPadY = 2;
  49   Int_t nPtBins=6 ;
  50   Double_t ptBinEdges[21]={1., 2., 3., 4., 5., 7., 10.} ;
  51
  52   PPRstyle();
  53   gStyle->SetPadLeftMargin(0.14);
  54   gStyle->SetPadRightMargin(0.01);
  55   //gStyle->SetPadTopMargin(0.01);
  56   gStyle->SetPadBottomMargin(0.08);
  57
  58   //Fit real only
  59   //Linear Bg
  60   char kkey[55];
  61   sprintf(kkey,outputKey.Data()) ;
  62   char key2[155];
  63   sprintf(key,"Mix%s",kkey) ;
  64   sprintf(key2,"%s_mr1",key) ;
  65   TH1D * mr1 = new TH1D(key2,"Mass",nPtBins,ptBinEdges) ;
  66   sprintf(key2,"%s_sr1",key) ;
  67   TH1D * sr1 = new TH1D(key2,"Width",nPtBins,ptBinEdges) ;
  68   sprintf(key2,"%s_ar1",key) ;
  69   TH1D * ar1 = new TH1D(key2,"a",nPtBins,ptBinEdges) ;
  70   sprintf(key2,"%s_br1",key) ;
  71   TH1D * br1 = new TH1D(key2,"a",nPtBins,ptBinEdges) ;
  72   sprintf(key2,"%s_yr1",key) ;
  73   TH1D * nr1 = new TH1D(key2,"Raw yield",nPtBins,ptBinEdges) ;
  74   sprintf(key2,"%s_yr1int",key) ;
  75   TH1D * nr1int = new TH1D(key2,"Raw yield, integrated",nPtBins,ptBinEdges) ;
  76
  77   //Quadratic Bg
  78   sprintf(key2,"%s_mr2",key) ;
  79   TH1D * mr2 = new TH1D(key2,"Mass",nPtBins,ptBinEdges) ;
  80   sprintf(key2,"%s_sr2",key) ;
  81   TH1D * sr2 = new TH1D(key2,"Width",nPtBins,ptBinEdges) ;
  82   sprintf(key2,"%s_ar2",key) ;
  83   TH1D * ar2 = new TH1D(key2,"a",nPtBins,ptBinEdges) ;
  84   sprintf(key2,"%s_br2",key) ;
  85   TH1D * br2 = new TH1D(key2,"a",nPtBins,ptBinEdges) ;
  86   sprintf(key2,"%s_cr2",key) ;
  87   TH1D * cr2 = new TH1D(key2,"a",nPtBins,ptBinEdges) ;
  88   sprintf(key2,"%s_yr2",key) ;
  89   TH1D * nr2 = new TH1D(key2,"Raw yield",nPtBins,ptBinEdges) ;
  90   sprintf(key2,"%s_yr2int",key) ;
  91   TH1D * nr2int = new TH1D(key2,"Raw yield, integrated",nPtBins,ptBinEdges) ;
  92
  93   TF1 * fit1 = new TF1("fit",CB,0.,1.,6) ;
  94   fit1->SetParName(0,"A") ;
  95   fit1->SetParName(1,"m_{0}") ;
  96   fit1->SetParName(2,"#sigma") ;
  97   fit1->SetParName(3,"a_{0}") ;
  98   fit1->SetParName(4,"a_{1}") ;
  99   fit1->SetParName(5,"a_{2}") ;
 100   fit1->SetLineWidth(2) ;
 101   fit1->SetLineColor(2) ;
 102   TF1 * fgs = new TF1("gs",CBs,0.,1.,4) ;
 103   fgs->SetParName(0,"A") ;
 104   fgs->SetParName(1,"m_{0}") ;
 105   fgs->SetParName(2,"#sigma") ;
 106   fgs->SetParName(3,"B") ;
 107   fgs->SetLineColor(2) ;
 108   fgs->SetLineWidth(1) ;
 109
 110   TF1 * fit2 = new TF1("fit2",CB2,0.,1.,7) ;
 111   fit2->SetParName(0,"A") ;
 112   fit2->SetParName(1,"m_{0}") ;
 113   fit2->SetParName(2,"#sigma") ;
 114   fit2->SetParName(3,"a_{0}") ;
 115   fit2->SetParName(4,"a_{1}") ;
 116   fit2->SetParName(5,"a_{2}") ;
 117   fit2->SetParName(6,"a_{3}") ;
 118   fit2->SetLineWidth(2) ;
 119   fit2->SetLineColor(4) ;
 120   fit2->SetLineStyle(2) ;
 121
 122   TF1 * fbg1 = new TF1("bg1",BG1,0.,1.,3) ;
 123   TF1 * fbg2 = new TF1("bg2",BG2,0.,1.,4) ;
 124
 125   TCanvas * c3 = new TCanvas("mggFit1_Signal","mggFitCB",10,10,1200,800) ;
 126   c3->Divide(nPadX,nPadY) ;
 127
 128   TCanvas * c1 = new TCanvas("mggFit1","mggFit1",10,10,1200,800) ;
 129   c1->Divide(nPadX,nPadY) ;
 130   c1->cd(0) ;
 131
 132   TCanvas * rawCanvas = new TCanvas("rawCanvas","rawCanvas",10,10,1200,800);
 133   rawCanvas->Divide(nPadX, nPadY);
 134
 135   TAxis * pta=h->GetYaxis() ;
 136   TAxis * ma=h->GetXaxis() ;
 137   for(Int_t ptBin=1; ptBin<=nPtBins; ptBin++){
 138     c1->cd(ptBin) ;
 139     printf("\n\t%.1f<pt<%.1f GeV/c\n",ptBinEdges[ptBin-1],ptBinEdges[ptBin]);
 140     Int_t imin=pta->FindBin(ptBinEdges[ptBin-1]+0.0001);
 141     Int_t imax=pta->FindBin(ptBinEdges[ptBin]-0.0001) ;
 142     Double_t pt=(ptBinEdges[ptBin]+ptBinEdges[ptBin-1])/2. ;
 143     TH1D * hp = h->ProjectionX(Form("re_%d",ptBin),imin,imax) ;
 144     hp->Sumw2() ;
 145     TH1D * hpm= hm->ProjectionX("mi",imin,imax) ;
 146     hpm->Sumw2() ;
 147     if(ptBin<=17){
 148       hp ->Rebin(4) ;
 149       hpm->Rebin(4) ;
 150     }
 151     else{
 152       hp ->Rebin(5) ;
 153       hpm->Rebin(5) ;
 154     }
 155     for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hp ->GetBinContent(ib)==0)hp ->SetBinError(ib,1.);}
 156     for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hpm->GetBinContent(ib)==0)hpm->SetBinError(ib,1.);}
 157     TH1D * hpm2   = (TH1D*)hpm->Clone("Bg1") ;
 158     TH1D * hpmScaled = (TH1D*)hpm->Clone("hpmScaled") ;
 159     TH1D * hpcopy = (TH1D*)hp ->Clone("hpcopy") ;
 160     TH1D * hp2    = (TH1D*)hp ->Clone("hp2") ;
 161     hpcopy->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})");
 162     hp2   ->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})");
 163     hpcopy->Divide(hpm) ;
 164     sprintf(key,"PID=%s, %3.1f<p_{T}<%3.1f GeV/c",pid,ptBinEdges[ptBin-1],ptBinEdges[ptBin]) ;
 165     hpcopy->SetTitle(key) ;
 166     hpcopy->SetMarkerStyle(20) ;
 167     hpcopy->SetMarkerSize(0.7) ;
 168
 169     fit1->SetParameters(0.001,0.136,0.0055,0.0002,-0.002,0.0) ;
 170     fit1->SetParLimits(0,0.000,1.000) ;
 171     fit1->SetParLimits(1,0.120,0.145) ;
 172     fit1->SetParLimits(2,0.005,0.012) ;
 173
 174     Double_t rangeMin=0.05 ;
 175     Double_t rangeMax=0.3 ;
 176     if(centrality==0) rangeMax=0.4 ;
 177     if(ptBin==1){
 178       rangeMin=0.06 ;
 179       rangeMax=0.25 ;
 180     }
 181     hpcopy->Fit(fit1,"Q" ,"",rangeMin,rangeMax) ;
 182     hpcopy->Fit(fit1,"MQ","",rangeMin,rangeMax) ;
 183
 184     ar1->SetBinContent(ptBin,fit1->GetParameter(3)) ;
 185     ar1->SetBinError  (ptBin,fit1->GetParError(3)) ;
 186     br1->SetBinContent(ptBin,fit1->GetParameter(4)) ;
 187     br1->SetBinError  (ptBin,fit1->GetParError(4)) ;
 188
 189     fit2->SetParameters(fit1->GetParameters()) ;
 190     fit2->SetParLimits(0,0.000,1.000) ;
 191     fit2->SetParLimits(1,0.120,0.145) ;
 192     fit2->SetParLimits(2,0.005,0.012) ;
 193
 194     hpcopy->Fit(fit2,"+NQ","",rangeMin,rangeMax) ;
 195     hpcopy->Fit(fit2,"+MQ","",rangeMin,rangeMax) ;
 196
 197     ar2->SetBinContent(ptBin,fit2->GetParameter(3)) ;
 198     ar2->SetBinError  (ptBin,fit2->GetParError(3)) ;
 199     br2->SetBinContent(ptBin,fit2->GetParameter(4)) ;
 200     br2->SetBinError  (ptBin,fit2->GetParError(4)) ;
 201     cr2->SetBinContent(ptBin,fit2->GetParameter(5)) ;
 202     cr2->SetBinError  (ptBin,fit2->GetParError(5)) ;
 203     hpcopy->GetXaxis()->SetRangeUser(0.05,0.30) ;
 204     hpcopy->Draw() ;
 205     c1->Update() ;
 206
 207     fbg1->SetParameters(fit1->GetParameter(3),fit1->GetParameter(4),fit1->GetParameter(5));
 208     fbg2->SetParameters(fit2->GetParameter(3),fit2->GetParameter(4),fit2->GetParameter(5),
 209                         fit2->GetParameter(6));
 210
 211     Double_t intRangeMin = PeakPosition(pt)-3.*PeakWidth(pt) ;
 212     Double_t intRangeMax = PeakPosition(pt)+3.*PeakWidth(pt) ;
 213     Int_t    intBinMin   = hp->GetXaxis()->FindBin(intRangeMin) ;
 214     Int_t    intBinMax   = hp->GetXaxis()->FindBin(intRangeMax) ;
 215     Double_t errStat     = hpm->Integral(intBinMin,intBinMax);
 216
 217     rawCanvas->cd(ptBin);
 218     hpmScaled->Scale(fbg1(0.1349));
 219     hp->SetTitle(key);
 220     hp->SetLineColor(kBlack);
 221     hp->SetAxisRange(0.05, 0.3);
 222     hp->DrawCopy();
 223     hpmScaled->SetLineColor(kRed);
 224     hpmScaled->DrawCopy("same");
 225     rawCanvas->Update();
 226
 227     hpm ->Multiply(fbg1) ;
 228     hpm2->Multiply(fbg2) ;
 229     hp  ->Add(hpm ,-1.) ;
 230     hp2 ->Add(hpm2,-1.) ;
 231
 232     c3->cd(ptBin) ;
 233
 234     Int_t binPi0 = hp->FindBin(kMean);
 235     fgs->SetParameters(hp->Integral(binPi0-1,binPi0+1)/3.,fit1->GetParameter(1),fit1->GetParameter(2)) ;
 236     fgs->SetParLimits(0,0.000,1.e+5) ;
 237     fgs->SetParLimits(1,0.120,0.145) ;
 238     fgs->SetParLimits(2,0.004,0.02) ;
 239     hp->Fit(fgs,"Q","",rangeMin,rangeMax) ;
 240     hp->SetMaximum(hp2->GetMaximum()*1.5) ;
 241     hp->SetMinimum(hp2->GetMinimum()*1.1) ;
 242     hp->SetMarkerStyle(20) ;
 243     hp->SetMarkerSize(0.7) ;
 244     mr1->SetBinContent(ptBin,fgs->GetParameter(1)) ;
 245     mr1->SetBinError  (ptBin,fgs->GetParError(1) ) ;
 246     sr1->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ;
 247     sr1->SetBinError  (ptBin,fgs->GetParError(2) ) ;
 248
 249     Double_t y=fgs->GetParameter(0)/hp->GetXaxis()->GetBinWidth(1) ;
 250     nr1->SetBinContent(ptBin,y) ;
 251     Double_t ey=fgs->GetParError(0)/hp->GetXaxis()->GetBinWidth(1) ;
 252     nr1->SetBinError(ptBin,ey) ;
 253
 254     Double_t npiInt = hp->Integral(intBinMin,intBinMax) ;
 255     Double_t norm   = fbg1->GetParameter(0) ;
 256     Double_t normErr= fbg1->GetParError(0) ;
 257     if(npiInt>0.){
 258       nr1int->SetBinContent(ptBin,npiInt) ;
 259       nr1int->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
 260     }
 261     hp2->GetXaxis()->SetRangeUser(0.05,0.3) ;
 262     hp2->SetMaximum(hp2->GetMaximum()*1.5) ;
 263     hp2->SetMinimum(hp2->GetMinimum()*1.1) ;
 264     hp2->SetMarkerStyle(20) ;
 265     hp2->SetMarkerSize(0.7) ;
 266     hp2->Fit(fgs,"Q","",rangeMin,rangeMax) ;
 267     mr2->SetBinContent(ptBin,fgs->GetParameter(1)) ;
 268     mr2->SetBinError  (ptBin,fgs->GetParError(1)) ;
 269     sr2->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ;
 270     sr2->SetBinError  (ptBin,fgs->GetParError(2)) ;
 271     y=fgs->GetParameter(0)/hp->GetXaxis()->GetBinWidth(1) ;
 272     nr2->SetBinContent(ptBin,y) ;
 273     ey= fgs->GetParError(0)/hp->GetXaxis()->GetBinWidth(1) ;
 274     nr2->SetBinError(ptBin,ey) ;
 275     npiInt=hp2->Integral(intBinMin,intBinMax) ;
 276     norm=fbg2->GetParameter(0) ;
 277     normErr=fbg2->GetParError(0) ;
 278     if(npiInt>0.){
 279       nr2int->SetBinContent(ptBin,npiInt) ;
 280       nr2int->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
 281     }
 282     hp2->SetTitle(key) ;
 283     hp2->Draw() ;
 284     c3->Update() ;
 285
 286     delete hp ;
 287     delete hpm ;
 288     delete hpm2 ;
 289   }
 290   char name[55] ;
 291
 292   sprintf(name,"Pi0_ratio_%s.eps",kkey) ;
 293   c1->Print(name) ;
 294   sprintf(name,"Pi0_signal_%s.eps",kkey) ;
 295   c3->Print(name) ;
 296
 297   //Normalize by the number of events
 298   Int_t cMin,cMax;
 299   if      (centrality == 0) {
 300     cMin=1;
 301     cMax=10;
 302   }
 303   else if (centrality == 1) {
 304     cMin=21;
 305     cMax=40;
 306   }
 307   else if (centrality == 2) {
 308     cMin=41;
 309     cMax=60;
 310   }
 311   else if (centrality == 3) {
 312     cMin=61;
 313     cMax=80;
 314   }
 315   Double_t nevents = hCentrality1->Integral(cMin,cMax);
 316   nr1   ->Scale(1./nevents) ;
 317   nr1int->Scale(1./nevents) ;
 318   nr2   ->Scale(1./nevents) ;
 319   nr2int->Scale(1./nevents) ;
 320   printf("Nevents=%f \n",nevents) ;
 321
 322
 323   TFile fout("LHC10h_Pi0_FitResult.root","update");
 324   mr1->Write(0,TObject::kOverwrite) ;
 325   sr1->Write(0,TObject::kOverwrite) ;
 326   ar1->Write(0,TObject::kOverwrite) ;
 327   br1->Write(0,TObject::kOverwrite) ;
 328   nr1->Write(0,TObject::kOverwrite) ;
 329   nr1int->Write(0,TObject::kOverwrite) ;
 330   ar2->Write(0,TObject::kOverwrite) ;
 331   br2->Write(0,TObject::kOverwrite) ;
 332   cr2->Write(0,TObject::kOverwrite) ;
 333   mr2->Write(0,TObject::kOverwrite) ;
 334   sr2->Write(0,TObject::kOverwrite) ;
 335   nr2->Write(0,TObject::kOverwrite) ;
 336   nr2int->Write(0,TObject::kOverwrite) ;
 337   fout.Close() ;
 338
 339 }
 340
 341 //-----------------------------------------------------------------------------
 342 Double_t PeakPosition(Double_t pt){
 343   //Fit to the measured peak position
 344   return 4.99292e-003*exp(-pt*9.32300e-001)+1.34944e-001 ;
 345 }
 346 //-----------------------------------------------------------------------------
 347 Double_t PeakWidth(Double_t pt){
 348   //fit to the measured peak width
 349   Double_t a=0.0068 ;
 350   Double_t b=0.0025 ;
 351   Double_t c=0.000319 ;
 352   return TMath::Sqrt(a*a+b*b/pt/pt+c*c*pt*pt) ;
 353 }
 354
 355 //-----------------------------------------------------------------------------
 356 Double_t CB(Double_t * x, Double_t * par){
 357   //Parameterization of Real/Mixed ratio
 358   Double_t m=par[1] ;
 359   Double_t s=par[2] ;
 360   Double_t dx=(x[0]-m)/s ;
 361   return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean) ;
 362 }
 363
 364 //-----------------------------------------------------------------------------
 365 Double_t CB2(Double_t * x, Double_t * par){
 366   //Another parameterization of Real/Mixed ratio7TeV/README
 367   Double_t m=par[1] ;
 368   Double_t s=par[2] ;
 369   Double_t dx=(x[0]-m)/s ;
 370   return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean)+par[5]*(x[0]-kMean)*(x[0]-kMean) ;
 371 }
 372 //-----------------------------------------------------------------------------
 373 Double_t CBs(Double_t * x, Double_t * par){
 374   //Parameterizatin of signal
 375   Double_t m=par[1] ;
 376   Double_t s=par[2] ;
 377   Double_t dx=(x[0]-m)/s ;
 378   return par[0]*exp(-dx*dx/2.)/TMath::Sqrt(TMath::TwoPi())/s+par[3] ;
 379 }
 380 //-----------------------------------------------------------------------------
 381 Double_t BG1(Double_t * x, Double_t * par){
 382   //Normalizatino of Mixed
 383   return par[0]+par[1]*(x[0]-kMean) ;
 384 }
 385 //-----------------------------------------------------------------------------
 386 Double_t BG2(Double_t * x, Double_t * par){
 387   //Another normalization of  Mixed
 388   return par[0]+par[1]*(x[0]-kMean)+par[2]*(x[0]-kMean)*(x[0]-kMean) ;
 389 }
 390
 391
 392 //-----------------------------------------------------------------------------
 393 PPRstyle()
 394 {
 395
 396   //////////////////////////////////////////////////////////////////////
 397   //
 398   // ROOT style macro for the TRD TDR
 399   //
 400   //////////////////////////////////////////////////////////////////////
 401
 402   gStyle->SetPalette(1);
 403   gStyle->SetCanvasBorderMode(-1);
 404   gStyle->SetCanvasBorderSize(1);
 405   gStyle->SetCanvasColor(10);
 406
 407   gStyle->SetFrameFillColor(10);
 408   gStyle->SetFrameBorderSize(1);
 409   gStyle->SetFrameBorderMode(-1);
 410   gStyle->SetFrameLineWidth(1.2);
 411   gStyle->SetFrameLineColor(1);
 412
 413   gStyle->SetHistFillColor(0);
 414   gStyle->SetHistLineWidth(1);
 415   gStyle->SetHistLineColor(1);
 416
 417   gStyle->SetPadColor(10);
 418   gStyle->SetPadBorderSize(1);
 419   gStyle->SetPadBorderMode(-1);
 420
 421   gStyle->SetStatColor(10);
 422   gStyle->SetTitleColor(kBlack,"X");
 423   gStyle->SetTitleColor(kBlack,"Y");
 424
 425   gStyle->SetLabelSize(0.04,"X");
 426   gStyle->SetLabelSize(0.04,"Y");
 427   gStyle->SetLabelSize(0.04,"Z");
 428   gStyle->SetTitleSize(0.04,"X");
 429   gStyle->SetTitleSize(0.04,"Y");
 430   gStyle->SetTitleSize(0.04,"Z");
 431   gStyle->SetTitleFont(42,"X");
 432   gStyle->SetTitleFont(42,"Y");
 433   gStyle->SetTitleFont(42,"X");
 434   gStyle->SetLabelFont(42,"X");
 435   gStyle->SetLabelFont(42,"Y");
 436   gStyle->SetLabelFont(42,"Z");
 437   gStyle->SetStatFont(42);
 438
 439   gStyle->SetTitleOffset(1.0,"X");
 440   gStyle->SetTitleOffset(1.4,"Y");
 441
 442   gStyle->SetFillColor(kWhite);
 443   gStyle->SetTitleFillColor(kWhite);
 444
 445   gStyle->SetOptDate(0);
 446   gStyle->SetOptTitle(1);
 447   gStyle->SetOptStat(0);
 448   gStyle->SetOptFit(111);
 449
 450 }
 451