updating makefile to include evgen headers to compile systematics selector
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / drawSystematics.C
CommitLineData
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
20extern TPad* gPad;
21
22void Track2Particle1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 10);
10ebe68d 23
24#endif
25
9f469bf5 26
10ebe68d 27void 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
39void SetRanges(TH1* hist)
40{
41 SetRanges(hist->GetXaxis());
42 SetRanges(hist->GetYaxis());
43 SetRanges(hist->GetZaxis());
44}
45
46void 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
55void Prepare2DPlot(TH2* hist)
56{
57 hist->SetStats(kFALSE);
58 hist->GetYaxis()->SetTitleOffset(1.4);
59
60 SetRanges(hist);
61}
62
7af955da 63void Prepare1DPlot(TH1* hist, Bool_t setRanges = kTRUE)
10ebe68d 64{
65 hist->SetLineWidth(2);
66 hist->SetStats(kFALSE);
67
7af955da 68 hist->GetXaxis()->SetTitleOffset(1.2);
69 hist->GetYaxis()->SetTitleOffset(1.2);
70
71 if (setRanges)
72 SetRanges(hist);
10ebe68d 73}
74
75void InitPad()
76{
9f469bf5 77 if (!gPad)
78 return;
79
10ebe68d 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
72e597d7 86 gPad->SetGridx();
87 gPad->SetGridy();
10ebe68d 88}
89
90void 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
100void Secondaries()
101{
102 TFile* file = TFile::Open("systematics.root");
103
7af955da 104 TH2F* secondaries = dynamic_cast<TH2F*> (file->Get("fSecondaries"));
10ebe68d 105 if (!secondaries)
106 {
107 printf("Could not read histogram\n");
108 return;
109 }
110
7af955da 111 TCanvas* canvas = new TCanvas("Secondaries", "Secondaries", 1000, 1000);
112 canvas->Divide(3, 3);
113 for (Int_t i=1; i<=8; i++)
10ebe68d 114 {
7af955da 115 TH1D* hist = secondaries->ProjectionY(Form("proj_%d", i), i, i);
116 hist->SetTitle(secondaries->GetXaxis()->GetBinLabel(i));
10ebe68d 117
7af955da 118 canvas->cd(i);
119 hist->Draw();
10ebe68d 120 }
121}
122
6b7fa615 123void Track2Particle1DComposition(const char** fileNames, Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9)
10ebe68d 124{
125 gSystem->Load("libPWG0base");
126
9f469bf5 127 TString canvasName;
128 canvasName.Form("Track2Particle1DComposition");
129 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
130 canvas->Divide(3, 1);
10ebe68d 131
9f469bf5 132 TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);
10ebe68d 133
9f469bf5 134 for (Int_t i=0; i<folderCount; ++i)
10ebe68d 135 {
6b7fa615 136 Track2Particle1DCreatePlots(fileNames[i], folderNames[i], upperPtLimit);
9f469bf5 137
6b7fa615 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])));
9f469bf5 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]);
10ebe68d 170 }
171
9f469bf5 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
6b7fa615 178TH1** 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
7af955da 192 TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
6b7fa615 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);
7af955da 233 ptDists[i]->SetStats(kFALSE);
234 ptDists[i]->SetTitle(";p_{T};Fraction of total");
6b7fa615 235 ptDists[i]->GetYaxis()->SetRangeUser(0, 1);
236 ptDists[i]->Draw((i == 0) ? "" : "SAME");
237 }
7af955da 238 legend->SetFillColor(0);
6b7fa615 239 legend->Draw();
240
241 canvas->SaveAs("DrawRatios.gif");
242
7af955da 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
6b7fa615 255 file->Close();
256
257 return ptDists;
258}
259
7af955da 260void 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
9f469bf5 270void DrawDifferentSpecies()
271{
272 gROOT->ProcessLine(".L drawPlots.C");
273
6b7fa615 274 const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "systematics.root" };
9f469bf5 275 const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
276
6b7fa615 277 Track2Particle1DComposition(fileNames, 4, folderNames);
9f469bf5 278}
10ebe68d 279
6b7fa615 280void 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
290void 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
300void 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
7af955da 310TH2F* HistogramDifferences(const char* filename, const char* folder1, const char* folder2)
6b7fa615 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
6b7fa615 324 TH3F* hist1 = fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
325 TH3F* hist2 = fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
326
7af955da 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
6b7fa615 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)
7af955da 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);
6b7fa615 337
338 printf("Over-/Underflow bins: %d %d\n", difference->GetBinContent(0), difference->GetBinContent(difference->GetNbinsX()+1));
339
340 return difference;
341}
342
343void HistogramDifferences()
344{
7af955da 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");
6b7fa615 349
7af955da 350 TCanvas* canvas = new TCanvas("HistogramDifferences", "HistogramDifferences", 1000, 1000);
351 canvas->Divide(2, 2);
6b7fa615 352
353 canvas->cd(1);
7af955da 354 KBoosted->GetXaxis()->SetRangeUser(-0.05, 0.05);
355 KBoosted->Draw("COLZ");
6b7fa615 356
357 canvas->cd(2);
7af955da 358 KReduced->GetXaxis()->SetRangeUser(-0.05, 0.05);
359 KReduced->Draw("COLZ");
6b7fa615 360
361 canvas->cd(3);
7af955da 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");
6b7fa615 368
369 canvas->SaveAs("HistogramDifferences.gif");
7af955da 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");
6b7fa615 417}
418
419
9f469bf5 420void ScalePtDependent(TH3F* hist, TF1* function)
421{
422 // assumes that pt is the third dimension of hist
423 // scales with function(pt)
10ebe68d 424
9f469bf5 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
6b7fa615 436void 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
452const char* ChangeComposition(void** correctionPointer, Int_t index)
9f469bf5 453{
454 AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer;
455
456 switch (index)
457 {
6b7fa615 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
7af955da 474 /*TH1** ptDists = DrawRatios("pythiaratios.root");
6b7fa615 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]);
7af955da 480 }*/
6b7fa615 481 return "PythiaRatios";
9f469bf5 482 break;
483
7af955da 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
6b7fa615 502 // each species rated with pythia ratios
7af955da 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
6b7fa615 509 TH1** ptDists = DrawRatios("pythiaratios.root");
510 Int_t functionIndex = (index - 2) / 2;
7af955da 511 Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
6b7fa615 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();
9f469bf5 522 break;
523
6b7fa615 524 case 999:
9f469bf5 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 }
6b7fa615 531
532 return "noname";
9f469bf5 533}
534
535void Composition()
536{
537 gSystem->Load("libPWG0base");
538
539 gSystem->Unlink("new_compositions.root");
540
6b7fa615 541 Int_t nCompositions = 8;
542 for (Int_t comp = 0; comp < nCompositions; ++comp)
9f469bf5 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
6b7fa615 556 const char* newName = ChangeComposition(fdNdEtaCorrection, comp);
9f469bf5 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 }
10ebe68d 573
9f469bf5 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]);
10ebe68d 575
9f469bf5 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]);
10ebe68d 577
9f469bf5 578 TList* collection = new TList;
10ebe68d 579
9f469bf5 580 for (Int_t i=0; i<4; ++i)
581 collection->Add(fdNdEtaCorrection[i]);
10ebe68d 582
6b7fa615 583 AlidNdEtaCorrection* newComposition = new AlidNdEtaCorrection(newName, newName);
9f469bf5 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 }
10ebe68d 594
595 gROOT->ProcessLine(".L drawPlots.C");
9f469bf5 596
6b7fa615 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" };
9f469bf5 599
6b7fa615 600 Track2Particle1DComposition(fileNames, nCompositions, folderNames);
9f469bf5 601}
602
603Double_t ConvSigma1To2D(Double_t sigma)
604{
605 return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
606}
607
608Double_t ConvDistance1DTo2D(Double_t distance)
609{
610 return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
611}
612
613Double_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
636TH2F* 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
4c6b34a8 653TH1F* Sigma2VertexGaussian()
9f469bf5 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
72e597d7 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)
9f469bf5 665 ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
666 ratio->SetMarkerStyle(21);
667
668 canvas->cd(2);
4c6b34a8 669 ratio->DrawCopy("P");
9f469bf5 670
72e597d7 671 TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 3 sigma / % included n sigma", 50, 0.05, 5.05);
9f469bf5 672 Double_t sigma3 = Sigma2VertexCount(tracks, 3);
72e597d7 673 for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
674 ratio2->Fill(nSigma, sigma3 / ratio->GetBinContent(ratio->FindBin(nSigma)));
9f469bf5 675 ratio2->SetMarkerStyle(21);
676
677 canvas->cd(3);
4c6b34a8 678 ratio2->DrawCopy("P");
9f469bf5 679
680 canvas->SaveAs("Sigma2Vertex.eps");
4c6b34a8 681
682 return ratio2;
9f469bf5 683}
684
4c6b34a8 685TH1F* Sigma2VertexSimulation()
9f469bf5 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
5c495d37 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());
9f469bf5 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);
4c6b34a8 710 ratio->DrawCopy("P");
711
712 return ratio;
10ebe68d 713}
714
4c6b34a8 715void 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
0bd1f8a0 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);
4c6b34a8 731 legend->AddEntry(ratio1, "Gaussian");
732 legend->AddEntry(ratio2, "Simulation");
733
0bd1f8a0 734 ratio2->SetTitle("");
735 ratio2->GetYaxis()->SetTitleOffset(1.5);
736 ratio2->GetXaxis()->SetRangeUser(2, 4);
4c6b34a8 737
738 TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
739 InitPad();
740
5c495d37 741 ratio2->SetMarkerStyle(21);
742 ratio1->SetMarkerStyle(22);
743
744 ratio2->GetYaxis()->SetRangeUser(0.8, 1.2);
4c6b34a8 745 ratio2->SetLineColor(kRed);
5c495d37 746 ratio2->SetMarkerColor(kRed);
747 ratio2->Draw("PL");
748 ratio1->Draw("SAMEPL");
4c6b34a8 749
750 legend->Draw();
0bd1f8a0 751
752 canvas->SaveAs("Sigma2VertexCompare.eps");
4c6b34a8 753}
9f469bf5 754
10ebe68d 755void drawSystematics()
756{
9f469bf5 757 //Secondaries();
758 //DrawDifferentSpecies();
759 //Composition();
760
761 Sigma2VertexSimulation();
10ebe68d 762}
7af955da 763
764void DrawdNdEtaDifferences()
765{
766 TH1* hists[5];
767
72e597d7 768 TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
f1076daa 769 legend->SetFillColor(0);
770
7af955da 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 {
72e597d7 778 hists[i] = 0;
7af955da 779 TFile* file = 0;
f1076daa 780 TString title;
7af955da 781
782 switch(i)
783 {
f1076daa 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;
7af955da 789 default: return;
790 }
791
72e597d7 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 }
7af955da 809 }
f1076daa 810 legend->Draw();
7af955da 811
812 canvas->cd(2);
f1076daa 813 gPad->SetLeftMargin(0.14);
814
815 TLegend* legend2 = new TLegend(0.73, 0.73, 0.98, 0.98);
816 legend2->SetFillColor(0);
7af955da 817
818 for (Int_t i=1; i<5; ++i)
819 {
72e597d7 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 }
7af955da 831 }
f1076daa 832
833 legend2->Draw();
834
835 canvas->SaveAs("particlecomposition_result.eps");
836 canvas->SaveAs("particlecomposition_result.gif");
7af955da 837}