removing checking of vertex resolution in getvertex
[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 void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char_t* correctionFileName="correction_mapprocess-types.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   // correctionTarget is of type AlidNdEtaCorrection::CorrectionType
742   //    kINEL = 3
743   //    kNSD = 4
744
745   loadlibs();
746
747   const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
748
749
750
751   const Char_t* changes[]  = { "pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25"};
752   Float_t scalesND[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0 };
753   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};
754   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};
755   Int_t nChanges = 17;
756
757   /*
758   const Char_t* changes[]  = { "pythia", "qgsm", "phojet"};
759   Float_t scalesND[] = {1.0, 1.10, 1.11};
760   Float_t scalesSD[] = {1.0, 0.69, 0.86};
761   Float_t scalesDD[] = {1.0, 0.98, 0.61};
762   Int_t nChanges = 3;
763   */
764   
765   // cross section from Pythia
766   // 14 TeV!
767   Float_t sigmaND = 55.2;
768   Float_t sigmaDD = 9.78;
769   Float_t sigmaSD = 14.30;
770
771   // standard correction
772   TFile::Open(correctionFileName);
773   AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
774   correctionStandard->LoadHistograms();
775
776   // dont take vertexreco from this one
777   correctionStandard->GetVertexRecoCorrection()->Reset();
778   // dont take triggerbias from this one
779   correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
780   correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
781   correctionStandard->GetTriggerBiasCorrectionND()->Reset();
782
783   AlidNdEtaCorrection* corrections[100];
784   TH1F* hRatios[100];
785
786   Int_t counter = 0;
787   for (Int_t j=0; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)
788
789     for (Int_t i=0; i<nChanges; i++) {
790       TFile::Open(correctionFileName);
791
792       TString name;
793       name.Form("dndeta_correction_syst_%s_%s", typeName[j], changes[i]);
794       AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
795       current->LoadHistograms("dndeta_correction");
796       current->Reset();
797
798       name.Form("dndeta_correction_ND");
799       AlidNdEtaCorrection* dNdEtaCorrectionND = new AlidNdEtaCorrection(name,name);
800       dNdEtaCorrectionND->LoadHistograms();
801       name.Form("dndeta_correction_DD");
802       AlidNdEtaCorrection* dNdEtaCorrectionDD = new AlidNdEtaCorrection(name,name);
803       dNdEtaCorrectionDD->LoadHistograms();
804       name.Form("dndeta_correction_SD");
805       AlidNdEtaCorrection* dNdEtaCorrectionSD = new AlidNdEtaCorrection(name,name);
806       dNdEtaCorrectionSD->LoadHistograms();
807
808       // calculating relative
809       Float_t nd = 100 * sigmaND/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
810       Float_t dd = 100 * (scalesDD[i]*sigmaDD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
811       Float_t sd = 100 * (scalesSD[i]*sigmaSD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
812
813       printf(Form("%s : ND=%.2f\%, DD=%.2f\%, SD=%.2f\% \n",changes[i],nd,dd,sd));
814       current->SetTitle(Form("ND=%.2f\%,DD=%.2f\%,SD=%.2f\%",nd,dd,sd));
815       current->SetTitle(name);
816
817       // scale
818       if (j == 0 || j == 2)
819       {
820         dNdEtaCorrectionND->GetVertexRecoCorrection()->Scale(scalesND[i]);
821         dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scalesDD[i]);
822         dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scalesSD[i]);
823       }
824       if (j == 1 || j == 2)
825       {
826         dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->Scale(scalesND[i]);
827         dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scalesDD[i]);
828         dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scalesSD[i]);
829
830         dNdEtaCorrectionND->GetTriggerBiasCorrectionNSD()->Scale(scalesND[i]);
831         dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scalesDD[i]);
832         dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scalesSD[i]);
833
834         dNdEtaCorrectionND->GetTriggerBiasCorrectionND()->Scale(scalesND[i]);
835         dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scalesDD[i]);
836         dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scalesSD[i]);
837       }
838
839       //clear track in correction
840       dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
841       dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
842       dNdEtaCorrectionSD->GetTrack2ParticleCorrection()->Reset();
843
844       TList collection;
845       collection.Add(correctionStandard);
846       collection.Add(dNdEtaCorrectionND);
847       collection.Add(dNdEtaCorrectionDD);
848       collection.Add(dNdEtaCorrectionSD);
849
850       current->Merge(&collection);
851       current->Finish();
852
853       corrections[counter] = current;
854
855       // now correct dNdeta distribution with modified correction map
856       TFile::Open(analysisFileName);
857
858       dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
859       fdNdEtaAnalysis->LoadHistograms();
860
861       fdNdEtaAnalysis->Finish(current, 0.3, correctionTarget, Form("%d %d", j, i));
862
863       name = "ratio";
864       if (j==0) name.Append("_vetexReco_");
865       if (j==1) name.Append("_triggerBias_");
866       if (j==2) name.Append("_vertexReco_triggerBias_");
867       name.Append(changes[i]);
868
869       hRatios[counter] = (TH1F*) fdNdEtaAnalysis->GetdNdEtaHistogram()->Clone(name);
870
871       name.Append(Form(" (DD #times %0.2f, SD #times %0.2f)",scalesDD[i],scalesSD[i]));
872       hRatios[counter]->SetTitle(name.Data());
873       hRatios[counter]->SetYTitle("ratio (Pythia x-sections)/(changed x-sections)");
874
875       if (counter > 0)
876         hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);
877
878       delete fdNdEtaAnalysis;
879
880       counter++;
881     }
882   }
883
884   TFile* fout = new TFile(outputFileName,"RECREATE");
885
886   // to make everything consistent
887   hRatios[0]->Divide(hRatios[0],hRatios[0],1,1);
888
889   for (Int_t i=0; i<nChanges * 3; i++)
890   {
891     corrections[i]->SaveHistograms();
892     hRatios[i]->Write();
893   }
894
895   fout->Write();
896   fout->Close();
897 }
898
899
900 DrawVertexRecoSyst() {
901   // Draws the ratio of the dN/dEta obtained with changed SD and DD
902   // cross-sections vertex reco corrections to the dN/dEta obtained
903   // from the standard pythia cross-sections 
904   //
905   // The files with the vertex reco corrections for different
906   // processes (and with the other standard corrections) are expected
907   // to have the names "analysis_esd_X.root", where the Xs are defined
908   // in the array changes below.
909
910   Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless"};
911   Char_t* descr[]  =   {"",
912                         "#sigma_{DD}' = 1.5#times#sigma_{DD}",
913                         "#sigma_{DD}' = 0.5#times#sigma_{DD}",
914                         "#sigma_{SD}' = 1.5#times#sigma_{SD}",
915                         "#sigma_{SD}' = 0.5#times#sigma_{SD}",
916                         "#sigma_{D}' = 1.5#times#sigma_{D}",
917                         "#sigma_{D}' = 0.5#times#sigma_{D}"};
918
919   Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.5, 0.5, 1.5, 0.5};
920   Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.5, 0.5};
921
922   Int_t markers[] = {20,20,21,22,23,28,29};
923   Int_t colors[]  = {1,2,3,4,6,8,102};
924
925   // cross section from Pythia
926   Float_t sigmaND = 55.2;
927   Float_t sigmaDD = 9.78;
928   Float_t sigmaSD = 14.30;
929
930   TH1F* dNdEta[7];
931   TH1F* ratios[7];
932
933   TFile* fin;
934
935   for (Int_t i=0; i<7; i++) {
936     // calculating relative
937     fin = TFile::Open(Form("analysis_esd_%s.root",changes[i]));
938
939     dNdEta[i] = (TH1F*)(fin->Get("dndeta/dndeta_dNdEta_corrected_2"))->Clone();
940
941     for (Int_t b=0; b<dNdEta[i]->GetNbinsX(); b++) {
942       if (TMath::Abs(dNdEta[i]->GetBinCenter(b))>0.9) {
943         dNdEta[i]->SetBinContent(b,0);
944         dNdEta[i]->SetBinError(b,0);
945       }
946     }
947
948     dNdEta[i]->Rebin(4);
949
950     dNdEta[i]->SetLineWidth(2);
951     dNdEta[i]->SetLineColor(colors[i]);
952     dNdEta[i]->SetMarkerStyle(markers[i]);
953     dNdEta[i]->SetMarkerSize(0.9);
954     dNdEta[i]->SetMarkerColor(colors[i]);
955
956     ratios[i] = (TH1F*)dNdEta[i]->Clone("ratio");
957     ratios[i]->Divide(ratios[i],dNdEta[0],1,1,"B");
958     
959     ratios[i]->SetName(changes[i]);
960     ratios[i]->SetTitle(changes[i]);
961   }
962   
963   //##########################################################
964   
965   gStyle->SetOptStat(0);
966   gStyle->SetOptTitle(0);
967   gStyle->SetOptFit(0);
968
969   gStyle->SetTextSize(0.2);
970   gStyle->SetTitleSize(0.05,"xyz");
971   gStyle->SetLabelOffset(0.01, "xyz");
972
973
974   gStyle->SetTitleOffset(1.2, "y");
975   gStyle->SetTitleOffset(1.2, "x");
976   gStyle->SetEndErrorSize(0.0);
977
978   //##############################################
979
980   //making canvas and pads
981   TCanvas *c = new TCanvas(Form("vertex_reco_syst_%s", plot), Form("vertex_reco_syst_%s", plot),600,500);
982
983   TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
984
985   p1->SetBottomMargin(0.15);
986   p1->SetTopMargin(0.03);
987   p1->SetLeftMargin(0.15);
988   p1->SetRightMargin(0.03);
989
990   p1->SetGridx();
991   p1->SetGridy();
992
993   p1->Draw();
994   p1->cd();
995   
996   
997   TH2F* null = new TH2F("","",100,-1.5,1.5,100,0.9601,1.099);
998   null->SetXTitle("#eta");
999   null->GetXaxis()->CenterTitle(kTRUE);
1000   null->SetYTitle("dN/d#eta / dN/d#eta_{pythia}");
1001   null->GetYaxis()->CenterTitle(kTRUE);
1002   null->Draw();
1003
1004   for (Int_t i=1; i<7; i++) 
1005     ratios[i]->Draw("same");
1006
1007   TLegend* legend = new TLegend(0.6, 0.6, 0.95, 0.95);
1008   legend->SetFillColor(0);
1009
1010   TLatex* text[7];
1011   for (Int_t i=1; i<7; i++) {
1012     legend->AddEntry(dNdEta[i], descr[i]);
1013   }
1014
1015   legend->Draw();
1016   //text(0.2,0.88,"Effect of changing",0.045,1,kTRUE);
1017   //text(0.2,0.83,"relative cross-sections",0.045,1,kTRUE);
1018   //text(0.2,0.78,"(vertex reconstruction corr.)",0.043,13,kTRUE);
1019
1020   c->SaveAs(Form("%s.gif", c->GetName()));
1021   c->SaveAs(Form("%s.eps", c->GetName()));
1022 }
1023
1024
1025
1026 DrawTriggerEfficiency(Char_t* fileName) {
1027
1028   gStyle->SetOptStat(0);
1029   gStyle->SetOptTitle(0);
1030   gStyle->SetOptFit(0);
1031
1032   gStyle->SetTextSize(0.04);
1033   gStyle->SetTitleSize(0.05,"xyz");
1034   //gStyle->SetTitleFont(133, "xyz");
1035   //gStyle->SetLabelFont(133, "xyz");
1036   //gStyle->SetLabelSize(17, "xyz");
1037   gStyle->SetLabelOffset(0.01, "xyz");
1038
1039   gStyle->SetTitleOffset(1.1, "y");
1040   gStyle->SetTitleOffset(1.1, "x");
1041   gStyle->SetEndErrorSize(0.0);
1042
1043   //##############################################
1044
1045   //making canvas and pads
1046   TCanvas *c = new TCanvas("trigger_eff", "",600,500);
1047
1048   TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
1049
1050   p1->SetBottomMargin(0.15);
1051   p1->SetTopMargin(0.03);
1052   p1->SetLeftMargin(0.15);
1053   p1->SetRightMargin(0.03);
1054   
1055   p1->SetGridx();
1056   p1->SetGridy();
1057
1058   p1->Draw();
1059   p1->cd();
1060
1061   gSystem->Load("libPWG0base");
1062
1063   AlidNdEtaCorrection* corrections[4];
1064   AliCorrectionMatrix2D* triggerBiasCorrection[4];
1065
1066   TH1F* hTriggerEffInv[4];
1067   TH1F* hTriggerEff[4];
1068
1069   Char_t* names[] = {"triggerBiasND", "triggerBiasDD", "triggerBiasSD", "triggerBiasINEL"};
1070   
1071   for (Int_t i=0; i<4; i++) {
1072     corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1073     corrections[i]->LoadHistograms(fileName, names[i]);    
1074
1075     triggerBiasCorrection[i] = corrections[i]->GetTriggerBiasCorrectionINEL();
1076
1077     
1078     hTriggerEffInv[i] = triggerBiasCorrection[i]->Get1DCorrection();
1079     hTriggerEff[i]    = (TH1F*)hTriggerEffInv[i]->Clone();
1080     
1081     for (Int_t b=0; b<=hTriggerEff[i]->GetNbinsX(); b++) {
1082       hTriggerEff[i]->SetBinContent(b,1);
1083       hTriggerEff[i]->SetBinError(b,0);
1084     }
1085     hTriggerEff[i]->Divide(hTriggerEff[i],hTriggerEffInv[i]);
1086     hTriggerEff[i]->Scale(100);
1087   }
1088
1089   Int_t colors[] = {2,3,4,1};
1090   Float_t effs[4];
1091   for (Int_t i=0; i<4; i++) {
1092     hTriggerEff[i]->Fit("pol0","","",-20,20);
1093     effs[i] = ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->GetParameter(0);
1094     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineWidth(1);
1095     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineStyle(2);
1096     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineColor(colors[i]);
1097     cout << effs[i] << endl;
1098   }
1099
1100
1101   Char_t* text[] = {"ND", "DD", "SD", "INEL"};
1102   TLatex* latex[4];
1103
1104   TH2F* null = new TH2F("","",100,-25,35,100,0,110);
1105   null->SetXTitle("Vertex z [cm]");
1106   null->GetXaxis()->CenterTitle(kTRUE);
1107   null->SetYTitle("Trigger efficiency [%]");
1108   null->GetYaxis()->CenterTitle(kTRUE);
1109   null->Draw();
1110
1111
1112   for (Int_t i=0; i<4; i++) {
1113     hTriggerEff[i]->SetLineWidth(2);
1114     hTriggerEff[i]->SetLineColor(colors[i]);
1115
1116     hTriggerEff[i]->Draw("same");
1117
1118     latex[i] = new TLatex(22,effs[i]-1.5, Form("%s (%0.1f)",text[i],effs[i]));
1119     latex[i]->SetTextColor(colors[i]);
1120     latex[i]->Draw();
1121   }
1122   
1123 }
1124
1125
1126 DrawSpectraPID(Char_t* fileName) {
1127
1128   gSystem->Load("libPWG0base");
1129
1130   Char_t* names[] = {"correction_0", "correction_1", "correction_2", "correction_3"};
1131   AlidNdEtaCorrection* corrections[4];
1132   AliCorrectionMatrix3D* trackToPartCorrection[4];
1133
1134   TH1F* measuredPt[4];
1135   TH1F* generatedPt[4];
1136   TH1F* ratioPt[4];
1137
1138   for (Int_t i=0; i<4; i++) {
1139     corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1140     corrections[i]->LoadHistograms(fileName, names[i]);    
1141
1142     trackToPartCorrection[i] = corrections[i]->GetTrack2ParticleCorrection();
1143
1144     Int_t binX1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(-10);
1145     Int_t binX2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(10);
1146     Int_t binY1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(-1);
1147     Int_t binY2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(1);
1148
1149     measuredPt[i]  = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->ProjectionZ(Form("m_%d",i),binX1,binX2,binY1,binY2);
1150     generatedPt[i] = (TH1F*)trackToPartCorrection[i]->GetGeneratedHistogram()->ProjectionZ(Form("g_%d",i),binX1,binX2,binY1,binY2);
1151     ratioPt[i] = (TH1F*)generatedPt[i]->Clone(Form("r_%d",i));
1152     ratioPt[i]->Divide(measuredPt[i], generatedPt[i], 1,1,"B");
1153   }
1154   
1155   ratioPt[0]->Draw();
1156   
1157   for (Int_t i=0; i<3; i++) {
1158     ratioPt[i]->SetLineColor(i+1);
1159     ratioPt[i]->SetLineWidth(2);
1160     
1161     ratioPt[i]->Draw("same");
1162     
1163   }
1164
1165   return;
1166   measuredPt[0]->SetLineColor(2);
1167   measuredPt[0]->SetLineWidth(5);
1168
1169   measuredPt[0]->Draw();
1170   generatedPt[0]->Draw("same");
1171 }
1172
1173 void changePtSpectrum(const char* fileName = "analysis_mc.root")
1174 {
1175   TFile* file = TFile::Open(fileName);
1176   TH1F* hist = dynamic_cast<TH1F*> (file->Get("dndeta_check_pt"));
1177
1178   //hist->Rebin(3);
1179   //hist->Scale(1.0/3);
1180
1181   TH1F* clone1 = dynamic_cast<TH1F*> (hist->Clone("clone1"));
1182   TH1F* clone2 = dynamic_cast<TH1F*> (hist->Clone("clone2"));
1183
1184   TH1F* scale1 =  dynamic_cast<TH1F*> (hist->Clone("scale1"));
1185   TH1F* scale2 =  dynamic_cast<TH1F*> (hist->Clone("scale2"));
1186
1187   Float_t ptCutOff = 0.3;
1188
1189   for (Int_t i=1; i <= hist->GetNbinsX(); ++i)
1190   {
1191     if (hist->GetBinCenter(i) > ptCutOff)
1192     {
1193       scale1->SetBinContent(i, 1);
1194       scale2->SetBinContent(i, 1);
1195     }
1196     else
1197     {
1198       // 90 % at pt = 0, 0% at pt = ptcutoff
1199       scale1->SetBinContent(i, 1 - (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * 0.3);
1200
1201       // 110% at pt = 0, ...
1202       scale2->SetBinContent(i, 1 + (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * 0.3);
1203     }
1204     scale1->SetBinError(i, 0);
1205     scale2->SetBinError(i, 0);
1206   }
1207
1208   new TCanvas;
1209
1210   scale1->Draw();
1211   scale2->SetLineColor(kRed);
1212   scale2->Draw("SAME");
1213
1214   clone1->Multiply(scale1);
1215   clone2->Multiply(scale2);
1216
1217   Prepare1DPlot(hist);
1218   Prepare1DPlot(clone1);
1219   Prepare1DPlot(clone2);
1220
1221   /*hist->SetMarkerStyle(markers[0]);
1222   clone1->SetMarkerStyle(markers[0]);
1223   clone2->SetMarkerStyle(markers[0]);*/
1224
1225   hist->SetTitle(";p_{T} in GeV/c;dN_{ch}/dp_{T} in c/GeV");
1226   hist->GetXaxis()->SetRangeUser(0, 0.7);
1227   hist->GetYaxis()->SetRangeUser(0.01, clone2->GetMaximum() * 1.1);
1228   hist->GetYaxis()->SetTitleOffset(1);
1229
1230   TCanvas* canvas = new TCanvas;
1231   gPad->SetBottomMargin(0.12);
1232   hist->Draw("H");
1233   clone1->SetLineColor(kRed);
1234   clone1->Draw("HSAME");
1235   clone2->SetLineColor(kBlue);
1236   clone2->Draw("HSAME");
1237   hist->Draw("HSAME");
1238
1239   Float_t fraction =  hist->Integral(hist->GetXaxis()->FindBin(ptCutOff), hist->GetNbinsX()) / hist->Integral(1, hist->GetNbinsX());
1240   Float_t fraction1 = clone1->Integral(clone1->GetXaxis()->FindBin(ptCutOff), clone1->GetNbinsX()) / clone1->Integral(1, clone1->GetNbinsX());
1241   Float_t fraction2 = clone2->Integral(clone2->GetXaxis()->FindBin(ptCutOff), clone2->GetNbinsX()) / clone2->Integral(1, clone2->GetNbinsX());
1242
1243   printf("%f %f %f\n", fraction, fraction1, fraction2);
1244   printf("Rel. %f %f\n", fraction1 / fraction, fraction2 / fraction);
1245
1246   canvas->SaveAs("changePtSpectrum.gif");
1247   canvas->SaveAs("changePtSpectrum.eps");
1248 }
1249
1250 void FractionBelowPt()
1251 {
1252   gSystem->Load("libPWG0base");
1253
1254   AlidNdEtaCorrection* fdNdEtaCorrection[4];
1255
1256   TFile::Open("systematics.root");
1257
1258   for (Int_t i=0; i<4; ++i)
1259   {
1260     TString name;
1261     name.Form("correction_%d", i);
1262     fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
1263     fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
1264   }
1265
1266   Double_t geneCount[5];
1267   Double_t measCount[5];
1268   geneCount[4] = 0;
1269   measCount[4] = 0;
1270
1271   for (Int_t i=0; i<4; ++i)
1272   {
1273     TH3F* hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1274     geneCount[i] = hist->Integral(hist->GetXaxis()->FindBin(-10), hist->GetXaxis()->FindBin(10),
1275                                   hist->GetYaxis()->FindBin(-0.8), hist->GetYaxis()->FindBin(0.8),
1276                                   1, hist->GetZaxis()->FindBin(0.3));
1277
1278     hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1279     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));
1280
1281     geneCount[4] += geneCount[i];
1282     measCount[4] += measCount[i];
1283
1284     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]);
1285   }
1286
1287   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]);
1288
1289   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]);
1290 }
1291
1292
1293 mergeCorrectionsMisalignment(Char_t* alignedFile = "correction_map_aligned.root",
1294                                              Char_t* misalignedFile = "correction_map_misaligned.root",
1295                                              Char_t* outputFileName="correction_map_misaligned_single.root")
1296 {
1297   //
1298   // from the aligned and misaligned corrections, 3 new corrections are created
1299   // in these new corrections only one of the corrections (track2particle, vertex, trigger)
1300   // is taken from the misaligned input to allow study of the effect on the different
1301   // corrections
1302
1303   gSystem->Load("libPWG0base");
1304
1305   const Char_t* typeName[] = { "track2particle", "vertex", "trigger" };
1306
1307   AlidNdEtaCorrection* corrections[3];
1308   for (Int_t j=0; j<3; j++) { // j = 0 (track2particle), j = 1 (vertex), j = 2 (trigger)
1309     AlidNdEtaCorrection* alignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1310     alignedCorrection->LoadHistograms(alignedFile);
1311
1312     AlidNdEtaCorrection* misalignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1313     misalignedCorrection->LoadHistograms(misalignedFile);
1314
1315     TString name;
1316     name.Form("dndeta_correction_alignment_%s", typeName[j]);
1317     AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
1318
1319     switch (j)
1320     {
1321       case 0:
1322         alignedCorrection->GetTrack2ParticleCorrection()->Reset();
1323         misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1324         misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1325         misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1326         misalignedCorrection->GetVertexRecoCorrection()->Reset();
1327         break;
1328
1329       case 1:
1330         misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1331         misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1332         misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1333         misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1334         alignedCorrection->GetVertexRecoCorrection()->Reset();
1335         break;
1336
1337       case 2:
1338         misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1339         alignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1340         alignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1341         alignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1342         misalignedCorrection->GetVertexRecoCorrection()->Reset();
1343         break;
1344
1345       default:
1346         return;
1347     }
1348
1349     TList collection;
1350     collection.Add(misalignedCorrection);
1351     collection.Add(alignedCorrection);
1352
1353     current->Merge(&collection);
1354     current->Finish();
1355
1356     corrections[j] = current;
1357   }
1358
1359   TFile* fout = new TFile(outputFileName, "RECREATE");
1360
1361   for (Int_t i=0; i<3; i++)
1362     corrections[i]->SaveHistograms();
1363
1364   fout->Write();
1365   fout->Close();
1366 }
1367
1368
1369 void DrawdNdEtaDifferencesAlignment()
1370 {
1371   TH1* hists[5];
1372
1373   TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
1374   legend->SetFillColor(0);
1375
1376   TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
1377   canvas->Divide(2, 1);
1378
1379   canvas->cd(1);
1380
1381   for (Int_t i=0; i<5; ++i)
1382   {
1383     hists[i] = 0;
1384     TFile* file = 0;
1385     TString title;
1386
1387     switch(i)
1388     {
1389       case 0 : file = TFile::Open("systematics_misalignment_aligned.root"); title = "aligned"; break;
1390       case 1 : file = TFile::Open("systematics_misalignment_misaligned.root"); title = "fully misaligned"; break;
1391       case 2 : file = TFile::Open("systematics_misalignment_track2particle.root"); title = "only track2particle"; break;
1392       case 3 : file = TFile::Open("systematics_misalignment_vertex.root"); title = "only vertex rec."; break;
1393       case 4 : file = TFile::Open("systematics_misalignment_trigger.root"); title = "only trigger bias"; break;
1394       default: return;
1395     }
1396
1397     if (file)
1398     {
1399       hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
1400       hists[i]->SetTitle("");
1401
1402       Prepare1DPlot(hists[i], kFALSE);
1403       hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
1404       hists[i]->GetYaxis()->SetRangeUser(6, 7);
1405       hists[i]->SetLineWidth(1);
1406       hists[i]->SetLineColor(colors[i]);
1407       hists[i]->SetMarkerColor(colors[i]);
1408       hists[i]->SetMarkerStyle(markers[i]);
1409       hists[i]->GetXaxis()->SetLabelOffset(0.015);
1410       hists[i]->GetYaxis()->SetTitleOffset(1.5);
1411       gPad->SetLeftMargin(0.12);
1412       hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
1413
1414       legend->AddEntry(hists[i], title);
1415       hists[i]->SetTitle(title);
1416     }
1417   }
1418   legend->Draw();
1419
1420   canvas->cd(2);
1421   gPad->SetLeftMargin(0.14);
1422
1423   TLegend* legend2 = new TLegend(0.63, 0.73, 0.98, 0.98);
1424   legend2->SetFillColor(0);
1425
1426   for (Int_t i=1; i<5; ++i)
1427   {
1428     if (hists[i])
1429     {
1430       legend2->AddEntry(hists[i]);
1431
1432       hists[i]->Divide(hists[0]);
1433       hists[i]->SetTitle("b)");
1434       hists[i]->GetYaxis()->SetRangeUser(0.9, 1.1);
1435       hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
1436       hists[i]->GetYaxis()->SetTitleOffset(1.8);
1437       hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
1438     }
1439   }
1440
1441   legend2->Draw();
1442
1443   canvas->SaveAs("misalignment_result.eps");
1444   canvas->SaveAs("misalignment_result.gif");
1445 }
1446
1447 void EvaluateMultiplicityEffect()
1448 {
1449   gSystem->Load("libPWG0base");
1450
1451   AlidNdEtaCorrection* fdNdEtaCorrectionLow[4];
1452   AlidNdEtaCorrection* fdNdEtaCorrectionHigh[4];
1453
1454   TFile::Open("systematics-low-multiplicity.root");
1455
1456   for (Int_t i=0; i<4; ++i)
1457   {
1458     TString name;
1459     name.Form("correction_%d", i);
1460     fdNdEtaCorrectionLow[i] = new AlidNdEtaCorrection(name, name);
1461     fdNdEtaCorrectionLow[i]->LoadHistograms("systematics-low-multiplicity.root", name);
1462   }
1463
1464   TList list;
1465   for (Int_t i=1; i<4; ++i)
1466     list.Add(fdNdEtaCorrectionLow[i]);
1467
1468   fdNdEtaCorrectionLow[0]->Merge(&list);
1469   fdNdEtaCorrectionLow[0]->Finish();
1470
1471   TFile::Open("systematics-high-multiplicity.root");
1472
1473   for (Int_t i=0; i<4; ++i)
1474   {
1475     TString name;
1476     name.Form("correction_%d", i);
1477     fdNdEtaCorrectionHigh[i] = new AlidNdEtaCorrection(name, name);
1478     fdNdEtaCorrectionHigh[i]->LoadHistograms("systematics-high-multiplicity.root", name);
1479   }
1480
1481   TList list2;
1482   for (Int_t i=1; i<4; ++i)
1483     list2.Add(fdNdEtaCorrectionHigh[i]);
1484
1485   fdNdEtaCorrectionHigh[0]->Merge(&list2);
1486   fdNdEtaCorrectionHigh[0]->Finish();
1487
1488   TH1F* outputLow = new TH1F("Track2ParticleLow", "Track2Particle at low multiplicity", 200, 0, 2);
1489   TH1F* outputHigh = new TH1F("Track2ParticleHigh", "Track2Particle at high multiplicity", 200, 0, 2);
1490
1491   TH3F* hist = fdNdEtaCorrectionLow[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1492   TH3F* hist2 = fdNdEtaCorrectionHigh[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1493   for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
1494     for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
1495       for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
1496       //for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
1497       {
1498         if (hist->GetBinContent(x, y, z) > 0)
1499           outputLow->Fill(hist->GetBinContent(x, y, z));
1500         //if (hist->GetBinContent(x, y, z) == 1)
1501         //  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));
1502
1503         if (hist2->GetBinContent(x, y, z) > 0)
1504           outputHigh->Fill(hist2->GetBinContent(x, y, z));
1505       }
1506
1507   TCanvas* canvas = new TCanvas("EvaluateMultiplicityEffect", "EvaluateMultiplicityEffect", 1000, 500);
1508   canvas->Divide(2, 1);
1509
1510   canvas->cd(1);
1511   outputLow->Draw();
1512   outputLow->Fit("gaus", "0");
1513   outputLow->GetFunction("gaus")->SetLineColor(2);
1514   outputLow->GetFunction("gaus")->DrawCopy("SAME");
1515
1516   canvas->cd(2);
1517   outputHigh->Draw();
1518   outputHigh->Fit("gaus", "0");
1519   outputHigh->GetFunction("gaus")->DrawCopy("SAME");
1520
1521   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1522 }
1523
1524 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")
1525 {
1526     gSystem->Load("libPWG0base");
1527         
1528         AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1529     dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
1530     
1531     dNdEtaCorrection->GetTrack2ParticleCorrection()->PlotBinErrors(xmin, xmax, ymin, ymax, zmin, zmax)->Draw();
1532 }
1533
1534
1535 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")
1536 {
1537   gSystem->Load("libPWG0base");
1538
1539   TFile* file = TFile::Open(output, "RECREATE");
1540
1541   const Int_t max = 5;
1542   dNdEtaAnalysis* fdNdEtaAnalysis[5];
1543
1544   new TCanvas;
1545   TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
1546   legend->SetFillColor(0);
1547
1548   for (Int_t i = 0; i < max; ++i)
1549   {
1550     TFile::Open(baseCorrectionMapFile);
1551     AlidNdEtaCorrection* baseCorrection = new AlidNdEtaCorrection(baseCorrectionMapFolder, baseCorrectionMapFolder);
1552     baseCorrection->LoadHistograms();
1553
1554     AlidNdEtaCorrection::CorrectionType correctionType = AlidNdEtaCorrection::kNone;
1555     const char* name = 0;
1556
1557     TFile::Open(changedCorrectionMapFile);
1558     switch (i)
1559     {
1560       case 0 :
1561         name = "default";
1562         break;
1563
1564       case 1 :
1565         baseCorrection->GetTrack2ParticleCorrection()->LoadHistograms(Form("%s/Track2Particle", changedCorrectionMapFolder));
1566         name = "Track2Particle";
1567         break;
1568
1569       case 2 :
1570         baseCorrection->GetVertexRecoCorrection()->LoadHistograms(Form("%s/VertexReconstruction", changedCorrectionMapFolder));
1571         name = "VertexReco";
1572         break;
1573
1574       case 3 :
1575         baseCorrection->GetTriggerBiasCorrectionINEL()->LoadHistograms(Form("%s/TriggerBias_MBToINEL", changedCorrectionMapFolder));
1576         name = "TriggerBias_MBToINEL";
1577         break;
1578
1579       case 4 :
1580         baseCorrection->LoadHistograms(changedCorrectionMapFolder);
1581         name = "all";
1582         break;
1583
1584       default: return;
1585     }
1586
1587     TFile::Open(dataFile);
1588     fdNdEtaAnalysis[i] = new dNdEtaAnalysis(name, name);
1589     fdNdEtaAnalysis[i]->LoadHistograms("dndeta");
1590
1591     fdNdEtaAnalysis[i]->Finish(baseCorrection, 0.3, AlidNdEtaCorrection::kINEL);
1592     file->cd();
1593     fdNdEtaAnalysis[i]->SaveHistograms();
1594
1595     TH1* hist = fdNdEtaAnalysis[i]->GetdNdEtaHistogram(0);
1596     hist->SetLineColor(colors[i]);
1597     hist->SetMarkerColor(colors[i]);
1598     hist->SetMarkerStyle(markers[i]+1);
1599     hist->DrawCopy((i == 0) ? "" : "SAME");
1600     legend->AddEntry(hist, name);
1601   }
1602
1603   legend->Draw();
1604 }