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 #include "TLegend.h"
  15
  16 using namespace std;
  17
  18 AliAnalysisMultPbTrackHistoManager * hManData = 0;
  19 AliAnalysisMultPbTrackHistoManager * hManCorr = 0;
  20 TH2D * hEvStatData = 0;
  21 TH2D * hEvStatCorr = 0;
  22
  23 const Int_t kHistoFitCompoments = 3;
  24 TH1D * gHistoCompoments[kHistoFitCompoments];
  25
  26 void LoadLibs(  Bool_t debug=0);
  27 void LoadData(TString dataFolder, TString correctionFolder);
  28 void SetStyle();
  29 void CheckSecondaries(Double_t & fracWeak, Double_t &fracMaterial);
  30 void CheckSanity();
  31 void ShowAcceptanceInVzSlices() ;
  32 TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) ;
  33 TH1D * GetCumulativeHisto (TH1 * h) ;
  34 static Double_t HistoSum(const double * x, const double* p);
  35 TF1 * GetFunctionHistoSum() ;
  36 TF1 * GetMTExp(Float_t norm=68, Float_t t=25) ;
  37 TF1 * GetHagedorn(Float_t norm=68, Float_t pt0=25, Float_t n=13) ;
  38 TF1 * GetLevy(Double_t temp=0.1, Double_t n=7, Double_t norm=10, const char * name="fLevy") ;
  39 void PrintCanvas(TCanvas* c,const TString formats) ;
  40
  41 // global switches
  42 Bool_t doPrint=kFALSE;// disable PrintCanvas
  43 Float_t zmin = -10;
  44 Float_t zmax = 10;
  45 Float_t etaMin = -0.5;
  46 Float_t etaMax = 0.5;
  47
  48 #define CORRECT_2D
  49
  50 void correct(TString dataFolder = "./output/LHC10g2d_130844_V0M_bin_10/", TString correctionFolder = "./output/LHC10g2a_130844_V0M_bin_10/",
  51              Float_t zminl=-10, Float_t zmaxl=10, Float_t etaMinl = -0.5, Float_t etaMaxl=0.5, Float_t npart=381.3) {
  52
  53   // Set vertex
  54   zmin = zminl;
  55   zmax = zmaxl;
  56   etaMin = etaMinl;
  57   etaMax = etaMaxl;
  58
  59   // Load stuff and set some styles
  60   LoadLibs();
  61   LoadData(dataFolder,correctionFolder);
  62   SetStyle();
  63   // ShowAcceptanceInVzSlices();
  64   // return;
  65
  66   // TODO add some cool printout for cuts and centrality selection
  67
  68   CheckSanity();
  69
  70   Double_t fractionWeak = 1, fractionMaterial=1;
  71   CheckSecondaries(fractionWeak, fractionMaterial);
  72   cout << "Rescaling secondaries correction, weak: " << fractionWeak << ", material: " << fractionMaterial <<endl;
  73
  74
  75   // Some shorthands
  76   // FIXME: Gen should be projected including overflow in z?
  77 #if defined CORRECT_1D
  78   TH1D * hDataPt   = (TH1D*) hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec, etaMin,etaMax,zmin,zmax)->Clone("hDataPt");
  79   TH1D * hMCPtGen  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,         etaMin,etaMax,zmin,zmax); //FIXME: che si fa qui?
  80   //zTH1D * hMCPtGen  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,         etaMin,etaMax,zmin,zmax);
  81   TH1D * hMCPtRec  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRec,         etaMin,etaMax,zmin,zmax);
  82   TH1D * hMCPtPri  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim,     etaMin,etaMax,zmin,zmax);
  83   TH1D * hMCPtSeM  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat,   etaMin,etaMax,zmin,zmax);
  84   TH1D * hMCPtSeW  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak,  etaMin,etaMax,zmin,zmax);
  85   TH1D * hMCPtFak  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecFake,     etaMin,etaMax,zmin,zmax);
  86 #elif defined CORRECT_2D
  87   TH1  * hDataPt   = (TH2D*) hManData->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRec,         etaMin,etaMax)->Clone("hDataPt");
  88   TH1  * hMCPtGen  = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoGen,         etaMin,etaMax);
  89   TH1  * hMCPtRec  = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRec,         etaMin,etaMax);
  90   TH1  * hMCPtPri  = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim,     etaMin,etaMax);
  91   TH1  * hMCPtSeM  = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat,   etaMin,etaMax);
  92   TH1  * hMCPtSeW  = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak,  etaMin,etaMax);
  93   TH1  * hMCPtFak  = (TH2D*) hManCorr->GetHistoPtVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake,     etaMin,etaMax);
  94 #endif
  95
  96   //  hMCPtRec->Draw("same");
  97
  98   TCanvas * cMC = new TCanvas ("cMC", "Monte Carlo");
  99   cMC->SetLogy();
 100   cMC->cd();
 101   // hMCPtSeW->Add(hMCPtSeM);
 102   // hMCPtSeW->Divide(hMCPtRec);
 103
 104
 105   hMCPtGen ->Draw();
 106   hMCPtRec ->Draw("same");
 107   hMCPtPri ->Draw("same");
 108   hMCPtSeM ->Draw("same");
 109   hMCPtSeW ->Draw("same");
 110   hMCPtFak ->Draw("same");
 111 #ifdef CORRECT_1D
 112   hMCPtGen ->GetXaxis()->SetRangeUser(0,4.5);
 113   hMCPtGen ->GetYaxis()->SetRangeUser(0.1,1e4);
 114 #endif
 115   TLegend * lMC = new TLegend(0.505034, 0.59965, 0.877517, 0.926573,"Monte Carlo");
 116   lMC->AddEntry( hMCPtGen, "Generated");
 117   lMC->AddEntry( hMCPtRec, "All Rec");
 118   lMC->AddEntry( hMCPtPri, "Rec Primaries");
 119   lMC->AddEntry( hMCPtSeM, "Rec Sec. Material");
 120   lMC->AddEntry( hMCPtSeW, "Rec Sec. Weak");
 121   lMC->AddEntry( hMCPtFak, "Rec Fakes");
 122   lMC->Draw();
 123
 124
 125   cout << "Fake/All Rec  = " << hMCPtFak->Integral()/hMCPtRec->Integral()  << endl;
 126   cout << "SM/All   Rec  = " << hMCPtSeM->Integral()/hMCPtRec->Integral()  << endl;
 127   cout << "SW/All   Rec  = " << hMCPtSeW->Integral()/hMCPtRec->Integral()  << endl;
 128
 129
 130   // Compute efficiency and subtract secondaries and fakes after rescaling to data
 131   // PRIM_DATA = ALL_DATA - SEC_MC/ALL_MC*ALL_DATA - FAK_MC/ALL_MC*ALL_DATA
 132   // TRUE_DATA = PRIM_DATA * GEN_MC/PRIM_MC
 133
 134
 135   TH1D * hEffPt = (TH1D*) hMCPtPri->Clone("hEffPt");
 136   hEffPt->Divide(hMCPtPri,hMCPtGen,1,1,"B");
 137
 138   TH1D * hCorSeM = (TH1D*) hMCPtSeM->Clone("hCorSeM");
 139   hCorSeM->Divide(hMCPtSeM,hMCPtRec,1,1,"B");
 140   hCorSeM->Scale(fractionMaterial);// rescale material correction
 141   hCorSeM->Multiply(hDataPt);
 142
 143   TH1D * hCorSeW = (TH1D*) hMCPtSeW->Clone("hCorSeW");
 144   hCorSeW->Divide(hMCPtSeW,hMCPtRec,1,1,"B");
 145   hCorSeW->Scale(fractionWeak);// rescale weak correction
 146   hCorSeW->Multiply(hDataPt);
 147
 148   TH1D * hCorFak = (TH1D*) hMCPtFak->Clone("hCorFak");
 149   hCorFak->Divide(hMCPtFak,hMCPtRec,1,1,"B");
 150   hCorFak->Multiply(hDataPt);
 151
 152   TH1D * hCorrected = (TH1D*) hDataPt->Clone("hCorrected");
 153   hCorrected->Add(hCorSeM,-1);
 154   hCorrected->Add(hCorSeW,-1);
 155   hCorrected->Add(hCorFak,-1);
 156   hCorrected->Divide(hEffPt);
 157   hCorrected->SetMarkerStyle(kOpenStar);
 158
 159
 160   TCanvas * cCorrections = new TCanvas("cCorrections", "cCorrections");
 161
 162   hEffPt->Draw();
 163   // hCorSeM->Draw();
 164   // hCorSeM->SetLineColor(kRed);
 165   // hCorSeM->SetMarkerColor(kRed);
 166   // hMCPtSeM->Draw("same");
 167   // hCorSeW->Draw("same");
 168   // hCorSeW->SetLineColor(kRed);
 169   // hCorSeW->SetMarkerColor(kRed);
 170   // hMCPtSeW->Draw("same");
 171
 172
 173   // hDataPt->Draw();
 174   // return;
 175   TCanvas * cdata = new TCanvas ("cData", "Data");
 176   cdata->SetLogy();
 177   cdata->cd();
 178 #ifdef CORRECT_2D
 179   Int_t minProj = hDataPt->GetYaxis()->FindBin(zmin);
 180   Int_t maxProj = hDataPt->GetYaxis()->FindBin(zmax-0.0001);
 181
 182   cout << minProj << "-" << maxProj << endl;
 183
 184   // This correction accounts for the vertex cut;
 185
 186   hDataPt    = ((TH2D*)hDataPt   )->ProjectionX("_px",minProj,maxProj,"E");
 187   hCorrected = ((TH2D*)hCorrected)->ProjectionX("_px",minProj,maxProj,"E");
 188   hMCPtGen   = ((TH2D*)hMCPtGen  )->ProjectionX("_px",minProj,maxProj,"E");
 189
 190   Float_t vertexCutCorrection =
 191     hManCorr->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Integral(-1,-1) /
 192     hManCorr->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Integral(minProj,maxProj) ;
 193   // cout << vertexCutCorrection << " " << hMCPtGen->Integral(-1,-1)  << " " << hMCPtPri->Integral() << endl;
 194   // vertexCutCorrection /= hMCPtGen->Integral(-1,-1);
 195   vertexCutCorrection = 1; // FIXME
 196   cout << "Vertex cut correction " << vertexCutCorrection << " (Efficiency " << 1./vertexCutCorrection << ")" << endl;
 197
 198   hDataPt    ->Scale(1.,"width");
 199   hCorrected ->Scale(vertexCutCorrection,"width");
 200   hMCPtGen   ->Scale(1.,"width");
 201 #endif
 202
 203   hDataPt->Draw();
 204   hDataPt ->GetXaxis()->SetRangeUser(0,4.5);
 205   hDataPt ->GetYaxis()->SetRangeUser(0.1,1e4);
 206   hCorrected->SetLineColor(kBlack);
 207   hCorrected->SetMarkerColor(kBlack);
 208   hCorrected->Draw("same");
 209   hMCPtGen->DrawClone("same");
 210   TF1 * f = GetLevy();
 211   //TF1 * f = GetHagedorn();
 212   hCorrected->Fit(f,"", "same");
 213   hCorrected->Fit(f,"IME", "same",0,2);
 214   cout << "dN/deta (function)  = " << f->Integral(0,100) << " +- " << f->IntegralError(0,100) << endl;
 215   cout << "dN/deta (data>0.15) = " <<  hCorrected->Integral(3,-1,"width") << endl;//
 216   cout << "dN/deta (func+data) = " << f->Integral(0,0.1) + hCorrected->Integral(3,-1,"width") << endl;//
 217   cout << "dN/deta (func 0.15+data) = " << f->Integral(0,0.15) + hCorrected->Integral(4,-1,"width") << endl;//
 218   cout << "Generated dN/deta (correction) = " << hMCPtGen->Integral("width") << endl;
 219   cout << "Generated dN/deta (correction, <0.13) = " << hMCPtGen->Integral(1,2,"width") << endl;
 220   cout << "-------" << endl;
 221   Double_t errorData = 0;
 222   Float_t dNdeta  =  (f->Integral(0,0.1) + hCorrected->IntegralAndError(3,-1,errorData, "width")) / (etaMax-etaMin);
 223   Float_t dNdetaE = TMath::Sqrt(f->IntegralError(0,0.1)*f->IntegralError(0,0.1) + errorData*errorData) / (etaMax-etaMin);
 224   cout << "dN/deta          " << dNdeta << " +- " << dNdetaE << endl;
 225   cout << "(dN/deta)/Npart  " << dNdeta/npart << " +- " << dNdetaE/npart << endl;
 226
 227   // FIXME: comment this out
 228   TH1D * hDataGen  = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        etaMin,etaMax,zmin,zmax);
 229   cout << "Generated dN/deta (data) =       " << hDataGen->Integral("width") << endl;
 230   hDataGen->Draw("same");
 231   TLegend * l = new TLegend(0.520134, 0.676573, 0.885906, 0.923077,"137161, p1+++");
 232   l->AddEntry(hDataPt, "Raw data");
 233   l->AddEntry(hCorrected, "Corrected data");
 234   l->AddEntry(hMCPtGen, "Monte Carlo (generated)");
 235   l->AddEntry(f, "Levy Fit");
 236   l->Draw();
 237 }
 238
 239 void CheckSecondaries(Double_t &fracWeak, Double_t &fracMaterial) {
 240   // Returns the fraction you need to rescale the secondaries from weak decays for.
 241
 242   // Some shorthands
 243   TH1D * hDataDCA   = hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 244   //  TH1D * hMCDCAGen  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
 245   TH1D * hMCDCARec  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 246   TH1D * hMCDCAPri  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
 247   TH1D * hMCDCASW  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak    );
 248   TH1D * hMCDCASM  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat    );
 249   TH1D * hMCDCAFak  = hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
 250
 251
 252   TCanvas * cCumulative = new TCanvas("cDCAculumative","DCA cumulative distributions");
 253   cCumulative->cd();
 254   GetCumulativeHisto(hMCDCAPri )->Draw();
 255   GetRatioIntegratedFractions(hMCDCASW, hMCDCARec  )->Draw("same");
 256   GetRatioIntegratedFractions(hMCDCASM, hMCDCARec  )->Draw("same");
 257   GetRatioIntegratedFractions(hMCDCAPri,hMCDCARec  )->Draw("same");
 258
 259
 260   TCanvas * c1 = new TCanvas("cDCAFit", "Fit to the DCA distributions");
 261   c1->SetLogy();
 262   // Draw all
 263   //  hDataDCA->Draw();
 264   //  hMCDCARec ->Draw("same");
 265   // hMCDCAPri ->Draw("same");
 266   // hMCDCASW ->Draw("same");
 267   // hMCDCASM ->Draw("same");
 268   // hMCDCAFak ->Draw("same");
 269   // return;
 270
 271   // Fit the DCA distribution. Uses a TF1 made by summing histograms
 272   TH1D * hMCPrimSMFak = (TH1D*) hMCDCAPri->Clone("hMCPrimSMFak");
 273   //  hMCPrimSMFak->Add(hMCDCASM);
 274   hMCPrimSMFak->Add(hMCDCAFak);
 275
 276   // Set the components which are used in HistoSum, the static
 277   // function for GetFunctionHistoSum
 278   gHistoCompoments[0] = (TH1D*) hMCPrimSMFak->Clone();
 279   gHistoCompoments[1] = (TH1D*) hMCDCASW->Clone();
 280   gHistoCompoments[2] = (TH1D*) hMCDCASM->Clone();
 281   TF1 * fHistos = GetFunctionHistoSum();
 282   fHistos->SetParameters(1,1,1);
 283   //fHistos->FixParameter(2,1);
 284   fHistos->SetLineColor(kRed);
 285   hDataDCA->Draw();
 286   // Fit!
 287   //  hDataDCA->Fit(fHistos,"","",0,200);
 288   // Rescale the components and draw to see how it looks like
 289   hMCPrimSMFak->Scale(fHistos->GetParameter(0));
 290   hMCDCASW    ->Scale(fHistos->GetParameter(1));
 291   hMCDCASM    ->Scale(fHistos->GetParameter(2));
 292   hMCPrimSMFak->Draw("same");
 293   hMCDCASW    ->Draw("same");
 294   hMCDCASM    ->Draw("same");
 295   fHistos->Draw("same");
 296   // compute scaling factors
 297   fracWeak     = fHistos->GetParameter(1)/fHistos->GetParameter(0);
 298   fracMaterial = fHistos->GetParameter(2)/fHistos->GetParameter(0);
 299
 300
 301 }
 302
 303 void CheckSanity() {
 304   // compares various distributions in data and in MC
 305   TCanvas * c1 = new TCanvas("cVz", "Vertex Z distribution");
 306   c1->cd();
 307   TH1D * hData  = hManData->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 308   TH1D * hCorr  = hManCorr->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 309   hData->SetLineColor  (kBlack);
 310   hData->SetMarkerColor(kBlack);
 311   hData->SetMarkerStyle(kFullCircle);
 312   hCorr->SetLineColor  (kRed);
 313   hCorr->SetMarkerColor(kRed);
 314   hCorr->SetMarkerStyle(kOpenSquare);
 315   TLegend * l = new TLegend(0.575503, 0.70979, 0.86745, 0.917832);
 316   l->AddEntry(hData, "Data");
 317   l->AddEntry(hCorr, "MC");
 318   hCorr->Draw("");
 319   hData->Draw("same");
 320   l->Draw();
 321
 322
 323   TCanvas * c2 = new TCanvas("cEta", "Eta distribution");
 324   c2->cd();
 325   hData  = hManData->GetHistoEta(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 326   hCorr  = hManCorr->GetHistoEta(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 327   hData->SetLineColor  (kBlack);
 328   hData->SetMarkerColor(kBlack);
 329   hData->SetMarkerStyle(kFullCircle);
 330   hCorr->SetLineColor  (kRed);
 331   hCorr->SetMarkerColor(kRed);
 332   hCorr->SetMarkerStyle(kOpenSquare);
 333   l = new TLegend(0.575503, 0.70979, 0.86745, 0.917832);
 334   l->AddEntry(hData, "Data");
 335   l->AddEntry(hCorr, "MC");
 336   hData->Draw("");
 337   hCorr->Draw("same");
 338   l->Draw();
 339
 340   TCanvas * c3 = new TCanvas("cMult", "Multiplicity distribution");
 341   c3->cd();
 342   hData  = hManData->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 343   hCorr  = hManCorr->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
 344   hData->SetLineColor  (kBlack);
 345   hData->SetMarkerColor(kBlack);
 346   hData->SetMarkerStyle(kFullCircle);
 347   hCorr->SetLineColor  (kRed);
 348   hCorr->SetMarkerColor(kRed);
 349   hCorr->SetMarkerStyle(kOpenSquare);
 350   l = new TLegend(0.575503, 0.70979, 0.86745, 0.917832);
 351   l->AddEntry(hData, "Data");
 352   l->AddEntry(hCorr, "MC");
 353   hData->Draw("");
 354   hCorr->Draw("same");
 355   l->Draw();
 356
 357 }
 358
 359 void LoadLibs(  Bool_t debug) {
 360
 361   gSystem->Load("libVMC");
 362   gSystem->Load("libTree");
 363   gSystem->Load("libSTEERBase");
 364   gSystem->Load("libESD");
 365   gSystem->Load("libAOD");
 366   gSystem->Load("libANALYSIS");
 367   gSystem->Load("libANALYSISalice");
 368   gSystem->Load("libCORRFW");
 369   gSystem->Load("libMinuit");
 370   gSystem->Load("libPWG2Spectra");
 371   gSystem->Load("libPWG0base");
 372
 373   gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG0/multPbPb"));
 374   gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG1/background"));
 375   // Load helper classes
 376   // TODO: replace this by a list of TOBJStrings
 377   TString taskName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisTaskMultPbTracks.cxx+");
 378   TString histoManName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx+");
 379   TString centrName("$ALICE_ROOT/PWG0/multPbPb/AliAnalysisMultPbCentralitySelector.cxx+");
 380   TString listName("$ALICE_ROOT/PWG1/background/AliHistoListWrapper.cxx+");
 381
 382   gSystem->ExpandPathName(taskName);
 383   gSystem->ExpandPathName(histoManName);
 384   gSystem->ExpandPathName(centrName);
 385   gSystem->ExpandPathName(listName);
 386
 387
 388   gROOT->LoadMacro(listName    +(debug?"+g":""));
 389   gROOT->LoadMacro(histoManName+(debug?"+g":""));
 390   gROOT->LoadMacro(centrName   +(debug?"+g":""));
 391   gROOT->LoadMacro(taskName    +(debug?"+g":""));
 392
 393   // Histo fitter
 394   gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/fcn.cxx+g");
 395   gROOT->LoadMacro("/Users/mfloris/Work/ALICE/ANALYSIS/HistoFitter/AliHistoFitter.cxx+g");
 396
 397
 398 }
 399
 400
 401 void SetStyle() {
 402
 403   hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kBlack);
 404   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetLineColor(kRed  );
 405   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kRed  );
 406   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
 407   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
 408   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue );
 409   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
 410
 411   hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kBlack);
 412   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerColor(kRed  );
 413   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kRed  );
 414   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
 415   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
 416   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue );
 417   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
 418
 419   hManData->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kFullCircle);
 420   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerStyle(kFullSquare);
 421   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kOpenSquare);
 422   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
 423   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
 424   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenCircle);
 425   hManCorr->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
 426
 427  hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kBlack);
 428   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetLineColor(kRed  );
 429   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetLineColor(kRed  );
 430   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetLineColor(kGreen);
 431   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetLineColor(kBlue );
 432   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetLineColor(kBlue-7 );
 433   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetLineColor(kCyan );
 434
 435   hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kBlack);
 436   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerColor(kRed  );
 437   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerColor(kRed  );
 438   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerColor(kGreen);
 439   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerColor(kBlue );
 440   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerColor(kBlue-7 );
 441   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerColor(kCyan );
 442
 443   hManData->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kFullCircle);
 444   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen)    ->SetMarkerStyle(kFullSquare);
 445   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec)    ->SetMarkerStyle(kOpenSquare);
 446   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->SetMarkerStyle(kOpenSquare);
 447   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat) ->SetMarkerStyle(kOpenSquare);
 448   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak) ->SetMarkerStyle(kOpenSquare);
 449   hManCorr->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake)->SetMarkerStyle(kOpenSquare);
 450
 451
 452 }
 453
 454 void LoadData(TString dataFolder, TString correctionFolder){
 455
 456   // Get histo manager for data and MC + stat histos
 457   TFile * fData = new TFile(dataFolder+"/multPbPbtracks.root");
 458   TFile * fCorr = new TFile(correctionFolder+"/multPbPbtracks.root");
 459   TFile * fStatData = new TFile(dataFolder+"/event_stat.root");
 460   TFile * fStatCorr = new TFile(correctionFolder+"/event_stat.root");
 461
 462   hManData = (AliAnalysisMultPbTrackHistoManager*) fData->Get("histoManager");
 463   hManCorr = (AliAnalysisMultPbTrackHistoManager*) fCorr->Get("histoManager");
 464   AliESDtrackCuts * cutsData = (AliESDtrackCuts*) fData->Get("AliESDtrackCuts");
 465   AliESDtrackCuts * cutsCorr = (AliESDtrackCuts*) fCorr->Get("AliESDtrackCuts");
 466   if (cutsData != cutsCorr) {
 467     cout << "ERROR: MC and data do not have the same cuts" << endl;
 468     // FIXME: exit here
 469   }
 470   cutsData->Print();
 471   hEvStatData = (TH2D*) fStatData->Get("fHistStatistics");
 472   hEvStatCorr = (TH2D*) fStatCorr->Get("fHistStatistics");
 473
 474   AliAnalysisMultPbCentralitySelector * centrData = (AliAnalysisMultPbCentralitySelector*) fData->Get("Cuts");
 475   if(!centrData) {
 476     cout << "ERROR:  cannot read centrality data" << endl;
 477   }
 478   centrData->Print();
 479   // Normalize
 480   Int_t irowGoodTrigger = 1;
 481   if (hEvStatCorr && hEvStatData) {
 482     //  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
 483     // hManData->ScaleHistos(hEvStatData->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
 484     // hManCorr->ScaleHistos(hEvStatCorr->GetBinContent(AliPhysicsSelection::kStatAccepted,irowGoodTrigger));
 485     // hManData->ScaleHistos(hManData->GetHistoStats()->GetBinContent(AliAnalysisMultPbTrackHistoManager::kStatVtx));
 486     // hManCorr->ScaleHistos(hManCorr->GetHistoStats()->GetBinContent(AliAnalysisMultPbTrackHistoManager::kStatVtx));
 487     TH1D* hvzData = hManData->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec);
 488     TH1D* hvzCorr = hManCorr->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec);
 489     hManData->ScaleHistos(hvzData->Integral(hvzData->FindBin(zmin),hvzData->FindBin(zmax-0.001)));
 490     hManCorr->ScaleHistos(hvzCorr->Integral(hvzCorr->FindBin(zmin),hvzCorr->FindBin(zmax-0.001)));
 491   } else {
 492     cout << "WARNING!!! ARBITRARY SCALING" << endl;
 493     hManData->ScaleHistos(1000);
 494     hManCorr->ScaleHistos(1000);
 495   }
 496 }
 497
 498 TF1 * GetHagedorn(Float_t norm, Float_t pt0, Float_t n) {
 499
 500   TF1 * f =0;
 501   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 502
 503   f=new TF1("fHagedorn",Form("(x/sqrt(x*x+%f*%f))*x*[0]*( 1 + x/[1] )^(-[2])",mass,mass),0,10);
 504   f->SetParameters(norm, pt0, n);
 505   f->SetParLimits(1, 0.01, 10);
 506   f->SetParNames("norm", "pt0", "n");
 507   f->SetLineWidth(1);
 508   return f;
 509
 510
 511 }
 512
 513 TF1 * GetMTExp(Float_t norm, Float_t t) {
 514
 515   TF1 * f =0;
 516   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 517
 518   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);
 519   f->SetParameters(norm, t);
 520   //  f->SetParLimits(1, 0.01);
 521   f->SetParNames("norm", "T");
 522   f->SetLineWidth(1);
 523   return f;
 524
 525
 526 }
 527
 528 TF1 * GetLevy(Double_t temp, Double_t n, Double_t norm, const char * name) {
 529
 530   char formula[500];
 531   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 532
 533
 534   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])");
 535   TF1* f=new TF1(name,formula,0,10);
 536   f->SetParameters(norm, n, temp,mass);
 537   f->SetParLimits(2, 0.01, 10);
 538   f->SetParNames("norm (dN/dy)", "n", "T", "mass");
 539   f->FixParameter(3,mass);
 540   f->SetLineWidth(1);
 541   return f;
 542
 543
 544 }
 545
 546
 547 TH1D * GetCumulativeHisto (TH1 * h) {
 548   // Returns a cumulative histogram
 549   // FIXME: put this method in a tools class
 550   TH1D * hInt = (TH1D*) h->Clone(TString(h->GetName())+"_Int");
 551   hInt->Reset();
 552   Float_t integral = h->Integral(-1,-1); // Consider under/over flow!
 553   Int_t nbin = h->GetNbinsX();
 554   for(Int_t ibin = 1; ibin <= nbin; ibin++){
 555     Double_t content = h->Integral(-1,ibin) / integral;
 556     hInt->SetBinContent(ibin, content);
 557   }
 558   return hInt;
 559 }
 560
 561 TH1D * GetRatioIntegratedFractions (TH1 * hNum, TH1 * hDenum) {
 562   // Returns the the ratio of integrated histograms
 563   // FIXME: put this method in a tools class
 564   TH1D * hRatio = (TH1D*) hNum->Clone(TString(hNum->GetName())+hDenum->GetName()+"_RatioIntegrated");
 565   hRatio->Reset();
 566   Int_t nbin = hNum->GetNbinsX();
 567   for(Int_t ibin = 1; ibin <= nbin; ibin++){
 568     Double_t content = hNum->Integral(-1,ibin) / hDenum->Integral(-1,ibin);// consider underflow
 569     hRatio->SetBinContent(ibin, content);
 570   }
 571   return hRatio;
 572 }
 573
 574 void ShowAcceptanceInVzSlices() {
 575   TCanvas * cvz = new TCanvas("cvz","acc #times eff vs vz");
 576   for(Int_t ivz = -10; ivz < -6; ivz+=2){
 577     ivz=0;//FIXME
 578     Bool_t first = kTRUE;
 579     TH1D * hMCPtPri  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim ,   etaMin,etaMax,ivz,ivz+2);
 580     TH1D * hMCPtGen  = hManCorr->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        etaMin,etaMax,ivz,ivz+2);
 581     TH1D * hMCPtPriD  = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim ,   etaMin,etaMax,ivz,ivz+2);
 582     TH1D * hMCPtGenD  = hManData->GetHistoPt(AliAnalysisMultPbTrackHistoManager::kHistoGen,        etaMin,etaMax,ivz,ivz+2);
 583     //    hEff= hMCPtGen;
 584     TH1D * hEff = (TH1D*)hMCPtPri->Clone(Form("hEff_vz_%d_%d",ivz,ivz+2));
 585     hEff->Divide(hMCPtPri,hMCPtGen,1,1,"B");
 586     TH1D * hEffD = (TH1D*)hMCPtPriD->Clone(Form("hEffD_vz_%d_%d",ivz,ivz+2));
 587     hEffD->Divide(hMCPtPriD,hMCPtGenD,1,1,"B");
 588     hEffD->SetLineColor(kRed);
 589     cout << "ivz " << ivz << endl;
 590     if(first) {
 591       first=kFALSE;
 592       cout << "First" << endl;
 593       hEff->Draw();
 594       hEffD->Draw("same");
 595       // hMCPtGen->Draw();
 596       // hMCPtPri->Draw("same");
 597     }
 598     else {
 599       cout << "Same" << endl;
 600       hEff->Draw("same");
 601       hEffD->Draw("same");
 602       // hMCPtGen->Draw("");
 603       // hMCPtPri->Draw("same");
 604     }
 605     cvz->Update();
 606     //    cvz->WaitPrimitive();
 607   }
 608
 609   //hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim )->Draw();
 610   //  hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoRec )->Draw("");
 611   // hManCorr->GetHistoVz(AliAnalysisMultPbTrackHistoManager::kHistoGen     )->Draw("same");
 612
 613
 614 }
 615
 616 TF1 * GetFunctionHistoSum() {
 617
 618   TF1 * f = new TF1 ("fHistoSum",HistoSum,0,50,kHistoFitCompoments);
 619   f->SetParNames("Primaries+Fakes", "Sec. Weak decays","Sec. Material");
 620   f->SetNpx(1000);
 621   return f;
 622
 623 }
 624
 625 static Double_t HistoSum(const double * x, const double* p){
 626   // This function uses a global array of histograms, rescaled by some
 627   // parameters, to define a function
 628   // The array is called gHistoCompoments
 629   // The size of the arrays is given by the global constant kHistoFitCompoments
 630
 631   Double_t xx = x[0];
 632   Double_t value = 0;
 633   for(Int_t icomp = 0; icomp < kHistoFitCompoments; icomp++){
 634     //    value += gHistoCompoments[icomp]-Interpolate(xx) * p[icomp];
 635     Int_t ibin = gHistoCompoments[icomp]->FindBin(xx);
 636     value += gHistoCompoments[icomp]->GetBinContent(ibin) * p[icomp];
 637   }
 638
 639   return value;
 640
 641 }
 642
 643 void PrintCanvas(TCanvas* c,const TString formats) {
 644   // print a canvas in every of the given comma-separated formats
 645   // ensure the canvas is updated
 646   if(!doPrint) return;
 647   c->Update();
 648   gSystem->ProcessEvents();
 649   TObjArray * arr = formats.Tokenize(",");
 650   TIterator * iter = arr->MakeIterator();
 651   TObjString * element = 0;
 652   TString name  =c ->GetName();
 653   name.ReplaceAll(" ","_");
 654   name.ReplaceAll("+","Plus");
 655   name.ReplaceAll("-","");
 656   while ((element = (TObjString*) iter->Next())) {
 657     c->Print(name+ "."+element->GetString());
 658   }
 659 }