Adding the macros and the class for 2-dim unfolding (Ydalia)
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 7 Jan 2009 12:56:48 +0000 (12:56 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 7 Jan 2009 12:56:48 +0000 (12:56 +0000)
PWG4/CMake_libPWG4JetTasks.txt
PWG4/JetTasks/AliJetSpectrumUnfolding.cxx [new file with mode: 0644]
PWG4/JetTasks/AliJetSpectrumUnfolding.h [new file with mode: 0644]
PWG4/PWG4JetTasksLinkDef.h
PWG4/libPWG4JetTasks.pkg
PWG4/macros/AnalysisTrainCAF.C
PWG4/macros/runJetSpectrumUnfolding.C [new file with mode: 0644]

index b76e6cc..2ad98b5 100644 (file)
@@ -1,7 +1,11 @@
 # -*- mode: cmake -*-
 
 set(SRCS
-       JetTasks/AliAnalysisTaskUE.cxx JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnaESDSpectraQA.cxx
+       JetTasks/AliAnalysisTaskUE.cxx
+       JetTasks/AliAnalysisTaskJetSpectrum.cxx
+       JetTasks/AliAnalysisHelperJetTasks.cxx
+       JetTasks/AliAnaESDSpectraQA.cxx
+       JetTasks/AliJetSpectrumUnfolding.cxx
 )
 
 # fill list of header files from list of source files
diff --git a/PWG4/JetTasks/AliJetSpectrumUnfolding.cxx b/PWG4/JetTasks/AliJetSpectrumUnfolding.cxx
new file mode 100644 (file)
index 0000000..b1b7fa6
--- /dev/null
@@ -0,0 +1,1204 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+//
+// This class is used to store correlation maps, generated and reconstructed data of the jet spectrum
+// It also contains functions to correct the spectrum using the bayesian unfolding
+//
+
+#include "AliJetSpectrumUnfolding.h"
+
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TDirectory.h>
+#include <TVirtualFitter.h>
+#include <TCanvas.h>
+#include <TString.h>
+#include <TF1.h>
+#include <TF2.h>
+#include <TMath.h>
+#include <TCollection.h>
+#include <TLegend.h>
+#include <TLine.h>
+#include <TRandom.h>
+#include <TProfile.h>
+#include <TProfile2D.h>
+#include <TStyle.h>
+#include <TColor.h>
+
+#include <THnSparse.h>
+
+ClassImp(AliJetSpectrumUnfolding)
+
+const Int_t NBINSE  = 50;
+const Int_t NBINSZ  = 50;
+const Int_t NEVENTS = 500000;
+const Double_t axisLowerLimitE = 0.;
+const Double_t axisLowerLimitZ = 0.;
+const Double_t axisUpperLimitE = 250.;
+const Double_t axisUpperLimitZ = 1.;
+
+Float_t AliJetSpectrumUnfolding::fgBayesianSmoothing  = 1;           // smoothing parameter (0 = no smoothing)
+Int_t   AliJetSpectrumUnfolding::fgBayesianIterations = 100;         // number of iterations in Bayesian method
+
+//____________________________________________________________________
+
+AliJetSpectrumUnfolding::AliJetSpectrumUnfolding() :
+  TNamed(), fCurrentRec(0), fCurrentCorrelation(0), fGenSpectrum(0),
+  fRecSpectrum(0), fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
+{
+  //
+  // default constructor
+  //
+
+  fGenSpectrum = 0;
+  fRecSpectrum = 0;
+  fUnfSpectrum = 0;
+  fCorrelation = 0;
+}
+
+//____________________________________________________________________
+AliJetSpectrumUnfolding::AliJetSpectrumUnfolding(const Char_t* name, const Char_t* title) :
+  TNamed(name, title), fCurrentRec(0), fCurrentCorrelation(0), fGenSpectrum(0),
+  fRecSpectrum(0), fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
+{
+  //
+  // named constructor
+  //
+
+  // do not add this hists to the directory
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  #define ZBINNING NBINSZ, axisLowerLimitZ, axisUpperLimitZ       // fragmentation of leading particle
+  #define EBINNING NBINSE, axisLowerLimitE, axisUpperLimitE       // energy of the jet
+
+  fRecSpectrum = new TH2F("fRecSpectrum", "Reconstructed Spectrum;E^{jet}_{rec} [GeV];z^{lp}_{rec}", EBINNING, ZBINNING);
+  fGenSpectrum = new TH2F("fGenSpectrum", "Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}_{gen}", EBINNING, ZBINNING);
+  fUnfSpectrum = new TH2F("fUnfSpectrum", "Unfolded Spectrum;E^{jet} [GeV];z^{lp}", EBINNING, ZBINNING);
+
+  const Int_t nbin[4]={NBINSE, NBINSE, NBINSZ, NBINSZ};
+  //arrays for bin limits
+  Double_t LowEdge[4] = {axisLowerLimitE, axisLowerLimitE, axisLowerLimitZ, axisLowerLimitZ};
+  Double_t UpEdge[4]  = {axisUpperLimitE, axisUpperLimitE, axisUpperLimitZ, axisUpperLimitZ};
+
+  fCorrelation = new THnSparseF("fCorrelation", "Correlation Function", 4, nbin, LowEdge, UpEdge);
+
+  TH1::AddDirectory(oldStatus);
+}
+
+//____________________________________________________________________
+AliJetSpectrumUnfolding::~AliJetSpectrumUnfolding()
+{
+  //
+  // Destructor
+  //
+
+  if (fGenSpectrum)
+    delete fGenSpectrum;
+  fGenSpectrum = 0;
+
+  if (fRecSpectrum)
+    delete fRecSpectrum;
+  fRecSpectrum = 0;
+
+  if (fUnfSpectrum)
+    delete fUnfSpectrum;
+  fUnfSpectrum = 0;
+  if (fCorrelation)
+    delete fCorrelation;
+  fCorrelation = 0;
+
+}
+
+//____________________________________________________________________
+Long64_t AliJetSpectrumUnfolding::Merge(TCollection* list)
+{
+  // Merge a list of AliJetSpectrumUnfolding objects with this (needed for
+  // PROOF).
+  // Returns the number of merged objects (including this).
+
+  if (!list)
+    return 0;
+
+  if (list->IsEmpty())
+    return 1;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj;
+
+  // collections of all histograms
+  TList collections[4];
+
+  Int_t count = 0;
+  while ((obj = iter->Next())) {
+
+    AliJetSpectrumUnfolding* entry = dynamic_cast<AliJetSpectrumUnfolding*> (obj);
+    if (entry == 0)
+      continue;
+
+    collections[0].Add(entry->fGenSpectrum);
+    collections[1].Add(entry->fRecSpectrum);
+    collections[2].Add(entry->fUnfSpectrum);
+    collections[3].Add(entry->fCorrelation);
+
+    count++;
+  }
+
+  fGenSpectrum->Merge(&collections[0]);
+  fRecSpectrum->Merge(&collections[1]);
+  fUnfSpectrum->Merge(&collections[2]);
+  fCorrelation->Merge(&collections[3]);
+
+  delete iter;
+
+  return count+1;
+}
+
+//____________________________________________________________________
+Bool_t AliJetSpectrumUnfolding::LoadHistograms(const Char_t* dir)
+{
+  //
+  // loads the histograms from a file
+  // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
+  //
+
+  if (!dir)
+    dir = GetName();
+
+  if (!gDirectory->cd(dir))
+    return kFALSE;
+
+  Bool_t success = kTRUE;
+
+  // store old histograms to delete them later
+  TList oldHistograms;
+  oldHistograms.SetOwner(1);
+
+  if (fGenSpectrum)  oldHistograms.Add(fGenSpectrum);
+  if (fRecSpectrum)  oldHistograms.Add(fRecSpectrum);
+  if (fUnfSpectrum)  oldHistograms.Add(fUnfSpectrum);
+  if (fCorrelation)  oldHistograms.Add(fCorrelation);
+
+  // load new histograms
+  fGenSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fGenSpectrum->GetName()));
+  if (!fGenSpectrum)
+    success = kFALSE;
+
+  fRecSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fRecSpectrum->GetName()));
+  if (!fRecSpectrum)
+    success = kFALSE;
+
+  fUnfSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fUnfSpectrum->GetName()));
+  if (!fUnfSpectrum)
+    success = kFALSE;
+
+  fCorrelation = dynamic_cast<THnSparseF*> (gDirectory->Get(fCorrelation->GetName()));
+  if (!fCorrelation)
+    success = kFALSE;
+
+  gDirectory->cd("..");
+
+  // delete old histograms
+  oldHistograms.Delete();
+
+  return success;
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::SaveHistograms()
+{
+  //
+  // saves the histograms
+  //
+
+  gDirectory->mkdir(GetName());
+  gDirectory->cd(GetName());
+
+  if (fGenSpectrum)
+    fGenSpectrum->Write();
+
+  if (fRecSpectrum)
+    fRecSpectrum->Write();
+
+  if (fUnfSpectrum)
+    fUnfSpectrum->Write();
+
+  if (fCorrelation)
+    fCorrelation->Write();
+
+  gDirectory->cd("..");
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::SetupCurrentHists(Bool_t createBigBin)
+{
+  //
+  // resets fUnfSpectrum
+  //
+
+  fUnfSpectrum->Reset();
+  fUnfSpectrum->Sumw2();
+
+  fCurrentRec = (TH2F*)fRecSpectrum->Clone("fCurrentRec");
+  fCurrentRec->Sumw2();
+
+  fCurrentCorrelation = (THnSparseF*)fCorrelation->Clone("fCurrentCorrelation");  
+  fCurrentCorrelation->Sumw2();
+
+  if (createBigBin)
+  {
+    Int_t maxBinE = 0, maxBinZ = 0;
+    Float_t maxE = 0, maxZ = 0;
+    for (Int_t me=1; me<=fCurrentRec->GetNbinsX(); me++)
+      for (Int_t mz=1; mz<=fCurrentRec->GetNbinsY(); mz++)
+      {
+        if (fCurrentRec->GetBinContent(me,mz) <= 5 && me>NBINSE/2 && mz>NBINSZ/2)
+        {
+          maxBinE = me;
+         maxBinZ = mz;
+         maxE = fCurrentRec->GetXaxis()->GetBinCenter(me);
+         maxZ = fCurrentRec->GetYaxis()->GetBinCenter(mz);
+          break;
+        }
+      }
+
+    if (maxBinE > 0 || maxBinZ > 0)
+    {
+      printf("Bin limit in measured spectrum is e = %d and z = %d.\n", maxBinE, maxBinZ);
+      fCurrentRec->SetBinContent(maxBinE, maxBinZ, fCurrentRec->Integral(maxBinE, fCurrentRec->GetNbinsX(), maxBinZ, fCurrentRec->GetNbinsY()));
+      for (Int_t me=maxBinE+1; me<=fCurrentRec->GetNbinsX(); me++)
+        for (Int_t mz=maxBinZ+1; mz<=fCurrentRec->GetNbinsY(); mz++)
+        {
+          fCurrentRec->SetBinContent(me, mz, 0);
+          fCurrentRec->SetBinError(me, mz, 0);
+        }
+      // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
+      fCurrentRec->SetBinError(maxBinE, maxBinZ, TMath::Sqrt(fCurrentRec->GetBinContent(maxBinE, maxBinZ)));
+
+      printf("This bin has now %f +- %f entries\n", fCurrentRec->GetBinContent(maxBinE, maxBinZ), fCurrentRec->GetBinError(maxBinE, maxBinZ));
+
+    /*  for (Int_t te=1; te<=NBINSE; te++)
+      {
+        for (Int_t tz=1; tz<=NBINSZ; tz++)
+       {
+          Int_t binMin[4] = {te, maxBinE, tz, maxBinZ};
+         Int_t binMax[4] = {NBINSE, NBINSE, NBINSZ, NBINSZ};
+         Float_t sum=0;
+         for (Int_t ite=te; ite<=NBINSE; ite++)
+           for (Int_t itz=tz; itz<=NBINSZ; itz++)
+             for (Int_t ime=maxBinE; ime<=NBINSE; ime++)
+               for (Int_t imz=maxBinZ; imz<=NBINSZ; imz++)
+               {
+                 Int_t bin[4] = {ite, ime, itz, imz};
+                 sum += fCurrentCorrelation->GetBinContent(bin);
+               }
+         fCurrentCorrelation->SetBinContent(binMin, sum);
+          fCurrentCorrelation->SetBinError(binMin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin)));
+          printf("create big bin1, nbins = %d, te  = %d, tz = %d\n", NBINSE, te, tz);
+          for (Int_t me=maxBinE; me<=NBINSE; me++)
+          {
+            for (Int_t mz=maxBinZ; mz<=NBINSZ; mz++)
+           {
+             Int_t bin[4] = {te, me, tz, mz};
+              fCurrentCorrelation->SetBinContent(bin, 0.);
+              fCurrentCorrelation->SetBinError(bin, 0.);
+              printf("create big bin2\n");
+           }
+         }
+        }
+      }*/
+      
+      for(Int_t idx = 0; idx<=fCurrentCorrelation->GetNbins(); idx++)
+      {
+        Int_t bin[4];
+        Float_t binContent = fCurrentCorrelation->GetBinContent(idx,bin);
+        Float_t binError   = fCurrentCorrelation->GetBinError(idx);
+        Int_t binMin[4] = {bin[0], maxBinE, bin[2], maxBinZ};
+        if ( (bin[1]>maxBinE && bin[1]<=NBINSE) && (bin[3]>maxBinZ && bin[3]<=NBINSZ) )
+        {
+         fCurrentCorrelation->SetBinContent(binMin, binContent + fCurrentCorrelation->GetBinContent(binMin));
+          fCurrentCorrelation->SetBinError(binMin, binError + TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin)));
+          fCurrentCorrelation->SetBinContent(bin, 0.);
+          fCurrentCorrelation->SetBinError(bin, 0.);         
+        } 
+        printf("create big bin1, nbins = %d, te  = %d, tz = %d\n", NBINSE, bin[0], bin[1]);
+      }
+
+      printf("Adjusted correlation matrix!\n");
+    }
+  } // end Create Big Bin
+
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
+{
+  //
+  // sets the parameters for Bayesian unfolding
+  //
+
+  fgBayesianSmoothing = smoothing;
+  fgBayesianIterations = nIterations;
+
+  printf("AliJetSpectrumUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::NormalizeToBinWidth(TH2* hist)
+{
+  //
+  // normalizes a 2-d histogram to its bin width (x width * y width)
+  //
+
+  for (Int_t i=1; i<=hist->GetNbinsX(); i++)
+    for (Int_t j=1; j<=hist->GetNbinsY(); j++)
+    {
+      Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
+      hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
+      hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
+    }
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::DrawHistograms()
+{
+  //
+  // draws the histograms of this class
+  //
+
+  gStyle->SetPalette(1);
+
+  TCanvas* canvas1 = new TCanvas("fRecSpectrum", "fRecSpectrum", 900, 600);
+  gPad->SetLogz();    
+  fRecSpectrum->DrawCopy("COLZ");
+
+  TCanvas* canvas2 = new TCanvas("fGenSpectrum", "fGenSpectrum", 900, 600);
+  gPad->SetLogz();    
+  fGenSpectrum->DrawCopy("COLZ");
+
+  TCanvas* canvas3 = new TCanvas("fUnfSpectrum", "fUnfSpectrum", 900, 600);
+  gPad->SetLogz();    
+  fUnfSpectrum->DrawCopy("COLZ");
+
+  TCanvas* canvas4 = new TCanvas("fCorrelation", "fCorrelation", 500, 500);
+  canvas1->Divide(2);  
+
+  canvas4->cd(1);
+  gPad->SetLogz();
+  TH2D* h0 = fCorrelation->Projection(1,0);
+  h0->SetXTitle("E^{jet}_{gen} [GeV]");
+  h0->SetYTitle("E^{jet}_{rec} [GeV]");
+  h0->SetTitle("Projection: Jet Energy");    
+  h0->DrawCopy("colz");
+
+  canvas1->cd(2);
+  gPad->SetLogz();  
+  TH2D* h1 = fCorrelation->Projection(3,2);
+  h1->SetXTitle("z^{lp}_{gen}");
+  h1->SetYTitle("z^{lp}_{rec}");
+  h1->SetTitle("Projection: Leading Particle Fragmentation");        
+  h1->DrawCopy("colz");
+
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::DrawComparison(const char* name, TH2* genHist)
+{
+
+  if (fUnfSpectrum->Integral() == 0)
+  {
+    printf("ERROR. Unfolded histogram is empty\n");
+    return;
+  }
+
+  //regain measured distribution used for unfolding, because the bins were modified in SetupCurrentHists
+  //in create big bin
+  fCurrentRec = (TH2F*)fRecSpectrum->Clone();
+  fCurrentRec->Sumw2();
+  fCurrentRec->Scale(1.0 / fCurrentRec->Integral());
+
+  // normalize unfolded result to 1
+  fUnfSpectrum->Scale(1.0 / fUnfSpectrum->Integral());
+
+  // find bin with <= 100 entries. this is used as limit for the chi2 calculation
+  Int_t mcBinLimitE = 0, mcBinLimitZ = 0;
+  for (Int_t i=0; i<genHist->GetNbinsX(); ++i)
+    for (Int_t j=0; j<genHist->GetNbinsY(); ++j)
+    {
+      if (genHist->GetBinContent(i,j) > 100)
+      {
+        mcBinLimitE = i;
+       mcBinLimitZ = j;
+      }
+      else
+        break;
+    }
+  Printf("AliJetSpectrumUnfolding::DrawComparison: Gen bin limit is (x,y) = (%d,%d)", mcBinLimitE,mcBinLimitZ);
+
+  // scale to 1 true spectrum
+  genHist->Sumw2();
+  genHist->Scale(1.0 / genHist->Integral());
+
+  // calculate residual
+  // for that we convolute the response matrix with the gathered result
+  TH2* tmpRecRecorrected = (TH2*) fUnfSpectrum->Clone("tmpRecRecorrected");
+  TH2* convoluted = CalculateRecSpectrum(tmpRecRecorrected);
+  if (convoluted->Integral() > 0)
+    convoluted->Scale(1.0 / convoluted->Integral());
+  else
+    printf("ERROR: convoluted is empty. Something went wrong calculating the convoluted histogram.\n");
+
+  TH2* residual = (TH2*) convoluted->Clone("residual");
+  residual->SetTitle("(R#otimesUnfolded - Reconstructed)/Reconstructed;E^{jet} [GeV]; z^{lp}");
+
+  fCurrentRec->Scale(1./fCurrentRec->Integral());
+  residual->Add(fCurrentRec, -1);
+  //residual->Divide(residual, fCurrentRec, 1, 1, "B");
+
+  // draw canvas
+  TCanvas* canvas1 = new TCanvas(name, name, 1000, 1000);
+  canvas1->Divide(2, 3);
+  
+  Int_t style = 1;
+  const Int_t NRGBs = 5;
+  const Int_t NCont = 500;
+
+  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+  Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
+  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
+  Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
+  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
+  gStyle->SetNumberContours(NCont);
+
+  canvas1->cd(1);
+  gStyle->SetPalette(style);
+  gPad->SetLogz();
+  genHist->SetTitle("Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}");
+  genHist->SetStats(0);
+  genHist->DrawCopy("colz");
+
+  canvas1->cd(2);
+  gStyle->SetPalette(style);
+  gPad->SetLogz();
+  fUnfSpectrum->SetStats(0);
+  fUnfSpectrum->DrawCopy("colz");
+
+  canvas1->cd(3);
+  gStyle->SetPalette(style);
+  gPad->SetLogz();
+  fCurrentRec->SetTitle(fRecSpectrum->GetTitle());
+  fCurrentRec->SetStats(0);
+  fCurrentRec->DrawCopy("colz");
+
+  canvas1->cd(4);
+  gStyle->SetPalette(style);
+  gPad->SetLogy();  
+  TH1D* projGenX = genHist->ProjectionX();
+  projGenX->SetTitle("Projection: Jet Energy; E^{jet} [GeV]");
+  TH1D* projUnfX = fUnfSpectrum->ProjectionX();
+  TH1D* projRecX = fCurrentRec->ProjectionX();  
+  projGenX->SetStats(0);
+  projRecX->SetStats(0);
+  projUnfX->SetStats(0);  
+  projGenX->SetLineColor(8);
+  projRecX->SetLineColor(2);  
+  projGenX->DrawCopy();
+  projUnfX->DrawCopy("same");
+  projRecX->DrawCopy("same");  
+
+  TLegend* legend = new TLegend(0.6, 0.85, 0.98, 0.98);
+  legend->AddEntry(projGenX, "Generated Spectrum");
+  legend->AddEntry(projUnfX, "Unfolded Spectrum");
+  legend->AddEntry(projRecX, "Reconstructed Spectrum");  
+  //legend->SetFillColor(0);
+  legend->Draw("same");
+
+  canvas1->cd(5);
+  gPad->SetLogy();
+  gStyle->SetPalette(style);
+  TH1D* projGenY = genHist->ProjectionY();
+  projGenY->SetTitle("Projection: Leading Particle Fragmentation; z^{lp}");
+  TH1D* projUnfY = fUnfSpectrum->ProjectionY();
+  TH1D* projRecY = fCurrentRec->ProjectionY();  
+  projGenY->SetStats(0);
+  projRecY->SetStats(0);
+  projUnfY->SetStats(0);  
+  projGenY->SetLineColor(8);
+  projRecY->SetLineColor(2);  
+  projGenY->DrawCopy();
+  projUnfY->DrawCopy("same");
+  projRecY->DrawCopy("same");  
+
+  TLegend* legend1 = new TLegend(0.6, 0.85, 0.98, 0.98);
+  legend1->AddEntry(projGenY, "Generated Spectrum");
+  legend1->AddEntry(projUnfY, "Unfolded Spectrum");
+  legend1->AddEntry(projRecY, "Recontructed Spectrum");  
+  //legend1->SetFillColor(0);
+  legend1->Draw("same");
+  
+  // Draw residuals
+  canvas1->cd(6);
+  gStyle->SetPalette(style);
+  gPad->SetLogz();
+  residual->SetStats(0);
+  residual->DrawCopy("colz");
+  
+  canvas1->SaveAs(Form("%s.png", canvas1->GetName()));
+}
+
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::GetComparisonResults(Float_t* gen, Int_t* genLimit, Float_t* residuals, Float_t* ratioAverage) const
+{
+  // Returns the chi2 between the Generated and the unfolded Reconstructed spectrum as well as between the Reconstructed and the folded unfolded
+  // These values are computed during DrawComparison, thus this function picks up the
+  // last calculation
+
+  if (gen)
+    *gen = fLastChi2MC;
+  if (genLimit)
+    *genLimit = fLastChi2MCLimit;
+  if (residuals)
+    *residuals = fLastChi2Residuals;
+  if (ratioAverage)
+    *ratioAverage = fRatioAverage;
+}
+
+//____________________________________________________________________
+void AliJetSpectrumUnfolding::ApplyBayesianMethod(Float_t regPar, Int_t nIterations, TH2* initialConditions, Bool_t determineError)
+{
+  //
+  // correct spectrum using bayesian unfolding
+  //
+
+  // initialize seed with current time 
+  gRandom->SetSeed(0);
+
+  printf("seting up current arrays and histograms...\n");
+  SetupCurrentHists(kFALSE); // kFALSE to not create big bin
+
+  // normalize Correlation Map to convert number of events into probabilities
+  /*for (Int_t te=1; te<=NBINSE; te++)
+    for (Int_t tz=1; tz<=NBINSZ; tz++)
+    {
+       Int_t bin[4];
+       Float_t sum=0.;
+       for (Int_t me = 1; me<=NBINSE; me++)
+         for (Int_t mz = 1; mz<=NBINSZ; mz++)
+         {
+           bin[0] = te; bin[1] = me; 
+           bin[2] = tz; bin[3] = mz;           
+           sum += fCurrentCorrelation->GetBinContent(bin);
+         }
+       if (sum > 0.)
+         for (Int_t me = 1; me<=NBINSE; me++)
+           for (Int_t mz = 1; mz<=NBINSZ; mz++)
+           {
+             bin[0] = te; bin[1] = me; 
+             bin[2] = tz; bin[3] = mz;           
+             fCurrentCorrelation->SetBinContent(bin, fCurrentCorrelation->GetBinContent(bin)/sum);
+            fCurrentCorrelation->SetBinError(bin, fCurrentCorrelation->GetBinError(bin)/sum);
+          }
+    }*/
+  Float_t sum[NBINSE+2][NBINSZ+2];
+  memset(sum,0,sizeof(Float_t)*(NBINSE+2)*(NBINSZ+2));
+
+  for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++)
+  {
+    Int_t bin[4];
+    Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin);
+    if ( (bin[1]>0 && bin[1]<=NBINSE) && (bin[3]>0 && bin[3]<=NBINSZ) )
+      sum[bin[0]][bin[2]] += binContent; 
+  }
+  
+  for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++)
+  {
+    Int_t bin[4];
+    Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin);
+    Float_t binError   = fCurrentCorrelation->GetBinError(bin);
+    if (sum[bin[0]][bin[2]]>0 && (bin[1]>0 && bin[1]<=NBINSE) &&
+        (bin[3]>0 && bin[3]<=NBINSZ) && (bin[0]>0 && bin[0]<=NBINSE) && (bin[2]>0 && bin[2]<=NBINSZ) )
+    {
+      fCurrentCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]);
+      fCurrentCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]);    
+    }  
+  }
+    
+  printf("calling UnfoldWithBayesian\n");
+  Int_t success = UnfoldWithBayesian(fCurrentCorrelation, fCurrentRec, initialConditions, fUnfSpectrum, regPar, nIterations, kFALSE); 
+  
+  if ( success != 0)
+    return;
+
+  if (!determineError)
+  {
+    Printf("AliJetSpectrumUnfolding::ApplyBayesianMethod: WARNING: No errors calculated.");
+    return;
+  }
+
+  // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
+  // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
+  // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic
+  // error of the unfolding method itself.
+
+  TH2* maxError = (TH2*) fUnfSpectrum->Clone("maxError");
+  maxError->Reset();
+
+  TH2* normalizedResult = (TH2*) fUnfSpectrum->Clone("normalizedResult");
+  normalizedResult->Scale(1.0 / normalizedResult->Integral());
+
+  const Int_t kErrorIterations = 20;
+
+  printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
+
+  TH2* randomized = (TH2*) fCurrentRec->Clone("randomized");
+  TH2* result2 = (TH2*) fUnfSpectrum->Clone("result2");
+  for (Int_t n=0; n<kErrorIterations; ++n)
+  {
+    // randomize the content of clone following a poisson with the mean = the value of that bin
+    for (Int_t x=1; x<=randomized->GetNbinsX(); x++)
+      for (Int_t y=1; y<=randomized->GetNbinsY(); y++)
+      {
+        Float_t randomValue = fCurrentRec->GetBinContent(x,y);
+        TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",randomValue*0.25, randomValue*1.25);
+        poisson->SetParameters(randomValue,0.);
+        randomValue = poisson->GetRandom();   
+        //printf("%e --> %e\n", fCurrentRec->GetBinContent(x,y), (Double_t)randomValue);
+        randomized->SetBinContent(x, y, randomValue);
+        delete poisson;
+      }
+
+    result2->Reset();
+    if (UnfoldWithBayesian(fCurrentCorrelation, randomized, initialConditions, result2, regPar, nIterations) != 0)
+      return;
+
+    result2->Scale(1.0 / result2->Integral());
+
+    // calculate ratio
+    result2->Divide(normalizedResult);
+
+    //new TCanvas; result2->DrawCopy("HIST");
+
+    // find max. deviation
+    for (Int_t i=1; i<=result2->GetNbinsX(); i++)
+      for (Int_t j=1; j<=result2->GetNbinsY(); j++)
+        maxError->SetBinContent(i, j, TMath::Max(maxError->GetBinContent(i,j), TMath::Abs(1 - result2->GetBinContent(i,j))));
+  }
+  delete randomized;
+  delete result2;
+
+  for (Int_t i=1; i<=fUnfSpectrum->GetNbinsX(); i++)
+    for (Int_t j=1; j<=fUnfSpectrum->GetNbinsY(); j++)
+      fUnfSpectrum->SetBinError(i, j, fUnfSpectrum->GetBinError(i,j) + maxError->GetBinContent(i,j)*fUnfSpectrum->GetBinContent(i,j));
+
+  delete maxError;
+  delete normalizedResult;
+}
+
+//____________________________________________________________________
+Int_t AliJetSpectrumUnfolding::UnfoldWithBayesian(THnSparseF* correlation, TH2* measured, TH2* initialConditions, TH2* aResult, Float_t regPar, Int_t nIterations, Bool_t calculateErrors)
+{
+  //
+  // unfolds a spectrum
+  //
+
+  if (measured->Integral() <= 0)
+  {
+    Printf("AliJetSpectrumUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
+    return 1;
+  }
+  const Int_t NFilledBins = correlation->GetNbins();  
+  const Int_t kStartBin = 1;
+
+  const Int_t kMaxTZ = NBINSZ; // max true axis fragmentation function
+  const Int_t kMaxMZ = NBINSZ; // max measured axis fragmentation function
+  const Int_t kMaxTE = NBINSE; // max true axis energy
+  const Int_t kMaxME = NBINSE; // max measured axis energy
+  
+  printf("NbinsE=%d - NbinsZ=%d\n", NBINSE, NBINSZ);
+
+  // store information in arrays, to increase processing speed 
+  Double_t measuredCopy[kMaxME+1][kMaxMZ+1];
+  Double_t prior[kMaxTE+1][kMaxTZ+1];
+  Double_t errors[kMaxTE+1][kMaxTZ+1];
+  Double_t result[kMaxTE+1][kMaxTZ+1];
+
+  THnSparseF *inverseCorrelation;
+  inverseCorrelation = (THnSparseF*)correlation->Clone("inverseCorrelation");
+  inverseCorrelation->Reset();
+  
+  Float_t inputDistIntegral = 1;
+  if (initialConditions)
+  {
+    printf("Using different starting conditions...\n");   
+    inputDistIntegral = initialConditions->Integral();
+  }
+  Float_t measuredIntegral = measured->Integral();  
+  for (Int_t me=1; me<=kMaxME; me++)
+    for (Int_t mz=1; mz<=kMaxMZ; mz++)
+    {
+      // normalization of the measured spectrum
+      measuredCopy[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral;
+      errors[me][mz] = measured->GetBinError(me, mz) / measuredIntegral;
+      // pick prior distribution and normalize it
+      if (initialConditions)
+        prior[me][mz] = initialConditions->GetBinContent(me,mz) / inputDistIntegral;
+      else
+        prior[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral;
+    }
+
+  // unfold...
+  for (Int_t i=0; i<nIterations; i++)
+  {
+   // calculate Inverse Correlation Map from Bayes theorem:
+   // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
+   /*Float_t norm = 0;
+   for (Int_t me=1; me<=kMaxME; me++)
+      for (Int_t mz=1; mz<=kMaxMZ; mz++)
+      {
+        norm = 0;
+        for (Int_t te=kStartBin; te<=kMaxTE; te++)
+         for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++)
+         {
+           Int_t bin[4] = {te, me, tz, mz};
+            norm += correlation->GetBinContent(bin)*prior[te][tz];
+          }
+        if (norm > 0)
+          for (Int_t te = kStartBin; te <= kMaxTE; te++)
+           for (Int_t tz = kStartBin; tz <= kMaxTZ; tz++)
+           {
+             Int_t bin[4] = {te, me, tz, mz};
+             inverseCorrelation->SetBinContent(bin, correlation->GetBinContent(bin)*prior[te][tz]/norm );
+            }
+        //else
+          // inverse response set to '0' wich has been already done in line 2069
+      }*/
+    inverseCorrelation->Reset();     
+    Float_t norm[kMaxTE+2][kMaxTZ+2];
+    for (Int_t te=0; te<(kMaxTE+2); te++)
+      for (Int_t tz=0; tz<(kMaxTZ+2); tz++)
+        norm[te][tz]=0;
+    for (Int_t idx=0; idx<=correlation->GetNbins(); idx++)
+    {
+      Int_t bin[4];
+      Float_t binContent = correlation->GetBinContent(idx, bin);
+      if (bin[1]>0 && bin[1]<=NBINSE && bin[3]>0 && bin[3]<=NBINSZ &&
+          bin[0]>0 && bin[0]<=NBINSE && bin[2]>0 && bin[2]<=NBINSZ)
+        norm[bin[1]][bin[3]] += binContent*prior[bin[0]][bin[2]];
+    }
+    Float_t chi2Measured=0, diff;    
+    for (Int_t idx=0; idx<=correlation->GetNbins(); idx++)
+    {
+      Int_t bin[4];
+      Float_t binContent = correlation->GetBinContent(idx, bin);
+      if (norm[bin[1]][bin[3]]>0 && bin[1]>0 && bin[1]<=NBINSE && 
+          bin[3]>0 && bin[3]<=NBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=NBINSE && bin[2]<=NBINSZ)
+      {    
+        inverseCorrelation->SetBinContent(bin, binContent*prior[bin[0]][bin[2]]/norm[bin[1]][bin[3]]);
+        if (errors[bin[1]][bin[3]]>0)
+        {
+          diff = ((measuredCopy[bin[1]][bin[3]]-norm[bin[1]][bin[3]])/(errors[bin[1]][bin[3]]));
+          chi2Measured += diff*diff;
+        }   
+      }
+    }
+    
+    // calculate "generated" spectrum
+    for (Int_t te = kStartBin; te<=kMaxTE; te++)
+      for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++)
+      {
+        Float_t value = 0;
+        for (Int_t me=1; me<=kMaxME; me++)
+         for (Int_t mz=1; mz<=kMaxMZ; mz++)
+         {
+           Int_t bin[4] = {te, me, tz, mz};
+           value += inverseCorrelation->GetBinContent(bin)*measuredCopy[me][mz];
+          }
+        result[te][tz] = value;
+        //printf("%e\n", result[te][tz]);
+      }
+
+    // regularization (simple smoothing)
+    Float_t chi2LastIter = 0;
+    for (Int_t te=kStartBin; te<=kMaxTE; te++)
+      for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++)
+      {
+        Float_t newValue = 0;
+        // 0 bin excluded from smoothing
+        if (( te >(kStartBin+1) && te<(kMaxTE-1) ) && ( tz > (kStartBin+1) && tz<(kMaxTZ-1) ))
+        {
+          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.;
+
+          // weight the average with the regularization parameter
+          newValue = (1 - regPar) * result[te][tz] + regPar * average;
+        }
+        else
+          newValue = result[te][tz];
+        if (prior[te][tz]>1.e-5)
+        { 
+          diff = ((prior[te][tz]-newValue)/prior[te][tz]); 
+          chi2LastIter = diff*diff;
+        }  
+        prior[te][tz] = newValue;
+      }
+    //printf(" iteration %d - chi2LastIter = %e - chi2Measured = %e \n", i, chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ), chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ)); 
+    if (chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-6 && chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-3)
+      break;
+  } // end of iterations
+  
+  // propagate errors of the reconstructed distribution through the unfolding
+  for (Int_t te = kStartBin; te<=kMaxTE; te++)
+      for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++)
+      {
+        Float_t valueError = 0;
+        Float_t binError = 0;
+        for (Int_t me=1; me<=kMaxME; me++)
+         for (Int_t mz=1; mz<=kMaxMZ; mz++)
+         {
+           Int_t bin[4] = {te, me, tz, mz};
+           valueError += inverseCorrelation->GetBinContent(bin)*inverseCorrelation->GetBinContent(bin)*errors[me][mz]*errors[me][mz];
+          }
+        //if (errors[te][tz]!=0)printf("errors[%d][%d]=%e\n", te, tz, valueError);
+        aResult->SetBinContent(te, tz, prior[te][tz]);
+        aResult->SetBinError(te, tz, TMath::Sqrt(valueError));   
+      }
+
+  // ***********************************************************************************************************
+  // Calculate the covariance matrix, all arguments are taken from G. D'Agostini (p.6-8)
+  if (calculateErrors)
+  {
+    printf("Covariance matrix will be calculated... this will take a lot of time (>1 day) ;)\n");
+    
+    //Variables and Matrices that will be use along the calculation    
+    const Int_t binsV[4] = {NBINSE,NBINSE, NBINSZ, NBINSZ};
+    const Double_t LowEdgeV[4] = {axisLowerLimitE, axisLowerLimitE, axisLowerLimitZ, axisLowerLimitZ};
+    const Double_t UpEdgeV[4] = {axisUpperLimitE, axisUpperLimitE, axisUpperLimitZ, axisUpperLimitZ};
+    
+    const Double_t Ntrue = (Double_t)measured->Integral();
+    
+    THnSparseF *V = new THnSparseF("V","",4, binsV, LowEdgeV, UpEdgeV);
+    V->Reset();            
+    Double_t invCorrContent1, Nt;
+    Double_t invCorrContent2, v11, v12, v2;        
+    // calculate V1 and V2
+    for (Int_t idx1=0; idx1<=NFilledBins; idx1++)
+    {
+      printf("Covariance Matrix calculation: iteration idx1=%d of %d\n", idx1, NFilledBins);
+      for (Int_t idx2=0; idx2<=NFilledBins; idx2++)
+      {      
+        Int_t bin1[4];
+        Int_t bin2[4];
+        invCorrContent1 = inverseCorrelation->GetBinContent(idx1, bin1);
+        invCorrContent2 = inverseCorrelation->GetBinContent(idx2, bin2);        
+        v11=0; v12=0; v2=0;
+        if(bin1[0]>0 && bin1[0]<=NBINSE && bin1[1]>0 && bin1[1]<=NBINSE && 
+           bin1[2]>0 && bin1[2]<=NBINSZ && bin1[3]>0 && bin1[3]<=NBINSZ &&
+           bin2[0]>0 && bin2[0]<=NBINSE && bin2[1]>0 && bin2[1]<=NBINSE && 
+           bin2[2]>0 && bin2[2]<=NBINSZ && bin2[3]>0 && bin2[3]<=NBINSZ)
+        {   
+          if (bin1[1]==bin2[1] && bin1[3]==bin2[3])    
+            v11 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]
+                  *(1. - measuredCopy[bin2[1]][bin2[3]]/Ntrue);                       
+          else
+            v12 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]*
+                  measuredCopy[bin2[1]][bin2[3]]/Ntrue;
+          Nt = (Double_t)prior[bin2[0]][bin2[2]];
+          v2 = measuredCopy[bin1[1]][bin1[3]]*measuredCopy[bin2[1]][bin2[3]]*
+               invCorrContent1*invCorrContent2*
+               BayesUncertaintyTerms(inverseCorrelation, correlation, bin1, bin2, Nt);   
+          Int_t binV[4] = {bin1[0],bin2[0],bin1[2],bin2[2]};         
+          V->SetBinContent(binV,v11-v12 + v2);
+        }  
+      }
+    }    
+
+    for(Int_t te = 1; te<=NBINSE; te++)
+      for(Int_t tz = 1; tz<=NBINSZ; tz++)
+      {
+        Int_t binV[4] = {te,te,tz,tz}; 
+        aResult->SetBinError( te, tz, V->GetBinContent(binV) );
+      }            
+      
+    TFile* f = new TFile("Covariance_UnfSpectrum.root");
+    f->Open("RECREATE");
+    V->Write();
+    f->Close();    
+  }  
+  
+  return 0;
+
+}
+
+//____________________________________________________________________
+Double_t AliJetSpectrumUnfolding::BayesUncertaintyTerms(THnSparseF *M, THnSparseF *C, Int_t* binTM, Int_t* binTM1, Double_t Nt)
+{
+  //
+  // helper function for the covariance matrix of the bayesian method
+  //
+
+  Double_t result = 0;
+  Float_t term[9];
+  Int_t tmpBin[4], tmpBin1[4];
+  const Int_t nFilledBins = C->GetNbins();
+  if (Nt==0)
+    return 0;
+    
+  Float_t CorrContent;
+  Float_t InvCorrContent;
+
+  tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1];  tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3];
+  tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3];      
+  if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
+  {
+    if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
+      term[0] = BayesCov(M, C, tmpBin, tmpBin1)/
+                (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
+    term[2] = term[0]*M->GetBinContent(tmpBin1);
+  }            
+  else
+  {
+    term[0] = 0;
+    term[2] = 0;    
+  }
+              
+  tmpBin[0]=binTM1[0]; tmpBin[1]=binTM[1]; tmpBin[2]=binTM1[2]; tmpBin[3]=binTM[3];
+  tmpBin1[0]=binTM1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM1[2]; tmpBin1[3]=binTM1[3];      
+  if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
+    term[6] = BayesCov(M, C, tmpBin, tmpBin1)*
+              M->GetBinContent(tmpBin)/
+              (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));                    
+  else 
+    term[6] = 0;
+  
+  for(Int_t idx1=0; idx1<=nFilledBins; idx1++)
+  { 
+    Int_t bin1[4];
+    CorrContent    = C->GetBinContent(idx1, bin1); 
+    InvCorrContent = M->GetBinContent(idx1, bin1); 
+    if(bin1[0]>0 && bin1[0]<=NBINSE && bin1[1]>0 && bin1[1]<=NBINSE && 
+       bin1[2]>0 && bin1[2]<=NBINSZ && bin1[3]>0 && bin1[3]<=NBINSZ)
+    {
+      tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3];
+      tmpBin1[0]=binTM[0]; tmpBin1[1]=bin1[1];  tmpBin1[2]=binTM[2]; tmpBin1[3]=bin1[3];      
+      if (C->GetBinContent(tmpBin)!=0 &&
+          binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
+        term[1] = BayesCov(M, C, tmpBin, tmpBin1)/C->GetBinContent(tmpBin);
+      else
+        term[1] = 0;
+
+      tmpBin[0] =binTM[0]; tmpBin[1] =bin1[1];   tmpBin[2] =binTM[2]; tmpBin[3] =bin1[3];
+      tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3];      
+      if (C->GetBinContent(tmpBin1)!=0)
+      {
+        if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2])
+          term[3] = BayesCov(M, C, tmpBin, tmpBin1)/
+                    C->GetBinContent(tmpBin1);
+        term[5] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin1)/
+                  C->GetBinContent(tmpBin1);
+      }            
+      else
+      {
+        term[3] = 0;
+        term[5] = 0;
+      }  
+   
+      tmpBin[0] =binTM1[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM1[2]; tmpBin[3] =binTM[3];
+      tmpBin1[0]=binTM1[0]; tmpBin1[1]=bin1[1];  tmpBin1[2]=binTM1[2]; tmpBin1[3]=bin1[3];      
+      if (C->GetBinContent(tmpBin)!=0)
+        term[7] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin)/
+                  C->GetBinContent(tmpBin);
+      else
+        term[7] = 0;
+
+      tmpBin[0] =bin1[0]; tmpBin[1] =binTM[1];  tmpBin[2] =bin1[2]; tmpBin[3] =binTM[3];
+      tmpBin1[0]=bin1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=bin1[2]; tmpBin1[3]=binTM1[3];      
+      if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0)
+        term[8] = BayesCov(M, C, tmpBin, tmpBin1)*
+                  M->GetBinContent(tmpBin)*M->GetBinContent(tmpBin)/
+                  (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1));
+      else 
+        term[8] = 0;
+                       
+      for (Int_t i=0; i<9; i++)
+        result += term[i]/Nt;                    
+    }          
+  }
+   
+  return result;
+}
+
+//____________________________________________________________________
+Double_t AliJetSpectrumUnfolding::BayesCov(THnSparseF *M, THnSparseF *correlation, Int_t* binTM, Int_t* bin1)
+{
+  Double_t result, result1, result2, result3;
+  
+  if (binTM[0]==bin1[0] && binTM[2]==bin1[2])
+  {
+    if (correlation->GetBinContent(bin1)!=0) 
+      result1 = 1./correlation->GetBinContent(bin1);
+    else 
+      result1 = 0;
+    result2 = 1.;
+  }
+  else
+  {
+    result1 = 0;
+    result2 = 0;
+  }  
+    
+  if (binTM[1]==bin1[1] && binTM[3]==bin1[3])
+  {
+    Int_t tmpbin[4] = {bin1[0], binTM[1], bin1[2], binTM[3]};
+    if(correlation->GetBinContent(tmpbin)!=0)
+      result3 = M->GetBinContent(tmpbin)/correlation->GetBinContent(tmpbin);
+    else 
+      result3 = 0;
+  }
+  else
+  {
+    result1 = 0;
+    result3 = 0;
+  }
+    
+  return result = result1 + result2 + result3;
+}
+
+//____________________________________________________________________
+TH2F* AliJetSpectrumUnfolding::CalculateRecSpectrum(TH2* inputGen)
+{
+  // runs the distribution given in inputGen through the correlation histogram identified by
+  // fCorrelation and produces a reconstructed spectrum
+
+  if (!inputGen)
+    return 0;
+
+  // normalize to convert number of events into probability
+  /*for (Int_t te=1; te<=NBINSE; te++)
+    for (Int_t tz=1; tz<=NBINSZ; tz++)
+    {
+       Int_t bin[4];
+       Float_t sum=0.;
+       for (Int_t me = 1; me<=NBINSE; me++)
+         for (Int_t mz = 1; mz<=NBINSZ; mz++)
+         {
+           bin[0] = te; bin[1] = me; 
+           bin[2] = tz; bin[3] = mz;           
+           sum += fCorrelation[correlationMap]->GetBinContent(bin);
+         }
+       if (sum > 0.)
+         for (Int_t me = 1; me<=NBINSE; me++)
+           for (Int_t mz = 1; mz<=NBINSZ; mz++)
+           {
+             bin[0] = te; bin[1] = me; 
+             bin[2] = tz; bin[3] = mz;           
+             fCorrelation[correlationMap]->SetBinContent(bin, fCorrelation[correlationMap]->GetBinContent(bin)/sum);
+            fCorrelation[correlationMap]->SetBinError(bin, fCorrelation[correlationMap]->GetBinError(bin)/sum);
+          }
+    }*/  
+  // normalize to convert number of events into probability (the following loop is much faster)
+  Float_t sum[NBINSE+2][NBINSZ+2];
+  memset(sum,0,sizeof(Float_t)*(NBINSE+2)*(NBINSZ+2));
+
+  for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++)
+  {
+    Int_t bin[4];
+    Float_t binContent = fCorrelation->GetBinContent(idx, bin);
+    if (bin[1]>0 && bin[1]<=NBINSE && bin[3]>0 && bin[3]<=NBINSZ){
+      sum[bin[0]][bin[2]] += binContent; 
+    }
+  }
+  
+  for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++)
+  {
+    Int_t bin[4];
+    Float_t binContent = fCorrelation->GetBinContent(idx, bin);
+    Float_t binError   = fCorrelation->GetBinError(bin);
+    if (sum[bin[0]][bin[2]]>0 && bin[1]>0 && bin[1]<=NBINSE && 
+        bin[3]>0 && bin[3]<=NBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=NBINSE && bin[2]<=NBINSZ) 
+    {
+      fCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]);
+      fCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]); 
+    }  
+  }
+
+  TH2F* target = dynamic_cast<TH2F*> (fRecSpectrum->Clone(Form("reconstructed_%s", inputGen->GetName())));
+  target->Reset();
+
+  for (Int_t me=1; me<=NBINSE; ++me)
+    for (Int_t mz=1; mz<=NBINSZ; ++mz)
+    {
+      Float_t measured = 0;
+      Float_t error = 0;
+
+      for (Int_t te=1; te<=NBINSE; ++te)
+        for (Int_t tz=1; tz<=NBINSZ; ++tz)
+        {
+         Int_t bin[4] = {te, me, tz, mz};
+          measured += inputGen->GetBinContent(te,tz) * fCorrelation->GetBinContent(bin);
+          error += inputGen->GetBinError(te,tz) * fCorrelation->GetBinContent(bin);
+        }
+      target->SetBinContent(me, mz, measured);
+      target->SetBinError(me, mz, error);
+    }
+
+  return target;
+}
+
+//__________________________________________________________________________________________________
+void AliJetSpectrumUnfolding::SetGenRecFromFunc(TF2* inputGen)
+{
+  // uses the given function to fill the input Generated histogram and generates from that
+  // the reconstructed histogram by applying the response histogram
+  // this can be used to evaluate if the methods work indepedently of the input
+  // distribution
+
+  if (!inputGen)
+    return;
+
+  TH2F* histtmp = new TH2F("histtmp", "tmp", EBINNING, ZBINNING);
+  TH2F* gen  = fGenSpectrum;
+
+  histtmp->Reset();
+  gen->Reset();
+
+  histtmp->FillRandom(inputGen->GetName(), NEVENTS);
+
+  for (Int_t i=1; i<=gen->GetNbinsX(); ++i)
+    for (Int_t j=1; j<=gen->GetNbinsY(); ++j)
+    {
+      gen->SetBinContent(i, j, histtmp->GetBinContent(i,j));
+      gen->SetBinError(i, j, histtmp->GetBinError(i,j));
+    }
+
+  delete histtmp;
+
+  //new TCanvas;
+  //gStyle->SetPalette(1);
+  //gPad->SetLogz();
+  //gen->Draw("COLZ");
+
+
+  TH2 *recsave = fRecSpectrum;
+
+  fRecSpectrum = CalculateRecSpectrum(gen);
+  fRecSpectrum->SetName(recsave->GetName());
+  delete recsave;
+
+  return;
+}
+//________________________________________________________________________________________
diff --git a/PWG4/JetTasks/AliJetSpectrumUnfolding.h b/PWG4/JetTasks/AliJetSpectrumUnfolding.h
new file mode 100644 (file)
index 0000000..e30015d
--- /dev/null
@@ -0,0 +1,89 @@
+#ifndef ALIJETSPECTRUMUNFOLDING_H
+#define ALIJETSPECTRUMUNFOLDING_H
+
+#include "TNamed.h"
+//
+// class that contains the correction matrix and the functions for
+// correction the jet spectrum
+// implements 2-dim bayesian method
+//
+
+class TH1;
+class TH2;
+class TH3;
+class TH1F;
+class TH2F;
+class TH3F;
+class TF1;
+class TF2;
+class TCollection;
+
+#include <TMatrixD.h>
+#include <TVectorD.h>
+#include <THnSparse.h>
+
+class AliJetSpectrumUnfolding : public TNamed {
+  public:
+    AliJetSpectrumUnfolding();
+    AliJetSpectrumUnfolding(const Char_t* name, const Char_t* title);
+    virtual ~AliJetSpectrumUnfolding();
+
+    virtual Long64_t Merge(TCollection* list);
+
+    Bool_t LoadHistograms(const Char_t* dir = 0);
+    void SaveHistograms();
+    void DrawHistograms();
+    void DrawComparison(const char* name, TH2* genSpec);
+
+    void SetBayesianParameters(Float_t smoothing, Int_t nIterations);
+    void ApplyBayesianMethod(Float_t regPar = 1, Int_t nIterations = 100, TH2* initialConditions = 0, Bool_t determineError = kTRUE);
+
+    TH2F* GetRecSpectrum()      { return fRecSpectrum; }
+    TH2F* GetGenSpectrum()      { return fGenSpectrum; }
+    TH2F* GetUnfSpectrum()      { return fUnfSpectrum; }
+    THnSparseF* GetCorrelation(){ return fCorrelation; }
+
+    void SetRecSpectrum(TH2F* hist)      { fRecSpectrum  = hist; }
+    void SetGenSpectrum(TH2F* hist)      { fGenSpectrum  = hist; }
+    void SetUnfSpectrum(TH2F* hist)      { fUnfSpectrum  = hist; }
+    void SetCorrelation(THnSparseF* hist){ fCorrelation  = hist; }
+
+    void SetGenRecFromFunc(TF2* inputGen);
+    TH2F* CalculateRecSpectrum(TH2* inputGen);
+
+    static void NormalizeToBinWidth(TH2* hist);
+
+    void GetComparisonResults(Float_t* gen = 0, Int_t* genLimit = 0, Float_t* residuals = 0, Float_t* ratioAverage = 0) const;
+
+  protected:
+    void SetupCurrentHists(Bool_t createBigBin);
+
+    static Double_t BayesCov(THnSparseF* M, THnSparseF* correlation, Int_t* binTM, Int_t* binTM1);
+    static Double_t BayesUncertaintyTerms(THnSparseF *M, THnSparseF *C, Int_t* binTM, Int_t* binTM1, Double_t Nt);
+    static Int_t UnfoldWithBayesian(THnSparseF* correlation, TH2* measured, TH2* initialConditions, TH2* aResult, Float_t regPar, Int_t nIterations, Bool_t calculateErrors = kFALSE);
+
+    static Float_t fgBayesianSmoothing;             //! smoothing parameter (0 = no smoothing)
+    static Int_t   fgBayesianIterations;            //! number of iterations in Bayesian method
+
+    TH2F* fCurrentRec;
+    THnSparseF* fCurrentCorrelation;
+
+    TH2F* fRecSpectrum;
+    TH2F* fGenSpectrum;
+    TH2F* fUnfSpectrum;
+    THnSparseF* fCorrelation;
+
+    Float_t fLastChi2MC;        //! last Chi2 between MC and unfolded ESD (calculated in DrawComparison)
+    Int_t   fLastChi2MCLimit;   //! bin where the last chi2 breached a certain threshold
+    Float_t fLastChi2Residuals; //! last Chi2 of the ESD and the folded unfolded ESD (calculated in DrawComparison)
+    Float_t fRatioAverage;      //! last average of |ratio-1| where ratio = unfolded / mc (bin 2..150)
+
+ private:
+    AliJetSpectrumUnfolding(const AliJetSpectrumUnfolding&);
+    AliJetSpectrumUnfolding& operator=(const AliJetSpectrumUnfolding&);
+
+  ClassDef(AliJetSpectrumUnfolding, 2);
+};
+
+#endif
+
index 5599fbd..4ef44c8 100644 (file)
@@ -8,6 +8,7 @@
 #pragma link C++ class AliAnalysisTaskJetSpectrum+;
 #pragma link C++ class AliAnalysisHelperJetTasks+;
 #pragma link C++ class AliAnaESDSpectraQA+;
