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