New book-keeping class
[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 #include <TVectorD.h>
42
43 #include <THnSparse.h>
44
45 ClassImp(AliJetSpectrumUnfolding)
46
47 const Int_t AliJetSpectrumUnfolding::fgkNBINSE  = 50;
48 const Int_t AliJetSpectrumUnfolding::fgkNBINSZ  = 50;
49 const Int_t AliJetSpectrumUnfolding::fgkNEVENTS = 500000;
50 const Double_t AliJetSpectrumUnfolding::fgkaxisLowerLimitE = 0.;
51 const Double_t AliJetSpectrumUnfolding::fgkaxisLowerLimitZ = 0.;
52 const Double_t AliJetSpectrumUnfolding::fgkaxisUpperLimitE = 250.;
53 const Double_t AliJetSpectrumUnfolding::fgkaxisUpperLimitZ = 1.;
54
55 Float_t AliJetSpectrumUnfolding::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
56 Int_t   AliJetSpectrumUnfolding::fgBayesianIterations = 100;         // number of iterations in Bayesian method
57
58 //____________________________________________________________________
59
60 AliJetSpectrumUnfolding::AliJetSpectrumUnfolding() :
61   TNamed(), fCurrentRec(0), fCurrentCorrelation(0), fRecSpectrum(0),  fGenSpectrum(0),
62   fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
63 {
64   //
65   // default constructor
66   //
67
68   fGenSpectrum = 0;
69   fRecSpectrum = 0;
70   fUnfSpectrum = 0;
71   fCorrelation = 0;
72 }
73
74 //____________________________________________________________________
75 AliJetSpectrumUnfolding::AliJetSpectrumUnfolding(const Char_t* name, const Char_t* title) :
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)
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);
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};
97   //arrays for bin limits
98   Double_t lowEdge[4] = {fgkaxisLowerLimitE, fgkaxisLowerLimitE, fgkaxisLowerLimitZ, fgkaxisLowerLimitZ};
99   Double_t upEdge[4]  = {fgkaxisUpperLimitE, fgkaxisUpperLimitE, fgkaxisUpperLimitZ, fgkaxisUpperLimitZ};
100
101   fCorrelation = new THnSparseF("fCorrelation", "Correlation Function", 4, nbin, lowEdge, upEdge);
102
103   TH1::AddDirectory(oldStatus);
104 }
105
106 //____________________________________________________________________
107 AliJetSpectrumUnfolding::~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 //____________________________________________________________________
132 Long64_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 //____________________________________________________________________
176 Bool_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 //____________________________________________________________________
226 void 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 //____________________________________________________________________
251 void 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
266   Printf("Correlation Matrix has %.0E filled bins", fCurrentCorrelation->GetNbins());
267
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       {
275         if (fCurrentRec->GetBinContent(me,mz) <= 5 && me>fgkNBINSE/2 && mz>fgkNBINSZ/2)
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};
337         if ( (bin[1]>maxBinE && bin[1]<=fgkNBINSE) && (bin[3]>maxBinZ && bin[3]<=fgkNBINSZ) )
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         } 
344         printf("create big bin1, nbins = %d, te  = %d, tz = %d\n", fgkNBINSE, bin[0], bin[1]);
345       }
346
347       printf("Adjusted correlation matrix!\n");
348     }
349   } // end Create Big Bin
350
351 }
352
353 //____________________________________________________________________
354 void 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 //____________________________________________________________________
367 void AliJetSpectrumUnfolding::NormalizeToBinWidth(TH2* const hist)
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 //____________________________________________________________________
383 void 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);
396   canvas2->cd();
397   gPad->SetLogz();    
398   fGenSpectrum->DrawCopy("COLZ");
399
400   TCanvas* canvas3 = new TCanvas("fUnfSpectrum", "fUnfSpectrum", 900, 600);
401   canvas3->cd();
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 //____________________________________________________________________
427 void AliJetSpectrumUnfolding::DrawComparison(const char* name, TH2* const genHist)
428 {
429   // 
430   // Draws the copmparison plot (gen,rec and unfolded distributions
431   //
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;
488   const Int_t nRGBs = 5;
489   const Int_t nCont = 500;
490
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);
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 //____________________________________________________________________
576 void AliJetSpectrumUnfolding::GetComparisonResults(Float_t* const gen, Int_t*  const genLimit, Float_t* const residuals, Float_t* const ratioAverage) const
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 //____________________________________________________________________
593 void AliJetSpectrumUnfolding::ApplyBayesianMethod(Float_t regPar, Int_t nIterations, TH2* const initialConditions, Bool_t determineError)
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     }*/
628   Float_t sum[fgkNBINSE+2][fgkNBINSZ+2];
629   memset(sum,0,sizeof(Float_t)*(fgkNBINSE+2)*(fgkNBINSZ+2));
630
631   for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++)
632   {
633     Int_t bin[4];
634     Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin);
635     if ( (bin[1]>0 && bin[1]<=fgkNBINSE) && (bin[3]>0 && bin[3]<=fgkNBINSZ) )
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);
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) )
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 //____________________________________________________________________
724 Int_t AliJetSpectrumUnfolding::UnfoldWithBayesian(THnSparseF* const correlation, TH2* const measured, TH2* const initialConditions, TH2* const aResult, Float_t regPar, Int_t nIterations, Bool_t calculateErrors)
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   }
735   const Int_t nFillesBins = correlation->GetNbins();  
736   const Int_t kStartBin = 1;
737
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
742   
743   printf("NbinsE=%d - NbinsZ=%d\n", fgkNBINSE, fgkNBINSZ);
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);
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)
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);
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)
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;
879         //        Float_t binError = 0;
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    
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};
901     
902     const Double_t nTrue = (Double_t)measured->Integral();
903     
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;        
908     // calculate V1 and V2
909     for (Int_t idx1=0; idx1<=nFillesBins; idx1++)
910     {
911       printf("Covariance Matrix calculation: iteration idx1=%d of %d\n", idx1, nFillesBins);
912       for (Int_t idx2=0; idx2<=nFillesBins; idx2++)
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;
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)
923         {   
924           if (bin1[1]==bin2[1] && bin1[3]==bin2[3])    
925             v11 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]
926                   *(1. - measuredCopy[bin2[1]][bin2[3]]/nTrue);                       
927           else
928             v12 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]*
929                   measuredCopy[bin2[1]][bin2[3]]/nTrue;
930           nt = (Double_t)prior[bin2[0]][bin2[2]];
931           v2 = measuredCopy[bin1[1]][bin1[3]]*measuredCopy[bin2[1]][bin2[3]]*
932                invCorrContent1*invCorrContent2*
933                BayesUncertaintyTerms(inverseCorrelation, correlation, bin1, bin2, nt);   
934           Int_t binV[4] = {bin1[0],bin2[0],bin1[2],bin2[2]};         
935           v->SetBinContent(binV,v11-v12 + v2);
936         }  
937       }
938     }    
939
940     for(Int_t te = 1; te<=fgkNBINSE; te++)
941       for(Int_t tz = 1; tz<=fgkNBINSZ; tz++)
942       {
943         Int_t binV[4] = {te,te,tz,tz}; 
944         aResult->SetBinError( te, tz, v->GetBinContent(binV) );
945       }            
946       
947     TFile* f = new TFile("Covariance_UnfSpectrum.root");
948     f->Open("RECREATE");
949     v->Write();
950     f->Close();    
951   }  
952   
953   return 0;
954
955 }
956
957 //____________________________________________________________________
958 Double_t AliJetSpectrumUnfolding::BayesUncertaintyTerms(THnSparseF* const M, THnSparseF* const C, Int_t* const binTM, Int_t* const binTM1, Double_t nt)
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();
968   if (nt==0)
969     return 0;
970     
971   Float_t corrContent;
972   Float_t invCorrContent;
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];
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)
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++)
1048         result += term[i]/nt;                    
1049     }          
1050   }
1051    
1052   return result;
1053 }
1054
1055 //____________________________________________________________________
1056 Double_t AliJetSpectrumUnfolding::BayesCov(THnSparseF* const M, THnSparseF* const correlation, Int_t* const binTM, Int_t* const bin1)
1057 {
1058
1059   //
1060   // get the covariance matrix  
1061   //
1062
1063
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 //____________________________________________________________________
1098 TH2F* AliJetSpectrumUnfolding::CalculateRecSpectrum(TH2* const inputGen)
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)
1130   Float_t sum[fgkNBINSE+2][fgkNBINSZ+2];
1131   memset(sum,0,sizeof(Float_t)*(fgkNBINSE+2)*(fgkNBINSZ+2));
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);
1137     if (bin[1]>0 && bin[1]<=fgkNBINSE && bin[3]>0 && bin[3]<=fgkNBINSZ){
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);
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) 
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
1158   for (Int_t me=1; me<=fgkNBINSE; ++me)
1159     for (Int_t mz=1; mz<=fgkNBINSZ; ++mz)
1160     {
1161       Float_t measured = 0;
1162       Float_t error = 0;
1163
1164       for (Int_t te=1; te<=fgkNBINSE; ++te)
1165         for (Int_t tz=1; tz<=fgkNBINSZ; ++tz)
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 //__________________________________________________________________________________________________
1179 void AliJetSpectrumUnfolding::SetGenRecFromFunc(TF2* const inputGen)
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
1189   TH2F* histtmp = new TH2F("histtmp", "tmp", 
1190                           fgkNBINSE, fgkaxisLowerLimitE, fgkaxisUpperLimitE,
1191                           fgkNBINSZ, fgkaxisLowerLimitZ, fgkaxisUpperLimitZ);
1192
1193   TH2F* gen  = fGenSpectrum;
1194
1195   histtmp->Reset();
1196   gen->Reset();
1197
1198   histtmp->FillRandom(inputGen->GetName(), fgkNEVENTS);
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 //________________________________________________________________________________________