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