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