New class to handle multi-dimensional unfolding + user macros (test/)
authorrvernet <rvernet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 8 Apr 2009 18:10:09 +0000 (18:10 +0000)
committerrvernet <rvernet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 8 Apr 2009 18:10:09 +0000 (18:10 +0000)
CORRFW/AliCFUnfolding.cxx [new file with mode: 0644]
CORRFW/AliCFUnfolding.h [new file with mode: 0644]
CORRFW/CORRFWLinkDef.h
CORRFW/libCORRFW.pkg
CORRFW/test/AliCFTaskForUnfolding.C [new file with mode: 0644]
CORRFW/test/AliCFTaskForUnfolding.cxx [new file with mode: 0644]
CORRFW/test/AliCFTaskForUnfolding.h [new file with mode: 0644]
CORRFW/test/testUnfolding.C [new file with mode: 0644]

diff --git a/CORRFW/AliCFUnfolding.cxx b/CORRFW/AliCFUnfolding.cxx
new file mode 100644 (file)
index 0000000..88cec41
--- /dev/null
@@ -0,0 +1,474 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//--------------------------------------------------------------------//
+//                                                                    //
+// AliCFUnfolding Class                                               //
+// Class to handle general unfolding procedure                        // 
+// For the moment only bayesian unfolding is supported                //
+// The next steps are to add chi2 minimisation and weighting methods  //
+//                                                                    //
+// Author : renaud.vernet@cern.ch                                     //
+//--------------------------------------------------------------------//
+
+
+#include "AliCFUnfolding.h"
+#include "TMath.h"
+#include "TAxis.h"
+#include "AliLog.h"
+
+ClassImp(AliCFUnfolding)
+
+//______________________________________________________________
+
+AliCFUnfolding::AliCFUnfolding() :
+  TNamed(),
+  fResponse(0x0),
+  fPrior(0x0),
+  fOriginalPrior(0x0),
+  fEfficiency(0x0),
+  fMeasured(0x0),
+  fMaxNumIterations(0),
+  fNVariables(0),
+  fMaxChi2(0),
+  fUseSmoothing(kFALSE),
+  fInverseResponse(0x0),
+  fMeasuredEstimate(0x0),
+  fConditional(0x0),
+  fProjResponseInT(0x0),
+  fUnfolded(0x0),
+  fCoordinates2N(0x0),
+  fCoordinatesN_M(0x0),
+  fCoordinatesN_T(0x0)
+{
+  //
+  // default constructor
+  //
+}
+
+//______________________________________________________________
+
+AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar, 
+                              const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior) :
+  TNamed(name,title),
+  fResponse((THnSparse*)response->Clone()),
+  fPrior(0x0),
+  fOriginalPrior(0x0),
+  fEfficiency((THnSparse*)efficiency->Clone()),
+  fMeasured((THnSparse*)measured->Clone()),
+  fMaxNumIterations(0),
+  fNVariables(nVar),
+  fMaxChi2(0),
+  fUseSmoothing(kFALSE),
+  fInverseResponse(0x0),
+  fMeasuredEstimate(0x0),
+  fConditional(0x0),
+  fProjResponseInT(0x0),
+  fUnfolded(0x0),
+  fCoordinates2N(0x0),
+  fCoordinatesN_M(0x0),
+  fCoordinatesN_T(0x0)
+{
+  //
+  // named constructor
+  //
+
+  AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));
+  
+  if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
+  else {
+    fPrior = (THnSparse*) prior->Clone();
+    fOriginalPrior = (THnSparse*)fPrior->Clone();
+  }
+  
+  for (Int_t iVar=0; iVar<fNVariables; iVar++) {
+    AliInfo(Form("prior      matrix has %d bins in dimension %d",fPrior     ->GetAxis(iVar)->GetNbins(),iVar));
+    AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar));
+    AliInfo(Form("measured   matrix has %d bins in dimension %d",fMeasured  ->GetAxis(iVar)->GetNbins(),iVar));
+  }
+  Init();
+}
+
+
+//______________________________________________________________
+
+AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) :
+  TNamed(c),
+  fResponse((THnSparse*)c.fResponse->Clone()),
+  fPrior((THnSparse*)c.fPrior->Clone()),
+  fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()),
+  fEfficiency((THnSparse*)c.fEfficiency->Clone()),
+  fMeasured((THnSparse*)c.fMeasured->Clone()),
+  fMaxNumIterations(c.fMaxNumIterations),
+  fNVariables(c.fNVariables),
+  fMaxChi2(c.fMaxChi2),
+  fUseSmoothing(c.fUseSmoothing),
+  fInverseResponse((THnSparse*)c.fInverseResponse->Clone()),
+  fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()),
+  fConditional((THnSparse*)c.fConditional->Clone()),
+  fProjResponseInT((THnSparse*)c.fProjResponseInT->Clone()),
+  fUnfolded((THnSparse*)c.fUnfolded->Clone()),
+  fCoordinates2N(new Int_t(*c.fCoordinates2N)),
+  fCoordinatesN_M(new Int_t(*c.fCoordinatesN_M)),
+  fCoordinatesN_T(new Int_t(*c.fCoordinatesN_T))
+{
+  //
+  // copy constructor
+  //
+}
+
+//______________________________________________________________
+
+AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) {
+  //
+  // assignment operator
+  //
+  
+  if (this!=&c) {
+    TNamed::operator=(c);
+    fResponse = (THnSparse*)c.fResponse->Clone() ;
+    fPrior = (THnSparse*)c.fPrior->Clone() ;
+    fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ;
+    fEfficiency = (THnSparse*)c.fEfficiency->Clone() ;
+    fMeasured = (THnSparse*)c.fMeasured->Clone() ;
+    fMaxNumIterations = c.fMaxNumIterations ;
+    fNVariables = c.fNVariables ;
+    fMaxChi2 = c.fMaxChi2 ;
+    fUseSmoothing = c.fUseSmoothing ;
+    fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ;
+    fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ;
+    fConditional = (THnSparse*)c.fConditional->Clone() ;
+    fProjResponseInT = (THnSparse*)c.fProjResponseInT->Clone() ;
+    fUnfolded = (THnSparse*)c.fUnfolded->Clone() ;
+    fCoordinates2N  = new Int_t(*c.fCoordinates2N)  ;
+    fCoordinatesN_M = new Int_t(*c.fCoordinatesN_M) ;
+    fCoordinatesN_T = new Int_t(*c.fCoordinatesN_T) ;
+  }
+  return *this;
+}
+
+//______________________________________________________________
+
+AliCFUnfolding::~AliCFUnfolding() {
+  //
+  // destructor
+  //
+  if (fResponse)           delete fResponse;
+  if (fPrior)              delete fPrior;
+  if (fOriginalPrior)      delete fOriginalPrior;
+  if (fEfficiency)         delete fEfficiency;
+  if (fMeasured)           delete fMeasured;
+  if (fInverseResponse)    delete fInverseResponse;
+  if (fMeasuredEstimate)   delete fMeasuredEstimate;
+  if (fConditional)        delete fConditional;
+  if (fProjResponseInT)    delete fProjResponseInT;
+  if (fCoordinates2N)      delete [] fCoordinates2N; 
+  if (fCoordinatesN_M)     delete [] fCoordinatesN_M; 
+  if (fCoordinatesN_T)     delete [] fCoordinatesN_T; 
+}
+
+//______________________________________________________________
+
+void AliCFUnfolding::Init() {
+  //
+  // initialisation function : creates internal settings
+  //
+
+  fCoordinates2N  = new Int_t[2*fNVariables];
+  fCoordinatesN_M = new Int_t[fNVariables];
+  fCoordinatesN_T = new Int_t[fNVariables];
+
+  // create the matrix of conditional probabilities P(M|T)
+  CreateConditional();
+  
+  // create the frame of the inverse response matrix
+  fInverseResponse  = (THnSparse*) fResponse->Clone();
+  // create the frame of the unfolded spectrum
+  fUnfolded = (THnSparse*) fPrior->Clone();
+  // create the frame of the measurement estimate spectrum
+  fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
+}
+
+//______________________________________________________________
+
+void AliCFUnfolding::CreateEstMeasured() {
+  //
+  // This function creates a estimate (M) of the reconstructed spectrum 
+  // given the a priori distribution (T) and the conditional matrix (COND)
+  //
+  // --> P(M) = SUM   { P(M|T)    * P(T) }
+  // --> M(i) = SUM_k { COND(i,k) * T(k) }
+  //
+  // This is needed to calculate the inverse response matrix
+  //
+
+
+  // clean the measured estimate spectrum
+  for (Long64_t i=0; i<fMeasuredEstimate->GetNbins(); i++) {
+    fMeasuredEstimate->GetBinContent(i,fCoordinatesN_M);
+    fMeasuredEstimate->SetBinContent(fCoordinatesN_M,0.);
+  }
+  
+  // fill it
+  for (Int_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
+    Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
+    GetCoordinates();
+    Double_t priorValue = fPrior->GetBinContent(fCoordinatesN_T);
+    Double_t fill = fMeasuredEstimate->GetBinContent(fCoordinatesN_M) + conditionalValue * priorValue * fEfficiency->GetBinContent(fCoordinatesN_T);
+    if (fill>0.) fMeasuredEstimate->SetBinContent(fCoordinatesN_M,fill);
+  }
+}
+
+//______________________________________________________________
+
+void AliCFUnfolding::CreateInvResponse() {
+  //
+  // Creates the inverse response matrix (INV) with Bayesian method
+  //  : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E)
+  //
+  // --> P(T|M)   = P(M|T)    * P(T) * eff(T) / SUM   { P(M|T)    * P(T) }
+  // --> INV(i,j) = COND(i,j) * T(j) * E(j)   / SUM_k { COND(i,k) * T(k) }
+  //
+
+  for (Int_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
+    Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
+    GetCoordinates();
+    Double_t priorValue = fPrior->GetBinContent(fCoordinatesN_T);
+    Double_t estimatedMeasured = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
+    Double_t fill = (estimatedMeasured>0. ? conditionalValue * priorValue * fEfficiency->GetBinContent(fCoordinatesN_T) / estimatedMeasured : 0. ) ;
+    if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) fInverseResponse->SetBinContent(fCoordinates2N,fill);
+  }
+}
+
+//______________________________________________________________
+
+void AliCFUnfolding::Unfold() {
+  //
+  // Main routine called by the user : 
+  // it calculates the unfolded spectrum from the response matrix and the measured spectrum
+  // several iterations are performed until a reasonable chi2 is reached
+  //
+
+  Int_t iIterBayes=0 ;
+  Double_t chi2=0 ;
+
+  for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations
+    CreateEstMeasured();
+    CreateInvResponse();
+    CreateUnfolded();
+    chi2 = GetChi2();
+    //printf("chi2 = %e\n",chi2);
+    if (fMaxChi2>0. && chi2<fMaxChi2) {
+      break;
+    }
+    // update the prior distribution
+    if (fUseSmoothing) Smooth();
+    fPrior = (THnSparse*)fUnfolded->Clone() ; // this should be changed (memory)
+  }
+  AliInfo(Form("Finished at iteration %d : Chi2 is %e and you required it to be < %e",iIterBayes,chi2,fMaxChi2));
+}
+
+//______________________________________________________________
+
+void AliCFUnfolding::CreateUnfolded() {
+  //
+  // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV)
+  // We have P(T) = SUM   { P(T|M)   * P(M) } 
+  //   -->   T(i) = SUM_k { INV(i,k) * M(k) }
+  //
+
+
+  // clear the unfolded spectrum
+  for (Long64_t i=0; i<fUnfolded->GetNbins(); i++) {
+    fUnfolded->GetBinContent(i,fCoordinatesN_T);
+    fUnfolded->SetBinContent(fCoordinatesN_T,0.);
+  }
+  
+  for (Int_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
+    Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
+    GetCoordinates();
+    Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T);
+    Double_t fill = fUnfolded->GetBinContent(fCoordinatesN_T) + (effValue>0. ? invResponseValue*fMeasured->GetBinContent(fCoordinatesN_M)/effValue : 0.) ;
+    if (fill>0.) fUnfolded->SetBinContent(fCoordinatesN_T,fill);
+  }
+}
+    
+void AliCFUnfolding::GetCoordinates() {
+  //
+  // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
+  //
+  for (Int_t i = 0; i<fNVariables ; i++) {
+    fCoordinatesN_M[i] = fCoordinates2N[i];
+    fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
+  }
+}
+
+//______________________________________________________________
+
+void AliCFUnfolding::CreateConditional() {
+  //
+  // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
+  //
+  //  --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
+  //
+
+  fConditional     = (THnSparse*) fResponse->Clone(); // output of this function
+  fProjResponseInT = (THnSparse*) fPrior->Clone();    // output denominator : 
+                                                      // projection of the response matrix on the TRUE axis
+  
+  // set in fProjResponseInT zero everywhere
+  for (Int_t iBin=0; iBin<fProjResponseInT->GetNbins(); iBin++) {
+    fProjResponseInT->GetBinContent(iBin,fCoordinatesN_T);
+    fProjResponseInT->SetBinContent(fCoordinatesN_T,0.);
+  }
+
+  // calculate the response projection on T axis
+  for (Int_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
+    Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
+    GetCoordinates();
+    Double_t fill = fProjResponseInT->GetBinContent(fCoordinatesN_T) + responseValue ;
+    if (fill>0.) fProjResponseInT->SetBinContent(fCoordinatesN_T,fill);
+  }
+  
+  // fill the conditional probability matrix
+  for (Int_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
+    Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
+    GetCoordinates();
+    Double_t fill = responseValue / fProjResponseInT->GetBinContent(fCoordinatesN_T) ;
+    if (fill>0. || fConditional->GetBinContent(fCoordinates2N)) fConditional->SetBinContent(fCoordinates2N,fill);
+  }
+}
+
+//______________________________________________________________
+
+Double_t AliCFUnfolding::GetChi2() {
+  //
+  // Returns the chi2 between unfolded and a priori spectrum
+  //
+
+  Double_t chi2 = 0. ;
+  for (Int_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
+    Double_t priorValue = fPrior->GetBinContent(iBin);
+    chi2 += (priorValue>0. ? TMath::Power(fUnfolded->GetBinContent(iBin) - priorValue,2) / priorValue : 0.) ;
+  }
+  return chi2;
+}
+
+//______________________________________________________________
+
+void AliCFUnfolding::SetMaxChi2PerDOF(Double_t val) {
+  //
+  // Max. chi2 per degree of freedom : user setting
+  //
+
+  Int_t nDOF = 1 ;
+  for (Int_t iDim=0; iDim<fNVariables; iDim++) {
+    nDOF *= fPrior->GetAxis(iDim)->GetNbins();
+  }
+  AliInfo(Form("Number of degrees of freedom = %d",nDOF));
+  fMaxChi2 = val * nDOF ;
+}
+
+//______________________________________________________________
+
+void AliCFUnfolding::Smooth() {
+  //
+  // Smoothes the unfolded spectrum
+  // Each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
+  //
+  
+  Int_t* numBins = new Int_t[fNVariables];
+  for (Int_t iVar=0; iVar<fNVariables; iVar++) numBins[iVar]=fUnfolded->GetAxis(iVar)->GetNbins();
+  
+  //need a copy because fUnfolded will be updated during the loop, and this creates problems
+  THnSparse* copy = (THnSparse*)fUnfolded->Clone();
+
+  for (Int_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
+    Double_t content = copy->GetBinContent(iBin,fCoordinatesN_T);
+
+    // skip the under/overflow bins...
+    Bool_t isOutside = kFALSE ;
+    for (Int_t iVar=0; iVar<fNVariables; iVar++) {
+      if (fCoordinatesN_T[iVar]<1 || fCoordinatesN_T[iVar]>numBins[iVar]) {
+       isOutside=kTRUE;
+       break;
+      }
+    }
+    if (isOutside) continue;
+    
+    Int_t neighbours = 0; // number of neighbours to average with
+
+    for (Int_t iVar=0; iVar<fNVariables; iVar++) {
+      if (fCoordinatesN_T[iVar] > 1) { // must not be on low edge border
+       fCoordinatesN_T[iVar]-- ; //get lower neighbouring bin 
+       Double_t contentNeighbour = copy->GetBinContent(fCoordinatesN_T);
+       content += contentNeighbour;
+       neighbours++;
+       fCoordinatesN_T[iVar]++ ; //back to initial coordinate
+      }
+      if (fCoordinatesN_T[iVar] < numBins[iVar]) { // must not be on up edge border
+       fCoordinatesN_T[iVar]++ ; //get upper neighbouring bin
+       Double_t contentNeighbour = copy->GetBinContent(fCoordinatesN_T);
+       content += contentNeighbour ;
+       neighbours++;
+       fCoordinatesN_T[iVar]-- ; //back to initial coordinate
+      }
+    }
+    content /= (1+neighbours) ; // make an average
+    fUnfolded->SetBinContent(fCoordinatesN_T,content);
+  }
+  delete [] numBins;
+  delete copy;
+}
+
+
+//______________________________________________________________
+
+void AliCFUnfolding::CreateFlatPrior() {
+  //
+  // Creates a flat prior distribution
+  // 
+
+  AliInfo("Creating a flat a priori distribution");
+  
+  // create the frame of the THnSparse given (for example) the one from the efficiency map
+  fPrior = (THnSparse*) fEfficiency->Clone();
+
+  Int_t nDim = fNVariables;
+  Int_t* bins = new Int_t[nDim]; // number of bins for each variable
+  Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
+
+  for (Int_t iVar=0; iVar<nDim; iVar++) {
+    bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
+    nBins *= bins[iVar];
+  }
+
+  Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
+
+  // loop that sets 1 in each bin
+  for (Long_t iBin=0; iBin<nBins; iBin++) {
+    Long_t bin_tmp = iBin ;
+    for (Int_t iVar=0; iVar<nDim; iVar++) {
+      bin[iVar] = 1 + bin_tmp % bins[iVar] ;
+      bin_tmp /= bins[iVar] ;
+    }
+    fPrior->SetBinContent(bin,1.); // put 1 everywhere
+  }
+  
+  fOriginalPrior = (THnSparse*)fPrior->Clone();
+
+  delete [] bin;
+  delete [] bins;
+}
diff --git a/CORRFW/AliCFUnfolding.h b/CORRFW/AliCFUnfolding.h
new file mode 100644 (file)
index 0000000..e7e211a
--- /dev/null
@@ -0,0 +1,84 @@
+
+#ifndef ALICFUNFOLDING_H
+#define ALICFUNFOLDING_H
+
+//--------------------------------------------------------------------//
+//                                                                    //
+// AliCFUnfolding Class                                               //
+// Class to handle general unfolding procedure                        // 
+// For the moment only bayesian unfolding is supported                //
+// The next steps are to add chi2 minimisation and weighting methods  //
+//                                                                    //
+// Author : renaud.vernet@cern.ch                                     //
+//--------------------------------------------------------------------//
+
+#include "TNamed.h"
+#include "THnSparse.h"
+
+class AliCFUnfolding : public TNamed {
+
+ public :
+  AliCFUnfolding();
+  AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar, 
+                const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior=0x0);
+  AliCFUnfolding(const AliCFUnfolding& c);
+  AliCFUnfolding& operator= (const AliCFUnfolding& c);
+  ~AliCFUnfolding();
+
+  void SetMaxNumberOfIterations(Int_t n)  {fMaxNumIterations=n;}
+  void SetMaxChi2(Double_t val)           {fMaxChi2=val;}
+  void SetMaxChi2PerDOF(Double_t val);
+  void UseSmoothing(Bool_t b=kTRUE)       {fUseSmoothing=b;}
+  void Unfold();
+
+  THnSparse* GetResponse()        const {return fResponse;}
+  THnSparse* GetInverseResponse() const {return fInverseResponse;}
+  THnSparse* GetPrior()           const {return fPrior;}
+  THnSparse* GetOriginalPrior()   const {return fOriginalPrior;}
+  THnSparse* GetEfficiency()      const {return fEfficiency;}
+  THnSparse* GetUnfolded()        const {return fUnfolded;}
+
+ private :
+  
+  // user-related settings
+  THnSparse     *fResponse;           // Response matrix : dimensions must be 2N = 2 x (number of variables)
+                                      // first N dimensions must be filled with reconstructed values
+                                      // last  N dimensions must be filled with generated values
+  THnSparse     *fPrior;              // This is the assumed generated distribution : dimensions must be N = number of variables
+                                      // it will be used at the first step 
+                                      // then will be updated automatically at each iteration
+  THnSparse     *fOriginalPrior;      // This is the original prior distribution : will not be modified
+  THnSparse     *fEfficiency;         // Efficiency map : dimensions must be N = number of variables
+                                      // this map must be filled only with "true" values of the variables (should not include resolution effects)
+  THnSparse     *fMeasured;           // Measured spectrum to be unfolded : dimensions must be N = number of variables
+  Int_t          fMaxNumIterations;   // Maximum  number of iterations to be performed
+  Int_t          fNVariables;         // Number of variables used in analysis spectra (pt, y, ...)
+  Double_t       fMaxChi2;            // Maximum Chi2 between unfolded and prior distributions. 
+  Bool_t         fUseSmoothing;       // Smooth the unfolded sectrum at each iteration
+
+  // internal settings
+  THnSparse     *fInverseResponse;   // Inverse response matrix
+  THnSparse     *fMeasuredEstimate;  // Estimation of the measured (M) spectrum given the a priori (T) distribution
+  THnSparse     *fConditional;       // Matrix holding the conditional probabilities P(M|T)
+  THnSparse     *fProjResponseInT;   // Projection of the response matrix on TRUE axis
+  THnSparse     *fUnfolded;          // Unfolded spectrum
+  Int_t         *fCoordinates2N;     // Coordinates in 2N (measured,true) space
+  Int_t         *fCoordinatesN_M;    // Coordinates in measured space
+  Int_t         *fCoordinatesN_T;    // Coordinates in true space
+  
+
+  // functions
+  void     Init();                // initialisation of the internal settings
+  void     GetCoordinates();      // gets a cell coordinates in Measured and True space
+  void     CreateConditional();   // creates the conditional matrix from the response matrix
+  void     CreateEstMeasured();   // creates the measured spectrum estimation from the conditional matrix and the prior distribution
+  void     CreateInvResponse();   // creates the inverse response function (Bayes Theorem) from the conditional matrix and the prior distribution
+  void     CreateUnfolded();      // creates the unfolded spectrum from the inverse response matrix and the measured distribution
+  void     CreateFlatPrior();     // creates a flat a priori distribution in case the one given in the constructor is null
+  Double_t GetChi2();             // returns the chi2 between unfolded and prior spectra
+  void     Smooth();              // smooth the unfolded spectrum using the neighbouring cells
+
+  ClassDef(AliCFUnfolding,0);
+};
+
+#endif
index df74690..75c4ee4 100644 (file)
@@ -29,5 +29,6 @@
 #pragma link C++ class  AliCFPairIsPrimaryCuts+;
 #pragma link C++ class  AliCFPairPidCut+;
 #pragma link C++ class  AliCFV0TopoCuts+;
