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   // //  hMCDCAGen ->Draw("same");
 154   //  hMCDCARec ->Draw("same");
 155   // hMCDCAPri ->Draw("same");
 156   // hMCDCASW ->Draw("same");
 157   // hMCDCASM ->Draw("same");
 158   // hMCDCAFak ->Draw("same");
 159   //  return;
 160
 161   TH1D * hMCPrimSMFak = (TH1D*) hMCDCAPri->Clone("hMCPrimSMFak");
 162   hMCPrimSMFak->Add(hMCDCASM);
 163   hMCPrimSMFak->Add(hMCDCAFak);
 164
 165   gHistoCompoments[0] = (TH1D*) hMCPrimSMFak->Clone();
 166   gHistoCompoments[1] = (TH1D*) hMCDCASW->Clone();
 167   TF1 * fHistos = GetFunctionHistoSum();
 168   fHistos->SetParameters(1,1);
 169   hDataDCA->Fit(fHistos,"","",0,50);
 170   hMCPrimSMFak->Scale(fHistos->GetParameter(0));
 171   hMCDCASW    ->Scale(fHistos->GetParameter(1));
 172   hMCPrimSMFak->Draw("same");
 173   hMCDCASW    ->Draw("same");
 174   return fHistos->GetParameter(1)/fHistos->GetParameter(0);
 175
 176
 177 }
 178
 179 void LoadLibs() {
 180
 181   gSystem->Load("libVMC");
 182   gSystem->Load("libTree");
 183   gSystem->Load("libSTEERBase");
 184   gSystem->Load("libESD");
 185   gSystem->Load("libAOD");
 186   gSystem->Load("libANALYSIS");
 187   gSystem->Load("libANALYSISalice");
 188   gSystem->Load("libCORRFW");
 189   gSystem->Load("libMinuit");
 190   gSystem->Load("libPWG2Spectra");
 191   gSystem->Load("libPWG0base");
 192
 193   gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG0"));
 194   gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG1/background"));
 195   // Load helper classes
 196   // TODO: replace this by a list of TOBJStrings
 197   TString taskName("AliAnalysisTaskMultPbTracks.cxx+");
 198   TString histoManName("AliAnalysisMultPbTrackHistoManager.cxx+");
 199   TString centrName("AliAnalysisMultPbCentralitySelector.cxx+");
 200   TString listName("$ALICE_ROOT/PWG1/background/AliHistoListWrapper.cxx+");
 201
 202   Bool_t debug=0;
 203   gROOT->LoadMacro(listName    +(debug?"+g":""));
 204   gROOT->LoadMacro(histoManName+(debug?"+g":""));
 205   gROOT->LoadMacro(centrName   +(debug?"+g":""));
 206   gROOT->LoadMacro(taskName    +(debug?"+g":""));
 207
 208   // Histo fitter
 209   gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/fcn.cxx+g");
 210   gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/AliHistoFitter.cxx+g");
 211
 212
 213 }
 214
 215
 216 void SetStyle() {
 217
 218   hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kBlack);
 219   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetLineColor(kRed  );
 220   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kRed  );
 221   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
 222   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
 223   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue );
 224   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
 225
 226   hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kBlack);
 227   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerColor(kRed  );
 228   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kRed  );
 229   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
 230   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
 231   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue );
 232   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
 233
 234   hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kFullCircle);
 235   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerStyle(kFullSquare);
 236   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kOpenSquare);
 237   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
 238   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
 239   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenCircle);
 240   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
 241
 242  hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kBlack);
 243   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetLineColor(kRed  );
 244   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kRed  );
 245   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
 246   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
 247   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue );
 248   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
 249
 250   hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kBlack);
 251   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerColor(kRed  );
 252   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kRed  );
 253   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
 254   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
 255   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue );
 256   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
 257
 258   hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kFullCircle);
 259   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerStyle(kFullSquare);
 260   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kOpenSquare);
 261   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
 262   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
 263   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenSquare);
 264   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
 265
 266
 267 }
 268
 269 void LoadData(TString dataFolder, TString correctionFolder){
 270
 271   // Get histo manager for data and MC + stat histos
 272   TFile * fData = new TFile(dataFolder+"multPbPbtracks.root");
 273   TFile * fCorr = new TFile(correctionFolder+"multPbPbtracks.root");
 274   TFile * fStatData = new TFile(dataFolder+"event_stat.root");
 275   TFile * fStatCorr = new TFile(correctionFolder+"event_stat.root");
 276
 277   hManData = (AliAnalysisMultPbTrackHistoManager*) fData->Get("histoManager");
 278   hManCorr = (AliAnalysisMultPbTrackHistoManager*) fCorr->Get("histoManager");
 279   AliESDtrackCuts * cutsData = (AliESDtrackCuts*) fData->Get("AliESDtrackCuts");
 280   AliESDtrackCuts * cutsCorr = (AliESDtrackCuts*) fCorr->Get("AliESDtrackCuts");
 281   if (cutsData != cutsCorr) {
 282     cout << "ERROR: MC and data do not have the same cuts" << endl;
 283     // FIXME: exit here
 284   }
 285   cutsData->Print();
 286   hEvStatData = (TH2D*) fStatData->Get("fHistStatistics");
 287   hEvStatCorr = (TH2D*) fStatCorr->Get("fHistStatistics");
 288
 289   // Normalize
 290   Int_t irowGoodTrigger = 1;
 291   if (hEvStatCorr && hEvStatData) {
 292     //  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
 293     hManData->ScaleHistos(hEvStatData->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
 294     hManCorr->ScaleHistos(hEvStatCorr->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
 295   } else {
 296     cout << "WARNING!!! ARBITRARY SCALING" << endl;
 297     hManData->ScaleHistos(1000);
 298     hManCorr->ScaleHistos(1000);
 299   }
 300 }
 301
 302 TF1 * GetHagedorn(Float_t norm, Float_t pt0, Float_t n) {
 303
 304   TF1 * f =0;
 305   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 306
 307   f=new TF1("fHagedorn",Form("(x/sqrt(x*x+%f*%f))*x*[0]*( 1 + x/[1] )^(-[2])",mass,mass),0,10);
 308   f->SetParameters(norm, pt0, n);
 309   f->SetParLimits(1, 0.01, 10);
 310   f->SetParNames("norm", "pt0", "n");
 311   f->SetLineWidth(1);
 312   return f;
 313
 314
 315 }
 316
 317 TF1 * GetMTExp(Float_t norm, Float_t t) {
 318
 319   TF1 * f =0;
 320   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 321
 322   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);
 323   f->SetParameters(norm, t);
 324   //  f->SetParLimits(1, 0.01);
 325   f->SetParNames("norm", "T");
 326   f->SetLineWidth(1);
 327   return f;
 328
 329
 330 }
 331
 332 TF1 * GetLevy(Double_t temp, Double_t n, Double_t norm, const char * name) {
 333
 334   char formula[500];
 335   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 336
 337
 338   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])");
 339   TF1* f=new TF1(name,formula,0,10);
 340   f->SetParameters(norm, n, temp,mass);
 341   f->SetParLimits(2, 0.01, 10);
 342   f->SetParNames("norm (dN/dy)", "n", "T", "mass");
 343   f->FixParameter(3,mass);
 344   f->SetLineWidth(1);
 345   return f;
 346
 347
 348 }
 349
 350
 351 TH1D * GetCumulativeHisto (TH1 * h) {
 352   // Returns a cumulative histogram
 353   // FIXME: put this method in a tools class
 354   TH1D * hInt = (TH1D*) h->Clone(TString(h->GetName())+"_Int");
 355   hInt->Reset();
 356   Float_t integral = h->Integral(-1,-1); // Consider under/over flow!
 357   Int_t nbin = h->GetNbinsX();
 358   for(Int_t ibin = 1; ibin <= nbin; ibin++){
 359     Double_t content = h->Integral(-1,ibin) / integral;
 360     hInt->SetBinContent(ibin, content);
 361   }
 362   return hInt;
 363 }
 364
 365 TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) {
 366   // Returns the the ratio of integrated histograms
 367   // FIXME: put this method in a tools class
 368   TH1D * hRatio = (TH1D*) hNum->Clone(TString(hNum->GetName())+hDenum->GetName()+"_RatioIntegrated");
 369   hRatio->Reset();
 370   Int_t nbin = hNum->GetNbinsX();
 371   for(Int_t ibin = 1; ibin <= nbin; ibin++){
 372     Double_t content = hNum->Integral(-1,ibin) / hDenum->Integral(-1,ibin);// consider underflow
 373     hRatio->SetBinContent(ibin, content);
 374   }
 375   return hRatio;
 376 }
 377
 378 void ShowAcceptanceInVzSlices() {
 379   TCanvas * cvz = new TCanvas("cvz","acc #times eff vs vz");
 380   for(Int_t ivz = -10; ivz < 10; ivz+=4){
 381     TH1D * hMCPtPri  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim ,   -0.5,0.5,ivz,ivz+4);
 382     TH1D * hMCPtGen  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        -0.5,0.5,ivz,ivz+4);
 383     //    hEff= hMCPtGen;
 384     TH1D * hEff = (TH1D*)hMCPtPri->Clone(Form("hEff_vz_%d_%d",ivz,ivz+4));
 385     hEff->Divide(hMCPtPri,hMCPtGen,1,1,"B");
 386     cout << "ivz " << ivz << endl;
 387     if(ivz < -9) {
 388       cout << "First" << endl;
 389       hEff->Draw();
 390       // hMCPtGen->Draw();
 391       // hMCPtPri->Draw("same");
 392     }
 393     else {
 394       cout << "Same" << endl;
 395       hEff->Draw("same");
 396       // hMCPtGen->Draw("");
 397       // hMCPtPri->Draw("same");
 398     }
 399     cvz->Update();
 400     //    cvz->WaitPrimitive();
 401   }
 402
 403   //hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim )->Draw();
 404    hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec )->Draw("");
 405   hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoGen     )->Draw("same");
 406
 407
 408 }
 409
 410 TF1 * GetFunctionHistoSum() {
 411
 412   TF1 * f = new TF1 ("fHistoSum",HistoSum,0,50,kHistoFitCompoments);
 413   f->SetParNames("All", "Sec. Weak decays");
 414   f->SetNpx(1000);
 415   return f;
 416
 417 }
 418
 419 static Double_t HistoSum(const double * x, const double* p){
 420   // This function uses a global array of histograms, rescaled by some
 421   // parameters, to define a function
 422   // The array is called gHistoCompoments
 423   // The size of the arrays is given by the global constant kHistoFitCompoments
 424
 425   Double_t xx = x[0];
 426   Double_t value = 0;
 427   for(Int_t icomp = 0; icomp < kHistoFitCompoments; icomp++){
 428     //    value += gHistoCompoments[icomp]-Interpolate(xx) * p[icomp];
 429     Int_t ibin = gHistoCompoments[icomp]->FindBin(xx);
 430     value += gHistoCompoments[icomp]->GetBinContent(ibin) * p[icomp];
 431   }
 432
 433   return value;
 434
 435 }