30ecad6463a4fc9c7beb940df776e63a4ed09b88
[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 Int_t markers[] = {20,20,21,22,23,28,29};
27 Int_t colors[]  = {1,2,3,4,6,8,102};
28
29 void SetRanges(TAxis* axis)
30 {
31   if (strcmp(axis->GetTitle(), "#eta") == 0)
32     axis->SetRangeUser(-1.7999, 1.7999);
33   if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0)
34     axis->SetRangeUser(0, 9.9999);
35   if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0)
36     axis->SetRangeUser(-15, 14.9999);
37   if (strcmp(axis->GetTitle(), "Ntracks") == 0)
38     axis->SetRangeUser(0, 99.9999);
39 }
40
41 void SetRanges(TH1* hist)
42 {
43   SetRanges(hist->GetXaxis());
44   SetRanges(hist->GetYaxis());
45   SetRanges(hist->GetZaxis());
46 }
47
48 void Prepare3DPlot(TH3* hist)
49 {
50   hist->GetXaxis()->SetTitleOffset(1.5);
51   hist->GetYaxis()->SetTitleOffset(1.5);
52   hist->GetZaxis()->SetTitleOffset(1.5);
53
54   hist->SetStats(kFALSE);
55 }
56
57 void Prepare2DPlot(TH2* hist)
58 {
59   hist->SetStats(kFALSE);
60   hist->GetYaxis()->SetTitleOffset(1.4);
61
62   SetRanges(hist);
63 }
64
65 void Prepare1DPlot(TH1* hist, Bool_t setRanges = kTRUE)
66 {
67   hist->SetLineWidth(2);
68   hist->SetStats(kFALSE);
69
70   hist->GetXaxis()->SetTitleOffset(1.2);
71   hist->GetYaxis()->SetTitleOffset(1.2);
72
73   if (setRanges)
74     SetRanges(hist);
75 }
76
77 void InitPad()
78 {
79   if (!gPad)
80     return;
81
82   gPad->Range(0, 0, 1, 1);
83   gPad->SetLeftMargin(0.15);
84   //gPad->SetRightMargin(0.05);
85   //gPad->SetTopMargin(0.13);
86   //gPad->SetBottomMargin(0.1);
87
88   gPad->SetGridx();
89   gPad->SetGridy();
90 }
91
92 void InitPadCOLZ()
93 {
94   gPad->Range(0, 0, 1, 1);
95   gPad->SetRightMargin(0.15);
96   gPad->SetLeftMargin(0.12);
97
98   gPad->SetGridx();
99   gPad->SetGridy();
100 }
101
102 void Secondaries()
103 {
104   TFile* file = TFile::Open("systematics.root");
105
106   TH2F* secondaries = dynamic_cast<TH2F*> (file->Get("fSecondaries"));
107   if (!secondaries)
108   {
109     printf("Could not read histogram\n");
110     return;
111   }
112
113   TCanvas* canvas = new TCanvas("Secondaries", "Secondaries", 1000, 1000);
114   canvas->Divide(3, 3);
115   for (Int_t i=1; i<=8; i++)
116   {
117     TH1D* hist = secondaries->ProjectionY(Form("proj_%d", i), i, i);
118     hist->SetTitle(secondaries->GetXaxis()->GetBinLabel(i));
119
120     canvas->cd(i);
121     hist->Draw();
122   }
123 }
124
125 void Track2Particle1DComposition(const char** fileNames, Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9)
126 {
127   gSystem->Load("libPWG0base");
128
129   TString canvasName;
130   canvasName.Form("Track2Particle1DComposition");
131   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
132   canvas->Divide(3, 1);
133
134   TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);
135
136   for (Int_t i=0; i<folderCount; ++i)
137   {
138     Track2Particle1DCreatePlots(fileNames[i], folderNames[i], upperPtLimit);
139
140     TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", folderNames[i], folderNames[i])));
141     TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", folderNames[i], folderNames[i])));
142     TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", folderNames[i], folderNames[i])));
143
144     Prepare1DPlot(corrX);
145     Prepare1DPlot(corrY);
146     Prepare1DPlot(corrZ);
147
148     const char* title = "Track2Particle Correction";
149     corrX->SetTitle(title);
150     corrY->SetTitle(title);
151     corrZ->SetTitle(title);
152
153     corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
154
155     corrX->SetLineColor(i+1);
156     corrY->SetLineColor(i+1);
157     corrZ->SetLineColor(i+1);
158
159     canvas->cd(1);
160     InitPad();
161     corrX->DrawCopy(((i>0) ? "SAME" : ""));
162
163     canvas->cd(2);
164     InitPad();
165     corrY->DrawCopy(((i>0) ? "SAME" : ""));
166
167     canvas->cd(3);
168     InitPad();
169     corrZ->DrawCopy(((i>0) ? "SAME" : ""));
170
171     legend->AddEntry(corrZ, folderNames[i]);
172   }
173
174   legend->Draw();
175
176   //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.gif", fileName, gMax, upperPtLimit));
177   //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.eps", fileName, gMax, upperPtLimit));
178 }
179
180 TH1** DrawRatios(const char* fileName = "systematics.root")
181 {
182   gSystem->Load("libPWG0base");
183
184   TFile* file = TFile::Open(fileName);
185
186   TString canvasName;
187   canvasName.Form("DrawRatios");
188   TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
189   canvas->Divide(2, 1);
190   canvas->cd(1);
191
192   TH1** ptDists = new TH1*[4];
193
194   TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
195
196   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
197   const char* particleNames[] = { "#pi", "K", "p", "other" };
198   for (Int_t i=0; i<4; ++i)
199   {
200     AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderNames[i], folderNames[i]);
201     dNdEtaCorrection->LoadHistograms(fileName, folderNames[i]);
202
203     TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
204
205     gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
206     gene->GetXaxis()->SetRangeUser(-10, 10);
207
208     ptDists[i] = dynamic_cast<TH1*> (gene->Project3D("z")->Clone(Form("%s_z", folderNames[i])));
209     ptDists[i]->SetTitle(";p_{T};Count");
210     if (!ptDists[i])
211     {
212       printf("Problem getting distribution %d.\n", i);
213       return 0;
214     }
215
216     ptDists[i]->GetYaxis()->SetRangeUser(1, ptDists[i]->GetMaximum()*1.1);
217     ptDists[i]->GetXaxis()->SetRangeUser(0, 9.9);
218     ptDists[i]->SetLineColor(i+1);
219     ptDists[i]->DrawCopy((i == 0) ? "" : "SAME");
220     ptDists[i]->GetYaxis()->SetRange(0, 0);
221
222     legend->AddEntry(ptDists[i], particleNames[i]);
223   }
224   gPad->SetLogy();
225
226   TH1* total = dynamic_cast<TH1*> (ptDists[0]->Clone("total"));
227
228   for (Int_t i=1; i<4; ++i)
229     total->Add(ptDists[i]);
230
231   canvas->cd(2);
232   for (Int_t i=0; i<4; ++i)
233   {
234     ptDists[i]->Divide(total);
235     ptDists[i]->SetStats(kFALSE);
236     ptDists[i]->SetTitle(";p_{T};Fraction of total");
237     ptDists[i]->GetYaxis()->SetRangeUser(0, 1);
238     ptDists[i]->Draw((i == 0) ? "" : "SAME");
239   }
240   legend->SetFillColor(0);
241   legend->Draw();
242
243   canvas->SaveAs("DrawRatios.gif");
244
245
246   canvas = new TCanvas("PythiaRatios", "PythiaRatios", 500, 500);
247   for (Int_t i=0; i<4; ++i)
248   {
249     TH1* hist = ptDists[i]->Clone();
250     hist->GetXaxis()->SetRangeUser(0, 1.9);
251     hist->Draw((i == 0) ? "" : "SAME");
252   }
253   legend->Draw();
254
255   canvas->SaveAs("pythiaratios.eps");
256
257   file->Close();
258
259   return ptDists;
260 }
261
262 void DrawpiKpAndCombinedZOnly(Float_t upperPtLimit=0.99)
263 {
264   gROOT->ProcessLine(".L drawPlots.C");
265   gSystem->Load("libPWG0base");
266
267   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "correction_map.root" };
268   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "dndeta_correction" };
269   const char* legendNames[] = { "#pi", "K", "p", "standard" };
270   Int_t folderCount = 3;
271
272   /*const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "systematics.root", "correction_map.root" };
273   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3", "dndeta_correction" };
274   const char* legendNames[] = { "#pi", "K", "p", "others", "standard" };
275   Int_t folderCount = 5;*/
276
277   TString canvasName;
278   canvasName.Form("Track2Particle1DComposition");
279   TCanvas* canvas = new TCanvas(canvasName, canvasName, 700, 500);
280   canvas->SetGridx();
281   canvas->SetGridy();
282   canvas->SetBottomMargin(0.12);
283   //InitPad();
284
285   TLegend* legend = new TLegend(0.8, 0.7, 0.95, 0.95);
286   legend->SetFillColor(0);
287
288   Int_t mycolors[] = {1, 2, 4};
289
290   for (Int_t i=0; i<folderCount; ++i)
291   {
292     Track2Particle1DCreatePlots(fileNames[i], folderNames[i], upperPtLimit);
293
294     TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", folderNames[i], folderNames[i])));
295
296     Prepare1DPlot(corrZ);
297
298     corrZ->SetTitle("");
299     corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
300     corrZ->GetYaxis()->SetRangeUser(0.51, 6);
301     corrZ->SetMarkerColor(mycolors[i]);
302     corrZ->SetLineColor(mycolors[i]);
303     corrZ->SetMarkerStyle(markers[i+1]);
304     corrZ->GetYaxis()->SetTitle("correction factor");
305
306     corrZ->DrawCopy(((i>0) ? "SAMEP" : "P"));
307
308     legend->AddEntry(corrZ, legendNames[i]);
309   }
310
311   legend->Draw();
312
313   canvas->SaveAs("ptcutoff_species.eps");
314 }
315
316 void DrawCompareToReal()
317 {
318   gROOT->ProcessLine(".L drawPlots.C");
319
320   const char* fileNames[] = { "correction_map.root", "new_compositions.root" };
321   const char* folderNames[] = { "dndeta_correction", "PythiaRatios" };
322
323   Track2Particle1DComposition(fileNames, 2, folderNames);
324 }
325
326 void DrawDifferentSpecies()
327 {
328   gROOT->ProcessLine(".L drawPlots.C");
329
330   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "systematics.root" };
331   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
332
333   Track2Particle1DComposition(fileNames, 4, folderNames);
334 }
335
336 void DrawpiKpAndCombined()
337 {
338   gROOT->ProcessLine(".L drawPlots.C");
339
340   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "correction_map.root" };
341   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "dndeta_correction" };
342
343   Track2Particle1DComposition(fileNames, 4, folderNames);
344 }
345
346 void DrawSpeciesAndCombination()
347 {
348   gROOT->ProcessLine(".L drawPlots.C");
349
350   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "new_compositions.root" };
351   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "PythiaRatios" };
352
353   Track2Particle1DComposition(fileNames, 4, folderNames);
354 }
355
356 void DrawSimulatedVsCombined()
357 {
358   gROOT->ProcessLine(".L drawPlots.C");
359
360   const char* fileNames[] = { "new_compositions.root", "new_compositions.root" };
361   const char* folderNames[] = { "Pythia", "PythiaRatios" };
362
363   Track2Particle1DComposition(fileNames, 2, folderNames);
364 }
365
366 void DrawBoosts()
367 {
368   gROOT->ProcessLine(".L drawPlots.C");
369
370   const char* fileNames[] = { "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root" };
371   const char* folderNames[] = { "PythiaRatios", "PiBoosted", "KBoosted", "pBoosted" };
372
373   Track2Particle1DComposition(fileNames, 4, folderNames);
374 }
375
376 TH2F* HistogramDifferences(const char* filename, const char* folder1, const char* folder2)
377 {
378   gSystem->Load("libPWG0base");
379
380   AlidNdEtaCorrection* fdNdEtaCorrection[2];
381
382   TFile::Open(filename);
383
384   fdNdEtaCorrection[0] = new AlidNdEtaCorrection(folder1, folder1);
385   fdNdEtaCorrection[0]->LoadHistograms(filename, folder1);
386
387   fdNdEtaCorrection[1] = new AlidNdEtaCorrection(folder2, folder2);
388   fdNdEtaCorrection[1]->LoadHistograms(filename, folder2);
389
390   TH3F* hist1 = fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
391   TH3F* hist2 = fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
392
393   //TH1F* difference = new TH1F("difference", Form(";#DeltaC_{pT, z, #eta} %s / %s;Count", folder2, folder1), 1000, 0.9, 1.1);
394   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());
395
396   for (Int_t x=hist1->GetXaxis()->FindBin(-10); x<=hist1->GetXaxis()->FindBin(10); ++x)
397     for (Int_t y=hist1->GetYaxis()->FindBin(-0.8); y<=hist1->GetYaxis()->FindBin(0.8); ++y)
398       for (Int_t z=hist1->GetZaxis()->FindBin(0.3); z<=hist1->GetZaxis()->FindBin(9.9); ++z)
399         if (hist1->GetBinContent(x, y, z) != 0)
400           difference->Fill(hist2->GetBinContent(x, y, z) / hist1->GetBinContent(x, y, z), hist1->GetYaxis()->GetBinCenter(y));
401
402   difference->GetYaxis()->SetRangeUser(-0.8, 0.8);
403
404   printf("Over-/Underflow bins: %d %d\n", difference->GetBinContent(0), difference->GetBinContent(difference->GetNbinsX()+1));
405
406   return difference;
407 }
408
409 void HistogramDifferences()
410 {
411   TH2F* KBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "KBoosted");
412   TH2F* pBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "pBoosted");
413   TH2F* KReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "KReduced");
414   TH2F* pReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "pReduced");
415
416   TCanvas* canvas = new TCanvas("HistogramDifferences", "HistogramDifferences", 1000, 1000);
417   canvas->Divide(2, 2);
418
419   canvas->cd(1);
420   KBoosted->GetXaxis()->SetRangeUser(-0.05, 0.05);
421   KBoosted->Draw("COLZ");
422
423   canvas->cd(2);
424   KReduced->GetXaxis()->SetRangeUser(-0.05, 0.05);
425   KReduced->Draw("COLZ");
426
427   canvas->cd(3);
428   pBoosted->GetXaxis()->SetRangeUser(-0.02, 0.02);
429   pBoosted->Draw("COLZ");
430
431   canvas->cd(4);
432   pReduced->GetXaxis()->SetRangeUser(-0.02, 0.02);
433   pReduced->Draw("COLZ");
434
435   canvas->SaveAs("HistogramDifferences.gif");
436
437   canvas = new TCanvas("HistogramDifferencesProfile", "HistogramDifferencesProfile", 1000, 500);
438   canvas->Divide(2, 1);
439
440   canvas->cd(1);
441   gPad->SetBottomMargin(0.13);
442   KBoosted->SetTitle("a) Ratio of correction maps");
443   KBoosted->SetStats(kFALSE);
444   KBoosted->GetXaxis()->SetTitleOffset(1.4);
445   KBoosted->GetXaxis()->SetLabelOffset(0.02);
446   KBoosted->Draw("COLZ");
447
448   canvas->cd(2);
449   gPad->SetGridx();
450   gPad->SetGridy();
451
452   TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
453
454   for (Int_t i=0; i<4; ++i)
455   {
456     TH2F* hist = 0;
457     TString name;
458     switch (i)
459     {
460       case 0: hist = KBoosted; name = "K enhanced"; break;
461       case 1: hist = KReduced; name = "K reduced"; break;
462       case 2: hist = pBoosted; name = "p enhanced"; break;
463       case 3: hist = pReduced; name = "p reduced"; break;
464     }
465
466     TProfile* profile = hist->ProfileY();
467     profile->SetTitle("b) Mean and RMS");
468     profile->GetXaxis()->SetRange(hist->GetYaxis()->GetFirst(), hist->GetYaxis()->GetLast());
469     profile->GetXaxis()->SetTitleOffset(1.2);
470     profile->GetXaxis()->SetLabelOffset(0.02);
471     profile->GetYaxis()->SetRangeUser(0.98, 1.02);
472     profile->SetStats(kFALSE);
473     profile->SetLineColor(i+1);
474     profile->SetMarkerColor(i+1);
475     profile->DrawCopy(((i > 0) ? "SAME" : ""));
476
477
478     legend->AddEntry(profile, name);
479   }
480
481   legend->Draw();
482   canvas->SaveAs("particlecomposition_result.eps");
483 }
484
485
486 void ScalePtDependent(TH3F* hist, TF1* function)
487 {
488   // assumes that pt is the third dimension of hist
489   // scales with function(pt)
490
491   for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
492   {
493     Double_t factor = function->Eval(hist->GetZaxis()->GetBinCenter(z));
494     printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor);
495
496     for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
497       for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
498         hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor);
499   }
500 }
501
502 void ScalePtDependent(TH3F* hist, TH1* function)
503 {
504   // assumes that pt is the third dimension of hist
505   // scales with histogram(pt)
506
507   for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
508   {
509     Double_t factor = function->GetBinContent(function->GetXaxis()->FindBin(hist->GetZaxis()->GetBinCenter(z)));
510     printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor);
511
512     for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
513       for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
514         hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor);
515   }
516 }
517
518 const char* ChangeComposition(void** correctionPointer, Int_t index)
519 {
520   AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer;
521
522   switch (index)
523   {
524     case 0: // result from pp events
525       {
526         TFile::Open("pythiaratios.root");
527
528         for (Int_t i=0; i<4; ++i)
529         {
530           TString name;
531           name.Form("correction_%d", i);
532           fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
533           fdNdEtaCorrection[i]->LoadHistograms("pythiaratios.root", name);
534         }
535       }
536       return "Pythia";
537       break;
538
539     case 1: // each species rated with pythia ratios
540       /*TH1** ptDists = DrawRatios("pythiaratios.root");
541
542       for (Int_t i=0; i<3; ++i)
543       {
544         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]);
545         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]);
546       }*/
547       return "PythiaRatios";
548       break;
549
550       // one species enhanced / reduced
551     case 2: // + 50% pions
552     case 3: // - 50% pions
553     case 4: // + 50% kaons
554     case 5: // - 50% kaons
555     case 6: // + 50% protons
556     case 7: // - 50% protons
557       Int_t correctionIndex = (index - 2) / 2;
558       Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
559
560       fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Scale(scaleFactor);
561       fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Scale(scaleFactor);
562
563       TString* str = new TString;
564       str->Form("%s%s", (correctionIndex == 0) ? "Pi" : ((correctionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced");
565       return str->Data();
566       break;
567
568       // each species rated with pythia ratios
569     case 12: // + 50% pions
570     case 13: // - 50% pions
571     case 14: // + 50% kaons
572     case 15: // - 50% kaons
573     case 16: // + 50% protons
574     case 17: // - 50% protons
575       TH1** ptDists = DrawRatios("pythiaratios.root");
576       Int_t functionIndex = (index - 2) / 2;
577       Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
578       ptDists[functionIndex]->Scale(scaleFactor);
579
580       for (Int_t i=0; i<3; ++i)
581       {
582         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]);
583         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]);
584       }
585       TString* str = new TString;
586       str->Form("%s%s", (functionIndex == 0) ? "Pi" : ((functionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced");
587       return str->Data();
588       break;
589
590     case 999:
591       TF1* ptDependence = new TF1("simple", "x", 0, 100);
592       ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDependence);
593       ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDependence);
594       break;
595
596   }
597
598   return "noname";
599 }
600
601 void Composition()
602 {
603   gSystem->Load("libPWG0base");
604
605   gSystem->Unlink("new_compositions.root");
606
607   Int_t nCompositions = 8;
608   for (Int_t comp = 0; comp < nCompositions; ++comp)
609   {
610     AlidNdEtaCorrection* fdNdEtaCorrection[4];
611
612     TFile::Open("systematics.root");
613
614     for (Int_t i=0; i<4; ++i)
615     {
616       TString name;
617       name.Form("correction_%d", i);
618       fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
619       fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
620     }
621
622     const char* newName = ChangeComposition(fdNdEtaCorrection, comp);
623
624     Double_t geneCount[5];
625     Double_t measCount[5];
626     geneCount[4] = 0;
627     measCount[4] = 0;
628
629     for (Int_t i=0; i<4; ++i)
630     {
631       geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Integral();
632       measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Integral();
633
634       geneCount[4] += geneCount[i];
635       measCount[4] += measCount[i];
636
637       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]);
638     }
639
640     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]);
641
642     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]);
643
644     TList* collection = new TList;
645
646     for (Int_t i=0; i<4; ++i)
647       collection->Add(fdNdEtaCorrection[i]);
648
649     AlidNdEtaCorrection* newComposition = new AlidNdEtaCorrection(newName, newName);
650     newComposition->Merge(collection);
651     newComposition->Finish();
652
653     delete collection;
654
655     TFile* file = TFile::Open("new_compositions.root", "UPDATE");
656     newComposition->SaveHistograms();
657     //file->Write();
658     file->Close();
659   }
660
661   gROOT->ProcessLine(".L drawPlots.C");
662
663   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" };
664   const char* folderNames[] = { "Pythia", "PythiaRatios", "PiBoosted", "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced" };
665
666   Track2Particle1DComposition(fileNames, nCompositions, folderNames);
667 }
668
669 Double_t ConvSigma1To2D(Double_t sigma)
670 {
671   return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
672 }
673
674 Double_t ConvDistance1DTo2D(Double_t distance)
675 {
676   return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
677 }
678
679 Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
680 {
681   Double_t count = 0;
682
683   //nSigma = ConvSigma1To2D(nSigma);
684
685   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
686     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
687     {
688       Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
689       Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
690
691       Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
692
693       d = ConvDistance1DTo2D(d);
694
695       if (d < nSigma)
696         count += tracks->GetBinContent(x, y);
697     }
698
699   return count;
700 }
701
702 TH2F* Sigma2VertexGaussianTracksHist()
703 {
704   TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
705
706   TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
707   gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
708
709   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
710     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
711       tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
712
713   //normalize
714   tracks->Scale(1.0 / tracks->Integral());
715
716   return tracks;
717 }
718
719 TH1F* Sigma2VertexGaussian()
720 {
721   TH2F* tracks = Sigma2VertexGaussianTracksHist();
722
723   TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
724   canvas->Divide(2, 2);
725
726   canvas->cd(1);
727   tracks->Draw("COLZ");
728
729   TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 50, 0.05, 5.05);
730   for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
731     ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
732   ratio->SetMarkerStyle(21);
733
734   canvas->cd(2);
735   ratio->DrawCopy("P");
736
737   TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 3 sigma / % included n sigma", 50, 0.05, 5.05);
738   Double_t sigma3 = Sigma2VertexCount(tracks, 3);
739   for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
740     ratio2->Fill(nSigma, sigma3 / ratio->GetBinContent(ratio->FindBin(nSigma)));
741   ratio2->SetMarkerStyle(21);
742
743   canvas->cd(3);
744   ratio2->DrawCopy("P");
745
746   canvas->SaveAs("Sigma2Vertex.eps");
747
748   return ratio2;
749 }
750
751 TH1F* Sigma2VertexSimulation()
752 {
753   TFile* file = TFile::Open("systematics.root");
754
755   TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertex"));
756   if (!sigmavertex)
757   {
758     printf("Could not read histogram\n");
759     return;
760   }
761
762   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());
763
764   for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
765     ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(3)) / sigmavertex->GetBinContent(i));
766
767   TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
768   canvas->Divide(2, 1);
769
770   canvas->cd(1);
771   sigmavertex->SetMarkerStyle(21);
772   sigmavertex->Draw("P");
773
774   canvas->cd(2);
775   ratio->SetMarkerStyle(21);
776   ratio->DrawCopy("P");
777
778   return ratio;
779 }
780
781 void Sigma2VertexCompare()
782 {
783   TH1F* ratio1 = Sigma2VertexGaussian();
784   TH1F* ratio2 = Sigma2VertexSimulation();
785
786   ratio1->SetStats(kFALSE);
787   ratio2->SetStats(kFALSE);
788
789   ratio1->SetMarkerStyle(0);
790   ratio2->SetMarkerStyle(0);
791
792   ratio1->SetLineWidth(2);
793   ratio2->SetLineWidth(2);
794
795   TLegend* legend = new TLegend(0.7, 0.8, 0.95, 0.95);
796   legend->SetFillColor(0);
797   legend->AddEntry(ratio1, "Gaussian");
798   legend->AddEntry(ratio2, "Simulation");
799
800   ratio2->SetTitle("");
801   ratio2->GetYaxis()->SetTitleOffset(1.5);
802   ratio2->GetXaxis()->SetRangeUser(2, 4);
803
804   TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
805   InitPad();
806
807   ratio2->SetMarkerStyle(21);
808   ratio1->SetMarkerStyle(22);
809
810   ratio2->GetYaxis()->SetRangeUser(0.8, 1.2);
811   ratio2->SetLineColor(kRed);
812   ratio2->SetMarkerColor(kRed);
813   ratio2->Draw("PL");
814   ratio1->Draw("SAMEPL");
815
816   legend->Draw();
817
818   canvas->SaveAs("Sigma2VertexCompare.eps");
819 }
820
821 void drawSystematics()
822 {
823   //Secondaries();
824   //DrawDifferentSpecies();
825   //Composition();
826
827   Sigma2VertexSimulation();
828
829 }
830
831 void DrawdNdEtaDifferences()
832 {
833   TH1* hists[5];
834
835   TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
836   legend->SetFillColor(0);
837
838   TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
839   canvas->Divide(2, 1);
840
841   canvas->cd(1);
842
843   for (Int_t i=0; i<5; ++i)
844   {
845     hists[i] = 0;
846     TFile* file = 0;
847     TString title;
848
849     switch(i)
850     {
851       case 0 : file = TFile::Open("systematics_dndeta_reference.root"); title = "standard composition"; break;
852       case 1 : file = TFile::Open("systematics_dndeta_KBoosted.root"); title = "+ 50% kaons"; break;
853       case 2 : file = TFile::Open("systematics_dndeta_KReduced.root"); title = "- 50% kaons"; break;
854       case 3 : file = TFile::Open("systematics_dndeta_pBoosted.root"); title = "+ 50% protons"; break;
855       case 4 : file = TFile::Open("systematics_dndeta_pReduced.root"); title = "- 50% protons"; break;
856       default: return;
857     }
858
859     if (file)
860     {
861       hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
862       hists[i]->SetTitle("a)");
863
864       Prepare1DPlot(hists[i], kFALSE);
865       hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
866       hists[i]->GetYaxis()->SetRangeUser(6, 7);
867       hists[i]->SetLineColor(colors[i]);
868       hists[i]->SetMarkerColor(colors[i]);
869       hists[i]->SetMarkerStyle(markers[i]);
870       hists[i]->GetXaxis()->SetLabelOffset(0.015);
871       hists[i]->GetYaxis()->SetTitleOffset(1.5);
872       gPad->SetLeftMargin(0.12);
873       hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
874
875       legend->AddEntry(hists[i], title);
876       hists[i]->SetTitle(title);
877     }
878   }
879   legend->Draw();
880
881   canvas->cd(2);
882   gPad->SetLeftMargin(0.14);
883
884   TLegend* legend2 = new TLegend(0.73, 0.73, 0.98, 0.98);
885   legend2->SetFillColor(0);
886
887   for (Int_t i=1; i<5; ++i)
888   {
889     if (hists[i])
890     {
891       legend2->AddEntry(hists[i]);
892
893       hists[i]->Divide(hists[0]);
894       hists[i]->SetTitle("b)");
895       hists[i]->SetLineColor(colors[i-1]);
896       hists[i]->SetMarkerColor(colors[i-1]);
897       hists[i]->GetYaxis()->SetRangeUser(0.95, 1.05);
898       hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
899       hists[i]->GetYaxis()->SetTitleOffset(1.8);
900       hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
901     }
902   }
903
904   legend2->Draw();
905
906   canvas->SaveAs("particlecomposition_result_detail.gif");
907
908   TCanvas* canvas2 = new TCanvas("DrawdNdEtaDifferences2", "DrawdNdEtaDifferences2", 700, 500);
909
910   for (Int_t i=1; i<5; ++i)
911   {
912     if (hists[i])
913     {
914       hists[i]->SetTitle("");
915       hists[i]->GetYaxis()->SetTitleOffset(1.1);
916       hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
917     }
918   }
919
920   legend2->Draw();
921
922   canvas2->SaveAs("particlecomposition_result.gif");
923   canvas2->SaveAs("particlecomposition_result.eps");
924 }
925
926 mergeCorrectionsWithDifferentCrosssections(Char_t* standardCorrectionFileName="correction_map.root",
927                                              Char_t* systematicCorrectionFileName="systematics.root",
928                                              Char_t* outputFileName="systematics_vtxtrigger_compositions.root") {
929   //
930   // Function used to merge standard corrections with vertex
931   // reconstruction corrections obtained by a certain mix of ND, DD
932   // and SD events.
933   // 
934
935   gSystem->Load("libPWG0base");
936
937   const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
938   const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless"};
939   Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.0, 1.0, 1.5, 0.5};
940   Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5};
941
942   // cross section from Pythia
943   Float_t sigmaND = 55.2;
944   Float_t sigmaDD = 9.78;
945   Float_t sigmaSD = 14.30;
946
947   AlidNdEtaCorrection* corrections[21];
948   for (Int_t j=0; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)
949     AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
950     correctionStandard->LoadHistograms(standardCorrectionFileName);
951
952     // dont take vertexreco from this one
953     correctionStandard->GetVertexRecoCorrection()->Reset();
954     // dont take triggerbias from this one
955     correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
956
957     for (Int_t i=0; i<7; i++) {
958       TString name;
959       name.Form("dndeta_correction_syst_%s_%s", typeName[j], changes[i]);
960       AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
961
962       name.Form("vertexRecoND");
963       AlidNdEtaCorrection* dNdEtaCorrectionVtxND = new AlidNdEtaCorrection(name,name);
964       dNdEtaCorrectionVtxND->LoadHistograms(systematicCorrectionFileName, name);
965       name.Form("vertexRecoDD");
966       AlidNdEtaCorrection* dNdEtaCorrectionVtxDD = new AlidNdEtaCorrection(name,name);
967       dNdEtaCorrectionVtxDD->LoadHistograms(systematicCorrectionFileName, name);
968       name.Form("vertexRecoSD");
969       AlidNdEtaCorrection* dNdEtaCorrectionVtxSD = new AlidNdEtaCorrection(name,name);
970       dNdEtaCorrectionVtxSD->LoadHistograms(systematicCorrectionFileName, name);
971
972       name.Form("triggerBiasND");
973       AlidNdEtaCorrection* dNdEtaCorrectionTriggerND = new AlidNdEtaCorrection(name,name);
974       dNdEtaCorrectionTriggerND->LoadHistograms(systematicCorrectionFileName, name);
975       name.Form("triggerBiasDD");
976       AlidNdEtaCorrection* dNdEtaCorrectionTriggerDD = new AlidNdEtaCorrection(name,name);
977       dNdEtaCorrectionTriggerDD->LoadHistograms(systematicCorrectionFileName, name);
978       name.Form("triggerBiasSD");
979       AlidNdEtaCorrection* dNdEtaCorrectionTriggerSD = new AlidNdEtaCorrection(name,name);
980       dNdEtaCorrectionTriggerSD->LoadHistograms(systematicCorrectionFileName, name);
981
982       // calculating relative
983       Float_t nd = 100 * sigmaND/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
984       Float_t dd = 100 * (scalesDD[i]*sigmaDD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
985       Float_t sd = 100 * (scalesSD[i]*sigmaSD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
986
987       printf(Form("%s : ND=%.1f\%, DD=%.1f\%, SD=%.1f\% \n",changes[i],nd,dd,sd));
988       current->SetTitle(Form("ND=%.2f\%,DD=%.2f\%,SD=%.2f\%",nd,dd,sd));
989       current->SetTitle(name);
990
991       // scale
992       if (j == 0 || j == 2)
993       {
994         dNdEtaCorrectionVtxDD->GetVertexRecoCorrection()->GetMeasuredHistogram()->Scale(scalesDD[i]);
995         dNdEtaCorrectionVtxSD->GetVertexRecoCorrection()->GetMeasuredHistogram()->Scale(scalesSD[i]);
996         dNdEtaCorrectionVtxDD->GetVertexRecoCorrection()->GetGeneratedHistogram()->Scale(scalesDD[i]);
997         dNdEtaCorrectionVtxSD->GetVertexRecoCorrection()->GetGeneratedHistogram()->Scale(scalesSD[i]);
998       }
999       if (j == 1 || j == 2)
1000       {
1001         dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionINEL()->GetMeasuredHistogram()->Scale(scalesDD[i]);
1002         dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionINEL()->GetMeasuredHistogram()->Scale(scalesSD[i]);
1003         dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionINEL()->GetGeneratedHistogram()->Scale(scalesDD[i]);
1004         dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionINEL()->GetGeneratedHistogram()->Scale(scalesSD[i]);
1005       }
1006
1007       //clear track, trigger in Vtx correction
1008       dNdEtaCorrectionVtxND->GetTrack2ParticleCorrection()->Reset();
1009       dNdEtaCorrectionVtxND->GetTriggerBiasCorrectionNSD()->Reset();
1010       dNdEtaCorrectionVtxND->GetTriggerBiasCorrectionND()->Reset();
1011       dNdEtaCorrectionVtxND->GetTriggerBiasCorrectionINEL()->Reset();
1012       dNdEtaCorrectionVtxSD->GetTrack2ParticleCorrection()->Reset();
1013       dNdEtaCorrectionVtxSD->GetTriggerBiasCorrectionNSD()->Reset();
1014       dNdEtaCorrectionVtxSD->GetTriggerBiasCorrectionND()->Reset();
1015       dNdEtaCorrectionVtxSD->GetTriggerBiasCorrectionINEL()->Reset();
1016       dNdEtaCorrectionVtxDD->GetTrack2ParticleCorrection()->Reset();
1017       dNdEtaCorrectionVtxDD->GetTriggerBiasCorrectionNSD()->Reset();
1018       dNdEtaCorrectionVtxDD->GetTriggerBiasCorrectionND()->Reset();
1019       dNdEtaCorrectionVtxDD->GetTriggerBiasCorrectionINEL()->Reset();
1020
1021       //clear track, vertexreco in trigger correction
1022       dNdEtaCorrectionTriggerND->GetTrack2ParticleCorrection()->Reset();
1023       dNdEtaCorrectionTriggerND->GetTriggerBiasCorrectionNSD()->Reset();
1024       dNdEtaCorrectionTriggerND->GetTriggerBiasCorrectionND()->Reset();
1025       dNdEtaCorrectionTriggerND->GetVertexRecoCorrection()->Reset();
1026       dNdEtaCorrectionTriggerSD->GetTrack2ParticleCorrection()->Reset();
1027       dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionNSD()->Reset();
1028       dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionND()->Reset();
1029       dNdEtaCorrectionTriggerSD->GetVertexRecoCorrection()->Reset();
1030       dNdEtaCorrectionTriggerDD->GetTrack2ParticleCorrection()->Reset();
1031       dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionNSD()->Reset();
1032       dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionND()->Reset();
1033       dNdEtaCorrectionTriggerDD->GetVertexRecoCorrection()->Reset();
1034
1035       TList collection;
1036       collection.Add(correctionStandard);
1037       collection.Add(dNdEtaCorrectionVtxND);
1038       collection.Add(dNdEtaCorrectionVtxDD);
1039       collection.Add(dNdEtaCorrectionVtxSD);
1040       collection.Add(dNdEtaCorrectionTriggerND);
1041       collection.Add(dNdEtaCorrectionTriggerDD);
1042       collection.Add(dNdEtaCorrectionTriggerSD);
1043
1044       current->Merge(&collection);
1045       current->Finish();
1046
1047       corrections[i+j*7] = current;
1048     }
1049   }
1050
1051   TFile* fout = new TFile(outputFileName,"RECREATE");
1052
1053   for (Int_t i=0; i<21; i++)
1054     corrections[i]->SaveHistograms();
1055
1056   fout->Write();
1057   fout->Close();
1058 }
1059
1060
1061 DrawVertexRecoSyst() {
1062   // Draws the ratio of the dN/dEta obtained with changed SD and DD
1063   // cross-sections vertex reco corrections to the dN/dEta obtained
1064   // from the standard pythia cross-sections 
1065   //
1066   // The files with the vertex reco corrections for different
1067   // processes (and with the other standard corrections) are expected
1068   // to have the names "analysis_esd_X.root", where the Xs are defined
1069   // in the array changes below.
1070
1071   Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless"};
1072   Char_t* descr[]  =   {"",
1073                         "#sigma_{DD}' = 1.5#times#sigma_{DD}",
1074                         "#sigma_{DD}' = 0.5#times#sigma_{DD}",
1075                         "#sigma_{SD}' = 1.5#times#sigma_{SD}",
1076                         "#sigma_{SD}' = 0.5#times#sigma_{SD}",
1077                         "#sigma_{D}' = 1.5#times#sigma_{D}",
1078                         "#sigma_{D}' = 0.5#times#sigma_{D}"};
1079
1080   Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.5, 0.5, 1.5, 0.5};
1081   Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.5, 0.5};
1082
1083   Int_t markers[] = {20,20,21,22,23,28,29};
1084   Int_t colors[]  = {1,2,3,4,6,8,102};
1085
1086   // cross section from Pythia
1087   Float_t sigmaND = 55.2;
1088   Float_t sigmaDD = 9.78;
1089   Float_t sigmaSD = 14.30;
1090
1091   TH1F* dNdEta[7];
1092   TH1F* ratios[7];
1093
1094   TFile* fin;
1095
1096   for (Int_t i=0; i<7; i++) {
1097     // calculating relative
1098     fin = TFile::Open(Form("analysis_esd_%s.root",changes[i]));
1099
1100     dNdEta[i] = (TH1F*)(fin->Get("dndeta/dndeta_dNdEta_corrected_2"))->Clone();
1101
1102     for (Int_t b=0; b<dNdEta[i]->GetNbinsX(); b++) {
1103       if (TMath::Abs(dNdEta[i]->GetBinCenter(b))>0.9) {
1104         dNdEta[i]->SetBinContent(b,0);
1105         dNdEta[i]->SetBinError(b,0);
1106       }
1107     }
1108
1109     dNdEta[i]->Rebin(4);
1110
1111     dNdEta[i]->SetLineWidth(2);
1112     dNdEta[i]->SetLineColor(colors[i]);
1113     dNdEta[i]->SetMarkerStyle(markers[i]);
1114     dNdEta[i]->SetMarkerSize(0.9);
1115     dNdEta[i]->SetMarkerColor(colors[i]);
1116
1117     ratios[i] = (TH1F*)dNdEta[i]->Clone("ratio");
1118     ratios[i]->Divide(ratios[i],dNdEta[0],1,1,"B");
1119     
1120     ratios[i]->SetName(changes[i]);
1121     ratios[i]->SetTitle(changes[i]);
1122   }
1123   
1124   //##########################################################
1125   
1126   gStyle->SetOptStat(0);
1127   gStyle->SetOptTitle(0);
1128   gStyle->SetOptFit(0);
1129
1130   gStyle->SetTextSize(0.2);
1131   gStyle->SetTitleSize(0.05,"xyz");
1132   gStyle->SetLabelOffset(0.01, "xyz");
1133
1134
1135   gStyle->SetTitleOffset(1.2, "y");
1136   gStyle->SetTitleOffset(1.2, "x");
1137   gStyle->SetEndErrorSize(0.0);
1138
1139   //##############################################
1140
1141   //making canvas and pads
1142   TCanvas *c = new TCanvas(Form("vertex_reco_syst_%s", plot), Form("vertex_reco_syst_%s", plot),600,500);
1143
1144   TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
1145
1146   p1->SetBottomMargin(0.15);
1147   p1->SetTopMargin(0.03);
1148   p1->SetLeftMargin(0.15);
1149   p1->SetRightMargin(0.03);
1150
1151   p1->SetGridx();
1152   p1->SetGridy();
1153
1154   p1->Draw();
1155   p1->cd();
1156   
1157   
1158   TH2F* null = new TH2F("","",100,-1.5,1.5,100,0.9601,1.099);
1159   null->SetXTitle("#eta");
1160   null->GetXaxis()->CenterTitle(kTRUE);
1161   null->SetYTitle("dN/d#eta / dN/d#eta_{pythia}");
1162   null->GetYaxis()->CenterTitle(kTRUE);
1163   null->Draw();
1164
1165   for (Int_t i=1; i<7; i++) 
1166     ratios[i]->Draw("same");
1167
1168   TLegend* legend = new TLegend(0.6, 0.6, 0.95, 0.95);
1169   legend->SetFillColor(0);
1170
1171   TLatex* text[7];
1172   for (Int_t i=1; i<7; i++) {
1173     legend->AddEntry(dNdEta[i], descr[i]);
1174   }
1175
1176   legend->Draw();
1177   //text(0.2,0.88,"Effect of changing",0.045,1,kTRUE);
1178   //text(0.2,0.83,"relative cross-sections",0.045,1,kTRUE);
1179   //text(0.2,0.78,"(vertex reconstruction corr.)",0.043,13,kTRUE);
1180
1181   c->SaveAs(Form("%s.gif", c->GetName()));
1182   c->SaveAs(Form("%s.eps", c->GetName()));
1183 }
1184
1185
1186
1187 DrawTriggerEfficiency(Char_t* fileName) {
1188
1189   gStyle->SetOptStat(0);
1190   gStyle->SetOptTitle(0);
1191   gStyle->SetOptFit(0);
1192
1193   gStyle->SetTextSize(0.04);
1194   gStyle->SetTitleSize(0.05,"xyz");
1195   //gStyle->SetTitleFont(133, "xyz");
1196   //gStyle->SetLabelFont(133, "xyz");
1197   //gStyle->SetLabelSize(17, "xyz");
1198   gStyle->SetLabelOffset(0.01, "xyz");
1199
1200   gStyle->SetTitleOffset(1.1, "y");
1201   gStyle->SetTitleOffset(1.1, "x");
1202   gStyle->SetEndErrorSize(0.0);
1203
1204   //##############################################
1205
1206   //making canvas and pads
1207   TCanvas *c = new TCanvas("trigger_eff", "",600,500);
1208
1209   TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
1210
1211   p1->SetBottomMargin(0.15);
1212   p1->SetTopMargin(0.03);
1213   p1->SetLeftMargin(0.15);
1214   p1->SetRightMargin(0.03);
1215   
1216   p1->SetGridx();
1217   p1->SetGridy();
1218
1219   p1->Draw();
1220   p1->cd();
1221
1222   gSystem->Load("libPWG0base");
1223
1224   AlidNdEtaCorrection* corrections[4];
1225   AliCorrectionMatrix2D* triggerBiasCorrection[4];
1226
1227   TH1F* hTriggerEffInv[4];
1228   TH1F* hTriggerEff[4];
1229
1230   Char_t* names[] = {"triggerBiasND", "triggerBiasDD", "triggerBiasSD", "triggerBiasINEL"};
1231   
1232   for (Int_t i=0; i<4; i++) {
1233     corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1234     corrections[i]->LoadHistograms(fileName, names[i]);    
1235
1236     triggerBiasCorrection[i] = corrections[i]->GetTriggerBiasCorrectionINEL();
1237
1238     
1239     hTriggerEffInv[i] = triggerBiasCorrection[i]->Get1DCorrection();
1240     hTriggerEff[i]    = (TH1F*)hTriggerEffInv[i]->Clone();
1241     
1242     for (Int_t b=0; b<=hTriggerEff[i]->GetNbinsX(); b++) {
1243       hTriggerEff[i]->SetBinContent(b,1);
1244       hTriggerEff[i]->SetBinError(b,0);
1245     }
1246     hTriggerEff[i]->Divide(hTriggerEff[i],hTriggerEffInv[i]);
1247     hTriggerEff[i]->Scale(100);
1248   }
1249
1250   Int_t colors[] = {2,3,4,1};
1251   Float_t effs[4];
1252   for (Int_t i=0; i<4; i++) {
1253     hTriggerEff[i]->Fit("pol0","","",-20,20);
1254     effs[i] = ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->GetParameter(0);
1255     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineWidth(1);
1256     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineStyle(2);
1257     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineColor(colors[i]);
1258     cout << effs[i] << endl;
1259   }
1260
1261
1262   Char_t* text[] = {"ND", "DD", "SD", "INEL"};
1263   TLatex* latex[4];
1264
1265   TH2F* null = new TH2F("","",100,-25,35,100,0,110);
1266   null->SetXTitle("Vertex z [cm]");
1267   null->GetXaxis()->CenterTitle(kTRUE);
1268   null->SetYTitle("Trigger efficiency [%]");
1269   null->GetYaxis()->CenterTitle(kTRUE);
1270   null->Draw();
1271
1272
1273   for (Int_t i=0; i<4; i++) {
1274     hTriggerEff[i]->SetLineWidth(2);
1275     hTriggerEff[i]->SetLineColor(colors[i]);
1276
1277     hTriggerEff[i]->Draw("same");
1278
1279     latex[i] = new TLatex(22,effs[i]-1.5, Form("%s (%0.1f)",text[i],effs[i]));
1280     latex[i]->SetTextColor(colors[i]);
1281     latex[i]->Draw();
1282   }
1283   
1284 }
1285
1286
1287 DrawSpectraPID(Char_t* fileName) {
1288
1289   gSystem->Load("libPWG0base");
1290
1291   Char_t* names[] = {"correction_0", "correction_1", "correction_2", "correction_3"};
1292   AlidNdEtaCorrection* corrections[4];
1293   AliCorrectionMatrix3D* trackToPartCorrection[4];
1294
1295   TH1F* measuredPt[4];
1296   TH1F* generatedPt[4];
1297   TH1F* ratioPt[4];
1298
1299   for (Int_t i=0; i<4; i++) {
1300     corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1301     corrections[i]->LoadHistograms(fileName, names[i]);    
1302
1303     trackToPartCorrection[i] = corrections[i]->GetTrack2ParticleCorrection();
1304
1305     Int_t binX1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(-10);
1306     Int_t binX2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(10);
1307     Int_t binY1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(-1);
1308     Int_t binY2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(1);
1309
1310     measuredPt[i]  = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->ProjectionZ(Form("m_%d",i),binX1,binX2,binY1,binY2);
1311     generatedPt[i] = (TH1F*)trackToPartCorrection[i]->GetGeneratedHistogram()->ProjectionZ(Form("g_%d",i),binX1,binX2,binY1,binY2);
1312     ratioPt[i] = (TH1F*)generatedPt[i]->Clone(Form("r_%d",i));
1313     ratioPt[i]->Divide(measuredPt[i], generatedPt[i], 1,1,"B");
1314   }
1315   
1316   ratioPt[0]->Draw();
1317   
1318   for (Int_t i=0; i<3; i++) {
1319     ratioPt[i]->SetLineColor(i+1);
1320     ratioPt[i]->SetLineWidth(2);
1321     
1322     ratioPt[i]->Draw("same");
1323     
1324   }
1325
1326   return;
1327   measuredPt[0]->SetLineColor(2);
1328   measuredPt[0]->SetLineWidth(5);
1329
1330   measuredPt[0]->Draw();
1331   generatedPt[0]->Draw("same");
1332 }
1333
1334 void changePtSpectrum(const char* fileName = "analysis_mc.root")
1335 {
1336   TFile* file = TFile::Open(fileName);
1337   TH1F* hist = dynamic_cast<TH1F*> (file->Get("dndeta_check_pt"));
1338
1339   //hist->Rebin(3);
1340   //hist->Scale(1.0/3);
1341
1342   TH1F* clone1 = dynamic_cast<TH1F*> (hist->Clone("clone1"));
1343   TH1F* clone2 = dynamic_cast<TH1F*> (hist->Clone("clone2"));
1344
1345   TH1F* scale1 =  dynamic_cast<TH1F*> (hist->Clone("scale1"));
1346   TH1F* scale2 =  dynamic_cast<TH1F*> (hist->Clone("scale2"));
1347
1348   Float_t ptCutOff = 0.3;
1349
1350   for (Int_t i=1; i <= hist->GetNbinsX(); ++i)
1351   {
1352     if (hist->GetBinCenter(i) > ptCutOff)
1353     {
1354       scale1->SetBinContent(i, 1);
1355       scale2->SetBinContent(i, 1);
1356     }
1357     else
1358     {
1359       // 90 % at pt = 0, 0% at pt = ptcutoff
1360       scale1->SetBinContent(i, 1 - (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * 0.3);
1361
1362       // 110% at pt = 0, ...
1363       scale2->SetBinContent(i, 1 + (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * 0.3);
1364     }
1365     scale1->SetBinError(i, 0);
1366     scale2->SetBinError(i, 0);
1367   }
1368
1369   new TCanvas;
1370
1371   scale1->Draw();
1372   scale2->SetLineColor(kRed);
1373   scale2->Draw("SAME");
1374
1375   clone1->Multiply(scale1);
1376   clone2->Multiply(scale2);
1377
1378   Prepare1DPlot(hist);
1379   Prepare1DPlot(clone1);
1380   Prepare1DPlot(clone2);
1381
1382   /*hist->SetMarkerStyle(markers[0]);
1383   clone1->SetMarkerStyle(markers[0]);
1384   clone2->SetMarkerStyle(markers[0]);*/
1385
1386   hist->SetTitle(";p_{T} in GeV/c;dN_{ch}/dp_{T} in c/GeV");
1387   hist->GetXaxis()->SetRangeUser(0, 0.7);
1388   hist->GetYaxis()->SetRangeUser(0.01, clone2->GetMaximum() * 1.1);
1389   hist->GetYaxis()->SetTitleOffset(1);
1390
1391   TCanvas* canvas = new TCanvas;
1392   gPad->SetBottomMargin(0.12);
1393   hist->Draw("H");
1394   clone1->SetLineColor(kRed);
1395   clone1->Draw("HSAME");
1396   clone2->SetLineColor(kBlue);
1397   clone2->Draw("HSAME");
1398   hist->Draw("HSAME");
1399
1400   Float_t fraction =  hist->Integral(hist->GetXaxis()->FindBin(ptCutOff), hist->GetNbinsX()) / hist->Integral(1, hist->GetNbinsX());
1401   Float_t fraction1 = clone1->Integral(clone1->GetXaxis()->FindBin(ptCutOff), clone1->GetNbinsX()) / clone1->Integral(1, clone1->GetNbinsX());
1402   Float_t fraction2 = clone2->Integral(clone2->GetXaxis()->FindBin(ptCutOff), clone2->GetNbinsX()) / clone2->Integral(1, clone2->GetNbinsX());
1403
1404   printf("%f %f %f\n", fraction, fraction1, fraction2);
1405   printf("Rel. %f %f\n", fraction1 / fraction, fraction2 / fraction);
1406
1407   canvas->SaveAs("changePtSpectrum.gif");
1408   canvas->SaveAs("changePtSpectrum.eps");
1409 }
1410
1411 void FractionBelowPt()
1412 {
1413   gSystem->Load("libPWG0base");
1414
1415   AlidNdEtaCorrection* fdNdEtaCorrection[4];
1416
1417   TFile::Open("systematics.root");
1418
1419   for (Int_t i=0; i<4; ++i)
1420   {
1421     TString name;
1422     name.Form("correction_%d", i);
1423     fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
1424     fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
1425   }
1426
1427   Double_t geneCount[5];
1428   Double_t measCount[5];
1429   geneCount[4] = 0;
1430   measCount[4] = 0;
1431
1432   for (Int_t i=0; i<4; ++i)
1433   {
1434     TH3F* hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1435     geneCount[i] = hist->Integral(hist->GetXaxis()->FindBin(-10), hist->GetXaxis()->FindBin(10),
1436                                   hist->GetYaxis()->FindBin(-0.8), hist->GetYaxis()->FindBin(0.8),
1437                                   1, hist->GetZaxis()->FindBin(0.3));
1438
1439     hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1440     measCount[i] = hist->Integral(hist->GetXaxis()->FindBin(-10), hist->GetXaxis()->FindBin(10), hist->GetYaxis()->FindBin(-0.8), hist->GetYaxis()->FindBin(0.8), 1, hist->GetZaxis()->FindBin(0.3));
1441
1442     geneCount[4] += geneCount[i];
1443     measCount[4] += measCount[i];
1444
1445     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]);
1446   }
1447
1448   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]);
1449
1450   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]);
1451 }
1452
1453
1454 mergeCorrectionsMisalignment(Char_t* alignedFile = "correction_map_aligned.root",
1455                                              Char_t* misalignedFile = "correction_map_misaligned.root",
1456                                              Char_t* outputFileName="correction_map_misaligned_single.root")
1457 {
1458   //
1459   // from the aligned and misaligned corrections, 3 new corrections are created
1460   // in these new corrections only one of the corrections (track2particle, vertex, trigger)
1461   // is taken from the misaligned input to allow study of the effect on the different
1462   // corrections
1463
1464   gSystem->Load("libPWG0base");
1465
1466   const Char_t* typeName[] = { "track2particle", "vertex", "trigger" };
1467
1468   AlidNdEtaCorrection* corrections[3];
1469   for (Int_t j=0; j<3; j++) { // j = 0 (track2particle), j = 1 (vertex), j = 2 (trigger)
1470     AlidNdEtaCorrection* alignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1471     alignedCorrection->LoadHistograms(alignedFile);
1472
1473     AlidNdEtaCorrection* misalignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1474     misalignedCorrection->LoadHistograms(misalignedFile);
1475
1476     TString name;
1477     name.Form("dndeta_correction_alignment_%s", typeName[j]);
1478     AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
1479
1480     switch (j)
1481     {
1482       case 0:
1483         alignedCorrection->GetTrack2ParticleCorrection()->Reset();
1484         misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1485         misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1486         misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1487         misalignedCorrection->GetVertexRecoCorrection()->Reset();
1488         break;
1489
1490       case 1:
1491         misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1492         misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1493         misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1494         misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1495         alignedCorrection->GetVertexRecoCorrection()->Reset();
1496         break;
1497
1498       case 2:
1499         misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1500         alignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1501         alignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1502         alignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1503         misalignedCorrection->GetVertexRecoCorrection()->Reset();
1504         break;
1505
1506       default:
1507         return;
1508     }
1509
1510     TList collection;
1511     collection.Add(misalignedCorrection);
1512     collection.Add(alignedCorrection);
1513
1514     current->Merge(&collection);
1515     current->Finish();
1516
1517     corrections[j] = current;
1518   }
1519
1520   TFile* fout = new TFile(outputFileName, "RECREATE");
1521
1522   for (Int_t i=0; i<3; i++)
1523     corrections[i]->SaveHistograms();
1524
1525   fout->Write();
1526   fout->Close();
1527 }
1528
1529
1530 void DrawdNdEtaDifferencesAlignment()
1531 {
1532   TH1* hists[5];
1533
1534   TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
1535   legend->SetFillColor(0);
1536
1537   TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
1538   canvas->Divide(2, 1);
1539
1540   canvas->cd(1);
1541
1542   for (Int_t i=0; i<5; ++i)
1543   {
1544     hists[i] = 0;
1545     TFile* file = 0;
1546     TString title;
1547
1548     switch(i)
1549     {
1550       case 0 : file = TFile::Open("systematics_misalignment_aligned.root"); title = "aligned"; break;
1551       case 1 : file = TFile::Open("systematics_misalignment_misaligned.root"); title = "fully misaligned"; break;
1552       case 2 : file = TFile::Open("systematics_misalignment_track2particle.root"); title = "only track2particle"; break;
1553       case 3 : file = TFile::Open("systematics_misalignment_vertex.root"); title = "only vertex rec."; break;
1554       case 4 : file = TFile::Open("systematics_misalignment_trigger.root"); title = "only trigger bias"; break;
1555       default: return;
1556     }
1557
1558     if (file)
1559     {
1560       hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
1561       hists[i]->SetTitle("");
1562
1563       Prepare1DPlot(hists[i], kFALSE);
1564       hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
1565       hists[i]->GetYaxis()->SetRangeUser(6, 7);
1566       hists[i]->SetLineWidth(1);
1567       hists[i]->SetLineColor(colors[i]);
1568       hists[i]->SetMarkerColor(colors[i]);
1569       hists[i]->SetMarkerStyle(markers[i]);
1570       hists[i]->GetXaxis()->SetLabelOffset(0.015);
1571       hists[i]->GetYaxis()->SetTitleOffset(1.5);
1572       gPad->SetLeftMargin(0.12);
1573       hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
1574
1575       legend->AddEntry(hists[i], title);
1576       hists[i]->SetTitle(title);
1577     }
1578   }
1579   legend->Draw();
1580
1581   canvas->cd(2);
1582   gPad->SetLeftMargin(0.14);
1583
1584   TLegend* legend2 = new TLegend(0.63, 0.73, 0.98, 0.98);
1585   legend2->SetFillColor(0);
1586
1587   for (Int_t i=1; i<5; ++i)
1588   {
1589     if (hists[i])
1590     {
1591       legend2->AddEntry(hists[i]);
1592
1593       hists[i]->Divide(hists[0]);
1594       hists[i]->SetTitle("b)");
1595       hists[i]->GetYaxis()->SetRangeUser(0.9, 1.1);
1596       hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
1597       hists[i]->GetYaxis()->SetTitleOffset(1.8);
1598       hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
1599     }
1600   }
1601
1602   legend2->Draw();
1603
1604   canvas->SaveAs("misalignment_result.eps");
1605   canvas->SaveAs("misalignment_result.gif");
1606 }
1607
1608 void EvaluateMultiplicityEffect()
1609 {
1610   gSystem->Load("libPWG0base");
1611
1612   AlidNdEtaCorrection* fdNdEtaCorrectionLow[4];
1613   AlidNdEtaCorrection* fdNdEtaCorrectionHigh[4];
1614
1615   TFile::Open("systematics-low-multiplicity.root");
1616
1617   for (Int_t i=0; i<4; ++i)
1618   {
1619     TString name;
1620     name.Form("correction_%d", i);
1621     fdNdEtaCorrectionLow[i] = new AlidNdEtaCorrection(name, name);
1622     fdNdEtaCorrectionLow[i]->LoadHistograms("systematics-low-multiplicity.root", name);
1623   }
1624
1625   TList list;
1626   for (Int_t i=1; i<4; ++i)
1627     list.Add(fdNdEtaCorrectionLow[i]);
1628
1629   fdNdEtaCorrectionLow[0]->Merge(&list);
1630   fdNdEtaCorrectionLow[0]->Finish();
1631
1632   TFile::Open("systematics-high-multiplicity.root");
1633
1634   for (Int_t i=0; i<4; ++i)
1635   {
1636     TString name;
1637     name.Form("correction_%d", i);
1638     fdNdEtaCorrectionHigh[i] = new AlidNdEtaCorrection(name, name);
1639     fdNdEtaCorrectionHigh[i]->LoadHistograms("systematics-high-multiplicity.root", name);
1640   }
1641
1642   TList list2;
1643   for (Int_t i=1; i<4; ++i)
1644     list2.Add(fdNdEtaCorrectionHigh[i]);
1645
1646   fdNdEtaCorrectionHigh[0]->Merge(&list2);
1647   fdNdEtaCorrectionHigh[0]->Finish();
1648
1649   TH1F* outputLow = new TH1F("Track2ParticleLow", "Track2Particle at low multiplicity", 200, 0, 2);
1650   TH1F* outputHigh = new TH1F("Track2ParticleHigh", "Track2Particle at high multiplicity", 200, 0, 2);
1651
1652   TH3F* hist = fdNdEtaCorrectionLow[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1653   TH3F* hist2 = fdNdEtaCorrectionHigh[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1654   for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
1655     for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
1656       for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
1657       //for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
1658       {
1659         if (hist->GetBinContent(x, y, z) > 0)
1660           outputLow->Fill(hist->GetBinContent(x, y, z));
1661         //if (hist->GetBinContent(x, y, z) == 1)
1662         //  printf("z = %f, eta = %f, pt = %f: %f %f %f\n", hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z), hist->GetBinContent(x, y, z), fdNdEtaCorrectionLow[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->GetBinContent(x, y, z), fdNdEtaCorrectionLow[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->GetBinContent(x, y, z));
1663
1664         if (hist2->GetBinContent(x, y, z) > 0)
1665           outputHigh->Fill(hist2->GetBinContent(x, y, z));
1666       }
1667
1668   TCanvas* canvas = new TCanvas("EvaluateMultiplicityEffect", "EvaluateMultiplicityEffect", 1000, 500);
1669   canvas->Divide(2, 1);
1670
1671   canvas->cd(1);
1672   outputLow->Draw();
1673   outputLow->Fit("gaus", "0");
1674   outputLow->GetFunction("gaus")->SetLineColor(2);
1675   outputLow->GetFunction("gaus")->DrawCopy("SAME");
1676
1677   canvas->cd(2);
1678   outputHigh->Draw();
1679   outputHigh->Fit("gaus", "0");
1680   outputHigh->GetFunction("gaus")->DrawCopy("SAME");
1681
1682   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1683 }
1684
1685 void PlotErrors(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
1686 {
1687     gSystem->Load("libPWG0base");
1688         
1689         AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1690     dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
1691     
1692     dNdEtaCorrection->GetTrack2ParticleCorrection()->PlotBinErrors(xmin, xmax, ymin, ymax, zmin, zmax)->Draw();
1693 }
1694
1695
1696 void runStudy(const char* baseCorrectionMapFile = "correction_map.root", const char* baseCorrectionMapFolder = "dndeta_correction", const char* changedCorrectionMapFile = "correction_map.root", const char* changedCorrectionMapFolder = "dndeta_correction", const char* dataFile = "analysis_esd_raw.root", const char* output = "analysis_esd_syst.root")
1697 {
1698   gSystem->Load("libPWG0base");
1699
1700   TFile* file = TFile::Open(output, "RECREATE");
1701
1702   const Int_t max = 5;
1703   dNdEtaAnalysis* fdNdEtaAnalysis[5];
1704
1705   new TCanvas;
1706   TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
1707   legend->SetFillColor(0);
1708
1709   for (Int_t i = 0; i < max; ++i)
1710   {
1711     TFile::Open(baseCorrectionMapFile);
1712     AlidNdEtaCorrection* baseCorrection = new AlidNdEtaCorrection(baseCorrectionMapFolder, baseCorrectionMapFolder);
1713     baseCorrection->LoadHistograms();
1714
1715     AlidNdEtaCorrection::CorrectionType correctionType = AlidNdEtaCorrection::kNone;
1716     const char* name = 0;
1717
1718     TFile::Open(changedCorrectionMapFile);
1719     switch (i)
1720     {
1721       case 0 :
1722         name = "default";
1723         break;
1724
1725       case 1 :
1726         baseCorrection->GetTrack2ParticleCorrection()->LoadHistograms(Form("%s/Track2Particle", changedCorrectionMapFolder));
1727         name = "Track2Particle";
1728         break;
1729
1730       case 2 :
1731         baseCorrection->GetVertexRecoCorrection()->LoadHistograms(Form("%s/VertexReconstruction", changedCorrectionMapFolder));
1732         name = "VertexReco";
1733         break;
1734
1735       case 3 :
1736         baseCorrection->GetTriggerBiasCorrectionINEL()->LoadHistograms(Form("%s/TriggerBias_MBToINEL", changedCorrectionMapFolder));
1737         name = "TriggerBias_MBToINEL";
1738         break;
1739
1740       case 4 :
1741         baseCorrection->LoadHistograms(changedCorrectionMapFolder);
1742         name = "all";
1743         break;
1744
1745       default: return;
1746     }
1747
1748     TFile::Open(dataFile);
1749     fdNdEtaAnalysis[i] = new dNdEtaAnalysis(name, name);
1750     fdNdEtaAnalysis[i]->LoadHistograms("dndeta");
1751
1752     fdNdEtaAnalysis[i]->Finish(baseCorrection, 0.3, AlidNdEtaCorrection::kINEL);
1753     file->cd();
1754     fdNdEtaAnalysis[i]->SaveHistograms();
1755
1756     TH1* hist = fdNdEtaAnalysis[i]->GetdNdEtaHistogram(0);
1757     hist->SetLineColor(colors[i]);
1758     hist->SetMarkerColor(colors[i]);
1759     hist->SetMarkerStyle(markers[i]+1);
1760     hist->DrawCopy((i == 0) ? "" : "SAME");
1761     legend->AddEntry(hist, name);
1762   }
1763
1764   legend->Draw();
1765 }