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