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