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