]>
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> | |
6b7fa615 | 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); | |
10ebe68d | 23 | |
24 | #endif | |
25 | ||
9f469bf5 | 26 | |
10ebe68d | 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) | |
64 | { | |
65 | hist->SetLineWidth(2); | |
66 | hist->SetStats(kFALSE); | |
67 | ||
68 | SetRanges(hist); | |
69 | } | |
70 | ||
71 | void InitPad() | |
72 | { | |
9f469bf5 | 73 | if (!gPad) |
74 | return; | |
75 | ||
10ebe68d | 76 | gPad->Range(0, 0, 1, 1); |
77 | gPad->SetLeftMargin(0.15); | |
78 | //gPad->SetRightMargin(0.05); | |
79 | //gPad->SetTopMargin(0.13); | |
80 | //gPad->SetBottomMargin(0.1); | |
81 | ||
82 | //gPad->SetGridx(); | |
83 | //gPad->SetGridy(); | |
84 | } | |
85 | ||
86 | void InitPadCOLZ() | |
87 | { | |
88 | gPad->Range(0, 0, 1, 1); | |
89 | gPad->SetRightMargin(0.15); | |
90 | gPad->SetLeftMargin(0.12); | |
91 | ||
92 | gPad->SetGridx(); | |
93 | gPad->SetGridy(); | |
94 | } | |
95 | ||
96 | void Secondaries() | |
97 | { | |
98 | TFile* file = TFile::Open("systematics.root"); | |
99 | ||
100 | TH3F* secondaries = dynamic_cast<TH3F*> (file->Get("fSecondaries")); | |
101 | if (!secondaries) | |
102 | { | |
103 | printf("Could not read histogram\n"); | |
104 | return; | |
105 | } | |
106 | ||
107 | for (Int_t ptBin=1; ptBin<=secondaries->GetNbinsZ(); ptBin++) | |
108 | //for (Int_t ptBin = 1; ptBin<=2; ptBin++) | |
109 | { | |
9f469bf5 | 110 | TH1F* hist = new TH1F(Form("secondaries_%d", ptBin), Form("secondaries_%d", ptBin), secondaries->GetNbinsY(), secondaries->GetYaxis()->GetXmin(), secondaries->GetYaxis()->GetXmax()); |
111 | ||
112 | hist->SetTitle(Form("%f < p_{T} < %f", secondaries->GetZaxis()->GetBinLowEdge(ptBin), secondaries->GetZaxis()->GetBinUpEdge(ptBin))); | |
10ebe68d | 113 | |
114 | for (Int_t cBin=1; cBin<=secondaries->GetNbinsY(); ++cBin) | |
115 | { | |
116 | if (secondaries->GetBinContent(0, cBin, ptBin) > 0) | |
117 | printf("WARNING: Underflow bin not empty!"); | |
118 | if (secondaries->GetBinContent(secondaries->GetNbinsX()+1, cBin, ptBin) > 0) | |
119 | printf("WARNING: Overflow bin not empty!"); | |
120 | ||
121 | Double_t sum = 0; | |
122 | Double_t count = 0; | |
123 | for (Int_t nBin=1; nBin<=secondaries->GetNbinsX(); ++nBin) | |
124 | { | |
125 | //printf("%f %f\n", secondaries->GetXaxis()->GetBinCenter(nBin), secondaries->GetBinContent(nBin, cBin, ptBin)); | |
126 | sum += secondaries->GetXaxis()->GetBinCenter(nBin) * secondaries->GetBinContent(nBin, cBin, ptBin); | |
127 | count += secondaries->GetBinContent(nBin, cBin, ptBin); | |
128 | } | |
129 | ||
130 | printf("%f %f\n", sum, count); | |
131 | ||
132 | if (count > 0) | |
9f469bf5 | 133 | { |
134 | hist->SetBinContent(cBin, sum / count); | |
135 | hist->SetBinError(cBin, TMath::Sqrt(sum) / count); | |
136 | } | |
10ebe68d | 137 | } |
138 | ||
9f469bf5 | 139 | hist->Scale(1.0 / hist->GetBinContent(hist->GetXaxis()->FindBin(1))); |
140 | hist->Add(new TF1("flat", "-1", 0, 2)); | |
141 | ||
10ebe68d | 142 | new TCanvas; |
9f469bf5 | 143 | hist->SetMarkerStyle(21); |
144 | hist->Draw(""); | |
10ebe68d | 145 | } |
146 | } | |
147 | ||
6b7fa615 | 148 | void Track2Particle1DComposition(const char** fileNames, Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9) |
10ebe68d | 149 | { |
150 | gSystem->Load("libPWG0base"); | |
151 | ||
9f469bf5 | 152 | TString canvasName; |
153 | canvasName.Form("Track2Particle1DComposition"); | |
154 | TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400); | |
155 | canvas->Divide(3, 1); | |
10ebe68d | 156 | |
9f469bf5 | 157 | TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95); |
10ebe68d | 158 | |
9f469bf5 | 159 | for (Int_t i=0; i<folderCount; ++i) |
10ebe68d | 160 | { |
6b7fa615 | 161 | Track2Particle1DCreatePlots(fileNames[i], folderNames[i], upperPtLimit); |
9f469bf5 | 162 | |
6b7fa615 | 163 | TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", folderNames[i], folderNames[i]))); |
164 | TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", folderNames[i], folderNames[i]))); | |
165 | TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", folderNames[i], folderNames[i]))); | |
9f469bf5 | 166 | |
167 | Prepare1DPlot(corrX); | |
168 | Prepare1DPlot(corrY); | |
169 | Prepare1DPlot(corrZ); | |
170 | ||
171 | const char* title = "Track2Particle Correction"; | |
172 | corrX->SetTitle(title); | |
173 | corrY->SetTitle(title); | |
174 | corrZ->SetTitle(title); | |
175 | ||
176 | corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit); | |
177 | ||
178 | corrX->SetLineColor(i+1); | |
179 | corrY->SetLineColor(i+1); | |
180 | corrZ->SetLineColor(i+1); | |
181 | ||
182 | canvas->cd(1); | |
183 | InitPad(); | |
184 | corrX->DrawCopy(((i>0) ? "SAME" : "")); | |
185 | ||
186 | canvas->cd(2); | |
187 | InitPad(); | |
188 | corrY->DrawCopy(((i>0) ? "SAME" : "")); | |
189 | ||
190 | canvas->cd(3); | |
191 | InitPad(); | |
192 | corrZ->DrawCopy(((i>0) ? "SAME" : "")); | |
193 | ||
194 | legend->AddEntry(corrZ, folderNames[i]); | |
10ebe68d | 195 | } |
196 | ||
9f469bf5 | 197 | legend->Draw(); |
198 | ||
199 | //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.gif", fileName, gMax, upperPtLimit)); | |
200 | //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.eps", fileName, gMax, upperPtLimit)); | |
201 | } | |
202 | ||
6b7fa615 | 203 | TH1** DrawRatios(const char* fileName = "systematics.root") |
204 | { | |
205 | gSystem->Load("libPWG0base"); | |
206 | ||
207 | TFile* file = TFile::Open(fileName); | |
208 | ||
209 | TString canvasName; | |
210 | canvasName.Form("DrawRatios"); | |
211 | TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400); | |
212 | canvas->Divide(2, 1); | |
213 | canvas->cd(1); | |
214 | ||
215 | TH1** ptDists = new TH1*[4]; | |
216 | ||
217 | TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95); | |
218 | ||
219 | const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" }; | |
220 | const char* particleNames[] = { "#pi", "K", "p", "other" }; | |
221 | for (Int_t i=0; i<4; ++i) | |
222 | { | |
223 | AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderNames[i], folderNames[i]); | |
224 | dNdEtaCorrection->LoadHistograms(fileName, folderNames[i]); | |
225 | ||
226 | TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram(); | |
227 | ||
228 | gene->GetYaxis()->SetRangeUser(-0.8, 0.8); | |
229 | gene->GetXaxis()->SetRangeUser(-10, 10); | |
230 | ||
231 | ptDists[i] = dynamic_cast<TH1*> (gene->Project3D("z")->Clone(Form("%s_z", folderNames[i]))); | |
232 | ptDists[i]->SetTitle(";p_{T};Count"); | |
233 | if (!ptDists[i]) | |
234 | { | |
235 | printf("Problem getting distribution %d.\n", i); | |
236 | return 0; | |
237 | } | |
238 | ||
239 | ptDists[i]->GetYaxis()->SetRangeUser(1, ptDists[i]->GetMaximum()*1.1); | |
240 | ptDists[i]->GetXaxis()->SetRangeUser(0, 9.9); | |
241 | ptDists[i]->SetLineColor(i+1); | |
242 | ptDists[i]->DrawCopy((i == 0) ? "" : "SAME"); | |
243 | ptDists[i]->GetYaxis()->SetRange(0, 0); | |
244 | ||
245 | legend->AddEntry(ptDists[i], particleNames[i]); | |
246 | } | |
247 | gPad->SetLogy(); | |
248 | ||
249 | TH1* total = dynamic_cast<TH1*> (ptDists[0]->Clone("total")); | |
250 | ||
251 | for (Int_t i=1; i<4; ++i) | |
252 | total->Add(ptDists[i]); | |
253 | ||
254 | canvas->cd(2); | |
255 | for (Int_t i=0; i<4; ++i) | |
256 | { | |
257 | ptDists[i]->Divide(total); | |
258 | ptDists[i]->SetTitle(";p_{T};Ratio"); | |
259 | ptDists[i]->GetYaxis()->SetRangeUser(0, 1); | |
260 | ptDists[i]->Draw((i == 0) ? "" : "SAME"); | |
261 | } | |
262 | legend->Draw(); | |
263 | ||
264 | canvas->SaveAs("DrawRatios.gif"); | |
265 | ||
266 | file->Close(); | |
267 | ||
268 | return ptDists; | |
269 | } | |
270 | ||
9f469bf5 | 271 | void DrawDifferentSpecies() |
272 | { | |
273 | gROOT->ProcessLine(".L drawPlots.C"); | |
274 | ||
6b7fa615 | 275 | const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "systematics.root" }; |
9f469bf5 | 276 | const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" }; |
277 | ||
6b7fa615 | 278 | Track2Particle1DComposition(fileNames, 4, folderNames); |
9f469bf5 | 279 | } |
10ebe68d | 280 | |
6b7fa615 | 281 | void DrawSpeciesAndCombination() |
282 | { | |
283 | gROOT->ProcessLine(".L drawPlots.C"); | |
284 | ||
285 | const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "new_compositions.root" }; | |
286 | const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "PythiaRatios" }; | |
287 | ||
288 | Track2Particle1DComposition(fileNames, 4, folderNames); | |
289 | } | |
290 | ||
291 | void DrawSimulatedVsCombined() | |
292 | { | |
293 | gROOT->ProcessLine(".L drawPlots.C"); | |
294 | ||
295 | const char* fileNames[] = { "new_compositions.root", "new_compositions.root" }; | |
296 | const char* folderNames[] = { "Pythia", "PythiaRatios" }; | |
297 | ||
298 | Track2Particle1DComposition(fileNames, 2, folderNames); | |
299 | } | |
300 | ||
301 | void DrawBoosts() | |
302 | { | |
303 | gROOT->ProcessLine(".L drawPlots.C"); | |
304 | ||
305 | const char* fileNames[] = { "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root" }; | |
306 | const char* folderNames[] = { "PythiaRatios", "PiBoosted", "KBoosted", "pBoosted" }; | |
307 | ||
308 | Track2Particle1DComposition(fileNames, 4, folderNames); | |
309 | } | |
310 | ||
311 | TH1F* HistogramDifferences(const char* filename, const char* folder1, const char* folder2) | |
312 | { | |
313 | gSystem->Load("libPWG0base"); | |
314 | ||
315 | AlidNdEtaCorrection* fdNdEtaCorrection[2]; | |
316 | ||
317 | TFile::Open(filename); | |
318 | ||
319 | fdNdEtaCorrection[0] = new AlidNdEtaCorrection(folder1, folder1); | |
320 | fdNdEtaCorrection[0]->LoadHistograms(filename, folder1); | |
321 | ||
322 | fdNdEtaCorrection[1] = new AlidNdEtaCorrection(folder2, folder2); | |
323 | fdNdEtaCorrection[1]->LoadHistograms(filename, folder2); | |
324 | ||
325 | TH1F* difference = new TH1F("difference", Form(";#DeltaC_{pT, z, #eta} %s - %s;Count", folder2, folder1), 1000, -0.5, 0.5); | |
326 | ||
327 | TH3F* hist1 = fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram(); | |
328 | TH3F* hist2 = fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetCorrectionHistogram(); | |
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 | difference->Fill(hist2->GetBinContent(x, y, z) - hist1->GetBinContent(x, y, z)); | |
334 | ||
335 | printf("Over-/Underflow bins: %d %d\n", difference->GetBinContent(0), difference->GetBinContent(difference->GetNbinsX()+1)); | |
336 | ||
337 | return difference; | |
338 | } | |
339 | ||
340 | void HistogramDifferences() | |
341 | { | |
342 | TH1F* PiBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "PiBoosted"); | |
343 | TH1F* KBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "KBoosted"); | |
344 | TH1F* pBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "pBoosted"); | |
345 | ||
346 | TCanvas* canvas = new TCanvas("HistogramDifferences", "HistogramDifferences", 1200, 400); | |
347 | canvas->Divide(3, 1); | |
348 | ||
349 | canvas->cd(1); | |
350 | PiBoosted->GetXaxis()->SetRangeUser(-0.01, 0.01); | |
351 | PiBoosted->Draw(); | |
352 | ||
353 | canvas->cd(2); | |
354 | KBoosted->GetXaxis()->SetRangeUser(-0.01, 0.01); | |
355 | KBoosted->Draw(); | |
356 | ||
357 | canvas->cd(3); | |
358 | pBoosted->GetXaxis()->SetRangeUser(-0.01, 0.01); | |
359 | pBoosted->Draw(); | |
360 | ||
361 | canvas->SaveAs("HistogramDifferences.gif"); | |
362 | } | |
363 | ||
364 | ||
9f469bf5 | 365 | void ScalePtDependent(TH3F* hist, TF1* function) |
366 | { | |
367 | // assumes that pt is the third dimension of hist | |
368 | // scales with function(pt) | |
10ebe68d | 369 | |
9f469bf5 | 370 | for (Int_t z=1; z<=hist->GetNbinsZ(); ++z) |
371 | { | |
372 | Double_t factor = function->Eval(hist->GetZaxis()->GetBinCenter(z)); | |
373 | printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor); | |
374 | ||
375 | for (Int_t x=1; x<=hist->GetNbinsX(); ++x) | |
376 | for (Int_t y=1; y<=hist->GetNbinsY(); ++y) | |
377 | hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor); | |
378 | } | |
379 | } | |
380 | ||
6b7fa615 | 381 | void ScalePtDependent(TH3F* hist, TH1* function) |
382 | { | |
383 | // assumes that pt is the third dimension of hist | |
384 | // scales with histogram(pt) | |
385 | ||
386 | for (Int_t z=1; z<=hist->GetNbinsZ(); ++z) | |
387 | { | |
388 | Double_t factor = function->GetBinContent(function->GetXaxis()->FindBin(hist->GetZaxis()->GetBinCenter(z))); | |
389 | printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor); | |
390 | ||
391 | for (Int_t x=1; x<=hist->GetNbinsX(); ++x) | |
392 | for (Int_t y=1; y<=hist->GetNbinsY(); ++y) | |
393 | hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor); | |
394 | } | |
395 | } | |
396 | ||
397 | const char* ChangeComposition(void** correctionPointer, Int_t index) | |
9f469bf5 | 398 | { |
399 | AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer; | |
400 | ||
401 | switch (index) | |
402 | { | |
6b7fa615 | 403 | case 0: // result from pp events |
404 | { | |
405 | TFile::Open("pythiaratios.root"); | |
406 | ||
407 | for (Int_t i=0; i<4; ++i) | |
408 | { | |
409 | TString name; | |
410 | name.Form("correction_%d", i); | |
411 | fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name); | |
412 | fdNdEtaCorrection[i]->LoadHistograms("pythiaratios.root", name); | |
413 | } | |
414 | } | |
415 | return "Pythia"; | |
416 | break; | |
417 | ||
418 | case 1: // each species rated with pythia ratios | |
419 | TH1** ptDists = DrawRatios("pythiaratios.root"); | |
420 | ||
421 | for (Int_t i=0; i<3; ++i) | |
422 | { | |
423 | ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]); | |
424 | ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]); | |
425 | } | |
426 | return "PythiaRatios"; | |
9f469bf5 | 427 | break; |
428 | ||
6b7fa615 | 429 | // each species rated with pythia ratios |
430 | case 2: // + 10% pions | |
431 | case 3: // - 10% pions | |
432 | case 4: // + 10% kaons | |
433 | case 5: // - 10% kaons | |
434 | case 6: // + 10% protons | |
435 | case 7: // - 10% protons | |
436 | TH1** ptDists = DrawRatios("pythiaratios.root"); | |
437 | Int_t functionIndex = (index - 2) / 2; | |
438 | Double_t scaleFactor = (index % 2 == 0) ? 1.1 : 0.9; | |
439 | ptDists[functionIndex]->Scale(scaleFactor); | |
440 | ||
441 | for (Int_t i=0; i<3; ++i) | |
442 | { | |
443 | ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]); | |
444 | ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]); | |
445 | } | |
446 | TString* str = new TString; | |
447 | str->Form("%s%s", (functionIndex == 0) ? "Pi" : ((functionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced"); | |
448 | return str->Data(); | |
9f469bf5 | 449 | break; |
450 | ||
6b7fa615 | 451 | case 999: |
9f469bf5 | 452 | TF1* ptDependence = new TF1("simple", "x", 0, 100); |
453 | ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDependence); | |
454 | ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDependence); | |
455 | break; | |
456 | ||
457 | } | |
6b7fa615 | 458 | |
459 | return "noname"; | |
9f469bf5 | 460 | } |
461 | ||
462 | void Composition() | |
463 | { | |
464 | gSystem->Load("libPWG0base"); | |
465 | ||
466 | gSystem->Unlink("new_compositions.root"); | |
467 | ||
6b7fa615 | 468 | Int_t nCompositions = 8; |
469 | for (Int_t comp = 0; comp < nCompositions; ++comp) | |
9f469bf5 | 470 | { |
471 | AlidNdEtaCorrection* fdNdEtaCorrection[4]; | |
472 | ||
473 | TFile::Open("systematics.root"); | |
474 | ||
475 | for (Int_t i=0; i<4; ++i) | |
476 | { | |
477 | TString name; | |
478 | name.Form("correction_%d", i); | |
479 | fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name); | |
480 | fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name); | |
481 | } | |
482 | ||
6b7fa615 | 483 | const char* newName = ChangeComposition(fdNdEtaCorrection, comp); |
9f469bf5 | 484 | |
485 | Double_t geneCount[5]; | |
486 | Double_t measCount[5]; | |
487 | geneCount[4] = 0; | |
488 | measCount[4] = 0; | |
489 | ||
490 | for (Int_t i=0; i<4; ++i) | |
491 | { | |
492 | geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Integral(); | |
493 | measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Integral(); | |
494 | ||
495 | geneCount[4] += geneCount[i]; | |
496 | measCount[4] += measCount[i]; | |
497 | ||
498 | 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]); | |
499 | } | |
10ebe68d | 500 | |
9f469bf5 | 501 | 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 | 502 | |
9f469bf5 | 503 | 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 | 504 | |
9f469bf5 | 505 | TList* collection = new TList; |
10ebe68d | 506 | |
9f469bf5 | 507 | for (Int_t i=0; i<4; ++i) |
508 | collection->Add(fdNdEtaCorrection[i]); | |
10ebe68d | 509 | |
6b7fa615 | 510 | AlidNdEtaCorrection* newComposition = new AlidNdEtaCorrection(newName, newName); |
9f469bf5 | 511 | newComposition->Merge(collection); |
512 | newComposition->Finish(); | |
513 | ||
514 | delete collection; | |
515 | ||
516 | TFile* file = TFile::Open("new_compositions.root", "UPDATE"); | |
517 | newComposition->SaveHistograms(); | |
518 | //file->Write(); | |
519 | file->Close(); | |
520 | } | |
10ebe68d | 521 | |
522 | gROOT->ProcessLine(".L drawPlots.C"); | |
9f469bf5 | 523 | |
6b7fa615 | 524 | 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" }; |
525 | const char* folderNames[] = { "Pythia", "PythiaRatios", "PiBoosted", "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced" }; | |
9f469bf5 | 526 | |
6b7fa615 | 527 | Track2Particle1DComposition(fileNames, nCompositions, folderNames); |
9f469bf5 | 528 | } |
529 | ||
530 | Double_t ConvSigma1To2D(Double_t sigma) | |
531 | { | |
532 | return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2); | |
533 | } | |
534 | ||
535 | Double_t ConvDistance1DTo2D(Double_t distance) | |
536 | { | |
537 | return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2); | |
538 | } | |
539 | ||
540 | Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma) | |
541 | { | |
542 | Double_t count = 0; | |
543 | ||
544 | //nSigma = ConvSigma1To2D(nSigma); | |
545 | ||
546 | for (Int_t x=1; x<=tracks->GetNbinsX(); ++x) | |
547 | for (Int_t y=1; y<=tracks->GetNbinsY(); ++y) | |
548 | { | |
549 | Double_t impactX = tracks->GetXaxis()->GetBinCenter(x); | |
550 | Double_t impactY = tracks->GetYaxis()->GetBinCenter(y); | |
551 | ||
552 | Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY); | |
553 | ||
554 | d = ConvDistance1DTo2D(d); | |
555 | ||
556 | if (d < nSigma) | |
557 | count += tracks->GetBinContent(x, y); | |
558 | } | |
559 | ||
560 | return count; | |
561 | } | |
562 | ||
563 | TH2F* Sigma2VertexGaussianTracksHist() | |
564 | { | |
565 | TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5); | |
566 | ||
567 | TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5); | |
568 | gaussian2D->SetParameters(1, 0, 1, 1, 0, 1); | |
569 | ||
570 | for (Int_t x=1; x<=tracks->GetNbinsX(); ++x) | |
571 | for (Int_t y=1; y<=tracks->GetNbinsY(); ++y) | |
572 | tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y))); | |
573 | ||
574 | //normalize | |
575 | tracks->Scale(1.0 / tracks->Integral()); | |
576 | ||
577 | return tracks; | |
578 | } | |
579 | ||
4c6b34a8 | 580 | TH1F* Sigma2VertexGaussian() |
9f469bf5 | 581 | { |
582 | TH2F* tracks = Sigma2VertexGaussianTracksHist(); | |
583 | ||
584 | TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000); | |
585 | canvas->Divide(2, 2); | |
586 | ||
587 | canvas->cd(1); | |
588 | tracks->Draw("COLZ"); | |
589 | ||
590 | TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 10, 0.25, 5.25); | |
591 | for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5) | |
592 | ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma)); | |
593 | ratio->SetMarkerStyle(21); | |
594 | ||
595 | canvas->cd(2); | |
4c6b34a8 | 596 | ratio->DrawCopy("P"); |
9f469bf5 | 597 | |
598 | TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 3 sigma / % included n sigma", 10, 0.25, 5.25); | |
599 | Double_t sigma3 = Sigma2VertexCount(tracks, 3); | |
600 | for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5) | |
601 | ratio2->Fill(nSigma, sigma3 / Sigma2VertexCount(tracks, nSigma)); | |
602 | ratio2->SetMarkerStyle(21); | |
603 | ||
604 | canvas->cd(3); | |
4c6b34a8 | 605 | ratio2->DrawCopy("P"); |
9f469bf5 | 606 | |
607 | canvas->SaveAs("Sigma2Vertex.eps"); | |
4c6b34a8 | 608 | |
609 | return ratio2; | |
9f469bf5 | 610 | } |
611 | ||
4c6b34a8 | 612 | TH1F* Sigma2VertexSimulation() |
9f469bf5 | 613 | { |
614 | TFile* file = TFile::Open("systematics.root"); | |
615 | ||
616 | TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertex")); | |
617 | if (!sigmavertex) | |
618 | { | |
619 | printf("Could not read histogram\n"); | |
620 | return; | |
621 | } | |
622 | ||
4c6b34a8 | 623 | TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;Nsigma;% included 3 sigma / % included n sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax()); |
9f469bf5 | 624 | |
625 | for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i) | |
626 | ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(3)) / sigmavertex->GetBinContent(i)); | |
627 | ||
628 | TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500); | |
629 | canvas->Divide(2, 1); | |
630 | ||
631 | canvas->cd(1); | |
632 | sigmavertex->SetMarkerStyle(21); | |
633 | sigmavertex->Draw("P"); | |
634 | ||
635 | canvas->cd(2); | |
636 | ratio->SetMarkerStyle(21); | |
4c6b34a8 | 637 | ratio->DrawCopy("P"); |
638 | ||
639 | return ratio; | |
10ebe68d | 640 | } |
641 | ||
4c6b34a8 | 642 | void Sigma2VertexCompare() |
643 | { | |
644 | TH1F* ratio1 = Sigma2VertexGaussian(); | |
645 | TH1F* ratio2 = Sigma2VertexSimulation(); | |
646 | ||
647 | ratio1->SetStats(kFALSE); | |
648 | ratio2->SetStats(kFALSE); | |
649 | ||
650 | ratio1->SetMarkerStyle(0); | |
651 | ratio2->SetMarkerStyle(0); | |
652 | ||
653 | TLegend* legend = new TLegend(0.647177,0.775424,0.961694,0.966102); | |
654 | legend->AddEntry(ratio1, "Gaussian"); | |
655 | legend->AddEntry(ratio2, "Simulation"); | |
656 | ||
657 | ratio1->GetXaxis()->SetTitleOffset(1.5); | |
658 | ||
659 | TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500); | |
660 | InitPad(); | |
661 | ||
662 | ratio1->Draw(); | |
663 | ratio2->SetLineColor(kRed); | |
664 | ratio2->Draw("SAME"); | |
665 | ||
666 | legend->Draw(); | |
667 | } | |
9f469bf5 | 668 | |
10ebe68d | 669 | void drawSystematics() |
670 | { | |
9f469bf5 | 671 | //Secondaries(); |
672 | //DrawDifferentSpecies(); | |
673 | //Composition(); | |
674 | ||
675 | Sigma2VertexSimulation(); | |
10ebe68d | 676 | } |