1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 // This class is used to store correlation maps, generated and reconstructed data of the jet spectrum
17 // It also contains functions to correct the spectrum using the bayesian unfolding
20 #include "AliJetSpectrumUnfolding.h"
26 #include <THnSparse.h>
27 #include <TDirectory.h>
28 #include <TVirtualFitter.h>
34 #include <TCollection.h>
39 #include <TProfile2D.h>
44 #include <THnSparse.h>
46 ClassImp(AliJetSpectrumUnfolding)
48 const Int_t AliJetSpectrumUnfolding::fgkNBINSE = 50;
49 const Int_t AliJetSpectrumUnfolding::fgkNBINSZ = 50;
50 const Int_t AliJetSpectrumUnfolding::fgkNEVENTS = 500000;
51 const Double_t AliJetSpectrumUnfolding::fgkaxisLowerLimitE = 0.;
52 const Double_t AliJetSpectrumUnfolding::fgkaxisLowerLimitZ = 0.;
53 const Double_t AliJetSpectrumUnfolding::fgkaxisUpperLimitE = 250.;
54 const Double_t AliJetSpectrumUnfolding::fgkaxisUpperLimitZ = 1.;
56 Float_t AliJetSpectrumUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing)
57 Int_t AliJetSpectrumUnfolding::fgBayesianIterations = 100; // number of iterations in Bayesian method
59 //____________________________________________________________________
61 AliJetSpectrumUnfolding::AliJetSpectrumUnfolding() :
62 TNamed(), fCurrentRec(0), fCurrentCorrelation(0), fRecSpectrum(0), fGenSpectrum(0),
63 fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
66 // default constructor
75 //____________________________________________________________________
76 AliJetSpectrumUnfolding::AliJetSpectrumUnfolding(const Char_t* name, const Char_t* title) :
77 TNamed(name, title), fCurrentRec(0), fCurrentCorrelation(0), fRecSpectrum(0),
78 fGenSpectrum(0), fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
84 // do not add this hists to the directory
85 Bool_t oldStatus = TH1::AddDirectoryStatus();
86 TH1::AddDirectory(kFALSE);
87 fRecSpectrum = new TH2F("fRecSpectrum", "Reconstructed Spectrum;E^{jet}_{rec} [GeV];z^{lp}_{rec}",
88 fgkNBINSE, fgkaxisLowerLimitE, fgkaxisUpperLimitE,
89 fgkNBINSZ, fgkaxisLowerLimitZ, fgkaxisUpperLimitZ);
90 fGenSpectrum = new TH2F("fGenSpectrum", "Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}_{gen}",
91 fgkNBINSE, fgkaxisLowerLimitE, fgkaxisUpperLimitE,
92 fgkNBINSZ, fgkaxisLowerLimitZ, fgkaxisUpperLimitZ);
93 fUnfSpectrum = new TH2F("fUnfSpectrum", "Unfolded Spectrum;E^{jet} [GeV];z^{lp}",
94 fgkNBINSE, fgkaxisLowerLimitE, fgkaxisUpperLimitE,
95 fgkNBINSZ, fgkaxisLowerLimitZ, fgkaxisUpperLimitZ);
97 const Int_t nbin[4]={fgkNBINSE, fgkNBINSE, fgkNBINSZ, fgkNBINSZ};
98 //arrays for bin limits
99 Double_t lowEdge[4] = {fgkaxisLowerLimitE, fgkaxisLowerLimitE, fgkaxisLowerLimitZ, fgkaxisLowerLimitZ};
100 Double_t upEdge[4] = {fgkaxisUpperLimitE, fgkaxisUpperLimitE, fgkaxisUpperLimitZ, fgkaxisUpperLimitZ};
102 fCorrelation = new THnSparseF("fCorrelation", "Correlation Function", 4, nbin, lowEdge, upEdge);
104 TH1::AddDirectory(oldStatus);
107 //____________________________________________________________________
108 AliJetSpectrumUnfolding::~AliJetSpectrumUnfolding()
132 //____________________________________________________________________
133 Long64_t AliJetSpectrumUnfolding::Merge(TCollection* list)
135 // Merge a list of AliJetSpectrumUnfolding objects with this (needed for
137 // Returns the number of merged objects (including this).
145 TIterator* iter = list->MakeIterator();
148 // collections of all histograms
149 TList collections[4];
152 while ((obj = iter->Next())) {
154 AliJetSpectrumUnfolding* entry = dynamic_cast<AliJetSpectrumUnfolding*> (obj);
158 collections[0].Add(entry->fGenSpectrum);
159 collections[1].Add(entry->fRecSpectrum);
160 collections[2].Add(entry->fUnfSpectrum);
161 collections[3].Add(entry->fCorrelation);
166 fGenSpectrum->Merge(&collections[0]);
167 fRecSpectrum->Merge(&collections[1]);
168 fUnfSpectrum->Merge(&collections[2]);
169 fCorrelation->Merge(&collections[3]);
176 //____________________________________________________________________
177 Bool_t AliJetSpectrumUnfolding::LoadHistograms(const Char_t* dir)
180 // loads the histograms from a file
181 // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
187 if (!gDirectory->cd(dir))
190 Bool_t success = kTRUE;
192 // store old histograms to delete them later
194 oldHistograms.SetOwner(1);
196 if (fGenSpectrum) oldHistograms.Add(fGenSpectrum);
198 if (fRecSpectrum) oldHistograms.Add(fRecSpectrum);
200 if (fUnfSpectrum) oldHistograms.Add(fUnfSpectrum);
202 if (fCorrelation) oldHistograms.Add(fCorrelation);
205 // load new histograms
206 fGenSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fGenSpectrum->GetName()));
210 fRecSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fRecSpectrum->GetName()));
214 fUnfSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fUnfSpectrum->GetName()));
218 fCorrelation = dynamic_cast<THnSparseF*> (gDirectory->Get(fCorrelation->GetName()));
222 gDirectory->cd("..");
224 // delete old histograms
225 oldHistograms.Delete();
230 //____________________________________________________________________
231 void AliJetSpectrumUnfolding::SaveHistograms()
234 // saves the histograms
237 gDirectory->mkdir(GetName());
238 gDirectory->cd(GetName());
241 fGenSpectrum->Write();
244 fRecSpectrum->Write();
247 fUnfSpectrum->Write();
250 fCorrelation->Write();
252 gDirectory->cd("..");
255 //____________________________________________________________________
256 void AliJetSpectrumUnfolding::SetupCurrentHists(Bool_t createBigBin)
259 // resets fUnfSpectrum
262 fUnfSpectrum->Reset();
263 fUnfSpectrum->Sumw2();
265 fCurrentRec = (TH2F*)fRecSpectrum->Clone("fCurrentRec");
266 fCurrentRec->Sumw2();
268 fCurrentCorrelation = (THnSparseF*)fCorrelation->Clone("fCurrentCorrelation");
269 fCurrentCorrelation->Sumw2();
271 Printf("Correlation Matrix has %d filled bins", (Int_t)fCurrentCorrelation->GetNbins());
275 Int_t maxBinE = 0, maxBinZ = 0;
276 Float_t maxE = 0, maxZ = 0;
277 for (Int_t me=1; me<=fCurrentRec->GetNbinsX(); me++)
278 for (Int_t mz=1; mz<=fCurrentRec->GetNbinsY(); mz++)
280 if (fCurrentRec->GetBinContent(me,mz) <= 5 && me>fgkNBINSE/2 && mz>fgkNBINSZ/2)
284 maxE = fCurrentRec->GetXaxis()->GetBinCenter(me);
285 maxZ = fCurrentRec->GetYaxis()->GetBinCenter(mz);
290 if (maxBinE > 0 || maxBinZ > 0)
292 printf("Bin limit in measured spectrum is e = %d and z = %d.\n", maxBinE, maxBinZ);
293 fCurrentRec->SetBinContent(maxBinE, maxBinZ, fCurrentRec->Integral(maxBinE, fCurrentRec->GetNbinsX(), maxBinZ, fCurrentRec->GetNbinsY()));
294 for (Int_t me=maxBinE+1; me<=fCurrentRec->GetNbinsX(); me++)
295 for (Int_t mz=maxBinZ+1; mz<=fCurrentRec->GetNbinsY(); mz++)
297 fCurrentRec->SetBinContent(me, mz, 0);
298 fCurrentRec->SetBinError(me, mz, 0);
300 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
301 fCurrentRec->SetBinError(maxBinE, maxBinZ, TMath::Sqrt(fCurrentRec->GetBinContent(maxBinE, maxBinZ)));
303 printf("This bin has now %f +- %f entries\n", fCurrentRec->GetBinContent(maxBinE, maxBinZ), fCurrentRec->GetBinError(maxBinE, maxBinZ));
305 /* for (Int_t te=1; te<=NBINSE; te++)
307 for (Int_t tz=1; tz<=NBINSZ; tz++)
309 Int_t binMin[4] = {te, maxBinE, tz, maxBinZ};
310 Int_t binMax[4] = {NBINSE, NBINSE, NBINSZ, NBINSZ};
312 for (Int_t ite=te; ite<=NBINSE; ite++)
313 for (Int_t itz=tz; itz<=NBINSZ; itz++)
314 for (Int_t ime=maxBinE; ime<=NBINSE; ime++)
315 for (Int_t imz=maxBinZ; imz<=NBINSZ; imz++)
317 Int_t bin[4] = {ite, ime, itz, imz};
318 sum += fCurrentCorrelation->GetBinContent(bin);
320 fCurrentCorrelation->SetBinContent(binMin, sum);
321 fCurrentCorrelation->SetBinError(binMin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin)));
322 printf("create big bin1, nbins = %d, te = %d, tz = %d\n", NBINSE, te, tz);
323 for (Int_t me=maxBinE; me<=NBINSE; me++)
325 for (Int_t mz=maxBinZ; mz<=NBINSZ; mz++)
327 Int_t bin[4] = {te, me, tz, mz};
328 fCurrentCorrelation->SetBinContent(bin, 0.);
329 fCurrentCorrelation->SetBinError(bin, 0.);
330 printf("create big bin2\n");
336 for(Int_t idx = 0; idx<=fCurrentCorrelation->GetNbins(); idx++)
339 Float_t binContent = fCurrentCorrelation->GetBinContent(idx,bin);
340 Float_t binError = fCurrentCorrelation->GetBinError(idx);
341 Int_t binMin[4] = {bin[0], maxBinE, bin[2], maxBinZ};
342 if ( (bin[1]>maxBinE && bin[1]<=fgkNBINSE) && (bin[3]>maxBinZ && bin[3]<=fgkNBINSZ) )
344 fCurrentCorrelation->SetBinContent(binMin, binContent + fCurrentCorrelation->GetBinContent(binMin));
345 fCurrentCorrelation->SetBinError(binMin, binError + TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin)));
346 fCurrentCorrelation->SetBinContent(bin, 0.);
347 fCurrentCorrelation->SetBinError(bin, 0.);
349 printf("create big bin1, nbins = %d, te = %d, tz = %d\n", fgkNBINSE, bin[0], bin[1]);
352 printf("Adjusted correlation matrix!\n");
354 } // end Create Big Bin
358 //____________________________________________________________________
359 void AliJetSpectrumUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
362 // sets the parameters for Bayesian unfolding
365 fgBayesianSmoothing = smoothing;
366 fgBayesianIterations = nIterations;
368 printf("AliJetSpectrumUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
371 //____________________________________________________________________
372 void AliJetSpectrumUnfolding::NormalizeToBinWidth(TH2* const hist)
375 // normalizes a 2-d histogram to its bin width (x width * y width)
378 for (Int_t i=1; i<=hist->GetNbinsX(); i++)
379 for (Int_t j=1; j<=hist->GetNbinsY(); j++)
381 Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
382 hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
383 hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
387 //____________________________________________________________________
388 void AliJetSpectrumUnfolding::DrawHistograms()
391 // draws the histograms of this class
394 gStyle->SetPalette(1);
396 TCanvas* canvas1 = new TCanvas("fRecSpectrum", "fRecSpectrum", 900, 600);
398 fRecSpectrum->DrawCopy("COLZ");
400 TCanvas* canvas2 = new TCanvas("fGenSpectrum", "fGenSpectrum", 900, 600);
403 fGenSpectrum->DrawCopy("COLZ");
405 TCanvas* canvas3 = new TCanvas("fUnfSpectrum", "fUnfSpectrum", 900, 600);
408 fUnfSpectrum->DrawCopy("COLZ");
410 TCanvas* canvas4 = new TCanvas("fCorrelation", "fCorrelation", 500, 500);
415 TH2D* h0 = fCorrelation->Projection(1,0);
416 h0->SetXTitle("E^{jet}_{gen} [GeV]");
417 h0->SetYTitle("E^{jet}_{rec} [GeV]");
418 h0->SetTitle("Projection: Jet Energy");
419 h0->DrawCopy("colz");
423 TH2D* h1 = fCorrelation->Projection(3,2);
424 h1->SetXTitle("z^{lp}_{gen}");
425 h1->SetYTitle("z^{lp}_{rec}");
426 h1->SetTitle("Projection: Leading Particle Fragmentation");
427 h1->DrawCopy("colz");
431 //____________________________________________________________________
432 void AliJetSpectrumUnfolding::DrawComparison(const char* name, TH2* const genHist)
435 // Draws the copmparison plot (gen,rec and unfolded distributions
438 if (fUnfSpectrum->Integral() == 0)
440 printf("ERROR. Unfolded histogram is empty\n");
444 //regain measured distribution used for unfolding, because the bins were modified in SetupCurrentHists
446 fCurrentRec = (TH2F*)fRecSpectrum->Clone();
447 fCurrentRec->Sumw2();
448 fCurrentRec->Scale(1.0 / fCurrentRec->Integral());
450 // normalize unfolded result to 1
451 fUnfSpectrum->Scale(1.0 / fUnfSpectrum->Integral());
453 // find bin with <= 100 entries. this is used as limit for the chi2 calculation
454 Int_t mcBinLimitE = 0, mcBinLimitZ = 0;
455 for (Int_t i=0; i<genHist->GetNbinsX(); ++i)
456 for (Int_t j=0; j<genHist->GetNbinsY(); ++j)
458 if (genHist->GetBinContent(i,j) > 100)
466 Printf("AliJetSpectrumUnfolding::DrawComparison: Gen bin limit is (x,y) = (%d,%d)", mcBinLimitE,mcBinLimitZ);
468 // scale to 1 true spectrum
470 genHist->Scale(1.0 / genHist->Integral());
472 // calculate residual
473 // for that we convolute the response matrix with the gathered result
474 TH2* tmpRecRecorrected = (TH2*) fUnfSpectrum->Clone("tmpRecRecorrected");
475 TH2* convoluted = CalculateRecSpectrum(tmpRecRecorrected);
476 if (convoluted->Integral() > 0)
477 convoluted->Scale(1.0 / convoluted->Integral());
479 printf("ERROR: convoluted is empty. Something went wrong calculating the convoluted histogram.\n");
481 TH2* residual = (TH2*) convoluted->Clone("residual");
482 residual->SetTitle("(R#otimesUnfolded - Reconstructed)/Reconstructed;E^{jet} [GeV]; z^{lp}");
484 fCurrentRec->Scale(1./fCurrentRec->Integral());
485 residual->Add(fCurrentRec, -1);
486 //residual->Divide(residual, fCurrentRec, 1, 1, "B");
489 TCanvas* canvas1 = new TCanvas(name, name, 1000, 1000);
490 canvas1->Divide(2, 3);
493 const Int_t nRGBs = 5;
494 const Int_t nCont = 500;
496 Double_t stops[nRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
497 Double_t red[nRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
498 Double_t green[nRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
499 Double_t blue[nRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
500 TColor::CreateGradientColorTable(nRGBs, stops, red, green, blue, nCont);
501 gStyle->SetNumberContours(nCont);
504 gStyle->SetPalette(style);
506 genHist->SetTitle("Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}");
507 genHist->SetStats(0);
508 genHist->DrawCopy("colz");
511 gStyle->SetPalette(style);
513 fUnfSpectrum->SetStats(0);
514 fUnfSpectrum->DrawCopy("colz");
517 gStyle->SetPalette(style);
519 fCurrentRec->SetTitle(fRecSpectrum->GetTitle());
520 fCurrentRec->SetStats(0);
521 fCurrentRec->DrawCopy("colz");
524 gStyle->SetPalette(style);
526 TH1D* projGenX = genHist->ProjectionX();
527 projGenX->SetTitle("Projection: Jet Energy; E^{jet} [GeV]");
528 TH1D* projUnfX = fUnfSpectrum->ProjectionX();
529 TH1D* projRecX = fCurrentRec->ProjectionX();
530 projGenX->SetStats(0);
531 projRecX->SetStats(0);
532 projUnfX->SetStats(0);
533 projGenX->SetLineColor(8);
534 projRecX->SetLineColor(2);
535 projGenX->DrawCopy();
536 projUnfX->DrawCopy("same");
537 projRecX->DrawCopy("same");
539 TLegend* legend = new TLegend(0.6, 0.85, 0.98, 0.98);
540 legend->AddEntry(projGenX, "Generated Spectrum");
541 legend->AddEntry(projUnfX, "Unfolded Spectrum");
542 legend->AddEntry(projRecX, "Reconstructed Spectrum");
543 //legend->SetFillColor(0);
544 legend->Draw("same");
548 gStyle->SetPalette(style);
549 TH1D* projGenY = genHist->ProjectionY();
550 projGenY->SetTitle("Projection: Leading Particle Fragmentation; z^{lp}");
551 TH1D* projUnfY = fUnfSpectrum->ProjectionY();
552 TH1D* projRecY = fCurrentRec->ProjectionY();
553 projGenY->SetStats(0);
554 projRecY->SetStats(0);
555 projUnfY->SetStats(0);
556 projGenY->SetLineColor(8);
557 projRecY->SetLineColor(2);
558 projGenY->DrawCopy();
559 projUnfY->DrawCopy("same");
560 projRecY->DrawCopy("same");
562 TLegend* legend1 = new TLegend(0.6, 0.85, 0.98, 0.98);
563 legend1->AddEntry(projGenY, "Generated Spectrum");
564 legend1->AddEntry(projUnfY, "Unfolded Spectrum");
565 legend1->AddEntry(projRecY, "Recontructed Spectrum");
566 //legend1->SetFillColor(0);
567 legend1->Draw("same");
571 gStyle->SetPalette(style);
573 residual->SetStats(0);
574 residual->DrawCopy("colz");
576 canvas1->SaveAs(Form("%s.png", canvas1->GetName()));
580 //____________________________________________________________________
581 void AliJetSpectrumUnfolding::ApplyBayesianMethod(Float_t regPar, Int_t nIterations, TH2* const initialConditions, Bool_t determineError)
584 // correct spectrum using bayesian unfolding
587 // initialize seed with current time
590 printf("seting up current arrays and histograms...\n");
591 SetupCurrentHists(kFALSE); // kFALSE to not create big bin
593 // normalize Correlation Map to convert number of events into probabilities
594 /*for (Int_t te=1; te<=NBINSE; te++)
595 for (Int_t tz=1; tz<=NBINSZ; tz++)
599 for (Int_t me = 1; me<=NBINSE; me++)
600 for (Int_t mz = 1; mz<=NBINSZ; mz++)
602 bin[0] = te; bin[1] = me;
603 bin[2] = tz; bin[3] = mz;
604 sum += fCurrentCorrelation->GetBinContent(bin);
607 for (Int_t me = 1; me<=NBINSE; me++)
608 for (Int_t mz = 1; mz<=NBINSZ; mz++)
610 bin[0] = te; bin[1] = me;
611 bin[2] = tz; bin[3] = mz;
612 fCurrentCorrelation->SetBinContent(bin, fCurrentCorrelation->GetBinContent(bin)/sum);
613 fCurrentCorrelation->SetBinError(bin, fCurrentCorrelation->GetBinError(bin)/sum);
616 Float_t sum[fgkNBINSE+2][fgkNBINSZ+2];
617 memset(sum,0,sizeof(Float_t)*(fgkNBINSE+2)*(fgkNBINSZ+2));
619 for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++)
622 Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin);
623 if ( (bin[1]>0 && bin[1]<=fgkNBINSE) && (bin[3]>0 && bin[3]<=fgkNBINSZ) )
624 sum[bin[0]][bin[2]] += binContent;
627 for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++)
630 Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin);
631 Float_t binError = fCurrentCorrelation->GetBinError(bin);
632 if (sum[bin[0]][bin[2]]>0 && (bin[1]>0 && bin[1]<=fgkNBINSE) &&
633 (bin[3]>0 && bin[3]<=fgkNBINSZ) && (bin[0]>0 && bin[0]<=fgkNBINSE) && (bin[2]>0 && bin[2]<=fgkNBINSZ) )
635 fCurrentCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]);
636 fCurrentCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]);
640 printf("calling UnfoldWithBayesian\n");
641 Int_t success = UnfoldWithBayesian(fCurrentCorrelation, fCurrentRec, initialConditions, fUnfSpectrum, regPar, nIterations, kFALSE);
648 Printf("AliJetSpectrumUnfolding::ApplyBayesianMethod: WARNING: No errors calculated.");
652 // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
653 // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
654 // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic
655 // error of the unfolding method itself.
657 TH2* maxError = (TH2*) fUnfSpectrum->Clone("maxError");
660 TH2* normalizedResult = (TH2*) fUnfSpectrum->Clone("normalizedResult");
661 normalizedResult->Scale(1.0 / normalizedResult->Integral());
663 const Int_t kErrorIterations = 20;
665 printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
667 TH2* randomized = (TH2*) fCurrentRec->Clone("randomized");
668 TH2* result2 = (TH2*) fUnfSpectrum->Clone("result2");
669 for (Int_t n=0; n<kErrorIterations; ++n)
671 // randomize the content of clone following a poisson with the mean = the value of that bin
672 for (Int_t x=1; x<=randomized->GetNbinsX(); x++)
673 for (Int_t y=1; y<=randomized->GetNbinsY(); y++)
675 Float_t randomValue = fCurrentRec->GetBinContent(x,y);
676 TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",randomValue*0.25, randomValue*1.25);
677 poisson->SetParameters(randomValue,0.);
678 randomValue = poisson->GetRandom();
679 //printf("%e --> %e\n", fCurrentRec->GetBinContent(x,y), (Double_t)randomValue);
680 randomized->SetBinContent(x, y, randomValue);
685 if (UnfoldWithBayesian(fCurrentCorrelation, randomized, initialConditions, result2, regPar, nIterations) != 0)
688 result2->Scale(1.0 / result2->Integral());
691 result2->Divide(normalizedResult);
693 //new TCanvas; result2->DrawCopy("HIST");
695 // find max. deviation
696 for (Int_t i=1; i<=result2->GetNbinsX(); i++)
697 for (Int_t j=1; j<=result2->GetNbinsY(); j++)
698 maxError->SetBinContent(i, j, TMath::Max(maxError->GetBinContent(i,j), TMath::Abs(1 - result2->GetBinContent(i,j))));
703 for (Int_t i=1; i<=fUnfSpectrum->GetNbinsX(); i++)
704 for (Int_t j=1; j<=fUnfSpectrum->GetNbinsY(); j++)
705 fUnfSpectrum->SetBinError(i, j, fUnfSpectrum->GetBinError(i,j) + maxError->GetBinContent(i,j)*fUnfSpectrum->GetBinContent(i,j));
708 delete normalizedResult;
711 //____________________________________________________________________
712 Int_t AliJetSpectrumUnfolding::UnfoldWithBayesian(THnSparseF* const correlation, TH2* const measured, TH2* const initialConditions, TH2* const aResult, Float_t regPar, Int_t nIterations, Bool_t calculateErrors)
715 // unfolds a spectrum
718 if (measured->Integral() <= 0)
720 Printf("AliJetSpectrumUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
723 const Int_t nFillesBins = correlation->GetNbins();
724 const Int_t kStartBin = 1;
726 const Int_t kMaxTZ = fgkNBINSZ; // max true axis fragmentation function
727 const Int_t kMaxMZ = fgkNBINSZ; // max measured axis fragmentation function
728 const Int_t kMaxTE = fgkNBINSE; // max true axis energy
729 const Int_t kMaxME = fgkNBINSE; // max measured axis energy
731 printf("NbinsE=%d - NbinsZ=%d\n", fgkNBINSE, fgkNBINSZ);
733 // store information in arrays, to increase processing speed
734 Double_t measuredCopy[kMaxME+1][kMaxMZ+1];
735 Double_t prior[kMaxTE+1][kMaxTZ+1];
736 Double_t errors[kMaxTE+1][kMaxTZ+1];
737 Double_t result[kMaxTE+1][kMaxTZ+1];
739 THnSparseF *inverseCorrelation;
740 inverseCorrelation = (THnSparseF*)correlation->Clone("inverseCorrelation");
741 inverseCorrelation->Reset();
743 Float_t inputDistIntegral = 1;
744 if (initialConditions)
746 printf("Using different starting conditions...\n");
747 inputDistIntegral = initialConditions->Integral();
749 Float_t measuredIntegral = measured->Integral();
750 for (Int_t me=1; me<=kMaxME; me++)
751 for (Int_t mz=1; mz<=kMaxMZ; mz++)
753 // normalization of the measured spectrum
754 measuredCopy[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral;
755 errors[me][mz] = measured->GetBinError(me, mz) / measuredIntegral;
756 // pick prior distribution and normalize it
757 if (initialConditions)
758 prior[me][mz] = initialConditions->GetBinContent(me,mz) / inputDistIntegral;
760 prior[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral;
764 for (Int_t i=0; i<nIterations; i++)
766 // calculate Inverse Correlation Map from Bayes theorem:
767 // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
769 for (Int_t me=1; me<=kMaxME; me++)
770 for (Int_t mz=1; mz<=kMaxMZ; mz++)
773 for (Int_t te=kStartBin; te<=kMaxTE; te++)
774 for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++)
776 Int_t bin[4] = {te, me, tz, mz};
777 norm += correlation->GetBinContent(bin)*prior[te][tz];
780 for (Int_t te = kStartBin; te <= kMaxTE; te++)
781 for (Int_t tz = kStartBin; tz <= kMaxTZ; tz++)
783 Int_t bin[4] = {te, me, tz, mz};
784 inverseCorrelation->SetBinContent(bin, correlation->GetBinContent(bin)*prior[te][tz]/norm );
787 // inverse response set to '0' wich has been already done in line 2069
789 inverseCorrelation->Reset();
790 Float_t norm[kMaxTE+2][kMaxTZ+2];
791 for (Int_t te=0; te<(kMaxTE+2); te++)
792 for (Int_t tz=0; tz<(kMaxTZ+2); tz++)
794 for (Int_t idx=0; idx<=correlation->GetNbins(); idx++)
797 Float_t binContent = correlation->GetBinContent(idx, bin);
798 if (bin[1]>0 && bin[1]<=fgkNBINSE && bin[3]>0 && bin[3]<=fgkNBINSZ &&
799 bin[0]>0 && bin[0]<=fgkNBINSE && bin[2]>0 && bin[2]<=fgkNBINSZ)
800 norm[bin[1]][bin[3]] += binContent*prior[bin[0]][bin[2]];
802 Float_t chi2Measured=0, diff;
803 for (Int_t idx=0; idx<=correlation->GetNbins(); idx++)
806 Float_t binContent = correlation->GetBinContent(idx, bin);
807 if (norm[bin[1]][bin[3]]>0 && bin[1]>0 && bin[1]<=fgkNBINSE &&
808 bin[3]>0 && bin[3]<=fgkNBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=fgkNBINSE && bin[2]<=fgkNBINSZ)
810 inverseCorrelation->SetBinContent(bin, binContent*prior[bin[0]][bin[2]]/norm[bin[1]][bin[3]]);
811 if (errors[bin[1]][bin[3]]>0)
813 diff = ((measuredCopy[bin[1]][bin[3]]-norm[bin[1]][bin[3]])/(errors[bin[1]][bin[3]]));
814 chi2Measured += diff*diff;
819 // calculate "generated" spectrum
820 for (Int_t te = kStartBin; te<=kMaxTE; te++)
821 for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++)
824 for (Int_t me=1; me<=kMaxME; me++)
825 for (Int_t mz=1; mz<=kMaxMZ; mz++)
827 Int_t bin[4] = {te, me, tz, mz};
828 value += inverseCorrelation->GetBinContent(bin)*measuredCopy[me][mz];
830 result[te][tz] = value;
831 //printf("%e\n", result[te][tz]);
834 // regularization (simple smoothing)
835 Float_t chi2LastIter = 0;
836 for (Int_t te=kStartBin; te<=kMaxTE; te++)
837 for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++)
839 Float_t newValue = 0;
840 // 0 bin excluded from smoothing
841 if (( te >(kStartBin+1) && te<(kMaxTE-1) ) && ( tz > (kStartBin+1) && tz<(kMaxTZ-1) ))
843 Float_t average = ((result[te-1][tz-1] + result[te-1][tz] + result[te-1][tz+1])+(result[te][tz-1] + result[te][tz] + result[te][tz+1])+(result[te+1][tz-1] + result[te+1][tz] + result[te+1][tz+1]))/9.;
845 // weight the average with the regularization parameter
846 newValue = (1 - regPar) * result[te][tz] + regPar * average;
849 newValue = result[te][tz];
850 if (prior[te][tz]>1.e-5)
852 diff = ((prior[te][tz]-newValue)/prior[te][tz]);
853 chi2LastIter = diff*diff;
855 prior[te][tz] = newValue;
857 //printf(" iteration %d - chi2LastIter = %e - chi2Measured = %e \n", i, chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ), chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ));
858 if (chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-6 && chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-3)
860 } // end of iterations
862 // propagate errors of the reconstructed distribution through the unfolding
863 for (Int_t te = kStartBin; te<=kMaxTE; te++)
864 for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++)
866 Float_t valueError = 0;
867 // Float_t binError = 0;
868 for (Int_t me=1; me<=kMaxME; me++)
869 for (Int_t mz=1; mz<=kMaxMZ; mz++)
871 Int_t bin[4] = {te, me, tz, mz};
872 valueError += inverseCorrelation->GetBinContent(bin)*inverseCorrelation->GetBinContent(bin)*errors[me][mz]*errors[me][mz];
874 //if (errors[te][tz]!=0)printf("errors[%d][%d]=%e\n", te, tz, valueError);
875 aResult->SetBinContent(te, tz, prior[te][tz]);
876 aResult->SetBinError(te, tz, TMath::Sqrt(valueError));
879 // ***********************************************************************************************************
880 // Calculate the covariance matrix, all arguments are taken from G. D'Agostini (p.6-8)
883 printf("Covariance matrix will be calculated... this will take a lot of time (>1 day) ;)\n");
885 //Variables and Matrices that will be use along the calculation
886 const Int_t binsV[4] = {fgkNBINSE,fgkNBINSE, fgkNBINSZ, fgkNBINSZ};
887 const Double_t lowEdgeV[4] = {fgkaxisLowerLimitE, fgkaxisLowerLimitE, fgkaxisLowerLimitZ, fgkaxisLowerLimitZ};
888 const Double_t upEdgeV[4] = {fgkaxisUpperLimitE, fgkaxisUpperLimitE, fgkaxisUpperLimitZ, fgkaxisUpperLimitZ};
890 const Double_t nTrue = (Double_t)measured->Integral();
892 THnSparseF *v = new THnSparseF("V","",4, binsV, lowEdgeV, upEdgeV);
894 Double_t invCorrContent1, nt;
895 Double_t invCorrContent2, v11,v12, v2;
896 // calculate V1 and V2
897 for (Int_t idx1=0; idx1<=nFillesBins; idx1++)
899 printf("Covariance Matrix calculation: iteration idx1=%d of %d\n", idx1, nFillesBins);
900 for (Int_t idx2=0; idx2<=nFillesBins; idx2++)
904 invCorrContent1 = inverseCorrelation->GetBinContent(idx1, bin1);
905 invCorrContent2 = inverseCorrelation->GetBinContent(idx2, bin2);
907 if(bin1[0]>0 && bin1[0]<=fgkNBINSE && bin1[1]>0 && bin1[1]<=fgkNBINSE &&
908 bin1[2]>0 && bin1[2]<=fgkNBINSZ && bin1[3]>0 && bin1[3]<=fgkNBINSZ &&
909 bin2[0]>0 && bin2[0]<=fgkNBINSE && bin2[1]>0 && bin2[1]<=fgkNBINSE &&
910 bin2[2]>0 && bin2[2]<=fgkNBINSZ && bin2[3]>0 && bin2[3]<=fgkNBINSZ)
912 if (bin1[1]==bin2[1] && bin1[3]==bin2[3])
913 v11 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]
914 *(1. - measuredCopy[bin2[1]][bin2[3]]/nTrue);
916 v12 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]*
917 measuredCopy[bin2[1]][bin2[3]]/nTrue;
918 nt = (Double_t)prior[bin2[0]][bin2[2]];
919 v2 = measuredCopy[bin1[1]][bin1[3]]*measuredCopy[bin2[1]][bin2[3]]*
920 invCorrContent1*invCorrContent2*
921 BayesUncertaintyTerms(inverseCorrelation, correlation, bin1, bin2, nt);
922 Int_t binV[4] = {bin1[0],bin2[0],bin1[2],bin2[2]};
923 v->SetBinContent(binV,v11-v12 + v2);
928 for(Int_t te = 1; te<=fgkNBINSE; te++)
929 for(Int_t tz = 1; tz<=fgkNBINSZ; tz++)
931 Int_t binV[4] = {te,te,tz,tz};
932 aResult->SetBinError( te, tz, v->GetBinContent(binV) );
935 TFile* f = new TFile("Covariance_UnfSpectrum.root");
945 //____________________________________________________________________
946 Double_t AliJetSpectrumUnfolding::BayesUncertaintyTerms(THnSparseF* const M, THnSparseF* const C,const Int_t* const binTM, const Int_t* const binTM1, Double_t nt)
949 // helper function for the covariance matrix of the bayesian method
953 Float_t term[9] = {0.};
954 Int_t tmpBin[4] = {0}, tmpBin1[4] = {0};
955 const Int_t nFilledBins = C->GetNbins();
960 Float_t invCorrContent;
962 tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3];
963 tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3];
964 if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
966 if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
967 term[0] = BayesCov(M, C, tmpBin, tmpBin1)/
968 (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
969 term[2] = term[0]*M->GetBinContent(tmpBin1);
977 tmpBin[0]=binTM1[0]; tmpBin[1]=binTM[1]; tmpBin[2]=binTM1[2]; tmpBin[3]=binTM[3];
978 tmpBin1[0]=binTM1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM1[2]; tmpBin1[3]=binTM1[3];
979 if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
980 term[6] = BayesCov(M, C, tmpBin, tmpBin1)*
981 M->GetBinContent(tmpBin)/
982 (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
986 for(Int_t idx1=0; idx1<=nFilledBins; idx1++)
989 corrContent = C->GetBinContent(idx1, bin1);
990 invCorrContent = M->GetBinContent(idx1, bin1);
991 if(bin1[0]>0 && bin1[0]<=fgkNBINSE && bin1[1]>0 && bin1[1]<=fgkNBINSE &&
992 bin1[2]>0 && bin1[2]<=fgkNBINSZ && bin1[3]>0 && bin1[3]<=fgkNBINSZ)
994 tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3];
995 tmpBin1[0]=binTM[0]; tmpBin1[1]=bin1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=bin1[3];
996 if (C->GetBinContent(tmpBin)!=0 &&
997 binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
998 term[1] = BayesCov(M, C, tmpBin, tmpBin1)/C->GetBinContent(tmpBin);
1002 tmpBin[0] =binTM[0]; tmpBin[1] =bin1[1]; tmpBin[2] =binTM[2]; tmpBin[3] =bin1[3];
1003 tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3];
1004 if (C->GetBinContent(tmpBin1)!=0)
1006 if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
1007 term[3] = BayesCov(M, C, tmpBin, tmpBin1)/
1008 C->GetBinContent(tmpBin1);
1009 term[5] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin1)/
1010 C->GetBinContent(tmpBin1);
1018 tmpBin[0] =binTM1[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM1[2]; tmpBin[3] =binTM[3];
1019 tmpBin1[0]=binTM1[0]; tmpBin1[1]=bin1[1]; tmpBin1[2]=binTM1[2]; tmpBin1[3]=bin1[3];
1020 if (C->GetBinContent(tmpBin)!=0)
1021 term[7] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin)/
1022 C->GetBinContent(tmpBin);
1026 tmpBin[0] =bin1[0]; tmpBin[1] =binTM[1]; tmpBin[2] =bin1[2]; tmpBin[3] =binTM[3];
1027 tmpBin1[0]=bin1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=bin1[2]; tmpBin1[3]=binTM1[3];
1028 if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
1029 term[8] = BayesCov(M, C, tmpBin, tmpBin1)*
1030 M->GetBinContent(tmpBin)*M->GetBinContent(tmpBin)/
1031 (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
1035 for (Int_t i=0; i<9; i++)
1036 result += term[i]/nt;
1043 //____________________________________________________________________
1044 Double_t AliJetSpectrumUnfolding::BayesCov(THnSparseF* const M, THnSparseF* const correlation,const Int_t* const binTM,const Int_t* const bin1)
1048 // get the covariance matrix
1052 Double_t result, result1, result2, result3;
1054 if (binTM[0]==bin1[0] && binTM[2]==bin1[2])
1056 if (correlation->GetBinContent(bin1)!=0)
1057 result1 = 1./correlation->GetBinContent(bin1);
1068 if (binTM[1]==bin1[1] && binTM[3]==bin1[3])
1070 Int_t tmpbin[4] = {bin1[0], binTM[1], bin1[2], binTM[3]};
1071 if(correlation->GetBinContent(tmpbin)!=0)
1072 result3 = M->GetBinContent(tmpbin)/correlation->GetBinContent(tmpbin);
1082 return result = result1 + result2 + result3;
1085 //____________________________________________________________________
1086 TH2F* AliJetSpectrumUnfolding::CalculateRecSpectrum(TH2* const inputGen)
1088 // runs the distribution given in inputGen through the correlation histogram identified by
1089 // fCorrelation and produces a reconstructed spectrum
1094 // normalize to convert number of events into probability
1095 /*for (Int_t te=1; te<=NBINSE; te++)
1096 for (Int_t tz=1; tz<=NBINSZ; tz++)
1100 for (Int_t me = 1; me<=NBINSE; me++)
1101 for (Int_t mz = 1; mz<=NBINSZ; mz++)
1103 bin[0] = te; bin[1] = me;
1104 bin[2] = tz; bin[3] = mz;
1105 sum += fCorrelation[correlationMap]->GetBinContent(bin);
1108 for (Int_t me = 1; me<=NBINSE; me++)
1109 for (Int_t mz = 1; mz<=NBINSZ; mz++)
1111 bin[0] = te; bin[1] = me;
1112 bin[2] = tz; bin[3] = mz;
1113 fCorrelation[correlationMap]->SetBinContent(bin, fCorrelation[correlationMap]->GetBinContent(bin)/sum);
1114 fCorrelation[correlationMap]->SetBinError(bin, fCorrelation[correlationMap]->GetBinError(bin)/sum);
1117 // normalize to convert number of events into probability (the following loop is much faster)
1118 Float_t sum[fgkNBINSE+2][fgkNBINSZ+2];
1119 memset(sum,0,sizeof(Float_t)*(fgkNBINSE+2)*(fgkNBINSZ+2));
1121 for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++)
1124 Float_t binContent = fCorrelation->GetBinContent(idx, bin);
1125 if (bin[1]>0 && bin[1]<=fgkNBINSE && bin[3]>0 && bin[3]<=fgkNBINSZ){
1126 sum[bin[0]][bin[2]] += binContent;
1130 for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++)
1133 Float_t binContent = fCorrelation->GetBinContent(idx, bin);
1134 Float_t binError = fCorrelation->GetBinError(bin);
1135 if (sum[bin[0]][bin[2]]>0 && bin[1]>0 && bin[1]<=fgkNBINSE &&
1136 bin[3]>0 && bin[3]<=fgkNBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=fgkNBINSE && bin[2]<=fgkNBINSZ)
1138 fCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]);
1139 fCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]);
1143 TH2F* target = dynamic_cast<TH2F*> (fRecSpectrum->Clone(Form("reconstructed_%s", inputGen->GetName())));
1144 if(!target)return 0;
1148 for (Int_t me=1; me<=fgkNBINSE; ++me)
1149 for (Int_t mz=1; mz<=fgkNBINSZ; ++mz)
1151 Float_t measured = 0;
1154 for (Int_t te=1; te<=fgkNBINSE; ++te)
1155 for (Int_t tz=1; tz<=fgkNBINSZ; ++tz)
1157 Int_t bin[4] = {te, me, tz, mz};
1158 measured += inputGen->GetBinContent(te,tz) * fCorrelation->GetBinContent(bin);
1159 error += inputGen->GetBinError(te,tz) * fCorrelation->GetBinContent(bin);
1161 target->SetBinContent(me, mz, measured);
1162 target->SetBinError(me, mz, error);
1168 //__________________________________________________________________________________________________
1169 void AliJetSpectrumUnfolding::SetGenRecFromFunc(const TF2* const inputGen)
1171 // uses the given function to fill the input Generated histogram and generates from that
1172 // the reconstructed histogram by applying the response histogram
1173 // this can be used to evaluate if the methods work indepedently of the input
1179 TH2F* histtmp = new TH2F("histtmp", "tmp",
1180 fgkNBINSE, fgkaxisLowerLimitE, fgkaxisUpperLimitE,
1181 fgkNBINSZ, fgkaxisLowerLimitZ, fgkaxisUpperLimitZ);
1183 TH2F* gen = fGenSpectrum;
1188 histtmp->FillRandom(inputGen->GetName(), fgkNEVENTS);
1190 for (Int_t i=1; i<=gen->GetNbinsX(); ++i)
1191 for (Int_t j=1; j<=gen->GetNbinsY(); ++j)
1193 gen->SetBinContent(i, j, histtmp->GetBinContent(i,j));
1194 gen->SetBinError(i, j, histtmp->GetBinError(i,j));
1200 //gStyle->SetPalette(1);
1202 //gen->Draw("COLZ");
1205 TH2 *recsave = fRecSpectrum;
1207 fRecSpectrum = CalculateRecSpectrum(gen);
1208 fRecSpectrum->SetName(recsave->GetName());
1213 //________________________________________________________________________________________