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