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