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