]>
Commit | Line | Data |
---|---|---|
10ebe68d | 1 | /* $Id$ */ |
2 | ||
3 | #if !defined(__CINT__) || defined(__MAKECINT__) | |
4 | ||
5 | #include "AliPWG0Helper.h" | |
6 | #include "dNdEtaAnalysis.h" | |
7 | #include "AlidNdEtaCorrection.h" | |
8 | ||
9 | #include <TCanvas.h> | |
10 | #include <TFile.h> | |
11 | #include <TH1.h> | |
12 | #include <TH2F.h> | |
13 | #include <TH3F.h> | |
14 | #include <TLine.h> | |
15 | #include <TSystem.h> | |
16 | ||
17 | #endif | |
18 | ||
9f469bf5 | 19 | extern TPad* gPad; |
20 | ||
10ebe68d | 21 | void SetRanges(TAxis* axis) |
22 | { | |
23 | if (strcmp(axis->GetTitle(), "#eta") == 0) | |
24 | axis->SetRangeUser(-1.7999, 1.7999); | |
25 | if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0) | |
26 | axis->SetRangeUser(0, 9.9999); | |
27 | if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0) | |
28 | axis->SetRangeUser(-15, 14.9999); | |
29 | if (strcmp(axis->GetTitle(), "Ntracks") == 0) | |
30 | axis->SetRangeUser(0, 99.9999); | |
31 | } | |
32 | ||
33 | void SetRanges(TH1* hist) | |
34 | { | |
35 | SetRanges(hist->GetXaxis()); | |
36 | SetRanges(hist->GetYaxis()); | |
37 | SetRanges(hist->GetZaxis()); | |
38 | } | |
39 | ||
40 | void Prepare3DPlot(TH3* hist) | |
41 | { | |
42 | hist->GetXaxis()->SetTitleOffset(1.5); | |
43 | hist->GetYaxis()->SetTitleOffset(1.5); | |
44 | hist->GetZaxis()->SetTitleOffset(1.5); | |
45 | ||
46 | hist->SetStats(kFALSE); | |
47 | } | |
48 | ||
49 | void Prepare2DPlot(TH2* hist) | |
50 | { | |
51 | hist->SetStats(kFALSE); | |
52 | hist->GetYaxis()->SetTitleOffset(1.4); | |
53 | ||
54 | SetRanges(hist); | |
55 | } | |
56 | ||
57 | void Prepare1DPlot(TH1* hist) | |
58 | { | |
59 | hist->SetLineWidth(2); | |
60 | hist->SetStats(kFALSE); | |
61 | ||
62 | SetRanges(hist); | |
63 | } | |
64 | ||
65 | void InitPad() | |
66 | { | |
9f469bf5 | 67 | if (!gPad) |
68 | return; | |
69 | ||
10ebe68d | 70 | gPad->Range(0, 0, 1, 1); |
71 | gPad->SetLeftMargin(0.15); | |
72 | //gPad->SetRightMargin(0.05); | |
73 | //gPad->SetTopMargin(0.13); | |
74 | //gPad->SetBottomMargin(0.1); | |
75 | ||
76 | //gPad->SetGridx(); | |
77 | //gPad->SetGridy(); | |
78 | } | |
79 | ||
80 | void InitPadCOLZ() | |
81 | { | |
82 | gPad->Range(0, 0, 1, 1); | |
83 | gPad->SetRightMargin(0.15); | |
84 | gPad->SetLeftMargin(0.12); | |
85 | ||
86 | gPad->SetGridx(); | |
87 | gPad->SetGridy(); | |
88 | } | |
89 | ||
90 | void Secondaries() | |
91 | { | |
92 | TFile* file = TFile::Open("systematics.root"); | |
93 | ||
94 | TH3F* secondaries = dynamic_cast<TH3F*> (file->Get("fSecondaries")); | |
95 | if (!secondaries) | |
96 | { | |
97 | printf("Could not read histogram\n"); | |
98 | return; | |
99 | } | |
100 | ||
101 | for (Int_t ptBin=1; ptBin<=secondaries->GetNbinsZ(); ptBin++) | |
102 | //for (Int_t ptBin = 1; ptBin<=2; ptBin++) | |
103 | { | |
9f469bf5 | 104 | TH1F* hist = new TH1F(Form("secondaries_%d", ptBin), Form("secondaries_%d", ptBin), secondaries->GetNbinsY(), secondaries->GetYaxis()->GetXmin(), secondaries->GetYaxis()->GetXmax()); |
105 | ||
106 | hist->SetTitle(Form("%f < p_{T} < %f", secondaries->GetZaxis()->GetBinLowEdge(ptBin), secondaries->GetZaxis()->GetBinUpEdge(ptBin))); | |
10ebe68d | 107 | |
108 | for (Int_t cBin=1; cBin<=secondaries->GetNbinsY(); ++cBin) | |
109 | { | |
110 | if (secondaries->GetBinContent(0, cBin, ptBin) > 0) | |
111 | printf("WARNING: Underflow bin not empty!"); | |
112 | if (secondaries->GetBinContent(secondaries->GetNbinsX()+1, cBin, ptBin) > 0) | |
113 | printf("WARNING: Overflow bin not empty!"); | |
114 | ||
115 | Double_t sum = 0; | |
116 | Double_t count = 0; | |
117 | for (Int_t nBin=1; nBin<=secondaries->GetNbinsX(); ++nBin) | |
118 | { | |
119 | //printf("%f %f\n", secondaries->GetXaxis()->GetBinCenter(nBin), secondaries->GetBinContent(nBin, cBin, ptBin)); | |
120 | sum += secondaries->GetXaxis()->GetBinCenter(nBin) * secondaries->GetBinContent(nBin, cBin, ptBin); | |
121 | count += secondaries->GetBinContent(nBin, cBin, ptBin); | |
122 | } | |
123 | ||
124 | printf("%f %f\n", sum, count); | |
125 | ||
126 | if (count > 0) | |
9f469bf5 | 127 | { |
128 | hist->SetBinContent(cBin, sum / count); | |
129 | hist->SetBinError(cBin, TMath::Sqrt(sum) / count); | |
130 | } | |
10ebe68d | 131 | } |
132 | ||
9f469bf5 | 133 | hist->Scale(1.0 / hist->GetBinContent(hist->GetXaxis()->FindBin(1))); |
134 | hist->Add(new TF1("flat", "-1", 0, 2)); | |
135 | ||
10ebe68d | 136 | new TCanvas; |
9f469bf5 | 137 | hist->SetMarkerStyle(21); |
138 | hist->Draw(""); | |
10ebe68d | 139 | } |
140 | } | |
141 | ||
9f469bf5 | 142 | void Track2Particle1DComposition(const char* fileName = "correction_map.root", Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9) |
10ebe68d | 143 | { |
144 | gSystem->Load("libPWG0base"); | |
145 | ||
9f469bf5 | 146 | TString canvasName; |
147 | canvasName.Form("Track2Particle1DComposition"); | |
148 | TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400); | |
149 | canvas->Divide(3, 1); | |
10ebe68d | 150 | |
9f469bf5 | 151 | TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95); |
10ebe68d | 152 | |
9f469bf5 | 153 | for (Int_t i=0; i<folderCount; ++i) |
10ebe68d | 154 | { |
9f469bf5 | 155 | Track2Particle1DCreatePlots(fileName, folderNames[i], upperPtLimit); |
156 | ||
157 | TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_x_div_meas_nTrackToNPart_x")); | |
158 | TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_y_div_meas_nTrackToNPart_y")); | |
159 | TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_z_div_meas_nTrackToNPart_z")); | |
160 | ||
161 | Prepare1DPlot(corrX); | |
162 | Prepare1DPlot(corrY); | |
163 | Prepare1DPlot(corrZ); | |
164 | ||
165 | const char* title = "Track2Particle Correction"; | |
166 | corrX->SetTitle(title); | |
167 | corrY->SetTitle(title); | |
168 | corrZ->SetTitle(title); | |
169 | ||
170 | corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit); | |
171 | ||
172 | corrX->SetLineColor(i+1); | |
173 | corrY->SetLineColor(i+1); | |
174 | corrZ->SetLineColor(i+1); | |
175 | ||
176 | canvas->cd(1); | |
177 | InitPad(); | |
178 | corrX->DrawCopy(((i>0) ? "SAME" : "")); | |
179 | ||
180 | canvas->cd(2); | |
181 | InitPad(); | |
182 | corrY->DrawCopy(((i>0) ? "SAME" : "")); | |
183 | ||
184 | canvas->cd(3); | |
185 | InitPad(); | |
186 | corrZ->DrawCopy(((i>0) ? "SAME" : "")); | |
187 | ||
188 | legend->AddEntry(corrZ, folderNames[i]); | |
10ebe68d | 189 | } |
190 | ||
9f469bf5 | 191 | legend->Draw(); |
192 | ||
193 | //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.gif", fileName, gMax, upperPtLimit)); | |
194 | //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.eps", fileName, gMax, upperPtLimit)); | |
195 | } | |
196 | ||
197 | void DrawDifferentSpecies() | |
198 | { | |
199 | gROOT->ProcessLine(".L drawPlots.C"); | |
200 | ||
201 | const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" }; | |
202 | ||
203 | Track2Particle1DComposition("systematics.root", 4, folderNames); | |
204 | } | |
10ebe68d | 205 | |
9f469bf5 | 206 | void ScalePtDependent(TH3F* hist, TF1* function) |
207 | { | |
208 | // assumes that pt is the third dimension of hist | |
209 | // scales with function(pt) | |
10ebe68d | 210 | |
9f469bf5 | 211 | for (Int_t z=1; z<=hist->GetNbinsZ(); ++z) |
212 | { | |
213 | Double_t factor = function->Eval(hist->GetZaxis()->GetBinCenter(z)); | |
214 | printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor); | |
215 | ||
216 | for (Int_t x=1; x<=hist->GetNbinsX(); ++x) | |
217 | for (Int_t y=1; y<=hist->GetNbinsY(); ++y) | |
218 | hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor); | |
219 | } | |
220 | } | |
221 | ||
222 | void ChangeComposition(void** correctionPointer, Int_t index) | |
223 | { | |
224 | AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer; | |
225 | ||
226 | switch (index) | |
227 | { | |
228 | case 0: | |
229 | break; | |
230 | ||
231 | case 1: | |
232 | fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Scale(10); | |
233 | fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Scale(10); | |
234 | break; | |
235 | ||
236 | case 2: | |
237 | TF1* ptDependence = new TF1("simple", "x", 0, 100); | |
238 | ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDependence); | |
239 | ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDependence); | |
240 | break; | |
241 | ||
242 | } | |
243 | } | |
244 | ||
245 | void Composition() | |
246 | { | |
247 | gSystem->Load("libPWG0base"); | |
248 | ||
249 | gSystem->Unlink("new_compositions.root"); | |
250 | ||
251 | for (Int_t comp = 0; comp < 3; ++comp) | |
252 | { | |
253 | AlidNdEtaCorrection* fdNdEtaCorrection[4]; | |
254 | ||
255 | TFile::Open("systematics.root"); | |
256 | ||
257 | for (Int_t i=0; i<4; ++i) | |
258 | { | |
259 | TString name; | |
260 | name.Form("correction_%d", i); | |
261 | fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name); | |
262 | fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name); | |
263 | } | |
264 | ||
265 | ChangeComposition(fdNdEtaCorrection, comp); | |
266 | ||
267 | Double_t geneCount[5]; | |
268 | Double_t measCount[5]; | |
269 | geneCount[4] = 0; | |
270 | measCount[4] = 0; | |
271 | ||
272 | for (Int_t i=0; i<4; ++i) | |
273 | { | |
274 | geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Integral(); | |
275 | measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Integral(); | |
276 | ||
277 | geneCount[4] += geneCount[i]; | |
278 | measCount[4] += measCount[i]; | |
279 | ||
280 | printf("Particle %s: %d gene, %d meas\n", ((i == 0) ? "pi" : (i == 1) ? "K" : (i == 2) ? "p" : "others"), (Int_t) geneCount[i], (Int_t) measCount[i]); | |
281 | } | |
10ebe68d | 282 | |
9f469bf5 | 283 | printf("Generated ratios are: %f pi, %f K, %f p, %f others\n", geneCount[0] / geneCount[4], geneCount[1] / geneCount[4], geneCount[2] / geneCount[4], geneCount[3] / geneCount[4]); |
10ebe68d | 284 | |
9f469bf5 | 285 | printf("Reconstructed ratios are: %f pi, %f K, %f p, %f others\n", measCount[0] / measCount[4], measCount[1] / measCount[4], measCount[2] / measCount[4], measCount[3] / measCount[4]); |
10ebe68d | 286 | |
9f469bf5 | 287 | TList* collection = new TList; |
10ebe68d | 288 | |
9f469bf5 | 289 | for (Int_t i=0; i<4; ++i) |
290 | collection->Add(fdNdEtaCorrection[i]); | |
10ebe68d | 291 | |
9f469bf5 | 292 | TString correctionName; |
293 | correctionName.Form("new_composition_%d", comp); | |
294 | ||
295 | AlidNdEtaCorrection* newComposition = new AlidNdEtaCorrection(correctionName, correctionName); | |
296 | newComposition->Merge(collection); | |
297 | newComposition->Finish(); | |
298 | ||
299 | delete collection; | |
300 | ||
301 | TFile* file = TFile::Open("new_compositions.root", "UPDATE"); | |
302 | newComposition->SaveHistograms(); | |
303 | //file->Write(); | |
304 | file->Close(); | |
305 | } | |
10ebe68d | 306 | |
307 | gROOT->ProcessLine(".L drawPlots.C"); | |
9f469bf5 | 308 | |
309 | const char* folderNames[] = { "new_composition_0", "new_composition_1", "new_composition_2" }; | |
310 | ||
311 | Track2Particle1DComposition("new_compositions.root", 3, folderNames); | |
312 | } | |
313 | ||
314 | Double_t ConvSigma1To2D(Double_t sigma) | |
315 | { | |
316 | return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2); | |
317 | } | |
318 | ||
319 | Double_t ConvDistance1DTo2D(Double_t distance) | |
320 | { | |
321 | return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2); | |
322 | } | |
323 | ||
324 | Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma) | |
325 | { | |
326 | Double_t count = 0; | |
327 | ||
328 | //nSigma = ConvSigma1To2D(nSigma); | |
329 | ||
330 | for (Int_t x=1; x<=tracks->GetNbinsX(); ++x) | |
331 | for (Int_t y=1; y<=tracks->GetNbinsY(); ++y) | |
332 | { | |
333 | Double_t impactX = tracks->GetXaxis()->GetBinCenter(x); | |
334 | Double_t impactY = tracks->GetYaxis()->GetBinCenter(y); | |
335 | ||
336 | Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY); | |
337 | ||
338 | d = ConvDistance1DTo2D(d); | |
339 | ||
340 | if (d < nSigma) | |
341 | count += tracks->GetBinContent(x, y); | |
342 | } | |
343 | ||
344 | return count; | |
345 | } | |
346 | ||
347 | TH2F* Sigma2VertexGaussianTracksHist() | |
348 | { | |
349 | TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5); | |
350 | ||
351 | TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5); | |
352 | gaussian2D->SetParameters(1, 0, 1, 1, 0, 1); | |
353 | ||
354 | for (Int_t x=1; x<=tracks->GetNbinsX(); ++x) | |
355 | for (Int_t y=1; y<=tracks->GetNbinsY(); ++y) | |
356 | tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y))); | |
357 | ||
358 | //normalize | |
359 | tracks->Scale(1.0 / tracks->Integral()); | |
360 | ||
361 | return tracks; | |
362 | } | |
363 | ||
364 | void Sigma2VertexGaussian() | |
365 | { | |
366 | TH2F* tracks = Sigma2VertexGaussianTracksHist(); | |
367 | ||
368 | TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000); | |
369 | canvas->Divide(2, 2); | |
370 | ||
371 | canvas->cd(1); | |
372 | tracks->Draw("COLZ"); | |
373 | ||
374 | TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 10, 0.25, 5.25); | |
375 | for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5) | |
376 | ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma)); | |
377 | ratio->SetMarkerStyle(21); | |
378 | ||
379 | canvas->cd(2); | |
380 | ratio->Draw("P"); | |
381 | ||
382 | TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 3 sigma / % included n sigma", 10, 0.25, 5.25); | |
383 | Double_t sigma3 = Sigma2VertexCount(tracks, 3); | |
384 | for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5) | |
385 | ratio2->Fill(nSigma, sigma3 / Sigma2VertexCount(tracks, nSigma)); | |
386 | ratio2->SetMarkerStyle(21); | |
387 | ||
388 | canvas->cd(3); | |
389 | ratio2->Draw("P"); | |
390 | ||
391 | canvas->SaveAs("Sigma2Vertex.eps"); | |
392 | } | |
393 | ||
394 | void Sigma2VertexSimulation() | |
395 | { | |
396 | TFile* file = TFile::Open("systematics.root"); | |
397 | ||
398 | TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertex")); | |
399 | if (!sigmavertex) | |
400 | { | |
401 | printf("Could not read histogram\n"); | |
402 | return; | |
403 | } | |
404 | ||
405 | TH1F* ratio = new TH1F("sigmavertex_ratio", "sigmavertex_ratio;Nsigma;% included 3 sigma / % included n sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax()); | |
406 | ||
407 | for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i) | |
408 | ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(3)) / sigmavertex->GetBinContent(i)); | |
409 | ||
410 | TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500); | |
411 | canvas->Divide(2, 1); | |
412 | ||
413 | canvas->cd(1); | |
414 | sigmavertex->SetMarkerStyle(21); | |
415 | sigmavertex->Draw("P"); | |
416 | ||
417 | canvas->cd(2); | |
418 | ratio->SetMarkerStyle(21); | |
419 | ratio->Draw("P"); | |
10ebe68d | 420 | } |
421 | ||
9f469bf5 | 422 | |
10ebe68d | 423 | void drawSystematics() |
424 | { | |
9f469bf5 | 425 | //Secondaries(); |
426 | //DrawDifferentSpecies(); | |
427 | //Composition(); | |
428 | ||
429 | Sigma2VertexSimulation(); | |
10ebe68d | 430 | } |