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