]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG0/dNdEta/drawSystematics.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
adding trigger classes
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / drawSystematics.C
CommitLineData
10ebe68d 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>
6b7fa615 16#include <TLegend.h>
17#include <TPad.h>
18#include <TF1.h>
19
20extern TPad* gPad;
21
22void Track2Particle1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 10);
10ebe68d 23
24#endif
25
8ca1a6d9 26Int_t markers[] = {20,20,21,22,23,28,29};
27Int_t colors[] = {1,2,3,4,6,8,102};
9f469bf5 28
dd367a14 29void 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
10ebe68d 39void 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
51void SetRanges(TH1* hist)
52{
53 SetRanges(hist->GetXaxis());
54 SetRanges(hist->GetYaxis());
55 SetRanges(hist->GetZaxis());
56}
57
58void 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
67void Prepare2DPlot(TH2* hist)
68{
69 hist->SetStats(kFALSE);
70 hist->GetYaxis()->SetTitleOffset(1.4);
71
72 SetRanges(hist);
73}
74
7af955da 75void Prepare1DPlot(TH1* hist, Bool_t setRanges = kTRUE)
10ebe68d 76{
77 hist->SetLineWidth(2);
78 hist->SetStats(kFALSE);
79
7af955da 80 hist->GetXaxis()->SetTitleOffset(1.2);
81 hist->GetYaxis()->SetTitleOffset(1.2);
82
83 if (setRanges)
84 SetRanges(hist);
10ebe68d 85}
86
87void InitPad()
88{
9f469bf5 89 if (!gPad)
90 return;
91
10ebe68d 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
72e597d7 98 gPad->SetGridx();
99 gPad->SetGridy();
10ebe68d 100}
101
102void 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
112void Secondaries()
113{
114 TFile* file = TFile::Open("systematics.root");
115
7af955da 116 TH2F* secondaries = dynamic_cast<TH2F*> (file->Get("fSecondaries"));
10ebe68d 117 if (!secondaries)
118 {
119 printf("Could not read histogram\n");
120 return;
121 }
122
7af955da 123 TCanvas* canvas = new TCanvas("Secondaries", "Secondaries", 1000, 1000);
124 canvas->Divide(3, 3);
125 for (Int_t i=1; i<=8; i++)
10ebe68d 126 {
7af955da 127 TH1D* hist = secondaries->ProjectionY(Form("proj_%d", i), i, i);
128 hist->SetTitle(secondaries->GetXaxis()->GetBinLabel(i));
10ebe68d 129
7af955da 130 canvas->cd(i);
131 hist->Draw();
10ebe68d 132 }
133}
134
6b7fa615 135void Track2Particle1DComposition(const char** fileNames, Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9)
10ebe68d 136{
137 gSystem->Load("libPWG0base");
138
9f469bf5 139 TString canvasName;
140 canvasName.Form("Track2Particle1DComposition");
141 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
142 canvas->Divide(3, 1);
10ebe68d 143
9f469bf5 144 TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);
10ebe68d 145
9f469bf5 146 for (Int_t i=0; i<folderCount; ++i)
10ebe68d 147 {
69b09e3b 148 Correction1DCreatePlots(fileNames[i], folderNames[i], upperPtLimit);
9f469bf5 149
69b09e3b 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"));
9f469bf5 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]);
10ebe68d 182 }
183
9f469bf5 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
6b7fa615 190TH1** 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
7af955da 204 TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
6b7fa615 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);
7af955da 245 ptDists[i]->SetStats(kFALSE);
246 ptDists[i]->SetTitle(";p_{T};Fraction of total");
6b7fa615 247 ptDists[i]->GetYaxis()->SetRangeUser(0, 1);
248 ptDists[i]->Draw((i == 0) ? "" : "SAME");
249 }
7af955da 250 legend->SetFillColor(0);
6b7fa615 251 legend->Draw();
252
253 canvas->SaveAs("DrawRatios.gif");
254
7af955da 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
6b7fa615 267 file->Close();
268
269 return ptDists;
270}
271
7af955da 272void 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
9f469bf5 282void DrawDifferentSpecies()
283{
284 gROOT->ProcessLine(".L drawPlots.C");
285
6b7fa615 286 const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "systematics.root" };
9f469bf5 287 const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
288
6b7fa615 289 Track2Particle1DComposition(fileNames, 4, folderNames);
9f469bf5 290}
10ebe68d 291
8ca1a6d9 292void 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
6b7fa615 302void 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
312void 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
322void 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
7af955da 332TH2F* HistogramDifferences(const char* filename, const char* folder1, const char* folder2)
6b7fa615 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
6b7fa615 346 TH3F* hist1 = fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
347 TH3F* hist2 = fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
348
7af955da 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
6b7fa615 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)
7af955da 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);
6b7fa615 359
360 printf("Over-/Underflow bins: %d %d\n", difference->GetBinContent(0), difference->GetBinContent(difference->GetNbinsX()+1));
361
362 return difference;
363}
364
365void HistogramDifferences()
366{
7af955da 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");
6b7fa615 371
7af955da 372 TCanvas* canvas = new TCanvas("HistogramDifferences", "HistogramDifferences", 1000, 1000);
373 canvas->Divide(2, 2);
6b7fa615 374
375 canvas->cd(1);
7af955da 376 KBoosted->GetXaxis()->SetRangeUser(-0.05, 0.05);
377 KBoosted->Draw("COLZ");
6b7fa615 378
379 canvas->cd(2);
7af955da 380 KReduced->GetXaxis()->SetRangeUser(-0.05, 0.05);
381 KReduced->Draw("COLZ");
6b7fa615 382
383 canvas->cd(3);
7af955da 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");
6b7fa615 390
391 canvas->SaveAs("HistogramDifferences.gif");
7af955da 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");
6b7fa615 439}
440
441
9f469bf5 442void ScalePtDependent(TH3F* hist, TF1* function)
443{
444 // assumes that pt is the third dimension of hist
445 // scales with function(pt)
10ebe68d 446
9f469bf5 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
6b7fa615 458void 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
474const char* ChangeComposition(void** correctionPointer, Int_t index)
9f469bf5 475{
476 AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer;
477
478 switch (index)
479 {
6b7fa615 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
7af955da 496 /*TH1** ptDists = DrawRatios("pythiaratios.root");
6b7fa615 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]);
7af955da 502 }*/
6b7fa615 503 return "PythiaRatios";
9f469bf5 504 break;
505
7af955da 506 // one species enhanced / reduced
69b09e3b 507 case 2: // + 50% kaons
508 case 3: // - 50% kaons
509 case 4: // + 50% protons
510 case 5: // - 50% protons
511 case 6: // + 50% kaons + 50% protons
512 case 7: // - 50% kaons - 50% protons
513 case 8: // + 50% kaons - 50% protons
514 case 9: // - 50% kaons + 50% protons
7af955da 515 TString* str = new TString;
69b09e3b 516 if (index < 6)
517 {
518 Int_t correctionIndex = index / 2;
519 Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
520
521 fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
522 str->Form("%s%s", (correctionIndex == 0) ? "Pi" : ((correctionIndex == 1) ? "K" : (correctionIndex == 2) ? "p" : "others"), (index % 2 == 0) ? "Boosted" : "Reduced");
523 }
524 else
525 {
526 Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
527 fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
528 str->Form("%s%s", "K", (scaleFactor > 1) ? "Boosted" : "Reduced");
529
530 if (index >= 8)
531 scaleFactor = (index % 2 == 0) ? 0.5 : 1.5;
532 fdNdEtaCorrection[2]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
533 *str += Form("%s%s", "p", (scaleFactor > 1) ? "Boosted" : "Reduced");
534 }
535
7af955da 536 return str->Data();
537 break;
538
6b7fa615 539 // each species rated with pythia ratios
7af955da 540 case 12: // + 50% pions
541 case 13: // - 50% pions
542 case 14: // + 50% kaons
543 case 15: // - 50% kaons
544 case 16: // + 50% protons
545 case 17: // - 50% protons
6b7fa615 546 TH1** ptDists = DrawRatios("pythiaratios.root");
547 Int_t functionIndex = (index - 2) / 2;
7af955da 548 Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
6b7fa615 549 ptDists[functionIndex]->Scale(scaleFactor);
550
551 for (Int_t i=0; i<3; ++i)
552 {
69b09e3b 553 ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram(), ptDists[i]);
554 ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram(), ptDists[i]);
6b7fa615 555 }
556 TString* str = new TString;
557 str->Form("%s%s", (functionIndex == 0) ? "Pi" : ((functionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced");
558 return str->Data();
9f469bf5 559 break;
560
6b7fa615 561 case 999:
9f469bf5 562 TF1* ptDependence = new TF1("simple", "x", 0, 100);
69b09e3b 563 ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram(), ptDependence);
564 ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram(), ptDependence);
9f469bf5 565 break;
566
567 }
6b7fa615 568
569 return "noname";
9f469bf5 570}
571
572void Composition()
573{
69b09e3b 574 loadlibs();
9f469bf5 575
576 gSystem->Unlink("new_compositions.root");
69b09e3b 577 gSystem->Unlink("new_compositions_analysis.root");
578
579 const char* names[] = { "pi", "K", "p", "other" };
580 TH1* hRatios[20];
9f469bf5 581
69b09e3b 582 Int_t nCompositions = 10;
583 Int_t counter = 0;
584 for (Int_t comp = 1; comp < nCompositions; ++comp)
9f469bf5 585 {
586 AlidNdEtaCorrection* fdNdEtaCorrection[4];
587
69b09e3b 588 TFile::Open("correction_mapparticle-species.root");
9f469bf5 589
590 for (Int_t i=0; i<4; ++i)
591 {
592 TString name;
69b09e3b 593 name.Form("dndeta_correction_%s", names[i]);
9f469bf5 594 fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
69b09e3b 595 fdNdEtaCorrection[i]->LoadHistograms();
9f469bf5 596 }
597
6b7fa615 598 const char* newName = ChangeComposition(fdNdEtaCorrection, comp);
9f469bf5 599
600 Double_t geneCount[5];
601 Double_t measCount[5];
602 geneCount[4] = 0;
603 measCount[4] = 0;
604
605 for (Int_t i=0; i<4; ++i)
606 {
69b09e3b 607 geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram()->Integral();
608 measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Integral();
9f469bf5 609
610 geneCount[4] += geneCount[i];
611 measCount[4] += measCount[i];
612
613 printf("Particle %s: %d gene, %d meas\n", ((i == 0) ? "pi" : (i == 1) ? "K" : (i == 2) ? "p" : "others"), (Int_t) geneCount[i], (Int_t) measCount[i]);
614 }
10ebe68d 615
9f469bf5 616 printf("Generated ratios are: %f pi, %f K, %f p, %f others\n", geneCount[0] / geneCount[4], geneCount[1] / geneCount[4], geneCount[2] / geneCount[4], geneCount[3] / geneCount[4]);
10ebe68d 617
9f469bf5 618 printf("Reconstructed ratios are: %f pi, %f K, %f p, %f others\n", measCount[0] / measCount[4], measCount[1] / measCount[4], measCount[2] / measCount[4], measCount[3] / measCount[4]);
10ebe68d 619
9f469bf5 620 TList* collection = new TList;
10ebe68d 621
69b09e3b 622 // skip "other" particle correction here
623 // with them has to be dealt differently, maybe just increasing the neutral particles...
624 for (Int_t i=1; i<3; ++i)
9f469bf5 625 collection->Add(fdNdEtaCorrection[i]);
10ebe68d 626
69b09e3b 627 fdNdEtaCorrection[0]->Merge(collection);
628 fdNdEtaCorrection[0]->Finish();
9f469bf5 629
630 delete collection;
631
69b09e3b 632 // save everything
9f469bf5 633 TFile* file = TFile::Open("new_compositions.root", "UPDATE");
69b09e3b 634 fdNdEtaCorrection[0]->SetName(newName);
635 fdNdEtaCorrection[0]->SaveHistograms();
9f469bf5 636 //file->Write();
637 file->Close();
69b09e3b 638
639 // correct dNdeta distribution with modified correction map
640 TFile::Open("analysis_esd_raw.root");
641
642 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
643 fdNdEtaAnalysis->LoadHistograms();
644
645 fdNdEtaAnalysis->Finish(fdNdEtaCorrection[0], 0.2, 3, newName);
646
647 hRatios[counter] = (TH1F*) fdNdEtaAnalysis->GetdNdEtaHistogram()->Clone(newName);
648 hRatios[counter]->SetTitle(newName);
649 hRatios[counter]->SetYTitle("dN_{ch}/d#eta ratio #frac{default composition}{modified composition}");
650
651 if (counter > 0)
652 hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);
653
654 file = TFile::Open("new_compositions_analysis.root", "UPDATE");
655 hRatios[counter]->Write();
656 file->Close();
657
658 delete fdNdEtaAnalysis;
659
660 counter++;
9f469bf5 661 }
10ebe68d 662
69b09e3b 663 /*
10ebe68d 664 gROOT->ProcessLine(".L drawPlots.C");
9f469bf5 665
6b7fa615 666 const char* fileNames[] = {"new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root" };
69b09e3b 667 const char* folderNames[] = { "PythiaRatios", "PiBoosted", "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced" };
9f469bf5 668
6b7fa615 669 Track2Particle1DComposition(fileNames, nCompositions, folderNames);
69b09e3b 670 */
9f469bf5 671}
672
9f469bf5 673
10ebe68d 674void drawSystematics()
675{
9f469bf5 676 //Secondaries();
677 //DrawDifferentSpecies();
678 //Composition();
679
680 Sigma2VertexSimulation();
27c04e01 681
10ebe68d 682}
7af955da 683
684void DrawdNdEtaDifferences()
685{
686 TH1* hists[5];
687
72e597d7 688 TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
f1076daa 689 legend->SetFillColor(0);
690
7af955da 691 TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
692 canvas->Divide(2, 1);
693
694 canvas->cd(1);
695
696 for (Int_t i=0; i<5; ++i)
697 {
72e597d7 698 hists[i] = 0;
7af955da 699 TFile* file = 0;
f1076daa 700 TString title;
7af955da 701
702 switch(i)
703 {
f1076daa 704 case 0 : file = TFile::Open("systematics_dndeta_reference.root"); title = "standard composition"; break;
705 case 1 : file = TFile::Open("systematics_dndeta_KBoosted.root"); title = "+ 50% kaons"; break;
706 case 2 : file = TFile::Open("systematics_dndeta_KReduced.root"); title = "- 50% kaons"; break;
707 case 3 : file = TFile::Open("systematics_dndeta_pBoosted.root"); title = "+ 50% protons"; break;
708 case 4 : file = TFile::Open("systematics_dndeta_pReduced.root"); title = "- 50% protons"; break;
7af955da 709 default: return;
710 }
711
72e597d7 712 if (file)
713 {
714 hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
715 hists[i]->SetTitle("a)");
716
717 Prepare1DPlot(hists[i], kFALSE);
718 hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
da7acf97 719 hists[i]->GetYaxis()->SetRangeUser(6, 7);
8ca1a6d9 720 hists[i]->SetLineColor(colors[i]);
721 hists[i]->SetMarkerColor(colors[i]);
722 hists[i]->SetMarkerStyle(markers[i]);
72e597d7 723 hists[i]->GetXaxis()->SetLabelOffset(0.015);
724 hists[i]->GetYaxis()->SetTitleOffset(1.5);
725 gPad->SetLeftMargin(0.12);
726 hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
727
728 legend->AddEntry(hists[i], title);
729 hists[i]->SetTitle(title);
730 }
7af955da 731 }
f1076daa 732 legend->Draw();
7af955da 733
734 canvas->cd(2);
f1076daa 735 gPad->SetLeftMargin(0.14);
736
737 TLegend* legend2 = new TLegend(0.73, 0.73, 0.98, 0.98);
738 legend2->SetFillColor(0);
7af955da 739
740 for (Int_t i=1; i<5; ++i)
741 {
72e597d7 742 if (hists[i])
743 {
744 legend2->AddEntry(hists[i]);
745
746 hists[i]->Divide(hists[0]);
747 hists[i]->SetTitle("b)");
8ca1a6d9 748 hists[i]->SetLineColor(colors[i-1]);
749 hists[i]->SetMarkerColor(colors[i-1]);
da7acf97 750 hists[i]->GetYaxis()->SetRangeUser(0.95, 1.05);
72e597d7 751 hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
752 hists[i]->GetYaxis()->SetTitleOffset(1.8);
753 hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
754 }
7af955da 755 }
f1076daa 756
757 legend2->Draw();
758
8ca1a6d9 759 canvas->SaveAs("particlecomposition_result_detail.gif");
760
761 TCanvas* canvas2 = new TCanvas("DrawdNdEtaDifferences2", "DrawdNdEtaDifferences2", 700, 500);
762
763 for (Int_t i=1; i<5; ++i)
764 {
765 if (hists[i])
766 {
767 hists[i]->SetTitle("");
768 hists[i]->GetYaxis()->SetTitleOffset(1.1);
769 hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
770 }
771 }
772
773 legend2->Draw();
774
775 canvas2->SaveAs("particlecomposition_result.gif");
776 canvas2->SaveAs("particlecomposition_result.eps");
7af955da 777}
27c04e01 778
51f6de65 779void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char_t* correctionFileName="correction_mapprocess-types.root", const char* analysisFileName = "analysis_esd_raw.root", const Char_t* outputFileName="systematics_vtxtrigger_compositions.root") {
27c04e01 780 //
781 // Function used to merge standard corrections with vertex
782 // reconstruction corrections obtained by a certain mix of ND, DD
783 // and SD events.
dd367a14 784 //
785 // the dn/deta spectrum is corrected and the ratios
786 // (standard to changed x-section) of the different dN/deta
787 // distributions are saved to a file.
788 //
51f6de65 789 // correctionTarget is of type AlidNdEtaCorrection::CorrectionType
790 // kINEL = 3
791 // kNSD = 4
27c04e01 792
dd367a14 793 loadlibs();
27c04e01 794
da7acf97 795 const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
eaa3702a 796
1d107532 797 //Karel:
798// fsd = 0.153 +- 0.031 (0.050 to take into account SD definition) --> change
799// fdd = 0.080 +- 0.050 --> change
800// fnd = 0.767 +- 0.059 --> keep (error small)
51f6de65 801
1d107532 802// const Char_t* changes[] = { "pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdlessddmore25", "sdmoreddless25"};
803 //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 };
a7f69e56 804 //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};
805 //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};
1d107532 806 //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};
807 //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};
808 Int_t nChanges = 9;
809
810 const Char_t* changes[] = { "ua5","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless" };
811 Float_t scalesND[] = {0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767};
812 Float_t scalesDD[] = {0.080, 0.130, 0.030, 0.080, 0.080, 0.130, 0.030, 0.130, 0.030};
813 Float_t scalesSD[] = {0.153, 0.153, 0.153, 0.203, 0.103, 0.203, 0.103, 0.103, 0.203};
814
815 for (Int_t i=0; i<9; i++)
816 scalesND[i] = 1.0 - scalesDD[i] - scalesSD[i];
eaa3702a 817
51f6de65 818 /*
eaa3702a 819 const Char_t* changes[] = { "pythia", "qgsm", "phojet"};
820 Float_t scalesND[] = {1.0, 1.10, 1.11};
821 Float_t scalesSD[] = {1.0, 0.69, 0.86};
822 Float_t scalesDD[] = {1.0, 0.98, 0.61};
823 Int_t nChanges = 3;
51f6de65 824 */
825
27c04e01 826 // cross section from Pythia
51f6de65 827 // 14 TeV!
27c04e01 828 Float_t sigmaND = 55.2;
829 Float_t sigmaDD = 9.78;
830 Float_t sigmaSD = 14.30;
831
dd367a14 832 // standard correction
833 TFile::Open(correctionFileName);
834 AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
835 correctionStandard->LoadHistograms();
836
837 // dont take vertexreco from this one
838 correctionStandard->GetVertexRecoCorrection()->Reset();
839 // dont take triggerbias from this one
840 correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
841 correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
842 correctionStandard->GetTriggerBiasCorrectionND()->Reset();
843
844 AlidNdEtaCorrection* corrections[100];
845 TH1F* hRatios[100];
846
847 Int_t counter = 0;
1d107532 848 for (Int_t j=2; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)
27c04e01 849
eaa3702a 850 for (Int_t i=0; i<nChanges; i++) {
dd367a14 851 TFile::Open(correctionFileName);
27c04e01 852
da7acf97 853 TString name;
854 name.Form("dndeta_correction_syst_%s_%s", typeName[j], changes[i]);
855 AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
770a1f1d 856 current->LoadHistograms("dndeta_correction");
857 current->Reset();
da7acf97 858
dd367a14 859 name.Form("dndeta_correction_ND");
860 AlidNdEtaCorrection* dNdEtaCorrectionND = new AlidNdEtaCorrection(name,name);
861 dNdEtaCorrectionND->LoadHistograms();
862 name.Form("dndeta_correction_DD");
863 AlidNdEtaCorrection* dNdEtaCorrectionDD = new AlidNdEtaCorrection(name,name);
864 dNdEtaCorrectionDD->LoadHistograms();
865 name.Form("dndeta_correction_SD");
866 AlidNdEtaCorrection* dNdEtaCorrectionSD = new AlidNdEtaCorrection(name,name);
867 dNdEtaCorrectionSD->LoadHistograms();
da7acf97 868
869 // calculating relative
1d107532 870 Float_t nd = dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
871 Float_t dd = dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
872 Float_t sd = dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
873 Float_t total = nd + dd + sd;
874
875 nd /= total;
876 sd /= total;
877 dd /= total;
878
879 Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
880
881 Float_t scaleND = scalesND[i] / nd;
882 Float_t scaleDD = scalesDD[i] / dd;
883 Float_t scaleSD = scalesSD[i] / sd;
884
885 Printf("ND=%.2f, DD=%.2f, SD=%.2f",scaleND, scaleDD, scaleSD);
886
887/* Float_t nd = 100 * sigmaND/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
da7acf97 888 Float_t dd = 100 * (scalesDD[i]*sigmaDD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
889 Float_t sd = 100 * (scalesSD[i]*sigmaSD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
890
1d107532 891 printf(Form("%s : ND=%.2f\%, DD=%.2f\%, SD=%.2f\% \n",changes[i],nd,dd,sd));*/
892 current->SetTitle(Form("ND=%.2f\%,DD=%.2f\%,SD=%.2f\%",scaleND,scaleDD,scaleSD));
da7acf97 893 current->SetTitle(name);
894
895 // scale
896 if (j == 0 || j == 2)
897 {
1d107532 898 dNdEtaCorrectionND->GetVertexRecoCorrection()->Scale(scaleND);
899 dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scaleDD);
900 dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scaleSD);
da7acf97 901 }
902 if (j == 1 || j == 2)
903 {
1d107532 904 dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->Scale(scaleND);
dd367a14 905 dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scalesDD[i]);
906 dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scalesSD[i]);
907
1d107532 908 dNdEtaCorrectionND->GetTriggerBiasCorrectionNSD()->Scale(scaleND);
909 dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scaleDD);
910 dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scaleSD);
dd367a14 911
1d107532 912 dNdEtaCorrectionND->GetTriggerBiasCorrectionND()->Scale(scaleND);
913 dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scaleDD);
914 dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scaleSD);
da7acf97 915 }
916
dd367a14 917 //clear track in correction
918 dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
919 dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
920 dNdEtaCorrectionSD->GetTrack2ParticleCorrection()->Reset();
da7acf97 921
922 TList collection;
923 collection.Add(correctionStandard);
dd367a14 924 collection.Add(dNdEtaCorrectionND);
925 collection.Add(dNdEtaCorrectionDD);
926 collection.Add(dNdEtaCorrectionSD);
da7acf97 927
928 current->Merge(&collection);
929 current->Finish();
930
dd367a14 931 corrections[counter] = current;
932
933 // now correct dNdeta distribution with modified correction map
934 TFile::Open(analysisFileName);
935
936 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
937 fdNdEtaAnalysis->LoadHistograms();
938
51f6de65 939 fdNdEtaAnalysis->Finish(current, 0.3, correctionTarget, Form("%d %d", j, i));
dd367a14 940
941 name = "ratio";
942 if (j==0) name.Append("_vetexReco_");
943 if (j==1) name.Append("_triggerBias_");
944 if (j==2) name.Append("_vertexReco_triggerBias_");
945 name.Append(changes[i]);
946
947 hRatios[counter] = (TH1F*) fdNdEtaAnalysis->GetdNdEtaHistogram()->Clone(name);
948
1d107532 949 name.Form("ND #times %0.2f DD #times %0.2f, SD #times %0.2f",scaleND,scaleDD,scaleSD);
dd367a14 950 hRatios[counter]->SetTitle(name.Data());
69b09e3b 951 hRatios[counter]->SetYTitle("dN_{ch}/d#eta ratio #frac{default cross-section}{modified cross-sections}");
dd367a14 952
953 if (counter > 0)
954 hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);
955
956 delete fdNdEtaAnalysis;
957
958 counter++;
da7acf97 959 }
27c04e01 960 }
da7acf97 961
27c04e01 962 TFile* fout = new TFile(outputFileName,"RECREATE");
963
dd367a14 964 // to make everything consistent
965 hRatios[0]->Divide(hRatios[0],hRatios[0],1,1);
966
1d107532 967 for (Int_t i=0; i<counter; i++)
dd367a14 968 {
27c04e01 969 corrections[i]->SaveHistograms();
dd367a14 970 hRatios[i]->Write();
971 }
27c04e01 972
973 fout->Write();
974 fout->Close();
975}
5e08578b 976
12bb57f1 977void CreateCorrectionsWithUA5CrossSections(Int_t origin, const Char_t* correctionFileName="correction_mapprocess-types.root", const Char_t* outputFileName="correction_map2.root") {
1d107532 978 //
979 // Function used to merge standard corrections with vertex
980 // reconstruction corrections obtained by a certain mix of ND, DD
981 // and SD events.
982 //
12bb57f1 983 // origin:
984 // -1 = Pythia (test)
985 // 0 = UA5
986 // 1 = Data 1.8 TeV
987 // 2 = Tel-Aviv
988 // 3 = Durham
989 //
1d107532 990
991 loadlibs();
992
993 const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
994
12bb57f1 995 Float_t ref_SD = -1;
996 Float_t ref_DD = -1;
997 Float_t ref_ND = -1;
998
999 switch (origin)
1000 {
1001 case -1: // Pythia, as test
1002 ref_SD = 0.223788;
1003 ref_DD = 0.123315;
1004 ref_ND = 0.652897;
1005 break;
1006
1007 case 0: // UA5
1008 ref_SD = 0.153;
1009 ref_DD = 0.080;
1010 ref_ND = 0.767;
1011 break;
1012
1013 case 1: // data 1.8 TeV
1014 ref_SD = 0.152;
1015 ref_DD = 0.092;
1016 ref_ND = 1 - ref_SD - ref_DD;
1017 break;
1018
1019 case 2: // tel-aviv model
1020 ref_SD = 0.171;
1021 ref_DD = 0.094;
1022 ref_ND = 1 - ref_SD - ref_DD;
1023 break;
1024
1025 case 3: // durham model
1026 ref_SD = 0.190;
1027 ref_DD = 0.125;
1028 ref_ND = 1 - ref_SD - ref_DD;
1029 break;
1030
1031 default:
1032 return;
1033 }
1034
1035 //Karel (UA5):
1d107532 1036// fsd = 0.153 +- 0.031
1037// fdd = 0.080 +- 0.050
1038// fnd = 0.767 +- 0.059
1039
12bb57f1 1040// Karel (1.8 TeV):
1041//
1042// Tel-Aviv model Sd/Inel = 0.171 Dd/Inel = 0.094
1043// Durham model Sd/Inel = 0.190 Dd/Inel = 0.125
1044// Data Sd/Inel = 0.152 +- 0.030 Dd/Inel = 0.092 +- 0.45
1045
1046
1d107532 1047
1048 // standard correction
1049 TFile::Open(correctionFileName);
1050 AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1051 correctionStandard->LoadHistograms();
1052
1053 // dont take vertexreco from this one
1054 correctionStandard->GetVertexRecoCorrection()->Reset();
1055 // dont take triggerbias from this one
1056 correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
1057 correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
1058 correctionStandard->GetTriggerBiasCorrectionND()->Reset();
1059
1060 AlidNdEtaCorrection* corrections[100];
1061 TH1F* hRatios[100];
1062
1063 Int_t counter = 0;
1064
1065 TFile::Open(correctionFileName);
1066
1067 AlidNdEtaCorrection* current = new AlidNdEtaCorrection("dndeta_correction_ua5", "dndeta_correction_ua5");
1068 current->LoadHistograms("dndeta_correction");
1069 current->Reset();
1070
1071 TString name;
1072 name.Form("dndeta_correction_ND");
1073 AlidNdEtaCorrection* dNdEtaCorrectionND = new AlidNdEtaCorrection(name,name);
1074 dNdEtaCorrectionND->LoadHistograms();
1075 name.Form("dndeta_correction_DD");
1076 AlidNdEtaCorrection* dNdEtaCorrectionDD = new AlidNdEtaCorrection(name,name);
1077 dNdEtaCorrectionDD->LoadHistograms();
1078 name.Form("dndeta_correction_SD");
1079 AlidNdEtaCorrection* dNdEtaCorrectionSD = new AlidNdEtaCorrection(name,name);
1080 dNdEtaCorrectionSD->LoadHistograms();
1081
1082 // calculating relative
1083 Float_t nd = dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1084 Float_t dd = dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1085 Float_t sd = dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1086 Float_t total = nd + dd + sd;
1087
1088 nd /= total;
1089 sd /= total;
1090 dd /= total;
1091
1092 Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
1093
12bb57f1 1094 Float_t scaleND = ref_ND / nd;
1095 Float_t scaleDD = ref_DD / dd;
1096 Float_t scaleSD = ref_SD / sd;
1d107532 1097
1098 Printf("ND=%.2f, DD=%.2f, SD=%.2f",scaleND, scaleDD, scaleSD);
1099
1100 // scale
1101 dNdEtaCorrectionND->GetVertexRecoCorrection()->Scale(scaleND);
1102 dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scaleDD);
1103 dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scaleSD);
1104
1105 dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->Scale(scaleND);
1106 dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scaleDD);
1107 dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scaleSD);
1108
1109 dNdEtaCorrectionND->GetTriggerBiasCorrectionNSD()->Scale(scaleND);
1110 dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scaleDD);
1111 dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scaleSD);
1112
1113 dNdEtaCorrectionND->GetTriggerBiasCorrectionND()->Scale(scaleND);
1114 dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scaleDD);
1115 dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scaleSD);
1116
1117 //clear track in correction
1118 dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
1119 dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
1120 dNdEtaCorrectionSD->GetTrack2ParticleCorrection()->Reset();
1121
1122 TList collection;
1123 collection.Add(correctionStandard);
1124 collection.Add(dNdEtaCorrectionND);
1125 collection.Add(dNdEtaCorrectionDD);
1126 collection.Add(dNdEtaCorrectionSD);
1127
1128 current->Merge(&collection);
1129 current->Finish();
1130
1131 TFile* fout = new TFile(outputFileName,"RECREATE");
1132 current->SaveHistograms();
1133
1134 fout->Write();
1135 fout->Close();
1136
1137 Printf("Trigger efficiencies:");
1138 Printf("ND: %.2f %%", 100.0 * dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1139 Printf("SD: %.2f %%", 100.0 * dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1140 Printf("DD: %.2f %%", 100.0 * dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1141 Printf("INEL: %.2f %%", 100.0 * current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1142 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()));
1143}
1144
5e08578b 1145DrawTriggerEfficiency(Char_t* fileName) {
1146
1147 gStyle->SetOptStat(0);
1148 gStyle->SetOptTitle(0);
1149 gStyle->SetOptFit(0);
1150
1151 gStyle->SetTextSize(0.04);
1152 gStyle->SetTitleSize(0.05,"xyz");
1153 //gStyle->SetTitleFont(133, "xyz");
1154 //gStyle->SetLabelFont(133, "xyz");
1155 //gStyle->SetLabelSize(17, "xyz");
1156 gStyle->SetLabelOffset(0.01, "xyz");
1157
1158 gStyle->SetTitleOffset(1.1, "y");
1159 gStyle->SetTitleOffset(1.1, "x");
1160 gStyle->SetEndErrorSize(0.0);
1161
1162 //##############################################
1163
1164 //making canvas and pads
1165 TCanvas *c = new TCanvas("trigger_eff", "",600,500);
1166
1167 TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
1168
1169 p1->SetBottomMargin(0.15);
1170 p1->SetTopMargin(0.03);
1171 p1->SetLeftMargin(0.15);
1172 p1->SetRightMargin(0.03);
1173
1174 p1->SetGridx();
1175 p1->SetGridy();
1176
1177 p1->Draw();
1178 p1->cd();
1179
1180 gSystem->Load("libPWG0base");
1181
1182 AlidNdEtaCorrection* corrections[4];
1183 AliCorrectionMatrix2D* triggerBiasCorrection[4];
1184
1185 TH1F* hTriggerEffInv[4];
1186 TH1F* hTriggerEff[4];
1187
1188 Char_t* names[] = {"triggerBiasND", "triggerBiasDD", "triggerBiasSD", "triggerBiasINEL"};
1189
1190 for (Int_t i=0; i<4; i++) {
1191 corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1192 corrections[i]->LoadHistograms(fileName, names[i]);
1193
1194 triggerBiasCorrection[i] = corrections[i]->GetTriggerBiasCorrectionINEL();
1195
1196
1197 hTriggerEffInv[i] = triggerBiasCorrection[i]->Get1DCorrection();
1198 hTriggerEff[i] = (TH1F*)hTriggerEffInv[i]->Clone();
1199
1200 for (Int_t b=0; b<=hTriggerEff[i]->GetNbinsX(); b++) {
1201 hTriggerEff[i]->SetBinContent(b,1);
1202 hTriggerEff[i]->SetBinError(b,0);
1203 }
1204 hTriggerEff[i]->Divide(hTriggerEff[i],hTriggerEffInv[i]);
1205 hTriggerEff[i]->Scale(100);
1206 }
1207
1208 Int_t colors[] = {2,3,4,1};
1209 Float_t effs[4];
1210 for (Int_t i=0; i<4; i++) {
1211 hTriggerEff[i]->Fit("pol0","","",-20,20);
1212 effs[i] = ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->GetParameter(0);
1213 ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineWidth(1);
1214 ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineStyle(2);
1215 ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineColor(colors[i]);
1216 cout << effs[i] << endl;
1217 }
1218
1219
1220 Char_t* text[] = {"ND", "DD", "SD", "INEL"};
1221 TLatex* latex[4];
1222
1223 TH2F* null = new TH2F("","",100,-25,35,100,0,110);
1224 null->SetXTitle("Vertex z [cm]");
1225 null->GetXaxis()->CenterTitle(kTRUE);
1226 null->SetYTitle("Trigger efficiency [%]");
1227 null->GetYaxis()->CenterTitle(kTRUE);
1228 null->Draw();
1229
1230
1231 for (Int_t i=0; i<4; i++) {
1232 hTriggerEff[i]->SetLineWidth(2);
1233 hTriggerEff[i]->SetLineColor(colors[i]);
1234
1235 hTriggerEff[i]->Draw("same");
1236
1237 latex[i] = new TLatex(22,effs[i]-1.5, Form("%s (%0.1f)",text[i],effs[i]));
1238 latex[i]->SetTextColor(colors[i]);
1239 latex[i]->Draw();
1240 }
1241
1242}
1243
1244
1245DrawSpectraPID(Char_t* fileName) {
1246
1247 gSystem->Load("libPWG0base");
1248
1249 Char_t* names[] = {"correction_0", "correction_1", "correction_2", "correction_3"};
1250 AlidNdEtaCorrection* corrections[4];
1251 AliCorrectionMatrix3D* trackToPartCorrection[4];
1252
1253 TH1F* measuredPt[4];
1254 TH1F* generatedPt[4];
1255 TH1F* ratioPt[4];
1256
1257 for (Int_t i=0; i<4; i++) {
1258 corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1259 corrections[i]->LoadHistograms(fileName, names[i]);
1260
1261 trackToPartCorrection[i] = corrections[i]->GetTrack2ParticleCorrection();
1262
1263 Int_t binX1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(-10);
1264 Int_t binX2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(10);
1265 Int_t binY1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(-1);
1266 Int_t binY2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(1);
1267
1268 measuredPt[i] = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->ProjectionZ(Form("m_%d",i),binX1,binX2,binY1,binY2);
1269 generatedPt[i] = (TH1F*)trackToPartCorrection[i]->GetGeneratedHistogram()->ProjectionZ(Form("g_%d",i),binX1,binX2,binY1,binY2);
1270 ratioPt[i] = (TH1F*)generatedPt[i]->Clone(Form("r_%d",i));
1271 ratioPt[i]->Divide(measuredPt[i], generatedPt[i], 1,1,"B");
1272 }
1273
1274 ratioPt[0]->Draw();
1275
1276 for (Int_t i=0; i<3; i++) {
1277 ratioPt[i]->SetLineColor(i+1);
1278 ratioPt[i]->SetLineWidth(2);
1279
1280 ratioPt[i]->Draw("same");
1281
1282 }
1283
1284 return;
1285 measuredPt[0]->SetLineColor(2);
1286 measuredPt[0]->SetLineWidth(5);
1287
1288 measuredPt[0]->Draw();
1289 generatedPt[0]->Draw("same");
1290}
da7acf97 1291
69b09e3b 1292void changePtSpectrum(const char* fileName = "analysis_mc.root", Float_t ptCutOff = 0.2, const char* fileName2 = 0)
da7acf97 1293{
69b09e3b 1294 Float_t factor = 0.5;
1295
8277513e 1296 TFile* file = TFile::Open(fileName);
69b09e3b 1297 TH1F* hist = dynamic_cast<TH1F*> (file->Get("dndeta_check_pt")->Clone());
1298
1299 TH1* hist2 = 0;
1300 if (fileName2)
1301 {
1302 file2 = TFile::Open(fileName2);
1303 hist2 = dynamic_cast<TH1*> (file2->Get("dndeta_check_pt")->Clone());
1304 hist2->Scale(hist->GetBinContent(hist->FindBin(ptCutOff)) / hist2->GetBinContent(hist2->FindBin(ptCutOff)));
1305 }
1306
1307 //hist->Scale(1.0 / hist->Integral());
da7acf97 1308
8ca1a6d9 1309 //hist->Rebin(3);
1310 //hist->Scale(1.0/3);
da7acf97 1311
1312 TH1F* clone1 = dynamic_cast<TH1F*> (hist->Clone("clone1"));
1313 TH1F* clone2 = dynamic_cast<TH1F*> (hist->Clone("clone2"));
1314
1315 TH1F* scale1 = dynamic_cast<TH1F*> (hist->Clone("scale1"));
1316 TH1F* scale2 = dynamic_cast<TH1F*> (hist->Clone("scale2"));
1317
da7acf97 1318 for (Int_t i=1; i <= hist->GetNbinsX(); ++i)
1319 {
1320 if (hist->GetBinCenter(i) > ptCutOff)
1321 {
1322 scale1->SetBinContent(i, 1);
1323 scale2->SetBinContent(i, 1);
1324 }
1325 else
1326 {
1327 // 90 % at pt = 0, 0% at pt = ptcutoff
69b09e3b 1328 scale1->SetBinContent(i, 1 - (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * factor);
da7acf97 1329
1330 // 110% at pt = 0, ...
69b09e3b 1331 scale2->SetBinContent(i, 1 + (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * factor);
da7acf97 1332 }
8ca1a6d9 1333 scale1->SetBinError(i, 0);
1334 scale2->SetBinError(i, 0);
da7acf97 1335 }
1336
69b09e3b 1337 /*
da7acf97 1338 new TCanvas;
1339 scale1->Draw();
1340 scale2->SetLineColor(kRed);
1341 scale2->Draw("SAME");
69b09e3b 1342 */
da7acf97 1343
1344 clone1->Multiply(scale1);
1345 clone2->Multiply(scale2);
1346
1347 Prepare1DPlot(hist);
1348 Prepare1DPlot(clone1);
1349 Prepare1DPlot(clone2);
1350
8ca1a6d9 1351 /*hist->SetMarkerStyle(markers[0]);
1352 clone1->SetMarkerStyle(markers[0]);
1353 clone2->SetMarkerStyle(markers[0]);*/
1354
8277513e 1355 hist->SetTitle(";p_{T} in GeV/c;dN_{ch}/dp_{T} in c/GeV");
69b09e3b 1356 hist->GetXaxis()->SetRangeUser(0, 0.5);
8ca1a6d9 1357 hist->GetYaxis()->SetRangeUser(0.01, clone2->GetMaximum() * 1.1);
da7acf97 1358 hist->GetYaxis()->SetTitleOffset(1);
1359
69b09e3b 1360 TCanvas* canvas = new TCanvas("c", "c", 600, 600);
1361 gPad->SetGridx();
1362 gPad->SetGridy();
1363 gPad->SetTopMargin(0.05);
1364 gPad->SetRightMargin(0.05);
8ca1a6d9 1365 gPad->SetBottomMargin(0.12);
1366 hist->Draw("H");
da7acf97 1367 clone1->SetLineColor(kRed);
8ca1a6d9 1368 clone1->Draw("HSAME");
da7acf97 1369 clone2->SetLineColor(kBlue);
8ca1a6d9 1370 clone2->Draw("HSAME");
1371 hist->Draw("HSAME");
69b09e3b 1372
1373 if (hist2)
1374 {
1375 Prepare1DPlot(hist2);
1376 hist2->SetLineStyle(2);
1377 hist2->Draw("HSAME");
1378 }
da7acf97 1379
1380 Float_t fraction = hist->Integral(hist->GetXaxis()->FindBin(ptCutOff), hist->GetNbinsX()) / hist->Integral(1, hist->GetNbinsX());
1381 Float_t fraction1 = clone1->Integral(clone1->GetXaxis()->FindBin(ptCutOff), clone1->GetNbinsX()) / clone1->Integral(1, clone1->GetNbinsX());
1382 Float_t fraction2 = clone2->Integral(clone2->GetXaxis()->FindBin(ptCutOff), clone2->GetNbinsX()) / clone2->Integral(1, clone2->GetNbinsX());
1383
1384 printf("%f %f %f\n", fraction, fraction1, fraction2);
1385 printf("Rel. %f %f\n", fraction1 / fraction, fraction2 / fraction);
1386
69b09e3b 1387 //canvas->SaveAs("changePtSpectrum.gif");
8ca1a6d9 1388 canvas->SaveAs("changePtSpectrum.eps");
da7acf97 1389}
fb30f6e3 1390
1391void FractionBelowPt()
1392{
1393 gSystem->Load("libPWG0base");
1394
1395 AlidNdEtaCorrection* fdNdEtaCorrection[4];
1396
1397 TFile::Open("systematics.root");
1398
1399 for (Int_t i=0; i<4; ++i)
1400 {
1401 TString name;
1402 name.Form("correction_%d", i);
1403 fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
1404 fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
1405 }
1406
1407 Double_t geneCount[5];
1408 Double_t measCount[5];
1409 geneCount[4] = 0;
1410 measCount[4] = 0;
1411
1412 for (Int_t i=0; i<4; ++i)
1413 {
1414 TH3F* hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1415 geneCount[i] = hist->Integral(hist->GetXaxis()->FindBin(-10), hist->GetXaxis()->FindBin(10),
1416 hist->GetYaxis()->FindBin(-0.8), hist->GetYaxis()->FindBin(0.8),
1417 1, hist->GetZaxis()->FindBin(0.3));
1418
1419 hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1420 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));
1421
1422 geneCount[4] += geneCount[i];
1423 measCount[4] += measCount[i];
1424
1425 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]);
1426 }
1427
1428 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]);
1429
1430 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]);
7f0efb86 1431}
8ca1a6d9 1432
1433
1434mergeCorrectionsMisalignment(Char_t* alignedFile = "correction_map_aligned.root",
1435 Char_t* misalignedFile = "correction_map_misaligned.root",
1436 Char_t* outputFileName="correction_map_misaligned_single.root")
1437{
1438 //
1439 // from the aligned and misaligned corrections, 3 new corrections are created
1440 // in these new corrections only one of the corrections (track2particle, vertex, trigger)
1441 // is taken from the misaligned input to allow study of the effect on the different
1442 // corrections
1443
1444 gSystem->Load("libPWG0base");
1445
1446 const Char_t* typeName[] = { "track2particle", "vertex", "trigger" };
1447
1448 AlidNdEtaCorrection* corrections[3];
1449 for (Int_t j=0; j<3; j++) { // j = 0 (track2particle), j = 1 (vertex), j = 2 (trigger)
1450 AlidNdEtaCorrection* alignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1451 alignedCorrection->LoadHistograms(alignedFile);
1452
1453 AlidNdEtaCorrection* misalignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1454 misalignedCorrection->LoadHistograms(misalignedFile);
1455
1456 TString name;
1457 name.Form("dndeta_correction_alignment_%s", typeName[j]);
1458 AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
1459
1460 switch (j)
1461 {
1462 case 0:
1463 alignedCorrection->GetTrack2ParticleCorrection()->Reset();
1464 misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1465 misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1466 misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1467 misalignedCorrection->GetVertexRecoCorrection()->Reset();
1468 break;
1469
1470 case 1:
1471 misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1472 misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1473 misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1474 misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1475 alignedCorrection->GetVertexRecoCorrection()->Reset();
1476 break;
1477
1478 case 2:
1479 misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1480 alignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1481 alignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1482 alignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1483 misalignedCorrection->GetVertexRecoCorrection()->Reset();
1484 break;
1485
1486 default:
1487 return;
1488 }
1489
1490 TList collection;
1491 collection.Add(misalignedCorrection);
1492 collection.Add(alignedCorrection);
1493
1494 current->Merge(&collection);
1495 current->Finish();
1496
1497 corrections[j] = current;
1498 }
1499
1500 TFile* fout = new TFile(outputFileName, "RECREATE");
1501
1502 for (Int_t i=0; i<3; i++)
1503 corrections[i]->SaveHistograms();
1504
1505 fout->Write();
1506 fout->Close();
1507}
1508
1509
1510void DrawdNdEtaDifferencesAlignment()
1511{
1512 TH1* hists[5];
1513
1514 TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
1515 legend->SetFillColor(0);
1516
1517 TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
1518 canvas->Divide(2, 1);
1519
1520 canvas->cd(1);
1521
1522 for (Int_t i=0; i<5; ++i)
1523 {
1524 hists[i] = 0;
1525 TFile* file = 0;
1526 TString title;
1527
1528 switch(i)
1529 {
1530 case 0 : file = TFile::Open("systematics_misalignment_aligned.root"); title = "aligned"; break;
1531 case 1 : file = TFile::Open("systematics_misalignment_misaligned.root"); title = "fully misaligned"; break;
1532 case 2 : file = TFile::Open("systematics_misalignment_track2particle.root"); title = "only track2particle"; break;
1533 case 3 : file = TFile::Open("systematics_misalignment_vertex.root"); title = "only vertex rec."; break;
1534 case 4 : file = TFile::Open("systematics_misalignment_trigger.root"); title = "only trigger bias"; break;
1535 default: return;
1536 }
1537
1538 if (file)
1539 {
1540 hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
1541 hists[i]->SetTitle("");
1542
1543 Prepare1DPlot(hists[i], kFALSE);
1544 hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
1545 hists[i]->GetYaxis()->SetRangeUser(6, 7);
1546 hists[i]->SetLineWidth(1);
1547 hists[i]->SetLineColor(colors[i]);
1548 hists[i]->SetMarkerColor(colors[i]);
1549 hists[i]->SetMarkerStyle(markers[i]);
1550 hists[i]->GetXaxis()->SetLabelOffset(0.015);
1551 hists[i]->GetYaxis()->SetTitleOffset(1.5);
1552 gPad->SetLeftMargin(0.12);
1553 hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
1554
1555 legend->AddEntry(hists[i], title);
1556 hists[i]->SetTitle(title);
1557 }
1558 }
1559 legend->Draw();
1560
1561 canvas->cd(2);
1562 gPad->SetLeftMargin(0.14);
1563
1564 TLegend* legend2 = new TLegend(0.63, 0.73, 0.98, 0.98);
1565 legend2->SetFillColor(0);
1566
1567 for (Int_t i=1; i<5; ++i)
1568 {
1569 if (hists[i])
1570 {
1571 legend2->AddEntry(hists[i]);
1572
1573 hists[i]->Divide(hists[0]);
1574 hists[i]->SetTitle("b)");
1575 hists[i]->GetYaxis()->SetRangeUser(0.9, 1.1);
1576 hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
1577 hists[i]->GetYaxis()->SetTitleOffset(1.8);
1578 hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
1579 }
1580 }
1581
1582 legend2->Draw();
1583
1584 canvas->SaveAs("misalignment_result.eps");
1585 canvas->SaveAs("misalignment_result.gif");
1586}
1587
1588void EvaluateMultiplicityEffect()
1589{
1590 gSystem->Load("libPWG0base");
1591
1592 AlidNdEtaCorrection* fdNdEtaCorrectionLow[4];
1593 AlidNdEtaCorrection* fdNdEtaCorrectionHigh[4];
1594
1595 TFile::Open("systematics-low-multiplicity.root");
1596
1597 for (Int_t i=0; i<4; ++i)
1598 {
1599 TString name;
1600 name.Form("correction_%d", i);
1601 fdNdEtaCorrectionLow[i] = new AlidNdEtaCorrection(name, name);
1602 fdNdEtaCorrectionLow[i]->LoadHistograms("systematics-low-multiplicity.root", name);
1603 }
1604
1605 TList list;
1606 for (Int_t i=1; i<4; ++i)
1607 list.Add(fdNdEtaCorrectionLow[i]);
1608
1609 fdNdEtaCorrectionLow[0]->Merge(&list);
1610 fdNdEtaCorrectionLow[0]->Finish();
1611
1612 TFile::Open("systematics-high-multiplicity.root");
1613
1614 for (Int_t i=0; i<4; ++i)
1615 {
1616 TString name;
1617 name.Form("correction_%d", i);
1618 fdNdEtaCorrectionHigh[i] = new AlidNdEtaCorrection(name, name);
1619 fdNdEtaCorrectionHigh[i]->LoadHistograms("systematics-high-multiplicity.root", name);
1620 }
1621
1622 TList list2;
1623 for (Int_t i=1; i<4; ++i)
1624 list2.Add(fdNdEtaCorrectionHigh[i]);
1625
1626 fdNdEtaCorrectionHigh[0]->Merge(&list2);
1627 fdNdEtaCorrectionHigh[0]->Finish();
1628
1629 TH1F* outputLow = new TH1F("Track2ParticleLow", "Track2Particle at low multiplicity", 200, 0, 2);
1630 TH1F* outputHigh = new TH1F("Track2ParticleHigh", "Track2Particle at high multiplicity", 200, 0, 2);
1631
1632 TH3F* hist = fdNdEtaCorrectionLow[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1633 TH3F* hist2 = fdNdEtaCorrectionHigh[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1634 for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
1635 for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
1636 for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
1637 //for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
1638 {
1639 if (hist->GetBinContent(x, y, z) > 0)
1640 outputLow->Fill(hist->GetBinContent(x, y, z));
1641 //if (hist->GetBinContent(x, y, z) == 1)
1642 // 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));
1643
1644 if (hist2->GetBinContent(x, y, z) > 0)
1645 outputHigh->Fill(hist2->GetBinContent(x, y, z));
1646 }
1647
1648 TCanvas* canvas = new TCanvas("EvaluateMultiplicityEffect", "EvaluateMultiplicityEffect", 1000, 500);
1649 canvas->Divide(2, 1);
1650
1651 canvas->cd(1);
1652 outputLow->Draw();
1653 outputLow->Fit("gaus", "0");
1654 outputLow->GetFunction("gaus")->SetLineColor(2);
1655 outputLow->GetFunction("gaus")->DrawCopy("SAME");
1656
1657 canvas->cd(2);
1658 outputHigh->Draw();
1659 outputHigh->Fit("gaus", "0");
1660 outputHigh->GetFunction("gaus")->DrawCopy("SAME");
1661
1662 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1663}
1664
1665void 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")
1666{
1667 gSystem->Load("libPWG0base");
1668
1669 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1670 dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
1671
1672 dNdEtaCorrection->GetTrack2ParticleCorrection()->PlotBinErrors(xmin, xmax, ymin, ymax, zmin, zmax)->Draw();
1673}
9e952c39 1674
1675
1676void 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")
1677{
1678 gSystem->Load("libPWG0base");
1679
1680 TFile* file = TFile::Open(output, "RECREATE");
1681
1682 const Int_t max = 5;
1683 dNdEtaAnalysis* fdNdEtaAnalysis[5];
1684
1685 new TCanvas;
1686 TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
1687 legend->SetFillColor(0);
1688
1689 for (Int_t i = 0; i < max; ++i)
1690 {
1691 TFile::Open(baseCorrectionMapFile);
1692 AlidNdEtaCorrection* baseCorrection = new AlidNdEtaCorrection(baseCorrectionMapFolder, baseCorrectionMapFolder);
1693 baseCorrection->LoadHistograms();
1694
1695 AlidNdEtaCorrection::CorrectionType correctionType = AlidNdEtaCorrection::kNone;
1696 const char* name = 0;
1697
1698 TFile::Open(changedCorrectionMapFile);
1699 switch (i)
1700 {
1701 case 0 :
1702 name = "default";
1703 break;
1704
1705 case 1 :
1706 baseCorrection->GetTrack2ParticleCorrection()->LoadHistograms(Form("%s/Track2Particle", changedCorrectionMapFolder));
1707 name = "Track2Particle";
1708 break;
1709
1710 case 2 :
1711 baseCorrection->GetVertexRecoCorrection()->LoadHistograms(Form("%s/VertexReconstruction", changedCorrectionMapFolder));
1712 name = "VertexReco";
1713 break;
1714
1715 case 3 :
1716 baseCorrection->GetTriggerBiasCorrectionINEL()->LoadHistograms(Form("%s/TriggerBias_MBToINEL", changedCorrectionMapFolder));
1717 name = "TriggerBias_MBToINEL";
1718 break;
1719
1720 case 4 :
1721 baseCorrection->LoadHistograms(changedCorrectionMapFolder);
1722 name = "all";
1723 break;
1724
1725 default: return;
1726 }
1727
1728 TFile::Open(dataFile);
1729 fdNdEtaAnalysis[i] = new dNdEtaAnalysis(name, name);
1730 fdNdEtaAnalysis[i]->LoadHistograms("dndeta");
1731
1732 fdNdEtaAnalysis[i]->Finish(baseCorrection, 0.3, AlidNdEtaCorrection::kINEL);
1733 file->cd();
1734 fdNdEtaAnalysis[i]->SaveHistograms();
1735
1736 TH1* hist = fdNdEtaAnalysis[i]->GetdNdEtaHistogram(0);
1737 hist->SetLineColor(colors[i]);
1738 hist->SetMarkerColor(colors[i]);
1739 hist->SetMarkerStyle(markers[i]+1);
1740 hist->DrawCopy((i == 0) ? "" : "SAME");
1741 legend->AddEntry(hist, name);
1742 }
1743
1744 legend->Draw();
1745}
a7f69e56 1746
1747void ChangePtInCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
1748{
1749 loadlibs();
1750 if (!TFile::Open(fileName))
1751 return;
1752
1753 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
1754 if (!dNdEtaCorrection->LoadHistograms())
1755 return;
1756
1757 dNdEtaCorrection->Finish();
1758
1759 AliCorrection* corr = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kTrack2Particle);
1760
1761 Printf(">>>> Before");
1762 corr->PrintInfo(0);
1763
1764 Float_t factor = 0.5;
1765 Float_t ptCutOff = 0.2;
1766
1767 TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram();
1768 TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram();
1769
1770 for (Int_t z = 1; z <= gene->GetZaxis()->FindBin(ptCutOff - 0.01); z++)
1771 {
1772 Float_t localFactor = 1 - (ptCutOff - gene->GetZaxis()->GetBinCenter(z)) / ptCutOff * factor;
1773 Printf("%f %f", gene->GetZaxis()->GetBinCenter(z), localFactor);
1774 for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1775 for (Int_t y = 1; y <= gene->GetNbinsY(); y++)
1776 {
1777 gene->SetBinContent(x, y, z, gene->GetBinContent(x, y, z) * localFactor);
1778 meas->SetBinContent(x, y, z, meas->GetBinContent(x, y, z) * localFactor);
1779 }
1780 }
1781
1782 dNdEtaCorrection->Finish();
1783
1784 Printf(">>>> After");
1785 corr->PrintInfo(0);
1786}