+#pragma link C++ class  AliCFUnfolding+;
 
 #endif
index c650d31..3535b14 100644 (file)
@@ -23,7 +23,8 @@ SRCS = AliCFFrame.cxx \
        AliCFPairQualityCuts.cxx \
        AliCFPairIsPrimaryCuts.cxx \
        AliCFPairPidCut.cxx \
-       AliCFV0TopoCuts.cxx  
+       AliCFV0TopoCuts.cxx \
+       AliCFUnfolding.cxx
 
 CHECKALIEN = $(shell root-config --has-alien)
 ifeq (yes,$(CHECKALIEN))
diff --git a/CORRFW/test/AliCFTaskForUnfolding.C b/CORRFW/test/AliCFTaskForUnfolding.C
new file mode 100644 (file)
index 0000000..446097b
--- /dev/null
@@ -0,0 +1,187 @@
+//DEFINITION OF A FEW CONSTANTS
+const Double_t ptmin  =   0.0 ;
+const Double_t ptmax  =   1.0 ;
+const Double_t etamin =  -1.5 ;
+const Double_t etamax =   1.5 ;
+const Int_t    PDG    =    211; 
+const Int_t    minclustersTPC = 40 ;
+//----------------------------------------------------
+
+Bool_t AliCFTaskForUnfolding()
+{
+  
+  TBenchmark benchmark;
+  benchmark.Start("AliSingleTrackTask");
+
+  AliLog::SetGlobalDebugLevel(0);
+
+  Load() ; //load the required libraries
+
+  TChain * analysisChain ;
+  printf("\n\nRunning on local file, please check the path\n\n");
+
+  analysisChain = new TChain("esdTree");
+  analysisChain->Add("AliESDs.root");
+  
+  Info("AliCFTaskForUnfolding",Form("CHAIN HAS %d ENTRIES",(Int_t)analysisChain->GetEntries()));
+
+  //CONTAINER DEFINITION
+  Info("AliCFTaskForUnfolding","SETUP CONTAINER");
+  //the sensitive variables (2 in this example), their indices
+  UInt_t ipt = 0;
+  UInt_t ieta  = 1;
+  //Setting up the container grid... 
+  UInt_t nstep = 3 ; //number of selection steps MC 
+  const Int_t nvar   = 2 ; //number of variables on the grid:pt,eta
+  const Int_t nbin[nvar] = {20,20} ;
+
+  //arrays for the number of bins in each dimension
+  Int_t iBin[nvar];
+  iBin[0]=nbin[0];
+  iBin[1]=nbin[1];
+
+  //arrays for lower bounds :
+  Double_t *binLim0=new Double_t[nbin[0]+1];
+  Double_t *binLim1=new Double_t[nbin[1]+1];
+
+  //values for bin lower bounds
+  for(Int_t i=0; i<=nbin[0]; i++) binLim0[i]=(Double_t)ptmin  + ( ptmax- ptmin)/nbin[0]*(Double_t)i ; 
+  for(Int_t i=0; i<=nbin[1]; i++) binLim1[i]=(Double_t)etamin + (etamax-etamin)/nbin[1]*(Double_t)i ; 
+
+  //one "container" for MC
+  AliCFContainer* container = new AliCFContainer("container","container for tracks",nstep,nvar,iBin);
+  //setting the bin limits
+  container -> SetBinLimits(ipt,binLim0);
+  container -> SetBinLimits(ieta ,binLim1);
+
+  // Gen-Level kinematic cuts
+  AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts");
+  mcKineCuts->SetPtRange (ptmin ,ptmax );
+  mcKineCuts->SetEtaRange(etamin,etamax);
+
+  //Particle-Level cuts:  
+  AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts");
+  mcGenCuts->SetRequireIsPrimary();
+  mcGenCuts->SetRequirePdgCode(PDG);
+
+  // Rec-Level kinematic cuts
+  AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts");
+  recKineCuts->SetPtRange( ptmin, ptmax);
+  recKineCuts->SetPtRange(etamin,etamax);
+
+  AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
+  recQualityCuts->SetMinNClusterTPC(minclustersTPC);
+
+  AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
+  recIsPrimaryCuts->SetMaxNSigmaToVertex(3);
+
+  printf("CREATE MC KINE CUTS\n");
+  TObjArray* mcList = new TObjArray(0) ;
+  mcList->AddLast(mcKineCuts);
+  mcList->AddLast(mcGenCuts);
+
+  printf("CREATE RECONSTRUCTION CUTS\n");
+  TObjArray* recList = new TObjArray(0) ;
+  recList->AddLast(recKineCuts);
+  recList->AddLast(recQualityCuts);
+  recList->AddLast(recIsPrimaryCuts);
+
+  TObjArray* emptyList = new TObjArray(0);
+
+  //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
+  printf("CREATE INTERFACE AND CUTS\n");
+  AliCFManager* man = new AliCFManager() ;
+  man->SetNStepEvent(0);
+  man->SetParticleContainer(container);
+  man->SetParticleCutsList(0,mcList);
+  man->SetParticleCutsList(1,recList);
+  man->SetParticleCutsList(2,emptyList); // this is special step for monte carlo spectrum
+
+  //CREATE THE TASK
+  printf("CREATE TASK\n");
+  // create the task
+  AliCFTaskForUnfolding *task = new AliCFTaskForUnfolding("AliCFTaskForUnfolding");
+  task->SetCFManager(man); //here is set the CF manager
+
+  //create correlation matrix for unfolding
+  Int_t thnDim[2*nvar];
+  for (int k=0; k<nvar; k++) {
+    //first half  : reconstructed 
+    //second half : MC
+    thnDim[k]      = nbin[k];
+    thnDim[k+nvar] = nbin[k];
+  }
+  THnSparseD* correlation = new THnSparseD("correlation","THnSparse with correlations",2*nvar,thnDim);
+  Double_t** binEdges = new Double_t[nvar];
+  for (int k=0; k<nvar; k++) {
+    binEdges[k]=new Double_t[nbin[k]+1];
+    container->GetBinLimits(k,binEdges[k]);
+    correlation->SetBinEdges(k,binEdges[k]);
+    correlation->SetBinEdges(k+nvar,binEdges[k]);
+  }
+  task->SetCorrelationMatrix(correlation);
+  //---
+
+  //SETUP THE ANALYSIS MANAGER TO READ INPUT CHAIN AND WRITE DESIRED OUTPUTS
+  printf("CREATE ANALYSIS MANAGER\n");
+  // Make the analysis manager
+  AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
+  mgr->SetAnalysisType(AliAnalysisManager::kLocalAnalysis);
+
+  AliMCEventHandler*  mcHandler = new AliMCEventHandler();
+  mgr->SetMCtruthEventHandler(mcHandler);
+  AliInputEventHandler* dataHandler = new AliESDInputHandler();
+  mgr->SetInputEventHandler(dataHandler);
+
+  // Create and connect containers for input/output
+
+  //------ input data ------
+  AliAnalysisDataContainer *cinput0  = mgr->CreateContainer("cchain0",TChain::Class(),AliAnalysisManager::kInputContainer);
+
+  // ----- output data -----
+  
+  //slot 0 : default output tree (by default handled by AliAnalysisTaskSE)
+  AliAnalysisDataContainer *coutput0 = mgr->CreateContainer("ctree0", TTree::Class(),AliAnalysisManager::kOutputContainer,"output.root");
+
+  //now comes user's output objects :
+  
+  // output TH1I for event counting
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist0", TH1I::Class(),AliAnalysisManager::kOutputContainer,"output.root");
+  // output Correction Framework Container (for acceptance & efficiency calculations)
+  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("ccontainer0", AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,"output.root");
+  AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("corr0", THnSparseD::Class(),AliAnalysisManager::kOutputContainer,"output.root");
+
+  cinput0->SetData(analysisChain);
+
+  mgr->AddTask(task);
+  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(task,0,coutput0);
+  mgr->ConnectOutput(task,1,coutput1);
+  mgr->ConnectOutput(task,2,coutput2);
+  mgr->ConnectOutput(task,3,coutput3);
+  printf("READY TO RUN\n");
+  //RUN !!!
+  if (mgr->InitAnalysis()) {
+    mgr->PrintStatus();
+    mgr->StartAnalysis("local",analysisChain);
+  }
+
+  benchmark.Stop("AliSingleTrackTask");
+  benchmark.Show("AliSingleTrackTask");
+
+  return kTRUE ;
+}
+
+void Load() {
+
+  //load the required aliroot libraries
+  gSystem->Load("libANALYSIS") ;
+  gSystem->Load("libANALYSISalice") ;
+  gSystem->Load("libCORRFW.so") ;
+
+  //compile online the task class
+  gSystem->SetIncludePath("-I. -I$ALICE_ROOT/include -I$ROOTSYS/include");
+  gROOT->LoadMacro("./AliCFTaskForUnfolding.cxx+");
+}
diff --git a/CORRFW/test/AliCFTaskForUnfolding.cxx b/CORRFW/test/AliCFTaskForUnfolding.cxx
new file mode 100644 (file)
index 0000000..44b71e5
--- /dev/null
@@ -0,0 +1,213 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------------
+// Task to prepare efficiency and response matrix for unfolding
+//-----------------------------------------------------------------------
+// Author : R. Vernet, INFN - Catania (it)
+//-----------------------------------------------------------------------
+
+
+#include "TStyle.h"
+#include "AliCFTaskForUnfolding.h"
+#include "TCanvas.h"
+#include "AliStack.h"
+#include "TParticle.h"
+#include "TH1I.h"
+#include "AliMCEvent.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliCFManager.h"
+#include "AliCFCutBase.h"
+#include "AliCFContainer.h"
+#include "TChain.h"
+#include "AliESDtrack.h"
+#include "AliLog.h"
+#include "THnSparse.h"
+#include "TH2D.h"
+
+ClassImp(AliCFTaskForUnfolding)
+
+//__________________________________________________________________________
+AliCFTaskForUnfolding::AliCFTaskForUnfolding() :
+  fCFManager(0x0),
+  fHistEventsProcessed(0x0),
+  fCorrelation(0x0)
+{
+  //
+  //Default ctor
+  //
+}
+//___________________________________________________________________________
+AliCFTaskForUnfolding::AliCFTaskForUnfolding(const Char_t* name) :
+  AliAnalysisTaskSE(name),
+  fCFManager(0x0),
+  fHistEventsProcessed(0x0),
+  fCorrelation(0x0)
+{
+  //
+  // Constructor. Initialization of Inputs and Outputs
+  //
+  Info("AliCFTaskForUnfolding","Calling Constructor");
+
+  /*
+    DefineInput(0) and DefineOutput(0)
+    are taken care of by AliAnalysisTaskSE constructor
+  */
+  DefineOutput(1,TH1I::Class());
+  DefineOutput(2,AliCFContainer::Class());
+  DefineOutput(3,THnSparseD::Class());
+}
+
+//___________________________________________________________________________
+AliCFTaskForUnfolding::~AliCFTaskForUnfolding() {
+  //
+  //destructor
+  //
+  Info("~AliCFTaskForUnfolding","Calling Destructor");
+  if (fCFManager)           delete fCFManager ;
+  if (fHistEventsProcessed) delete fHistEventsProcessed ;
+  if (fCorrelation)         delete fCorrelation ;
+}
+
+//_________________________________________________
+void AliCFTaskForUnfolding::UserExec(Option_t *)
+{
+  //
+  // Main loop function
+  //
+  Info("UserExec","") ;
+
+  AliVEvent*    fEvent = fInputEvent ;
+  AliVParticle* track ;
+  
+  if (!fEvent) {
+    Error("UserExec","NO EVENT FOUND!");
+    return;
+  }
+
+  if (!fMCEvent) Error("UserExec","NO MC INFO FOUND");
+  
+  //pass the MC evt handler to the cuts that need it 
+  fCFManager->SetEventInfo(fMCEvent);
+
+  // MC-event selection
+  Double_t containerInput[2] ;
+        
+  //loop on the MC event
+  for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) { 
+    AliMCParticle *mcPart  = fMCEvent->GetTrack(ipart);
+    
+    if (!fCFManager->CheckParticleCuts(0,mcPart)) continue;
+    containerInput[0] = (Float_t)mcPart->Pt();
+    containerInput[1] = (Float_t)mcPart->Eta();
+    fCFManager->GetParticleContainer()->Fill(containerInput,0);
+  }    
+
+  //Now go to rec level
+  for (Int_t iTrack = 0; iTrack<fEvent->GetNumberOfTracks(); iTrack++) {
+    
+    track = fEvent->GetTrack(iTrack);
+    if (!fCFManager->CheckParticleCuts(1,track)) continue;
+    
+    Int_t label = track->GetLabel();
+    if (label<0) continue;
+    AliMCParticle* mcPart = fMCEvent->GetTrack(label);
+    // check if this track was part of the signal
+    if (!fCFManager->CheckParticleCuts(0,mcPart)) continue;
+
+    //fill the container
+    Double_t mom[3];
+    track->PxPyPz(mom);
+    Double_t pt=TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
+    containerInput[0] = pt;
+    containerInput[1] = track->Eta();
+    fCFManager->GetParticleContainer()->Fill(containerInput,1) ;
+    containerInput[0] = mcPart->Pt();
+    containerInput[1] = mcPart->Eta();
+    fCFManager->GetParticleContainer()->Fill(containerInput,2);
+    
+    Double_t fill[4]; //fill response matrix
+    // dimensions 0&1 : pt,eta (Rec)
+    fill[0] = pt ;
+    fill[1] = track->Eta();
+    // dimensions 2&3 : pt,eta (MC)
+    fill[2] = mcPart->Pt();
+    fill[3] = mcPart->Eta();
+    fCorrelation->Fill(fill);
+  }
+  
+  fHistEventsProcessed->Fill(0);
+
+  /* PostData(0) is taken care of by AliAnalysisTaskSE */
+  PostData(1,fHistEventsProcessed) ;
+  PostData(2,fCFManager->GetParticleContainer()) ;
+  PostData(3,fCorrelation) ;
+}
+
+
+//___________________________________________________________________________
+void AliCFTaskForUnfolding::Terminate(Option_t*)
+{
+  // The Terminate() function is the last function to be called during
+  // a query. It always runs on the client, it can be used to present
+  // the results graphically or save the results to file.
+
+  Info("Terminate","");
+  AliAnalysisTaskSE::Terminate();
+
+  gStyle->SetPalette(1);
+
+  //draw some example plots....
+
+  AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
+  TH2D* h00 =   cont->ShowProjection(0,1,0) ;
+  TH2D* h01 =   cont->ShowProjection(0,1,1) ;
+  THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
+
+  TCanvas * c =new TCanvas("c","",800,400);
+  c->Divide(2,1);
+  c->cd(1);
+  h00->Draw("text");
+  c->cd(2);
+  h01->Draw("text");
+  c->SaveAs("spectra.eps");
+
+  TCanvas * c2 =new TCanvas("c2","",800,400);
+  c2->Divide(2,1);
+  c2->cd(1);
+  hcorr->Projection(0,2)->Draw("text");
+  c2->cd(2);
+  hcorr->Projection(1,3)->Draw("text");
+  c2->SaveAs("correlation.eps");
+}
+
+
+//___________________________________________________________________________
+void AliCFTaskForUnfolding::UserCreateOutputObjects() {
+  //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
+  //TO BE SET BEFORE THE EXECUTION OF THE TASK
+  //
+  Info("CreateOutputObjects","CreateOutputObjects of task %s", GetName());
+
+  //slot #1
+  OpenFile(1);
+  fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
+
+//   OpenFile(2);
+//   OpenFile(3);
+}
+
diff --git a/CORRFW/test/AliCFTaskForUnfolding.h b/CORRFW/test/AliCFTaskForUnfolding.h
new file mode 100644 (file)
index 0000000..8507d87
--- /dev/null
@@ -0,0 +1,62 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------------
+// Author : R. Vernet, INFN - Catania (it)
+//-----------------------------------------------------------------------
+
+#ifndef ALICFTASKFORUNFOLDING_H
+#define ALICFTASKFORUNFOLDING_H
+
+#include "AliAnalysisTaskSE.h"
+
+class TH1I;
+class TParticle ;
+class TFile ;
+class AliStack ;
+class AliCFManager;
+class AliESDtrack;
+class AliVParticle;
+class THnSparse;
+
+class AliCFTaskForUnfolding : public AliAnalysisTaskSE {
+  public:
+
+  AliCFTaskForUnfolding();
+  AliCFTaskForUnfolding(const Char_t* name);
+  virtual ~AliCFTaskForUnfolding();
+
+  // ANALYSIS FRAMEWORK STUFF to loop on data and fill output objects
+  void     UserCreateOutputObjects();
+  void     UserExec(Option_t *option);
+  void     Terminate(Option_t *);
+  
+  // CORRECTION FRAMEWORK RELATED FUNCTIONS
+  void           SetCFManager(AliCFManager* io) {fCFManager = io;}   // global correction manager
+  AliCFManager * GetCFManager() const {return fCFManager;}           // get corr manager
+  void           SetCorrelationMatrix(THnSparse* h) {fCorrelation=h;}
+
+ protected:
+  
+  AliCFManager   *fCFManager    ;  // pointer to the CF manager
+
+  // Histograms
+  //Number of events
+  TH1I  *fHistEventsProcessed; // simple histo for monitoring the number of events processed
+  THnSparse* fCorrelation;          //response matrix for unfolding  
+  ClassDef(AliCFTaskForUnfolding,0);
+};
+
+#endif
diff --git a/CORRFW/test/testUnfolding.C b/CORRFW/test/testUnfolding.C
new file mode 100644 (file)
index 0000000..36880dc
--- /dev/null
@@ -0,0 +1,119 @@
+
+void testUnfolding() {
+  Load();
+
+  // get the essential
+  AliCFDataGrid *measured    = (AliCFDataGrid*) GetMeasuredSpectrum();
+  AliCFDataGrid *generated   = (AliCFDataGrid*) GetGeneratedSpectrum();
+  AliCFEffGrid  *efficiency  = (AliCFEffGrid*)  GetEfficiency();
+  THnSparse     *response    = (THnSparse*)     GetResponseMatrix();
+  
+  // create a guessed "a priori" distribution using binning of MC
+  THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ;
+  //----
+  AliCFUnfolding unfolding("unfolding","",2,response,efficiency->GetGrid(),((AliCFGridSparse*)measured->GetData())->GetGrid(),guessed);
+  unfolding.SetMaxNumberOfIterations(100);
+  unfolding.SetMaxChi2PerDOF(1.e-07);
+  //unfolding.UseSmoothing();
+  unfolding.Unfold();
+
+  THnSparse* result = unfolding.GetUnfolded();
+  //----
+  
+  TCanvas * canvas = new TCanvas("canvas","",1000,700);
+  canvas->Divide(3,3);
+
+  canvas->cd(1);
+  TH2D* h_gen = generated->Project(0,1);
+  h_gen->SetTitle("generated");
+  h_gen->Draw("lego2");
+
+  canvas->cd(2);
+  TH2D* h_meas = measured->Project(0,1);
+  h_meas->SetTitle("measured");
+  h_meas->Draw("lego2");
+  
+  canvas->cd(3);
+  TH2D* h_guessed = guessed->Projection(1,0);
+  h_guessed->SetTitle("a priori");
+  h_guessed->Draw("lego2");
+
+  canvas->cd(4);
+  TH2D* h_eff = efficiency->Project(0,1);
+  h_eff->SetTitle("efficiency");
+  h_eff->Draw("lego2");
+
+  canvas->cd(5);
+  TH2D* h_unf = result->Projection(1,0);
+  h_unf->SetTitle("unfolded");
+  h_unf->Draw("lego2");
+
+  canvas->cd(6);
+  TH2D* ratio = (TH2D*)h_unf->Clone();
+  ratio->SetTitle("unfolded / generated");
+  ratio->Divide(h_unf,h_gen,1,1);
+//   ratio->Draw("cont4z");
+//   ratio->Draw("surf2");
+  ratio->Draw("lego2");
+
+  canvas->cd(7);
+  TH2* orig = unfolding.GetOriginalPrior()->Projection(1,0);
+  orig->SetTitle("original prior");
+  orig->Draw("lego2");
+
+  canvas->cd(8);
+  AliCFDataGrid* corrected = (AliCFDataGrid*)measured->Clone();
+  corrected->ApplyEffCorrection(*efficiency);
+  TH2D* corr = corrected->Project(0,1);
+  corr->SetTitle("simple correction");
+  corr->Draw("lego2");
+
+  canvas->cd(9);
+  TH2D* ratio2 = (TH2D*) corr->Clone();
+  ratio2->Divide(corr,h_gen,1,1);
+  ratio2->SetTitle("simple correction / generated");
+  ratio2->Draw("cont4z");
+
+  return;
+}
+
+// ====================================================================
+
+
+void Load(){
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libCORRFW");
+  gStyle->SetPalette(1);
+  gStyle->SetOptStat(0);
+}
+
+TObject* GetMeasuredSpectrum() {
+  TFile * f = TFile::Open("test/output.root","read");
+  AliCFContainer* c = (AliCFContainer*)f->Get("container");
+  AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
+  data->SetMeasured(1);
+  return data;
+}
+TObject* GetGeneratedSpectrum() {
+  TFile * f = TFile::Open("test/output.root","read");
+  AliCFContainer* c = (AliCFContainer*)f->Get("container");
+  AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
+  data->SetMeasured(0);
+  return data;
+}
+TObject* GetEfficiency() {
+  TFile * f = TFile::Open("test/output.root","read");
+  AliCFContainer* c = (AliCFContainer*)f->Get("container");
+  AliCFEffGrid* eff = new AliCFEffGrid("eff","",*c);
+  eff->CalculateEfficiency(2,0);
+  return eff;
+}
+THnSparse* GetResponseMatrix() {
+  TFile * f = TFile::Open("test/output.root","read");
+  THnSparse* h = (THnSparse*)f->Get("correlation");
+  return h;
+}
+THnSparse* CreateGuessed(const THnSparse* h) {
+  THnSparse* guessed = (THnSparse*) h->Clone();
+  return guessed ;
+}