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