+#pragma link C++ class AliJetSpectrumUnfolding+;
 
 
 #endif
index 8f11dd7..ebf8ecb 100644 (file)
@@ -1,6 +1,6 @@
 #-*- Mode: Makefile -*-
 
-SRCS = JetTasks/AliAnalysisTaskUE.cxx JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnaESDSpectraQA.cxx
+SRCS = JetTasks/AliAnalysisTaskUE.cxx JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnaESDSpectraQA.cxx JetTasks/AliJetSpectrumUnfolding.cxx
 
 HDRS:= $(SRCS:.cxx=.h) 
 
index 3505bb9..167907e 100644 (file)
@@ -20,7 +20,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
     Int_t iESDfilter     = 1;  // Only active if iAODanalysis=0
     Int_t iJETAN         = 1;
     Int_t iJETANAOD      = 0;
-    Int_t iJETANMC       = 0;
+    Int_t iJETANMC       = 1;
     Int_t iDIJETAN       = 0;
     Int_t iPWG4SPECTRUM  = 1;
     Int_t iPWG4UE        = 0;
@@ -35,10 +35,11 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
     if (iESDfilter) iAODhandler=1;
 
     // Dataset from CAF
-    TString dataset = "/PWG4/arian/jetjet15-50";
-    // CKB quick hack
+    //    TString dataset = "/PWG4/kleinb/LHC08v_jetjet15-50";
+    TString dataset = "/PWG4/kleinb/LHC08r_jetjet50";
+    // CKB quick hack for local analysis
     gROOT->LoadMacro("CreateESDChain.C");
