3 #if !defined(__CINT__) || defined(__MAKECINT__)
5 #include "AliPWG0Helper.h"
6 #include "dNdEtaAnalysis.h"
7 #include "AlidNdEtaCorrection.h"
22 void Track2Particle1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 10);
26 Int_t markers[] = {20,20,21,22,23,28,29};
27 Int_t colors[] = {1,2,3,4,6,8,102};
31 gSystem->Load("libTree");
32 gSystem->Load("libVMC");
34 gSystem->Load("libSTEERBase");
35 gSystem->Load("libANALYSIS");
36 gSystem->Load("libPWG0base");
39 void SetRanges(TAxis* axis)
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);
51 void SetRanges(TH1* hist)
53 SetRanges(hist->GetXaxis());
54 SetRanges(hist->GetYaxis());
55 SetRanges(hist->GetZaxis());
58 void Prepare3DPlot(TH3* hist)
60 hist->GetXaxis()->SetTitleOffset(1.5);
61 hist->GetYaxis()->SetTitleOffset(1.5);
62 hist->GetZaxis()->SetTitleOffset(1.5);
64 hist->SetStats(kFALSE);
67 void Prepare2DPlot(TH2* hist)
69 hist->SetStats(kFALSE);
70 hist->GetYaxis()->SetTitleOffset(1.4);
75 void Prepare1DPlot(TH1* hist, Bool_t setRanges = kTRUE)
77 hist->SetLineWidth(2);
78 hist->SetStats(kFALSE);
80 hist->GetXaxis()->SetTitleOffset(1.2);
81 hist->GetYaxis()->SetTitleOffset(1.2);
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);
104 gPad->Range(0, 0, 1, 1);
105 gPad->SetRightMargin(0.15);
106 gPad->SetLeftMargin(0.12);
114 TFile* file = TFile::Open("systematics.root");
116 TH2F* secondaries = dynamic_cast<TH2F*> (file->Get("fSecondaries"));
119 printf("Could not read histogram\n");
123 TCanvas* canvas = new TCanvas("Secondaries", "Secondaries", 1000, 1000);
124 canvas->Divide(3, 3);
125 for (Int_t i=1; i<=8; i++)
127 TH1D* hist = secondaries->ProjectionY(Form("proj_%d", i), i, i);
128 hist->SetTitle(secondaries->GetXaxis()->GetBinLabel(i));
135 void Track2Particle1DComposition(const char** fileNames, Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9)
137 gSystem->Load("libPWG0base");
140 canvasName.Form("Track2Particle1DComposition");
141 TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
142 canvas->Divide(3, 1);
144 TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);
146 for (Int_t i=0; i<folderCount; ++i)
148 Correction1DCreatePlots(fileNames[i], folderNames[i], upperPtLimit);
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"));
154 Prepare1DPlot(corrX);
155 Prepare1DPlot(corrY);
156 Prepare1DPlot(corrZ);
158 const char* title = "Track2Particle Correction";
159 corrX->SetTitle(title);
160 corrY->SetTitle(title);
161 corrZ->SetTitle(title);
163 corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
165 corrX->SetLineColor(i+1);
166 corrY->SetLineColor(i+1);
167 corrZ->SetLineColor(i+1);
171 corrX->DrawCopy(((i>0) ? "SAME" : ""));
175 corrY->DrawCopy(((i>0) ? "SAME" : ""));
179 corrZ->DrawCopy(((i>0) ? "SAME" : ""));
181 legend->AddEntry(corrZ, folderNames[i]);
186 //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.gif", fileName, gMax, upperPtLimit));
187 //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.eps", fileName, gMax, upperPtLimit));
190 TH1** DrawRatios(const char* fileName = "systematics.root")
192 gSystem->Load("libPWG0base");
194 TFile* file = TFile::Open(fileName);
197 canvasName.Form("DrawRatios");
198 TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
199 canvas->Divide(2, 1);
202 TH1** ptDists = new TH1*[4];
204 TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
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)
210 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderNames[i], folderNames[i]);
211 dNdEtaCorrection->LoadHistograms(fileName, folderNames[i]);
213 TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
215 gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
216 gene->GetXaxis()->SetRangeUser(-10, 10);
218 ptDists[i] = dynamic_cast<TH1*> (gene->Project3D("z")->Clone(Form("%s_z", folderNames[i])));
219 ptDists[i]->SetTitle(";p_{T};Count");
222 printf("Problem getting distribution %d.\n", i);
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);
232 legend->AddEntry(ptDists[i], particleNames[i]);
236 TH1* total = dynamic_cast<TH1*> (ptDists[0]->Clone("total"));
238 for (Int_t i=1; i<4; ++i)
239 total->Add(ptDists[i]);
242 for (Int_t i=0; i<4; ++i)
244 ptDists[i]->Divide(total);
245 ptDists[i]->SetStats(kFALSE);
246 ptDists[i]->SetTitle(";p_{T};Fraction of total");
247 ptDists[i]->GetYaxis()->SetRangeUser(0, 1);
248 ptDists[i]->Draw((i == 0) ? "" : "SAME");
250 legend->SetFillColor(0);
253 canvas->SaveAs("DrawRatios.gif");
256 canvas = new TCanvas("PythiaRatios", "PythiaRatios", 500, 500);
257 for (Int_t i=0; i<4; ++i)
259 TH1* hist = ptDists[i]->Clone();
260 hist->GetXaxis()->SetRangeUser(0, 1.9);
261 hist->Draw((i == 0) ? "" : "SAME");
265 canvas->SaveAs("pythiaratios.eps");
272 void DrawCompareToReal()
274 gROOT->ProcessLine(".L drawPlots.C");
276 const char* fileNames[] = { "correction_map.root", "new_compositions.root" };
277 const char* folderNames[] = { "dndeta_correction", "PythiaRatios" };
279 Track2Particle1DComposition(fileNames, 2, folderNames);
282 void DrawDifferentSpecies()
284 gROOT->ProcessLine(".L drawPlots.C");
286 const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "systematics.root" };
287 const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
289 Track2Particle1DComposition(fileNames, 4, folderNames);
292 void DrawpiKpAndCombined()
294 gROOT->ProcessLine(".L drawPlots.C");
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" };
299 Track2Particle1DComposition(fileNames, 4, folderNames);
302 void DrawSpeciesAndCombination()
304 gROOT->ProcessLine(".L drawPlots.C");
306 const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "new_compositions.root" };
307 const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "PythiaRatios" };
309 Track2Particle1DComposition(fileNames, 4, folderNames);
312 void DrawSimulatedVsCombined()
314 gROOT->ProcessLine(".L drawPlots.C");
316 const char* fileNames[] = { "new_compositions.root", "new_compositions.root" };
317 const char* folderNames[] = { "Pythia", "PythiaRatios" };
319 Track2Particle1DComposition(fileNames, 2, folderNames);
324 gROOT->ProcessLine(".L drawPlots.C");
326 const char* fileNames[] = { "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root" };
327 const char* folderNames[] = { "PythiaRatios", "PiBoosted", "KBoosted", "pBoosted" };
329 Track2Particle1DComposition(fileNames, 4, folderNames);
332 TH2F* HistogramDifferences(const char* filename, const char* folder1, const char* folder2)
334 gSystem->Load("libPWG0base");
336 AlidNdEtaCorrection* fdNdEtaCorrection[2];
338 TFile::Open(filename);
340 fdNdEtaCorrection[0] = new AlidNdEtaCorrection(folder1, folder1);
341 fdNdEtaCorrection[0]->LoadHistograms(filename, folder1);
343 fdNdEtaCorrection[1] = new AlidNdEtaCorrection(folder2, folder2);
344 fdNdEtaCorrection[1]->LoadHistograms(filename, folder2);
346 TH3F* hist1 = fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
347 TH3F* hist2 = fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
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());
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)
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));
358 difference->GetYaxis()->SetRangeUser(-0.8, 0.8);
360 printf("Over-/Underflow bins: %d %d\n", difference->GetBinContent(0), difference->GetBinContent(difference->GetNbinsX()+1));
365 void HistogramDifferences()
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");
372 TCanvas* canvas = new TCanvas("HistogramDifferences", "HistogramDifferences", 1000, 1000);
373 canvas->Divide(2, 2);
376 KBoosted->GetXaxis()->SetRangeUser(-0.05, 0.05);
377 KBoosted->Draw("COLZ");
380 KReduced->GetXaxis()->SetRangeUser(-0.05, 0.05);
381 KReduced->Draw("COLZ");
384 pBoosted->GetXaxis()->SetRangeUser(-0.02, 0.02);
385 pBoosted->Draw("COLZ");
388 pReduced->GetXaxis()->SetRangeUser(-0.02, 0.02);
389 pReduced->Draw("COLZ");
391 canvas->SaveAs("HistogramDifferences.gif");
393 canvas = new TCanvas("HistogramDifferencesProfile", "HistogramDifferencesProfile", 1000, 500);
394 canvas->Divide(2, 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");
408 TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
410 for (Int_t i=0; i<4; ++i)
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;
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" : ""));
434 legend->AddEntry(profile, name);
438 canvas->SaveAs("particlecomposition_result.eps");
442 void ScalePtDependent(TH3F* hist, TF1* function)
444 // assumes that pt is the third dimension of hist
445 // scales with function(pt)
447 for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
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);
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);
458 void ScalePtDependent(TH3F* hist, TH1* function)
460 // assumes that pt is the third dimension of hist
461 // scales with histogram(pt)
463 for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
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);
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);
474 const char* ChangeComposition(void** correctionPointer, Int_t index)
476 AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer;
480 case 0: // result from pp events
482 TFile::Open("pythiaratios.root");
484 for (Int_t i=0; i<4; ++i)
487 name.Form("correction_%d", i);
488 fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
489 fdNdEtaCorrection[i]->LoadHistograms("pythiaratios.root", name);
495 case 1: // each species rated with pythia ratios
496 /*TH1** ptDists = DrawRatios("pythiaratios.root");
498 for (Int_t i=0; i<3; ++i)
500 ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]);
501 ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]);
503 return "PythiaRatios";
506 // one species enhanced / reduced
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
515 TString* str = new TString;
518 Int_t correctionIndex = index / 2;
519 Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
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");
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");
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");
539 // each species rated with pythia ratios
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
546 TH1** ptDists = DrawRatios("pythiaratios.root");
547 Int_t functionIndex = (index - 2) / 2;
548 Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
549 ptDists[functionIndex]->Scale(scaleFactor);
551 for (Int_t i=0; i<3; ++i)
553 ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram(), ptDists[i]);
554 ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram(), ptDists[i]);
556 TString* str = new TString;
557 str->Form("%s%s", (functionIndex == 0) ? "Pi" : ((functionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced");
562 TF1* ptDependence = new TF1("simple", "x", 0, 100);
563 ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram(), ptDependence);
564 ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram(), ptDependence);
576 gSystem->Unlink("new_compositions.root");
577 gSystem->Unlink("new_compositions_analysis.root");
579 const char* names[] = { "pi", "K", "p", "other" };
582 Int_t nCompositions = 10;
584 for (Int_t comp = 1; comp < nCompositions; ++comp)
586 AlidNdEtaCorrection* fdNdEtaCorrection[4];
588 TFile::Open("correction_mapparticle-species.root");
590 for (Int_t i=0; i<4; ++i)
593 name.Form("dndeta_correction_%s", names[i]);
594 fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
595 fdNdEtaCorrection[i]->LoadHistograms();
598 const char* newName = ChangeComposition(fdNdEtaCorrection, comp);
600 Double_t geneCount[5];
601 Double_t measCount[5];
605 for (Int_t i=0; i<4; ++i)
607 geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram()->Integral();
608 measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Integral();
610 geneCount[4] += geneCount[i];
611 measCount[4] += measCount[i];
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]);
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]);
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]);
620 TList* collection = new TList;
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)
625 collection->Add(fdNdEtaCorrection[i]);
627 fdNdEtaCorrection[0]->Merge(collection);
628 fdNdEtaCorrection[0]->Finish();
633 TFile* file = TFile::Open("new_compositions.root", "UPDATE");
634 fdNdEtaCorrection[0]->SetName(newName);
635 fdNdEtaCorrection[0]->SaveHistograms();
639 // correct dNdeta distribution with modified correction map
640 TFile::Open("analysis_esd_raw.root");
642 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
643 fdNdEtaAnalysis->LoadHistograms();
645 fdNdEtaAnalysis->Finish(fdNdEtaCorrection[0], 0.2, 3, newName);
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}");
652 hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);
654 file = TFile::Open("new_compositions_analysis.root", "UPDATE");
655 hRatios[counter]->Write();
658 delete fdNdEtaAnalysis;
664 gROOT->ProcessLine(".L drawPlots.C");
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" };
667 const char* folderNames[] = { "PythiaRatios", "PiBoosted", "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced" };
669 Track2Particle1DComposition(fileNames, nCompositions, folderNames);
674 void drawSystematics()
677 //DrawDifferentSpecies();
680 Sigma2VertexSimulation();
684 void DrawdNdEtaDifferences()
688 TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
689 legend->SetFillColor(0);
691 TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
692 canvas->Divide(2, 1);
696 for (Int_t i=0; i<5; ++i)
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;
714 hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
715 hists[i]->SetTitle("a)");
717 Prepare1DPlot(hists[i], kFALSE);
718 hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
719 hists[i]->GetYaxis()->SetRangeUser(6, 7);
720 hists[i]->SetLineColor(colors[i]);
721 hists[i]->SetMarkerColor(colors[i]);
722 hists[i]->SetMarkerStyle(markers[i]);
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" : ""));
728 legend->AddEntry(hists[i], title);
729 hists[i]->SetTitle(title);
735 gPad->SetLeftMargin(0.14);
737 TLegend* legend2 = new TLegend(0.73, 0.73, 0.98, 0.98);
738 legend2->SetFillColor(0);
740 for (Int_t i=1; i<5; ++i)
744 legend2->AddEntry(hists[i]);
746 hists[i]->Divide(hists[0]);
747 hists[i]->SetTitle("b)");
748 hists[i]->SetLineColor(colors[i-1]);
749 hists[i]->SetMarkerColor(colors[i-1]);
750 hists[i]->GetYaxis()->SetRangeUser(0.95, 1.05);
751 hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
752 hists[i]->GetYaxis()->SetTitleOffset(1.8);
753 hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
759 canvas->SaveAs("particlecomposition_result_detail.gif");
761 TCanvas* canvas2 = new TCanvas("DrawdNdEtaDifferences2", "DrawdNdEtaDifferences2", 700, 500);
763 for (Int_t i=1; i<5; ++i)
767 hists[i]->SetTitle("");
768 hists[i]->GetYaxis()->SetTitleOffset(1.1);
769 hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
775 canvas2->SaveAs("particlecomposition_result.gif");
776 canvas2->SaveAs("particlecomposition_result.eps");
779 void 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") {
781 // Function used to merge standard corrections with vertex
782 // reconstruction corrections obtained by a certain mix of ND, DD
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.
789 // correctionTarget is of type AlidNdEtaCorrection::CorrectionType
795 const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
799 const Char_t* changes[] = { "pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdlessddmore25", "sdmoreddless25"};
800 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 };
801 //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};
802 //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};
803 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};
804 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 const Char_t* changes[] = { "pythia", "qgsm", "phojet"};
809 Float_t scalesND[] = {1.0, 1.10, 1.11};
810 Float_t scalesSD[] = {1.0, 0.69, 0.86};
811 Float_t scalesDD[] = {1.0, 0.98, 0.61};
815 // cross section from Pythia
817 Float_t sigmaND = 55.2;
818 Float_t sigmaDD = 9.78;
819 Float_t sigmaSD = 14.30;
821 // standard correction
822 TFile::Open(correctionFileName);
823 AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
824 correctionStandard->LoadHistograms();
826 // dont take vertexreco from this one
827 correctionStandard->GetVertexRecoCorrection()->Reset();
828 // dont take triggerbias from this one
829 correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
830 correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
831 correctionStandard->GetTriggerBiasCorrectionND()->Reset();
833 AlidNdEtaCorrection* corrections[100];
837 for (Int_t j=0; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)
839 for (Int_t i=0; i<nChanges; i++) {
840 TFile::Open(correctionFileName);
843 name.Form("dndeta_correction_syst_%s_%s", typeName[j], changes[i]);
844 AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
845 current->LoadHistograms("dndeta_correction");
848 name.Form("dndeta_correction_ND");
849 AlidNdEtaCorrection* dNdEtaCorrectionND = new AlidNdEtaCorrection(name,name);
850 dNdEtaCorrectionND->LoadHistograms();
851 name.Form("dndeta_correction_DD");
852 AlidNdEtaCorrection* dNdEtaCorrectionDD = new AlidNdEtaCorrection(name,name);
853 dNdEtaCorrectionDD->LoadHistograms();
854 name.Form("dndeta_correction_SD");
855 AlidNdEtaCorrection* dNdEtaCorrectionSD = new AlidNdEtaCorrection(name,name);
856 dNdEtaCorrectionSD->LoadHistograms();
858 // calculating relative
859 Float_t nd = 100 * sigmaND/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
860 Float_t dd = 100 * (scalesDD[i]*sigmaDD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
861 Float_t sd = 100 * (scalesSD[i]*sigmaSD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
863 printf(Form("%s : ND=%.2f\%, DD=%.2f\%, SD=%.2f\% \n",changes[i],nd,dd,sd));
864 current->SetTitle(Form("ND=%.2f\%,DD=%.2f\%,SD=%.2f\%",nd,dd,sd));
865 current->SetTitle(name);
868 if (j == 0 || j == 2)
870 dNdEtaCorrectionND->GetVertexRecoCorrection()->Scale(scalesND[i]);
871 dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scalesDD[i]);
872 dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scalesSD[i]);
874 if (j == 1 || j == 2)
876 dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->Scale(scalesND[i]);
877 dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scalesDD[i]);
878 dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scalesSD[i]);
880 dNdEtaCorrectionND->GetTriggerBiasCorrectionNSD()->Scale(scalesND[i]);
881 dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scalesDD[i]);
882 dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scalesSD[i]);
884 dNdEtaCorrectionND->GetTriggerBiasCorrectionND()->Scale(scalesND[i]);
885 dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scalesDD[i]);
886 dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scalesSD[i]);
889 //clear track in correction
890 dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
891 dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
892 dNdEtaCorrectionSD->GetTrack2ParticleCorrection()->Reset();
895 collection.Add(correctionStandard);
896 collection.Add(dNdEtaCorrectionND);
897 collection.Add(dNdEtaCorrectionDD);
898 collection.Add(dNdEtaCorrectionSD);
900 current->Merge(&collection);
903 corrections[counter] = current;
905 // now correct dNdeta distribution with modified correction map
906 TFile::Open(analysisFileName);
908 dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
909 fdNdEtaAnalysis->LoadHistograms();
911 fdNdEtaAnalysis->Finish(current, 0.3, correctionTarget, Form("%d %d", j, i));
914 if (j==0) name.Append("_vetexReco_");
915 if (j==1) name.Append("_triggerBias_");
916 if (j==2) name.Append("_vertexReco_triggerBias_");
917 name.Append(changes[i]);
919 hRatios[counter] = (TH1F*) fdNdEtaAnalysis->GetdNdEtaHistogram()->Clone(name);
921 name.Form("DD #times %0.2f, SD #times %0.2f",scalesDD[i],scalesSD[i]);
922 hRatios[counter]->SetTitle(name.Data());
923 hRatios[counter]->SetYTitle("dN_{ch}/d#eta ratio #frac{default cross-section}{modified cross-sections}");
926 hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);
928 delete fdNdEtaAnalysis;
934 TFile* fout = new TFile(outputFileName,"RECREATE");
936 // to make everything consistent
937 hRatios[0]->Divide(hRatios[0],hRatios[0],1,1);
939 for (Int_t i=0; i<nChanges * 3; i++)
941 corrections[i]->SaveHistograms();
949 DrawTriggerEfficiency(Char_t* fileName) {
951 gStyle->SetOptStat(0);
952 gStyle->SetOptTitle(0);
953 gStyle->SetOptFit(0);
955 gStyle->SetTextSize(0.04);
956 gStyle->SetTitleSize(0.05,"xyz");
957 //gStyle->SetTitleFont(133, "xyz");
958 //gStyle->SetLabelFont(133, "xyz");
959 //gStyle->SetLabelSize(17, "xyz");
960 gStyle->SetLabelOffset(0.01, "xyz");
962 gStyle->SetTitleOffset(1.1, "y");
963 gStyle->SetTitleOffset(1.1, "x");
964 gStyle->SetEndErrorSize(0.0);
966 //##############################################
968 //making canvas and pads
969 TCanvas *c = new TCanvas("trigger_eff", "",600,500);
971 TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
973 p1->SetBottomMargin(0.15);
974 p1->SetTopMargin(0.03);
975 p1->SetLeftMargin(0.15);
976 p1->SetRightMargin(0.03);
984 gSystem->Load("libPWG0base");
986 AlidNdEtaCorrection* corrections[4];
987 AliCorrectionMatrix2D* triggerBiasCorrection[4];
989 TH1F* hTriggerEffInv[4];
990 TH1F* hTriggerEff[4];
992 Char_t* names[] = {"triggerBiasND", "triggerBiasDD", "triggerBiasSD", "triggerBiasINEL"};
994 for (Int_t i=0; i<4; i++) {
995 corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
996 corrections[i]->LoadHistograms(fileName, names[i]);
998 triggerBiasCorrection[i] = corrections[i]->GetTriggerBiasCorrectionINEL();
1001 hTriggerEffInv[i] = triggerBiasCorrection[i]->Get1DCorrection();
1002 hTriggerEff[i] = (TH1F*)hTriggerEffInv[i]->Clone();
1004 for (Int_t b=0; b<=hTriggerEff[i]->GetNbinsX(); b++) {
1005 hTriggerEff[i]->SetBinContent(b,1);
1006 hTriggerEff[i]->SetBinError(b,0);
1008 hTriggerEff[i]->Divide(hTriggerEff[i],hTriggerEffInv[i]);
1009 hTriggerEff[i]->Scale(100);
1012 Int_t colors[] = {2,3,4,1};
1014 for (Int_t i=0; i<4; i++) {
1015 hTriggerEff[i]->Fit("pol0","","",-20,20);
1016 effs[i] = ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->GetParameter(0);
1017 ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineWidth(1);
1018 ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineStyle(2);
1019 ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineColor(colors[i]);
1020 cout << effs[i] << endl;
1024 Char_t* text[] = {"ND", "DD", "SD", "INEL"};
1027 TH2F* null = new TH2F("","",100,-25,35,100,0,110);
1028 null->SetXTitle("Vertex z [cm]");
1029 null->GetXaxis()->CenterTitle(kTRUE);
1030 null->SetYTitle("Trigger efficiency [%]");
1031 null->GetYaxis()->CenterTitle(kTRUE);
1035 for (Int_t i=0; i<4; i++) {
1036 hTriggerEff[i]->SetLineWidth(2);
1037 hTriggerEff[i]->SetLineColor(colors[i]);
1039 hTriggerEff[i]->Draw("same");
1041 latex[i] = new TLatex(22,effs[i]-1.5, Form("%s (%0.1f)",text[i],effs[i]));
1042 latex[i]->SetTextColor(colors[i]);
1049 DrawSpectraPID(Char_t* fileName) {
1051 gSystem->Load("libPWG0base");
1053 Char_t* names[] = {"correction_0", "correction_1", "correction_2", "correction_3"};
1054 AlidNdEtaCorrection* corrections[4];
1055 AliCorrectionMatrix3D* trackToPartCorrection[4];
1057 TH1F* measuredPt[4];
1058 TH1F* generatedPt[4];
1061 for (Int_t i=0; i<4; i++) {
1062 corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1063 corrections[i]->LoadHistograms(fileName, names[i]);
1065 trackToPartCorrection[i] = corrections[i]->GetTrack2ParticleCorrection();
1067 Int_t binX1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(-10);
1068 Int_t binX2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(10);
1069 Int_t binY1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(-1);
1070 Int_t binY2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(1);
1072 measuredPt[i] = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->ProjectionZ(Form("m_%d",i),binX1,binX2,binY1,binY2);
1073 generatedPt[i] = (TH1F*)trackToPartCorrection[i]->GetGeneratedHistogram()->ProjectionZ(Form("g_%d",i),binX1,binX2,binY1,binY2);
1074 ratioPt[i] = (TH1F*)generatedPt[i]->Clone(Form("r_%d",i));
1075 ratioPt[i]->Divide(measuredPt[i], generatedPt[i], 1,1,"B");
1080 for (Int_t i=0; i<3; i++) {
1081 ratioPt[i]->SetLineColor(i+1);
1082 ratioPt[i]->SetLineWidth(2);
1084 ratioPt[i]->Draw("same");
1089 measuredPt[0]->SetLineColor(2);
1090 measuredPt[0]->SetLineWidth(5);
1092 measuredPt[0]->Draw();
1093 generatedPt[0]->Draw("same");
1096 void changePtSpectrum(const char* fileName = "analysis_mc.root", Float_t ptCutOff = 0.2, const char* fileName2 = 0)
1098 Float_t factor = 0.5;
1100 TFile* file = TFile::Open(fileName);
1101 TH1F* hist = dynamic_cast<TH1F*> (file->Get("dndeta_check_pt")->Clone());
1106 file2 = TFile::Open(fileName2);
1107 hist2 = dynamic_cast<TH1*> (file2->Get("dndeta_check_pt")->Clone());
1108 hist2->Scale(hist->GetBinContent(hist->FindBin(ptCutOff)) / hist2->GetBinContent(hist2->FindBin(ptCutOff)));
1111 //hist->Scale(1.0 / hist->Integral());
1114 //hist->Scale(1.0/3);
1116 TH1F* clone1 = dynamic_cast<TH1F*> (hist->Clone("clone1"));
1117 TH1F* clone2 = dynamic_cast<TH1F*> (hist->Clone("clone2"));
1119 TH1F* scale1 = dynamic_cast<TH1F*> (hist->Clone("scale1"));
1120 TH1F* scale2 = dynamic_cast<TH1F*> (hist->Clone("scale2"));
1122 for (Int_t i=1; i <= hist->GetNbinsX(); ++i)
1124 if (hist->GetBinCenter(i) > ptCutOff)
1126 scale1->SetBinContent(i, 1);
1127 scale2->SetBinContent(i, 1);
1131 // 90 % at pt = 0, 0% at pt = ptcutoff
1132 scale1->SetBinContent(i, 1 - (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * factor);
1134 // 110% at pt = 0, ...
1135 scale2->SetBinContent(i, 1 + (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * factor);
1137 scale1->SetBinError(i, 0);
1138 scale2->SetBinError(i, 0);
1144 scale2->SetLineColor(kRed);
1145 scale2->Draw("SAME");
1148 clone1->Multiply(scale1);
1149 clone2->Multiply(scale2);
1151 Prepare1DPlot(hist);
1152 Prepare1DPlot(clone1);
1153 Prepare1DPlot(clone2);
1155 /*hist->SetMarkerStyle(markers[0]);
1156 clone1->SetMarkerStyle(markers[0]);
1157 clone2->SetMarkerStyle(markers[0]);*/
1159 hist->SetTitle(";p_{T} in GeV/c;dN_{ch}/dp_{T} in c/GeV");
1160 hist->GetXaxis()->SetRangeUser(0, 0.5);
1161 hist->GetYaxis()->SetRangeUser(0.01, clone2->GetMaximum() * 1.1);
1162 hist->GetYaxis()->SetTitleOffset(1);
1164 TCanvas* canvas = new TCanvas("c", "c", 600, 600);
1167 gPad->SetTopMargin(0.05);
1168 gPad->SetRightMargin(0.05);
1169 gPad->SetBottomMargin(0.12);
1171 clone1->SetLineColor(kRed);
1172 clone1->Draw("HSAME");
1173 clone2->SetLineColor(kBlue);
1174 clone2->Draw("HSAME");
1175 hist->Draw("HSAME");
1179 Prepare1DPlot(hist2);
1180 hist2->SetLineStyle(2);
1181 hist2->Draw("HSAME");
1184 Float_t fraction = hist->Integral(hist->GetXaxis()->FindBin(ptCutOff), hist->GetNbinsX()) / hist->Integral(1, hist->GetNbinsX());
1185 Float_t fraction1 = clone1->Integral(clone1->GetXaxis()->FindBin(ptCutOff), clone1->GetNbinsX()) / clone1->Integral(1, clone1->GetNbinsX());
1186 Float_t fraction2 = clone2->Integral(clone2->GetXaxis()->FindBin(ptCutOff), clone2->GetNbinsX()) / clone2->Integral(1, clone2->GetNbinsX());
1188 printf("%f %f %f\n", fraction, fraction1, fraction2);
1189 printf("Rel. %f %f\n", fraction1 / fraction, fraction2 / fraction);
1191 //canvas->SaveAs("changePtSpectrum.gif");
1192 canvas->SaveAs("changePtSpectrum.eps");
1195 void FractionBelowPt()
1197 gSystem->Load("libPWG0base");
1199 AlidNdEtaCorrection* fdNdEtaCorrection[4];
1201 TFile::Open("systematics.root");
1203 for (Int_t i=0; i<4; ++i)
1206 name.Form("correction_%d", i);
1207 fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
1208 fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
1211 Double_t geneCount[5];
1212 Double_t measCount[5];
1216 for (Int_t i=0; i<4; ++i)
1218 TH3F* hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1219 geneCount[i] = hist->Integral(hist->GetXaxis()->FindBin(-10), hist->GetXaxis()->FindBin(10),
1220 hist->GetYaxis()->FindBin(-0.8), hist->GetYaxis()->FindBin(0.8),
1221 1, hist->GetZaxis()->FindBin(0.3));
1223 hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1224 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));
1226 geneCount[4] += geneCount[i];
1227 measCount[4] += measCount[i];
1229 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]);
1232 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]);
1234 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]);
1238 mergeCorrectionsMisalignment(Char_t* alignedFile = "correction_map_aligned.root",
1239 Char_t* misalignedFile = "correction_map_misaligned.root",
1240 Char_t* outputFileName="correction_map_misaligned_single.root")
1243 // from the aligned and misaligned corrections, 3 new corrections are created
1244 // in these new corrections only one of the corrections (track2particle, vertex, trigger)
1245 // is taken from the misaligned input to allow study of the effect on the different
1248 gSystem->Load("libPWG0base");
1250 const Char_t* typeName[] = { "track2particle", "vertex", "trigger" };
1252 AlidNdEtaCorrection* corrections[3];
1253 for (Int_t j=0; j<3; j++) { // j = 0 (track2particle), j = 1 (vertex), j = 2 (trigger)
1254 AlidNdEtaCorrection* alignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1255 alignedCorrection->LoadHistograms(alignedFile);
1257 AlidNdEtaCorrection* misalignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1258 misalignedCorrection->LoadHistograms(misalignedFile);
1261 name.Form("dndeta_correction_alignment_%s", typeName[j]);
1262 AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
1267 alignedCorrection->GetTrack2ParticleCorrection()->Reset();
1268 misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1269 misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1270 misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1271 misalignedCorrection->GetVertexRecoCorrection()->Reset();
1275 misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1276 misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1277 misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1278 misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1279 alignedCorrection->GetVertexRecoCorrection()->Reset();
1283 misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1284 alignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1285 alignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1286 alignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1287 misalignedCorrection->GetVertexRecoCorrection()->Reset();
1295 collection.Add(misalignedCorrection);
1296 collection.Add(alignedCorrection);
1298 current->Merge(&collection);
1301 corrections[j] = current;
1304 TFile* fout = new TFile(outputFileName, "RECREATE");
1306 for (Int_t i=0; i<3; i++)
1307 corrections[i]->SaveHistograms();
1314 void DrawdNdEtaDifferencesAlignment()
1318 TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
1319 legend->SetFillColor(0);
1321 TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
1322 canvas->Divide(2, 1);
1326 for (Int_t i=0; i<5; ++i)
1334 case 0 : file = TFile::Open("systematics_misalignment_aligned.root"); title = "aligned"; break;
1335 case 1 : file = TFile::Open("systematics_misalignment_misaligned.root"); title = "fully misaligned"; break;
1336 case 2 : file = TFile::Open("systematics_misalignment_track2particle.root"); title = "only track2particle"; break;
1337 case 3 : file = TFile::Open("systematics_misalignment_vertex.root"); title = "only vertex rec."; break;
1338 case 4 : file = TFile::Open("systematics_misalignment_trigger.root"); title = "only trigger bias"; break;
1344 hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
1345 hists[i]->SetTitle("");
1347 Prepare1DPlot(hists[i], kFALSE);
1348 hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
1349 hists[i]->GetYaxis()->SetRangeUser(6, 7);
1350 hists[i]->SetLineWidth(1);
1351 hists[i]->SetLineColor(colors[i]);
1352 hists[i]->SetMarkerColor(colors[i]);
1353 hists[i]->SetMarkerStyle(markers[i]);
1354 hists[i]->GetXaxis()->SetLabelOffset(0.015);
1355 hists[i]->GetYaxis()->SetTitleOffset(1.5);
1356 gPad->SetLeftMargin(0.12);
1357 hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
1359 legend->AddEntry(hists[i], title);
1360 hists[i]->SetTitle(title);
1366 gPad->SetLeftMargin(0.14);
1368 TLegend* legend2 = new TLegend(0.63, 0.73, 0.98, 0.98);
1369 legend2->SetFillColor(0);
1371 for (Int_t i=1; i<5; ++i)
1375 legend2->AddEntry(hists[i]);
1377 hists[i]->Divide(hists[0]);
1378 hists[i]->SetTitle("b)");
1379 hists[i]->GetYaxis()->SetRangeUser(0.9, 1.1);
1380 hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
1381 hists[i]->GetYaxis()->SetTitleOffset(1.8);
1382 hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
1388 canvas->SaveAs("misalignment_result.eps");
1389 canvas->SaveAs("misalignment_result.gif");
1392 void EvaluateMultiplicityEffect()
1394 gSystem->Load("libPWG0base");
1396 AlidNdEtaCorrection* fdNdEtaCorrectionLow[4];
1397 AlidNdEtaCorrection* fdNdEtaCorrectionHigh[4];
1399 TFile::Open("systematics-low-multiplicity.root");
1401 for (Int_t i=0; i<4; ++i)
1404 name.Form("correction_%d", i);
1405 fdNdEtaCorrectionLow[i] = new AlidNdEtaCorrection(name, name);
1406 fdNdEtaCorrectionLow[i]->LoadHistograms("systematics-low-multiplicity.root", name);
1410 for (Int_t i=1; i<4; ++i)
1411 list.Add(fdNdEtaCorrectionLow[i]);
1413 fdNdEtaCorrectionLow[0]->Merge(&list);
1414 fdNdEtaCorrectionLow[0]->Finish();
1416 TFile::Open("systematics-high-multiplicity.root");
1418 for (Int_t i=0; i<4; ++i)
1421 name.Form("correction_%d", i);
1422 fdNdEtaCorrectionHigh[i] = new AlidNdEtaCorrection(name, name);
1423 fdNdEtaCorrectionHigh[i]->LoadHistograms("systematics-high-multiplicity.root", name);
1427 for (Int_t i=1; i<4; ++i)
1428 list2.Add(fdNdEtaCorrectionHigh[i]);
1430 fdNdEtaCorrectionHigh[0]->Merge(&list2);
1431 fdNdEtaCorrectionHigh[0]->Finish();
1433 TH1F* outputLow = new TH1F("Track2ParticleLow", "Track2Particle at low multiplicity", 200, 0, 2);
1434 TH1F* outputHigh = new TH1F("Track2ParticleHigh", "Track2Particle at high multiplicity", 200, 0, 2);
1436 TH3F* hist = fdNdEtaCorrectionLow[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1437 TH3F* hist2 = fdNdEtaCorrectionHigh[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1438 for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
1439 for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
1440 for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
1441 //for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
1443 if (hist->GetBinContent(x, y, z) > 0)
1444 outputLow->Fill(hist->GetBinContent(x, y, z));
1445 //if (hist->GetBinContent(x, y, z) == 1)
1446 // 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));
1448 if (hist2->GetBinContent(x, y, z) > 0)
1449 outputHigh->Fill(hist2->GetBinContent(x, y, z));
1452 TCanvas* canvas = new TCanvas("EvaluateMultiplicityEffect", "EvaluateMultiplicityEffect", 1000, 500);
1453 canvas->Divide(2, 1);
1457 outputLow->Fit("gaus", "0");
1458 outputLow->GetFunction("gaus")->SetLineColor(2);
1459 outputLow->GetFunction("gaus")->DrawCopy("SAME");
1463 outputHigh->Fit("gaus", "0");
1464 outputHigh->GetFunction("gaus")->DrawCopy("SAME");
1466 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1469 void 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")
1471 gSystem->Load("libPWG0base");
1473 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1474 dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
1476 dNdEtaCorrection->GetTrack2ParticleCorrection()->PlotBinErrors(xmin, xmax, ymin, ymax, zmin, zmax)->Draw();
1480 void 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")
1482 gSystem->Load("libPWG0base");
1484 TFile* file = TFile::Open(output, "RECREATE");
1486 const Int_t max = 5;
1487 dNdEtaAnalysis* fdNdEtaAnalysis[5];
1490 TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
1491 legend->SetFillColor(0);
1493 for (Int_t i = 0; i < max; ++i)
1495 TFile::Open(baseCorrectionMapFile);
1496 AlidNdEtaCorrection* baseCorrection = new AlidNdEtaCorrection(baseCorrectionMapFolder, baseCorrectionMapFolder);
1497 baseCorrection->LoadHistograms();
1499 AlidNdEtaCorrection::CorrectionType correctionType = AlidNdEtaCorrection::kNone;
1500 const char* name = 0;
1502 TFile::Open(changedCorrectionMapFile);
1510 baseCorrection->GetTrack2ParticleCorrection()->LoadHistograms(Form("%s/Track2Particle", changedCorrectionMapFolder));
1511 name = "Track2Particle";
1515 baseCorrection->GetVertexRecoCorrection()->LoadHistograms(Form("%s/VertexReconstruction", changedCorrectionMapFolder));
1516 name = "VertexReco";
1520 baseCorrection->GetTriggerBiasCorrectionINEL()->LoadHistograms(Form("%s/TriggerBias_MBToINEL", changedCorrectionMapFolder));
1521 name = "TriggerBias_MBToINEL";
1525 baseCorrection->LoadHistograms(changedCorrectionMapFolder);
1532 TFile::Open(dataFile);
1533 fdNdEtaAnalysis[i] = new dNdEtaAnalysis(name, name);
1534 fdNdEtaAnalysis[i]->LoadHistograms("dndeta");
1536 fdNdEtaAnalysis[i]->Finish(baseCorrection, 0.3, AlidNdEtaCorrection::kINEL);
1538 fdNdEtaAnalysis[i]->SaveHistograms();
1540 TH1* hist = fdNdEtaAnalysis[i]->GetdNdEtaHistogram(0);
1541 hist->SetLineColor(colors[i]);
1542 hist->SetMarkerColor(colors[i]);
1543 hist->SetMarkerStyle(markers[i]+1);
1544 hist->DrawCopy((i == 0) ? "" : "SAME");
1545 legend->AddEntry(hist, name);
1551 void ChangePtInCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
1554 if (!TFile::Open(fileName))
1557 AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
1558 if (!dNdEtaCorrection->LoadHistograms())
1561 dNdEtaCorrection->Finish();
1563 AliCorrection* corr = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kTrack2Particle);
1565 Printf(">>>> Before");
1568 Float_t factor = 0.5;
1569 Float_t ptCutOff = 0.2;
1571 TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram();
1572 TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram();
1574 for (Int_t z = 1; z <= gene->GetZaxis()->FindBin(ptCutOff - 0.01); z++)
1576 Float_t localFactor = 1 - (ptCutOff - gene->GetZaxis()->GetBinCenter(z)) / ptCutOff * factor;
1577 Printf("%f %f", gene->GetZaxis()->GetBinCenter(z), localFactor);
1578 for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1579 for (Int_t y = 1; y <= gene->GetNbinsY(); y++)
1581 gene->SetBinContent(x, y, z, gene->GetBinContent(x, y, z) * localFactor);
1582 meas->SetBinContent(x, y, z, meas->GetBinContent(x, y, z) * localFactor);
1586 dNdEtaCorrection->Finish();
1588 Printf(">>>> After");