adding flag for non-field data
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / drawSystematicsNew.C
1
2 Int_t markers[] = {20,21,22,23,28,29};
3 Int_t colors[]  = {1,2,3,4,6,8,102};
4
5 void loadlibs()
6 {
7   gSystem->Load("libTree");
8   gSystem->Load("libVMC");
9
10   gSystem->Load("libSTEERBase");
11   gSystem->Load("libANALYSIS");
12   gSystem->Load("libPWG0base");
13 }
14
15 void InitPad()
16 {
17   if (!gPad)
18     return;
19
20   gPad->Range(0, 0, 1, 1);
21   gPad->SetLeftMargin(0.15);
22   //gPad->SetRightMargin(0.05);
23   //gPad->SetTopMargin(0.13);
24   //gPad->SetBottomMargin(0.1);
25
26   gPad->SetGridx();
27   gPad->SetGridy();
28 }
29
30 void DrawpiKpAndCombinedZOnly(Float_t upperPtLimit=0.99)
31 {
32   //gROOT->ProcessLine(".L drawPlots.C");
33   gSystem->Load("libPWG0base");
34
35   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "correction_map.root" };
36   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "dndeta_correction" };
37   const char* legendNames[] = { "#pi", "K", "p", "standard" };
38   Int_t folderCount = 3;
39
40   TH2F* h2DCorrections[4];
41   TH1F* h1DCorrections[4];
42   for (Int_t i=0; i<4; i++) {
43     TFile::Open(fileNames[i]);
44     AlidNdEtaCorrection* correctionTmp = new AlidNdEtaCorrection(folderNames[i],folderNames[i]);
45     correctionTmp->LoadHistograms();
46     
47     //    h2DCorrections[i] = correctionTmp->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get2DCorrectionHistogram("yz",-1,1);
48
49     h1DCorrections[i] = correctionTmp->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get1DCorrectionHistogram("z",-10,10,0.8,0.8);
50   }
51
52   TH2F* null = new TH2F("","",100,0.1,10,100,0.5,9.99);
53   null->SetXTitle("p_{T} (GeV/c)");
54   null->SetXTitle("Correction");
55
56   null->Draw();
57
58   //h1DCorrections[0]->SetMaximum(5);
59   //h1DCorrections[0]->Draw();
60   //  h2DCorrections[0]->Draw("colz");
61   for (Int_t i=1; i<4; i++) {
62     h1DCorrections[i]->Draw("same");
63   }
64   
65 }
66
67 TPad* DrawChange(Bool_t spd, const char* basename, const char** changes, Int_t nChanges, Int_t nDraw, Int_t* colors, const char** names = 0, Float_t scale = 0.10)
68 {  
69   Float_t etaMax = 1.05;
70   if (spd)
71     etaMax = 1.79;
72
73   TH1F* hRatios[100];
74   for(Int_t i=0; i<nChanges; i++) {
75     hRatios[i] = (TH1F*)gFile->Get(Form("%s%s",basename,changes[i]));
76     hRatios[i]->SetLineWidth(1);
77     hRatios[i]->SetMarkerStyle(22);
78     hRatios[i]->SetMarkerSize(0.8);
79
80     Float_t average = hRatios[i]->Integral(hRatios[i]->FindBin(-1), hRatios[i]->FindBin(1)) / (hRatios[i]->FindBin(1) - hRatios[i]->FindBin(-1) + 1);
81     Printf("%s: %.2f %%" , hRatios[i]->GetTitle(), (average - 1) * 100);
82   }
83   
84   TPad* p = DrawCanvasAndPad("syst_changeInXsection",700,400);
85   p->SetRightMargin(0.2);
86   p->SetLeftMargin(0.13);
87   
88   TH2F* null = new TH2F("","",100,-etaMax, etaMax,100,1. - scale,1. + scale);
89   null->GetXaxis()->SetTitle("#eta");
90   null->GetYaxis()->SetTitle(hRatios[0]->GetYaxis()->GetTitle());
91   null->Draw();
92   
93   line = new TLine(-etaMax, 1, etaMax, 1);
94   line->Draw();
95
96   TLatex* text[100];
97
98   for(Int_t i=1; i<nDraw; i++) {
99     hRatios[i]->SetLineColor(colors[i]);
100     hRatios[i]->SetMarkerColor(colors[i]);
101     hRatios[i]->Draw("HISTPL SAME");
102     
103     TString str(hRatios[i]->GetTitle());
104     
105     if (names)
106       str = names[i];
107     
108     text[i] = new TLatex(etaMax + 0.03,hRatios[i]->GetBinContent(hRatios[i]->FindBin(etaMax-0.1))-0.002,str.Data());
109     text[i]->SetTextAlign(11);
110     text[i]->SetTextColor(colors[i]);
111     text[i]->Draw();
112   }
113   
114   return p;
115 }
116
117 TPad* DrawChangeRatios(Bool_t spd, TFile* file1, TFile* file2, const char* basename, const char** changes, Int_t nChanges, Int_t nDraw, Int_t* colors, const char** names = 0, Float_t scale = 0.10)
118 {  
119   Float_t etaMax = 1.05;
120   if (spd)
121     etaMax = 1.79;
122
123   TH1F* hRatios[100];
124   for(Int_t i=0; i<nChanges; i++) {
125     hRatios[i] = (TH1F*)file1->Get(Form("%s%s",basename,changes[i]));
126     hRatios[i]->SetLineWidth(1);
127     hRatios[i]->SetMarkerStyle(22);
128     hRatios[i]->SetMarkerSize(0.8);
129     
130     hRatio2 = (TH1F*)file2->Get(Form("%s%s",basename,changes[i]));
131     hRatios[i]->Divide(hRatio2);
132
133     Float_t average = hRatios[i]->Integral(hRatios[i]->FindBin(-1), hRatios[i]->FindBin(1)) / (hRatios[i]->FindBin(1) - hRatios[i]->FindBin(-1) + 1);
134     Printf("%s: %.2f %%" , hRatios[i]->GetTitle(), (average - 1) * 100);
135   }
136   
137   TPad* p = DrawCanvasAndPad("syst_changeInXsection",700,400);
138   p->SetRightMargin(0.2);
139   p->SetLeftMargin(0.13);
140   
141   TH2F* null = new TH2F("","",100,-etaMax, etaMax,100,1. - scale,1. + scale);
142   null->GetXaxis()->SetTitle("#eta");
143   null->GetYaxis()->SetTitle(hRatios[0]->GetYaxis()->GetTitle());
144   null->Draw();
145   
146   line = new TLine(-etaMax, 1, etaMax, 1);
147   line->Draw();
148
149   TLatex* text[100];
150
151   for(Int_t i=1; i<nDraw; i++) {
152     hRatios[i]->SetLineColor(colors[i]);
153     hRatios[i]->SetMarkerColor(colors[i]);
154     hRatios[i]->Draw("HISTPL SAME");
155     
156     TString str(hRatios[i]->GetTitle());
157     
158     if (names)
159       str = names[i];
160     
161     text[i] = new TLatex(etaMax + 0.03,hRatios[i]->GetBinContent(hRatios[i]->FindBin(etaMax-0.1))-0.002,str.Data());
162     text[i]->SetTextAlign(11);
163     text[i]->SetTextColor(colors[i]);
164     text[i]->Draw();
165   }
166   
167   return p;
168 }
169
170 void DrawEffectOfChangeInCrossSection(Bool_t spd = kFALSE, const char* fileName = "systematics_vtxtrigger_compositions.root") 
171 {
172   TFile::Open(fileName);
173
174   const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25" };
175   //const Char_t* changes[]  = { "pythia", "qgsm", "phojet" };
176   //const Int_t nChanges = 3;
177   Int_t colors[] = {1,1,4,1,2,2,4,2,1};
178
179   c = DrawChange(spd, "ratio_vertexReco_triggerBias_", changes, 17, 9, colors, 0);
180   c->SaveAs("cross_sections.eps");
181 }
182
183 void DrawEffectOfChangeInCrossSection2Files(Bool_t spd = kFALSE, const char* fileName1 = "systematics_vtxtrigger_compositions.root", const char* fileName2) 
184 {
185   file1 = TFile::Open(fileName1);
186   file2 = TFile::Open(fileName2);
187
188   const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25" };
189   //const Char_t* changes[]  = { "pythia", "qgsm", "phojet" };
190   //const Int_t nChanges = 3;
191   Int_t colors[] = {1,1,4,1,2,2,4,2,1};
192
193   c = DrawChangeRatios(spd, file1, file2, "ratio_vertexReco_triggerBias_", changes, 17, 9, colors, 0);
194   c->SaveAs("cross_sections.eps");
195 }
196
197 void DrawEffectOfChangeInComposition(Bool_t spd = kFALSE, const char* fileName = "new_compositions_analysis.root") 
198 {
199   TFile::Open(fileName);
200
201   const Char_t* changes[]  = { "PythiaRatios", "KBoosted", "KReduced", "pBoosted", "pReduced", "KBoostedpBoosted", "KReducedpReduced", "KBoostedpReduced", "KReducedpBoosted"};
202   const char*   names[]    = { "",             "K #times 1.5", "K #times 0.5", "p #times 1.5", "p #times 0.5", "K #times 1.5, p #times 1.5", "K #times 0.5, p #times 0.5", "K #times 1.5, p #times 0.5", "K #times 0.5, p #times 1.5" };
203   //const Char_t* changes[]  = { "PythiaRatios", "PiBoosted",      "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced", "othersBoosted", "othersReduced" };
204   //const char*   names[]    = { "",             "#pi #times 1.5", "#pi #times 0.5", "K #times 1.5", "K #times 0.5", "p #times 1.5", "p #times 0.5",  "others #times 1.5", "others #times 0.5" };
205   Int_t colors[] = {1,1,2,2,1,2,1,4,4};
206
207   c = DrawChange(spd, "", changes, 9, 9, colors, names, 0.03);
208   c->SaveAs("compositions.eps");
209 }
210
211 TPad* DrawCanvasAndPad(const Char_t* name, Int_t sizeX=600, Int_t sizeY=500) {
212
213   gStyle->SetOptStat(0);
214   gStyle->SetOptTitle(0);
215   gStyle->SetOptFit(0);
216
217   gStyle->SetTextSize(0.04);
218   gStyle->SetTitleSize(0.05,"xyz");
219   //gStyle->SetTitleFont(133, "xyz");
220   //gStyle->SetLabelFont(133, "xyz");
221   //gStyle->SetLabelSize(17, "xyz");
222   gStyle->SetLabelOffset(0.01, "xyz");
223
224   gStyle->SetTitleOffset(1.1, "y");
225   gStyle->SetTitleOffset(1.1, "x");
226   gStyle->SetEndErrorSize(0.0);
227
228   //##############################################
229
230   //making canvas and pads
231   TCanvas *c = new TCanvas(name,name,sizeX,sizeY);
232
233   TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
234
235   p1->SetBottomMargin(0.15);
236   p1->SetTopMargin(0.03);
237   p1->SetLeftMargin(0.15);
238   p1->SetRightMargin(0.03);
239   
240   p1->SetGridx();
241   p1->SetGridy();
242
243   p1->Draw();
244   p1->cd();
245
246   return p1;
247 }
248
249 void MisalignmentShowRawTrackPlots(const char* dirName = "fdNdEtaAnalysisESD")
250 {
251   loadlibs();
252
253   TFile* file = TFile::Open("fullA-simrec/MB2/analysis_esd_raw.root");
254   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
255   fdNdEtaAnalysis->LoadHistograms();
256
257   TFile* file2 = TFile::Open("fullA-sim/analysis_esd_raw.root");
258   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis(dirName, dirName);
259   fdNdEtaAnalysis2->LoadHistograms();
260
261   TH3* track1 = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track1");
262   TH3* track2 = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track2");
263
264   // normalize to number of events;
265   TH2* event1 = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
266   TH2* event2 = fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram();
267   Int_t event1Count = event1->Integral();
268   Int_t event2Count = event2->Integral();
269   track1->Scale(1.0 / event1Count);
270   track2->Scale(1.0 / event2Count);
271
272   const Float_t innerLimit = 0.49;
273   const Float_t outerLimit = 0.99;
274
275   track1->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
276   track2->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
277   AliPWG0Helper::CreateDividedProjections(track1, track2, "ze1");
278   TH1* fullRange = gROOT->FindObject("track1_ze1_div_track2_ze1");
279
280   track1->GetYaxis()->SetRangeUser(-innerLimit, innerLimit);
281   track2->GetYaxis()->SetRangeUser(-innerLimit, innerLimit);
282   AliPWG0Helper::CreateDividedProjections(track1, track2, "ze2");
283   TH1* central = gROOT->FindObject("track1_ze2_div_track2_ze2");
284   central->SetLineColor(1);
285   central->SetMarkerStyle(21);
286
287   for (Int_t x=1; x<track1->GetXaxis()->GetNbins(); ++x)
288     for (Int_t y=track1->GetYaxis()->FindBin(-innerLimit); y<track1->GetYaxis()->FindBin(innerLimit); ++y)
289       for (Int_t z=1; z<track1->GetZaxis()->GetNbins(); ++z)
290       {
291         track1->SetBinContent(x, y, z, 0);
292         track1->SetBinError(x, y, z, 0);
293         track2->SetBinContent(x, y, z, 0);
294         track2->SetBinError(x, y, z, 0);
295       }
296
297   track1->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
298   track2->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
299   AliPWG0Helper::CreateDividedProjections(track1, track2, "ze3");
300   TH1* peripheral = gROOT->FindObject("track1_ze3_div_track2_ze3");
301   peripheral->SetLineColor(2);
302   peripheral->SetMarkerStyle(22);
303   peripheral->SetMarkerColor(2);
304
305   TH2* tmp = new TH2F("tmp", ";p_{T} [GeV/c]      ;#frac{tracks full misalignment during rec.}{tracks ideal misalignment during rec.}", 1, 0.1, 10, 1, 0.9, 1.3);
306
307   tmp->SetStats(kFALSE);
308   //tmp->GetXaxis()->SetNoExponent();
309
310   Float_t ptStart = 0.1;
311
312   fullRange->GetXaxis()->SetRangeUser(ptStart, 9.9);
313   central->GetXaxis()->SetRangeUser(ptStart, 9.9);
314   peripheral->GetXaxis()->SetRangeUser(ptStart, 9.9);
315
316   TCanvas* canvas = new TCanvas("MisalignmentShowRawTrackPlots", "MisalignmentShowRawTrackPlots", 700, 400);
317   gPad->SetLogx();
318   gPad->SetGridx();
319   gPad->SetGridy();
320
321   TLegend* legend = new TLegend(0.2, 0.7, 0.4, 0.8);
322
323   legend->AddEntry(central, "|#eta| < 0.5");
324   legend->AddEntry(peripheral, "0.5 < |#eta| < 1.0");
325
326   legend->SetFillColor(0);
327
328   tmp->Draw();
329   //fullRange->Draw("SAME");
330   central->Draw("SAME");
331   peripheral->Draw("SAME");
332
333   legend->Draw();
334
335   canvas->SaveAs("syst_mis_ntracks.eps");
336 }
337
338
339 void drawdNdEtaRatios(const char* canvasName, Int_t n, const char** files, const char** dirs, const char** names, Int_t* histID)
340 {
341   loadlibs();
342   gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/dNdEta/drawPlots.C");
343
344   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
345   canvas->Divide(2, 1);
346   canvas->cd(2)->SetGridx();
347   canvas->cd(2)->SetGridy();
348
349   TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
350   legend->SetFillColor(0);
351
352   TH1* base = 0;
353
354   for (Int_t i = 0; i < n; ++i)
355   {
356     TFile::Open(files[i]);
357
358     dNdEtaAnalysis* tmp = new dNdEtaAnalysis(dirs[i], dirs[i]);
359     tmp->LoadHistograms();
360
361     TH1* hist = tmp->GetdNdEtaPtCutOffCorrectedHistogram(histID[i]);
362
363     if (i == 0)
364       base = hist;
365
366     legend->AddEntry(hist, names[i]);
367
368     hist->SetMarkerColor(colors[i]);
369     hist->SetMarkerStyle(markers[i]);
370
371     canvas->cd(1);
372     hist->DrawCopy((i == 0) ? "" : "SAME");
373
374     if (i != 0)
375     {
376       PrintIntegratedDeviation(hist, base, names[i]);
377
378       canvas->cd(2);
379       hist->Divide(hist, base, 1, 1);
380       hist->GetYaxis()->SetRangeUser(0.9, 1.1);
381       hist->DrawCopy((i == 1) ? "" : "SAME");
382     }
383   }
384
385   legend->Draw();
386
387   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
388 }
389
390 void drawdNdEtaRatios(const char* canvasName, Int_t n, const char** files1, const char** files2, const char** dirs1, const char** dirs2, const char** names, Int_t* histID)
391 {
392   gSystem->Load("libPWG0base");
393
394   TCanvas* canvas = new TCanvas(canvasName, canvasName, 700, 400);
395   canvas->SetLeftMargin(0.12);
396
397   TLegend* legend = new TLegend(0.35, 0.7, 0.65, 0.85);
398   legend->SetFillColor(0);
399
400   for (Int_t i = 0; i < n; ++i)
401   {
402     TFile::Open(files1[i]);
403
404     dNdEtaAnalysis* tmp1 = new dNdEtaAnalysis(dirs1[i], dirs1[i]);
405     tmp1->LoadHistograms();
406
407     TFile::Open(files2[i]);
408
409     dNdEtaAnalysis* tmp2 = new dNdEtaAnalysis(dirs2[i], dirs2[i]);
410     tmp2->LoadHistograms();
411
412     TH1* hist1 = tmp1->GetdNdEtaPtCutOffCorrectedHistogram(histID[i]);
413     TH1* hist2 = tmp2->GetdNdEtaPtCutOffCorrectedHistogram(histID[i]);
414
415     TH1* division = hist1->Clone();
416
417     division->Divide(hist1, hist2, 1, 1, "B");
418
419     division->SetMarkerColor(colors[i]);
420     division->SetMarkerStyle(markers[i]);
421
422     legend->AddEntry(division, names[i]);
423
424     division->SetTitle("");
425     division->GetYaxis()->SetTitle("#frac{dN_{ch}/d#eta using MC vtx}{dN_{ch}/d#eta using ESD vtx}");
426     division->SetStats(kFALSE);
427     division->GetYaxis()->SetTitleOffset(1.3);
428     division->GetXaxis()->SetRangeUser(-0.99, 0.99);
429     division->GetYaxis()->SetRangeUser(0.981, 1.02);
430     division->DrawCopy((i == 0) ? "" : "SAME");
431   }
432
433   gPad->SetGridx();
434   gPad->SetGridy();
435
436   legend->Draw();
437
438   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
439 }
440
441 void MaterialBudgetChange()
442 {
443   const char* files[] =
444     { "Material-normal-mcvtx/analysis_esd.root",
445       "Material-increased-mcvtx/analysis_esd.root",
446       "Material-decreased-mcvtx/analysis_esd.root" };
447
448   const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
449   const char* names[] = { "no change", "+ 10 % material", "- 10 % material" };
450   Int_t hist[] = { 0, 0, 0 };
451
452   drawdNdEtaRatios("MaterialBudgetChange", 3, files, dirs, names, hist);
453 }
454
455 void MisalignmentChange()
456 {
457   const char* files[] =
458     { "maps/v4-09-Release/tpc-only/fullC-simrec--fullA-sim-mcvtx/analysis_esd.root",
459       "maps/v4-09-Release/tpc-only/fullC-simrec--fullA-simrec-mcvtx/analysis_esd.root" };
460
461   const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
462   const char* names[] = { "no change", "increased in sim/rec", "increased in sim" };
463   Int_t hist[] = { 0, 0, 0 };
464
465   drawdNdEtaRatios("MisalignmentChange", 2, files, dirs, names, hist);
466 }
467
468 void dNdEtaVertexRanges()
469 {
470   const char* files[] =
471     { "analysis_esd.root",
472       "analysis_esd.root",
473       "analysis_esd.root" };
474
475   const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
476   const char* names[] = { "|vtx-z| < 10 cm", "-10 cm < vtx-z < 0 cm", "0 cm < vtx-z < 10 cm" };
477   Int_t hist[] = { 0, 1, 2 };
478
479   drawdNdEtaRatios("dNdEtaVertexRanges", 3, files, dirs, names, hist);
480 }
481
482 void vertexShiftStudy(Int_t histID)
483 {
484   const char* files[] = { "maps/idealA/mc-vertex/analysis_esd.root", "results/idealC-idealA/analysis_esd.root", "maps/idealA/mc-vertex-shift-0.05/analysis_esd.root", "maps/idealA/mc-vertex-shift-0.1/analysis_esd.root", "maps/idealA/mc-vertex-shift-dep/analysis_esd.root" };
485   const char* dirs[] = { "dndeta", "dndeta", "dndeta", "dndeta", "dndeta" };
486   const char* names[] = { "mc vtx", "esd vtx", "+ 0.05 cm", "+ 0.1 cm", "old vtx shift" };
487   Int_t hist[] = { histID, histID, histID, histID, histID };
488
489   drawdNdEtaRatios("syst_vertex_shift1", 5, files, dirs, names, hist);
490
491   const char* files1[] = { "maps/idealA/mc-vertex/analysis_esd.root", "maps/idealA/mc-vertex/analysis_esd.root", "maps/idealA/mc-vertex/analysis_esd.root" };
492   const char* files2[] = { "results/idealC-idealA/analysis_esd.root", "results/idealC-idealA/analysis_esd.root", "results/idealC-idealA/analysis_esd.root" };
493   const char* dirs1[] = { "dndeta", "dndeta", "dndeta"};
494   const char* names[] = { "|vtx-z| < 10 cm", "-10 cm < vtx-z < 0 cm", "0 cm < vtx-z < 10 cm" };
495   Int_t hist[] = { 0, 1, 2 };
496
497   drawdNdEtaRatios("syst_vertex_shift2", 3, files1, files2, dirs, dirs, names, hist);
498 }
499
500 void vertexShift()
501 {
502   TFile::Open("vertex.root");
503
504   TH2* hist = gFile->Get("fVertexCorr");
505   TProfile* prof = hist->ProfileX();
506
507   prof->SetStats(kFALSE);
508   prof->SetTitle(";MC vtx-z [cm];mean (ESD vtx-z - MC vtx-z) [cm]");
509   prof->GetYaxis()->SetTitleOffset(1.2);
510
511   prof->SetLineWidth(2);
512
513   TCanvas* canvas = new TCanvas("syst_vertex_shift", "syst_vertex_shift", 700, 400);
514
515   gPad->SetGridx();
516   gPad->SetGridy();
517
518   prof->Draw();
519
520   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
521 }
522
523 void CompareRawTrackPlots(const char* fileName1, const char* fileName2, Float_t ptCut = 0.0, Int_t multCut = 1)
524 {
525   loadlibs();
526
527   const char* dirName = "fdNdEtaAnalysisESD";
528
529   TFile* file = TFile::Open(fileName1);
530   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
531   fdNdEtaAnalysis->LoadHistograms();
532
533   TFile* file2 = TFile::Open(fileName2);
534   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis(dirName, dirName);
535   fdNdEtaAnalysis2->LoadHistograms();
536
537   TH3* track1 = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track1");
538   TH3* track2 = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track2");
539
540   TH2* event1 = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
541   TH2* event2 = fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram();
542   Int_t event1Count = event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), multCut, event1->GetNbinsY() + 1);
543   Int_t event2Count = event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), multCut, event1->GetNbinsY() + 1);
544
545   Float_t nTrack1 = track1->Integral(track1->GetXaxis()->FindBin(-9.9), track1->GetXaxis()->FindBin(9.9), track1->GetYaxis()->FindBin(-0.79), track1->GetYaxis()->FindBin(0.79), track1->GetZaxis()->FindBin(ptCut), track1->GetZaxis()->GetNbins());
546   Float_t nTrack2 = track2->Integral(track2->GetXaxis()->FindBin(-9.9), track2->GetXaxis()->FindBin(9.9), track2->GetYaxis()->FindBin(-0.79), track2->GetYaxis()->FindBin(0.79), track2->GetZaxis()->FindBin(ptCut), track2->GetZaxis()->GetNbins());
547
548   Printf("%d tracks in %d events in first sample; %d tracks in %d events in second sample", (Int_t) nTrack1, event1Count, (Int_t) nTrack2, event2Count);
549
550   // normalize to number of events;
551   nTrack1 /= event1Count;
552   nTrack2 /= event2Count;
553
554   Printf("There are %.2f tracks/event in the first sample and %.2f tracks/event in the second sample. %.2f %% difference (with pt cut at %.2f GeV/c)", nTrack1, nTrack2, 100.0 * (nTrack1 - nTrack2) / nTrack1, ptCut);
555
556   gROOT->cd();
557
558   AliPWG0Helper::CreateDividedProjections(track1, track2);
559
560   new TCanvas; gROOT->FindObject("track1_yx_div_track2_yx")->Draw("COLZ");
561   new TCanvas; gROOT->FindObject("track1_zx_div_track2_zx")->Draw("COLZ");
562   new TCanvas; gROOT->FindObject("track1_zy_div_track2_zy")->Draw("COLZ");
563
564   for (Int_t i=0; i<3; i++)
565   {
566     char c = 'x' + (char) i;
567
568     /*proj1 = track1->Project3D(Form("%ce2", c));
569     proj2 = track2->Project3D(Form("%ce2", c));
570     AliPWG0Helper::NormalizeToBinWidth(proj1);
571     AliPWG0Helper::NormalizeToBinWidth(proj2);
572
573     new TCanvas;
574     proj1->DrawCopy();
575     proj2->SetLineColor(2);
576     proj2->SetMarkerColor(2);
577     proj2->DrawCopy("SAME");*/
578
579     AliPWG0Helper::CreateDividedProjections(track1, track2, Form("%ce2", c));
580     TH1* pt = gROOT->FindObject(Form("track1_%ce2_div_track2_%ce2", c, c));
581     new TCanvas; pt->DrawCopy();
582     gPad->SetGridx(); gPad->SetGridy();
583   }
584
585   event1_x = event1->ProjectionX("event1_x");
586   event1_y = event1->ProjectionY("event1_y");
587   event2_x = event2->ProjectionX("event2_x");
588   event2_y = event2->ProjectionY("event2_y");
589
590   new TCanvas; event1_x->DrawCopy(); event2_x->SetLineColor(2); event2_x->DrawCopy("SAME"); 
591
592   event1_x->Divide(event2_x);
593   event1_y->Divide(event2_y);
594
595   new TCanvas; event1_x->DrawCopy();
596   new TCanvas; event1_y->DrawCopy();
597
598   event1->Divide(event2);
599   new TCanvas;
600   event1->Draw("COLZ");
601   event1->SetMinimum(0.5);
602   event1->SetMaximum(2);
603
604
605 }
606
607 void MagnitudeOfCorrection(const char* fileName, const char* dirName = "dndeta", Float_t ptCut = 0.3)
608 {
609   loadlibs();
610
611   TFile* file = TFile::Open(fileName);
612   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
613   fdNdEtaAnalysis->LoadHistograms();
614   fdNdEtaAnalysis->GetData()->PrintInfo(ptCut);
615 }
616
617 Double_t ConvSigma1To2D(Double_t sigma)
618 {
619   return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
620 }
621
622 Double_t ConvDistance1DTo2D(Double_t distance)
623 {
624   return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
625 }
626
627 Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
628 {
629   Double_t count = 0;
630
631   //nSigma = ConvSigma1To2D(nSigma);
632
633   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
634     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
635     {
636       Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
637       Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
638
639       Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
640
641       d = ConvDistance1DTo2D(d);
642
643       if (d < nSigma)
644         count += tracks->GetBinContent(x, y);
645     }
646
647   return count;
648 }
649
650 TH2F* Sigma2VertexGaussianTracksHist()
651 {
652   TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
653
654   TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
655   gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
656
657   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
658     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
659       tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
660
661   //normalize
662   tracks->Scale(1.0 / tracks->Integral());
663
664   return tracks;
665 }
666
667 TH1F* Sigma2VertexGaussian()
668 {
669   TH2F* tracks = Sigma2VertexGaussianTracksHist();
670
671   TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
672   canvas->Divide(2, 2);
673
674   canvas->cd(1);
675   tracks->Draw("COLZ");
676
677   TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 50, 0.05, 5.05);
678   for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
679     ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
680   ratio->SetMarkerStyle(21);
681
682   canvas->cd(2);
683   ratio->DrawCopy("P");
684
685   TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 4 sigma / % included n sigma", 50, 0.05, 5.05);
686   Double_t sigma3 = Sigma2VertexCount(tracks, 4);
687   for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
688     ratio2->Fill(nSigma, sigma3 / ratio->GetBinContent(ratio->FindBin(nSigma)));
689   ratio2->SetMarkerStyle(21);
690
691   canvas->cd(3);
692   ratio2->DrawCopy("P");
693
694   canvas->SaveAs("Sigma2Vertex.eps");
695
696   return ratio2;
697 }
698
699 TH1F** Sigma2VertexSimulation(const char* fileName = "systematics.root")
700 {
701   TFile* file = TFile::Open(fileName);
702
703   TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertexTracks"));
704   TH1F* sigmavertexPrim = dynamic_cast<TH1F*> (file->Get("fSigmaVertexPrim"));
705   if (!sigmavertex || !sigmavertexPrim)
706   {
707     printf("Could not read histogram(s)\n");
708     return;
709   }
710
711   // calculate ratio
712   TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;N#sigma;% included in 4 #sigma / % included in N#sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax());
713
714   // calculate contamination
715   TH1F* contamination = ratio->Clone("sigmavertexsimulation_contamination");
716   contamination->SetTitle("sigmavertexsimulation_contamination;N#sigma;1 + N_{secondaries} / N_{all}");
717
718   for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
719   {
720     ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(4)) / sigmavertex->GetBinContent(i));
721     contamination->SetBinContent(i, 1 + (sigmavertex->GetBinContent(i) - sigmavertexPrim->GetBinContent(i)) / sigmavertex->GetBinContent(i));
722   }
723
724   // print stats
725   for (Float_t sigma = 2.0; sigma < 5.25; sigma += 0.5)
726   {
727     Float_t error1 = 1 - ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma - 0.5));
728     Float_t error2 = -1 + ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma + 0.5));
729     Float_t cont = -1 + contamination->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma));
730     Printf("%.2f sigma --> syst. error = - %.2f %% + %.2f %%, cont. = %.2f %%", sigma, error1 * 100, error2 * 100, cont * 100);
731   }
732
733   TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
734   canvas->Divide(2, 1);
735
736   canvas->cd(1);
737   sigmavertex->SetMarkerStyle(21);
738   sigmavertex->Draw("P");
739
740   canvas->cd(2);
741   ratio->SetMarkerStyle(21);
742   ratio->DrawCopy("P");
743
744   contamination->DrawCopy("SAME");
745
746   TH1F** returnContainer = new TH1F*[2];
747   returnContainer[0] = ratio;
748   returnContainer[1] = contamination;
749
750   return returnContainer;
751 }
752
753 void Sigma2VertexCompare(const char* fileName = "systematics.root")
754 {
755   TH1F* ratio1 = Sigma2VertexGaussian();
756
757   TH1F** hists = Sigma2VertexSimulation(fileName);
758   TH1F* ratio2 = hists[0];
759   TH1F* contamination = hists[1];
760
761   ratio1->SetStats(kFALSE);
762   ratio2->SetStats(kFALSE);
763
764   ratio1->SetMarkerStyle(0);
765   ratio2->SetMarkerStyle(0);
766
767   ratio1->SetLineWidth(2);
768   ratio2->SetLineWidth(2);
769
770   TLegend* legend = new TLegend(0.6, 0.8, 0.95, 0.95);
771   legend->SetFillColor(0);
772   legend->AddEntry(ratio1, "Gaussian");
773   legend->AddEntry(ratio2, "Simulation");
774   legend->AddEntry(contamination, "1 + Contamination");
775
776   ratio2->SetTitle("");
777   ratio2->GetYaxis()->SetTitleOffset(1.5);
778   ratio2->GetXaxis()->SetRangeUser(2, 5);
779
780   TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
781   InitPad();
782
783   ratio2->SetMarkerStyle(21);
784   ratio1->SetMarkerStyle(22);
785
786   ratio2->GetYaxis()->SetRangeUser(0.8, 1.2);
787   ratio2->SetLineColor(kRed);
788   ratio2->SetMarkerColor(kRed);
789   ratio2->Draw("PL");
790   ratio1->Draw("SAMEPL");
791
792   contamination->Draw("SAME");
793
794   legend->Draw();
795
796   canvas->SaveAs("Sigma2VertexCompare.eps");
797 }