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