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