-    TChain *chain = CreateESDChain("jetjet15-50.txt",10);
+    TChain *chain = CreateESDChain("jetjet15-50.txt",1000);
 
 
     printf("==================================================================\n");
@@ -69,11 +70,11 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
 
 
     // Reset user processes if CAF if not responding anymore
-    //TProof::Reset("lxb6046"); 
+    // TProof::Reset("alicecaf"); 
     // One may enable a different ROOT version on CAF
 
-    const char* proofNode = "lxb6046";
-    //    const char* proofNode = "kleinb@localhost";
+    //    const char* proofNode = "alicecaf";
+    const char* proofNode = "localhost";
 
 
 
@@ -81,7 +82,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
     // Connect to proof
     if(bPROOF){
       TProof::Mgr(proofNode)->ShowROOTVersions();
-      TProof::Mgr(proofNode)->SetROOTVersion("v5-21-01-alice_dbg");
+      //      TProof::Mgr(proofNode)->SetROOTVersion("v5-21-01-alice_dbg");
       TProof::Open(proofNode); 
 
       // Clear packages if changing ROOT version on CAF or local
@@ -162,10 +163,9 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
     if (iAODhandler) {
        // AOD output handler
        AliAODHandler* aodHandler   = new AliAODHandler();
-       mgr->SetOutputEventHandler(aodHandler);
+       aodHandler->SetFillAOD(kFALSE);
+       mgr->SetOutputEventHandler(aodHandler);       
        aodHandler->SetOutputFileName(Form("AliAODs_CKB_%07d-%07d.root",nOffset,nOffset+nEvents));
-       //       aodHandler->SetSelectAll(kFALSE);
-       //       aodHandler->SetCreateNonStandardAOD();
        cout_aod = mgr->CreateContainer("cAOD", TTree::Class(),
                                       AliAnalysisManager::kOutputContainer, "default");
        cout_aod->SetSpecialOutput();
@@ -211,7 +211,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
        mgr->AddTask(jetana);
        // Output histograms list for jet analysis                       
        AliAnalysisDataContainer *cout_jet = mgr->CreateContainer("jethist", TList::Class(),
-                                                             AliAnalysisManager::kOutputContainer, "jethist.root");
+                                                             AliAnalysisManager::kOutputContainer,"jethist.root");
        // Dummy AOD output container for jet analysis (no client yet)
        c_aodjet = mgr->CreateContainer("cAODjet", TTree::Class(),
                            AliAnalysisManager::kExchangeContainer);
@@ -285,7 +285,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
       //      if(iAODanalysis)pwg4spec->SetAODInput(kTRUE);
       pwg4spec->SetDebugLevel(11); 
       //      pwg4spec->SetBranchRec("jetsMC"); 
-      //      pwg4spec->SetBranchGen("jetsAOD"); 
+      //      pwg4spec->SetBranchGen("jetsMC"); 
       mgr->AddTask(pwg4spec);
 
       AliAnalysisDataContainer *coutput1_Spec = mgr->CreateContainer("pwg4spec", TList::Class(),AliAnalysisManager::kOutputContainer, "histos_pwg4spec.root");
@@ -296,8 +296,6 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
       mgr->ConnectInput  (pwg4spec,  0, cinput);    
       // mgr->ConnectInput  (pwg4spec,  0, c_aodjet);    
       mgr->ConnectOutput (pwg4spec,  0, c_aodSpec );
-      //       mgr->ConnectOutput (pwg4spec,  1, coutput1_Spec );
-       
       mgr->ConnectOutput (pwg4spec,  1, coutput1_Spec );
     }   
 
