PWGGA/PHOSTasks/PHOS_pp_pi0/macros/makeMmixPU_GS.C

   1 //-----------------------------------------------------------------------------
   2 void makeMmixPU_GS(const TString histoFile="LHC11a_pass4_20130913.root",
   3                    const Int_t nSigma=2,
   4                    const char* module="A10")
   5 {
   6   //Fit Real/Mixed ratio, normalize Mixed and subtract it from Real.
   7   // The pi0 peak if fitted by the Gaussian function,
   8   // the background is fitted by pol1 or pol2
   9
  10   TString hMassName;
  11   TFile * file = new TFile(histoFile) ;
  12   THashList *hList = (THashList*)file->Get("histESD");
  13   char key[125] ;
  14   hMassName = "hMassPt";
  15   hMassName += module;
  16   TH2F * h   = (TH2F*)hList->FindObject(hMassName) ;
  17
  18   hMassName = "hMiMassPt";
  19   hMassName += module;
  20   TH2F * hm  = (TH2F*)hList->FindObject(hMassName) ;
  21
  22   TH1F * hev = (TH1F*)hList->FindObject("hSelEvents") ;
  23
  24   // Array of pt bins
  25   Int_t nbin=18 ;
  26   Double_t xa[]={0.6,0.8,1.,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,3.0,3.5,4.0,5.,6.,8.,10.,12.} ;
  27   PPRstyle();
  28   gStyle->SetPadLeftMargin(0.16);
  29   gStyle->SetPadRightMargin(0.01);
  30   gStyle->SetPadTopMargin(0.02);
  31   gStyle->SetPadBottomMargin(0.11);
  32   gStyle->SetTitleX(0.80);
  33   gStyle->SetTitleY(0.98);
  34
  35   //Fit real only
  36   //Linear Bg
  37   char key2[155];
  38   sprintf(key,"Mix") ;
  39   sprintf(key2,"%s_mr1",key) ;
  40   TH1D * mr1 = new TH1D(key2,"Mass",nbin,xa) ;
  41   sprintf(key2,"%s_sr1",key) ;
  42   TH1D * sr1 = new TH1D(key2,"Width",nbin,xa) ;
  43   sprintf(key2,"%s_ar1",key) ;
  44   TH1D * ar1 = new TH1D(key2,"a",nbin,xa) ;
  45   sprintf(key2,"%s_br1",key) ;
  46   TH1D * br1 = new TH1D(key2,"a",nbin,xa) ;
  47   sprintf(key2,"%s_yr1",key) ;
  48   TH1D * nr1 = new TH1D(key2,"Raw yield",nbin,xa) ;
  49   sprintf(key2,"%s_yr1int",key) ;
  50   TH1D * nr1int = new TH1D(key2,"Raw yield, integrated",nbin,xa) ;
  51
  52   //Quadratic Bg
  53   sprintf(key2,"%s_mr2",key) ;
  54   TH1D * mr2 = new TH1D(key2,"Mass",nbin,xa) ;
  55   sprintf(key2,"%s_sr2",key) ;
  56   TH1D * sr2 = new TH1D(key2,"Width",nbin,xa) ;
  57   sprintf(key2,"%s_ar2",key) ;
  58   TH1D * ar2 = new TH1D(key2,"a",nbin,xa) ;
  59   sprintf(key2,"%s_br2",key) ;
  60   TH1D * br2 = new TH1D(key2,"a",nbin,xa) ;
  61   sprintf(key2,"%s_cr2",key) ;
  62   TH1D * cr2 = new TH1D(key2,"a",nbin,xa) ;
  63   sprintf(key2,"%s_yr2",key) ;
  64   TH1D * nr2 = new TH1D(key2,"Raw yield",nbin,xa) ;
  65   sprintf(key2,"%s_yr2int",key) ;
  66   TH1D * nr2int = new TH1D(key2,"Raw yield, integrated",nbin,xa) ;
  67
  68   TF1 * fit1 = new TF1("fit",CB,0.,1.,5) ;
  69   fit->SetParName(0,"A") ;
  70   fit->SetParName(1,"m_{0}") ;
  71   fit->SetParName(2,"#sigma") ;
  72   fit->SetParName(3,"a_{0}") ;
  73   fit->SetParName(4,"a_{1}") ;
  74   fit->SetLineWidth(2) ;
  75   fit->SetLineColor(2) ;
  76   TF1 * fgs = new TF1("gs",CBs,0.,1.,3) ;
  77   fgs->SetLineColor(2) ;
  78   fgs->SetLineWidth(1) ;
  79
  80   TF1 * fit2 = new TF1("fit2",CB2,0.,1.,6) ;
  81   fit2->SetParName(0,"A") ;
  82   fit2->SetParName(1,"m_{0}") ;
  83   fit2->SetParName(2,"#sigma") ;
  84   fit2->SetParName(3,"a_{0}") ;
  85   fit2->SetParName(4,"a_{1}") ;
  86   fit2->SetLineWidth(2) ;
  87   fit2->SetLineColor(4) ;
  88   fit2->SetLineStyle(2) ;
  89   TF1 * fbg1 = new TF1("bg1",BG1,0.,1.,2) ;
  90   TF1 * fbg2 = new TF1("bg2",BG2,0.,1.,3) ;
  91
  92   TCanvas * c3 = new TCanvas("mggFit1_Signal","mggFitCB",10,10,1400,800) ;
  93   c3->Divide(6,3) ;
  94
  95   TCanvas * cReal = new TCanvas("cReal","Mgg real events",10,10,1400,800) ;
  96   cReal->Divide(6,3) ;
  97
  98   TCanvas * c1 = new TCanvas("mggFit1","mggFit1",10,10,1400,800) ;
  99   c1->Divide(6,3) ;
 100   c1->cd(0) ;
 101   TCanvas * c2=0,*c4=0,*c5=0,*c6=0 ;
 102
 103   TAxis * pta=h->GetYaxis() ;
 104   TAxis * ma=h->GetXaxis() ;
 105   for(Int_t i=1;i<=nbin;i++){
 106     c1->cd(i) ;
 107     Int_t imin=pta->FindBin(xa[i-1]+0.0001);
 108     Int_t imax=pta->FindBin(xa[i]-0.0001) ;
 109     Double_t pt=(xa[i]+xa[i-1])/2. ;
 110     TH1D * hp = h->ProjectionX(Form("re%d",i),imin,imax) ;
 111     hp->Sumw2() ;
 112     TH1D * hpm= hm->ProjectionX(Form("mi%d",i),imin,imax) ;
 113     hpm->Sumw2() ;
 114     if(i<15){
 115       hp ->Rebin(2) ;
 116       hpm->Rebin(2) ;
 117     }
 118     else{
 119       hp ->Rebin(5) ;
 120       hpm->Rebin(5) ;
 121     }
 122     hp ->SetNdivisions(506);
 123     hpm->SetNdivisions(506);
 124     hp ->SetLabelSize(0.05,"X");
 125     hpm->SetLabelSize(0.05,"X");
 126     hp ->SetTitleSize(0.00,"X");
 127     hpm->SetTitleSize(0.00,"X");
 128     hp ->SetLabelSize(0.05,"Y");
 129     hpm->SetLabelSize(0.05,"Y");
 130     hp ->SetTitleSize(0.05,"Y");
 131     hpm->SetTitleSize(0.05,"Y");
 132     hp ->SetTitleOffset(1.0,"X");
 133     hpm->SetTitleOffset(1.0,"X");
 134     if (i>12) {
 135       hp ->SetLabelSize(0.05,"X");
 136       hpm->SetLabelSize(0.05,"X");
 137       hp ->SetTitleSize(0.05,"X");
 138       hpm->SetTitleSize(0.05,"X");
 139       hp ->SetXTitle("m_{#gamma#gamma} (GeV/c^{2})");
 140       hpm->SetXTitle("m_{#gamma#gamma} (GeV/c^{2})");
 141     }
 142
 143     // for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hp ->GetBinContent(ib)==0)hp ->SetBinError(ib,1.);}
 144     // for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hpm->GetBinContent(ib)==0)hpm->SetBinError(ib,1.);}
 145     TH1D * hpm2   = (TH1D*)hpm->Clone("Bg1") ;
 146     TH1D * hpcopy = (TH1D*)hp ->Clone("hpcopy") ;
 147     TH1D * hp2    = (TH1D*)hp ->Clone("hp2") ;
 148     TH1D * hpReal = (TH1D*)hp ->Clone("hpReal") ;
 149     hpcopy->Divide(hpm) ;
 150     sprintf(key,"%3.1f<p_{t}<%3.1f GeV/c",xa[i-1],xa[i]) ;
 151     hpcopy->SetTitle(key) ;
 152     hpReal->SetTitle(key) ;
 153     hpReal->SetLineWidth(2);
 154     hpcopy->SetMarkerStyle(20) ;
 155     hpcopy->SetMarkerSize(0.7) ;
 156
 157     fit1->SetParameters(TMath::Min(0.3,0.001*i*i),0.136,0.0077,0.0013,-0.0007) ;
 158     fit1->SetParLimits(2,0.003,0.010) ;
 159     fit1->SetParLimits(1,0.132,0.139) ;
 160     fit1->SetParLimits(0,0.,1.) ;
 161
 162     Double_t rangeMin=0.06 ;
 163     Double_t rangeMax=0.22 ;
 164     if(i>=16)rangeMax=0.3 ; //More points to fix background
 165     hpcopy->Fit(fit1,"NQ","",rangeMin,rangeMax) ;
 166     hpcopy->Fit(fit1,"MQ","",rangeMin,rangeMax) ;
 167     // hpcopy->Fit(fit1,"LQ","",rangeMin,rangeMax) ;
 168     ar1->SetBinContent(i,fit1->GetParameter(3)) ;
 169     ar1->SetBinError  (i,fit1->GetParError (3)) ;
 170     br1->SetBinContent(i,fit1->GetParameter(4)) ;
 171     br1->SetBinError  (i,fit1->GetParError (4)) ;
 172
 173     fit2->SetParameters(fit1->GetParameters()) ;
 174     fit2->SetParLimits(2,0.003,0.010) ;
 175     fit2->SetParLimits(1,0.132,0.139) ;
 176     fit2->SetParLimits(0,0.,1.) ;
 177     fit2->SetParameter(5,0.) ;
 178
 179     hpcopy->Fit(fit2,"+NQ","",rangeMin,rangeMax) ;
 180     hpcopy->Fit(fit2,"+MQ","",rangeMin,rangeMax) ;
 181     // hpcopy->Fit(fit2,"+LQ","",rangeMin,rangeMax) ;
 182     ar2->SetBinContent(i,fit2->GetParameter(3)) ;
 183     ar2->SetBinError  (i,fit2->GetParError (3)) ;
 184     br2->SetBinContent(i,fit2->GetParameter(4)) ;
 185     br2->SetBinError  (i,fit2->GetParError (4)) ;
 186     cr2->SetBinContent(i,fit2->GetParameter(5)) ;
 187     cr2->SetBinError  (i,fit2->GetParError (5)) ;
 188
 189     c1->cd(i) ;
 190     hpcopy->SetAxisRange(0.065,0.229,"X") ;
 191     hpcopy->Draw() ;
 192
 193     cReal->cd(i) ;
 194     hpReal->SetAxisRange(0.065,0.229,"X") ;
 195     hpReal->Draw("ehist") ;
 196
 197     fbg1->SetParameters(fit1->GetParameter(3),fit1->GetParameter(4));
 198     fbg2->SetParameters(fit2->GetParameter(3),fit2->GetParameter(4),fit2->GetParameter(5));
 199
 200     Double_t intRangeMin = PeakPosition(pt)-nSigma*PeakWidth(pt) ;
 201     Double_t intRangeMax = PeakPosition(pt)+nSigma*PeakWidth(pt) ;
 202     Int_t    intBinMin   = hp->GetXaxis()->FindBin(intRangeMin) ;
 203     Int_t    intBinMax   = hp->GetXaxis()->FindBin(intRangeMax) ;
 204     Double_t errStat     = hpm->Integral(intBinMin,intBinMax);
 205
 206     hpm ->Multiply(fbg1) ;
 207     hpm2->Multiply(fbg2) ;
 208     hp  ->Add(hpm,-1.) ;
 209     hp2 ->Add(hpm2,-1.) ;
 210
 211     c3->cd(i) ;
 212
 213     if(i<15)
 214       fgs->SetParameters(hp->Integral(32,36)/5.,fit1->GetParameter(1),fit1->GetParameter(2)) ;
 215     else
 216       fgs->SetParameters(hp->Integral(13,15)/3.,0.136,0.005) ;
 217     fgs->SetParLimits(2,0.003,0.010) ;
 218     fgs->SetParLimits(1,0.132,0.142) ;
 219     fgs->SetParLimits(0,0.,1.e+6) ;
 220
 221     hp->Fit(fgs,"Q","",rangeMin,rangeMax) ;
 222     hp->SetMaximum(hp2->GetMaximum()*1.1) ;
 223     hp->SetMinimum(hp2->GetMinimum()*1.1) ;
 224     hp->SetMarkerStyle(20) ;
 225     hp->SetMarkerSize(0.7) ;
 226     mr1->SetBinContent(i,fgs->GetParameter(1)) ;
 227     mr1->SetBinError  (i,fgs->GetParError(1)) ;
 228     sr1->SetBinContent(i,TMath::Abs(fgs->GetParameter(2))) ;
 229     sr1->SetBinError  (i,fgs->GetParError(2)) ;
 230
 231     Double_t y=fgs->Integral(intRangeMin,intRangeMax)/hp->GetXaxis()->GetBinWidth(1) ;
 232     nr1->SetBinContent(i,y) ;
 233     Double_t ey=0 ;
 234     if(fgs->GetParameter(0)!=0. && fgs->GetParameter(2)!=0.){
 235       Double_t en=fgs->GetParError(0)/fgs->GetParameter(0) ;
 236       Double_t es=fgs->GetParError(2)/fgs->GetParameter(2) ;
 237       ey=y*TMath::Sqrt(en*en+es*es) ;
 238     }
 239     nr1->SetBinError(i,ey) ;
 240
 241     Double_t npiInt=hp->Integral(intBinMin,intBinMax) ;
 242     Double_t norm=fbg1->GetParameter(0) ;
 243     Double_t normErr=fbg1->GetParError(0) ;
 244     if(npiInt>0.){
 245       nr1int->SetBinContent(i,npiInt) ;
 246       nr1int->SetBinError(i,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
 247     }
 248     hp2->GetXaxis()->SetRangeUser(0.06,0.229) ;
 249     hp2->SetMaximum(hp2->GetMaximum()*1.1) ;
 250     hp2->SetMinimum(hp2->GetMinimum()*1.1) ;
 251     hp2->SetMarkerStyle(20) ;
 252     hp2->SetMarkerSize(0.7) ;
 253
 254     hp2->Fit(fgs,"Q","",rangeMin,rangeMax) ;
 255     mr2->SetBinContent(i,fgs->GetParameter(1)) ;
 256     mr2->SetBinError  (i,fgs->GetParError(1)) ;
 257     sr2->SetBinContent(i,TMath::Abs(fgs->GetParameter(2))) ;
 258     sr2->SetBinError  (i,fgs->GetParError(2)) ;
 259     y=fgs->Integral(intRangeMin,intRangeMax)/hp->GetXaxis()->GetBinWidth(1) ;
 260     nr2->SetBinContent(i,y) ;
 261     ey=0 ;
 262     if(fgs->GetParameter(0)!=0. && fgs->GetParameter(2)!=0.){
 263       Double_t en=fgs->GetParError(0)/fgs->GetParameter(0) ;
 264       Double_t es=fgs->GetParError(2)/fgs->GetParameter(2) ;
 265       ey=y*TMath::Sqrt(en*en+es*es) ;
 266     }
 267     nr2->SetBinError(i,ey) ;
 268     npiInt=hp2->Integral(intBinMin,intBinMax) ;
 269     norm=fbg2->GetParameter(0) ;
 270     normErr=fbg2->GetParError(0) ;
 271     if(npiInt>0.){
 272       nr2int->SetBinContent(i,npiInt) ;
 273       nr2int->SetBinError(i,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
 274     }
 275     hp2->SetTitle(key) ;
 276     hp2->SetAxisRange(0.065,0.229,"X") ;
 277     hp2->Draw() ;
 278     hp ->Draw("same") ;
 279     delete hpm ;
 280     delete hpm2 ;
 281     c1->Update() ;
 282     c3->Update() ;
 283     cReal->Update() ;
 284   }
 285
 286   if (c3) c3->Print("Pi0_Signal_Gaus.eps") ;
 287   if (c1) c1->Print("Pi0_Ratio_Gaus.eps") ;
 288   if (cReal) cReal->Print("Pi0_InvMass_Gaus.eps") ;
 289
 290   //Normalize by the number of events
 291   Double_t nMBOR   = hev->GetBinContent(2); // MBOR events pass4
 292   Double_t nevents = nMBOR;
 293   printf("==============\nN events = %d\n==============\n",nevents);
 294
 295   nr1   ->Scale(1./nevents) ;
 296   nr1int->Scale(1./nevents) ;
 297   nr2   ->Scale(1./nevents) ;
 298   nr2int->Scale(1./nevents) ;
 299
 300
 301   TFile fout("LHC11a_FitResult.root","update");
 302   mr1->Write(0,TObject::kOverwrite) ;
 303   sr1->Write(0,TObject::kOverwrite) ;
 304   ar1->Write(0,TObject::kOverwrite) ;
 305   br1->Write(0,TObject::kOverwrite) ;
 306   nr1->Write(0,TObject::kOverwrite) ;
 307   nr1int->Write(0,TObject::kOverwrite) ;
 308   ar2->Write(0,TObject::kOverwrite) ;
 309   br2->Write(0,TObject::kOverwrite) ;
 310   cr2->Write(0,TObject::kOverwrite) ;
 311   mr2->Write(0,TObject::kOverwrite) ;
 312   sr2->Write(0,TObject::kOverwrite) ;
 313   nr2->Write(0,TObject::kOverwrite) ;
 314   nr2int->Write(0,TObject::kOverwrite) ;
 315   fout.Close() ;
 316
 317 }
 318
 319 //-----------------------------------------------------------------------------
 320 const Double_t kMean=0.135 ; //Approximate peak position to facilitate error estimate
 321
 322 Double_t PeakPosition(Double_t pt){
 323   //Fit to the measured peak position
 324   return 0.454*exp(-pt*7.13)+0.13569 ;
 325 }
 326 Double_t PeakWidth(Double_t pt){
 327   //fit to the measured peak width
 328   return 1.69e-02*exp(-pt*2.82)+4.93e-03 ;
 329 }
 330
 331 Double_t CB(Double_t * x, Double_t * par){
 332   //Parameterization of Real/Mixed ratio
 333   Double_t m=par[1] ;
 334   Double_t s=par[2] ;
 335   Double_t dx=(x[0]-m)/s ;
 336   return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean) ;
 337 }
 338 Double_t CB2(Double_t * x, Double_t * par){
 339   //Another parameterization of Real/Mixed ratio
 340   Double_t m=par[1] ;
 341   Double_t s=par[2] ;
 342   Double_t dx=(x[0]-m)/s ;
 343   return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean)+par[5]*(x[0]-kMean)*(x[0]-kMean) ;
 344 }
 345 Double_t CBs(Double_t * x, Double_t * par){
 346   //Parameterizatin of signal
 347   Double_t m=par[1] ;
 348   Double_t s=par[2] ;
 349   Double_t dx=(x[0]-m)/s ;
 350   return par[0]*exp(-dx*dx/2.) ;
 351 }
 352 Double_t BG1(Double_t * x, Double_t * par){
 353   //Normalizatino of Mixed
 354   return par[0]+par[1]*(x[0]-kMean) ;
 355 }
 356 Double_t BG2(Double_t * x, Double_t * par){
 357   //Another normalization of  Mixed
 358   return par[0]+par[1]*(x[0]-kMean)+par[2]*(x[0]-kMean)*(x[0]-kMean) ;
 359 }
 360
 361
 362 //-----------------------------------------------------------------------------
 363 PPRstyle()
 364 {
 365
 366   //////////////////////////////////////////////////////////////////////
 367   //
 368   // ROOT style macro for the TRD TDR
 369   //
 370   //////////////////////////////////////////////////////////////////////
 371
 372   gStyle->SetPalette(1);
 373   gStyle->SetCanvasBorderMode(-1);
 374   gStyle->SetCanvasBorderSize(1);
 375   gStyle->SetCanvasColor(10);
 376
 377   gStyle->SetFrameFillColor(10);
 378   gStyle->SetFrameBorderSize(1);
 379   gStyle->SetFrameBorderMode(-1);
 380   gStyle->SetFrameLineWidth(1.2);
 381   gStyle->SetFrameLineColor(1);
 382
 383   gStyle->SetHistFillColor(0);
 384   gStyle->SetHistLineWidth(1);
 385   gStyle->SetHistLineColor(1);
 386
 387   gStyle->SetPadColor(10);
 388   gStyle->SetPadBorderSize(1);
 389   gStyle->SetPadBorderMode(-1);
 390
 391   gStyle->SetStatColor(10);
 392   gStyle->SetTitleColor(kBlack,"X");
 393   gStyle->SetTitleColor(kBlack,"Y");
 394
 395   gStyle->SetLabelSize(0.04,"X");
 396   gStyle->SetLabelSize(0.04,"Y");
 397   gStyle->SetLabelSize(0.04,"Z");
 398   gStyle->SetTitleSize(0.04,"X");
 399   gStyle->SetTitleSize(0.04,"Y");
 400   gStyle->SetTitleSize(0.04,"Z");
 401   gStyle->SetTitleFont(42,"X");
 402   gStyle->SetTitleFont(42,"Y");
 403   gStyle->SetTitleFont(42,"X");
 404   gStyle->SetLabelFont(42,"X");
 405   gStyle->SetLabelFont(42,"Y");
 406   gStyle->SetLabelFont(42,"Z");
 407   gStyle->SetStatFont(42);
 408
 409   gStyle->SetTitleOffset(1.0,"X");
 410   gStyle->SetTitleOffset(1.4,"Y");
 411
 412   gStyle->SetFillColor(kWhite);
 413   gStyle->SetTitleFillColor(kWhite);
 414
 415   gStyle->SetOptDate(0);
 416   gStyle->SetOptTitle(1);
 417   gStyle->SetOptStat(0);
 418   gStyle->SetOptFit(0);
 419
 420 }
 421