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