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