PWGGA/PHOSTasks/PHOS_PbPb/macros/production/MakeMmixPi0.C

   1 /* $Id$ */
   2
   3 #include "TFile.h"
   4 #include "TF1.h"
   5 #include "TH1F.h"
   6 #include "TH2F.h"
   7 #include "TSystem.h"
   8 #include "TStyle.h"
   9 #include "TCanvas.h"
  10 #include "TMath.h"
  11
  12 void PPRstyle();
  13
  14 Double_t PeakPosition(Double_t pt);
  15 Double_t PeakWidth(Double_t pt);
  16 Double_t CB(Double_t * x, Double_t * par);
  17 Double_t CB2(Double_t * x, Double_t * par);
  18 Double_t CBs(Double_t * x, Double_t * par);
  19 Double_t BG1(Double_t * x, Double_t * par);
  20 Double_t BG2(Double_t * x, Double_t * par);
  21
  22
  23 // const Int_t nPadX = 3, nPadY = 2;
  24 // const Int_t nPtBins=6 ;
  25 // const Double_t ptBinEdges[21]={1., 2., 3., 4., 5., 7., 10.} ;
  26
  27 const Int_t nPadX = 5, nPadY = 4;
  28 Int_t nPtBins=0;
  29 Double_t ptBinEdges[1000] = {0};
  30 double GetPtBin(int bin){
  31   if( bin==0 )
  32     return 1.;
  33
  34   // return GetPtBin(bin-1) * 1.1;
  35
  36   // if ( bin % 2 )
  37   //   return GetPtBin(bin-1) + 0.4;
  38   // else
  39   //   return GetPtBin(bin-1) + 0.2;
  40
  41   double previousBin = GetPtBin(bin-1);
  42   double linInc = 0.2;
  43   double threshold = 5.;
  44   double logFact = 1 + linInc/threshold;
  45   if ( previousBin < threshold ) // linear
  46     return previousBin + linInc;
  47   else { // logarithmic
  48     return previousBin * logFact;
  49   }
  50 }
  51 void MakePtBins() {
  52   int bin = -1;
  53   do {
  54     ++bin;
  55     ptBinEdges[bin] = GetPtBin(bin);
  56   } while(ptBinEdges[bin] < 40.);
  57   nPtBins = bin -2;
  58
  59   printf("Making Pt Bins:\n");
  60   for(int b=0; b < nPtBins+1; ++b)
  61     printf("%.1f, ", ptBinEdges[b]);
  62   printf("\n N. Bins: %d\n", nPtBins);
  63
  64
  65   // for(int bin = 0; bin <= nPtBins; ++bin){
  66   //   ptBinEdges[bin] = GetPtBin(bin);
  67   //   printf("%.1f, ", ptBinEdges[bin]);
  68   // }
  69   // printf("\n");
  70 }
  71
  72
  73 const int kNCentralityBins = 3;
  74
  75 const Double_t kMean=0.136 ; //Approximate peak position to facilitate error estimate
  76
  77 const char format[] = ".pdf"; // say if you want .pdf
  78
  79 TH2F* FindPi0(TList *histoList, bool mix, const char* pid, int centrality);
  80
  81 //-----------------------------------------------------------------------------
  82 void MakeMmixPi0(const TString filename,
  83                  const TString listPath = "PHOSPi0Flow/PHOSPi0FlowCoutput1", // lego train
  84                  const Int_t centrality=0,
  85                  const char* pid="CPV",
  86                  const char* trigger="kCentral",
  87                  const char* saveToDir="")
  88 {
  89   MakePtBins();
  90   Printf("\nMakeMmixPi0(%s, %s, %i, %s, %s)", filename.Data(), listPath.Data(), centrality, pid, saveToDir);
  91
  92   if( TString(saveToDir).Length() )
  93     gSystem->mkdir(saveToDir, true);
  94
  95   //Fit Real/Mixed ratio, normalize Mixed and subtract it from Real
  96
  97   TFile * file = new TFile(filename) ;
  98   TList *histoList = (TList*)file->Get(listPath);
  99
 100   char key[125] ;
 101
 102   TH1F * hev = (TH1F*)histoList->FindObject("hTotSelEvents") ;
 103   TH2F * hCentrality  = (TH2F*)histoList->FindObject("hCenPHOSCells") ;
 104   TH1D * hCentralityX = hCentrality->ProjectionX();
 105
 106   printf("TotSelEvents (4): %.0f \n",hev->GetBinContent(4)) ;
 107   printf("Centrality:   %.0f \n",hCentralityX->Integral()) ;
 108
 109   TH2F *hPi0 = FindPi0(histoList, false, pid, centrality);
 110   TH2F *hPi0Mix = FindPi0(histoList, true, pid, centrality);
 111   if( !hPi0 || !hPi0Mix || hPi0->GetEntries() < 10000) {
 112     Printf(Form("no histogram(0x%p, 0x%p) or to low number of entries(%.1f) in hPi0, skipping this cent", hPi0, hPi0Mix, hPi0->GetEntries()));
 113     file->Close();
 114     delete file;
 115     return;
 116   }
 117
 118
 119   TFile* outFile = new TFile(Form("%sPi0_FitResult.root", saveToDir, centrality),"update");
 120
 121   PPRstyle();
 122   gStyle->SetPadLeftMargin(0.14);
 123   gStyle->SetPadRightMargin(0.01);
 124   //gStyle->SetPadTopMargin(0.01);
 125   gStyle->SetPadBottomMargin(0.08);
 126
 127   //Fit real only
 128   //Linear Bg
 129   char kkey[55];
 130   sprintf(kkey, Form("%s_cent%d_%s",pid, centrality, trigger)) ;
 131   char key2[155];
 132   sprintf(key,"Mix%s",kkey) ;
 133   sprintf(key2,"%s_mr1",key) ;
 134   TH1D * mr1 = new TH1D(key2,"Mass",nPtBins,ptBinEdges) ;
 135   sprintf(key2,"%s_sr1",key) ;
 136   TH1D * sr1 = new TH1D(key2,"Width",nPtBins,ptBinEdges) ;
 137   sprintf(key2,"%s_ar1",key) ;
 138   TH1D * ar1 = new TH1D(key2,"a",nPtBins,ptBinEdges) ;
 139   sprintf(key2,"%s_br1",key) ;
 140   TH1D * br1 = new TH1D(key2,"a",nPtBins,ptBinEdges) ;
 141   sprintf(key2,"%s_yr1",key) ;
 142   TH1D * nr1 = new TH1D(key2,"Raw yield",nPtBins,ptBinEdges) ;
 143   sprintf(key2,"%s_yr1int",key) ;
 144   TH1D * nr1int = new TH1D(key2,"Raw yield, integrated",nPtBins,ptBinEdges) ;
 145
 146   //Quadratic Bg
 147   sprintf(key2,"%s_mr2",key) ;
 148   TH1D * mr2 = new TH1D(key2,"Mass",nPtBins,ptBinEdges) ;
 149   sprintf(key2,"%s_sr2",key) ;
 150   TH1D * sr2 = new TH1D(key2,"Width",nPtBins,ptBinEdges) ;
 151   sprintf(key2,"%s_ar2",key) ;
 152   TH1D * ar2 = new TH1D(key2,"a",nPtBins,ptBinEdges) ;
 153   sprintf(key2,"%s_br2",key) ;
 154   TH1D * br2 = new TH1D(key2,"a",nPtBins,ptBinEdges) ;
 155   sprintf(key2,"%s_cr2",key) ;
 156   TH1D * cr2 = new TH1D(key2,"a",nPtBins,ptBinEdges) ;
 157   sprintf(key2,"%s_yr2",key) ;
 158   TH1D * nr2 = new TH1D(key2,"Raw yield",nPtBins,ptBinEdges) ;
 159   sprintf(key2,"%s_yr2int",key) ;
 160   TH1D * nr2int = new TH1D(key2,"Raw yield, integrated",nPtBins,ptBinEdges) ;
 161
 162   // Functions:
 163   TF1 * funcRatioFit1 = new TF1("funcRatioFit1",CB,0.,1.,6) ;
 164   funcRatioFit1->SetParName(0,"A") ;
 165   funcRatioFit1->SetParName(1,"m_{0}") ;
 166   funcRatioFit1->SetParName(2,"#sigma") ;
 167   funcRatioFit1->SetParName(3,"a_{0}") ;
 168   funcRatioFit1->SetParName(4,"a_{1}") ;
 169   funcRatioFit1->SetParName(5,"a_{2}") ;
 170   funcRatioFit1->SetLineWidth(2) ;
 171   funcRatioFit1->SetLineColor(2) ;
 172   TF1 * fgs = new TF1("gs",CBs,0.,1.,4) ;
 173   fgs->SetParName(0,"A") ;
 174   fgs->SetParName(1,"m_{0}") ;
 175   fgs->SetParName(2,"#sigma") ;
 176   fgs->SetParName(3,"B") ;
 177   fgs->SetLineColor(2) ;
 178   fgs->SetLineWidth(1) ;
 179   TF1 * funcRatioFit2 = new TF1("funcRatioFit2",CB2,0.,1.,7) ;
 180   funcRatioFit2->SetParName(0,"A") ;
 181   funcRatioFit2->SetParName(1,"m_{0}") ;
 182   funcRatioFit2->SetParName(2,"#sigma") ;
 183   funcRatioFit2->SetParName(3,"a_{0}") ;
 184   funcRatioFit2->SetParName(4,"a_{1}") ;
 185   funcRatioFit2->SetParName(5,"a_{2}") ;
 186   funcRatioFit2->SetParName(6,"a_{3}") ;
 187   funcRatioFit2->SetLineWidth(2) ;
 188   funcRatioFit2->SetLineColor(4) ;
 189   funcRatioFit2->SetLineStyle(2) ;
 190   TF1 * fbg1 = new TF1("bg1",BG1,0.,1.,3) ;
 191   TF1 * fbg2 = new TF1("bg2",BG2,0.,1.,4) ;
 192
 193
 194   // Canvases
 195   TCanvas * rawCanvas = new TCanvas("rawCanvas","rawCanvas",10,10,1200,800);
 196   rawCanvas->Divide(nPadX, nPadY);
 197   TCanvas * ratioCanvas = new TCanvas("mggFit1","mggFit1",10,10,1200,800) ;
 198   ratioCanvas->Divide(nPadX,nPadY) ;
 199   TCanvas * sbsCanvas = new TCanvas("mggFit1_Signal","mggFitCB",10,10,1200,800) ;
 200   sbsCanvas->Divide(nPadX,nPadY) ;
 201
 202
 203   for(Int_t ptBin=1; ptBin<=nPtBins; ptBin++){
 204     ratioCanvas->cd(ptBin) ;
 205     TAxis * pta=hPi0->GetYaxis() ;
 206     Int_t imin=pta->FindBin(ptBinEdges[ptBin-1]+0.0001);
 207     Int_t imax=pta->FindBin(ptBinEdges[ptBin]-0.0001) ;
 208     Double_t pt=(ptBinEdges[ptBin]+ptBinEdges[ptBin-1])/2. ;
 209
 210     TH1D * hPi0Proj = hPi0->ProjectionX(Form("hPi0_ptBin%d",ptBin), imin, imax) ;
 211     TH1D * hPi0ProjMix= hPi0Mix->ProjectionX(Form("hPi0Mix_ptBin%d",ptBin), imin, imax) ;
 212
 213     hPi0Proj->SetTitle(Form("M_{#gamma#gamma}, PID=%s, %.1f<p_{T}<%.1f GeV/c",pid,ptBinEdges[ptBin-1],ptBinEdges[ptBin]));
 214     hPi0Proj->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})");
 215     hPi0ProjMix->SetTitle(Form("M_{#gamma#gamma}^{Mix}, PID=%s, %.1f<p_{T}<%.1f GeV/c",pid,ptBinEdges[ptBin-1],ptBinEdges[ptBin]));
 216     hPi0ProjMix->SetXTitle("M_{#gamma#gamma}^{Mix} (GeV/c^{2})");
 217
 218     Printf("\n\t%.1f<pt<%.1f GeV/c, entries: %.0f",ptBinEdges[ptBin-1],ptBinEdges[ptBin],  hPi0Proj->GetEntries());
 219     if( hPi0Proj->GetEntries() < 100. ) {
 220       Printf("skipping this bin as n. entries is to low");
 221       continue;
 222     }
 223
 224     hPi0Proj ->Rebin(4) ;
 225     hPi0ProjMix->Rebin(4) ;
 226     hPi0Proj->SetName(Form("%s_rebin", hPi0Proj->GetName()));
 227     hPi0ProjMix->SetName(Form("%s_rebin", hPi0ProjMix->GetName()));
 228
 229     // Error Fix
 230     for(Int_t ib=1; ib<=hPi0Proj->GetNbinsX();ib++){if(hPi0Proj ->GetBinContent(ib)==0)hPi0Proj ->SetBinError(ib,1.);}
 231     for(Int_t ib=1; ib<=hPi0Proj->GetNbinsX();ib++){if(hPi0ProjMix->GetBinContent(ib)==0)hPi0ProjMix->SetBinError(ib,1.);}
 232
 233     // Signal-Mix Ratio
 234     TH1D * hPi0Ratio = (TH1D*)hPi0Proj->Clone( Form("hPi0Ratio_ptBin%d",ptBin) ) ;
 235     hPi0Ratio->SetTitle(Form("#frac{M_{#gamma#gamma}}{M_{#gamma#gamma}^{Mix}}, %.1f<p_{T}<%.1f GeV/c", ptBinEdges[ptBin-1], ptBinEdges[ptBin]));
 236     hPi0Ratio->Divide(hPi0ProjMix) ;
 237     hPi0Ratio->SetMarkerStyle(20) ;
 238     hPi0Ratio->SetMarkerSize(0.7) ;
 239
 240     funcRatioFit1->SetParameters(0.001,0.136,0.0055,0.0002,-0.002,0.0) ;
 241     funcRatioFit1->SetParLimits(0,0.000,1.000) ;
 242     funcRatioFit1->SetParLimits(1,0.120,0.145) ;
 243     funcRatioFit1->SetParLimits(2,0.004,0.012) ;
 244
 245     Double_t rangeMin=0.05 ;
 246     Double_t rangeMax=0.3 ;
 247     if(centrality==0) rangeMax=0.4 ;
 248     if(ptBin==1){
 249       rangeMin=0.06 ;
 250       rangeMax=0.25 ;
 251     }
 252     int error = hPi0Ratio->Fit(funcRatioFit1,"Q" ,"",rangeMin,rangeMax) ;
 253     if( error ) {
 254       Printf("fit (ratio) error: %d ", error);
 255       continue;
 256     }
 257     error = hPi0Ratio->Fit(funcRatioFit1,"MQ","",rangeMin,rangeMax) ;
 258     if( error % 4000) {
 259       Printf("fit (ratio more) error: %d ", error);
 260       continue;
 261     }
 262
 263
 264     ar1->SetBinContent(ptBin,funcRatioFit1->GetParameter(3)) ;
 265     ar1->SetBinError  (ptBin,funcRatioFit1->GetParError(3)) ;
 266     br1->SetBinContent(ptBin,funcRatioFit1->GetParameter(4)) ;
 267     br1->SetBinError  (ptBin,funcRatioFit1->GetParError(4)) ;
 268
 269     funcRatioFit2->SetParameters(funcRatioFit1->GetParameters()) ;
 270     funcRatioFit2->SetParLimits(0,0.000,1.000) ;
 271     funcRatioFit2->SetParLimits(1,0.110,0.145) ;
 272     funcRatioFit2->SetParLimits(2,0.004,0.012) ;
 273
 274     error = hPi0Ratio->Fit(funcRatioFit2,"+NQ","",rangeMin,rangeMax) ;
 275     if( error )  {
 276       Printf("fit (funcRatioFit2) error: %d ", error);
 277       continue;
 278     }
 279     error =hPi0Ratio->Fit(funcRatioFit2,"+MQ","",rangeMin,rangeMax) ;
 280     if( error % 4000) {
 281       Printf("fit (funcRatioFit2 more) error: %d ", error);
 282       continue;
 283     }
 284
 285     ar2->SetBinContent(ptBin,funcRatioFit2->GetParameter(3)) ;
 286     ar2->SetBinError  (ptBin,funcRatioFit2->GetParError(3)) ;
 287     br2->SetBinContent(ptBin,funcRatioFit2->GetParameter(4)) ;
 288     br2->SetBinError  (ptBin,funcRatioFit2->GetParError(4)) ;
 289     cr2->SetBinContent(ptBin,funcRatioFit2->GetParameter(5)) ;
 290     cr2->SetBinError  (ptBin,funcRatioFit2->GetParError(5)) ;
 291     hPi0Ratio->GetXaxis()->SetRangeUser(0.05,0.30) ;
 292     hPi0Ratio->DrawCopy() ;
 293     hPi0Ratio->Write(0,TObject::kOverwrite) ;
 294     ratioCanvas->Update() ;
 295
 296     fbg1->SetParameters(funcRatioFit1->GetParameter(3),funcRatioFit1->GetParameter(4),funcRatioFit1->GetParameter(5));
 297     fbg2->SetParameters(funcRatioFit2->GetParameter(3),funcRatioFit2->GetParameter(4),funcRatioFit2->GetParameter(5),
 298                         funcRatioFit2->GetParameter(6));
 299
 300     Double_t intRangeMin = PeakPosition(pt)-3.*PeakWidth(pt) ;
 301     Double_t intRangeMax = PeakPosition(pt)+3.*PeakWidth(pt) ;
 302     Int_t    intBinMin   = hPi0Proj->GetXaxis()->FindBin(intRangeMin) ;
 303     Int_t    intBinMax   = hPi0Proj->GetXaxis()->FindBin(intRangeMax) ;
 304     Double_t errStat     = hPi0ProjMix->Integral(intBinMin,intBinMax);
 305
 306     rawCanvas->cd(ptBin);
 307     TH1D * hPi0MixScaled = (TH1D*)hPi0ProjMix ->Clone(Form("%s_clone2", hPi0Proj->GetName())) ;
 308     hPi0MixScaled->Scale(fbg1->Eval(0.1349));
 309     hPi0Proj->SetLineColor(kBlack);
 310     hPi0Proj->SetAxisRange(0.0, 0.3);
 311     hPi0Proj->DrawCopy();
 312     hPi0Proj->Write(0,TObject::kOverwrite) ;
 313     hPi0MixScaled->SetLineColor(kRed);
 314     hPi0MixScaled->SetTitle(Form("M_{#gamma#gamma}^{Mix,scaled}, PID=%s, %.1f<p_{T}<%.1f GeV/c",pid,ptBinEdges[ptBin-1],ptBinEdges[ptBin]));
 315     hPi0MixScaled->DrawCopy("same");
 316     rawCanvas->Update();
 317     hPi0MixScaled->Write(0,TObject::kOverwrite) ;
 318
 319     // Linear background subtraction
 320     TH1D * hPi0MixScaledPol1 = (TH1D*)hPi0ProjMix->Clone(Form("hPi0MixScaledPol1_ptBin%d", ptBin)) ;
 321     hPi0MixScaledPol1 ->Multiply(fbg1) ;
 322     TH1D * hPi0BSPol1 = (TH1D*)hPi0Proj->Clone(Form("hPi0BSPol1_ptBin%d", ptBin)) ;
 323     hPi0BSPol1->Add(hPi0MixScaledPol1 ,-1.) ;
 324
 325     // Quadratic background
 326     TH1D * hPi0MixScaledPol2 = (TH1D*)hPi0ProjMix->Clone(Form("hPi0MixScaledPol2_ptBin%d", ptBin)) ;
 327     hPi0MixScaledPol2->Multiply(fbg2) ;
 328     TH1D * hPi0BSPol2     = (TH1D*)hPi0Proj    ->Clone(Form("hPi0BSPol2_ptBin%d", ptBin)) ;
 329     hPi0BSPol2 ->Add(hPi0MixScaledPol2,-1.) ;
 330
 331     sbsCanvas->cd(ptBin) ;
 332
 333     Int_t binPi0 = hPi0BSPol1->FindBin(funcRatioFit1->GetParameter(1));
 334     Int_t nWidPi0 = 2 * (Int_t) (funcRatioFit1->GetParameter(2)/hPi0BSPol1->GetBinWidth(1));
 335     fgs->SetParameters(hPi0BSPol1->Integral(binPi0-nWidPi0,binPi0+nWidPi0)/5., funcRatioFit1->GetParameter(1), funcRatioFit1->GetParameter(2)) ;
 336     fgs->SetParLimits(0,0.,1.e+5) ;
 337     fgs->SetParLimits(1,0.110,0.145) ;
 338     fgs->SetParLimits(2,0.004,0.02) ;
 339     error = hPi0BSPol1->Fit(fgs,"Q","",rangeMin,rangeMax) ;
 340     if( error ) {
 341       Printf("fit (hPi0BSPol1 fgs) error: %d ", error);
 342       continue;
 343     }
 344     hPi0BSPol1->SetMaximum(hPi0BSPol2->GetMaximum()*1.5) ;
 345     hPi0BSPol1->SetMinimum(hPi0BSPol2->GetMinimum()*1.1) ;
 346     hPi0BSPol1->SetMarkerStyle(20) ;
 347     hPi0BSPol1->SetMarkerSize(0.7) ;
 348     mr1->SetBinContent(ptBin,fgs->GetParameter(1)) ;
 349     mr1->SetBinError  (ptBin,fgs->GetParError(1) ) ;
 350     sr1->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ;
 351     sr1->SetBinError  (ptBin,fgs->GetParError(2) ) ;
 352
 353     Double_t y=fgs->GetParameter(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ;
 354     nr1->SetBinContent(ptBin,y) ;
 355     Double_t ey=fgs->GetParError(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ;
 356     nr1->SetBinError(ptBin,ey) ;
 357
 358     Double_t npiInt = hPi0BSPol1->Integral(intBinMin,intBinMax) ;
 359     Double_t norm   = fbg1->GetParameter(0) ;
 360     Double_t normErr= fbg1->GetParError(0) ;
 361     if(npiInt>0.){
 362       nr1int->SetBinContent(ptBin,npiInt) ;
 363       nr1int->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
 364     }
 365     hPi0BSPol2->GetXaxis()->SetRangeUser(0.05,0.3) ;
 366     hPi0BSPol2->SetMaximum(hPi0BSPol2->GetMaximum()*1.5) ;
 367     hPi0BSPol2->SetMinimum(hPi0BSPol2->GetMinimum()*1.1) ;
 368     hPi0BSPol2->SetMarkerStyle(20) ;
 369     hPi0BSPol2->SetMarkerSize(0.7) ;
 370     hPi0BSPol2->Fit(fgs,"Q","",rangeMin,rangeMax) ;
 371     if( error ) {
 372       Printf("fit (hPi0BSPol2 fgs) error: %d ", error);
 373       continue;
 374     }
 375     mr2->SetBinContent(ptBin,fgs->GetParameter(1)) ;
 376     mr2->SetBinError  (ptBin,fgs->GetParError(1)) ;
 377     sr2->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ;
 378     sr2->SetBinError  (ptBin,fgs->GetParError(2)) ;
 379     y=fgs->GetParameter(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ;
 380     nr2->SetBinContent(ptBin,y) ;
 381     ey= fgs->GetParError(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ;
 382     nr2->SetBinError(ptBin,ey) ;
 383     Double_t npiInt2=hPi0BSPol2->Integral(intBinMin,intBinMax) ;
 384     norm=fbg2->GetParameter(0) ;
 385     normErr=fbg2->GetParError(0) ;
 386     if(npiInt2>0.){
 387       nr2int->SetBinContent(ptBin,npiInt2) ;
 388       nr2int->SetBinError(ptBin,TMath::Sqrt(npiInt2 + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
 389     }
 390     hPi0BSPol2->SetTitle(Form("M_{#gamma#gamma}^{BS_{2}}, PID=%s, %.1f<p_{T}<%.1f GeV/c",pid,ptBinEdges[ptBin-1],ptBinEdges[ptBin]));
 391     hPi0BSPol2->DrawCopy() ;
 392     sbsCanvas->Update() ;
 393     hPi0BSPol2->Write(0,TObject::kOverwrite) ;
 394   }
 395
 396   char name[55] ;
 397   sprintf(name,"%sPi0_ratio_%s%s", saveToDir, kkey, format) ;
 398   ratioCanvas->Print(name) ;
 399   sprintf(name,"%sPi0_signal_%s%s", saveToDir, kkey, format) ;
 400   sbsCanvas->Print(name) ;
 401   sprintf(name,"%sPi0_raw_%s%s", saveToDir, kkey, format) ;
 402   rawCanvas->Print(name);
 403
 404   //Normalize by the number of events
 405   Int_t cMin=0, cMax=100;
 406   if      (centrality == 0) {
 407     cMin=1;
 408     cMax=10;
 409   }
 410   else if (centrality == 1) {
 411     cMin=11;
 412     cMax=40;
 413   }
 414   else if (centrality == 2) {
 415     cMin=41;
 416     cMax=80;
 417   }
 418   else if (centrality == -1) {
 419     cMin=1;
 420     cMax=80;
 421   }
 422   else {
 423     Printf("ERROR: Centrality: %d not defined !!! ERROR", centrality);
 424     return;
 425   }
 426
 427   Double_t nevents = hCentralityX->Integral(cMin,cMax);
 428   if ( nevents > 0.9 ) {
 429     nr1   ->Scale(1./nevents) ;
 430     nr1int->Scale(1./nevents) ;
 431     nr2   ->Scale(1./nevents) ;
 432     nr2int->Scale(1./nevents) ;
 433   } else {
 434     Printf("WARNING: non positive nEvents in centrality range, cMin:%d, cMax:%d, nEvents:%f", cMin, cMax, nevents );
 435
 436   }
 437
 438   mr1->Write(0,TObject::kOverwrite) ;
 439   sr1->Write(0,TObject::kOverwrite) ;
 440   ar1->Write(0,TObject::kOverwrite) ;
 441   br1->Write(0,TObject::kOverwrite) ;
 442   nr1->Write(0,TObject::kOverwrite) ;
 443   nr1int->Write(0,TObject::kOverwrite) ;
 444   ar2->Write(0,TObject::kOverwrite) ;
 445   br2->Write(0,TObject::kOverwrite) ;
 446   cr2->Write(0,TObject::kOverwrite) ;
 447   mr2->Write(0,TObject::kOverwrite) ;
 448   sr2->Write(0,TObject::kOverwrite) ;
 449   nr2->Write(0,TObject::kOverwrite) ;
 450   nr2int->Write(0,TObject::kOverwrite) ;
 451
 452
 453   outFile->Close();
 454   delete outFile;
 455
 456   file->Close();
 457   delete file;
 458
 459   delete ratioCanvas;
 460   delete sbsCanvas;
 461   delete rawCanvas;
 462 }
 463
 464 //-----------------------------------------------------------------------------
 465 Double_t PeakPosition(Double_t pt){
 466   //Fit to the measured peak position
 467   return 4.99292e-003*exp(-pt*9.32300e-001)+1.34944e-001 ;
 468 }
 469 //-----------------------------------------------------------------------------
 470 Double_t PeakWidth(Double_t pt){
 471   //fit to the measured peak width
 472   Double_t a=0.0068 ;
 473   Double_t b=0.0025 ;
 474   Double_t c=0.000319 ;
 475   return TMath::Sqrt(a*a+b*b/pt/pt+c*c*pt*pt) ;
 476 }
 477
 478 //-----------------------------------------------------------------------------
 479 Double_t CB(Double_t * x, Double_t * par){
 480   //Parameterization of Real/Mixed ratio
 481   Double_t m=par[1] ;
 482   Double_t s=par[2] ;
 483   Double_t dx=(x[0]-m)/s ;
 484   return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean) ;
 485 }
 486
 487 //-----------------------------------------------------------------------------
 488 Double_t CB2(Double_t * x, Double_t * par){
 489   //Another parameterization of Real/Mixed ratio7TeV/README
 490   Double_t m=par[1] ;
 491   Double_t s=par[2] ;
 492   Double_t dx=(x[0]-m)/s ;
 493   return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean)+par[5]*(x[0]-kMean)*(x[0]-kMean) ;
 494 }
 495 //-----------------------------------------------------------------------------
 496 Double_t CBs(Double_t * x, Double_t * par){
 497   //Parameterizatin of signal
 498   Double_t m=par[1] ;
 499   Double_t s=par[2] ;
 500   Double_t dx=(x[0]-m)/s ;
 501   return par[0]*exp(-dx*dx/2.)/TMath::Sqrt(TMath::TwoPi())/s+par[3] ;
 502 }
 503 //-----------------------------------------------------------------------------
 504 Double_t BG1(Double_t * x, Double_t * par){
 505   //Normalizatino of Mixed
 506   return par[0]+par[1]*(x[0]-kMean) ;
 507 }
 508 //-----------------------------------------------------------------------------
 509 Double_t BG2(Double_t * x, Double_t * par){
 510   //Another normalization of  Mixed
 511   return par[0]+par[1]*(x[0]-kMean)+par[2]*(x[0]-kMean)*(x[0]-kMean) ;
 512 }
 513
 514
 515 //-----------------------------------------------------------------------------
 516 void PPRstyle()
 517 {
 518
 519   //////////////////////////////////////////////////////////////////////
 520   //
 521   // ROOT style macro for the TRD TDR
 522   //
 523   //////////////////////////////////////////////////////////////////////
 524
 525   gStyle->SetPalette(1);
 526   gStyle->SetCanvasBorderMode(-1);
 527   gStyle->SetCanvasBorderSize(1);
 528   gStyle->SetCanvasColor(10);
 529
 530   gStyle->SetFrameFillColor(10);
 531   gStyle->SetFrameBorderSize(1);
 532   gStyle->SetFrameBorderMode(-1);
 533   gStyle->SetFrameLineWidth(1.2);
 534   gStyle->SetFrameLineColor(1);
 535
 536   gStyle->SetHistFillColor(0);
 537   gStyle->SetHistLineWidth(1);
 538   gStyle->SetHistLineColor(1);
 539
 540   gStyle->SetPadColor(10);
 541   gStyle->SetPadBorderSize(1);
 542   gStyle->SetPadBorderMode(-1);
 543
 544   gStyle->SetStatColor(10);
 545   gStyle->SetTitleColor(kBlack,"X");
 546   gStyle->SetTitleColor(kBlack,"Y");
 547
 548   gStyle->SetLabelSize(0.04,"X");
 549   gStyle->SetLabelSize(0.04,"Y");
 550   gStyle->SetLabelSize(0.04,"Z");
 551   gStyle->SetTitleSize(0.04,"X");
 552   gStyle->SetTitleSize(0.04,"Y");
 553   gStyle->SetTitleSize(0.04,"Z");
 554   gStyle->SetTitleFont(42,"X");
 555   gStyle->SetTitleFont(42,"Y");
 556   gStyle->SetTitleFont(42,"X");
 557   gStyle->SetLabelFont(42,"X");
 558   gStyle->SetLabelFont(42,"Y");
 559   gStyle->SetLabelFont(42,"Z");
 560   gStyle->SetStatFont(42);
 561
 562   gStyle->SetTitleOffset(1.0,"X");
 563   gStyle->SetTitleOffset(1.4,"Y");
 564
 565   gStyle->SetFillColor(kWhite);
 566   gStyle->SetTitleFillColor(kWhite);
 567
 568   gStyle->SetOptDate(0);
 569   gStyle->SetOptTitle(1);
 570   gStyle->SetOptStat(0);
 571   gStyle->SetOptFit(111);
 572
 573 }
 574
 575 // For different tasks (e.g. triggers)
 576 void MakeMmixPi0()
 577 {
 578   // Take care when uncommenting, seems to leak memory consumption.
 579
 580   // // All Centrality
 581   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow/PHOSPi0FlowCoutput1", -1, "CPV", "out/kCentral/");
 582   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_SemiCentral/PHOSPi0Flow_SemiCentralCoutput1", -1, "CPV", "out/kSemiCentral/");
 583   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_MB/PHOSPi0Flow_MBCoutput1", -1, "CPV", "out/kMB/");
 584   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_PHOSPb/PHOSPi0Flow_PHOSPbCoutput1", -1, "CPV", "out/kPHOSPb/");
 585
 586   // 0-10%
 587   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow/PHOSPi0FlowCoutput1", 0, "CPV", "out/kCentral/");
 588   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_MB/PHOSPi0Flow_MBCoutput1", 0, "CPV", "out/kMB/");
 589   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_PHOSPb/PHOSPi0Flow_PHOSPbCoutput1", 0, "CPV", "out/kPHOSPb/");
 590
 591   // // 10-40%
 592   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_SemiCentral/PHOSPi0Flow_SemiCentralCoutput1", 1, "CPV", "out/kSemiCentral/");
 593   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_MB/PHOSPi0Flow_MBCoutput1", 1, "CPV", "out/kMB/");
 594   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_PHOSPb/PHOSPi0Flow_PHOSPbCoutput1", 1, "CPV", "out/kPHOSPb/");
 595
 596   // // 40-80%
 597   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_SemiCentral/PHOSPi0Flow_SemiCentralCoutput1", 2, "CPV", "out/kSemiCentral/");
 598   // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_MB/PHOSPi0Flow_MBCoutput1", 2, "CPV", "out/kMB/");
 599   MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_PHOSPb/PHOSPi0Flow_PHOSPbCoutput1", 2, "CPV", "out/kPHOSPb/");
 600 }
 601
 602
 603
 604 TH2F* FindPi0(TList *histoList, bool mix, const char* pid, int centrality)
 605 {
 606   // consider mixed or not
 607   TString mixStr("");
 608   if( mix )
 609     mixStr = "Mi";
 610
 611   // If centrality is some integer in range [0, kNCentralityBins] return Pi0 plot for that cent.
 612   if( centrality >= 0 && centrality < kNCentralityBins ) {
 613     TString name = Form("h%sPi0%s_cen%d", mixStr.Data(), pid, centrality);
 614     TH2F* hist = (TH2F*) histoList->FindObject(name.Data());
 615     hist->Sumw2();
 616     return hist;
 617   }
 618   // If centrality is '-1' Merge [0, kNCentralityBins)
 619   else if ( centrality == -1 ) {
 620     TString name = Form("h%sPi0%s_cen%d", mixStr.Data(), pid, 0);
 621     TH2F* histMerge = (TH2F*) histoList->FindObject(name.Data()) -> Clone(Form("h%sPi0%s_merged", mixStr.Data(), pid));
 622     histMerge->Sumw2();
 623
 624     for( int cent = 1; cent < kNCentralityBins; ++cent ) {
 625       name = Form("h%sPi0%s_cen%d", mixStr.Data(), pid, cent);
 626       TH2F* other = (TH2F*) histoList->FindObject(name.Data());
 627       other->Sumw2();
 628       histMerge->Add( other );
 629     }
 630     return histMerge;
 631   }
 632   else { // case not defined
 633     Printf("ERROR: Centrality must be in range: [%d,%d]", -1, kNCentralityBins - 1 );
 634     return 0;
 635   }
 636 }