PWG0/multPbPb/correct.C

   1 #include "AliAnalysisMultPbTrackHistoManager.h"
   2 #include "TH2D.h"
   3 #include "TH1D.h"
   4 #include "TF1.h"
   5 #include "TCanvas.h"
   6 #include "TFile.h"
   7 #include "TSystem.h"
   8 #include "TROOT.h"
   9 #include <iostream>
  10 #include "TDatabasePDG.h"
  11 #include "AliPhysicsSelection.h"
  12 #include "AliESDtrackCuts.h"
  13
  14 using namespace std;
  15
  16 AliAnalysisMultPbTrackHistoManager * hManData = 0;
  17 AliAnalysisMultPbTrackHistoManager * hManCorr = 0;
  18 TH2D * hEvStatData = 0;
  19 TH2D * hEvStatCorr = 0;
  20
  21 const Int_t kHistoFitCompoments = 2;
  22 TH1D * gHistoCompoments[kHistoFitCompoments];
  23
  24 void LoadLibs();
  25 void LoadData(TString dataFolder, TString correctionFolder);
  26 void SetStyle();
  27 Double_t CheckSecondaries();
  28 void ShowAcceptanceInVzSlices() ;
  29 TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) ;
  30   TH1D * GetCumulativeHisto (TH1 * h) ;
  31 static Double_t HistoSum(const double * x, const double* p);
  32 TF1 * GetFunctionHistoSum() ;
  33 TF1 * GetMTExp(Float_t norm=68, Float_t t=25) ;
  34 TF1 * GetHagedorn(Float_t norm=68, Float_t pt0=25, Float_t n=13) ;
  35 TF1 * GetLevy(Double_t temp=0.1, Double_t n=7, Double_t norm=10, const char * name="fLevy") ;
  36
  37
  38
  39
  40 void correct(TString dataFolder = "output/LHC09d_000104892_p4/", TString correctionFolder = "output/LHC10a8_104867/") {
  41
  42   // Load stuff and set some styles
  43   LoadLibs();
  44   LoadData(dataFolder,correctionFolder);
  45   SetStyle();
  46   // ShowAcceptanceInVzSlices();
  47   // return;
  48
  49   // TODO add some cool printout for cuts and centrality selection
  50
  51   Double_t fractionWeak = CheckSecondaries();
  52   cout << "Rescaling weak correction: " << fractionWeak << endl;
  53
  54
  55   // Some shorthands
  56   // FIXME: Gen should be projected including overflow in z?
  57   TH1D * hDataPt   = (TH1D*) hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec,        -0.5,0.5,-10,10)->Clone("hDataPt");
  58   TH1D * hMCPtGen  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        -0.5,0.5,-10,10);
  59   TH1D * hMCPtRec  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec,        -0.5,0.5,-10,10);
  60   TH1D * hMCPtPri  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim,    -0.5,0.5,-10,10);
  61   TH1D * hMCPtSeM  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat,  -0.5,0.5,-10,10);
  62   TH1D * hMCPtSeW  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak, -0.5,0.5,-10,10);
  63   TH1D * hMCPtFak  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecFake,    -0.5,0.5,-10,10);
  64
  65   TCanvas * cdata = new TCanvas ("cData", "Data");
  66
  67   hDataPt->Draw();
  68   //  hMCPtRec->Draw("same");
  69
  70   TCanvas * cMC = new TCanvas ("cMC", "Monte Carlo");
  71   cMC->cd();
  72   hMCPtGen ->Draw();
  73   hMCPtRec ->Draw("same");
  74   hMCPtPri ->Draw("same");
  75   hMCPtSeM ->Draw("same");
  76   hMCPtSeW ->Draw("same");
  77   hMCPtFak ->Draw("same");
  78
  79   cout << "Fake/All Rec  = " << hMCPtFak->Integral()/hMCPtRec->Integral()  << endl;
  80   cout << "SM/All   Rec  = " << hMCPtSeM->Integral()/hMCPtRec->Integral()  << endl;
  81   cout << "SW/All   Rec  = " << hMCPtSeW->Integral()/hMCPtRec->Integral()  << endl;
  82
  83
  84   // Compute efficiency and subtract secondaries and fakes after rescaling to data
  85   // PRIM_DATA = ALL_DATA - SEC_MC/ALL_MC*ALL_DATA - FAK_MC/ALL_MC*ALL_DATA
  86   // TRUE_DATA = PRIM_DATA * GEN_MC/PRIM_MC
  87
  88   TH1D * hEffPt = (TH1D*) hMCPtPri->Clone("hEffPt");
  89   hEffPt->Divide(hMCPtPri,hMCPtGen,1,1,"B");
  90
  91   TH1D * hCorSeM = (TH1D*) hMCPtSeM->Clone("hEffPt");
  92   hCorSeM->Divide(hMCPtSeM,hMCPtRec,1,1,"B");
  93   hCorSeM->Multiply(hDataPt);
  94
  95   TH1D * hCorSeW = (TH1D*) hMCPtSeW->Clone("hEffPt");
  96   hCorSeW->Divide(hMCPtSeW,hMCPtRec,1,1,"B");
  97   hCorSeW->Scale(fractionWeak);// rescale weak correction
  98   hCorSeW->Multiply(hDataPt);
  99
 100   TH1D * hCorFak = (TH1D*) hMCPtFak->Clone("hEffPt");
 101   hCorFak->Divide(hMCPtFak,hMCPtRec,1,1,"B");
 102   hCorFak->Multiply(hDataPt);
 103
 104   TH1D * hCorrected = (TH1D*) hDataPt->Clone("hCorrected");
 105   hCorrected->Add(hCorSeM,-1);
 106   hCorrected->Add(hCorSeW,-1);
 107   hCorrected->Add(hCorFak,-1);
 108   hCorrected->Divide(hEffPt);
 109   hCorrected->SetMarkerStyle(kOpenStar);
 110
 111   cdata->cd();
 112   hDataPt->Draw();
 113   hCorrected->SetLineColor(kBlack);
 114   hCorrected->SetMarkerColor(kBlack);
 115   hCorrected->Draw("same");
 116   hMCPtGen->DrawClone("same");
 117   TF1 * f = GetLevy();
 118   hCorrected->Fit(f,"", "same");
 119   hCorrected->Fit(f,"IME", "same");
 120   cout << "dN/deta = " << f->Integral(0,100) << " +- " << f->IntegralError(0,100) << endl;
 121   cout << "Generated dN/deta (correction) = " << hMCPtGen->Integral() << endl;
 122   // FIXME: comment this out
 123   TH1D * hDataGen  = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        -0.5,0.5,-10,10);
 124   cout << "Generated dN/deta (data) =       " << hDataGen->Integral() << endl;
 125   hDataGen->Draw("same");
 126 }
 127
 128 Double_t CheckSecondaries() {
 129   // Returns the fraction you need to rescale the secondaries from weak decays for.
 130
 131   // Some shorthands
 132   TH1D * hDataDCA   = hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 133   //  TH1D * hMCDCAGen  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
 134   TH1D * hMCDCARec  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 135   TH1D * hMCDCAPri  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
 136   TH1D * hMCDCASW  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak    );
 137   TH1D * hMCDCASM  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat    );
 138   TH1D * hMCDCAFak  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
 139
 140
 141   TCanvas * cCumulative = new TCanvas("cDCAculumative","DCA cumulative distributions");
 142   cCumulative->cd();
 143   GetCumulativeHisto(hMCDCAPri )->Draw();
 144   GetRatioIntegratedFractions(hMCDCASW, hMCDCARec  )->Draw("same");
 145   GetRatioIntegratedFractions(hMCDCASM, hMCDCARec  )->Draw("same");
 146   GetRatioIntegratedFractions(hMCDCAPri,hMCDCARec  )->Draw("same");
 147
 148
 149   TCanvas * c1 = new TCanvas("cDCAFit", "Fit to the DCA distributions");
 150   c1->SetLogy();
 151   // Draw all
 152   //  hDataDCA->Draw();
 153   //  hMCDCARec ->Draw("same");
 154   // hMCDCAPri ->Draw("same");
 155   // hMCDCASW ->Draw("same");
 156   // hMCDCASM ->Draw("same");
 157   // hMCDCAFak ->Draw("same");
 158   //  return 1;
 159
 160   TH1D * hMCPrimSMFak = (TH1D*) hMCDCAPri->Clone("hMCPrimSMFak");
 161   hMCPrimSMFak->Add(hMCDCASM);
 162   hMCPrimSMFak->Add(hMCDCAFak);
 163
 164   gHistoCompoments[0] = (TH1D*) hMCPrimSMFak->Clone();
 165   gHistoCompoments[1] = (TH1D*) hMCDCASW->Clone();
 166   TF1 * fHistos = GetFunctionHistoSum();
 167   fHistos->SetParameters(1,1);
 168   hDataDCA->Fit(fHistos,"","",0,50);
 169   hMCPrimSMFak->Scale(fHistos->GetParameter(0));
 170   hMCDCASW    ->Scale(fHistos->GetParameter(1));
 171   hMCPrimSMFak->Draw("same");
 172   hMCDCASW    ->Draw("same");
 173   return fHistos->GetParameter(1)/fHistos->GetParameter(0);
 174
 175
 176 }
 177
 178 void LoadLibs() {
 179
 180   gSystem->Load("libVMC");
 181   gSystem->Load("libTree");
 182   gSystem->Load("libSTEERBase");
 183   gSystem->Load("libESD");
 184   gSystem->Load("libAOD");
 185   gSystem->Load("libANALYSIS");
 186   gSystem->Load("libANALYSISalice");
 187   gSystem->Load("libCORRFW");
 188   gSystem->Load("libMinuit");
 189   gSystem->Load("libPWG2Spectra");
 190   gSystem->Load("libPWG0base");
 191
 192   gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG0"));
 193   gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG1/background"));
 194   // Load helper classes
 195   // TODO: replace this by a list of TOBJStrings
 196   TString taskName("AliAnalysisTaskMultPbTracks.cxx+");
 197   TString histoManName("AliAnalysisMultPbTrackHistoManager.cxx+");
 198   TString centrName("AliAnalysisMultPbCentralitySelector.cxx+");
 199   TString listName("$ALICE_ROOT/PWG1/background/AliHistoListWrapper.cxx+");
 200
 201   Bool_t debug=0;
 202   gROOT->LoadMacro(listName    +(debug?"+g":""));
 203   gROOT->LoadMacro(histoManName+(debug?"+g":""));
 204   gROOT->LoadMacro(centrName   +(debug?"+g":""));
 205   gROOT->LoadMacro(taskName    +(debug?"+g":""));
 206
 207   // Histo fitter
 208   gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/fcn.cxx+g");
 209   gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/AliHistoFitter.cxx+g");
 210
 211
 212 }
 213
 214
 215 void SetStyle() {
 216
 217   hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kBlack);
 218   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetLineColor(kRed  );
 219   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kRed  );
 220   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
 221   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
 222   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue );
 223   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
 224
 225   hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kBlack);
 226   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerColor(kRed  );
 227   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kRed  );
 228   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
 229   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
 230   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue );
 231   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
 232
 233   hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kFullCircle);
 234   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerStyle(kFullSquare);
 235   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kOpenSquare);
 236   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
 237   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
 238   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenCircle);
 239   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
 240
 241  hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kBlack);
 242   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetLineColor(kRed  );
 243   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kRed  );
 244   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
 245   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
 246   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue );
 247   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
 248
 249   hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kBlack);
 250   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerColor(kRed  );
 251   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kRed  );
 252   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
 253   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
 254   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue );
 255   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
 256
 257   hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kFullCircle);
 258   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerStyle(kFullSquare);
 259   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kOpenSquare);
 260   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
 261   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
 262   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenSquare);
 263   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
 264
 265
 266 }
 267
 268 void LoadData(TString dataFolder, TString correctionFolder){
 269
 270   // Get histo manager for data and MC + stat histos
 271   TFile * fData = new TFile(dataFolder+"multPbPbtracks.root");
 272   TFile * fCorr = new TFile(correctionFolder+"multPbPbtracks.root");
 273   TFile * fStatData = new TFile(dataFolder+"event_stat.root");
 274   TFile * fStatCorr = new TFile(correctionFolder+"event_stat.root");
 275
 276   hManData = (AliAnalysisMultPbTrackHistoManager*) fData->Get("histoManager");
 277   hManCorr = (AliAnalysisMultPbTrackHistoManager*) fCorr->Get("histoManager");
 278   AliESDtrackCuts * cutsData = (AliESDtrackCuts*) fData->Get("AliESDtrackCuts");
 279   AliESDtrackCuts * cutsCorr = (AliESDtrackCuts*) fCorr->Get("AliESDtrackCuts");
 280   if (cutsData != cutsCorr) {
 281     cout << "ERROR: MC and data do not have the same cuts" << endl;
 282     // FIXME: exit here
 283   }
 284   cutsData->Print();
 285   hEvStatData = (TH2D*) fStatData->Get("fHistStatistics");
 286   hEvStatCorr = (TH2D*) fStatCorr->Get("fHistStatistics");
 287
 288   // Normalize
 289   Int_t irowGoodTrigger = 1;
 290   if (hEvStatCorr && hEvStatData) {
 291     //  hManData->ScaleHistos(75351.36/1.015);// Nint for run 104892 estimated correcting for the trigger efficiency, multiplied for the physics selection efficiency which I'm not correcting for the time being
 292     hManData->ScaleHistos(hEvStatData->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
 293     hManCorr->ScaleHistos(hEvStatCorr->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
 294   } else {
 295     cout << "WARNING!!! ARBITRARY SCALING" << endl;
 296     hManData->ScaleHistos(1000);
 297     hManCorr->ScaleHistos(1000);
 298   }
 299 }
 300
 301 TF1 * GetHagedorn(Float_t norm, Float_t pt0, Float_t n) {
 302
 303   TF1 * f =0;
 304   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 305
 306   f=new TF1("fHagedorn",Form("(x/sqrt(x*x+%f*%f))*x*[0]*( 1 + x/[1] )^(-[2])",mass,mass),0,10);
 307   f->SetParameters(norm, pt0, n);
 308   f->SetParLimits(1, 0.01, 10);
 309   f->SetParNames("norm", "pt0", "n");
 310   f->SetLineWidth(1);
 311   return f;
 312
 313
 314 }
 315
 316 TF1 * GetMTExp(Float_t norm, Float_t t) {
 317
 318   TF1 * f =0;
 319   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 320
 321   f=new TF1("fMTExp",Form("(x/sqrt(x*x+%f*%f))*x*[0]*exp(-sqrt(x*x+%f*%f)/[1])",mass,mass,mass,mass),0,10);
 322   f->SetParameters(norm, t);
 323   //  f->SetParLimits(1, 0.01);
 324   f->SetParNames("norm", "T");
 325   f->SetLineWidth(1);
 326   return f;
 327
 328
 329 }
 330
 331 TF1 * GetLevy(Double_t temp, Double_t n, Double_t norm, const char * name) {
 332
 333   char formula[500];
 334   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 335
 336
 337   sprintf(formula,"(x/sqrt(x*x+[3]*[3]))*x*( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2])  )^(-[1])");
 338   TF1* f=new TF1(name,formula,0,10);
 339   f->SetParameters(norm, n, temp,mass);
 340   f->SetParLimits(2, 0.01, 10);
 341   f->SetParNames("norm (dN/dy)", "n", "T", "mass");
 342   f->FixParameter(3,mass);
 343   f->SetLineWidth(1);
 344   return f;
 345
 346
 347 }
 348
 349
 350 TH1D * GetCumulativeHisto (TH1 * h) {
 351   // Returns a cumulative histogram
 352   // FIXME: put this method in a tools class
 353   TH1D * hInt = (TH1D*) h->Clone(TString(h->GetName())+"_Int");
 354   hInt->Reset();
 355   Float_t integral = h->Integral(-1,-1); // Consider under/over flow!
 356   Int_t nbin = h->GetNbinsX();
 357   for(Int_t ibin = 1; ibin <= nbin; ibin++){
 358     Double_t content = h->Integral(-1,ibin) / integral;
 359     hInt->SetBinContent(ibin, content);
 360   }
 361   return hInt;
 362 }
 363
 364 TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) {
 365   // Returns the the ratio of integrated histograms
 366   // FIXME: put this method in a tools class
 367   TH1D * hRatio = (TH1D*) hNum->Clone(TString(hNum->GetName())+hDenum->GetName()+"_RatioIntegrated");
 368   hRatio->Reset();
 369   Int_t nbin = hNum->GetNbinsX();
 370   for(Int_t ibin = 1; ibin <= nbin; ibin++){
 371     Double_t content = hNum->Integral(-1,ibin) / hDenum->Integral(-1,ibin);// consider underflow
 372     hRatio->SetBinContent(ibin, content);
 373   }
 374   return hRatio;
 375 }
 376
 377 void ShowAcceptanceInVzSlices() {
 378   TCanvas * cvz = new TCanvas("cvz","acc #times eff vs vz");
 379   for(Int_t ivz = -10; ivz < 10; ivz+=4){
 380     TH1D * hMCPtPri  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim ,   -0.5,0.5,ivz,ivz+4);
 381     TH1D * hMCPtGen  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        -0.5,0.5,ivz,ivz+4);
 382     //    hEff= hMCPtGen;
 383     TH1D * hEff = (TH1D*)hMCPtPri->Clone(Form("hEff_vz_%d_%d",ivz,ivz+4));
 384     hEff->Divide(hMCPtPri,hMCPtGen,1,1,"B");
 385     cout << "ivz " << ivz << endl;
 386     if(ivz < -9) {
 387       cout << "First" << endl;
 388       hEff->Draw();
 389       // hMCPtGen->Draw();
 390       // hMCPtPri->Draw("same");
 391     }
 392     else {
 393       cout << "Same" << endl;
 394       hEff->Draw("same");
 395       // hMCPtGen->Draw("");
 396       // hMCPtPri->Draw("same");
 397     }
 398     cvz->Update();
 399     //    cvz->WaitPrimitive();
 400   }
 401
 402   //hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim )->Draw();
 403    hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec )->Draw("");
 404   hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoGen     )->Draw("same");
 405
 406
 407 }
 408
 409 TF1 * GetFunctionHistoSum() {
 410
 411   TF1 * f = new TF1 ("fHistoSum",HistoSum,0,50,kHistoFitCompoments);
 412   f->SetParNames("All", "Sec. Weak decays");
 413   f->SetNpx(1000);
 414   return f;
 415
 416 }
 417
 418 static Double_t HistoSum(const double * x, const double* p){
 419   // This function uses a global array of histograms, rescaled by some
 420   // parameters, to define a function
 421   // The array is called gHistoCompoments
 422   // The size of the arrays is given by the global constant kHistoFitCompoments
 423
 424   Double_t xx = x[0];
 425   Double_t value = 0;
 426   for(Int_t icomp = 0; icomp < kHistoFitCompoments; icomp++){
 427     //    value += gHistoCompoments[icomp]-Interpolate(xx) * p[icomp];
 428     Int_t ibin = gHistoCompoments[icomp]->FindBin(xx);
 429     value += gHistoCompoments[icomp]->GetBinContent(ibin) * p[icomp];
 430   }
 431
 432   return value;
 433
 434 }