Macros for post-processing histograms produced by AliAnalysisTaskPi0
authorkharlov <kharlov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 25 Nov 2013 17:51:20 +0000 (17:51 +0000)
committerkharlov <kharlov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 25 Nov 2013 17:51:20 +0000 (17:51 +0000)
PWGGA/PHOSTasks/PHOS_pp_pi0/macros/MakeFinalSpectrum.C [new file with mode: 0644]
PWGGA/PHOSTasks/PHOS_pp_pi0/macros/makeMmixPU_CB.C [new file with mode: 0644]
PWGGA/PHOSTasks/PHOS_pp_pi0/macros/makeMmixPU_GS.C [new file with mode: 0644]

diff --git a/PWGGA/PHOSTasks/PHOS_pp_pi0/macros/MakeFinalSpectrum.C b/PWGGA/PHOSTasks/PHOS_pp_pi0/macros/MakeFinalSpectrum.C
new file mode 100644 (file)
index 0000000..6251635
--- /dev/null
@@ -0,0 +1,384 @@
+void MakeFinalSpectrum()\r
+{\r
+  //-----------------------------------------------------------------------------\r
+  // This macro takes the raw pi0 spectrum from LHC11a_FitResult_20130314.root,\r
+  // correct it for feed down, then for efficiency,\r
+  // adds all systematic errors and produces the production invariant spectrum\r
+  //--\r
+  // Last modification: 06.07.2012 by Yuri Kharlov\r
+  //-----------------------------------------------------------------------------\r
+\r
+  TFile *f  = new TFile("LHC11a_FitResult_20130913.root");\r
+  TH1D * nr1CB    = (TH1D*)f->Get("Mix_CB_yr1")  ;\r
+  TH1D * nr1intCB = (TH1D*)f->Get("Mix_CB_yr1int")  ;\r
+  TH1D * nr2CB    = (TH1D*)f->Get("Mix_CB_yr2")  ;\r
+  TH1D * nr2intCB = (TH1D*)f->Get("Mix_CB_yr2int")  ;\r
+\r
+  TH1D * nr1GS    = (TH1D*)f->Get("Mix_yr1")  ;\r
+  TH1D * nr1intGS = (TH1D*)f->Get("Mix_yr1int")  ;\r
+  TH1D * nr2GS    = (TH1D*)f->Get("Mix_yr2")  ;\r
+  TH1D * nr2intGS = (TH1D*)f->Get("Mix_yr2int")  ;\r
+\r
+  //Divide by bin width\r
+  for(Int_t i=1;i<= nr1CB->GetNbinsX();i++){\r
+    nr1CB   ->SetBinContent(i,nr1CB->GetBinContent(i)/nr1CB->GetXaxis()->GetBinWidth(i)) ;\r
+    nr1CB   ->SetBinError  (i,nr1CB->GetBinError(i)/nr1CB->GetXaxis()->GetBinWidth(i)) ;\r
+    nr1intCB->SetBinContent(i,nr1intCB->GetBinContent(i)/nr1intCB->GetXaxis()->GetBinWidth(i)) ;\r
+    nr1intCB->SetBinError  (i,nr1intCB->GetBinError(i)/nr1intCB->GetXaxis()->GetBinWidth(i)) ;\r
+    nr2CB   ->SetBinContent(i,nr2CB->GetBinContent(i)/nr2CB->GetXaxis()->GetBinWidth(i)) ;\r
+    nr2CB   ->SetBinError  (i,nr2CB->GetBinError(i)/nr2CB->GetXaxis()->GetBinWidth(i)) ;\r
+    nr2intCB->SetBinContent(i,nr2intCB->GetBinContent(i)/nr2intCB->GetXaxis()->GetBinWidth(i)) ;\r
+    nr2intCB->SetBinError  (i,nr2intCB->GetBinError(i)/nr2intCB->GetXaxis()->GetBinWidth(i)) ;\r
+\r
+    nr1GS   ->SetBinContent(i,nr1GS->GetBinContent(i)/nr1GS->GetXaxis()->GetBinWidth(i)) ;\r
+    nr1GS   ->SetBinError  (i,nr1GS->GetBinError(i)/nr1GS->GetXaxis()->GetBinWidth(i)) ;\r
+    nr1intGS->SetBinContent(i,nr1intGS->GetBinContent(i)/nr1intGS->GetXaxis()->GetBinWidth(i)) ;\r
+    nr1intGS->SetBinError  (i,nr1intGS->GetBinError(i)/nr1intGS->GetXaxis()->GetBinWidth(i)) ;\r
+    nr2GS   ->SetBinContent(i,nr2GS->GetBinContent(i)/nr2GS->GetXaxis()->GetBinWidth(i)) ;\r
+    nr2GS   ->SetBinError  (i,nr2GS->GetBinError(i)/nr2GS->GetXaxis()->GetBinWidth(i)) ;\r
+    nr2intGS->SetBinContent(i,nr2intGS->GetBinContent(i)/nr2intGS->GetXaxis()->GetBinWidth(i)) ;\r
+    nr2intGS->SetBinError  (i,nr2intGS->GetBinError(i)/nr2intGS->GetXaxis()->GetBinWidth(i)) ;\r
+  }\r
+\r
+  // feed down correction\r
+\r
+  TF1 *fKaonContaminationToPi0 = new TF1("kaonCont","1./(1.-1.33*1.2*exp(-2.95-0.16*x))",0.,30.);\r
+  nr1GS   ->Divide(fKaonContaminationToPi0);\r
+  nr1intGS->Divide(fKaonContaminationToPi0);\r
+  nr2GS   ->Divide(fKaonContaminationToPi0);\r
+  nr2intGS->Divide(fKaonContaminationToPi0);\r
+  nr1CB   ->Divide(fKaonContaminationToPi0);\r
+  nr2CB   ->Divide(fKaonContaminationToPi0);\r
+  nr1intCB->Divide(fKaonContaminationToPi0);\r
+  nr2intCB->Divide(fKaonContaminationToPi0);\r
+\r
+  // SPD pileup correction\r
+\r
+  TF1 *fSPDpileup = new TF1("SPDpileup","0.988",0.,30.);\r
+  nr1GS   ->Multiply(fSPDpileup);\r
+  nr1intGS->Multiply(fSPDpileup);\r
+  nr2GS   ->Multiply(fSPDpileup);\r
+  nr2intGS->Multiply(fSPDpileup);\r
+  nr1CB   ->Multiply(fSPDpileup);\r
+  nr2CB   ->Multiply(fSPDpileup);\r
+  nr1intCB->Multiply(fSPDpileup);\r
+  nr2intCB->Multiply(fSPDpileup);\r
+\r
+  //correct for efficiency\r
+  TFile *fEff = new TFile("Pi0_efficiency_LHC11a__20131029_Mall.root") ;\r
+\r
+  TF1 * effGS=fEff->Get("eff_Pi0_Gaus_2760GeV") ;\r
+  TF1 * effCB=fEff->Get("eff_Pi0_CB_2760GeV") ;\r
+  effGS   ->SetRange(0.,25.) ;\r
+  effCB   ->SetRange(0.,25.) ;\r
+\r
+  nr1GS   ->Divide(effGS) ;\r
+  nr1intGS->Divide(effGS) ;\r
+  nr2GS   ->Divide(effGS) ;\r
+  nr2intGS->Divide(effGS) ;\r
+  nr1CB   ->Divide(effCB) ;\r
+  nr2CB   ->Divide(effCB) ;\r
+  nr1intCB->Divide(effCB) ;\r
+  nr2intCB->Divide(effCB) ;\r
+\r
+  //make 1/pt\r
+  for(Int_t i=1;i<=nr1CB->GetNbinsX();i++){\r
+    Double_t pt = TMath::TwoPi()*nr1CB->GetXaxis()->GetBinCenter(i);\r
+    nr1CB   ->SetBinContent(i,nr1CB   ->GetBinContent(i)/pt) ;\r
+    nr1CB   ->SetBinError  (i,nr1CB   ->GetBinError(i)  /pt) ;\r
+    nr1intCB->SetBinContent(i,nr1intCB->GetBinContent(i)/pt) ;\r
+    nr1intCB->SetBinError  (i,nr1intCB->GetBinError(i)  /pt) ;\r
+    nr2CB   ->SetBinContent(i,nr2CB   ->GetBinContent(i)/pt) ;\r
+    nr2CB   ->SetBinError  (i,nr2CB   ->GetBinError(i)  /pt) ;\r
+    nr2intCB->SetBinContent(i,nr2intCB->GetBinContent(i)/pt) ;\r
+    nr2intCB->SetBinError  (i,nr2intCB->GetBinError(i)  /pt) ;\r
+\r
+    nr1GS   ->SetBinContent(i,nr1GS   ->GetBinContent(i)/pt) ;\r
+    nr1GS   ->SetBinError  (i,nr1GS   ->GetBinError(i)  /pt) ;\r
+    nr1intGS->SetBinContent(i,nr1intGS->GetBinContent(i)/pt) ;\r
+    nr1intGS->SetBinError  (i,nr1intGS->GetBinError(i)  /pt) ;\r
+    nr2GS   ->SetBinContent(i,nr2GS   ->GetBinContent(i)/pt) ;\r
+    nr2GS   ->SetBinError  (i,nr2GS   ->GetBinError(i)  /pt) ;\r
+    nr2intGS->SetBinContent(i,nr2intGS->GetBinContent(i)/pt) ;\r
+    nr2intGS->SetBinError  (i,nr2intGS->GetBinError(i)  /pt) ;\r
+  }\r
+\r
+  //For the final spectrum we take average of fits\r
+  //with numerical integration of entries in signal\r
+  \r
+  TH1D * hStat = (TH1D*)nr2intCB->Clone("hPi02760GeVStat") ;\r
+  TH1D * hSys  = (TH1D*)hStat   ->Clone("hPi02760GeVSys") ;\r
+  TH1D * hSys2 = (TH1D*)hStat   ->Clone("hPi02760GeVSysTypeB") ;\r
+  TH1D * hSys3 = (TH1D*)hStat   ->Clone("hPi02760GeVSysTypeC") ;\r
+  hStat->SetAxisRange(0.,14.9,"X");\r
+  hSys ->SetAxisRange(0.,14.9,"X");\r
+\r
+  //For systematic error estimate take largest deviation\r
+  //of integrated yeilds (note, they are efficiency corrected)\r
+  for(Int_t i=1;i<=nr1CB->GetNbinsX();i++){\r
+    Double_t mean= hStat->GetBinContent(i) ;\r
+    Double_t dev = TMath::Max(\r
+                   TMath::Max(TMath::Abs(nr1intCB->GetBinContent(i)-mean),\r
+                              TMath::Abs(nr2intCB->GetBinContent(i)-mean)),\r
+                   TMath::Max(TMath::Abs(nr1intGS->GetBinContent(i)-mean),\r
+                              TMath::Abs(nr2intGS->GetBinContent(i)-mean))\r
+                             );\r
+    hSys ->SetBinError(i,dev) ;\r
+    hSys2->SetBinError(i,dev) ;\r
+  }\r
+\r
+  //Add other sys errors\r
+  TF1 * globalE = new TF1("globalE","1.-((x+1.354)/(x*1.002+1.354))^6.18 ",1.,30.) ; \r
+  TF1 * conv    = new TF1("conversion","0.035",0.,30.) ;\r
+  TF1 * accept  = new TF1("accept"    ,"0.01" ,0.,30.) ;\r
+  TF1 * pileup  = new TF1("pileup"    ,"0.004",0.,30.) ;\r
+  TF1 * calib   = new TF1("calib"     ,"0.005",0.,30.) ;\r
+  TF1 * modDiff = new TF1("modDiff"   ,"0.04",0.,30.) ;\r
+  // TF1 * modDiff = new TF1("modDiff"   ,"16.9*exp(-4.5*x)+0.033",0.,30.) ;\r
+  TF1 * tofCut  = new TF1("tofCut"    ,"0.0105" ,0.,30.) ;\r
+\r
+  //Borya's estimate of non-linearity (found for pp @ 7 TeV)\r
+  TF1 * nonlin= new TF1("nl","0.015+7.38*exp(-x/0.24)",0.,30.) ;\r
+\r
+  //Draw sys errors\r
+  TH1D * hRelSysRaw = (TH1D*)hSys->Clone("RelSysRaw") ;\r
+  hRelSysRaw->SetTitle("Summary of systematic errors");\r
+  for(Int_t i=1;i<=hSys->GetNbinsX();i++){\r
+    Double_t mean= hSys->GetBinContent(i) ;\r
+    Double_t a=hSys->GetBinError(i) ;\r
+    if(mean>0)\r
+      hRelSysRaw->SetBinContent(i,a/mean) ;\r
+    else\r
+      hRelSysRaw->SetBinContent(i,0.) ;\r
+      hRelSysRaw->SetBinError(i,0.) ;\r
+  }\r
+\r
+  //Add errors in sys errors\r
+\r
+  TH1D * hRelSysTot = (TH1D*)hSys->Clone("RelSys") ;\r
+  hRelSysTot->SetTitle("Summary of systematic uncertainties");\r
+\r
+  for(Int_t i=1;i<=hSys->GetNbinsX();i++){\r
+    Double_t pt   = hSys->GetXaxis()->GetBinCenter(i) ;\r
+    Double_t mean = hSys->GetBinContent(i) ;\r
+    Double_t a    = hSys->GetBinError(i) ;\r
+    // Double_t b    = mean * hSysErrModules->GetBinContent(i) ;\r
+    Double_t tot= mean*mean*nonlin ->Eval(pt)*nonlin ->Eval(pt) \r
+                 +mean*mean*conv   ->Eval(pt)*conv   ->Eval(pt)\r
+                 +mean*mean*accept ->Eval(pt)*accept ->Eval(pt)\r
+                 +mean*mean*pileup ->Eval(pt)*pileup ->Eval(pt)\r
+                 +mean*mean*calib  ->Eval(pt)*calib  ->Eval(pt)\r
+                 +mean*mean*modDiff->Eval(pt)*modDiff->Eval(pt)\r
+                 +mean*mean*tofCut ->Eval(pt)*tofCut ->Eval(pt)\r
+                 +mean*mean*globalE->Eval(pt)*globalE->Eval(pt); \r
+    Double_t raa= mean*mean*nonlin->Eval(pt)*nonlin->Eval(pt) \r
+                 +mean*mean*pileup->Eval(pt)*pileup->Eval(pt)\r
+                 +mean*mean*calib ->Eval(pt)*calib ->Eval(pt); \r
+    hSys3->SetBinError(i,TMath::Sqrt(tot)) ;\r
+    // hSys->SetBinError(i,TMath::Sqrt(tot + a*a + b*b)) ;\r
+    hSys2->SetBinError(i,TMath::Sqrt(raa + a*a)) ;\r
+    hSys ->SetBinError(i,TMath::Sqrt(tot + a*a)) ;\r
+    \r
+    a = hSys->GetBinError(i) ;\r
+    printf("i=%d, %g+-%g\n",i,mean,a);\r
+    if(mean>0)\r
+      hRelSysTot->SetBinContent(i,a/mean) ;\r
+    else {\r
+      hRelSysTot->SetBinContent(i,0.) ;\r
+      hRelSysTot->SetBinError(i,0.) ;\r
+    }\r
+  }\r
+\r
+  // TFile *fSysErrModules = TFile::Open("PHOS_sysErr_modules.root");\r
+  // TH1D* hSysErrModules = (TH1D*)fSysErrModules->Get("hSyserr");\r
+\r
+  gStyle->SetOptStat(0);\r
+  TCanvas * c = new TCanvas("c","SysErrors") ;\r
+  c->SetLogy();\r
+  c->cd() ;\r
+  gPad->SetRightMargin(0.02);\r
+  gPad->SetTopMargin(0.07);\r
+  hRelSysTot->SetAxisRange(0.8    ,11.9,"X");\r
+  hRelSysTot->SetAxisRange(0.0031,0.45,"Y");\r
+  hRelSysTot->GetYaxis()->SetMoreLogLabels();\r
+  hRelSysTot->GetYaxis()->SetNoExponent();\r
+  hRelSysTot->SetNdivisions(520,"X");\r
+  hRelSysTot->SetLineColor(1) ;\r
+  hRelSysTot->SetLineWidth(2) ;\r
+  hRelSysRaw->SetLineColor(2) ;\r
+  hRelSysRaw->SetLineWidth(2) ;\r
+  globalE->SetLineColor(kGreen+2) ;\r
+  nonlin ->SetLineColor(4) ;\r
+  conv   ->SetLineColor(6) ;\r
+  accept ->SetLineColor(kOrange) ;\r
+  pileup ->SetLineColor(42);\r
+  calib  ->SetLineColor(44);\r
+  calib  ->SetLineStyle(2);\r
+  modDiff->SetLineColor(52);\r
+  modDiff->SetLineStyle(2);\r
+  tofCut ->SetLineColor(53);\r
+  tofCut ->SetLineStyle(3);\r
+  // hSysErrModules->SetLineColor(kOrange+2);\r
+  // hSysErrModules->SetLineStyle(2);\r
+  hRelSysTot->SetXTitle("p_{T} (GeV/c)") ;\r
+  hRelSysTot->SetYTitle("Rel.syst.error") ;\r
+  hRelSysTot->GetYaxis()->SetTitleOffset(1.2) ;\r
+  hRelSysTot->Draw("hist") ;\r
+  hRelSysRaw->Draw("h same") ;\r
+  globalE->Draw("same") ; \r
+  nonlin ->Draw("same") ;\r
+  conv   ->Draw("same") ;\r
+  accept ->Draw("same") ;\r
+  pileup ->Draw("same");\r
+  calib  ->Draw("same");\r
+  modDiff->Draw("same");\r
+  tofCut ->Draw("same");\r
+  // hSysErrModules->Draw("same ][");\r
+  TLegend * l = new TLegend(0.57,0.62,0.85,0.925) ;\r
+  l->SetFillColor(kWhite);\r
+  l->SetBorderSize(0);\r
+  l->AddEntry(hRelSysTot,"Total uncertainty","l") ;\r
+  l->AddEntry(hRelSysRaw,"Raw extraction","l") ;\r
+  l->AddEntry(conv   ,"Conversion","l") ;\r
+  l->AddEntry(nonlin ,"Non-linearity","l") ;\r
+  l->AddEntry(accept ,"Acceptance","l");\r
+  l->AddEntry(pileup ,"Pileup","l");\r
+  l->AddEntry(calib  ,"Rel.calib.","l");\r
+  l->AddEntry(modDiff,"Per-module yield","l");\r
+  l->AddEntry(tofCut ,"Timing cut","l");\r
+  // l->AddEntry(hSysErrModules,"Intermodule spectra","l");\r
+  l->AddEntry(globalE,"Global E scale","l") ;\r
+  l->Draw() ;\r
+\r
+  c->Print("LHC11a_SysErrors.eps");\r
+\r
+  hStat->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, stat.err.");\r
+  hStat->SetXTitle("p_{T}, GeV/c");\r
+  hStat->SetYTitle("1/N_{ev} 1/(2#pi p_{T}) d^{2}N/dydp_{T} (GeV^{-2}c^{2})");\r
+  hSys ->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, total syst.err.");\r
+  hSys ->SetXTitle("p_{T}, GeV/c");\r
+  hSys ->SetYTitle("1/N_{ev} 1/(2#pi p_{T}) d^{2}N/dydp_{T} (GeV^{-2}c^{2})");\r
+  hSys2->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, syst.err. for R_{AA}");\r
+  hSys3->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, apparatus syst.err.");\r
+  hSys ->SetFillColor(kBlue-10) ;\r
+  hSys ->SetNdivisions(520,"X");\r
+\r
+  hSysFinal = (TH1F*)hSys->Clone("hSysFinal");\r
+  hSysFinal ->SetTitle("Production #pi^{0} yield, pp @ 2.76 TeV, 07.11.2013");\r
+  hSysFinal ->GetYaxis()->SetTitleOffset(1.2);\r
+  hStat->GetYaxis()->SetTitleOffset(1.2);\r
+  hSysFinal ->SetAxisRange(0.8,9.9,"X");\r
+\r
+  TCanvas *c2 = new TCanvas("c2","Production spectrum");\r
+  gPad->SetRightMargin(0.02);\r
+  gPad->SetTopMargin(0.07);\r
+  gPad->SetGridx();\r
+  gPad->SetGridy();\r
+  c2->SetLogy();\r
+  hSysFinal->Draw("E2") ;\r
+  hStat->SetMarkerStyle(20) ;\r
+  hStat->SetMarkerColor(4) ;\r
+  hStat->SetLineColor(4) ;\r
+  hStat->Draw("same") ;\r
+  \r
+  // //Apply bin width correction\r
+  // TH1D * hBWcorr = BinWidthCorrection(hStat) ;\r
+  // hSys->Divide(hBWcorr) ;\r
+  // hSys2->Divide(hBWcorr) ;\r
+  // hSys3->Divide(hBWcorr) ;\r
+  // hStat->Divide(hBWcorr) ;\r
+\r
+  c2->Print("LHC11a_pi0Spectrum.eps");\r
+\r
+  TF1 *tsallis = FitTsallis(hStat,hSys);\r
+  \r
+  TFile fout("PHOS_pi0_2760eV_noBinWidthCorr_20131112.root","recreate") ;\r
+  hSys   ->Write() ;\r
+  hSys2  ->Write() ;\r
+  hSys3  ->Write() ;\r
+  hStat  ->Write() ;\r
+  tsallis->Write() ;\r
+  fout.Close() ;\r
+\r
+\r
+}\r
+\r
+//-----------------------------------------------------------------------------\r
+TF1 *FitTsallis(TH1D* hStat, TH1D* hSys)\r
+{\r
+  TF1 * fit = new TF1("Tsalis","[0]/2./3.1415*([2]-1.)*([2]-2.)/([2]*[1]*([2]*[1]+0.135*([2]-2.)))*(1.+(sqrt(x*x+0.135*0.135)-0.135)/([2]*[1]))^-[2]",0.5,25.) ;\r
+  fit->SetParameters(10.,0.2,8.) ;\r
+  fit->SetLineColor(kBlack) ;\r
+  fit->SetLineWidth(2) ;\r
+  fit->SetParName(0,"dN/dy") ;\r
+  fit->SetParName(1,"T") ;\r
+  fit->SetParName(2,"n") ;\r
+  \r
+  TH1D * hSysRatio  = (TH1D*)hSys ->Clone("RatioSys") ;\r
+  TH1D * hStatRatio = (TH1D*)hStat->Clone("RatioStat") ;\r
+\r
+  TH1D * hsum = (TH1D*)hStat->Clone("sum") ;\r
+  for(Int_t i=1; i<=hsum->GetNbinsX();i++){\r
+    Double_t a=hSys->GetBinError(i) ;\r
+    Double_t b=hStat->GetBinError(i) ;\r
+    hsum->SetBinError(i,TMath::Sqrt(a*a+b*b)) ;\r
+  }\r
+  hsum->SetStats(1) ;\r
+  hsum->Fit(fit,"Q") ;\r
+  Double_t meanPt=fit->Moment(1,1.,10.) ;\r
+  printf("<pt>=%f \n",meanPt) ;\r
+\r
+  hSysRatio ->Divide(fit) ;\r
+  hStatRatio->Divide(fit) ;\r
+\r
+  return fit;\r
+}\r
+//-----------------------------------------------------------------------------\r
+TH1D * BinWidthCorrection(TH1D * h){\r
+  //We apply bin width a-la PHENIX \r
+  //Use Tsalis fit to calculate shift in y direction\r
\r
+  TF1 * fit = new TF1("hag","[0]*(([2]*[1]+sqrt(x*x+0.135*0.135))/([2]*[1]+0.135))^-[2]",0.5,25.) ;\r
+  fit->SetParameters(10.,0.2,8.) ;\r
+  TCanvas * corr = new TCanvas("BWcorr","Bin width correction") ;\r
+  Int_t col[6]={kRed,kOrange,kMagenta,kGreen,kCyan,kBlue} ;\r
+  TH1D * hcorr[20] ;\r
+  char key[55] ;\r
+  Double_t rMax=10 ;\r
+  Int_t iteration=0 ;\r
+  TH1D * htmp = (TH1D*)h->Clone("tmp") ;\r
+  while(iteration<6){\r
+    printf(" Iteration %d: rMax=%f \n",iteration, rMax) ;\r
+    htmp->Fit(fit,"N") ;\r
+    sprintf(key,"Ineration%d",iteration) ;\r
+    hcorr[iteration]=(TH1D*)h->Clone(key);\r
+    rMax= 0; \r
+    for(Int_t i=1;i<=h->GetNbinsX();i++){\r
+      Double_t a=h->GetXaxis()->GetBinLowEdge(i) ;\r
+      Double_t b=h->GetXaxis()->GetBinUpEdge(i) ;\r
+      Double_t r=fit->Integral(a,b)/(b-a)/fit->Eval(0.5*(a+b)) ;\r
+      hcorr[iteration]->SetBinContent(i,r) ;\r
+      hcorr[iteration]->SetBinError(i,0.) ;\r
+      if(rMax<r)rMax=r ;\r
+    }\r
+    delete htmp ;\r
+    htmp = (TH1D*)h->Clone("tmp") ;\r
+    htmp->Divide(hcorr[iteration]) ;\r
+    corr->cd() ;\r
+    hcorr[iteration]->SetLineColor(col[iteration]);\r
+    if(iteration==0)\r
+      hcorr[iteration]->Draw() ;\r
+    else\r
+      hcorr[iteration]->Draw("same") ;\r
+    corr->Update() ;\r
+    iteration++ ;\r
+  } \r
+\r
+  hcorr[5]->SetTitle("Bin-width correction for #pi^{0} spectrum");\r
+  hcorr[5]->SetYTitle("Bin-width corrected / uncorrected");\r
+  TFile fout("PHOS_pi0_7TeV_BinWidthCorrection.root","recreate") ;\r
+  hcorr[5]->Write();\r
+  fout.Close();\r
+\r
+  return hcorr[5] ;\r
+}\r
diff --git a/PWGGA/PHOSTasks/PHOS_pp_pi0/macros/makeMmixPU_CB.C b/PWGGA/PHOSTasks/PHOS_pp_pi0/macros/makeMmixPU_CB.C
new file mode 100644 (file)
index 0000000..69a165b
--- /dev/null
@@ -0,0 +1,484 @@
+void makeMmixPU_CB(const TString histoFile="LHC11a_pass4_20130913.root",\r
+                  const Int_t nSigma=2,\r
+                   const char* module="A10")\r
+{\r
+  //Fit Real/Mixed ratio, normalize Mixed and subtract it from Real.\r
+  // The pi0 peak if fitted by the Crystal Ball function, \r
+  // the background is fitted by pol1 or pol2\r
+\r
+  TString hMassName;\r
+  TFile * file = new TFile(histoFile) ;\r
+  THashList *hList = (THashList*)file->Get("histESD");\r
+  char key[125] ;\r
+  hMassName = "hMassPt";\r
+  hMassName += module;\r
+  TH2F * h   = (TH2F*)hList->FindObject(hMassName) ;\r
+\r
+  hMassName = "hMiMassPt";\r
+  hMassName += module;\r
+  TH2F * hm  = (TH2F*)hList->FindObject(hMassName) ;\r
+\r
+  TH1F * hev = (TH1F*)hList->FindObject("hSelEvents") ;\r
+\r
+  // Array of pt bins\r
+  Int_t nPadX = 6, nPadY = 3;\r
+  Int_t nbin=18 ;\r
+  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.} ;\r
+  PPRstyle();\r
+  gStyle->SetPadLeftMargin(0.16);\r
+  gStyle->SetPadRightMargin(0.01);\r
+  gStyle->SetPadTopMargin(0.02);\r
+  gStyle->SetPadBottomMargin(0.11);\r
+  gStyle->SetTitleX(0.80);\r
+  gStyle->SetTitleY(0.98);\r
+  \r
+  //Fit real only \r
+  //Linear Bg\r
+  char key2[155];\r
+  sprintf(key,"Mix_CB") ;\r
+  sprintf(key2,"%s_mr1",key) ;\r
+  TH1D * mr1 = new TH1D(key2,"Mass",nbin,xa) ;\r
+  sprintf(key2,"%s_sr1",key) ;\r
+  TH1D * sr1 = new TH1D(key2,"Width",nbin,xa) ;\r
+  sprintf(key2,"%s_rms1",key) ;\r
+  TH1D * rms1 = new TH1D(key2,"Width",nbin,xa) ;\r
+  sprintf(key2,"%s_nnr1",key) ;\r
+  TH1D * nnr1 = new TH1D(key2,"Mass",nbin,xa) ;\r
+  sprintf(key2,"%s_ar1",key) ;\r
+  TH1D * ar1 = new TH1D(key2,"Width",nbin,xa) ;\r
+  sprintf(key2,"%s_yr1",key) ;\r
+  TH1D * nr1 = new TH1D(key2,"Raw yield",nbin,xa) ;\r
+  sprintf(key2,"%s_yr1int",key) ;\r
+  TH1D * nr1int = new TH1D(key2,"Raw yield, integrated",nbin,xa) ;\r
+\r
+  //Quadratic Bg\r
+  sprintf(key2,"%s_mr2",key) ;\r
+  TH1D * mr2 = new TH1D(key2,"Mass",nbin,xa) ;\r
+  sprintf(key2,"%s_sr2",key) ;\r
+  TH1D * sr2 = new TH1D(key2,"Width",nbin,xa) ;\r
+  sprintf(key2,"%s_rms2",key) ;\r
+  TH1D * rms2 = new TH1D(key2,"Width",nbin,xa) ;\r
+  sprintf(key2,"%s_nnr2",key) ;\r
+  TH1D * nnr2 = new TH1D(key2,"Mass",nbin,xa) ;\r
+  sprintf(key2,"%s_ar2",key) ;\r
+  TH1D * ar2 = new TH1D(key2,"Width",nbin,xa) ;\r
+  sprintf(key2,"%s_yr2",key) ;\r
+  TH1D * nr2 = new TH1D(key2,"Raw yield",nbin,xa) ;\r
+  sprintf(key2,"%s_yr2int",key) ;\r
+  TH1D * nr2int = new TH1D(key2,"Raw yield, integrated",nbin,xa) ;\r
+\r
+\r
+  TF1 * fit1 = new TF1("fit",CB,0.,1.,7) ;\r
+  fit->SetParName(0,"A") ;\r
+  fit->SetParName(1,"m_{0}") ;\r
+  fit->SetParName(2,"#sigma") ;\r
+  fit->SetParName(3,"a_{0}") ;\r
+  fit->SetParName(4,"a_{1}") ;\r
+  fit->SetLineWidth(2) ;\r
+  fit->SetLineColor(2) ;\r
+\r
+  TF1 * fgs = new TF1("gs",CBs,0.,1.,5) ;\r
+  fgs->SetLineColor(2) ;\r
+  fgs->SetLineWidth(1) ;\r
+\r
+  TF1 * fit2 = new TF1("fit2",CB2,0.,1.,8) ;\r
+  fit2->SetParName(0,"A") ;\r
+  fit2->SetParName(1,"m_{0}") ;\r
+  fit2->SetParName(2,"#sigma") ;\r
+  fit2->SetParName(3,"a_{0}") ;\r
+  fit2->SetParName(4,"a_{1}") ;\r
+  fit2->SetLineWidth(2) ;\r
+  fit2->SetLineColor(4) ;\r
+  fit2->SetLineStyle(2) ;\r
+\r
+  TF1 * fbg1 = new TF1("bg1",BG1,0.,1.,2) ;\r
+  TF1 * fbg2 = new TF1("bg2",BG2,0.,1.,3) ;\r
+\r
+  TCanvas * c3 = new TCanvas("mggFit1_CB_Signal","mggFitCB",10,10,1400,800) ;\r
+  c3->Divide(nPadX,nPadY) ;\r
+\r
+  TCanvas * cReal = new TCanvas("cReal","Mgg real events",10,10,1400,800) ;\r
+  cReal->Divide(nPadX,nPadY) ;\r
+\r
+  TCanvas * c1 = new TCanvas("mggFit1_CB","mggFit1",10,10,1400,800) ;\r
+  c1->Divide(nPadX,nPadY) ;\r
+  c1->cd(0) ; \r
+  TCanvas * c2=0,*c4=0,*c5=0,*c6=0 ; \r
+\r
+  TAxis * pta=h->GetYaxis() ;\r
+  TAxis * ma=h->GetXaxis() ;\r
+  for(Int_t i=1;i<=nbin;i++){\r
+    c1->cd(i) ;\r
+    Int_t imin=pta->FindBin(xa[i-1]+0.0001);\r
+    Int_t imax=pta->FindBin(xa[i]-0.0001) ;\r
+    Double_t pt=(xa[i]+xa[i-1])/2. ;\r
+    TH1D * hp = h->ProjectionX("re",imin,imax) ;\r
+    hp->Sumw2() ;\r
+    hp->SetNdivisions(505,"X");\r
+    TH1D * hpm= hm->ProjectionX("mi",imin,imax) ;\r
+    hpm->Sumw2() ;\r
+    hpm->SetNdivisions(505,"X");\r
+    if(i<15){\r
+      hp ->Rebin(2) ;\r
+      hpm->Rebin(2) ;\r
+    }\r
+    else{\r
+      hp ->Rebin(5) ;\r
+      hpm->Rebin(5) ;\r
+    }\r
+    hp ->SetNdivisions(506);\r
+    hpm->SetNdivisions(506);\r
+    hp ->SetLabelSize(0.05,"X");\r
+    hpm->SetLabelSize(0.05,"X");\r
+    hp ->SetTitleSize(0.00,"X");\r
+    hpm->SetTitleSize(0.00,"X");\r
+    hp ->SetLabelSize(0.05,"Y");\r
+    hpm->SetLabelSize(0.05,"Y");\r
+    hp ->SetTitleSize(0.05,"Y");\r
+    hpm->SetTitleSize(0.05,"Y");\r
+    hp ->SetTitleOffset(1.0,"X");\r
+    hpm->SetTitleOffset(1.0,"X");\r
+    if (i>12) {\r
+      hp ->SetLabelSize(0.05,"X");\r
+      hpm->SetLabelSize(0.05,"X");\r
+      hp ->SetTitleSize(0.05,"X");\r
+      hpm->SetTitleSize(0.05,"X");\r
+      hp ->SetXTitle("m_{#gamma#gamma} (GeV/c^{2})");\r
+      hpm->SetXTitle("m_{#gamma#gamma} (GeV/c^{2})");\r
+    }\r
+    // //Asign errors to the zero bins\r
+    // for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hp->GetBinContent(ib)==0)hp->SetBinError(ib,1.);}\r
+    // for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hpm->GetBinContent(ib)==0)hpm->SetBinError(ib,1.);}\r
+\r
+    TH1D * hpm2   = (TH1D*)hpm->Clone("Bg1") ;\r
+    TH1D * hpcopy = (TH1D*)hp ->Clone("hpcopy") ;\r
+    TH1D * hp2    = (TH1D*)hp ->Clone("hp2") ;\r
+    TH1D * hpReal = (TH1D*)hp ->Clone("hpReal") ;\r
+    hpcopy->Divide(hpm) ;\r
+    sprintf(key,"%3.1f<p_{T}<%3.1f GeV/c",xa[i-1],xa[i]) ;\r
+    hpcopy->SetTitle(key) ;\r
+    hpReal->SetTitle(key) ;\r
+    hpReal->SetLineWidth(2);\r
+    hpcopy->SetMarkerStyle(20) ;\r
+    hpcopy->SetMarkerSize(0.7) ;\r
+\r
+    Double_t mInit = 0.136., wInit = 0.005;\r
+    if(i==1) {\r
+      mInit = 0.139;\r
+      wInit = 0.007;\r
+    }\r
+    if(i==2) {\r
+      mInit = 0.136;\r
+      wInit = 0.007;\r
+    }\r
+    fit1->SetParameters(0.001,mInit,wInit,9.,0.5,0.0013,0.) ;\r
+    fit1->SetParLimits(2,0.003,0.010) ;\r
+    fit1->SetParLimits(1,0.132,0.139) ;\r
+    fit1->SetParLimits(0,0.,1e+6) ;\r
+//    fit1->SetParLimits(3,0.,10.) ;\r
+//    fit1->SetParLimits(4,0.,1e+2) ;\r
+    fit1->FixParameter(3,1.60) ;\r
+    fit1->FixParameter(4,1.27) ;\r
+\r
+    //Select fitting range\r
+    Double_t rangeMin=0.06 ;\r
+    Double_t rangeMax=0.22 ;\r
+    if(i>=16)rangeMax=0.3 ; //More points to fix background\r
+    hpcopy->Fit(fit1,"NQ","",rangeMin,rangeMax) ;\r
+    hpcopy->Fit(fit1,"MQ","",rangeMin,rangeMax) ;\r
+\r
+    fit2->SetParameters(fit1->GetParameters()) ;\r
+    fit2->SetParameter(2,wInit);\r
+    fit2->SetParLimits(2,0.003,0.008) ;\r
+    fit2->SetParLimits(1,0.130,0.142) ;\r
+    fit2->SetParLimits(0,0.,1e+6) ;\r
+//    fit2->SetParLimits(3,0.,10.) ;\r
+//    fit2->SetParLimits(4,0.,1e+2) ;\r
+    fit2->FixParameter(3,1.60) ;\r
+    fit2->FixParameter(4,1.27) ;\r
+    fit2->SetParameter(7,0.) ;\r
+\r
+    hpcopy->Fit(fit2,"+NQ","",rangeMin,rangeMax) ;\r
+    hpcopy->Fit(fit2,"+MQ","",rangeMin,rangeMax) ;\r
+\r
+    c1->cd(i) ;\r
+    hpcopy->SetAxisRange(0.065,0.229,"X") ;\r
+    hpcopy->Draw() ;\r
+\r
+    cReal->cd(i) ;\r
+    hpReal->SetAxisRange(0.065,0.229,"X") ;\r
+    hpReal->Draw("ehist") ;\r
+\r
+    fbg1->SetParameters(fit1->GetParameter(5),fit1->GetParameter(6)); \r
+    fbg2->SetParameters(fit2->GetParameter(5),fit2->GetParameter(6),fit2->GetParameter(7)); \r
+    Double_t intRangeMin = PeakPosition(pt)-nSigma*PeakWidth(pt) ;\r
+    Double_t intRangeMax = PeakPosition(pt)+nSigma*PeakWidth(pt) ;\r
+    Int_t    intBinMin   = hp->GetXaxis()->FindBin(intRangeMin) ;\r
+    Int_t    intBinMax   = hp->GetXaxis()->FindBin(intRangeMax) ;\r
+    Double_t errStat     = hpm->Integral(intBinMin,intBinMax); \r
+\r
+    hpm ->Multiply(fbg1) ;\r
+    hpm2->Multiply(fbg2) ;\r
+    hp  ->Add(hpm,-1.) ;\r
+    hp2 ->Add(hpm2,-1.) ;\r
+\r
+    c3->cd(i) ;\r
+\r
+    if(i<15)\r
+      fgs->SetParameters(hp->Integral(32,36)/5.,fit1->GetParameter(1),fit1->GetParameter(2),fit1->GetParameter(3),fit1->GetParameter(4)) ;\r
+    else\r
+      fgs->SetParameters(hp->Integral(13,15)/3.,fit1->GetParameter(1),fit1->GetParameter(2),fit1->GetParameter(3),fit1->GetParameter(4)) ;\r
+//      fgs->SetParameters(hp->Integral(13,15)/3.,0.135,0.008,0.7,2.) ;\r
+\r
+    fgs->SetParLimits(0,0.,1e+06) ;\r
+    fgs->SetParLimits(1,0.130,0.142) ;\r
+    fgs->SetParLimits(2,0.003,0.008) ;\r
+    fgs->FixParameter(3,1.60) ;\r
+    fgs->FixParameter(4,1.27) ;\r
+    \r
+    hp->Fit(fgs,"Q","",rangeMin,rangeMax) ;   \r
+    hp->SetMaximum(hp2->GetMaximum()*1.1) ;\r
+    hp->SetMinimum(hp2->GetMinimum()*1.1) ;\r
+    hp->SetMarkerStyle(20) ;\r
+    hp->SetMarkerSize(0.7) ;\r
+    mr1->SetBinContent (i,fgs->GetParameter(1)) ;\r
+    mr1->SetBinError   (i,fgs->GetParError(1)) ;\r
+    sr1->SetBinContent (i,TMath::Abs(fgs->GetParameter(2))) ;\r
+    sr1->SetBinError   (i,fgs->GetParError(2)) ;\r
+    nnr1->SetBinContent(i,fgs->GetParameter(3)) ;\r
+    nnr1->SetBinError  (i,fgs->GetParError(3)) ;\r
+    ar1->SetBinContent (i,fgs->GetParameter(4)) ;\r
+    ar1->SetBinError   (i,fgs->GetParError(4)) ;\r
+\r
+    Double_t y=fgs->Integral(intRangeMin,intRangeMax)/hp->GetXaxis()->GetBinWidth(1) ;\r
+    nr1->SetBinContent(i,y) ;\r
+    if(y>0)\r
+      rms1->SetBinContent(i,fgs->CentralMoment(2.,0.05,0.2)) ;\r
+    Double_t ey=0 ;\r
+    if(fgs->GetParameter(0)!=0. && fgs->GetParameter(2)!=0.){\r
+      Double_t en=fgs->GetParError(0)/fgs->GetParameter(0) ;\r
+      Double_t es=fgs->GetParError(2)/fgs->GetParameter(2) ;\r
+      ey=y*TMath::Sqrt(en*en+es*es) ;\r
+    }\r
+    nr1->SetBinError(i,ey) ;\r
+\r
+    Double_t npiInt=hp->Integral(intBinMin,intBinMax) ;\r
+    Double_t norm=fbg1->GetParameter(0) ;\r
+    Double_t normErr=fbg1->GetParError(0) ;\r
+    if(npiInt>0.){\r
+      nr1int->SetBinContent(i,npiInt) ;\r
+      nr1int->SetBinError(i,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;\r
+    }\r
+\r
+    hp2->GetXaxis()->SetRangeUser(0.05,0.25) ;\r
+    hp2->SetMaximum(hp2->GetMaximum()*1.1) ;\r
+    hp2->SetMinimum(hp2->GetMinimum()*1.1) ;\r
+    hp2->SetMarkerStyle(20) ;\r
+    hp2->SetMarkerSize(0.7) ;\r
+\r
+    hp2->Fit(fgs,"Q","",rangeMin,rangeMax) ;\r
+    mr2->SetBinContent (i,fgs->GetParameter(1)) ;\r
+    mr2->SetBinError   (i,fgs->GetParError(1)) ;\r
+    sr2->SetBinContent (i,TMath::Abs(fgs->GetParameter(2))) ;\r
+    sr2->SetBinError   (i,fgs->GetParError(2)) ;\r
+    nnr2->SetBinContent(i,fgs->GetParameter(3)) ;\r
+    nnr2->SetBinError  (i,fgs->GetParError(3)) ;\r
+    ar2->SetBinContent (i,fgs->GetParameter(4)) ;\r
+    ar2->SetBinError   (i,fgs->GetParError(4)) ;\r
+\r
+    y=fgs->Integral(intRangeMin,intRangeMax)/hp->GetXaxis()->GetBinWidth(1) ;\r
+    nr2->SetBinContent(i,y) ;\r
+    ey=0 ;\r
+    if(fgs->GetParameter(0)!=0. && fgs->GetParameter(2)!=0.){\r
+      Double_t en=fgs->GetParError(0)/fgs->GetParameter(0) ;\r
+      Double_t es=fgs->GetParError(2)/fgs->GetParameter(2) ;\r
+      ey=y*TMath::Sqrt(en*en+es*es) ;\r
+    }\r
+    nr2->SetBinError(i,ey) ;\r
+    npiInt=hp2->Integral(intBinMin,intBinMax) ;\r
+    norm=fbg2->GetParameter(0) ;\r
+    normErr=fbg2->GetParError(0) ;\r
+    if(npiInt>0.){\r
+      nr2int->SetBinContent(i,npiInt) ;\r
+      nr2int->SetBinError(i,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;\r
+    } \r
+\r
+    hp2->SetTitle(key) ;\r
+    hp2->SetAxisRange(0.065,0.229,"X") ;\r
+    hp2->Draw() ;\r
+    delete hp ;\r
+//    delete hp2 ;\r
+//    delete hpcopy ;\r
+    delete hpm ;\r
+    delete hpm2 ;\r
+    c1->Update() ;\r
+    c3->Update() ;\r
+    cReal->Update() ;\r
+  }\r
+\r
+  if (c3) c3->Print("Pi0_Signal_CB.eps") ;\r
+  if (c1) c1->Print("Pi0_Ratio_CB.eps") ;\r
+  if (cReal) cReal->Print("Pi0_InvMass_CB.eps") ;\r
+\r
+  //Normalize by the number of non-pileup events\r
+  Double_t nMBOR   = hev->GetBinContent(2); // MBOR events pass4\r
+  // Double_t nPU     = hev->GetBinContent(7);\r
+  // Double_t nMBOR   = hev->GetBinContent(1); // MBOR events pass3\r
+  // Double_t nPU     = hev->GetBinContent(6);\r
+  Double_t nevents = nMBOR;\r
+  // Double_t nevents = nMBOR - nPU;\r
+  printf("==============\nN events = %d\n==============\n",nevents);\r
+\r
+  nr1   ->Scale(1./nevents) ;\r
+  nr1int->Scale(1./nevents) ;\r
+  nr2   ->Scale(1./nevents) ;\r
+  nr2int->Scale(1./nevents) ;\r
+\r
+\r
+  TFile fout("LHC11a_FitResult.root","update");\r
+  mr1   ->Write(0,TObject::kOverwrite) ;\r
+  sr1   ->Write(0,TObject::kOverwrite) ;\r
+  rms1  ->Write(0,TObject::kOverwrite) ;\r
+  nnr1  ->Write(0,TObject::kOverwrite) ;\r
+  ar1   ->Write(0,TObject::kOverwrite) ;\r
+  nr1   ->Write(0,TObject::kOverwrite) ;\r
+  nr1int->Write(0,TObject::kOverwrite) ;\r
+  mr2   ->Write(0,TObject::kOverwrite) ;\r
+  sr2   ->Write(0,TObject::kOverwrite) ;\r
+  rms2  ->Write(0,TObject::kOverwrite) ;\r
+  nnr2  ->Write(0,TObject::kOverwrite) ;\r
+  ar2   ->Write(0,TObject::kOverwrite) ;\r
+  nr2   ->Write(0,TObject::kOverwrite) ;\r
+  nr2int->Write(0,TObject::kOverwrite) ;\r
+  fout.Close() ;\r
+\r
+}\r
+\r
+\r
+//-----------------------------------------------------------------------------\r
+const Double_t kMean=0.135 ; //Approximate peak position to facilitate error estimate\r
+\r
+Double_t PeakPosition(Double_t pt){\r
+  //Fit to the measured peak position\r
+  // return 4.99292e-003*exp(-pt*9.32300e-001)+1.34944e-001 ;\r
+  return 0.57*exp(-pt*7.62)+0.13592 ;\r
+}\r
+Double_t PeakWidth(Double_t pt){\r
+  //fit to the measured peak width\r
+  return 1.60935e-02*exp(-pt*2.25609e+00)+4.65743e-03 ;\r
+  // return 0.0068 ;\r
+}\r
\r
+Double_t CB(Double_t * x, Double_t * par){\r
+  Double_t m=par[1] ;\r
+  Double_t s=par[2] ;\r
+  Double_t n=par[3] ;\r
+  Double_t a=par[4] ;\r
+  Double_t dx=(x[0]-m)/s ;\r
+  if(dx>-a)\r
+    return par[0]*exp(-dx*dx/2.)+par[5]+par[6]*(x[0]-kMean) ;\r
+  else{\r
+    Double_t A=TMath::Power((n/TMath::Abs(a)),n)*TMath::Exp(-a*a/2) ;\r
+    Double_t B=n/TMath::Abs(a)-TMath::Abs(a) ;\r
+    return par[0]*A*TMath::Power((B-dx),-n)+par[5]+par[6]*(x[0]-kMean) ;\r
+  }\r
+\r
+}\r
+Double_t CB2(Double_t * x, Double_t * par){\r
+  Double_t m=par[1] ;\r
+  Double_t s=par[2] ;\r
+  Double_t n=par[3] ;\r
+  Double_t a=par[4] ;\r
+  Double_t dx=(x[0]-m)/s ;\r
+  if(dx>-a)\r
+    return par[0]*exp(-dx*dx/2.)+par[5]+par[6]*(x[0]-kMean)+par[7]*(x[0]-kMean)*(x[0]-kMean) ;\r
+  else{\r
+    Double_t A=TMath::Power((n/TMath::Abs(a)),n)*TMath::Exp(-a*a/2) ;\r
+    Double_t B=n/TMath::Abs(a)-TMath::Abs(a) ;\r
+    return par[0]*A*TMath::Power((B-dx),-n)+par[5]+par[6]*(x[0]-kMean)+par[7]*(x[0]-kMean)*(x[0]-kMean) ;\r
+  }\r
+\r
+}\r
+Double_t CBs(Double_t * x, Double_t * par){\r
+  Double_t m=par[1] ;\r
+  Double_t s=par[2] ;\r
+  Double_t n=par[3] ;\r
+  Double_t a=par[4] ;\r
+  Double_t dx=(x[0]-m)/s ;\r
+  if(dx>-a)\r
+    return par[0]*exp(-dx*dx/2.) ;\r
+  else{\r
+    Double_t A=TMath::Power((n/TMath::Abs(a)),n)*TMath::Exp(-a*a/2) ;\r
+    Double_t B=n/TMath::Abs(a)-TMath::Abs(a) ;\r
+    return par[0]*A*TMath::Power((B-dx),-n) ;\r
+  }\r
+}\r
+Double_t BG1(Double_t * x, Double_t * par){\r
+  //Normalizatino of Mixed\r
+  return par[0]+par[1]*(x[0]-kMean) ;\r
+}\r
+Double_t BG2(Double_t * x, Double_t * par){\r
+  //Another normalization of  Mixed\r
+  return par[0]+par[1]*(x[0]-kMean)+par[2]*(x[0]-kMean)*(x[0]-kMean) ;\r
+}\r
+\r
+//-----------------------------------------------------------------------------\r
+PPRstyle()\r
+{\r
+\r
+  //////////////////////////////////////////////////////////////////////\r
+  //\r
+  // ROOT style macro for the TRD TDR\r
+  //\r
+  //////////////////////////////////////////////////////////////////////\r
+\r
+  gStyle->SetPalette(1);\r
+  gStyle->SetCanvasBorderMode(-1);\r
+  gStyle->SetCanvasBorderSize(1);\r
+  gStyle->SetCanvasColor(10);\r
+\r
+  gStyle->SetFrameFillColor(10);\r
+  gStyle->SetFrameBorderSize(1);\r
+  gStyle->SetFrameBorderMode(-1);\r
+  gStyle->SetFrameLineWidth(1.2);\r
+  gStyle->SetFrameLineColor(1);\r
+\r
+  gStyle->SetHistFillColor(0);\r
+  gStyle->SetHistLineWidth(1);\r
+  gStyle->SetHistLineColor(1);\r
+\r
+  gStyle->SetPadColor(10);\r
+  gStyle->SetPadBorderSize(1);\r
+  gStyle->SetPadBorderMode(-1);\r
+\r
+  gStyle->SetStatColor(10);\r
+  gStyle->SetTitleColor(kBlack,"X");\r
+  gStyle->SetTitleColor(kBlack,"Y");\r
+\r
+  gStyle->SetLabelSize(0.04,"X");\r
+  gStyle->SetLabelSize(0.04,"Y");\r
+  gStyle->SetLabelSize(0.04,"Z");\r
+  gStyle->SetTitleSize(0.04,"X");\r
+  gStyle->SetTitleSize(0.04,"Y");\r
+  gStyle->SetTitleSize(0.04,"Z");\r
+  gStyle->SetTitleFont(42,"X");\r
+  gStyle->SetTitleFont(42,"Y");\r
+  gStyle->SetTitleFont(42,"X");\r
+  gStyle->SetLabelFont(42,"X");\r
+  gStyle->SetLabelFont(42,"Y");\r
+  gStyle->SetLabelFont(42,"Z");\r
+  gStyle->SetStatFont(42);\r
+\r
+  gStyle->SetTitleOffset(1.0,"X");\r
+  gStyle->SetTitleOffset(1.4,"Y");\r
+\r
+  gStyle->SetFillColor(kWhite);\r
+  gStyle->SetTitleFillColor(kWhite);\r
+\r
+  gStyle->SetOptDate(0);\r
+  gStyle->SetOptTitle(1);\r
+  gStyle->SetOptStat(0);\r
+  gStyle->SetOptFit(0);\r
+\r
+}\r
+\r
diff --git a/PWGGA/PHOSTasks/PHOS_pp_pi0/macros/makeMmixPU_GS.C b/PWGGA/PHOSTasks/PHOS_pp_pi0/macros/makeMmixPU_GS.C
new file mode 100644 (file)
index 0000000..ea6c4b6
--- /dev/null
@@ -0,0 +1,421 @@
+//-----------------------------------------------------------------------------\r
+void makeMmixPU_GS(const TString histoFile="LHC11a_pass4_20130913.root",\r
+                  const Int_t nSigma=2,\r
+                  const char* module="A10")\r
+{\r
+  //Fit Real/Mixed ratio, normalize Mixed and subtract it from Real.\r
+  // The pi0 peak if fitted by the Gaussian function, \r
+  // the background is fitted by pol1 or pol2\r
+\r
+  TString hMassName;\r
+  TFile * file = new TFile(histoFile) ;\r
+  THashList *hList = (THashList*)file->Get("histESD");\r
+  char key[125] ;\r
+  hMassName = "hMassPt";\r
+  hMassName += module;\r
+  TH2F * h   = (TH2F*)hList->FindObject(hMassName) ;\r
+\r
+  hMassName = "hMiMassPt";\r
+  hMassName += module;\r
+  TH2F * hm  = (TH2F*)hList->FindObject(hMassName) ;\r
+\r
+  TH1F * hev = (TH1F*)hList->FindObject("hSelEvents") ;\r
+\r
+  // Array of pt bins\r
+  Int_t nbin=18 ;\r
+  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.} ;\r
+  PPRstyle();\r
+  gStyle->SetPadLeftMargin(0.16);\r
+  gStyle->SetPadRightMargin(0.01);\r
+  gStyle->SetPadTopMargin(0.02);\r
+  gStyle->SetPadBottomMargin(0.11);\r
+  gStyle->SetTitleX(0.80);\r
+  gStyle->SetTitleY(0.98);\r
+\r
+  //Fit real only \r
+  //Linear Bg\r
+  char key2[155];\r
+  sprintf(key,"Mix") ;\r
+  sprintf(key2,"%s_mr1",key) ;\r
+  TH1D * mr1 = new TH1D(key2,"Mass",nbin,xa) ;\r
+  sprintf(key2,"%s_sr1",key) ;\r
+  TH1D * sr1 = new TH1D(key2,"Width",nbin,xa) ;\r
+  sprintf(key2,"%s_ar1",key) ;\r
+  TH1D * ar1 = new TH1D(key2,"a",nbin,xa) ;\r
+  sprintf(key2,"%s_br1",key) ;\r
+  TH1D * br1 = new TH1D(key2,"a",nbin,xa) ;\r
+  sprintf(key2,"%s_yr1",key) ;\r
+  TH1D * nr1 = new TH1D(key2,"Raw yield",nbin,xa) ;\r
+  sprintf(key2,"%s_yr1int",key) ;\r
+  TH1D * nr1int = new TH1D(key2,"Raw yield, integrated",nbin,xa) ;\r
+\r
+  //Quadratic Bg\r
+  sprintf(key2,"%s_mr2",key) ;\r
+  TH1D * mr2 = new TH1D(key2,"Mass",nbin,xa) ;\r
+  sprintf(key2,"%s_sr2",key) ;\r
+  TH1D * sr2 = new TH1D(key2,"Width",nbin,xa) ;\r
+  sprintf(key2,"%s_ar2",key) ;\r
+  TH1D * ar2 = new TH1D(key2,"a",nbin,xa) ;\r
+  sprintf(key2,"%s_br2",key) ;\r
+  TH1D * br2 = new TH1D(key2,"a",nbin,xa) ;\r
+  sprintf(key2,"%s_cr2",key) ;\r
+  TH1D * cr2 = new TH1D(key2,"a",nbin,xa) ;\r
+  sprintf(key2,"%s_yr2",key) ;\r
+  TH1D * nr2 = new TH1D(key2,"Raw yield",nbin,xa) ;\r
+  sprintf(key2,"%s_yr2int",key) ;\r
+  TH1D * nr2int = new TH1D(key2,"Raw yield, integrated",nbin,xa) ;\r
+\r
+  TF1 * fit1 = new TF1("fit",CB,0.,1.,5) ;\r
+  fit->SetParName(0,"A") ;\r
+  fit->SetParName(1,"m_{0}") ;\r
+  fit->SetParName(2,"#sigma") ;\r
+  fit->SetParName(3,"a_{0}") ;\r
+  fit->SetParName(4,"a_{1}") ;\r
+  fit->SetLineWidth(2) ;\r
+  fit->SetLineColor(2) ;\r
+  TF1 * fgs = new TF1("gs",CBs,0.,1.,3) ;\r
+  fgs->SetLineColor(2) ;\r
+  fgs->SetLineWidth(1) ;\r
+\r
+  TF1 * fit2 = new TF1("fit2",CB2,0.,1.,6) ;\r
+  fit2->SetParName(0,"A") ;\r
+  fit2->SetParName(1,"m_{0}") ;\r
+  fit2->SetParName(2,"#sigma") ;\r
+  fit2->SetParName(3,"a_{0}") ;\r
+  fit2->SetParName(4,"a_{1}") ;\r
+  fit2->SetLineWidth(2) ;\r
+  fit2->SetLineColor(4) ;\r
+  fit2->SetLineStyle(2) ;\r
+  TF1 * fbg1 = new TF1("bg1",BG1,0.,1.,2) ;\r
+  TF1 * fbg2 = new TF1("bg2",BG2,0.,1.,3) ;\r
+\r
+  TCanvas * c3 = new TCanvas("mggFit1_Signal","mggFitCB",10,10,1400,800) ;\r
+  c3->Divide(6,3) ;\r
+\r
+  TCanvas * cReal = new TCanvas("cReal","Mgg real events",10,10,1400,800) ;\r
+  cReal->Divide(6,3) ;\r
+\r
+  TCanvas * c1 = new TCanvas("mggFit1","mggFit1",10,10,1400,800) ;\r
+  c1->Divide(6,3) ;\r
+  c1->cd(0) ;\r
+  TCanvas * c2=0,*c4=0,*c5=0,*c6=0 ; \r
+\r
+  TAxis * pta=h->GetYaxis() ;\r
+  TAxis * ma=h->GetXaxis() ;\r
+  for(Int_t i=1;i<=nbin;i++){\r
+    c1->cd(i) ;\r
+    Int_t imin=pta->FindBin(xa[i-1]+0.0001);\r
+    Int_t imax=pta->FindBin(xa[i]-0.0001) ;\r
+    Double_t pt=(xa[i]+xa[i-1])/2. ;\r
+    TH1D * hp = h->ProjectionX(Form("re%d",i),imin,imax) ;\r
+    hp->Sumw2() ;\r
+    TH1D * hpm= hm->ProjectionX(Form("mi%d",i),imin,imax) ;\r
+    hpm->Sumw2() ;\r
+    if(i<15){\r
+      hp ->Rebin(2) ;\r
+      hpm->Rebin(2) ;\r
+    }\r
+    else{\r
+      hp ->Rebin(5) ;\r
+      hpm->Rebin(5) ;\r
+    }\r
+    hp ->SetNdivisions(506);\r
+    hpm->SetNdivisions(506);\r
+    hp ->SetLabelSize(0.05,"X");\r
+    hpm->SetLabelSize(0.05,"X");\r
+    hp ->SetTitleSize(0.00,"X");\r
+    hpm->SetTitleSize(0.00,"X");\r
+    hp ->SetLabelSize(0.05,"Y");\r
+    hpm->SetLabelSize(0.05,"Y");\r
+    hp ->SetTitleSize(0.05,"Y");\r
+    hpm->SetTitleSize(0.05,"Y");\r
+    hp ->SetTitleOffset(1.0,"X");\r
+    hpm->SetTitleOffset(1.0,"X");\r
+    if (i>12) {\r
+      hp ->SetLabelSize(0.05,"X");\r
+      hpm->SetLabelSize(0.05,"X");\r
+      hp ->SetTitleSize(0.05,"X");\r
+      hpm->SetTitleSize(0.05,"X");\r
+      hp ->SetXTitle("m_{#gamma#gamma} (GeV/c^{2})");\r
+      hpm->SetXTitle("m_{#gamma#gamma} (GeV/c^{2})");\r
+    }\r
+\r
+    // for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hp ->GetBinContent(ib)==0)hp ->SetBinError(ib,1.);}\r
+    // for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hpm->GetBinContent(ib)==0)hpm->SetBinError(ib,1.);}\r
+    TH1D * hpm2   = (TH1D*)hpm->Clone("Bg1") ;\r
+    TH1D * hpcopy = (TH1D*)hp ->Clone("hpcopy") ;\r
+    TH1D * hp2    = (TH1D*)hp ->Clone("hp2") ;\r
+    TH1D * hpReal = (TH1D*)hp ->Clone("hpReal") ;\r
+    hpcopy->Divide(hpm) ;\r
+    sprintf(key,"%3.1f<p_{t}<%3.1f GeV/c",xa[i-1],xa[i]) ;\r
+    hpcopy->SetTitle(key) ;\r
+    hpReal->SetTitle(key) ;\r
+    hpReal->SetLineWidth(2);\r
+    hpcopy->SetMarkerStyle(20) ;\r
+    hpcopy->SetMarkerSize(0.7) ;\r
+    \r
+    fit1->SetParameters(TMath::Min(0.3,0.001*i*i),0.136,0.0077,0.0013,-0.0007) ;\r
+    fit1->SetParLimits(2,0.003,0.010) ;\r
+    fit1->SetParLimits(1,0.132,0.139) ;\r
+    fit1->SetParLimits(0,0.,1.) ;\r
+\r
+    Double_t rangeMin=0.06 ;\r
+    Double_t rangeMax=0.22 ;\r
+    if(i>=16)rangeMax=0.3 ; //More points to fix background\r
+    hpcopy->Fit(fit1,"NQ","",rangeMin,rangeMax) ;\r
+    hpcopy->Fit(fit1,"MQ","",rangeMin,rangeMax) ;\r
+    // hpcopy->Fit(fit1,"LQ","",rangeMin,rangeMax) ;\r
+    ar1->SetBinContent(i,fit1->GetParameter(3)) ;\r
+    ar1->SetBinError  (i,fit1->GetParError (3)) ;\r
+    br1->SetBinContent(i,fit1->GetParameter(4)) ;\r
+    br1->SetBinError  (i,fit1->GetParError (4)) ;\r
+\r
+    fit2->SetParameters(fit1->GetParameters()) ;\r
+    fit2->SetParLimits(2,0.003,0.010) ;\r
+    fit2->SetParLimits(1,0.132,0.139) ;\r
+    fit2->SetParLimits(0,0.,1.) ;\r
+    fit2->SetParameter(5,0.) ;\r
+\r
+    hpcopy->Fit(fit2,"+NQ","",rangeMin,rangeMax) ;\r
+    hpcopy->Fit(fit2,"+MQ","",rangeMin,rangeMax) ;\r
+    // hpcopy->Fit(fit2,"+LQ","",rangeMin,rangeMax) ;\r
+    ar2->SetBinContent(i,fit2->GetParameter(3)) ;\r
+    ar2->SetBinError  (i,fit2->GetParError (3)) ;\r
+    br2->SetBinContent(i,fit2->GetParameter(4)) ;\r
+    br2->SetBinError  (i,fit2->GetParError (4)) ;\r
+    cr2->SetBinContent(i,fit2->GetParameter(5)) ;\r
+    cr2->SetBinError  (i,fit2->GetParError (5)) ;\r
+\r
+    c1->cd(i) ;\r
+    hpcopy->SetAxisRange(0.065,0.229,"X") ;\r
+    hpcopy->Draw() ;\r
+\r
+    cReal->cd(i) ;\r
+    hpReal->SetAxisRange(0.065,0.229,"X") ;\r
+    hpReal->Draw("ehist") ;\r
+\r
+    fbg1->SetParameters(fit1->GetParameter(3),fit1->GetParameter(4)); \r
+    fbg2->SetParameters(fit2->GetParameter(3),fit2->GetParameter(4),fit2->GetParameter(5)); \r
+\r
+    Double_t intRangeMin = PeakPosition(pt)-nSigma*PeakWidth(pt) ;\r
+    Double_t intRangeMax = PeakPosition(pt)+nSigma*PeakWidth(pt) ;\r
+    Int_t    intBinMin   = hp->GetXaxis()->FindBin(intRangeMin) ;\r
+    Int_t    intBinMax   = hp->GetXaxis()->FindBin(intRangeMax) ;\r
+    Double_t errStat     = hpm->Integral(intBinMin,intBinMax); \r
+\r
+    hpm ->Multiply(fbg1) ;\r
+    hpm2->Multiply(fbg2) ;\r
+    hp  ->Add(hpm,-1.) ;\r
+    hp2 ->Add(hpm2,-1.) ;\r
+\r
+    c3->cd(i) ;\r
+    \r
+    if(i<15)\r
+      fgs->SetParameters(hp->Integral(32,36)/5.,fit1->GetParameter(1),fit1->GetParameter(2)) ;\r
+    else\r
+      fgs->SetParameters(hp->Integral(13,15)/3.,0.136,0.005) ;\r
+    fgs->SetParLimits(2,0.003,0.010) ;\r
+    fgs->SetParLimits(1,0.132,0.142) ;\r
+    fgs->SetParLimits(0,0.,1.e+6) ;\r
+    \r
+    hp->Fit(fgs,"Q","",rangeMin,rangeMax) ;   \r
+    hp->SetMaximum(hp2->GetMaximum()*1.1) ;\r
+    hp->SetMinimum(hp2->GetMinimum()*1.1) ;\r
+    hp->SetMarkerStyle(20) ;\r
+    hp->SetMarkerSize(0.7) ;\r
+    mr1->SetBinContent(i,fgs->GetParameter(1)) ;\r
+    mr1->SetBinError  (i,fgs->GetParError(1)) ;\r
+    sr1->SetBinContent(i,TMath::Abs(fgs->GetParameter(2))) ;\r
+    sr1->SetBinError  (i,fgs->GetParError(2)) ;\r
+\r
+    Double_t y=fgs->Integral(intRangeMin,intRangeMax)/hp->GetXaxis()->GetBinWidth(1) ;\r
+    nr1->SetBinContent(i,y) ;\r
+    Double_t ey=0 ;\r
+    if(fgs->GetParameter(0)!=0. && fgs->GetParameter(2)!=0.){\r
+      Double_t en=fgs->GetParError(0)/fgs->GetParameter(0) ;\r
+      Double_t es=fgs->GetParError(2)/fgs->GetParameter(2) ;\r
+      ey=y*TMath::Sqrt(en*en+es*es) ;\r
+    }\r
+    nr1->SetBinError(i,ey) ;\r
+\r
+    Double_t npiInt=hp->Integral(intBinMin,intBinMax) ;\r
+    Double_t norm=fbg1->GetParameter(0) ;\r
+    Double_t normErr=fbg1->GetParError(0) ;\r
+    if(npiInt>0.){\r
+      nr1int->SetBinContent(i,npiInt) ;\r
+      nr1int->SetBinError(i,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;\r
+    }\r
+    hp2->GetXaxis()->SetRangeUser(0.06,0.229) ;\r
+    hp2->SetMaximum(hp2->GetMaximum()*1.1) ;\r
+    hp2->SetMinimum(hp2->GetMinimum()*1.1) ;\r
+    hp2->SetMarkerStyle(20) ;\r
+    hp2->SetMarkerSize(0.7) ;\r
+\r
+    hp2->Fit(fgs,"Q","",rangeMin,rangeMax) ;\r
+    mr2->SetBinContent(i,fgs->GetParameter(1)) ;\r
+    mr2->SetBinError  (i,fgs->GetParError(1)) ;\r
+    sr2->SetBinContent(i,TMath::Abs(fgs->GetParameter(2))) ;\r
+    sr2->SetBinError  (i,fgs->GetParError(2)) ;\r
+    y=fgs->Integral(intRangeMin,intRangeMax)/hp->GetXaxis()->GetBinWidth(1) ;\r
+    nr2->SetBinContent(i,y) ;\r
+    ey=0 ;\r
+    if(fgs->GetParameter(0)!=0. && fgs->GetParameter(2)!=0.){\r
+      Double_t en=fgs->GetParError(0)/fgs->GetParameter(0) ;\r
+      Double_t es=fgs->GetParError(2)/fgs->GetParameter(2) ;\r
+      ey=y*TMath::Sqrt(en*en+es*es) ;\r
+    }\r
+    nr2->SetBinError(i,ey) ;\r
+    npiInt=hp2->Integral(intBinMin,intBinMax) ;\r
+    norm=fbg2->GetParameter(0) ;\r
+    normErr=fbg2->GetParError(0) ;\r
+    if(npiInt>0.){\r
+      nr2int->SetBinContent(i,npiInt) ;\r
+      nr2int->SetBinError(i,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;\r
+    } \r
+    hp2->SetTitle(key) ;\r
+    hp2->SetAxisRange(0.065,0.229,"X") ;\r
+    hp2->Draw() ;\r
+    hp ->Draw("same") ;\r
+    delete hpm ;\r
+    delete hpm2 ;\r
+    c1->Update() ;\r
+    c3->Update() ;\r
+    cReal->Update() ;\r
+  }\r
+\r
+  if (c3) c3->Print("Pi0_Signal_Gaus.eps") ;\r
+  if (c1) c1->Print("Pi0_Ratio_Gaus.eps") ;\r
+  if (cReal) cReal->Print("Pi0_InvMass_Gaus.eps") ;\r
+\r
+  //Normalize by the number of events\r
+  Double_t nMBOR   = hev->GetBinContent(2); // MBOR events pass4\r
+  Double_t nevents = nMBOR;\r
+  printf("==============\nN events = %d\n==============\n",nevents);\r
+\r
+  nr1   ->Scale(1./nevents) ;\r
+  nr1int->Scale(1./nevents) ;\r
+  nr2   ->Scale(1./nevents) ;\r
+  nr2int->Scale(1./nevents) ;\r
+\r
+\r
+  TFile fout("LHC11a_FitResult.root","update");\r
+  mr1->Write(0,TObject::kOverwrite) ;\r
+  sr1->Write(0,TObject::kOverwrite) ;\r
+  ar1->Write(0,TObject::kOverwrite) ;\r
+  br1->Write(0,TObject::kOverwrite) ;\r
+  nr1->Write(0,TObject::kOverwrite) ;\r
+  nr1int->Write(0,TObject::kOverwrite) ;\r
+  ar2->Write(0,TObject::kOverwrite) ;\r
+  br2->Write(0,TObject::kOverwrite) ;\r
+  cr2->Write(0,TObject::kOverwrite) ;\r
+  mr2->Write(0,TObject::kOverwrite) ;\r
+  sr2->Write(0,TObject::kOverwrite) ;\r
+  nr2->Write(0,TObject::kOverwrite) ;\r
+  nr2int->Write(0,TObject::kOverwrite) ;\r
+  fout.Close() ;\r
+\r
+}\r
+\r
+//-----------------------------------------------------------------------------\r
+const Double_t kMean=0.135 ; //Approximate peak position to facilitate error estimate\r
+\r
+Double_t PeakPosition(Double_t pt){\r
+  //Fit to the measured peak position\r
+  return 0.454*exp(-pt*7.13)+0.13569 ;\r
+}\r
+Double_t PeakWidth(Double_t pt){\r
+  //fit to the measured peak width\r
+  return 1.69e-02*exp(-pt*2.82)+4.93e-03 ;\r
+}\r
\r
+Double_t CB(Double_t * x, Double_t * par){\r
+  //Parameterization of Real/Mixed ratio\r
+  Double_t m=par[1] ;\r
+  Double_t s=par[2] ;\r
+  Double_t dx=(x[0]-m)/s ;\r
+  return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean) ;\r
+}\r
+Double_t CB2(Double_t * x, Double_t * par){\r
+  //Another parameterization of Real/Mixed ratio\r
+  Double_t m=par[1] ;\r
+  Double_t s=par[2] ;\r
+  Double_t dx=(x[0]-m)/s ;\r
+  return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean)+par[5]*(x[0]-kMean)*(x[0]-kMean) ;\r
+}\r
+Double_t CBs(Double_t * x, Double_t * par){\r
+  //Parameterizatin of signal\r
+  Double_t m=par[1] ;\r
+  Double_t s=par[2] ;\r
+  Double_t dx=(x[0]-m)/s ;\r
+  return par[0]*exp(-dx*dx/2.) ;\r
+}\r
+Double_t BG1(Double_t * x, Double_t * par){\r
+  //Normalizatino of Mixed\r
+  return par[0]+par[1]*(x[0]-kMean) ;\r
+}\r
+Double_t BG2(Double_t * x, Double_t * par){\r
+  //Another normalization of  Mixed\r
+  return par[0]+par[1]*(x[0]-kMean)+par[2]*(x[0]-kMean)*(x[0]-kMean) ;\r
+}\r
+\r
+\r
+//-----------------------------------------------------------------------------\r
+PPRstyle()\r
+{\r
+\r
+  //////////////////////////////////////////////////////////////////////\r
+  //\r
+  // ROOT style macro for the TRD TDR\r
+  //\r
+  //////////////////////////////////////////////////////////////////////\r
+\r
+  gStyle->SetPalette(1);\r
+  gStyle->SetCanvasBorderMode(-1);\r
+  gStyle->SetCanvasBorderSize(1);\r
+  gStyle->SetCanvasColor(10);\r
+\r
+  gStyle->SetFrameFillColor(10);\r
+  gStyle->SetFrameBorderSize(1);\r
+  gStyle->SetFrameBorderMode(-1);\r
+  gStyle->SetFrameLineWidth(1.2);\r
+  gStyle->SetFrameLineColor(1);\r
+\r
+  gStyle->SetHistFillColor(0);\r
+  gStyle->SetHistLineWidth(1);\r
+  gStyle->SetHistLineColor(1);\r
+\r
+  gStyle->SetPadColor(10);\r
+  gStyle->SetPadBorderSize(1);\r
+  gStyle->SetPadBorderMode(-1);\r
+\r
+  gStyle->SetStatColor(10);\r
+  gStyle->SetTitleColor(kBlack,"X");\r
+  gStyle->SetTitleColor(kBlack,"Y");\r
+\r
+  gStyle->SetLabelSize(0.04,"X");\r
+  gStyle->SetLabelSize(0.04,"Y");\r
+  gStyle->SetLabelSize(0.04,"Z");\r
+  gStyle->SetTitleSize(0.04,"X");\r
+  gStyle->SetTitleSize(0.04,"Y");\r
+  gStyle->SetTitleSize(0.04,"Z");\r
+  gStyle->SetTitleFont(42,"X");\r
+  gStyle->SetTitleFont(42,"Y");\r
+  gStyle->SetTitleFont(42,"X");\r
+  gStyle->SetLabelFont(42,"X");\r
+  gStyle->SetLabelFont(42,"Y");\r
+  gStyle->SetLabelFont(42,"Z");\r
+  gStyle->SetStatFont(42);\r
+\r
+  gStyle->SetTitleOffset(1.0,"X");\r
+  gStyle->SetTitleOffset(1.4,"Y");\r
+\r
+  gStyle->SetFillColor(kWhite);\r
+  gStyle->SetTitleFillColor(kWhite);\r
+\r
+  gStyle->SetOptDate(0);\r
+  gStyle->SetOptTitle(1);\r
+  gStyle->SetOptStat(0);\r
+  gStyle->SetOptFit(0);\r
+\r
+}\r
+\r