diff --git a/PWG4/macros/runJetSpectrumUnfolding.C b/PWG4/macros/runJetSpectrumUnfolding.C
new file mode 100644 (file)
index 0000000..84d1958
--- /dev/null
@@ -0,0 +1,411 @@
+//
+// script with functions to use AliJetSpectrumUnfolding class
+//
+
+
+void loadlibs(){
+  // load all the needed libs to run wit root only
+
+  gSystem->Load("libTree");
+  gSystem->Load("libPhysics");
+  gSystem->Load("libHist");
+  gSystem->Load("libVMC");
+  gSystem->Load("libSTEERBase");
+  gSystem->Load("libAOD");
+  gSystem->Load("libESD");
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+  gSystem->Load("libJETAN");
+  gSystem->Load("libPWG4JetTasks");
+
+
+}
+
+
+
+void draw(const char* fileName = "unfolded_pwg4spec.root", const char* folder = "unfolding", Bool_t proj = kTRUE)
+{
+
+  loadlibs();
+
+  AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
+
+  TFile::Open(fileName);
+  jetSpec->LoadHistograms();
+  
+  
+  if (proj)
+  {
+    canvas1 = new TCanvas("Response Map Projection", "Response Map Projection", 500, 500);
+    canvas1->Divide(2);
+  
+    Int_t style = 1;
+    const Int_t NRGBs = 5;
+    const Int_t NCont = 500;
+
+    Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+    Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
+    Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
+    Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
+    TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
+    gStyle->SetNumberContours(NCont);
+    gStyle->SetPalette(style);
+
+    canvas1->cd(1);
+    gPad->SetLogz();
+    h2 = (jetSpec->GetCorrelation())->Projection(2,3);
+    h2->SetXTitle("z^{lp}_{gen}");
+    h2->SetYTitle("z^{lp}_{rec}");    
+    h2->DrawCopy("colz");
+  
+    canvas1->cd(2);
+    gPad->SetLogz();  
+    h3 = (jetSpec->GetCorrelation())->Projection(0,1);
+    h3->SetXTitle("E^{jet}_{gen} [GeV]");
+    h3->SetYTitle("E^{jet}_{rec} [GeV]");    
+    h3->DrawCopy("colz");
+  }
+  jetSpec->DrawComparison("Draw_unfolded_pwg4spec", jetSpec->GetGenSpectrum());
+
+  return;
+}
+
+//________________________________________________________________________________________________________________
+void unfold(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameUnf = "unfolded_pwg4spec.root")
+{
+  // function to load jet spectra from the output file of the task AliAnalysisTaskJetSpectrum
+  // to do the unfolding
+
+  loadlibs();
+
+  AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
+
+  TFile::Open(fileNameRec);
+  jetSpec->LoadHistograms();
+
+  TFile::Open(fileNameGen);
+  TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
+  jetSpec->SetGenSpectrum(hist);
+
+  jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0);
+  // last parameter = calculateErrors  <- this method to calculate the errors takes a lot of time
+   
+  TFile* file = TFile::Open(fileNameUnf,"RECREATE");
+  jetSpec->SaveHistograms();
+  file->Close();
+}
+
+//___________________________________________________________________________
+void buildSpectra(Int_t caseNo, const char* inFile, const char* outFile)
+{
+  // function to test 2D Bayesian unfolding with other spectra
+  // build from a function
+
+  loadlibs();
+
+
+
+  AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
+
+  TFile::Open(inFile);
+  jetSpec->LoadHistograms("unfolding");
+  
+  TCanvas *c1 = new TCanvas();TH2 *h2 = (jetSpec->GetCorrelation())->Projection(0,3);
+  h2->DrawCopy("colz");
+  c1->Update();
+  if(getchar()=='q')return;
+  
+
+  switch (caseNo)
+  {
+    case 1: func = new TF2("func", "501-x-y"); break;
+    case 2: func = new TF2("func", "1000 * 1/(y+x+1)"); break;
+    case 3: func = new TF2("func", "1./(y*pow(x,5.7))"); break;
+    case 4: func = new TF2("func", "exp(-0.1*x - 0.1*y)"); break;
+    case 5: func = new TF2("func", "x*x + (y**3)/100."); break;
+    case 6: func = new TF2("func", "1000*y*exp(-0.1*x)"); break;
+    case 7: func = new TF2("func", "exp(-((x-100.)/(0.3*x))**2 - ((y-0.6)/(0.8*y))**2)"); break;
+    default: return;
+  }
+
+  //new TCanvas; func->Draw();
+
+  jetSpec->SetGenRecFromFunc(func);
+              
+  TFile* file = TFile::Open(outFile,"RECREATE");
+  jetSpec->SaveHistograms();
+  
+  h2 = (jetSpec->GetCorrelation())->Projection(0,3);
+  h2->DrawCopy("colz");
+  c1->Update();
+  file->Close();
+
+  //new TCanvas; jetSpec->GetRecSpectrum()->DrawCopy("colz");
+  //new TCanvas; jetSpec->GetGenSpectrum()->DrawCopy("colz");
+}
+
+//_____________________________________________________________________________________________
+void buildResponseMap(const char* fileName = "responseMap.root")
+{
+  // function to build a Response Map with a gaussian distribution
+  loadlibs();
+
+  AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
+
+  TF2* func = new TF2("func", "exp(-((x-[0])/[1])**2 - ((y-[2])/[3])**2)");
+  
+  bool bPerfect = false;
+  
+  Double_t var[4];
+  Float_t sigmax, sigmay;
+  for (Int_t tx=1; tx<=jetSpec->GetCorrelation()->GetAxis(0)->GetNbins(); tx++)
+    for (Int_t ty=1; ty<=jetSpec->GetCorrelation()->GetAxis(2)->GetNbins(); ty++)
+    {
+      var[0] = jetSpec->GetCorrelation()->GetAxis(0)->GetBinCenter(tx);
+      var[2] = jetSpec->GetCorrelation()->GetAxis(2)->GetBinCenter(ty);
+      sigmax = 0.2*var[0];
+      sigmay = 0.2*var[2];
+      func->SetParameters(var[0],sigmax,var[2],sigmay);
+      for (Int_t mx=1; mx<=jetSpec->GetCorrelation()->GetAxis(1)->GetNbins(); mx++)
+        for (Int_t my=1; my<=jetSpec->GetCorrelation()->GetAxis(3)->GetNbins(); my++)
+        {
+          var[1] = jetSpec->GetCorrelation()->GetAxis(1)->GetBinCenter(mx);
+         var[3] = jetSpec->GetCorrelation()->GetAxis(3)->GetBinCenter(my);
+
+
+         if(bPerfect){
+           if (var[1]==var[0] && var[3]==var[2])
+             jetSpec->GetCorrelation()->Fill(var,1);
+         }
+         else {
+           // cut at  sigma
+           if (TMath::Abs(var[1]-var[0]) < sigmax || TMath::Abs(var[3]-var[2]) < sigmay)
+             jetSpec->GetCorrelation()->Fill(var,func->Eval(var[1],var[3]));
+         }
+        }
+    }
+
+
+  TFile* file = TFile::Open(fileName,"RECREATE");
+  jetSpec->SaveHistograms();
+  file->Close();
+}
+
+//_____________________________________________________________________________
+void StatisticalUncertainties(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameOut = "Uncertainties2DBayesUnf.root")
+{
+  // This function gives the statistical uncertainties due to the 2D Bayesian Unfolding
+
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+  gSystem->Load("libJETAN");
+  gSystem->Load("libPWG4JetTasks");
+
+  AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
+
+  TFile::Open(fileNameRec);
+  jetSpec->LoadHistograms();
+
+  TFile::Open(fileNameGen);
+  TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
+  jetSpec->SetGenSpectrum(hist);
+    
+  // create sigma histogram
+  TH2F* sigma = (TH2F*)(jetSpec->GetGenSpectrum())->Clone("sigma");
+  sigma->Reset();
+    
+  THnSparseF* hcorr = (THnSparseF*)(jetSpec->GetCorrelation())->Clone();
+  TH2F*       htrue = (TH2F*)(jetSpec->GetGenSpectrum())->Clone();  
+  TH2F*       hmeas = (TH2F*)(jetSpec->GetRecSpectrum())->Clone();    
+  TH2F*       hunfo = (TH2F*)(jetSpec->GetUnfSpectrum())->Clone();      
+  
+  Int_t nIterations = 1000;
+  Float_t binContent;
+  for(Int_t i=0; i<nIterations; i++)
+  {
+    printf("iteration = %d\n", i);
+    // reset histograms
+    jetSpec->GetRecSpectrum()->Reset();
+    jetSpec->GetGenSpectrum()->Reset();
+    jetSpec->GetCorrelation()->Reset();
+    jetSpec->GetUnfSpectrum()->Reset();
+  
+    THnSparseF* tmpcorr = (THnSparseF*)hcorr->Clone("tmpcorr"); 
+    TH2F*       tmptrue = (TH2F*)htrue->Clone("tmptrue");  
+    
+    jetSpec->SetGenSpectrum(tmptrue);
+    jetSpec->SetCorrelation(tmpcorr);
+  
+    // randomize reconstructed distribution
+    for (Int_t me=1; me<=hmeas->GetNbinsX(); me++)
+      for (Int_t mz=1; mz<=hmeas->GetNbinsY(); mz++)
+      {
+        binContent = hmeas->GetBinContent(me,mz);
+        if (binContent>0)
+        {
+          TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",binContent*0.25, binContent*1.25);
+          poisson->SetParameters(binContent,0.);
+          binContent = poisson->GetRandom();
+          delete poisson;
+        }  
+        jetSpec->GetRecSpectrum()->SetBinContent(me,mz, binContent);
+      } 
+        
+    // unfold
+    jetSpec->ApplyBayesianMethod(0.2, 20, 0, 0);
+    
+    // calculate sigma^2
+    for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
+      for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
+      {
+        if (htrue->GetBinContent(te,tz)!=0)
+        {
+          binContent = (jetSpec->GetUnfSpectrum()->GetBinContent(te,tz) -
+                        htrue->GetBinContent(te,tz))/htrue->GetBinContent(te,tz);
+          binContent *= binContent;
+          sigma->SetBinContent(te,tz, binContent + sigma->GetBinContent(te,tz));
+        }  
+      } 
+  }
+  // calculate sigma   
+  for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
+    for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
+    {
+      binContent = sigma->GetBinContent(te,tz);
+      binContent = TMath::Sqrt(binContent/(Float_t)nIterations);
+      sigma->SetBinContent(te,tz, binContent);
+    } 
+          
+  const Int_t NRGBs = 5;
+  const Int_t NCont = 500;
+
+  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+  Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
+  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
+  Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
+  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
+  gStyle->SetNumberContours(NCont);   
+  gStyle->SetPalette(1);  
+
+  new TCanvas();
+  gPad->SetLogz();
+  sigma->SetTitle("#sigma((U_{R} - U)/U)");
+  sigma->SetXTitle("E^{jet} [GeV]");
+  sigma->SetYTitle("z^{lp}");  
+  sigma->DrawCopy();
+  
+  TFile* file = TFile::Open(fileNameOut,"RECREATE");
+  sigma->Write();
+  file->Close();   
+}
+
+//_______________________________________________________________________________________________________________
+void FillSpecFromFile(const char* fileNameSpec = "histos_pwg4spec.root")
+{
+  // This functions avoids problems when the number of bins or the bin limits
+  // in the histograms of the AliJetSpectrumUnfolding and AliAnalysisTaskJetSpectrum classes
+  // are different.
+
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+  gSystem->Load("libJETAN");
+  gSystem->Load("libPWG4JetTasks");
+
+  file = new TFile(fileNameSpec,"read");
+  tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
+  THnSparseF *fhCorrelation  = (THnSparseF*)(tlist->FindObject("fhCorrelation_less5tracks"));
+  THnSparseF *fhCorrelation2 = (THnSparseF*)(tlist->FindObject("fhCorrelation_5to10tracks"));
+  THnSparseF *fhCorrelation3 = (THnSparseF*)(tlist->FindObject("fhCorrelation_more10tracks"));
+  TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fhEGenZGen"));
+  TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fhERecZRec"));  
+
+  fhCorrelation->Add(fhCorrelation2, 1);
+  fhCorrelation->Add(fhCorrelation3, 1);
+
+  delete fhCorrelation2;
+  delete fhCorrelation3;
+
+  AliJetSpectrumUnfolding *jetSpec = new AliJetSpectrumUnfolding("unfolding","unfolding");
+
+  // generated jets (true distribution)
+  for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++)
+    for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++)
+    {
+       Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te);
+       Float_t  z = fhEGenZGen->GetYaxis()->GetBinCenter(tz);
+       Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej);
+       Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z);
+       Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz);
+       Float_t err  = jetSpec->GetGenSpectrum()->GetBinError(bine,binz);
+       jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
+       jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
+    }
+  file = TFile::Open("gen_pwg4spec.root", "RECREATE");
+  jetSpec->SaveHistograms();
+  file->Close();
+  jetSpec->GetGenSpectrum()->Reset();
+  printf("true distribution has been set\n");
+
+  // reconstructed jets (measured distribution)
+  for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++)
+    for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++)
+    {
+       Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me);
+       Float_t   zm = fhERecZRec->GetYaxis()->GetBinCenter(mz);
+       Int_t bine   = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec);
+       Int_t binz   = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm);
+       Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz);
+       Float_t err  = jetSpec->GetRecSpectrum()->GetBinError(bine, binz);
+       jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz));
+       jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
+    }
+
+  // Response function
+  jetSpec->GetCorrelation()->Reset();
+  jetSpec->GetCorrelation()->Sumw2();
+        
+  for (Int_t idx=1; idx<=fhCorrelation->GetNbins(); idx++)
+  {
+    //printf("idx: %d\n",idx);
+    Double_t var[4];
+    Int_t bin[4];
+    Float_t BinContent = fhCorrelation->GetBinContent(idx, bin);
+    var[0] = fhCorrelation->GetAxis(0)->GetBinCenter(bin[0]);
+    var[1] = fhCorrelation->GetAxis(1)->GetBinCenter(bin[1]);
+    var[2] = fhCorrelation->GetAxis(2)->GetBinCenter(bin[2]);
+    var[3] = fhCorrelation->GetAxis(3)->GetBinCenter(bin[3]);
+    bin[0] = jetSpec->GetCorrelation()->GetAxis(0)->FindBin(var[0]);
+    bin[1] = jetSpec->GetCorrelation()->GetAxis(1)->FindBin(var[1]);    
+    bin[2] = jetSpec->GetCorrelation()->GetAxis(2)->FindBin(var[2]);
+    bin[3] = jetSpec->GetCorrelation()->GetAxis(3)->FindBin(var[3]);        
+    Float_t cont = jetSpec->GetCorrelation()->GetBinContent(bin);
+    Float_t err  = jetSpec->GetCorrelation()->GetBinError(bin);
+    jetSpec->GetCorrelation()->SetBinContent(bin, cont + fhCorrelation->GetBinContent(idx));
+    jetSpec->GetCorrelation()->SetBinError(bin, err + fhCorrelation->GetBinError(idx));
+  }
+
+  file = TFile::Open("rec_pwg4spec.root", "RECREATE");
+  jetSpec->SaveHistograms();
+  file->Close();
+
+  printf("reconstructed distribution has been set\n");    
+  printf("response map has been set\n");
+  
+}
+
+
+
+
+void correct(){
+  // simple steering to correct a given distribution;
+  loadlibs();
+  FillSpecFromFile("/home/ydelgado/pcalice014.cern.ch/macros/jets/CAFdata/histos_pwg4spec.root");
+
+  char name[100];
+  sprintf(name, "unfolded_pwg4spec.root");
+
+  unfold("gen_pwg4spec.root", "unfolding", "rec_pwg4spec.root", name);
+  //draw(name, "unfolding", 1); 
+
+}