]>
Commit | Line | Data |
---|---|---|
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 | #include <TLegend.h> | |
17 | #include <TPad.h> | |
18 | #include <TF1.h> | |
19 | ||
20 | extern TPad* gPad; | |
21 | ||
22 | void Track2Particle1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 10); | |
23 | ||
24 | #endif | |
25 | ||
26 | ||
27 | void SetRanges(TAxis* axis) | |
28 | { | |
29 | if (strcmp(axis->GetTitle(), "#eta") == 0) | |
30 | axis->SetRangeUser(-1.7999, 1.7999); | |
31 | if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0) | |
32 | axis->SetRangeUser(0, 9.9999); | |
33 | if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0) | |
34 | axis->SetRangeUser(-15, 14.9999); | |
35 | if (strcmp(axis->GetTitle(), "Ntracks") == 0) | |
36 | axis->SetRangeUser(0, 99.9999); | |
37 | } | |
38 | ||
39 | void SetRanges(TH1* hist) | |
40 | { | |
41 | SetRanges(hist->GetXaxis()); | |
42 | SetRanges(hist->GetYaxis()); | |
43 | SetRanges(hist->GetZaxis()); | |
44 | } | |
45 | ||
46 | void Prepare3DPlot(TH3* hist) | |
47 | { | |
48 | hist->GetXaxis()->SetTitleOffset(1.5); | |
49 | hist->GetYaxis()->SetTitleOffset(1.5); | |
50 | hist->GetZaxis()->SetTitleOffset(1.5); | |
51 | ||
52 | hist->SetStats(kFALSE); | |
53 | } | |
54 | ||
55 | void Prepare2DPlot(TH2* hist) | |
56 | { | |
57 | hist->SetStats(kFALSE); | |
58 | hist->GetYaxis()->SetTitleOffset(1.4); | |
59 | ||
60 | SetRanges(hist); | |
61 | } | |
62 | ||
63 | void Prepare1DPlot(TH1* hist, Bool_t setRanges = kTRUE) | |
64 | { | |
65 | hist->SetLineWidth(2); | |
66 | hist->SetStats(kFALSE); | |
67 | ||
68 | hist->GetXaxis()->SetTitleOffset(1.2); | |
69 | hist->GetYaxis()->SetTitleOffset(1.2); | |
70 | ||
71 | if (setRanges) | |
72 | SetRanges(hist); | |
73 | } | |
74 | ||
75 | void InitPad() | |
76 | { | |
77 | if (!gPad) | |
78 | return; | |
79 | ||
80 | gPad->Range(0, 0, 1, 1); | |
81 | gPad->SetLeftMargin(0.15); | |
82 | //gPad->SetRightMargin(0.05); | |
83 | //gPad->SetTopMargin(0.13); | |
84 | //gPad->SetBottomMargin(0.1); | |
85 | ||
86 | gPad->SetGridx(); | |
87 | gPad->SetGridy(); | |
88 | } | |
89 | ||
90 | void InitPadCOLZ() | |
91 | { | |
92 | gPad->Range(0, 0, 1, 1); | |
93 | gPad->SetRightMargin(0.15); | |
94 | gPad->SetLeftMargin(0.12); | |
95 | ||
96 | gPad->SetGridx(); | |
97 | gPad->SetGridy(); | |
98 | } | |
99 | ||
100 | void Secondaries() | |
101 | { | |
102 | TFile* file = TFile::Open("systematics.root"); | |
103 | ||
104 | TH2F* secondaries = dynamic_cast<TH2F*> (file->Get("fSecondaries")); | |
105 | if (!secondaries) | |
106 | { | |
107 | printf("Could not read histogram\n"); | |
108 | return; | |
109 | } | |
110 | ||
111 | TCanvas* canvas = new TCanvas("Secondaries", "Secondaries", 1000, 1000); | |
112 | canvas->Divide(3, 3); | |
113 | for (Int_t i=1; i<=8; i++) | |
114 | { | |
115 | TH1D* hist = secondaries->ProjectionY(Form("proj_%d", i), i, i); | |
116 | hist->SetTitle(secondaries->GetXaxis()->GetBinLabel(i)); | |
117 | ||
118 | canvas->cd(i); | |
119 | hist->Draw(); | |
120 | } | |
121 | } | |
122 | ||
123 | void Track2Particle1DComposition(const char** fileNames, Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9) | |
124 | { | |
125 | gSystem->Load("libPWG0base"); | |
126 | ||
127 | TString canvasName; | |
128 | canvasName.Form("Track2Particle1DComposition"); | |
129 | TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400); | |
130 | canvas->Divide(3, 1); | |
131 | ||
132 | TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95); | |
133 | ||
134 | for (Int_t i=0; i<folderCount; ++i) | |
135 | { | |
136 | Track2Particle1DCreatePlots(fileNames[i], folderNames[i], upperPtLimit); | |
137 | ||
138 | TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", folderNames[i], folderNames[i]))); | |
139 | TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", folderNames[i], folderNames[i]))); | |
140 | TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", folderNames[i], folderNames[i]))); | |
141 | ||
142 | Prepare1DPlot(corrX); | |
143 | Prepare1DPlot(corrY); | |
144 | Prepare1DPlot(corrZ); | |
145 | ||
146 | const char* title = "Track2Particle Correction"; | |
147 | corrX->SetTitle(title); | |
148 | corrY->SetTitle(title); | |
149 | corrZ->SetTitle(title); | |
150 | ||
151 | corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit); | |
152 | ||
153 | corrX->SetLineColor(i+1); | |
154 | corrY->SetLineColor(i+1); | |
155 | corrZ->SetLineColor(i+1); | |
156 | ||
157 | canvas->cd(1); | |
158 | InitPad(); | |
159 | corrX->DrawCopy(((i>0) ? "SAME" : "")); | |
160 | ||
161 | canvas->cd(2); | |
162 | InitPad(); | |
163 | corrY->DrawCopy(((i>0) ? "SAME" : "")); | |
164 | ||
165 | canvas->cd(3); | |
166 | InitPad(); | |
167 | corrZ->DrawCopy(((i>0) ? "SAME" : "")); | |
168 | ||
169 | legend->AddEntry(corrZ, folderNames[i]); | |
170 | } | |
171 | ||
172 | legend->Draw(); | |
173 | ||
174 | //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.gif", fileName, gMax, upperPtLimit)); | |
175 | //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.eps", fileName, gMax, upperPtLimit)); | |
176 | } | |
177 | ||
178 | TH1** DrawRatios(const char* fileName = "systematics.root") | |
179 | { | |
180 | gSystem->Load("libPWG0base"); | |
181 | ||
182 | TFile* file = TFile::Open(fileName); | |
183 | ||
184 | TString canvasName; | |
185 | canvasName.Form("DrawRatios"); | |
186 | TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400); | |
187 | canvas->Divide(2, 1); | |
188 | canvas->cd(1); | |
189 | ||
190 | TH1** ptDists = new TH1*[4]; | |
191 | ||
192 | TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98); | |
193 | ||
194 | const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" }; | |
195 | const char* particleNames[] = { "#pi", "K", "p", "other" }; | |
196 | for (Int_t i=0; i<4; ++i) | |
197 | { | |
198 | AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderNames[i], folderNames[i]); | |
199 | dNdEtaCorrection->LoadHistograms(fileName, folderNames[i]); | |
200 | ||
201 | TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram(); | |
202 | ||
203 | gene->GetYaxis()->SetRangeUser(-0.8, 0.8); | |
204 | gene->GetXaxis()->SetRangeUser(-10, 10); | |
205 | ||
206 | ptDists[i] = dynamic_cast<TH1*> (gene->Project3D("z")->Clone(Form("%s_z", folderNames[i]))); | |
207 | ptDists[i]->SetTitle(";p_{T};Count"); | |
208 | if (!ptDists[i]) | |
209 | { | |
210 | printf("Problem getting distribution %d.\n", i); | |
211 | return 0; | |
212 | } | |
213 | ||
214 | ptDists[i]->GetYaxis()->SetRangeUser(1, ptDists[i]->GetMaximum()*1.1); | |
215 | ptDists[i]->GetXaxis()->SetRangeUser(0, 9.9); | |
216 | ptDists[i]->SetLineColor(i+1); | |
217 | ptDists[i]->DrawCopy((i == 0) ? "" : "SAME"); | |
218 | ptDists[i]->GetYaxis()->SetRange(0, 0); | |
219 | ||
220 | legend->AddEntry(ptDists[i], particleNames[i]); | |
221 | } | |
222 | gPad->SetLogy(); | |
223 | ||
224 | TH1* total = dynamic_cast<TH1*> (ptDists[0]->Clone("total")); | |
225 | ||
226 | for (Int_t i=1; i<4; ++i) | |
227 | total->Add(ptDists[i]); | |
228 | ||
229 | canvas->cd(2); | |
230 | for (Int_t i=0; i<4; ++i) | |
231 | { | |
232 | ptDists[i]->Divide(total); | |
233 | ptDists[i]->SetStats(kFALSE); | |
234 | ptDists[i]->SetTitle(";p_{T};Fraction of total"); | |
235 | ptDists[i]->GetYaxis()->SetRangeUser(0, 1); | |
236 | ptDists[i]->Draw((i == 0) ? "" : "SAME"); | |
237 | } | |
238 | legend->SetFillColor(0); | |
239 | legend->Draw(); | |
240 | ||
241 | canvas->SaveAs("DrawRatios.gif"); | |
242 | ||
243 | ||
244 | canvas = new TCanvas("PythiaRatios", "PythiaRatios", 500, 500); | |
245 | for (Int_t i=0; i<4; ++i) | |
246 | { | |
247 | TH1* hist = ptDists[i]->Clone(); | |
248 | hist->GetXaxis()->SetRangeUser(0, 1.9); | |
249 | hist->Draw((i == 0) ? "" : "SAME"); | |
250 | } | |
251 | legend->Draw(); | |
252 | ||
253 | canvas->SaveAs("pythiaratios.eps"); | |
254 | ||
255 | file->Close(); | |
256 | ||
257 | return ptDists; | |
258 | } | |
259 | ||
260 | void DrawCompareToReal() | |
261 | { | |
262 | gROOT->ProcessLine(".L drawPlots.C"); | |
263 | ||
264 | const char* fileNames[] = { "correction_map.root", "new_compositions.root" }; | |
265 | const char* folderNames[] = { "dndeta_correction", "PythiaRatios" }; | |
266 | ||
267 | Track2Particle1DComposition(fileNames, 2, folderNames); | |
268 | } | |
269 | ||
270 | void DrawDifferentSpecies() | |
271 | { | |
272 | gROOT->ProcessLine(".L drawPlots.C"); | |
273 | ||
274 | const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "systematics.root" }; | |
275 | const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" }; | |
276 | ||
277 | Track2Particle1DComposition(fileNames, 4, folderNames); | |
278 | } | |
279 | ||
280 | void DrawSpeciesAndCombination() | |
281 | { | |
282 | gROOT->ProcessLine(".L drawPlots.C"); | |
283 | ||
284 | const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "new_compositions.root" }; | |
285 | const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "PythiaRatios" }; | |
286 | ||
287 | Track2Particle1DComposition(fileNames, 4, folderNames); | |
288 | } | |
289 | ||
290 | void DrawSimulatedVsCombined() | |
291 | { | |
292 | gROOT->ProcessLine(".L drawPlots.C"); | |
293 | ||
294 | const char* fileNames[] = { "new_compositions.root", "new_compositions.root" }; | |
295 | const char* folderNames[] = { "Pythia", "PythiaRatios" }; | |
296 | ||
297 | Track2Particle1DComposition(fileNames, 2, folderNames); | |
298 | } | |
299 | ||
300 | void DrawBoosts() | |
301 | { | |
302 | gROOT->ProcessLine(".L drawPlots.C"); | |
303 | ||
304 | const char* fileNames[] = { "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root" }; | |
305 | const char* folderNames[] = { "PythiaRatios", "PiBoosted", "KBoosted", "pBoosted" }; | |
306 | ||
307 | Track2Particle1DComposition(fileNames, 4, folderNames); | |
308 | } | |
309 | ||
310 | TH2F* HistogramDifferences(const char* filename, const char* folder1, const char* folder2) | |
311 | { | |
312 | gSystem->Load("libPWG0base"); | |
313 | ||
314 | AlidNdEtaCorrection* fdNdEtaCorrection[2]; | |
315 | ||
316 | TFile::Open(filename); | |
317 | ||
318 | fdNdEtaCorrection[0] = new AlidNdEtaCorrection(folder1, folder1); | |
319 | fdNdEtaCorrection[0]->LoadHistograms(filename, folder1); | |
320 | ||
321 | fdNdEtaCorrection[1] = new AlidNdEtaCorrection(folder2, folder2); | |
322 | fdNdEtaCorrection[1]->LoadHistograms(filename, folder2); | |
323 | ||
324 | TH3F* hist1 = fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram(); | |
325 | TH3F* hist2 = fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetCorrectionHistogram(); | |
326 | ||
327 | //TH1F* difference = new TH1F("difference", Form(";#DeltaC_{pT, z, #eta} %s / %s;Count", folder2, folder1), 1000, 0.9, 1.1); | |
328 | TH2F* difference = new TH2F(Form("difference_%s_%s", folder1, folder2), Form(";#Sigma (C_{pT, z} %s / C_{pT, z} %s);#eta;Count", folder2, folder1), 100, 0.9, 1.1, hist1->GetYaxis()->GetNbins(), hist1->GetYaxis()->GetXmin(), hist1->GetYaxis()->GetXmax()); | |
329 | ||
330 | for (Int_t x=hist1->GetXaxis()->FindBin(-10); x<=hist1->GetXaxis()->FindBin(10); ++x) | |
331 | for (Int_t y=hist1->GetYaxis()->FindBin(-0.8); y<=hist1->GetYaxis()->FindBin(0.8); ++y) | |
332 | for (Int_t z=hist1->GetZaxis()->FindBin(0.3); z<=hist1->GetZaxis()->FindBin(9.9); ++z) | |
333 | if (hist1->GetBinContent(x, y, z) != 0) | |
334 | difference->Fill(hist2->GetBinContent(x, y, z) / hist1->GetBinContent(x, y, z), hist1->GetYaxis()->GetBinCenter(y)); | |
335 | ||
336 | difference->GetYaxis()->SetRangeUser(-0.8, 0.8); | |
337 | ||
338 | printf("Over-/Underflow bins: %d %d\n", difference->GetBinContent(0), difference->GetBinContent(difference->GetNbinsX()+1)); | |
339 | ||
340 | return difference; | |
341 | } | |
342 | ||
343 | void HistogramDifferences() | |
344 | { | |
345 | TH2F* KBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "KBoosted"); | |
346 | TH2F* pBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "pBoosted"); | |
347 | TH2F* KReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "KReduced"); | |
348 | TH2F* pReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "pReduced"); | |
349 | ||
350 | TCanvas* canvas = new TCanvas("HistogramDifferences", "HistogramDifferences", 1000, 1000); | |
351 | canvas->Divide(2, 2); | |
352 | ||
353 | canvas->cd(1); | |
354 | KBoosted->GetXaxis()->SetRangeUser(-0.05, 0.05); | |
355 | KBoosted->Draw("COLZ"); | |
356 | ||
357 | canvas->cd(2); | |
358 | KReduced->GetXaxis()->SetRangeUser(-0.05, 0.05); | |
359 | KReduced->Draw("COLZ"); | |
360 | ||
361 | canvas->cd(3); | |
362 | pBoosted->GetXaxis()->SetRangeUser(-0.02, 0.02); | |
363 | pBoosted->Draw("COLZ"); | |
364 | ||
365 | canvas->cd(4); | |
366 | pReduced->GetXaxis()->SetRangeUser(-0.02, 0.02); | |
367 | pReduced->Draw("COLZ"); | |
368 | ||
369 | canvas->SaveAs("HistogramDifferences.gif"); | |
370 | ||
371 | canvas = new TCanvas("HistogramDifferencesProfile", "HistogramDifferencesProfile", 1000, 500); | |
372 | canvas->Divide(2, 1); | |
373 | ||
374 | canvas->cd(1); | |
375 | gPad->SetBottomMargin(0.13); | |
376 | KBoosted->SetTitle("a) Ratio of correction maps"); | |
377 | KBoosted->SetStats(kFALSE); | |
378 | KBoosted->GetXaxis()->SetTitleOffset(1.4); | |
379 | KBoosted->GetXaxis()->SetLabelOffset(0.02); | |
380 | KBoosted->Draw("COLZ"); | |
381 | ||
382 | canvas->cd(2); | |
383 | gPad->SetGridx(); | |
384 | gPad->SetGridy(); | |
385 | ||
386 | TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98); | |
387 | ||
388 | for (Int_t i=0; i<4; ++i) | |
389 | { | |
390 | TH2F* hist = 0; | |
391 | TString name; | |
392 | switch (i) | |
393 | { | |
394 | case 0: hist = KBoosted; name = "K enhanced"; break; | |
395 | case 1: hist = KReduced; name = "K reduced"; break; | |
396 | case 2: hist = pBoosted; name = "p enhanced"; break; | |
397 | case 3: hist = pReduced; name = "p reduced"; break; | |
398 | } | |
399 | ||
400 | TProfile* profile = hist->ProfileY(); | |
401 | profile->SetTitle("b) Mean and RMS"); | |
402 | profile->GetXaxis()->SetRange(hist->GetYaxis()->GetFirst(), hist->GetYaxis()->GetLast()); | |
403 | profile->GetXaxis()->SetTitleOffset(1.2); | |
404 | profile->GetXaxis()->SetLabelOffset(0.02); | |
405 | profile->GetYaxis()->SetRangeUser(0.98, 1.02); | |
406 | profile->SetStats(kFALSE); | |
407 | profile->SetLineColor(i+1); | |
408 | profile->SetMarkerColor(i+1); | |
409 | profile->DrawCopy(((i > 0) ? "SAME" : "")); | |
410 | ||
411 | ||
412 | legend->AddEntry(profile, name); | |
413 | } | |
414 | ||
415 | legend->Draw(); | |
416 | canvas->SaveAs("particlecomposition_result.eps"); | |
417 | } | |
418 | ||
419 | ||
420 | void ScalePtDependent(TH3F* hist, TF1* function) | |
421 | { | |
422 | // assumes that pt is the third dimension of hist | |
423 | // scales with function(pt) | |
424 | ||
425 | for (Int_t z=1; z<=hist->GetNbinsZ(); ++z) | |
426 | { | |
427 | Double_t factor = function->Eval(hist->GetZaxis()->GetBinCenter(z)); | |
428 | printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor); | |
429 | ||
430 | for (Int_t x=1; x<=hist->GetNbinsX(); ++x) | |
431 | for (Int_t y=1; y<=hist->GetNbinsY(); ++y) | |
432 | hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor); | |
433 | } | |
434 | } | |
435 | ||
436 | void ScalePtDependent(TH3F* hist, TH1* function) | |
437 | { | |
438 | // assumes that pt is the third dimension of hist | |
439 | // scales with histogram(pt) | |
440 | ||
441 | for (Int_t z=1; z<=hist->GetNbinsZ(); ++z) | |
442 | { | |
443 | Double_t factor = function->GetBinContent(function->GetXaxis()->FindBin(hist->GetZaxis()->GetBinCenter(z))); | |
444 | printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor); | |
445 | ||
446 | for (Int_t x=1; x<=hist->GetNbinsX(); ++x) | |
447 | for (Int_t y=1; y<=hist->GetNbinsY(); ++y) | |
448 | hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor); | |
449 | } | |
450 | } | |
451 | ||
452 | const char* ChangeComposition(void** correctionPointer, Int_t index) | |
453 | { | |
454 | AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer; | |
455 | ||
456 | switch (index) | |
457 | { | |
458 | case 0: // result from pp events | |
459 | { | |
460 | TFile::Open("pythiaratios.root"); | |
461 | ||
462 | for (Int_t i=0; i<4; ++i) | |
463 | { | |
464 | TString name; | |
465 | name.Form("correction_%d", i); | |
466 | fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name); | |
467 | fdNdEtaCorrection[i]->LoadHistograms("pythiaratios.root", name); | |
468 | } | |
469 | } | |
470 | return "Pythia"; | |
471 | break; | |
472 | ||
473 | case 1: // each species rated with pythia ratios | |
474 | /*TH1** ptDists = DrawRatios("pythiaratios.root"); | |
475 | ||
476 | for (Int_t i=0; i<3; ++i) | |
477 | { | |
478 | ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]); | |
479 | ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]); | |
480 | }*/ | |
481 | return "PythiaRatios"; | |
482 | break; | |
483 | ||
484 | // one species enhanced / reduced | |
485 | case 2: // + 50% pions | |
486 | case 3: // - 50% pions | |
487 | case 4: // + 50% kaons | |
488 | case 5: // - 50% kaons | |
489 | case 6: // + 50% protons | |
490 | case 7: // - 50% protons | |
491 | Int_t correctionIndex = (index - 2) / 2; | |
492 | Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5; | |
493 | ||
494 | fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Scale(scaleFactor); | |
495 | fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Scale(scaleFactor); | |
496 | ||
497 | TString* str = new TString; | |
498 | str->Form("%s%s", (correctionIndex == 0) ? "Pi" : ((correctionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced"); | |
499 | return str->Data(); | |
500 | break; | |
501 | ||
502 | // each species rated with pythia ratios | |
503 | case 12: // + 50% pions | |
504 | case 13: // - 50% pions | |
505 | case 14: // + 50% kaons | |
506 | case 15: // - 50% kaons | |
507 | case 16: // + 50% protons | |
508 | case 17: // - 50% protons | |
509 | TH1** ptDists = DrawRatios("pythiaratios.root"); | |
510 | Int_t functionIndex = (index - 2) / 2; | |
511 | Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5; | |
512 | ptDists[functionIndex]->Scale(scaleFactor); | |
513 | ||
514 | for (Int_t i=0; i<3; ++i) | |
515 | { | |
516 | ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]); | |
517 | ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]); | |
518 | } | |
519 | TString* str = new TString; | |
520 | str->Form("%s%s", (functionIndex == 0) ? "Pi" : ((functionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced"); | |
521 | return str->Data(); | |
522 | break; | |
523 | ||
524 | case 999: | |
525 | TF1* ptDependence = new TF1("simple", "x", 0, 100); | |
526 | ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDependence); | |
527 | ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDependence); | |
528 | break; | |
529 | ||
530 | } | |
531 | ||
532 | return "noname"; | |
533 | } | |
534 | ||
535 | void Composition() | |
536 | { | |
537 | gSystem->Load("libPWG0base"); | |
538 | ||
539 | gSystem->Unlink("new_compositions.root"); | |
540 | ||
541 | Int_t nCompositions = 8; | |
542 | for (Int_t comp = 0; comp < nCompositions; ++comp) | |
543 | { | |
544 | AlidNdEtaCorrection* fdNdEtaCorrection[4]; | |
545 | ||
546 | TFile::Open("systematics.root"); | |
547 | ||
548 | for (Int_t i=0; i<4; ++i) | |
549 | { | |
550 | TString name; | |
551 | name.Form("correction_%d", i); | |
552 | fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name); | |
553 | fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name); | |
554 | } | |
555 | ||
556 | const char* newName = ChangeComposition(fdNdEtaCorrection, comp); | |
557 | ||
558 | Double_t geneCount[5]; | |
559 | Double_t measCount[5]; | |
560 | geneCount[4] = 0; | |
561 | measCount[4] = 0; | |
562 | ||
563 | for (Int_t i=0; i<4; ++i) | |
564 | { | |
565 | geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Integral(); | |
566 | measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Integral(); | |
567 | ||
568 | geneCount[4] += geneCount[i]; | |
569 | measCount[4] += measCount[i]; | |
570 | ||
571 | 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]); | |
572 | } | |
573 | ||
574 | 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]); | |
575 | ||
576 | 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]); | |
577 | ||
578 | TList* collection = new TList; | |
579 | ||
580 | for (Int_t i=0; i<4; ++i) | |
581 | collection->Add(fdNdEtaCorrection[i]); | |
582 | ||
583 | AlidNdEtaCorrection* newComposition = new AlidNdEtaCorrection(newName, newName); | |
584 | newComposition->Merge(collection); | |
585 | newComposition->Finish(); | |
586 | ||
587 | delete collection; | |
588 | ||
589 | TFile* file = TFile::Open("new_compositions.root", "UPDATE"); | |
590 | newComposition->SaveHistograms(); | |
591 | //file->Write(); | |
592 | file->Close(); | |
593 | } | |
594 | ||
595 | gROOT->ProcessLine(".L drawPlots.C"); | |
596 | ||
597 | const char* fileNames[] = {"new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root" }; | |
598 | const char* folderNames[] = { "Pythia", "PythiaRatios", "PiBoosted", "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced" }; | |
599 | ||
600 | Track2Particle1DComposition(fileNames, nCompositions, folderNames); | |
601 | } | |
602 | ||
603 | Double_t ConvSigma1To2D(Double_t sigma) | |
604 | { | |
605 | return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2); | |
606 | } | |
607 | ||
608 | Double_t ConvDistance1DTo2D(Double_t distance) | |
609 | { | |
610 | return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2); | |
611 | } | |
612 | ||
613 | Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma) | |
614 | { | |
615 | Double_t count = 0; | |
616 | ||
617 | //nSigma = ConvSigma1To2D(nSigma); | |
618 | ||
619 | for (Int_t x=1; x<=tracks->GetNbinsX(); ++x) | |
620 | for (Int_t y=1; y<=tracks->GetNbinsY(); ++y) | |
621 | { | |
622 | Double_t impactX = tracks->GetXaxis()->GetBinCenter(x); | |
623 | Double_t impactY = tracks->GetYaxis()->GetBinCenter(y); | |
624 | ||
625 | Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY); | |
626 | ||
627 | d = ConvDistance1DTo2D(d); | |
628 | ||
629 | if (d < nSigma) | |
630 | count += tracks->GetBinContent(x, y); | |
631 | } | |
632 | ||
633 | return count; | |
634 | } | |
635 | ||
636 | TH2F* Sigma2VertexGaussianTracksHist() | |
637 | { | |
638 | TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5); | |
639 | ||
640 | TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5); | |
641 | gaussian2D->SetParameters(1, 0, 1, 1, 0, 1); | |
642 | ||
643 | for (Int_t x=1; x<=tracks->GetNbinsX(); ++x) | |
644 | for (Int_t y=1; y<=tracks->GetNbinsY(); ++y) | |
645 | tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y))); | |
646 | ||
647 | //normalize | |
648 | tracks->Scale(1.0 / tracks->Integral()); | |
649 | ||
650 | return tracks; | |
651 | } | |
652 | ||
653 | TH1F* Sigma2VertexGaussian() | |
654 | { | |
655 | TH2F* tracks = Sigma2VertexGaussianTracksHist(); | |
656 | ||
657 | TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000); | |
658 | canvas->Divide(2, 2); | |
659 | ||
660 | canvas->cd(1); | |
661 | tracks->Draw("COLZ"); | |
662 | ||
663 | TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 50, 0.05, 5.05); | |
664 | for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1) | |
665 | ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma)); | |
666 | ratio->SetMarkerStyle(21); | |
667 | ||
668 | canvas->cd(2); | |
669 | ratio->DrawCopy("P"); | |
670 | ||
671 | TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 3 sigma / % included n sigma", 50, 0.05, 5.05); | |
672 | Double_t sigma3 = Sigma2VertexCount(tracks, 3); | |
673 | for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1) | |
674 | ratio2->Fill(nSigma, sigma3 / ratio->GetBinContent(ratio->FindBin(nSigma))); | |
675 | ratio2->SetMarkerStyle(21); | |
676 | ||
677 | canvas->cd(3); | |
678 | ratio2->DrawCopy("P"); | |
679 | ||
680 | canvas->SaveAs("Sigma2Vertex.eps"); | |
681 | ||
682 | return ratio2; | |
683 | } | |
684 | ||
685 | TH1F* Sigma2VertexSimulation() | |
686 | { | |
687 | TFile* file = TFile::Open("systematics.root"); | |
688 | ||
689 | TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertex")); | |
690 | if (!sigmavertex) | |
691 | { | |
692 | printf("Could not read histogram\n"); | |
693 | return; | |
694 | } | |
695 | ||
696 | TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;N#sigma;% included in 3 #sigma / % included in N#sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax()); | |
697 | ||
698 | for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i) | |
699 | ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(3)) / sigmavertex->GetBinContent(i)); | |
700 | ||
701 | TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500); | |
702 | canvas->Divide(2, 1); | |
703 | ||
704 | canvas->cd(1); | |
705 | sigmavertex->SetMarkerStyle(21); | |
706 | sigmavertex->Draw("P"); | |
707 | ||
708 | canvas->cd(2); | |
709 | ratio->SetMarkerStyle(21); | |
710 | ratio->DrawCopy("P"); | |
711 | ||
712 | return ratio; | |
713 | } | |
714 | ||
715 | void Sigma2VertexCompare() | |
716 | { | |
717 | TH1F* ratio1 = Sigma2VertexGaussian(); | |
718 | TH1F* ratio2 = Sigma2VertexSimulation(); | |
719 | ||
720 | ratio1->SetStats(kFALSE); | |
721 | ratio2->SetStats(kFALSE); | |
722 | ||
723 | ratio1->SetMarkerStyle(0); | |
724 | ratio2->SetMarkerStyle(0); | |
725 | ||
726 | ratio1->SetLineWidth(2); | |
727 | ratio2->SetLineWidth(2); | |
728 | ||
729 | TLegend* legend = new TLegend(0.7, 0.8, 0.95, 0.95); | |
730 | legend->SetFillColor(0); | |
731 | legend->AddEntry(ratio1, "Gaussian"); | |
732 | legend->AddEntry(ratio2, "Simulation"); | |
733 | ||
734 | ratio2->SetTitle(""); | |
735 | ratio2->GetYaxis()->SetTitleOffset(1.5); | |
736 | ratio2->GetXaxis()->SetRangeUser(2, 4); | |
737 | ||
738 | TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500); | |
739 | InitPad(); | |
740 | ||
741 | ratio2->SetMarkerStyle(21); | |
742 | ratio1->SetMarkerStyle(22); | |
743 | ||
744 | ratio2->GetYaxis()->SetRangeUser(0.8, 1.2); | |
745 | ratio2->SetLineColor(kRed); | |
746 | ratio2->SetMarkerColor(kRed); | |
747 | ratio2->Draw("PL"); | |
748 | ratio1->Draw("SAMEPL"); | |
749 | ||
750 | legend->Draw(); | |
751 | ||
752 | canvas->SaveAs("Sigma2VertexCompare.eps"); | |
753 | } | |
754 | ||
755 | void drawSystematics() | |
756 | { | |
757 | //Secondaries(); | |
758 | //DrawDifferentSpecies(); | |
759 | //Composition(); | |
760 | ||
761 | Sigma2VertexSimulation(); | |
762 | } | |
763 | ||
764 | void DrawdNdEtaDifferences() | |
765 | { | |
766 | TH1* hists[5]; | |
767 | ||
768 | TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98); | |
769 | legend->SetFillColor(0); | |
770 | ||
771 | TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500); | |
772 | canvas->Divide(2, 1); | |
773 | ||
774 | canvas->cd(1); | |
775 | ||
776 | for (Int_t i=0; i<5; ++i) | |
777 | { | |
778 | hists[i] = 0; | |
779 | TFile* file = 0; | |
780 | TString title; | |
781 | ||
782 | switch(i) | |
783 | { | |
784 | case 0 : file = TFile::Open("systematics_dndeta_reference.root"); title = "standard composition"; break; | |
785 | case 1 : file = TFile::Open("systematics_dndeta_KBoosted.root"); title = "+ 50% kaons"; break; | |
786 | case 2 : file = TFile::Open("systematics_dndeta_KReduced.root"); title = "- 50% kaons"; break; | |
787 | case 3 : file = TFile::Open("systematics_dndeta_pBoosted.root"); title = "+ 50% protons"; break; | |
788 | case 4 : file = TFile::Open("systematics_dndeta_pReduced.root"); title = "- 50% protons"; break; | |
789 | default: return; | |
790 | } | |
791 | ||
792 | if (file) | |
793 | { | |
794 | hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2"); | |
795 | hists[i]->SetTitle("a)"); | |
796 | ||
797 | Prepare1DPlot(hists[i], kFALSE); | |
798 | hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999); | |
799 | hists[i]->SetLineColor(i+1); | |
800 | hists[i]->SetMarkerColor(i+1); | |
801 | hists[i]->GetXaxis()->SetLabelOffset(0.015); | |
802 | hists[i]->GetYaxis()->SetTitleOffset(1.5); | |
803 | gPad->SetLeftMargin(0.12); | |
804 | hists[i]->DrawCopy(((i > 0) ? "SAME" : "")); | |
805 | ||
806 | legend->AddEntry(hists[i], title); | |
807 | hists[i]->SetTitle(title); | |
808 | } | |
809 | } | |
810 | legend->Draw(); | |
811 | ||
812 | canvas->cd(2); | |
813 | gPad->SetLeftMargin(0.14); | |
814 | ||
815 | TLegend* legend2 = new TLegend(0.73, 0.73, 0.98, 0.98); | |
816 | legend2->SetFillColor(0); | |
817 | ||
818 | for (Int_t i=1; i<5; ++i) | |
819 | { | |
820 | if (hists[i]) | |
821 | { | |
822 | legend2->AddEntry(hists[i]); | |
823 | ||
824 | hists[i]->Divide(hists[0]); | |
825 | hists[i]->SetTitle("b)"); | |
826 | hists[i]->GetYaxis()->SetRangeUser(0.98, 1.02); | |
827 | hists[i]->GetYaxis()->SetTitle("Ratio to standard composition"); | |
828 | hists[i]->GetYaxis()->SetTitleOffset(1.8); | |
829 | hists[i]->DrawCopy(((i > 1) ? "SAME" : "")); | |
830 | } | |
831 | } | |
832 | ||
833 | legend2->Draw(); | |
834 | ||
835 | canvas->SaveAs("particlecomposition_result.eps"); | |
836 | canvas->SaveAs("particlecomposition_result.gif"); | |
837 | } |