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