Removing dummy if
[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: // + 30% kaons
508     case 3: // - 30% kaons
509     case 4: // + 30% protons
510     case 5: // - 30% protons
511     case 6: // + 30% kaons + 30% protons
512     case 7: // - 30% kaons - 30% protons
513     case 8: // + 30% kaons - 30% protons
514     case 9: // - 30% kaons + 30% protons
515     case 10: // + 30% others
516     case 11: // - 30% others
517       TString* str = new TString;
518       if (index < 6)
519       {
520         Int_t correctionIndex = index / 2;
521         Double_t scaleFactor = (index % 2 == 0) ? 1.3 : 0.7;
522   
523         fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
524         str->Form("%s%s", (correctionIndex == 0) ? "Pi" : ((correctionIndex == 1) ? "K" : (correctionIndex == 2) ? "p" : "others"), (index % 2 == 0) ? "Boosted" : "Reduced");
525       }
526       else if (index < 10)
527       {
528         Double_t scaleFactor = (index % 2 == 0) ? 1.3 : 0.7;
529         fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
530         str->Form("%s%s", "K", (scaleFactor > 1) ? "Boosted" : "Reduced");
531         
532         if (index >= 8)
533           scaleFactor = (index % 2 == 0) ? 0.3 : 1.7;
534         fdNdEtaCorrection[2]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
535         *str += Form("%s%s", "p", (scaleFactor > 1) ? "Boosted" : "Reduced");
536       }
537       else
538       {
539         Double_t scaleFactor = (index % 2 == 0) ? 1.3 : 0.7;
540         fdNdEtaCorrection[3]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
541         str->Form("%s%s", "others", (scaleFactor > 1) ? "Boosted" : "Reduced");
542       }
543
544       return str->Data();
545       break;
546
547       // each species rated with pythia ratios
548     case 12: // + 50% pions
549     case 13: // - 50% pions
550     case 14: // + 50% kaons
551     case 15: // - 50% kaons
552     case 16: // + 50% protons
553     case 17: // - 50% protons
554       TH1** ptDists = DrawRatios("pythiaratios.root");
555       Int_t functionIndex = (index - 2) / 2;
556       Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
557       ptDists[functionIndex]->Scale(scaleFactor);
558
559       for (Int_t i=0; i<3; ++i)
560       {
561         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram(), ptDists[i]);
562         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram(), ptDists[i]);
563       }
564       TString* str = new TString;
565       str->Form("%s%s", (functionIndex == 0) ? "Pi" : ((functionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced");
566       return str->Data();
567       break;
568
569     case 999:
570       TF1* ptDependence = new TF1("simple", "x", 0, 100);
571       ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram(), ptDependence);
572       ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram(), ptDependence);
573       break;
574
575   }
576
577   return "noname";
578 }
579
580 void Composition()
581 {
582   loadlibs();
583
584   gSystem->Unlink("new_compositions.root");
585   gSystem->Unlink("new_compositions_analysis.root");
586   
587   const char* names[] = { "pi", "K", "p", "other" };
588   TH1* hRatios[20];
589
590   //backgroundEvents = 1162+434; // Michele for MB1, run 104892, 15.02.10
591   backgroundEvents = -1;    // use 0 bin from MC! for 2.36 TeV
592   
593   Printf("Subtracting %d background events!!!", backgroundEvents);
594   gSystem->Sleep(1000);
595   
596   Int_t nCompositions = 12;
597   Int_t counter = 0;
598   for (Int_t comp = 1; comp < nCompositions; ++comp)
599   {
600     AlidNdEtaCorrection* fdNdEtaCorrection[4];
601
602     TFile::Open("correction_mapparticle-species.root");
603
604     for (Int_t i=0; i<4; ++i)
605     {
606       TString name;
607       name.Form("dndeta_correction_%s", names[i]);
608       fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
609       fdNdEtaCorrection[i]->LoadHistograms();
610     }
611
612     const char* newName = ChangeComposition(fdNdEtaCorrection, comp);
613
614     Double_t geneCount[5];
615     Double_t measCount[5];
616     geneCount[4] = 0;
617     measCount[4] = 0;
618
619     for (Int_t i=0; i<4; ++i)
620     {
621       geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram()->Integral();
622       measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Integral();
623
624       geneCount[4] += geneCount[i];
625       measCount[4] += measCount[i];
626
627       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]);
628     }
629
630     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]);
631
632     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]);
633
634     TList* collection = new TList;
635
636     // skip "other" particle correction here
637     // with them has to be dealt differently, maybe just increasing the neutral particles...
638     for (Int_t i=1; i<4; ++i)
639       collection->Add(fdNdEtaCorrection[i]);
640
641     fdNdEtaCorrection[0]->Merge(collection);
642     fdNdEtaCorrection[0]->Finish();
643
644     delete collection;
645
646     // save everything
647     TFile* file = TFile::Open("new_compositions.root", "UPDATE");
648     fdNdEtaCorrection[0]->SetName(newName);
649     fdNdEtaCorrection[0]->SaveHistograms();
650     //file->Write();
651     file->Close();
652     
653     // correct dNdeta distribution with modified correction map
654     TFile::Open("analysis_esd_raw.root");
655
656     dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
657     fdNdEtaAnalysis->LoadHistograms();
658
659     fdNdEtaAnalysis->Finish(fdNdEtaCorrection[0], 0.2, 3, newName);
660     
661     hRatios[counter] = (TH1F*) fdNdEtaAnalysis->GetdNdEtaHistogram()->Clone(newName);
662     hRatios[counter]->SetTitle(newName);
663     hRatios[counter]->SetYTitle("dN_{ch}/d#eta ratio #frac{default composition}{modified composition}");
664
665     if (counter > 0)
666       hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);
667
668     file = TFile::Open("new_compositions_analysis.root", "UPDATE");
669     hRatios[counter]->Write();
670     file->Close();
671     
672     delete fdNdEtaAnalysis;
673
674     counter++;
675   }
676
677   /*
678   gROOT->ProcessLine(".L drawPlots.C");
679
680   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" };
681   const char* folderNames[] = { "PythiaRatios", "PiBoosted", "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced" };
682
683   Track2Particle1DComposition(fileNames, nCompositions, folderNames);
684   */
685 }
686
687
688 void drawSystematics()
689 {
690   //Secondaries();
691   //DrawDifferentSpecies();
692   //Composition();
693
694   Sigma2VertexSimulation();
695
696 }
697
698 void DrawdNdEtaDifferences()
699 {
700   TH1* hists[5];
701
702   TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
703   legend->SetFillColor(0);
704
705   TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
706   canvas->Divide(2, 1);
707
708   canvas->cd(1);
709
710   for (Int_t i=0; i<5; ++i)
711   {
712     hists[i] = 0;
713     TFile* file = 0;
714     TString title;
715
716     switch(i)
717     {
718       case 0 : file = TFile::Open("systematics_dndeta_reference.root"); title = "standard composition"; break;
719       case 1 : file = TFile::Open("systematics_dndeta_KBoosted.root"); title = "+ 50% kaons"; break;
720       case 2 : file = TFile::Open("systematics_dndeta_KReduced.root"); title = "- 50% kaons"; break;
721       case 3 : file = TFile::Open("systematics_dndeta_pBoosted.root"); title = "+ 50% protons"; break;
722       case 4 : file = TFile::Open("systematics_dndeta_pReduced.root"); title = "- 50% protons"; break;
723       default: return;
724     }
725
726     if (file)
727     {
728       hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
729       hists[i]->SetTitle("a)");
730
731       Prepare1DPlot(hists[i], kFALSE);
732       hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
733       hists[i]->GetYaxis()->SetRangeUser(6, 7);
734       hists[i]->SetLineColor(colors[i]);
735       hists[i]->SetMarkerColor(colors[i]);
736       hists[i]->SetMarkerStyle(markers[i]);
737       hists[i]->GetXaxis()->SetLabelOffset(0.015);
738       hists[i]->GetYaxis()->SetTitleOffset(1.5);
739       gPad->SetLeftMargin(0.12);
740       hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
741
742       legend->AddEntry(hists[i], title);
743       hists[i]->SetTitle(title);
744     }
745   }
746   legend->Draw();
747
748   canvas->cd(2);
749   gPad->SetLeftMargin(0.14);
750
751   TLegend* legend2 = new TLegend(0.73, 0.73, 0.98, 0.98);
752   legend2->SetFillColor(0);
753
754   for (Int_t i=1; i<5; ++i)
755   {
756     if (hists[i])
757     {
758       legend2->AddEntry(hists[i]);
759
760       hists[i]->Divide(hists[0]);
761       hists[i]->SetTitle("b)");
762       hists[i]->SetLineColor(colors[i-1]);
763       hists[i]->SetMarkerColor(colors[i-1]);
764       hists[i]->GetYaxis()->SetRangeUser(0.95, 1.05);
765       hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
766       hists[i]->GetYaxis()->SetTitleOffset(1.8);
767       hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
768     }
769   }
770
771   legend2->Draw();
772
773   canvas->SaveAs("particlecomposition_result_detail.gif");
774
775   TCanvas* canvas2 = new TCanvas("DrawdNdEtaDifferences2", "DrawdNdEtaDifferences2", 700, 500);
776
777   for (Int_t i=1; i<5; ++i)
778   {
779     if (hists[i])
780     {
781       hists[i]->SetTitle("");
782       hists[i]->GetYaxis()->SetTitleOffset(1.1);
783       hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
784     }
785   }
786
787   legend2->Draw();
788
789   canvas2->SaveAs("particlecomposition_result.gif");
790   canvas2->SaveAs("particlecomposition_result.eps");
791 }
792
793 void mergeCorrectionsWithDifferentCrosssections(Int_t origin, Int_t correctionTarget = 3, Char_t* correctionFileName="correction_mapprocess-types.root", const char* analysisFileName = "analysis_esd_raw.root", const Char_t* outputFileName=0) {
794   //
795   // Function used to merge standard corrections with vertex
796   // reconstruction corrections obtained by a certain mix of ND, DD
797   // and SD events.
798   //
799   // the dn/deta spectrum is corrected and the ratios
800   // (standard to changed x-section) of the different dN/deta
801   // distributions are saved to a file.
802   //
803   // correctionTarget is of type AlidNdEtaCorrection::CorrectionType
804   //    kINEL = 3
805   //    kNSD = 4
806   //    kOnePart = 6
807
808   if (outputFileName == 0)
809   {
810     if (correctionTarget == 3)
811       outputFileName = "systematics_vtxtrigger_compositions_inel.root";
812     if (correctionTarget == 4)
813       outputFileName = "systematics_vtxtrigger_compositions_nsd.root";
814     if (correctionTarget == 6)
815       outputFileName = "systematics_vtxtrigger_compositions_onepart.root";
816   }
817
818   loadlibs();
819
820   const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
821
822   //Karel:
823 //     fsd = 0.153 +- 0.031 (0.050 to take into account SD definition) --> change
824 //     fdd = 0.080 +- 0.050 --> change 
825 //     fnd = 0.767 +- 0.059 --> keep (error small)
826
827 //  const Char_t* changes[]  = { "pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdlessddmore25", "sdmoreddless25"};
828   //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 };
829   //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};
830   //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};
831   //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};
832   //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};
833 /*  Int_t nChanges = 9;
834
835   const Char_t* changes[]  = { "ua5","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless" };
836   Float_t scalesND[] = {0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767};
837   Float_t scalesDD[] = {0.080, 0.130, 0.030, 0.080, 0.080, 0.130, 0.030, 0.130, 0.030};
838   Float_t scalesSD[] = {0.153, 0.153, 0.153, 0.203, 0.103, 0.203, 0.103, 0.103, 0.203};*/
839   
840   Float_t ref_SD = -1;
841   Float_t ref_DD = -1;
842   Float_t ref_ND = -1;
843   
844   Float_t error_SD = -1;
845   Float_t error_DD = -1;
846   Float_t error_ND = -1;
847   
848   GetRelativeFractions(origin, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
849   
850   Printf("Using x-sections:\n SD: %f +- %f\n DD: %f +- %f\n ND: %f +- %f", ref_SD, error_SD, ref_DD, error_DD, ref_ND, error_ND);
851   
852   const Char_t* changes[]  = { "default","sdless","sdmore","ddless","ddmore", "dless", "dmore", "sdlessddmore", "sdmoreddless" };
853   Int_t nChanges = 9;
854   Float_t scalesSD[9];
855   Float_t scalesDD[9];
856   Float_t scalesND[9];
857   
858   if (1)
859   {
860     // sample 8 points on the error ellipse
861     for (Int_t i=0; i<9; i++)
862     {
863       Float_t factorSD = 0;
864       Float_t factorDD = 0;
865       
866       if (i > 0 && i < 3)
867         factorSD = (i % 2 == 0) ? 1 : -1;
868       else if (i >= 3 && i < 5)
869         factorDD = (i % 2 == 0) ? 1 : -1;
870       else if (i >= 5 && i < 9)
871       {
872         factorSD = ((i % 2 == 0) ? 1.0 : -1.0) / TMath::Sqrt(2);
873         if (i == 5 || i == 6)
874           factorDD = factorSD;
875         else
876           factorDD = -factorSD;
877       }
878       
879       scalesSD[i] = ref_SD + factorSD * error_SD;
880       scalesDD[i] = ref_DD + factorDD * error_DD;
881       scalesND[i] = 1.0 - scalesDD[i] - scalesSD[i];
882       
883       Printf("Case %d: SD: %f DD: %f ND: %f", i, scalesSD[i], scalesDD[i], scalesND[i]);
884     }
885   }
886   else
887   {
888     Printf("WARNING: Special treatment for ratios active");
889     gSystem->Sleep(1000);
890     
891     // constrained values by allowed changing of cross-sections
892     Float_t pythiaScaling = 0.224 / 0.189;
893
894     if (origin == 10)
895     {
896       // 900 GeV
897       for (Int_t i=0; i<9; i++)
898       {
899         scalesSD[i] = 15.3;
900         scalesDD[i] = 9.5;
901       }
902
903       scalesSD[1] = 15.7;
904       scalesSD[2] = 17.6;
905       scalesSD[3] = 13.5;
906       scalesSD[4] = 17.6;
907
908       scalesDD[5] = 15.5;
909       scalesDD[6] = 8.8;
910       scalesDD[7] = 13.8;
911       scalesDD[8] = 7.6;
912     }
913     else if (origin == 20)
914     {
915       // 2.36 TeV
916       pythiaScaling = 0.217 / 0.167;
917       
918       for (Int_t i=0; i<9; i++)
919       {
920         scalesSD[i] = 15.9;
921         scalesDD[i] = 10.7;
922       }
923
924       scalesSD[1] = 13.5;
925       scalesSD[2] = 15.2;
926       scalesSD[3] = 13.5;
927       scalesSD[4] = 17.6;
928
929       scalesDD[5] = 13.8;
930       scalesDD[6] = 7.6;
931       scalesDD[7] = 13.8;
932       scalesDD[8] = 7.6;
933     }
934     else
935       AliFatal("Not supported");
936
937     for (Int_t i=0; i<9; i++)
938     {
939       scalesSD[i] /= 100;
940       scalesSD[i] *= pythiaScaling;
941       scalesDD[i] /= 100;
942       scalesND[i] = 1.0 - scalesDD[i] - scalesSD[i];
943       Printf("Case %d: SD: %f DD: %f ND: %f", i, scalesSD[i], scalesDD[i], scalesND[i]);
944     }
945   }
946   
947   Int_t backgroundEvents = 0;
948   
949   //backgroundEvents = 1162+434; // Michele for MB1, run 104892, 15.02.10
950   //backgroundEvents = 6;          // Michele for V0AND, run 104892, 15.02.10
951   
952   //backgroundEvents = 4398+961;   // Michele for MB1, run 104824-52, 16.02.10
953   //backgroundEvents = 19;         // Michele for V0AND, run 104824-52, 16.02.10
954   
955   backgroundEvents = -1;    // use 0 bin from MC! for 2.36 TeV
956   
957   Printf("Subtracting %d background events!!!", backgroundEvents);
958   gSystem->Sleep(1000);
959   
960   /*
961   const Char_t* changes[]  = { "pythia", "qgsm", "phojet"};
962   Float_t scalesND[] = {1.0, 1.10, 1.11};
963   Float_t scalesSD[] = {1.0, 0.69, 0.86};
964   Float_t scalesDD[] = {1.0, 0.98, 0.61};
965   Int_t nChanges = 3;
966   */
967   
968   // cross section from Pythia
969   // 14 TeV!
970 //   Float_t sigmaND = 55.2;
971 //   Float_t sigmaDD = 9.78;
972 //   Float_t sigmaSD = 14.30;
973
974   // standard correction
975   TFile::Open(correctionFileName);
976   AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
977   correctionStandard->LoadHistograms();
978
979   // dont take vertexreco from this one
980   correctionStandard->GetVertexRecoCorrection()->Reset();
981   // dont take triggerbias from this one
982   correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
983   correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
984   correctionStandard->GetTriggerBiasCorrectionND()->Reset();
985   correctionStandard->GetTriggerBiasCorrectionOnePart()->Reset();
986
987   AlidNdEtaCorrection* corrections[100];
988   TH1F* hRatios[100];
989
990   Int_t counter = 0;
991   for (Int_t j=2; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)
992
993     for (Int_t i=0; i<nChanges; i++) {
994       TFile::Open(correctionFileName);
995
996       TString name;
997       name.Form("dndeta_correction_syst_%s_%s", typeName[j], changes[i]);
998       AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
999       current->LoadHistograms("dndeta_correction");
1000       current->Reset();
1001
1002       name.Form("dndeta_correction_ND");
1003       AlidNdEtaCorrection* dNdEtaCorrectionND = new AlidNdEtaCorrection(name,name);
1004       dNdEtaCorrectionND->LoadHistograms();
1005       name.Form("dndeta_correction_DD");
1006       AlidNdEtaCorrection* dNdEtaCorrectionDD = new AlidNdEtaCorrection(name,name);
1007       dNdEtaCorrectionDD->LoadHistograms();
1008       name.Form("dndeta_correction_SD");
1009       AlidNdEtaCorrection* dNdEtaCorrectionSD = new AlidNdEtaCorrection(name,name);
1010       dNdEtaCorrectionSD->LoadHistograms();
1011
1012       // calculating relative
1013       Float_t nd = dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1014       Float_t dd = dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1015       Float_t sd = dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1016       Float_t total = nd + dd + sd;
1017       
1018       nd /= total;
1019       sd /= total;
1020       dd /= total;
1021       
1022       Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
1023       
1024       Float_t scaleND = scalesND[i] / nd;
1025       Float_t scaleDD = scalesDD[i] / dd;
1026       Float_t scaleSD = scalesSD[i] / sd;
1027       
1028       Printf("ND=%.2f, DD=%.2f, SD=%.2f",scaleND, scaleDD, scaleSD);      
1029       
1030 /*      Float_t nd = 100 * sigmaND/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
1031       Float_t dd = 100 * (scalesDD[i]*sigmaDD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
1032       Float_t sd = 100 * (scalesSD[i]*sigmaSD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
1033
1034       printf(Form("%s : ND=%.2f\%, DD=%.2f\%, SD=%.2f\% \n",changes[i],nd,dd,sd));*/
1035       current->SetTitle(Form("ND=%.2f\%,DD=%.2f\%,SD=%.2f\%",scaleND,scaleDD,scaleSD));
1036       current->SetTitle(name);
1037
1038       // scale
1039       if (j == 0 || j == 2)
1040       {
1041         dNdEtaCorrectionND->GetVertexRecoCorrection()->Scale(scaleND);
1042         dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scaleDD);
1043         dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scaleSD);
1044       }
1045       if (j == 1 || j == 2)
1046       {
1047         dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->Scale(scaleND);
1048         dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scaleDD);
1049         dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scaleSD);
1050
1051         dNdEtaCorrectionND->GetTriggerBiasCorrectionNSD()->Scale(scaleND);
1052         dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scaleDD);
1053         dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scaleSD);
1054
1055         dNdEtaCorrectionND->GetTriggerBiasCorrectionND()->Scale(scaleND);
1056         dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scaleDD);
1057         dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scaleSD);
1058
1059         dNdEtaCorrectionND->GetTriggerBiasCorrectionOnePart()->Scale(scaleND);
1060         dNdEtaCorrectionDD->GetTriggerBiasCorrectionOnePart()->Scale(scaleDD);
1061         dNdEtaCorrectionSD->GetTriggerBiasCorrectionOnePart()->Scale(scaleSD);
1062       }
1063
1064       //clear track in correction
1065       dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
1066       dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
1067       dNdEtaCorrectionSD->GetTrack2ParticleCorrection()->Reset();
1068
1069       TList collection;
1070       collection.Add(correctionStandard);
1071       collection.Add(dNdEtaCorrectionND);
1072       collection.Add(dNdEtaCorrectionDD);
1073       collection.Add(dNdEtaCorrectionSD);
1074
1075       current->Merge(&collection);
1076       current->Finish();
1077       
1078       // print 0 bin efficiency
1079       TH1* hist2 = current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
1080       if (hist2->GetBinContent(1))
1081       {
1082         Float_t triggerEff = 100.0 / hist2->GetBinContent(1);
1083         Printf("trigger eff in 0 bin is: %.2f %%", triggerEff);
1084       }
1085
1086       corrections[counter] = current;
1087
1088       // now correct dNdeta distribution with modified correction map
1089       TFile::Open(analysisFileName);
1090
1091       dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
1092       fdNdEtaAnalysis->LoadHistograms();
1093
1094       fdNdEtaAnalysis->Finish(current, 0.151, correctionTarget, Form("%d %d", j, i), backgroundEvents);
1095
1096       name = "ratio";
1097       if (j==0) name.Append("_vetexReco_");
1098       if (j==1) name.Append("_triggerBias_");
1099       if (j==2) name.Append("_vertexReco_triggerBias_");
1100       name.Append(changes[i]);
1101
1102       hRatios[counter] = (TH1F*) fdNdEtaAnalysis->GetdNdEtaHistogram()->Clone(name);
1103
1104       name.Form("ND #times %0.2f DD #times %0.2f, SD #times %0.2f",scaleND,scaleDD,scaleSD);
1105       hRatios[counter]->SetTitle(name.Data());
1106       hRatios[counter]->SetYTitle("dN_{ch}/d#eta ratio #frac{default cross-section}{modified cross-sections}");
1107       
1108       TF1* pol0 = new TF1("pol0", "[0]", -0.49, 0.49);
1109       pol0->SetParameter(0, 0);
1110       hRatios[counter]->Fit(pol0, "RN");
1111       Printf("Case %d: %f", i, pol0->GetParameter(0));
1112       
1113       if (counter > 0)
1114         hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);
1115
1116       delete fdNdEtaAnalysis;
1117
1118       counter++;
1119     }
1120   }
1121
1122   TFile* fout = new TFile(outputFileName,"RECREATE");
1123
1124   // to make everything consistent
1125   hRatios[0]->Divide(hRatios[0],hRatios[0],1,1);
1126
1127   for (Int_t i=0; i<counter; i++)
1128   {
1129     corrections[i]->SaveHistograms();
1130     hRatios[i]->Write();
1131   }
1132
1133   fout->Write();
1134   fout->Close();
1135 }
1136
1137 void GetRelativeFractions(Int_t origin, Float_t& ref_SD, Float_t& ref_DD, Float_t& ref_ND, Float_t& error_SD, Float_t& error_DD, Float_t& error_ND)
1138 {
1139   // origin: 
1140   //   -1 = Pythia (test)
1141   //   0 = UA5
1142   //   1 = Data 1.8 TeV
1143   //   2 = Tel-Aviv
1144   //   3 = Durham
1145   //
1146
1147   switch (origin)
1148   {
1149     case -10: // Pythia default at 7 GeV, 50% error
1150       Printf("PYTHIA x-sections");
1151       ref_SD = 0.192637; error_SD = ref_SD * 0.5;
1152       ref_DD = 0.129877; error_DD = ref_DD * 0.5;
1153       ref_ND = 0.677486; error_ND = 0;
1154       break;
1155
1156     case -1: // Pythia default at 900 GeV, as test
1157       Printf("PYTHIA x-sections");
1158       ref_SD = 0.223788;
1159       ref_DD = 0.123315;
1160       ref_ND = 0.652897;
1161       break;
1162       
1163     case 0: // UA5
1164       Printf("UA5 x-sections a la first paper");
1165       ref_SD = 0.153; error_SD = 0.05;
1166       ref_DD = 0.080; error_DD = 0.05;
1167       ref_ND = 0.767; error_ND = 0;
1168       break;
1169       
1170     case 10: // UA5
1171       Printf("UA5 x-sections hadron level definition for Pythia"); 
1172       // Fractions in Pythia with UA5 cuts selection for SD
1173       // ND: 0.688662
1174       // SD: 0.188588 --> this should be 15.3
1175       // DD: 0.122750
1176       ref_SD = 0.224 * 0.153 / 0.189; error_SD = 0.023 * 0.224 / 0.189;
1177       ref_DD = 0.095;                 error_DD = 0.06; 
1178       ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
1179       break;
1180     
1181     case 11: // UA5
1182       Printf("UA5 x-sections hadron level definition for Phojet"); 
1183       // Fractions in Phojet with UA5 cuts selection for SD
1184       // ND: 0.783573
1185       // SD: 0.151601 --> this should be 15.3
1186       // DD: 0.064827
1187       ref_SD = 0.191 * 0.153 / 0.152; error_SD = 0.023 * 0.191 / 0.152;
1188       ref_DD = 0.095;                 error_DD = 0.06; 
1189       ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
1190       break;
1191       
1192     case 20: // E710, 1.8 TeV
1193       Printf("E710 x-sections hadron level definition for Pythia");
1194       // ND: 0.705709
1195       // SD: 0.166590 --> this should be 15.9
1196       // DD: 0.127701
1197       ref_SD = 0.217 * 0.159 / 0.167; error_SD = 0.024 * 0.217 / 0.167;
1198       ref_DD = 0.075 * 1.43;          error_DD = 0.02 * 1.43; 
1199       ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
1200       break;
1201     
1202     case 21: // E710, 1.8 TeV
1203       Printf("E710 x-sections hadron level definition for Phojet"); 
1204       // ND: 0.817462
1205       // SD: 0.125506 --> this should be 15.9
1206       // DD: 0.057032
1207       ref_SD = 0.161 * 0.159 / 0.126; error_SD = 0.024 * 0.161 / 0.126;
1208       ref_DD = 0.075 * 1.43;         error_DD = 0.02 * 1.43;
1209       ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
1210       break;
1211     
1212     case 1: // data 1.8 TeV
1213       Printf("??? x-sections");
1214       ref_SD = 0.152;
1215       ref_DD = 0.092;
1216       ref_ND = 1 - ref_SD - ref_DD;
1217       break;
1218       
1219     case 2: // tel-aviv model
1220       Printf("Tel-aviv model x-sections");
1221       ref_SD = 0.171;
1222       ref_DD = 0.094;
1223       ref_ND = 1 - ref_SD - ref_DD;
1224       break;
1225     
1226     case 3: // durham model
1227       Printf("Durham model x-sections");
1228       ref_SD = 0.190;
1229       ref_DD = 0.125;
1230       ref_ND = 1 - ref_SD - ref_DD;
1231       break;
1232     
1233     default:
1234       AliFatal(Form("Unknown origin %d", origin));
1235   }
1236 }
1237
1238 void CreateCorrectionsWithUA5CrossSections(Int_t origin, const Char_t* correctionFileName="correction_mapprocess-types.root", const Char_t* outputFileName="correction_map2.root") {
1239   //
1240   // Function used to merge standard corrections with vertex
1241   // reconstruction corrections obtained by a certain mix of ND, DD
1242   // and SD events.
1243   //
1244   loadlibs();
1245
1246   const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
1247   
1248   Float_t ref_SD = -1;
1249   Float_t ref_DD = -1;
1250   Float_t ref_ND = -1;
1251   
1252   Float_t error_SD = -1;
1253   Float_t error_DD = -1;
1254   Float_t error_ND = -1;
1255   
1256   GetRelativeFractions(origin, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
1257   
1258   Printf("Using x-sections:\n SD: %f +- %f\n DD: %f +- %f\n ND: %f +- %f", ref_SD, error_SD, ref_DD, error_DD, ref_ND, error_ND);
1259   
1260 //Karel (UA5):
1261 //     fsd = 0.153 +- 0.031
1262 //     fdd = 0.080 +- 0.050
1263 //     fnd = 0.767 +- 0.059
1264
1265 //       Karel (1.8 TeV):
1266 //       
1267 //       Tel-Aviv model Sd/Inel = 0.171           Dd/Inel = 0.094
1268 //       Durham model   Sd/Inel = 0.190           Dd/Inel = 0.125
1269 //       Data           Sd/Inel = 0.152 +- 0.030  Dd/Inel = 0.092 +- 0.45
1270
1271   // standard correction
1272   TFile::Open(correctionFileName);
1273   AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1274   correctionStandard->LoadHistograms();
1275
1276   // dont take vertexreco from this one
1277   correctionStandard->GetVertexRecoCorrection()->Reset();
1278   // dont take triggerbias from this one
1279   correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
1280   correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
1281   correctionStandard->GetTriggerBiasCorrectionND()->Reset();
1282   correctionStandard->GetTriggerBiasCorrectionOnePart()->Reset();
1283
1284   AlidNdEtaCorrection* corrections[100];
1285   TH1F* hRatios[100];
1286
1287   Int_t counter = 0;
1288       
1289   TFile::Open(correctionFileName);
1290
1291   AlidNdEtaCorrection* current = new AlidNdEtaCorrection("dndeta_correction_ua5", "dndeta_correction_ua5");
1292   current->LoadHistograms("dndeta_correction");
1293   current->Reset();
1294
1295   TString name;
1296   name.Form("dndeta_correction_ND");
1297   AlidNdEtaCorrection* dNdEtaCorrectionND = new AlidNdEtaCorrection(name,name);
1298   dNdEtaCorrectionND->LoadHistograms();
1299   name.Form("dndeta_correction_DD");
1300   AlidNdEtaCorrection* dNdEtaCorrectionDD = new AlidNdEtaCorrection(name,name);
1301   dNdEtaCorrectionDD->LoadHistograms();
1302   name.Form("dndeta_correction_SD");
1303   AlidNdEtaCorrection* dNdEtaCorrectionSD = new AlidNdEtaCorrection(name,name);
1304   dNdEtaCorrectionSD->LoadHistograms();
1305
1306   // calculating relative
1307   Float_t nd = dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1308   Float_t dd = dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1309   Float_t sd = dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1310   Float_t total = nd + dd + sd;
1311   
1312   nd /= total;
1313   sd /= total;
1314   dd /= total;
1315   
1316   Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
1317   
1318   Float_t scaleND = ref_ND / nd;
1319   Float_t scaleDD = ref_DD / dd;
1320   Float_t scaleSD = ref_SD / sd;
1321   
1322   Printf("ND=%.2f, DD=%.2f, SD=%.2f",scaleND, scaleDD, scaleSD);
1323
1324   // scale
1325   dNdEtaCorrectionND->GetVertexRecoCorrection()->Scale(scaleND);
1326   dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scaleDD);
1327   dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scaleSD);
1328     
1329   dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->Scale(scaleND);
1330   dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scaleDD);
1331   dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scaleSD);
1332
1333   dNdEtaCorrectionND->GetTriggerBiasCorrectionNSD()->Scale(scaleND);
1334   dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scaleDD);
1335   dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scaleSD);
1336
1337   dNdEtaCorrectionND->GetTriggerBiasCorrectionND()->Scale(scaleND);
1338   dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scaleDD);
1339   dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scaleSD);
1340
1341   dNdEtaCorrectionND->GetTriggerBiasCorrectionOnePart()->Scale(scaleND);
1342   dNdEtaCorrectionDD->GetTriggerBiasCorrectionOnePart()->Scale(scaleDD);
1343   dNdEtaCorrectionSD->GetTriggerBiasCorrectionOnePart()->Scale(scaleSD);
1344
1345   //clear track in correction
1346   dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
1347   dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
1348   dNdEtaCorrectionSD->GetTrack2ParticleCorrection()->Reset();
1349
1350   TList collection;
1351   collection.Add(correctionStandard);
1352   collection.Add(dNdEtaCorrectionND);
1353   collection.Add(dNdEtaCorrectionDD);
1354   collection.Add(dNdEtaCorrectionSD);
1355
1356   current->Merge(&collection);
1357   current->Finish();
1358
1359   // print 0 bin efficiency
1360   TH1* hist2 = current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
1361   if (hist2->GetBinContent(1) > 0)
1362   {
1363     Float_t triggerEff = 100.0 / hist2->GetBinContent(1);
1364     Printf("trigger eff in 0 bin is: %.2f %%", triggerEff);
1365   }
1366   
1367   TFile* fout = new TFile(outputFileName,"RECREATE");
1368   current->SaveHistograms();
1369
1370   fout->Write();
1371   fout->Close();
1372
1373   Printf("Trigger efficiencies:");
1374   Printf("ND: %.2f %%", 100.0 * dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1375   Printf("SD: %.2f %%", 100.0 * dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1376   Printf("DD: %.2f %%", 100.0 * dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1377   Printf("INEL: %.2f %%", 100.0 * current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1378   Printf("NSD: %.2f %%", 100.0 * (dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() + dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral()) / (dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral() + dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral()));
1379 }
1380
1381 DrawTriggerEfficiency(Char_t* fileName) {
1382
1383   gStyle->SetOptStat(0);
1384   gStyle->SetOptTitle(0);
1385   gStyle->SetOptFit(0);
1386
1387   gStyle->SetTextSize(0.04);
1388   gStyle->SetTitleSize(0.05,"xyz");
1389   //gStyle->SetTitleFont(133, "xyz");
1390   //gStyle->SetLabelFont(133, "xyz");
1391   //gStyle->SetLabelSize(17, "xyz");
1392   gStyle->SetLabelOffset(0.01, "xyz");
1393
1394   gStyle->SetTitleOffset(1.1, "y");
1395   gStyle->SetTitleOffset(1.1, "x");
1396   gStyle->SetEndErrorSize(0.0);
1397
1398   //##############################################
1399
1400   //making canvas and pads
1401   TCanvas *c = new TCanvas("trigger_eff", "",600,500);
1402
1403   TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
1404
1405   p1->SetBottomMargin(0.15);
1406   p1->SetTopMargin(0.03);
1407   p1->SetLeftMargin(0.15);
1408   p1->SetRightMargin(0.03);
1409   
1410   p1->SetGridx();
1411   p1->SetGridy();
1412
1413   p1->Draw();
1414   p1->cd();
1415
1416   gSystem->Load("libPWG0base");
1417
1418   AlidNdEtaCorrection* corrections[4];
1419   AliCorrectionMatrix2D* triggerBiasCorrection[4];
1420
1421   TH1F* hTriggerEffInv[4];
1422   TH1F* hTriggerEff[4];
1423
1424   Char_t* names[] = {"triggerBiasND", "triggerBiasDD", "triggerBiasSD", "triggerBiasINEL"};
1425   
1426   for (Int_t i=0; i<4; i++) {
1427     corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1428     corrections[i]->LoadHistograms(fileName, names[i]);    
1429
1430     triggerBiasCorrection[i] = corrections[i]->GetTriggerBiasCorrectionINEL();
1431
1432     
1433     hTriggerEffInv[i] = triggerBiasCorrection[i]->Get1DCorrection();
1434     hTriggerEff[i]    = (TH1F*)hTriggerEffInv[i]->Clone();
1435     
1436     for (Int_t b=0; b<=hTriggerEff[i]->GetNbinsX(); b++) {
1437       hTriggerEff[i]->SetBinContent(b,1);
1438       hTriggerEff[i]->SetBinError(b,0);
1439     }
1440     hTriggerEff[i]->Divide(hTriggerEff[i],hTriggerEffInv[i]);
1441     hTriggerEff[i]->Scale(100);
1442   }
1443
1444   Int_t colors[] = {2,3,4,1};
1445   Float_t effs[4];
1446   for (Int_t i=0; i<4; i++) {
1447     hTriggerEff[i]->Fit("pol0","","",-20,20);
1448     effs[i] = ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->GetParameter(0);
1449     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineWidth(1);
1450     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineStyle(2);
1451     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineColor(colors[i]);
1452     cout << effs[i] << endl;
1453   }
1454
1455
1456   Char_t* text[] = {"ND", "DD", "SD", "INEL"};
1457   TLatex* latex[4];
1458
1459   TH2F* null = new TH2F("","",100,-25,35,100,0,110);
1460   null->SetXTitle("Vertex z [cm]");
1461   null->GetXaxis()->CenterTitle(kTRUE);
1462   null->SetYTitle("Trigger efficiency [%]");
1463   null->GetYaxis()->CenterTitle(kTRUE);
1464   null->Draw();
1465
1466
1467   for (Int_t i=0; i<4; i++) {
1468     hTriggerEff[i]->SetLineWidth(2);
1469     hTriggerEff[i]->SetLineColor(colors[i]);
1470
1471     hTriggerEff[i]->Draw("same");
1472
1473     latex[i] = new TLatex(22,effs[i]-1.5, Form("%s (%0.1f)",text[i],effs[i]));
1474     latex[i]->SetTextColor(colors[i]);
1475     latex[i]->Draw();
1476   }
1477   
1478 }
1479
1480
1481 DrawSpectraPID(Char_t* fileName) {
1482
1483   gSystem->Load("libPWG0base");
1484
1485   Char_t* names[] = {"correction_0", "correction_1", "correction_2", "correction_3"};
1486   AlidNdEtaCorrection* corrections[4];
1487   AliCorrectionMatrix3D* trackToPartCorrection[4];
1488
1489   TH1F* measuredPt[4];
1490   TH1F* generatedPt[4];
1491   TH1F* ratioPt[4];
1492
1493   for (Int_t i=0; i<4; i++) {
1494     corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1495     corrections[i]->LoadHistograms(fileName, names[i]);    
1496
1497     trackToPartCorrection[i] = corrections[i]->GetTrack2ParticleCorrection();
1498
1499     Int_t binX1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(-10);
1500     Int_t binX2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(10);
1501     Int_t binY1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(-1);
1502     Int_t binY2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(1);
1503
1504     measuredPt[i]  = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->ProjectionZ(Form("m_%d",i),binX1,binX2,binY1,binY2);
1505     generatedPt[i] = (TH1F*)trackToPartCorrection[i]->GetGeneratedHistogram()->ProjectionZ(Form("g_%d",i),binX1,binX2,binY1,binY2);
1506     ratioPt[i] = (TH1F*)generatedPt[i]->Clone(Form("r_%d",i));
1507     ratioPt[i]->Divide(measuredPt[i], generatedPt[i], 1,1,"B");
1508   }
1509   
1510   ratioPt[0]->Draw();
1511   
1512   for (Int_t i=0; i<3; i++) {
1513     ratioPt[i]->SetLineColor(i+1);
1514     ratioPt[i]->SetLineWidth(2);
1515     
1516     ratioPt[i]->Draw("same");
1517     
1518   }
1519
1520   return;
1521   measuredPt[0]->SetLineColor(2);
1522   measuredPt[0]->SetLineWidth(5);
1523
1524   measuredPt[0]->Draw();
1525   generatedPt[0]->Draw("same");
1526 }
1527
1528 void changePtSpectrum(const char* fileName = "analysis_mc.root", Float_t ptCutOff = 0.2, const char* fileName2 = 0)
1529 {
1530   Float_t factor = 0.5;
1531
1532   TFile* file = TFile::Open(fileName);
1533   TH1F* hist = dynamic_cast<TH1F*> (file->Get("dndeta_check_pt")->Clone());
1534   
1535   TH1* hist2 = 0;
1536   if (fileName2)
1537   {
1538     file2 = TFile::Open(fileName2);
1539     hist2 = dynamic_cast<TH1*> (file2->Get("dndeta_check_pt")->Clone());
1540     hist2->Scale(hist->GetBinContent(hist->FindBin(ptCutOff)) / hist2->GetBinContent(hist2->FindBin(ptCutOff)));
1541   }
1542   
1543   //hist->Scale(1.0 / hist->Integral());
1544
1545   //hist->Rebin(3);
1546   //hist->Scale(1.0/3);
1547
1548   TH1F* clone1 = dynamic_cast<TH1F*> (hist->Clone("clone1"));
1549   TH1F* clone2 = dynamic_cast<TH1F*> (hist->Clone("clone2"));
1550
1551   TH1F* scale1 =  dynamic_cast<TH1F*> (hist->Clone("scale1"));
1552   TH1F* scale2 =  dynamic_cast<TH1F*> (hist->Clone("scale2"));
1553
1554   for (Int_t i=1; i <= hist->GetNbinsX(); ++i)
1555   {
1556     if (hist->GetBinCenter(i) > ptCutOff)
1557     {
1558       scale1->SetBinContent(i, 1);
1559       scale2->SetBinContent(i, 1);
1560     }
1561     else
1562     {
1563       // 90 % at pt = 0, 0% at pt = ptcutoff
1564       scale1->SetBinContent(i, 1 - (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * factor);
1565
1566       // 110% at pt = 0, ...
1567       scale2->SetBinContent(i, 1 + (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * factor);
1568     }
1569     scale1->SetBinError(i, 0);
1570     scale2->SetBinError(i, 0);
1571   }
1572
1573   /*
1574   new TCanvas;
1575   scale1->Draw();
1576   scale2->SetLineColor(kRed);
1577   scale2->Draw("SAME");
1578   */
1579
1580   clone1->Multiply(scale1);
1581   clone2->Multiply(scale2);
1582
1583   Prepare1DPlot(hist);
1584   Prepare1DPlot(clone1);
1585   Prepare1DPlot(clone2);
1586
1587   /*hist->SetMarkerStyle(markers[0]);
1588   clone1->SetMarkerStyle(markers[0]);
1589   clone2->SetMarkerStyle(markers[0]);*/
1590
1591   hist->SetTitle(";p_{T} in GeV/c;dN_{ch}/dp_{T} in c/GeV");
1592   hist->GetXaxis()->SetRangeUser(0, 0.5);
1593   hist->GetYaxis()->SetRangeUser(0.01, clone2->GetMaximum() * 1.1);
1594   hist->GetYaxis()->SetTitleOffset(1);
1595
1596   TCanvas* canvas = new TCanvas("c", "c", 600, 600);
1597   gPad->SetGridx();
1598   gPad->SetGridy();
1599   gPad->SetTopMargin(0.05);
1600   gPad->SetRightMargin(0.05);
1601   gPad->SetBottomMargin(0.12);
1602   hist->Draw("H");
1603   clone1->SetLineColor(kRed);
1604   clone1->Draw("HSAME");
1605   clone2->SetLineColor(kBlue);
1606   clone2->Draw("HSAME");
1607   hist->Draw("HSAME");
1608   
1609   if (hist2)
1610   {
1611     Prepare1DPlot(hist2);
1612     hist2->SetLineStyle(2);
1613     hist2->Draw("HSAME");
1614   }
1615
1616   Float_t fraction =  hist->Integral(hist->GetXaxis()->FindBin(ptCutOff), hist->GetNbinsX()) / hist->Integral(1, hist->GetNbinsX());
1617   Float_t fraction1 = clone1->Integral(clone1->GetXaxis()->FindBin(ptCutOff), clone1->GetNbinsX()) / clone1->Integral(1, clone1->GetNbinsX());
1618   Float_t fraction2 = clone2->Integral(clone2->GetXaxis()->FindBin(ptCutOff), clone2->GetNbinsX()) / clone2->Integral(1, clone2->GetNbinsX());
1619
1620   printf("%f %f %f\n", fraction, fraction1, fraction2);
1621   printf("Rel. %f %f\n", fraction1 / fraction, fraction2 / fraction);
1622
1623   //canvas->SaveAs("changePtSpectrum.gif");
1624   canvas->SaveAs("changePtSpectrum.eps");
1625 }
1626
1627 void FractionBelowPt()
1628 {
1629   gSystem->Load("libPWG0base");
1630
1631   AlidNdEtaCorrection* fdNdEtaCorrection[4];
1632
1633   TFile::Open("systematics.root");
1634
1635   for (Int_t i=0; i<4; ++i)
1636   {
1637     TString name;
1638     name.Form("correction_%d", i);
1639     fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
1640     fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
1641   }
1642
1643   Double_t geneCount[5];
1644   Double_t measCount[5];
1645   geneCount[4] = 0;
1646   measCount[4] = 0;
1647
1648   for (Int_t i=0; i<4; ++i)
1649   {
1650     TH3F* hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1651     geneCount[i] = hist->Integral(hist->GetXaxis()->FindBin(-10), hist->GetXaxis()->FindBin(10),
1652                                   hist->GetYaxis()->FindBin(-0.8), hist->GetYaxis()->FindBin(0.8),
1653                                   1, hist->GetZaxis()->FindBin(0.3));
1654
1655     hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1656     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));
1657
1658     geneCount[4] += geneCount[i];
1659     measCount[4] += measCount[i];
1660
1661     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]);
1662   }
1663
1664   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]);
1665
1666   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]);
1667 }
1668
1669
1670 mergeCorrectionsMisalignment(Char_t* alignedFile = "correction_map_aligned.root",
1671                                              Char_t* misalignedFile = "correction_map_misaligned.root",
1672                                              Char_t* outputFileName="correction_map_misaligned_single.root")
1673 {
1674   //
1675   // from the aligned and misaligned corrections, 3 new corrections are created
1676   // in these new corrections only one of the corrections (track2particle, vertex, trigger)
1677   // is taken from the misaligned input to allow study of the effect on the different
1678   // corrections
1679
1680   gSystem->Load("libPWG0base");
1681
1682   const Char_t* typeName[] = { "track2particle", "vertex", "trigger" };
1683
1684   AlidNdEtaCorrection* corrections[3];
1685   for (Int_t j=0; j<3; j++) { // j = 0 (track2particle), j = 1 (vertex), j = 2 (trigger)
1686     AlidNdEtaCorrection* alignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1687     alignedCorrection->LoadHistograms(alignedFile);
1688
1689     AlidNdEtaCorrection* misalignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1690     misalignedCorrection->LoadHistograms(misalignedFile);
1691
1692     TString name;
1693     name.Form("dndeta_correction_alignment_%s", typeName[j]);
1694     AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
1695
1696     switch (j)
1697     {
1698       case 0:
1699         alignedCorrection->GetTrack2ParticleCorrection()->Reset();
1700         misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1701         misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1702         misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1703         misalignedCorrection->GetVertexRecoCorrection()->Reset();
1704         break;
1705
1706       case 1:
1707         misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1708         misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1709         misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1710         misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1711         alignedCorrection->GetVertexRecoCorrection()->Reset();
1712         break;
1713
1714       case 2:
1715         misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1716         alignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1717         alignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1718         alignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1719         misalignedCorrection->GetVertexRecoCorrection()->Reset();
1720         break;
1721
1722       default:
1723         return;
1724     }
1725
1726     TList collection;
1727     collection.Add(misalignedCorrection);
1728     collection.Add(alignedCorrection);
1729
1730     current->Merge(&collection);
1731     current->Finish();
1732
1733     corrections[j] = current;
1734   }
1735
1736   TFile* fout = new TFile(outputFileName, "RECREATE");
1737
1738   for (Int_t i=0; i<3; i++)
1739     corrections[i]->SaveHistograms();
1740
1741   fout->Write();
1742   fout->Close();
1743 }
1744
1745
1746 void DrawdNdEtaDifferencesAlignment()
1747 {
1748   TH1* hists[5];
1749
1750   TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
1751   legend->SetFillColor(0);
1752
1753   TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
1754   canvas->Divide(2, 1);
1755
1756   canvas->cd(1);
1757
1758   for (Int_t i=0; i<5; ++i)
1759   {
1760     hists[i] = 0;
1761     TFile* file = 0;
1762     TString title;
1763
1764     switch(i)
1765     {
1766       case 0 : file = TFile::Open("systematics_misalignment_aligned.root"); title = "aligned"; break;
1767       case 1 : file = TFile::Open("systematics_misalignment_misaligned.root"); title = "fully misaligned"; break;
1768       case 2 : file = TFile::Open("systematics_misalignment_track2particle.root"); title = "only track2particle"; break;
1769       case 3 : file = TFile::Open("systematics_misalignment_vertex.root"); title = "only vertex rec."; break;
1770       case 4 : file = TFile::Open("systematics_misalignment_trigger.root"); title = "only trigger bias"; break;
1771       default: return;
1772     }
1773
1774     if (file)
1775     {
1776       hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
1777       hists[i]->SetTitle("");
1778
1779       Prepare1DPlot(hists[i], kFALSE);
1780       hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
1781       hists[i]->GetYaxis()->SetRangeUser(6, 7);
1782       hists[i]->SetLineWidth(1);
1783       hists[i]->SetLineColor(colors[i]);
1784       hists[i]->SetMarkerColor(colors[i]);
1785       hists[i]->SetMarkerStyle(markers[i]);
1786       hists[i]->GetXaxis()->SetLabelOffset(0.015);
1787       hists[i]->GetYaxis()->SetTitleOffset(1.5);
1788       gPad->SetLeftMargin(0.12);
1789       hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
1790
1791       legend->AddEntry(hists[i], title);
1792       hists[i]->SetTitle(title);
1793     }
1794   }
1795   legend->Draw();
1796
1797   canvas->cd(2);
1798   gPad->SetLeftMargin(0.14);
1799
1800   TLegend* legend2 = new TLegend(0.63, 0.73, 0.98, 0.98);
1801   legend2->SetFillColor(0);
1802
1803   for (Int_t i=1; i<5; ++i)
1804   {
1805     if (hists[i])
1806     {
1807       legend2->AddEntry(hists[i]);
1808
1809       hists[i]->Divide(hists[0]);
1810       hists[i]->SetTitle("b)");
1811       hists[i]->GetYaxis()->SetRangeUser(0.9, 1.1);
1812       hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
1813       hists[i]->GetYaxis()->SetTitleOffset(1.8);
1814       hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
1815     }
1816   }
1817
1818   legend2->Draw();
1819
1820   canvas->SaveAs("misalignment_result.eps");
1821   canvas->SaveAs("misalignment_result.gif");
1822 }
1823
1824 void EvaluateMultiplicityEffect()
1825 {
1826   gSystem->Load("libPWG0base");
1827
1828   AlidNdEtaCorrection* fdNdEtaCorrectionLow[4];
1829   AlidNdEtaCorrection* fdNdEtaCorrectionHigh[4];
1830
1831   TFile::Open("systematics-low-multiplicity.root");
1832
1833   for (Int_t i=0; i<4; ++i)
1834   {
1835     TString name;
1836     name.Form("correction_%d", i);
1837     fdNdEtaCorrectionLow[i] = new AlidNdEtaCorrection(name, name);
1838     fdNdEtaCorrectionLow[i]->LoadHistograms("systematics-low-multiplicity.root", name);
1839   }
1840
1841   TList list;
1842   for (Int_t i=1; i<4; ++i)
1843     list.Add(fdNdEtaCorrectionLow[i]);
1844
1845   fdNdEtaCorrectionLow[0]->Merge(&list);
1846   fdNdEtaCorrectionLow[0]->Finish();
1847
1848   TFile::Open("systematics-high-multiplicity.root");
1849
1850   for (Int_t i=0; i<4; ++i)
1851   {
1852     TString name;
1853     name.Form("correction_%d", i);
1854     fdNdEtaCorrectionHigh[i] = new AlidNdEtaCorrection(name, name);
1855     fdNdEtaCorrectionHigh[i]->LoadHistograms("systematics-high-multiplicity.root", name);
1856   }
1857
1858   TList list2;
1859   for (Int_t i=1; i<4; ++i)
1860     list2.Add(fdNdEtaCorrectionHigh[i]);
1861
1862   fdNdEtaCorrectionHigh[0]->Merge(&list2);
1863   fdNdEtaCorrectionHigh[0]->Finish();
1864
1865   TH1F* outputLow = new TH1F("Track2ParticleLow", "Track2Particle at low multiplicity", 200, 0, 2);
1866   TH1F* outputHigh = new TH1F("Track2ParticleHigh", "Track2Particle at high multiplicity", 200, 0, 2);
1867
1868   TH3F* hist = fdNdEtaCorrectionLow[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1869   TH3F* hist2 = fdNdEtaCorrectionHigh[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1870   for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
1871     for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
1872       for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
1873       //for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
1874       {
1875         if (hist->GetBinContent(x, y, z) > 0)
1876           outputLow->Fill(hist->GetBinContent(x, y, z));
1877         //if (hist->GetBinContent(x, y, z) == 1)
1878         //  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));
1879
1880         if (hist2->GetBinContent(x, y, z) > 0)
1881           outputHigh->Fill(hist2->GetBinContent(x, y, z));
1882       }
1883
1884   TCanvas* canvas = new TCanvas("EvaluateMultiplicityEffect", "EvaluateMultiplicityEffect", 1000, 500);
1885   canvas->Divide(2, 1);
1886
1887   canvas->cd(1);
1888   outputLow->Draw();
1889   outputLow->Fit("gaus", "0");
1890   outputLow->GetFunction("gaus")->SetLineColor(2);
1891   outputLow->GetFunction("gaus")->DrawCopy("SAME");
1892
1893   canvas->cd(2);
1894   outputHigh->Draw();
1895   outputHigh->Fit("gaus", "0");
1896   outputHigh->GetFunction("gaus")->DrawCopy("SAME");
1897
1898   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1899 }
1900
1901 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")
1902 {
1903     gSystem->Load("libPWG0base");
1904         
1905         AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1906     dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
1907     
1908     dNdEtaCorrection->GetTrack2ParticleCorrection()->PlotBinErrors(xmin, xmax, ymin, ymax, zmin, zmax)->Draw();
1909 }
1910
1911
1912 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")
1913 {
1914   gSystem->Load("libPWG0base");
1915
1916   TFile* file = TFile::Open(output, "RECREATE");
1917
1918   const Int_t max = 5;
1919   dNdEtaAnalysis* fdNdEtaAnalysis[5];
1920
1921   new TCanvas;
1922   TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
1923   legend->SetFillColor(0);
1924
1925   for (Int_t i = 0; i < max; ++i)
1926   {
1927     TFile::Open(baseCorrectionMapFile);
1928     AlidNdEtaCorrection* baseCorrection = new AlidNdEtaCorrection(baseCorrectionMapFolder, baseCorrectionMapFolder);
1929     baseCorrection->LoadHistograms();
1930
1931     AlidNdEtaCorrection::CorrectionType correctionType = AlidNdEtaCorrection::kNone;
1932     const char* name = 0;
1933
1934     TFile::Open(changedCorrectionMapFile);
1935     switch (i)
1936     {
1937       case 0 :
1938         name = "default";
1939         break;
1940
1941       case 1 :
1942         baseCorrection->GetTrack2ParticleCorrection()->LoadHistograms(Form("%s/Track2Particle", changedCorrectionMapFolder));
1943         name = "Track2Particle";
1944         break;
1945
1946       case 2 :
1947         baseCorrection->GetVertexRecoCorrection()->LoadHistograms(Form("%s/VertexReconstruction", changedCorrectionMapFolder));
1948         name = "VertexReco";
1949         break;
1950
1951       case 3 :
1952         baseCorrection->GetTriggerBiasCorrectionINEL()->LoadHistograms(Form("%s/TriggerBias_MBToINEL", changedCorrectionMapFolder));
1953         name = "TriggerBias_MBToINEL";
1954         break;
1955
1956       case 4 :
1957         baseCorrection->LoadHistograms(changedCorrectionMapFolder);
1958         name = "all";
1959         break;
1960
1961       default: return;
1962     }
1963
1964     TFile::Open(dataFile);
1965     fdNdEtaAnalysis[i] = new dNdEtaAnalysis(name, name);
1966     fdNdEtaAnalysis[i]->LoadHistograms("dndeta");
1967
1968     fdNdEtaAnalysis[i]->Finish(baseCorrection, 0.3, AlidNdEtaCorrection::kINEL);
1969     file->cd();
1970     fdNdEtaAnalysis[i]->SaveHistograms();
1971
1972     TH1* hist = fdNdEtaAnalysis[i]->GetdNdEtaHistogram(0);
1973     hist->SetLineColor(colors[i]);
1974     hist->SetMarkerColor(colors[i]);
1975     hist->SetMarkerStyle(markers[i]+1);
1976     hist->DrawCopy((i == 0) ? "" : "SAME");
1977     legend->AddEntry(hist, name);
1978   }
1979
1980   legend->Draw();
1981 }
1982
1983 void ChangePtInCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
1984 {
1985   loadlibs();
1986   if (!TFile::Open(fileName))
1987     return;
1988
1989   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
1990   if (!dNdEtaCorrection->LoadHistograms())
1991     return;
1992
1993   dNdEtaCorrection->Finish();
1994
1995   AliCorrection* corr = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kTrack2Particle);
1996  
1997   Printf(">>>> Before");
1998   corr->PrintInfo(0);
1999
2000   Float_t factor = 0.5;
2001   Float_t ptCutOff = 0.2;
2002   
2003   TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram();
2004   TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram();
2005   
2006   for (Int_t z = 1; z <= gene->GetZaxis()->FindBin(ptCutOff - 0.01); z++)
2007   {
2008     Float_t localFactor = 1 - (ptCutOff - gene->GetZaxis()->GetBinCenter(z)) / ptCutOff * factor;
2009     Printf("%f %f", gene->GetZaxis()->GetBinCenter(z), localFactor);
2010     for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
2011       for (Int_t y = 1; y <= gene->GetNbinsY(); y++)
2012       {
2013         gene->SetBinContent(x, y, z, gene->GetBinContent(x, y, z) * localFactor);
2014         meas->SetBinContent(x, y, z, meas->GetBinContent(x, y, z) * localFactor);
2015       }
2016   }
2017   
2018   dNdEtaCorrection->Finish();
2019   
2020   Printf(">>>> After");
2021   corr->PrintInfo(0);
2022 }
2023
2024 Float_t FitAverage(TH1* hist, Int_t n, Float_t* begin, Float_t *end, Int_t color, Int_t& totalBins)
2025 {
2026   Float_t average = 0;
2027   totalBins = 0;
2028   
2029   for (Int_t i=0; i<n; i++)
2030   {
2031     func = new TF1("func", "[0]", hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->FindBin(begin[i])), hist->GetXaxis()->GetBinUpEdge(hist->GetXaxis()->FindBin(end[i])));
2032     Int_t bins = hist->GetXaxis()->FindBin(end[i]) - hist->GetXaxis()->FindBin(begin[i]) + 1;
2033     func->SetParameter(0, 1);
2034     func->SetLineColor(color);
2035
2036     hist->Fit(func, "RNQ");
2037     func->Draw("SAME");
2038     
2039     average += func->GetParameter(0) * bins;
2040     totalBins += bins;
2041   }
2042   
2043   return average / totalBins;
2044 }
2045
2046 void SPDIntegrateGaps(Bool_t all, const char* mcFile = "../../../LHC10b2/v0or/spd/analysis_esd_raw.root")
2047 {
2048   Float_t eta = 1.29;
2049   Int_t binBegin = ((TH2*) gFile->Get("fEtaPhi"))->GetXaxis()->FindBin(-eta);
2050   Int_t binEnd   = ((TH2*) gFile->Get("fEtaPhi"))->GetXaxis()->FindBin(eta);
2051   
2052   Printf("eta range: %f bins: %d %d", eta, binBegin, binEnd);
2053   
2054   if (!all)
2055     Printf("Eta smaller than 0 side");
2056   
2057   c = new TCanvas;
2058   TFile::Open("analysis_esd_raw.root");
2059   hist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("hist", binBegin, (all) ? binEnd : 40);
2060   hist->Rebin(2);
2061   hist->SetStats(0);
2062   hist->Sumw2();
2063   hist->Draw("HIST");
2064   gPad->SetGridx();
2065   gPad->SetGridy();
2066   
2067   TFile::Open(mcFile);  
2068   mcHist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("mcHist", binBegin, (all) ? binEnd : 40);
2069   mcHist->Rebin(2);
2070   mcHist->SetLineColor(2);
2071   mcHist->Scale(hist->Integral() / mcHist->Integral());
2072   mcHist->Draw("SAME");
2073   
2074   Float_t add = 0;
2075   Int_t bins;
2076   
2077   Float_t okRangeBegin[] = { 0.04, 0.67, 1.34 };
2078   Float_t okRangeEnd[] =   { 0.55, 1.24, 1.63 };
2079   Float_t gapRangeBegin[] = { 0.6, 1.27  };
2080   Float_t gapRangeEnd[] =   { 0.65, 1.32 };
2081   Float_t averageOK  = FitAverage(hist, 3, okRangeBegin, okRangeEnd, 1, bins);
2082   Float_t averageGap = FitAverage(hist, 2, gapRangeBegin, gapRangeEnd, 2, bins);
2083   add += bins * (averageOK - averageGap);
2084   Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
2085   
2086   Float_t okRangeBegin2[] = { 2.4, 2.65 };
2087   Float_t okRangeEnd2[] =   { 2.55, 3.2 };
2088   Float_t gapRangeBegin2[] = { 2.59, 3.3 };
2089   Float_t gapRangeEnd2[] =   { 2.61, 3.3 };
2090   averageOK  = FitAverage(hist, 2, okRangeBegin2, okRangeEnd2, 1, bins);
2091   averageGap = FitAverage(hist, 2, gapRangeBegin2, gapRangeEnd2, 2, bins);
2092   add += bins * (averageOK - averageGap);
2093   Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
2094   
2095   Float_t okRangeBegin3[] = { 3.55, 3.9 };
2096   Float_t okRangeEnd3[] =   { 3.8, 4.15 };
2097   Float_t gapRangeBegin3[] = { 3.83  };
2098   Float_t gapRangeEnd3[] =   { 3.86 };
2099   averageOK  = FitAverage(hist, 2, okRangeBegin3, okRangeEnd3, 1, bins);
2100   averageGap = FitAverage(hist, 1, gapRangeBegin3, gapRangeEnd3, 2, bins);
2101   add += bins * (averageOK - averageGap);
2102   Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
2103   
2104   Float_t okRangeBegin4[] = { 4.2, 4.5 };
2105   Float_t okRangeEnd4[] =   { 4.4, 4.7 };
2106   Float_t gapRangeBegin4[] = { 4.45  };
2107   Float_t gapRangeEnd4[] =   { 4.45 };
2108   averageOK  = FitAverage(hist, 2, okRangeBegin4, okRangeEnd4, 1, bins);
2109   averageGap = FitAverage(hist, 1, gapRangeBegin4, gapRangeEnd4, 2, bins);
2110   add += bins * (averageOK - averageGap);
2111   Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
2112   
2113   Float_t okRangeBegin5[] = { 5.4, 5.7 };
2114   Float_t okRangeEnd5[] =   { 5.6, 5.8 };
2115   Float_t gapRangeBegin5[] = { 5.63  };
2116   Float_t gapRangeEnd5[] =   { 5.67 };
2117   averageOK  = FitAverage(hist, 2, okRangeBegin5, okRangeEnd5, 1, bins);
2118   averageGap = FitAverage(hist, 1, gapRangeBegin5, gapRangeEnd5, 2, bins);
2119   add += bins * (averageOK - averageGap);
2120   Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
2121   
2122   Printf("This adds %.2f %% to the total number of tracklets (%f)", 100.0 * add / hist->Integral(), hist->Integral());
2123   c->SaveAs("gap1.png");
2124   
2125   add1 = add;
2126   total1 = hist->Integral();
2127
2128   if (all)
2129     return;
2130
2131   Printf("\nEta larger than 0 side");
2132   
2133   c = new TCanvas;
2134   TFile::Open("analysis_esd_raw.root");
2135   hist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("hist2", 41, binEnd);
2136   hist->Rebin(2);
2137   hist->SetStats(0);
2138   hist->Sumw2();
2139   hist->Draw("HIST");
2140   gPad->SetGridx();
2141   gPad->SetGridy();
2142   
2143   TFile::Open(mcFile);  
2144   mcHist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("mcHist", 41, binEnd);
2145   mcHist->Rebin(2);
2146   mcHist->SetLineColor(2);
2147   mcHist->Scale(hist->Integral() / mcHist->Integral());
2148   mcHist->Draw("SAME");
2149   
2150   add = 0;
2151   
2152   Float_t okRangeBegin[] = { 0.04, 0.67, 1.34 };
2153   Float_t okRangeEnd[] =   { 0.55, 1.24, 1.63 };
2154   Float_t gapRangeBegin[] = { 0.6, 1.27  };
2155   Float_t gapRangeEnd[] =   { 0.65, 1.32 };
2156   Float_t averageOK  = FitAverage(hist, 3, okRangeBegin, okRangeEnd, 1, bins);
2157   Float_t averageGap = FitAverage(hist, 2, gapRangeBegin, gapRangeEnd, 2, bins);
2158   add += bins * (averageOK - averageGap);
2159   Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
2160   
2161   Float_t okRangeBegin2[] = { 2.32, 2.65 };
2162   Float_t okRangeEnd2[] =   { 2.55, 3.2 };
2163   Float_t gapRangeBegin2[] = { 2.59, 3.3 };
2164   Float_t gapRangeEnd2[] =   { 2.61, 3.3 };
2165   averageOK  = FitAverage(hist, 2, okRangeBegin2, okRangeEnd2, 1, bins);
2166   averageGap = FitAverage(hist, 2, gapRangeBegin2, gapRangeEnd2, 2, bins);
2167   add += bins * (averageOK - averageGap);
2168   Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
2169   
2170   Float_t okRangeBegin3[] = { 3.55, 3.9 };
2171   Float_t okRangeEnd3[] =   { 3.8, 4.15 };
2172   Float_t gapRangeBegin3[] = { 3.83  };
2173   Float_t gapRangeEnd3[] =   { 3.86 };
2174   averageOK  = FitAverage(hist, 2, okRangeBegin3, okRangeEnd3, 1, bins);
2175   averageGap = FitAverage(hist, 1, gapRangeBegin3, gapRangeEnd3, 2, bins);
2176   add += bins * (averageOK - averageGap);
2177   Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
2178   
2179   Float_t okRangeBegin4[] = { 4.2, 4.5 };
2180   Float_t okRangeEnd4[] =   { 4.4, 4.7 };
2181   Float_t gapRangeBegin4[] = { 4.45  };
2182   Float_t gapRangeEnd4[] =   { 4.45 };
2183   averageOK  = FitAverage(hist, 2, okRangeBegin4, okRangeEnd4, 1, bins);
2184   averageGap = FitAverage(hist, 1, gapRangeBegin4, gapRangeEnd4, 2, bins);
2185   add += bins * (averageOK - averageGap);
2186   Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
2187
2188   Printf("This adds %.2f %% to the total number of tracklets (%f)", 100.0 * add / hist->Integral(), hist->Integral());
2189   c->SaveAs("gap2.png");
2190   
2191   Printf("In total we add %.2f %%.", 100.0 * (add1 + add) / (total1 + hist->Integral()));
2192 }