git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adding the macros and the class for 2-dim unfolding (Ydalia)
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliJetSpectrumUnfolding.cxx
CommitLineData
6d75bdb8 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15//
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
18//
19
20#include "AliJetSpectrumUnfolding.h"
21
22#include <TFile.h>
23#include <TH1F.h>
24#include <TH2F.h>
25#include <TH3F.h>
26#include <TDirectory.h>
27#include <TVirtualFitter.h>
28#include <TCanvas.h>
29#include <TString.h>
30#include <TF1.h>
31#include <TF2.h>
32#include <TMath.h>
33#include <TCollection.h>
34#include <TLegend.h>
35#include <TLine.h>
36#include <TRandom.h>
37#include <TProfile.h>
38#include <TProfile2D.h>
39#include <TStyle.h>
40#include <TColor.h>
41
42#include <THnSparse.h>
43
44ClassImp(AliJetSpectrumUnfolding)
45
46const Int_t NBINSE = 50;
47const Int_t NBINSZ = 50;
48const Int_t NEVENTS = 500000;
49const Double_t axisLowerLimitE = 0.;
50const Double_t axisLowerLimitZ = 0.;
51const Double_t axisUpperLimitE = 250.;
52const Double_t axisUpperLimitZ = 1.;
53
54Float_t AliJetSpectrumUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing)
55Int_t AliJetSpectrumUnfolding::fgBayesianIterations = 100; // number of iterations in Bayesian method
56
57//____________________________________________________________________
58
59AliJetSpectrumUnfolding::AliJetSpectrumUnfolding() :
60 TNamed(), fCurrentRec(0), fCurrentCorrelation(0), fGenSpectrum(0),
61 fRecSpectrum(0), fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
62{
63 //
64 // default constructor
65 //
66
67 fGenSpectrum = 0;
68 fRecSpectrum = 0;
69 fUnfSpectrum = 0;
70 fCorrelation = 0;
71}
72
73//____________________________________________________________________
74AliJetSpectrumUnfolding::AliJetSpectrumUnfolding(const Char_t* name, const Char_t* title) :
75 TNamed(name, title), fCurrentRec(0), fCurrentCorrelation(0), fGenSpectrum(0),
76 fRecSpectrum(0), fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
77{
78 //
79 // named constructor
80 //
81
82 // do not add this hists to the directory
83 Bool_t oldStatus = TH1::AddDirectoryStatus();
84 TH1::AddDirectory(kFALSE);
85
86 #define ZBINNING NBINSZ, axisLowerLimitZ, axisUpperLimitZ // fragmentation of leading particle
87 #define EBINNING NBINSE, axisLowerLimitE, axisUpperLimitE // energy of the jet
88
89 fRecSpectrum = new TH2F("fRecSpectrum", "Reconstructed Spectrum;E^{jet}_{rec} [GeV];z^{lp}_{rec}", EBINNING, ZBINNING);
90 fGenSpectrum = new TH2F("fGenSpectrum", "Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}_{gen}", EBINNING, ZBINNING);
91 fUnfSpectrum = new TH2F("fUnfSpectrum", "Unfolded Spectrum;E^{jet} [GeV];z^{lp}", EBINNING, ZBINNING);
92
93 const Int_t nbin[4]={NBINSE, NBINSE, NBINSZ, NBINSZ};
94 //arrays for bin limits
95 Double_t LowEdge[4] = {axisLowerLimitE, axisLowerLimitE, axisLowerLimitZ, axisLowerLimitZ};
96 Double_t UpEdge[4] = {axisUpperLimitE, axisUpperLimitE, axisUpperLimitZ, axisUpperLimitZ};
97
98 fCorrelation = new THnSparseF("fCorrelation", "Correlation Function", 4, nbin, LowEdge, UpEdge);
99
100 TH1::AddDirectory(oldStatus);
101}
102
103//____________________________________________________________________
104AliJetSpectrumUnfolding::~AliJetSpectrumUnfolding()
105{
106 //
107 // Destructor
108 //
109
110 if (fGenSpectrum)
111 delete fGenSpectrum;
112 fGenSpectrum = 0;
113
114 if (fRecSpectrum)
115 delete fRecSpectrum;
116 fRecSpectrum = 0;
117
118 if (fUnfSpectrum)
119 delete fUnfSpectrum;
120 fUnfSpectrum = 0;
121
122 if (fCorrelation)
123 delete fCorrelation;
124 fCorrelation = 0;
125
126}
127
128//____________________________________________________________________
129Long64_t AliJetSpectrumUnfolding::Merge(TCollection* list)
130{
131 // Merge a list of AliJetSpectrumUnfolding objects with this (needed for
132 // PROOF).
133 // Returns the number of merged objects (including this).
134
135 if (!list)
136 return 0;
137
138 if (list->IsEmpty())
139 return 1;
140
141 TIterator* iter = list->MakeIterator();
142 TObject* obj;
143
144 // collections of all histograms
145 TList collections[4];
146
147 Int_t count = 0;
148 while ((obj = iter->Next())) {
149
150 AliJetSpectrumUnfolding* entry = dynamic_cast<AliJetSpectrumUnfolding*> (obj);
151 if (entry == 0)
152 continue;
153
154 collections[0].Add(entry->fGenSpectrum);
155 collections[1].Add(entry->fRecSpectrum);
156 collections[2].Add(entry->fUnfSpectrum);
157 collections[3].Add(entry->fCorrelation);
158
159 count++;
160 }
161
162 fGenSpectrum->Merge(&collections[0]);
163 fRecSpectrum->Merge(&collections[1]);
164 fUnfSpectrum->Merge(&collections[2]);
165 fCorrelation->Merge(&collections[3]);
166
167 delete iter;
168
169 return count+1;
170}
171
172//____________________________________________________________________
173Bool_t AliJetSpectrumUnfolding::LoadHistograms(const Char_t* dir)
174{
175 //
176 // loads the histograms from a file
177 // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
178 //
179
180 if (!dir)
181 dir = GetName();
182
183 if (!gDirectory->cd(dir))
184 return kFALSE;
185
186 Bool_t success = kTRUE;
187
188 // store old histograms to delete them later
189 TList oldHistograms;
190 oldHistograms.SetOwner(1);
191
192 if (fGenSpectrum) oldHistograms.Add(fGenSpectrum);
193 if (fRecSpectrum) oldHistograms.Add(fRecSpectrum);
194 if (fUnfSpectrum) oldHistograms.Add(fUnfSpectrum);
195 if (fCorrelation) oldHistograms.Add(fCorrelation);
196
197 // load new histograms
198 fGenSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fGenSpectrum->GetName()));
199 if (!fGenSpectrum)
200 success = kFALSE;
201
202 fRecSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fRecSpectrum->GetName()));
203 if (!fRecSpectrum)
204 success = kFALSE;
205
206 fUnfSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fUnfSpectrum->GetName()));
207 if (!fUnfSpectrum)
208 success = kFALSE;
209
210 fCorrelation = dynamic_cast<THnSparseF*> (gDirectory->Get(fCorrelation->GetName()));
211 if (!fCorrelation)
212 success = kFALSE;
213
214 gDirectory->cd("..");
215
216 // delete old histograms
217 oldHistograms.Delete();
218
219 return success;
220}
221
222//____________________________________________________________________
223void AliJetSpectrumUnfolding::SaveHistograms()
224{
225 //
226 // saves the histograms
227 //
228
229 gDirectory->mkdir(GetName());
230 gDirectory->cd(GetName());
231
232 if (fGenSpectrum)
233 fGenSpectrum->Write();
234
235 if (fRecSpectrum)
236 fRecSpectrum->Write();
237
238 if (fUnfSpectrum)
239 fUnfSpectrum->Write();
240
241 if (fCorrelation)
242 fCorrelation->Write();
243
244 gDirectory->cd("..");
245}
246
247//____________________________________________________________________
248void AliJetSpectrumUnfolding::SetupCurrentHists(Bool_t createBigBin)
249{
250 //
251 // resets fUnfSpectrum
252 //
253
254 fUnfSpectrum->Reset();
255 fUnfSpectrum->Sumw2();
256
257 fCurrentRec = (TH2F*)fRecSpectrum->Clone("fCurrentRec");
258 fCurrentRec->Sumw2();
259
260 fCurrentCorrelation = (THnSparseF*)fCorrelation->Clone("fCurrentCorrelation");
261 fCurrentCorrelation->Sumw2();
262
263 if (createBigBin)
264 {
265 Int_t maxBinE = 0, maxBinZ = 0;
266 Float_t maxE = 0, maxZ = 0;
267 for (Int_t me=1; me<=fCurrentRec->GetNbinsX(); me++)
268 for (Int_t mz=1; mz<=fCurrentRec->GetNbinsY(); mz++)
269 {
270 if (fCurrentRec->GetBinContent(me,mz) <= 5 && me>NBINSE/2 && mz>NBINSZ/2)
271 {
272 maxBinE = me;
273 maxBinZ = mz;
274 maxE = fCurrentRec->GetXaxis()->GetBinCenter(me);
275 maxZ = fCurrentRec->GetYaxis()->GetBinCenter(mz);
276 break;
277 }
278 }
279
280 if (maxBinE > 0 || maxBinZ > 0)
281 {
282 printf("Bin limit in measured spectrum is e = %d and z = %d.\n", maxBinE, maxBinZ);
283 fCurrentRec->SetBinContent(maxBinE, maxBinZ, fCurrentRec->Integral(maxBinE, fCurrentRec->GetNbinsX(), maxBinZ, fCurrentRec->GetNbinsY()));
284 for (Int_t me=maxBinE+1; me<=fCurrentRec->GetNbinsX(); me++)
285 for (Int_t mz=maxBinZ+1; mz<=fCurrentRec->GetNbinsY(); mz++)
286 {
287 fCurrentRec->SetBinContent(me, mz, 0);
288 fCurrentRec->SetBinError(me, mz, 0);
289 }
290 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
291 fCurrentRec->SetBinError(maxBinE, maxBinZ, TMath::Sqrt(fCurrentRec->GetBinContent(maxBinE, maxBinZ)));
292
293 printf("This bin has now %f +- %f entries\n", fCurrentRec->GetBinContent(maxBinE, maxBinZ), fCurrentRec->GetBinError(maxBinE, maxBinZ));
294
295 /* for (Int_t te=1; te<=NBINSE; te++)
296 {
297 for (Int_t tz=1; tz<=NBINSZ; tz++)
298 {
299 Int_t binMin[4] = {te, maxBinE, tz, maxBinZ};
300 Int_t binMax[4] = {NBINSE, NBINSE, NBINSZ, NBINSZ};
301 Float_t sum=0;
302 for (Int_t ite=te; ite<=NBINSE; ite++)
303 for (Int_t itz=tz; itz<=NBINSZ; itz++)
304 for (Int_t ime=maxBinE; ime<=NBINSE; ime++)
305 for (Int_t imz=maxBinZ; imz<=NBINSZ; imz++)
306 {
307 Int_t bin[4] = {ite, ime, itz, imz};
308 sum += fCurrentCorrelation->GetBinContent(bin);
309 }
310 fCurrentCorrelation->SetBinContent(binMin, sum);
311 fCurrentCorrelation->SetBinError(binMin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin)));
312 printf("create big bin1, nbins = %d, te = %d, tz = %d\n", NBINSE, te, tz);
313 for (Int_t me=maxBinE; me<=NBINSE; me++)
314 {
315 for (Int_t mz=maxBinZ; mz<=NBINSZ; mz++)
316 {
317 Int_t bin[4] = {te, me, tz, mz};
318 fCurrentCorrelation->SetBinContent(bin, 0.);
319 fCurrentCorrelation->SetBinError(bin, 0.);
320 printf("create big bin2\n");
321 }
322 }
323 }
324 }*/
325
326 for(Int_t idx = 0; idx<=fCurrentCorrelation->GetNbins(); idx++)
327 {
328 Int_t bin[4];
329 Float_t binContent = fCurrentCorrelation->GetBinContent(idx,bin);
330 Float_t binError = fCurrentCorrelation->GetBinError(idx);
331 Int_t binMin[4] = {bin[0], maxBinE, bin[2], maxBinZ};
332 if ( (bin[1]>maxBinE && bin[1]<=NBINSE) && (bin[3]>maxBinZ && bin[3]<=NBINSZ) )
333 {
334 fCurrentCorrelation->SetBinContent(binMin, binContent + fCurrentCorrelation->GetBinContent(binMin));
335 fCurrentCorrelation->SetBinError(binMin, binError + TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin)));
336 fCurrentCorrelation->SetBinContent(bin, 0.);
337 fCurrentCorrelation->SetBinError(bin, 0.);
338 }
339 printf("create big bin1, nbins = %d, te = %d, tz = %d\n", NBINSE, bin[0], bin[1]);
340 }
341
342 printf("Adjusted correlation matrix!\n");
343 }
344 } // end Create Big Bin
345
346}
347
348//____________________________________________________________________
349void AliJetSpectrumUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
350{
351 //
352 // sets the parameters for Bayesian unfolding
353 //
354
355 fgBayesianSmoothing = smoothing;
356 fgBayesianIterations = nIterations;
357
358 printf("AliJetSpectrumUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
359}
360
361//____________________________________________________________________
362void AliJetSpectrumUnfolding::NormalizeToBinWidth(TH2* hist)
363{
364 //
365 // normalizes a 2-d histogram to its bin width (x width * y width)
366 //
367
368 for (Int_t i=1; i<=hist->GetNbinsX(); i++)
369 for (Int_t j=1; j<=hist->GetNbinsY(); j++)
370 {
371 Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
372 hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
373 hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
374 }
375}
376
377//____________________________________________________________________
378void AliJetSpectrumUnfolding::DrawHistograms()
379{
380 //
381 // draws the histograms of this class
382 //
383
384 gStyle->SetPalette(1);
385
386 TCanvas* canvas1 = new TCanvas("fRecSpectrum", "fRecSpectrum", 900, 600);
387 gPad->SetLogz();
388 fRecSpectrum->DrawCopy("COLZ");
389
390 TCanvas* canvas2 = new TCanvas("fGenSpectrum", "fGenSpectrum", 900, 600);
391 gPad->SetLogz();
392 fGenSpectrum->DrawCopy("COLZ");
393
394 TCanvas* canvas3 = new TCanvas("fUnfSpectrum", "fUnfSpectrum", 900, 600);
395 gPad->SetLogz();
396 fUnfSpectrum->DrawCopy("COLZ");
397
398 TCanvas* canvas4 = new TCanvas("fCorrelation", "fCorrelation", 500, 500);
399 canvas1->Divide(2);
400
401 canvas4->cd(1);
402 gPad->SetLogz();
403 TH2D* h0 = fCorrelation->Projection(1,0);
404 h0->SetXTitle("E^{jet}_{gen} [GeV]");
405 h0->SetYTitle("E^{jet}_{rec} [GeV]");
406 h0->SetTitle("Projection: Jet Energy");
407 h0->DrawCopy("colz");
408
409 canvas1->cd(2);
410 gPad->SetLogz();
411 TH2D* h1 = fCorrelation->Projection(3,2);
412 h1->SetXTitle("z^{lp}_{gen}");
413 h1->SetYTitle("z^{lp}_{rec}");
414 h1->SetTitle("Projection: Leading Particle Fragmentation");
415 h1->DrawCopy("colz");
416
417}
418
419//____________________________________________________________________
420void AliJetSpectrumUnfolding::DrawComparison(const char* name, TH2* genHist)
421{
422
423 if (fUnfSpectrum->Integral() == 0)
424 {
425 printf("ERROR. Unfolded histogram is empty\n");
426 return;
427 }
428
429 //regain measured distribution used for unfolding, because the bins were modified in SetupCurrentHists
430 //in create big bin
431 fCurrentRec = (TH2F*)fRecSpectrum->Clone();
432 fCurrentRec->Sumw2();
433 fCurrentRec->Scale(1.0 / fCurrentRec->Integral());
434
435 // normalize unfolded result to 1
436 fUnfSpectrum->Scale(1.0 / fUnfSpectrum->Integral());
437
438 // find bin with <= 100 entries. this is used as limit for the chi2 calculation
439 Int_t mcBinLimitE = 0, mcBinLimitZ = 0;
440 for (Int_t i=0; i<genHist->GetNbinsX(); ++i)
441 for (Int_t j=0; j<genHist->GetNbinsY(); ++j)
442 {
443 if (genHist->GetBinContent(i,j) > 100)
444 {
445 mcBinLimitE = i;
446 mcBinLimitZ = j;
447 }
448 else
449 break;
450 }
451 Printf("AliJetSpectrumUnfolding::DrawComparison: Gen bin limit is (x,y) = (%d,%d)", mcBinLimitE,mcBinLimitZ);
452
453 // scale to 1 true spectrum
454 genHist->Sumw2();
455 genHist->Scale(1.0 / genHist->Integral());
456
457 // calculate residual
458 // for that we convolute the response matrix with the gathered result
459 TH2* tmpRecRecorrected = (TH2*) fUnfSpectrum->Clone("tmpRecRecorrected");
460 TH2* convoluted = CalculateRecSpectrum(tmpRecRecorrected);
461 if (convoluted->Integral() > 0)
462 convoluted->Scale(1.0 / convoluted->Integral());
463 else
464 printf("ERROR: convoluted is empty. Something went wrong calculating the convoluted histogram.\n");
465
466 TH2* residual = (TH2*) convoluted->Clone("residual");
467 residual->SetTitle("(R#otimesUnfolded - Reconstructed)/Reconstructed;E^{jet} [GeV]; z^{lp}");
468
469 fCurrentRec->Scale(1./fCurrentRec->Integral());
470 residual->Add(fCurrentRec, -1);
471 //residual->Divide(residual, fCurrentRec, 1, 1, "B");
472
473 // draw canvas
474 TCanvas* canvas1 = new TCanvas(name, name, 1000, 1000);
475 canvas1->Divide(2, 3);
476
477 Int_t style = 1;
478 const Int_t NRGBs = 5;
479 const Int_t NCont = 500;
480
481 Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
482 Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
483 Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
484 Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
485 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
486 gStyle->SetNumberContours(NCont);
487
488 canvas1->cd(1);
489 gStyle->SetPalette(style);
490 gPad->SetLogz();
491 genHist->SetTitle("Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}");
492 genHist->SetStats(0);
493 genHist->DrawCopy("colz");
494
495 canvas1->cd(2);
496 gStyle->SetPalette(style);
497 gPad->SetLogz();
498 fUnfSpectrum->SetStats(0);
499 fUnfSpectrum->DrawCopy("colz");
500
501 canvas1->cd(3);
502 gStyle->SetPalette(style);
503 gPad->SetLogz();
504 fCurrentRec->SetTitle(fRecSpectrum->GetTitle());
505 fCurrentRec->SetStats(0);
506 fCurrentRec->DrawCopy("colz");
507
508 canvas1->cd(4);
509 gStyle->SetPalette(style);
510 gPad->SetLogy();
511 TH1D* projGenX = genHist->ProjectionX();
512 projGenX->SetTitle("Projection: Jet Energy; E^{jet} [GeV]");
513 TH1D* projUnfX = fUnfSpectrum->ProjectionX();
514 TH1D* projRecX = fCurrentRec->ProjectionX();
515 projGenX->SetStats(0);
516 projRecX->SetStats(0);
517 projUnfX->SetStats(0);
518 projGenX->SetLineColor(8);
519 projRecX->SetLineColor(2);
520 projGenX->DrawCopy();
521 projUnfX->DrawCopy("same");
522 projRecX->DrawCopy("same");
523
524 TLegend* legend = new TLegend(0.6, 0.85, 0.98, 0.98);
525 legend->AddEntry(projGenX, "Generated Spectrum");
526 legend->AddEntry(projUnfX, "Unfolded Spectrum");
527 legend->AddEntry(projRecX, "Reconstructed Spectrum");
528 //legend->SetFillColor(0);
529 legend->Draw("same");
530
531 canvas1->cd(5);
532 gPad->SetLogy();
533 gStyle->SetPalette(style);
534 TH1D* projGenY = genHist->ProjectionY();
535 projGenY->SetTitle("Projection: Leading Particle Fragmentation; z^{lp}");
536 TH1D* projUnfY = fUnfSpectrum->ProjectionY();
537 TH1D* projRecY = fCurrentRec->ProjectionY();
538 projGenY->SetStats(0);
539 projRecY->SetStats(0);
540 projUnfY->SetStats(0);
541 projGenY->SetLineColor(8);
542 projRecY->SetLineColor(2);
543 projGenY->DrawCopy();
544 projUnfY->DrawCopy("same");
545 projRecY->DrawCopy("same");
546
547 TLegend* legend1 = new TLegend(0.6, 0.85, 0.98, 0.98);
548 legend1->AddEntry(projGenY, "Generated Spectrum");
549 legend1->AddEntry(projUnfY, "Unfolded Spectrum");
550 legend1->AddEntry(projRecY, "Recontructed Spectrum");
551 //legend1->SetFillColor(0);
552 legend1->Draw("same");
553
554 // Draw residuals
555 canvas1->cd(6);
556 gStyle->SetPalette(style);
557 gPad->SetLogz();
558 residual->SetStats(0);
559 residual->DrawCopy("colz");
560
561 canvas1->SaveAs(Form("%s.png", canvas1->GetName()));
562}
563
564
565//____________________________________________________________________
566void AliJetSpectrumUnfolding::GetComparisonResults(Float_t* gen, Int_t* genLimit, Float_t* residuals, Float_t* ratioAverage) const
567{
568 // Returns the chi2 between the Generated and the unfolded Reconstructed spectrum as well as between the Reconstructed and the folded unfolded
569 // These values are computed during DrawComparison, thus this function picks up the
570 // last calculation
571
572 if (gen)
573 *gen = fLastChi2MC;
574 if (genLimit)
575 *genLimit = fLastChi2MCLimit;
576 if (residuals)
577 *residuals = fLastChi2Residuals;
578 if (ratioAverage)
579 *ratioAverage = fRatioAverage;
580}
581
582//____________________________________________________________________
583void AliJetSpectrumUnfolding::ApplyBayesianMethod(Float_t regPar, Int_t nIterations, TH2* initialConditions, Bool_t determineError)
584{
585 //
586 // correct spectrum using bayesian unfolding
587 //
588
589 // initialize seed with current time
590 gRandom->SetSeed(0);
591
592 printf("seting up current arrays and histograms...\n");
593 SetupCurrentHists(kFALSE); // kFALSE to not create big bin
594
595 // normalize Correlation Map to convert number of events into probabilities
596 /*for (Int_t te=1; te<=NBINSE; te++)
597 for (Int_t tz=1; tz<=NBINSZ; tz++)
598 {
599 Int_t bin[4];
600 Float_t sum=0.;
601 for (Int_t me = 1; me<=NBINSE; me++)
602 for (Int_t mz = 1; mz<=NBINSZ; mz++)
603 {
604 bin[0] = te; bin[1] = me;
605 bin[2] = tz; bin[3] = mz;
606 sum += fCurrentCorrelation->GetBinContent(bin);
607 }
608 if (sum > 0.)
609 for (Int_t me = 1; me<=NBINSE; me++)
610 for (Int_t mz = 1; mz<=NBINSZ; mz++)
611 {
612 bin[0] = te; bin[1] = me;
613 bin[2] = tz; bin[3] = mz;
614 fCurrentCorrelation->SetBinContent(bin, fCurrentCorrelation->GetBinContent(bin)/sum);
615 fCurrentCorrelation->SetBinError(bin, fCurrentCorrelation->GetBinError(bin)/sum);
616 }
617 }*/
618 Float_t sum[NBINSE+2][NBINSZ+2];
619 memset(sum,0,sizeof(Float_t)*(NBINSE+2)*(NBINSZ+2));
620
621 for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++)
622 {
623 Int_t bin[4];
624 Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin);
625 if ( (bin[1]>0 && bin[1]<=NBINSE) && (bin[3]>0 && bin[3]<=NBINSZ) )
626 sum[bin[0]][bin[2]] += binContent;
627 }
628
629 for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++)
630 {
631 Int_t bin[4];
632 Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin);
633 Float_t binError = fCurrentCorrelation->GetBinError(bin);
634 if (sum[bin[0]][bin[2]]>0 && (bin[1]>0 && bin[1]<=NBINSE) &&
635 (bin[3]>0 && bin[3]<=NBINSZ) && (bin[0]>0 && bin[0]<=NBINSE) && (bin[2]>0 && bin[2]<=NBINSZ) )
636 {
637 fCurrentCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]);
638 fCurrentCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]);
639 }
640 }
641
642 printf("calling UnfoldWithBayesian\n");
643 Int_t success = UnfoldWithBayesian(fCurrentCorrelation, fCurrentRec, initialConditions, fUnfSpectrum, regPar, nIterations, kFALSE);
644
645 if ( success != 0)
646 return;
647
648 if (!determineError)
649 {
650 Printf("AliJetSpectrumUnfolding::ApplyBayesianMethod: WARNING: No errors calculated.");
651 return;
652 }
653
654 // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
655 // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
656 // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic
657 // error of the unfolding method itself.
658
659 TH2* maxError = (TH2*) fUnfSpectrum->Clone("maxError");
660 maxError->Reset();
661
662 TH2* normalizedResult = (TH2*) fUnfSpectrum->Clone("normalizedResult");
663 normalizedResult->Scale(1.0 / normalizedResult->Integral());
664
665 const Int_t kErrorIterations = 20;
666
667 printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
668
669 TH2* randomized = (TH2*) fCurrentRec->Clone("randomized");
670 TH2* result2 = (TH2*) fUnfSpectrum->Clone("result2");
671 for (Int_t n=0; n<kErrorIterations; ++n)
672 {
673 // randomize the content of clone following a poisson with the mean = the value of that bin
674 for (Int_t x=1; x<=randomized->GetNbinsX(); x++)
675 for (Int_t y=1; y<=randomized->GetNbinsY(); y++)
676 {
677 Float_t randomValue = fCurrentRec->GetBinContent(x,y);
678 TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",randomValue*0.25, randomValue*1.25);
679 poisson->SetParameters(randomValue,0.);
680 randomValue = poisson->GetRandom();
681 //printf("%e --> %e\n", fCurrentRec->GetBinContent(x,y), (Double_t)randomValue);
682 randomized->SetBinContent(x, y, randomValue);
683 delete poisson;
684 }
685
686 result2->Reset();
687 if (UnfoldWithBayesian(fCurrentCorrelation, randomized, initialConditions, result2, regPar, nIterations) != 0)
688 return;
689
690 result2->Scale(1.0 / result2->Integral());
691
692 // calculate ratio
693 result2->Divide(normalizedResult);
694
695 //new TCanvas; result2->DrawCopy("HIST");
696
697 // find max. deviation
698 for (Int_t i=1; i<=result2->GetNbinsX(); i++)
699 for (Int_t j=1; j<=result2->GetNbinsY(); j++)
700 maxError->SetBinContent(i, j, TMath::Max(maxError->GetBinContent(i,j), TMath::Abs(1 - result2->GetBinContent(i,j))));
701 }
702 delete randomized;
703 delete result2;
704
705 for (Int_t i=1; i<=fUnfSpectrum->GetNbinsX(); i++)
706 for (Int_t j=1; j<=fUnfSpectrum->GetNbinsY(); j++)
707 fUnfSpectrum->SetBinError(i, j, fUnfSpectrum->GetBinError(i,j) + maxError->GetBinContent(i,j)*fUnfSpectrum->GetBinContent(i,j));
708
709 delete maxError;
710 delete normalizedResult;
711}
712
713//____________________________________________________________________
714Int_t AliJetSpectrumUnfolding::UnfoldWithBayesian(THnSparseF* correlation, TH2* measured, TH2* initialConditions, TH2* aResult, Float_t regPar, Int_t nIterations, Bool_t calculateErrors)
715{
716 //
717 // unfolds a spectrum
718 //
719
720 if (measured->Integral() <= 0)
721 {
722 Printf("AliJetSpectrumUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
723 return 1;
724 }
725 const Int_t NFilledBins = correlation->GetNbins();
726 const Int_t kStartBin = 1;
727
728 const Int_t kMaxTZ = NBINSZ; // max true axis fragmentation function
729 const Int_t kMaxMZ = NBINSZ; // max measured axis fragmentation function
730 const Int_t kMaxTE = NBINSE; // max true axis energy
731 const Int_t kMaxME = NBINSE; // max measured axis energy
732
733 printf("NbinsE=%d - NbinsZ=%d\n", NBINSE, NBINSZ);
734
735 // store information in arrays, to increase processing speed
736 Double_t measuredCopy[kMaxME+1][kMaxMZ+1];
737 Double_t prior[kMaxTE+1][kMaxTZ+1];
738 Double_t errors[kMaxTE+1][kMaxTZ+1];
739 Double_t result[kMaxTE+1][kMaxTZ+1];
740
741 THnSparseF *inverseCorrelation;
742 inverseCorrelation = (THnSparseF*)correlation->Clone("inverseCorrelation");
743 inverseCorrelation->Reset();
744
745 Float_t inputDistIntegral = 1;
746 if (initialConditions)
747 {
748 printf("Using different starting conditions...\n");
749 inputDistIntegral = initialConditions->Integral();
750 }
751 Float_t measuredIntegral = measured->Integral();
752 for (Int_t me=1; me<=kMaxME; me++)
753 for (Int_t mz=1; mz<=kMaxMZ; mz++)
754 {
755 // normalization of the measured spectrum
756 measuredCopy[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral;
757 errors[me][mz] = measured->GetBinError(me, mz) / measuredIntegral;
758 // pick prior distribution and normalize it
759 if (initialConditions)
760 prior[me][mz] = initialConditions->GetBinContent(me,mz) / inputDistIntegral;
761 else
762 prior[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral;
763 }
764
765 // unfold...
766 for (Int_t i=0; i<nIterations; i++)
767 {
768 // calculate Inverse Correlation Map from Bayes theorem:
769 // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
770 /*Float_t norm = 0;
771 for (Int_t me=1; me<=kMaxME; me++)
772 for (Int_t mz=1; mz<=kMaxMZ; mz++)
773 {
774 norm = 0;
775 for (Int_t te=kStartBin; te<=kMaxTE; te++)
776 for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++)
777 {
778 Int_t bin[4] = {te, me, tz, mz};
779 norm += correlation->GetBinContent(bin)*prior[te][tz];
780 }
781 if (norm > 0)
782 for (Int_t te = kStartBin; te <= kMaxTE; te++)
783 for (Int_t tz = kStartBin; tz <= kMaxTZ; tz++)
784 {
785 Int_t bin[4] = {te, me, tz, mz};
786 inverseCorrelation->SetBinContent(bin, correlation->GetBinContent(bin)*prior[te][tz]/norm );
787 }
788 //else
789 // inverse response set to '0' wich has been already done in line 2069
790 }*/
791 inverseCorrelation->Reset();
792 Float_t norm[kMaxTE+2][kMaxTZ+2];
793 for (Int_t te=0; te<(kMaxTE+2); te++)
794 for (Int_t tz=0; tz<(kMaxTZ+2); tz++)
795 norm[te][tz]=0;
796 for (Int_t idx=0; idx<=correlation->GetNbins(); idx++)
797 {
798 Int_t bin[4];
799 Float_t binContent = correlation->GetBinContent(idx, bin);
800 if (bin[1]>0 && bin[1]<=NBINSE && bin[3]>0 && bin[3]<=NBINSZ &&
801 bin[0]>0 && bin[0]<=NBINSE && bin[2]>0 && bin[2]<=NBINSZ)
802 norm[bin[1]][bin[3]] += binContent*prior[bin[0]][bin[2]];
803 }
804 Float_t chi2Measured=0, diff;
805 for (Int_t idx=0; idx<=correlation->GetNbins(); idx++)
806 {
807 Int_t bin[4];
808 Float_t binContent = correlation->GetBinContent(idx, bin);
809 if (norm[bin[1]][bin[3]]>0 && bin[1]>0 && bin[1]<=NBINSE &&
810 bin[3]>0 && bin[3]<=NBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=NBINSE && bin[2]<=NBINSZ)
811 {
812 inverseCorrelation->SetBinContent(bin, binContent*prior[bin[0]][bin[2]]/norm[bin[1]][bin[3]]);
813 if (errors[bin[1]][bin[3]]>0)
814 {
815 diff = ((measuredCopy[bin[1]][bin[3]]-norm[bin[1]][bin[3]])/(errors[bin[1]][bin[3]]));
816 chi2Measured += diff*diff;
817 }
818 }
819 }
820
821 // calculate "generated" spectrum
822 for (Int_t te = kStartBin; te<=kMaxTE; te++)
823 for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++)
824 {
825 Float_t value = 0;
826 for (Int_t me=1; me<=kMaxME; me++)
827 for (Int_t mz=1; mz<=kMaxMZ; mz++)
828 {
829 Int_t bin[4] = {te, me, tz, mz};
830 value += inverseCorrelation->GetBinContent(bin)*measuredCopy[me][mz];
831 }
832 result[te][tz] = value;
833 //printf("%e\n", result[te][tz]);
834 }
835
836 // regularization (simple smoothing)
837 Float_t chi2LastIter = 0;
838 for (Int_t te=kStartBin; te<=kMaxTE; te++)
839 for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++)
840 {
841 Float_t newValue = 0;
842 // 0 bin excluded from smoothing
843 if (( te >(kStartBin+1) && te<(kMaxTE-1) ) && ( tz > (kStartBin+1) && tz<(kMaxTZ-1) ))
844 {
845 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.;
846
847 // weight the average with the regularization parameter
848 newValue = (1 - regPar) * result[te][tz] + regPar * average;
849 }
850 else
851 newValue = result[te][tz];
852 if (prior[te][tz]>1.e-5)
853 {
854 diff = ((prior[te][tz]-newValue)/prior[te][tz]);
855 chi2LastIter = diff*diff;
856 }
857 prior[te][tz] = newValue;
858 }
859 //printf(" iteration %d - chi2LastIter = %e - chi2Measured = %e \n", i, chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ), chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ));
860 if (chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-6 && chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-3)
861 break;
862 } // end of iterations
863
864 // propagate errors of the reconstructed distribution through the unfolding
865 for (Int_t te = kStartBin; te<=kMaxTE; te++)
866 for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++)
867 {
868 Float_t valueError = 0;
869 Float_t binError = 0;
870 for (Int_t me=1; me<=kMaxME; me++)
871 for (Int_t mz=1; mz<=kMaxMZ; mz++)
872 {
873 Int_t bin[4] = {te, me, tz, mz};
874 valueError += inverseCorrelation->GetBinContent(bin)*inverseCorrelation->GetBinContent(bin)*errors[me][mz]*errors[me][mz];
875 }
876 //if (errors[te][tz]!=0)printf("errors[%d][%d]=%e\n", te, tz, valueError);
877 aResult->SetBinContent(te, tz, prior[te][tz]);
878 aResult->SetBinError(te, tz, TMath::Sqrt(valueError));
879 }
880
881 // ***********************************************************************************************************
882 // Calculate the covariance matrix, all arguments are taken from G. D'Agostini (p.6-8)
883 if (calculateErrors)
884 {
885 printf("Covariance matrix will be calculated... this will take a lot of time (>1 day) ;)\n");
886
887 //Variables and Matrices that will be use along the calculation
888 const Int_t binsV[4] = {NBINSE,NBINSE, NBINSZ, NBINSZ};
889 const Double_t LowEdgeV[4] = {axisLowerLimitE, axisLowerLimitE, axisLowerLimitZ, axisLowerLimitZ};
890 const Double_t UpEdgeV[4] = {axisUpperLimitE, axisUpperLimitE, axisUpperLimitZ, axisUpperLimitZ};
891
892 const Double_t Ntrue = (Double_t)measured->Integral();
893
894 THnSparseF *V = new THnSparseF("V","",4, binsV, LowEdgeV, UpEdgeV);
895 V->Reset();
896 Double_t invCorrContent1, Nt;
897 Double_t invCorrContent2, v11, v12, v2;
898 // calculate V1 and V2
899 for (Int_t idx1=0; idx1<=NFilledBins; idx1++)
900 {
901 printf("Covariance Matrix calculation: iteration idx1=%d of %d\n", idx1, NFilledBins);
902 for (Int_t idx2=0; idx2<=NFilledBins; idx2++)
903 {
904 Int_t bin1[4];
905 Int_t bin2[4];
906 invCorrContent1 = inverseCorrelation->GetBinContent(idx1, bin1);
907 invCorrContent2 = inverseCorrelation->GetBinContent(idx2, bin2);
908 v11=0; v12=0; v2=0;
909 if(bin1[0]>0 && bin1[0]<=NBINSE && bin1[1]>0 && bin1[1]<=NBINSE &&
910 bin1[2]>0 && bin1[2]<=NBINSZ && bin1[3]>0 && bin1[3]<=NBINSZ &&
911 bin2[0]>0 && bin2[0]<=NBINSE && bin2[1]>0 && bin2[1]<=NBINSE &&
912 bin2[2]>0 && bin2[2]<=NBINSZ && bin2[3]>0 && bin2[3]<=NBINSZ)
913 {
914 if (bin1[1]==bin2[1] && bin1[3]==bin2[3])
915 v11 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]
916 *(1. - measuredCopy[bin2[1]][bin2[3]]/Ntrue);
917 else
918 v12 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]*
919 measuredCopy[bin2[1]][bin2[3]]/Ntrue;
920 Nt = (Double_t)prior[bin2[0]][bin2[2]];
921 v2 = measuredCopy[bin1[1]][bin1[3]]*measuredCopy[bin2[1]][bin2[3]]*
922 invCorrContent1*invCorrContent2*
923 BayesUncertaintyTerms(inverseCorrelation, correlation, bin1, bin2, Nt);
924 Int_t binV[4] = {bin1[0],bin2[0],bin1[2],bin2[2]};
925 V->SetBinContent(binV,v11-v12 + v2);
926 }
927 }
928 }
929
930 for(Int_t te = 1; te<=NBINSE; te++)
931 for(Int_t tz = 1; tz<=NBINSZ; tz++)
932 {
933 Int_t binV[4] = {te,te,tz,tz};
934 aResult->SetBinError( te, tz, V->GetBinContent(binV) );
935 }
936
937 TFile* f = new TFile("Covariance_UnfSpectrum.root");
938 f->Open("RECREATE");
939 V->Write();
940 f->Close();
941 }
942
943 return 0;
944
945}
946
947//____________________________________________________________________
948Double_t AliJetSpectrumUnfolding::BayesUncertaintyTerms(THnSparseF *M, THnSparseF *C, Int_t* binTM, Int_t* binTM1, Double_t Nt)
949{
950 //
951 // helper function for the covariance matrix of the bayesian method
952 //
953
954 Double_t result = 0;
955 Float_t term[9];
956 Int_t tmpBin[4], tmpBin1[4];
957 const Int_t nFilledBins = C->GetNbins();
958 if (Nt==0)
959 return 0;
960
961 Float_t CorrContent;
962 Float_t InvCorrContent;
963
964 tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3];
965 tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3];
966 if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
967 {
968 if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
969 term[0] = BayesCov(M, C, tmpBin, tmpBin1)/
970 (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
971 term[2] = term[0]*M->GetBinContent(tmpBin1);
972 }
973 else
974 {
975 term[0] = 0;
976 term[2] = 0;
977 }
978
979 tmpBin[0]=binTM1[0]; tmpBin[1]=binTM[1]; tmpBin[2]=binTM1[2]; tmpBin[3]=binTM[3];
980 tmpBin1[0]=binTM1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM1[2]; tmpBin1[3]=binTM1[3];
981 if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
982 term[6] = BayesCov(M, C, tmpBin, tmpBin1)*
983 M->GetBinContent(tmpBin)/
984 (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
985 else
986 term[6] = 0;
987
988 for(Int_t idx1=0; idx1<=nFilledBins; idx1++)
989 {
990 Int_t bin1[4];
991 CorrContent = C->GetBinContent(idx1, bin1);
992 InvCorrContent = M->GetBinContent(idx1, bin1);
993 if(bin1[0]>0 && bin1[0]<=NBINSE && bin1[1]>0 && bin1[1]<=NBINSE &&
994 bin1[2]>0 && bin1[2]<=NBINSZ && bin1[3]>0 && bin1[3]<=NBINSZ)
995 {
996 tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3];
997 tmpBin1[0]=binTM[0]; tmpBin1[1]=bin1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=bin1[3];
998 if (C->GetBinContent(tmpBin)!=0 &&
999 binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
1000 term[1] = BayesCov(M, C, tmpBin, tmpBin1)/C->GetBinContent(tmpBin);
1001 else
1002 term[1] = 0;
1003
1004 tmpBin[0] =binTM[0]; tmpBin[1] =bin1[1]; tmpBin[2] =binTM[2]; tmpBin[3] =bin1[3];
1005 tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3];
1006 if (C->GetBinContent(tmpBin1)!=0)
1007 {
1008 if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
1009 term[3] = BayesCov(M, C, tmpBin, tmpBin1)/
1010 C->GetBinContent(tmpBin1);
1011 term[5] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin1)/
1012 C->GetBinContent(tmpBin1);
1013 }
1014 else
1015 {
1016 term[3] = 0;
1017 term[5] = 0;
1018 }
1019
1020 tmpBin[0] =binTM1[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM1[2]; tmpBin[3] =binTM[3];
1021 tmpBin1[0]=binTM1[0]; tmpBin1[1]=bin1[1]; tmpBin1[2]=binTM1[2]; tmpBin1[3]=bin1[3];
1022 if (C->GetBinContent(tmpBin)!=0)
1023 term[7] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin)/
1024 C->GetBinContent(tmpBin);
1025 else
1026 term[7] = 0;
1027
1028 tmpBin[0] =bin1[0]; tmpBin[1] =binTM[1]; tmpBin[2] =bin1[2]; tmpBin[3] =binTM[3];
1029 tmpBin1[0]=bin1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=bin1[2]; tmpBin1[3]=binTM1[3];
1030 if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
1031 term[8] = BayesCov(M, C, tmpBin, tmpBin1)*
1032 M->GetBinContent(tmpBin)*M->GetBinContent(tmpBin)/
1033 (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
1034 else
1035 term[8] = 0;
1036
1037 for (Int_t i=0; i<9; i++)
1038 result += term[i]/Nt;
1039 }
1040 }
1041
1042 return result;
1043}
1044
1045//____________________________________________________________________
1046Double_t AliJetSpectrumUnfolding::BayesCov(THnSparseF *M, THnSparseF *correlation, Int_t* binTM, Int_t* bin1)
1047{
1048 Double_t result, result1, result2, result3;
1049
1050 if (binTM[0]==bin1[0] && binTM[2]==bin1[2])
1051 {
1052 if (correlation->GetBinContent(bin1)!=0)
1053 result1 = 1./correlation->GetBinContent(bin1);
1054 else
1055 result1 = 0;
1056 result2 = 1.;
1057 }
1058 else
1059 {
1060 result1 = 0;
1061 result2 = 0;
1062 }
1063
1064 if (binTM[1]==bin1[1] && binTM[3]==bin1[3])
1065 {
1066 Int_t tmpbin[4] = {bin1[0], binTM[1], bin1[2], binTM[3]};
1067 if(correlation->GetBinContent(tmpbin)!=0)
1068 result3 = M->GetBinContent(tmpbin)/correlation->GetBinContent(tmpbin);
1069 else
1070 result3 = 0;
1071 }
1072 else
1073 {
1074 result1 = 0;
1075 result3 = 0;
1076 }
1077
1078 return result = result1 + result2 + result3;
1079}
1080
1081//____________________________________________________________________
1082TH2F* AliJetSpectrumUnfolding::CalculateRecSpectrum(TH2* inputGen)
1083{
1084 // runs the distribution given in inputGen through the correlation histogram identified by
1085 // fCorrelation and produces a reconstructed spectrum
1086
1087 if (!inputGen)
1088 return 0;
1089
1090 // normalize to convert number of events into probability
1091 /*for (Int_t te=1; te<=NBINSE; te++)
1092 for (Int_t tz=1; tz<=NBINSZ; tz++)
1093 {
1094 Int_t bin[4];
1095 Float_t sum=0.;
1096 for (Int_t me = 1; me<=NBINSE; me++)
1097 for (Int_t mz = 1; mz<=NBINSZ; mz++)
1098 {
1099 bin[0] = te; bin[1] = me;
1100 bin[2] = tz; bin[3] = mz;
1101 sum += fCorrelation[correlationMap]->GetBinContent(bin);
1102 }
1103 if (sum > 0.)
1104 for (Int_t me = 1; me<=NBINSE; me++)
1105 for (Int_t mz = 1; mz<=NBINSZ; mz++)
1106 {
1107 bin[0] = te; bin[1] = me;
1108 bin[2] = tz; bin[3] = mz;
1109 fCorrelation[correlationMap]->SetBinContent(bin, fCorrelation[correlationMap]->GetBinContent(bin)/sum);
1110 fCorrelation[correlationMap]->SetBinError(bin, fCorrelation[correlationMap]->GetBinError(bin)/sum);
1111 }
1112 }*/
1113 // normalize to convert number of events into probability (the following loop is much faster)
1114 Float_t sum[NBINSE+2][NBINSZ+2];
1115 memset(sum,0,sizeof(Float_t)*(NBINSE+2)*(NBINSZ+2));
1116
1117 for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++)
1118 {
1119 Int_t bin[4];
1120 Float_t binContent = fCorrelation->GetBinContent(idx, bin);
1121 if (bin[1]>0 && bin[1]<=NBINSE && bin[3]>0 && bin[3]<=NBINSZ){
1122 sum[bin[0]][bin[2]] += binContent;
1123 }
1124 }
1125
1126 for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++)
1127 {
1128 Int_t bin[4];
1129 Float_t binContent = fCorrelation->GetBinContent(idx, bin);
1130 Float_t binError = fCorrelation->GetBinError(bin);
1131 if (sum[bin[0]][bin[2]]>0 && bin[1]>0 && bin[1]<=NBINSE &&
1132 bin[3]>0 && bin[3]<=NBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=NBINSE && bin[2]<=NBINSZ)
1133 {
1134 fCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]);
1135 fCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]);
1136 }
1137 }
1138
1139 TH2F* target = dynamic_cast<TH2F*> (fRecSpectrum->Clone(Form("reconstructed_%s", inputGen->GetName())));
1140 target->Reset();
1141
1142 for (Int_t me=1; me<=NBINSE; ++me)
1143 for (Int_t mz=1; mz<=NBINSZ; ++mz)
1144 {
1145 Float_t measured = 0;
1146 Float_t error = 0;
1147
1148 for (Int_t te=1; te<=NBINSE; ++te)
1149 for (Int_t tz=1; tz<=NBINSZ; ++tz)
1150 {
1151 Int_t bin[4] = {te, me, tz, mz};
1152 measured += inputGen->GetBinContent(te,tz) * fCorrelation->GetBinContent(bin);
1153 error += inputGen->GetBinError(te,tz) * fCorrelation->GetBinContent(bin);
1154 }
1155 target->SetBinContent(me, mz, measured);
1156 target->SetBinError(me, mz, error);
1157 }
1158
1159 return target;
1160}
1161
1162//__________________________________________________________________________________________________
1163void AliJetSpectrumUnfolding::SetGenRecFromFunc(TF2* inputGen)
1164{
1165 // uses the given function to fill the input Generated histogram and generates from that
1166 // the reconstructed histogram by applying the response histogram
1167 // this can be used to evaluate if the methods work indepedently of the input
1168 // distribution
1169
1170 if (!inputGen)
1171 return;
1172
1173 TH2F* histtmp = new TH2F("histtmp", "tmp", EBINNING, ZBINNING);
1174 TH2F* gen = fGenSpectrum;
1175
1176 histtmp->Reset();
1177 gen->Reset();
1178
1179 histtmp->FillRandom(inputGen->GetName(), NEVENTS);
1180
1181 for (Int_t i=1; i<=gen->GetNbinsX(); ++i)
1182 for (Int_t j=1; j<=gen->GetNbinsY(); ++j)
1183 {
1184 gen->SetBinContent(i, j, histtmp->GetBinContent(i,j));
1185 gen->SetBinError(i, j, histtmp->GetBinError(i,j));
1186 }
1187
1188 delete histtmp;
1189
1190 //new TCanvas;
1191 //gStyle->SetPalette(1);
1192 //gPad->SetLogz();
1193 //gen->Draw("COLZ");
1194
1195
1196 TH2 *recsave = fRecSpectrum;
1197
1198 fRecSpectrum = CalculateRecSpectrum(gen);
1199 fRecSpectrum->SetName(recsave->GetName());
1200 delete recsave;
1201
1202 return;
1203}
1204//________________________________________________________________________________________