updating makefile to include evgen headers to compile systematics selector
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / drawSystematics.C
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 }