bug fix - Adam
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 May 2012 16:16:09 +0000 (16:16 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 May 2012 16:16:09 +0000 (16:16 +0000)
EMCAL/AliEMCALUnfolding.cxx

index 37f878e..8abbc4c 100644 (file)
-/**************************************************************************
- * 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.                  *
- **************************************************************************/
-
-//_________________________________________________________________________
-//  Base class for the cluster unfolding algorithm 
-//*-- Author: Adam Matyja (SUBATECH)
-//  Based on unfolding in clusterizerv1 done by Cynthia Hadjidakis
-//-- Unfolding for eta~0: Cynthia Hadjidakis - still in AliEMCALCLusterizerv1
-//-- Unfolding extension for whole EMCAL: Adam Matyja (SUBATECH & INP PAN)
-//
-//  unfolds the clusters having several local maxima. 
-//////////////////////////////////////////////////////////////////////////////
-
-// --- ROOT system ---
-#include "TClonesArray.h"
-#include <TMath.h> 
-#include <TMinuit.h>
-
-// --- Standard library ---
-#include <cassert>
-
-// --- AliRoot header files ---
-#include "AliEMCALUnfolding.h"
-#include "AliEMCALGeometry.h"
-#include "AliRunLoader.h"
-#include "AliRun.h"
-#include "AliEMCAL.h"
-#include "AliEMCALRecParam.h"
-#include "AliEMCALRecPoint.h"
-#include "AliEMCALDigit.h"
-#include "AliEMCALReconstructor.h"
-
-#include "AliLog.h"
-#include "AliCDBManager.h"
-class AliCDBStorage;
-#include "AliCDBEntry.h"
-
-Double_t AliEMCALUnfolding::fgSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};
-Double_t AliEMCALUnfolding::fgPar5[3]={12.31,-0.007381,-0.06936};
-Double_t AliEMCALUnfolding::fgPar6[3]={0.05452,0.0001228,0.001361};
-
-ClassImp(AliEMCALUnfolding)
-  
-//____________________________________________________________________________
-AliEMCALUnfolding::AliEMCALUnfolding():
-  fNumberOfECAClusters(0),
-  fECALocMaxCut(0),
-  fThreshold(0.01),//10 MeV
-  fGeom(NULL),
-  fRecPoints(NULL),
-  fDigitsArr(NULL)
-{
-  // ctor with the indication of the file where header Tree and digits Tree are stored
-  Init() ;
-}
-
-//____________________________________________________________________________
-AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry):
-  fNumberOfECAClusters(0),
-  fECALocMaxCut(0),
-  fThreshold(0.01),//10 MeV
-  fGeom(geometry),
-  fRecPoints(NULL),
-  fDigitsArr(NULL)
-{
-  // ctor with the indication of the file where header Tree and digits Tree are stored
-  // use this contructor to avoid usage of Init() which uses runloader
-  // change needed by HLT - MP
-  if (!fGeom)
-  {
-    AliFatal("AliEMCALUnfolding: Geometry not initialized.");
-  }
-  
-}
-
-//____________________________________________________________________________
-AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMaxCut,Double_t *SSPars,Double_t *Par5,Double_t *Par6):
-  fNumberOfECAClusters(0),
-  fECALocMaxCut(ECALocMaxCut),
-  fThreshold(0.01),//10 MeV
-  fGeom(geometry),
-  fRecPoints(NULL),
-  fDigitsArr(NULL)
-{
-  // ctor with the indication of the file where header Tree and digits Tree are stored
-  // use this contructor to avoid usage of Init() which uses runloader
-  // change needed by HLT - MP
-  if (!fGeom)
-  {
-    AliFatal("AliEMCALUnfolding: Geometry not initialized.");
-  }
-  Int_t i=0;
-  for (i = 0; i < 8; i++) fgSSPars[i] = SSPars[i];
-  for (i = 0; i < 3; i++) {
-    fgPar5[i] = Par5[i];
-    fgPar6[i] = Par6[i];
-  }
-
-}
-
-//____________________________________________________________________________
-void AliEMCALUnfolding::Init()
-{
-  // Make all memory allocations which can not be done in default constructor.
-  // Attach the Clusterizer task to the list of EMCAL tasks
-
-  AliRunLoader *rl = AliRunLoader::Instance();
-  if (rl && rl->GetAliRun()){
-    AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
-    if(emcal)fGeom = emcal->GetGeometry();
-  }
-  
-  if(!fGeom)
-    fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
-  
-  AliDebug(1,Form("geom %p",fGeom));
-  
-  if(!gMinuit) 
-    gMinuit = new TMinuit(100) ;
-  
-}
-
-//____________________________________________________________________________
-  AliEMCALUnfolding::~AliEMCALUnfolding()
-{
-  // dtor
-}
-
-//____________________________________________________________________________
-void AliEMCALUnfolding::SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,TClonesArray *digitsArr)
-{
-  //
-  //Set input for unfolding purposes
-  SetNumberOfECAClusters(numberOfECAClusters);
-  SetRecPoints(recPoints);
-  SetDigitsArr(digitsArr);
-}
-
-//____________________________________________________________________________
-void AliEMCALUnfolding::MakeUnfolding()
-{
-  // Unfolds clusters using the shape of an ElectroMagnetic shower
-  // Performs unfolding of all clusters
-  
-  if(fNumberOfECAClusters > 0){
-    if (fGeom==0)
-      AliFatal("Did not get geometry from EMCALLoader") ;
-    //Int_t nModulesToUnfold = fGeom->GetNCells();
-    
-    Int_t numberofNotUnfolded = fNumberOfECAClusters ;
-    Int_t index ;
-    for(index = 0 ; index < numberofNotUnfolded ; index++){
-      AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
-      if(recPoint){
-        
-        Int_t nMultipl = recPoint->GetMultiplicity() ;
-        AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
-        Float_t * maxAtEnergy = new Float_t[nMultipl] ;
-        Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
-        
-        if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
-          if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){
-            fRecPoints->Remove(recPoint);
-            fRecPoints->Compress() ;//is it really needed
-            index-- ;
-            fNumberOfECAClusters-- ;
-            numberofNotUnfolded-- ;
-          }
-        }
-        else{
-          recPoint->SetNExMax(1) ; //Only one local maximum
-        }
-        
-        delete[] maxAt ;
-        delete[] maxAtEnergy ;
-      } else AliError("RecPoint NULL");
-    } // rec point loop
-  }
-  // End of Unfolding of clusters
-}
-
-//____________________________________________________________________________
-Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, 
-                                         Int_t nMax, 
-                                         AliEMCALDigit ** maxAt, 
-                                         Float_t * maxAtEnergy)
-{
-  // Extended to whole EMCAL 
-  
-  //**************************** part 1 *******************************************
-  // Performs the unfolding of a cluster with nMax overlapping showers 
-  
-  //printf("Original cluster E %f\n",iniTower->GetEnergy());
-  
-  Int_t nPar = 3 * nMax ;
-  Float_t * fitparameters = new Float_t[nPar] ;
-  
-  if (fGeom==0)
-    AliFatal("Did not get geometry from EMCALLoader") ;
-  
-  Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
-  if( !rv ) 
-  {
-    // Fit failed, return (and remove cluster? - why? I leave the cluster)
-    iniTower->SetNExMax(-1) ;
-    delete[] fitparameters ;
-    return kFALSE;
-  }
-  
-  //**************************** part 2 *******************************************
-  // create unfolded rec points and fill them with new energy lists
-  // First calculate energy deposited in each sell in accordance with
-  // fit (without fluctuations): efit[]
-  // and later correct this number in acordance with actual energy
-  // deposition
-  
-  Int_t nDigits = iniTower->GetMultiplicity() ;
-  Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells
-  Float_t xpar=0.,zpar=0.,epar=0.  ;//center of gravity in cell units
-  
-  AliEMCALDigit * digit = 0 ;
-  Int_t * digitsList = iniTower->GetDigitsList() ;
-  
-  Int_t iSupMod =  0 ;
-  Int_t iTower  =  0 ;
-  Int_t iIphi   =  0 ;
-  Int_t iIeta   =  0 ;
-  Int_t iphi    =  0 ;//x direction
-  Int_t ieta    =  0 ;//z direstion
-  
-  Int_t iparam = 0 ;
-  Int_t iDigit = 0 ;
-  
-  for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
-  {
-    digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
-    if(digit)
-    {
-      fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
-      fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
-                                         iIphi, iIeta,iphi,ieta);
-      EvalParsPhiDependence(digit->GetId(),fGeom);
-      
-      efit[iDigit] = 0.;
-      iparam = 0;
-      while(iparam < nPar )
-      {
-        xpar = fitparameters[iparam] ;
-        zpar = fitparameters[iparam+1] ;
-        epar = fitparameters[iparam+2] ;
-        iparam += 3 ;
-        
-        efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
-      }
-    } else AliDebug(1,"Digit NULL part 2!");
-    
-  }//digit loop
-  
-  //**************************** part 3 *******************************************
-  // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
-  // so that energy deposited in each cell is distributed between new clusters proportionally
-  // to its contribution to efit
-  
-  Float_t * energiesList = iniTower->GetEnergiesList() ;
-  Float_t ratio = 0 ;
-  Float_t eDigit = 0. ;
-  Int_t nSplittedClusters=(Int_t)nPar/3;
-  
-  Float_t * correctedEnergyList = new Float_t[nDigits*nSplittedClusters];
-  //above - temporary table with energies after unfolding.
-  //the orderis following: 
-  //first cluster <first cell - last cell>, 
-  //second cluster <first cell - last cell>, etc.
-  
-  //**************************** sub-part 3.1 *************************************
-  //here we check if energy of the cell in the cluster after unfolding is above threshold. 
-  //If not the energy from a given cell in the cluster is divided in correct proportions 
-  //in accordance to the other clusters and added to them and set to 0.
-  
-  iparam = 0 ;
-  while(iparam < nPar )
-  {
-    xpar = fitparameters[iparam] ;
-    zpar = fitparameters[iparam+1] ;
-    epar = fitparameters[iparam+2] ;
-    
-    for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
-    {
-      digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
-      if(digit)
-      {
-        fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
-        fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
-                                           iIphi, iIeta,iphi,ieta);
-        
-        EvalParsPhiDependence(digit->GetId(),fGeom);
-        
-        if(efit[iDigit]==0) 
-        {//just for sure
-          correctedEnergyList[iparam/3+iDigit] = 0;
-          continue;
-        }
-        
-        ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;
-        eDigit = energiesList[iDigit] * ratio ;
-        
-        //add energy to temporary matrix
-        correctedEnergyList[iparam/3+iDigit] = eDigit;
-        
-      } else AliDebug(1,"NULL digit part 3");
-    }//digit loop 
-    iparam += 3 ;
-  }//while
-  
-  //**************************** sub-part 3.2 *************************************
-  //here we correct energy for each cell and cluster
-  Float_t maximumEne=0;
-  Int_t maximumIndex=0;
-  Bool_t isAnyBelowThreshold=kFALSE;
-  //  Float_t Threshold=0.01;
-  Float_t * energyFraction = new Float_t[nSplittedClusters];
-  Int_t iparam2 = 0 ;
-  for(iDigit = 0 ; iDigit < nDigits ; iDigit++)
-  {
-    isAnyBelowThreshold=kFALSE;
-    maximumEne=0;
-    for(iparam = 0 ; iparam < nPar ; iparam+=3)
-    {
-
-      if(correctedEnergyList[iparam/3+iDigit] < fThreshold ) isAnyBelowThreshold = kTRUE;
-      if(correctedEnergyList[iparam/3+iDigit] > maximumEne) 
-      {
-        maximumEne = correctedEnergyList[iparam/3+iDigit];
-        maximumIndex = iparam;
-      }
-    }//end of loop over clusters after unfolding
-    
-    if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold 
-    
-    //printf("Correct E cell %f < %f, org Digit index %d, e = %f\n",correctedEnergyList[iparam/3+iDigit],fThreshold,iDigit, energiesList[iDigit]);
-    //if( energiesList[iDigit] < correctedEnergyList[iparam/3+iDigit]) printf("\t What? \n");
-
-    if(maximumEne < fThreshold) 
-    {//add all cluster cells and put energy into max index, other set to 0
-      maximumEne=0.;
-      for(iparam = 0 ; iparam < nPar ; iparam+=3)
-      {
-        maximumEne+=correctedEnergyList[iparam/3+iDigit];
-        correctedEnergyList[iparam/3+iDigit]=0;
-      }
-      correctedEnergyList[maximumIndex/3+iDigit]=maximumEne;
-      continue;
-    }//end if
-    
-    //divide energy of cell below threshold in the correct proportion and add to other cells
-    maximumEne=0;//not used any more so use it for the energy sum
-    for(iparam = 0 ; iparam < nPar ; iparam+=3)
-    {//calculate energy sum
-      if(correctedEnergyList[iparam/3+iDigit] < fThreshold) energyFraction[iparam/3]=0;
-      else 
-      {
-        energyFraction[iparam/3]=1;
-        maximumEne+=correctedEnergyList[iparam/3+iDigit];
-      }
-    }//end of loop over clusters after unfolding
-    if(maximumEne>0)
-    {
-      for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction
-        energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3+iDigit] / maximumEne;
-      }
-      
-      for(iparam = 0 ; iparam < nPar ; iparam+=3)
-      {//add energy from cells below threshold to others
-        if(energyFraction[iparam/3]>0) continue;
-        else
-        {
-          for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3)
-          {
-            correctedEnergyList[iparam2/3+iDigit] += (energyFraction[iparam2/3] * 
-                                                      correctedEnergyList[iparam/3+iDigit]) ;
-          }//inner loop
-          correctedEnergyList[iparam/3+iDigit] = 0;
-        }
-      }
-    }
-    else
-    {
-      //digit energy to be set to 0
-      for(iparam = 0 ; iparam < nPar ; iparam+=3)
-      {
-        correctedEnergyList[iparam/3+iDigit] = 0;
-      }
-    }//new adam correction for is energy>0
-    
-  }//end of loop over digits
-  delete[] energyFraction;
-  
-  //**************************** sub-part 3.3 *************************************
-  //here we add digits to recpoints with corrected energy
-  iparam = 0 ;
-  while(iparam < nPar )
-  {
-    AliEMCALRecPoint * recPoint = 0 ;
-    
-    if(fNumberOfECAClusters >= fRecPoints->GetSize())
-      fRecPoints->Expand(2*fNumberOfECAClusters) ;
-    
-    //add recpoint
-    (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
-    recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
-    
-    if(recPoint)
-    {
-      fNumberOfECAClusters++ ;
-      recPoint->SetNExMax(nSplittedClusters) ;
-      
-      for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
-      {
-        digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
-        
-        
-        if(digit && correctedEnergyList[iparam/3+iDigit]>0. )
-        {
-          //if(correctedEnergyList[iparam/3+iDigit]<fThreshold) printf("Final E cell %f < %f\n",correctedEnergyList[iparam/3+iDigit],fThreshold);
-          recPoint->AddDigit( *digit, correctedEnergyList[iparam/3+iDigit], kFALSE ) ; //FIXME, need to study the shared case
-        } else 
-        {
-          AliDebug(1,Form("NULL digit part3.3 or NULL energy=%f",correctedEnergyList[iparam/3+iDigit]));
-          //cout<<"nDigits "<<nDigits<<" iParam/3 "<<iparam/3<< endl;
-        }
-      }//digit loop 
-    } else AliError("NULL RecPoint");
-    
-    //protection from recpoint with no digits
-    //cout<<"multi rec "<<recPoint->GetMultiplicity()<<endl;
-    if(recPoint->GetMultiplicity()==0)
-    {
-      delete (*fRecPoints)[fNumberOfECAClusters];
-      //cout<<"size fRecPoints before "<<fRecPoints->GetSize()<<endl;
-      fRecPoints->RemoveAt(fNumberOfECAClusters);
-      //cout<<"size fRecPoints after "<<fRecPoints->GetSize()<<endl;
-      fNumberOfECAClusters--;
-      nSplittedClusters--;
-      
-    }
-    
-    iparam += 3 ;
-  }//while
-  
-  delete[] fitparameters ;
-  delete[] efit ;
-  delete[] correctedEnergyList ;
-  
-  return kTRUE;
-}
-
-
-//____________________________________________________________________________
-Bool_t AliEMCALUnfolding::UnfoldClusterV2old(AliEMCALRecPoint * iniTower, 
-                                         Int_t nMax, 
-                                         AliEMCALDigit ** maxAt, 
-                                         Float_t * maxAtEnergy)
-{
-  // Extended to whole EMCAL 
-  // Performs the unfolding of a cluster with nMax overlapping showers 
-  
-  Int_t nPar = 3 * nMax ;
-  Float_t * fitparameters = new Float_t[nPar] ;
-  
-  if (fGeom==0)
-    AliFatal("Did not get geometry from EMCALLoader") ;
-  
-  Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
-  if( !rv ) {
-    // Fit failed, return (and remove cluster? - why? I leave the cluster)
-    iniTower->SetNExMax(-1) ;
-    delete[] fitparameters ;
-    return kFALSE;
-  }
-  
-  // create unfolded rec points and fill them with new energy lists
-  // First calculate energy deposited in each sell in accordance with
-  // fit (without fluctuations): efit[]
-  // and later correct this number in acordance with actual energy
-  // deposition
-  
-  Int_t nDigits = iniTower->GetMultiplicity() ;
-  Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells
-  Float_t xpar=0.,zpar=0.,epar=0.  ;//center of gravity in cell units
-  
-  AliEMCALDigit * digit = 0 ;
-  Int_t * digitsList = iniTower->GetDigitsList() ;
-  
-  Int_t iSupMod =  0 ;
-  Int_t iTower  =  0 ;
-  Int_t iIphi   =  0 ;
-  Int_t iIeta   =  0 ;
-  Int_t iphi    =  0 ;//x direction
-  Int_t ieta    =  0 ;//z direstion
-  
-  Int_t iparam = 0 ;
-  Int_t iDigit = 0 ;
-  
-  for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
-    digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
-    if(digit){
-      fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
-      fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
-                                         iIphi, iIeta,iphi,ieta);
-      EvalParsPhiDependence(digit->GetId(),fGeom);
-      
-      efit[iDigit] = 0.;
-      iparam = 0;
-      while(iparam < nPar ){
-        xpar = fitparameters[iparam] ;
-        zpar = fitparameters[iparam+1] ;
-        epar = fitparameters[iparam+2] ;
-        iparam += 3 ;
-        
-        efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
-      }
-    } else AliError("Digit NULL!");
-    
-  }//digit loop
-  
-  // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
-  // so that energy deposited in each cell is distributed between new clusters proportionally
-  // to its contribution to efit
-  
-  Float_t * energiesList = iniTower->GetEnergiesList() ;
-  Float_t ratio = 0 ;
-  
-  iparam = 0 ;
-  while(iparam < nPar ){
-    xpar = fitparameters[iparam] ;
-    zpar = fitparameters[iparam+1] ;
-    epar = fitparameters[iparam+2] ;
-    iparam += 3 ;
-    
-    AliEMCALRecPoint * recPoint = 0 ;
-    
-    if(fNumberOfECAClusters >= fRecPoints->GetSize())
-      fRecPoints->Expand(2*fNumberOfECAClusters) ;
-    
-    //add recpoint
-    (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
-    recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
-    
-    if(recPoint){
-      
-      fNumberOfECAClusters++ ;
-      recPoint->SetNExMax((Int_t)nPar/3) ;
-      
-      Float_t eDigit = 0. ;
-      for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
-        digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
-        if(digit){
-          fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
-          fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
-                                             iIphi, iIeta,iphi,ieta);
-          EvalParsPhiDependence(digit->GetId(),fGeom);
-          if(efit[iDigit]==0) continue;//just for sure
-          ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;
-          eDigit = energiesList[iDigit] * ratio ;
-          recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case
-        } else AliError("NULL digit");
-      }//digit loop 
-    } else AliError("NULL RecPoint");
-  }//while
-  
-  delete[] fitparameters ;
-  delete[] efit ;
-  
-  return kTRUE;
-}
-
-
-//____________________________________________________________________________
-Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, 
-                                       const Float_t* maxAtEnergy,
-                                       Int_t nPar, Float_t * fitparameters) const
-{
-  // Calls TMinuit to fit the energy distribution of a cluster with several maxima
-  // The initial values for fitting procedure are set equal to the
-  // positions of local maxima.       
-  // Cluster will be fitted as a superposition of nPar/3
-  // electromagnetic showers
-
-  if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
-       
-  if(!gMinuit)
-    gMinuit = new TMinuit(100) ;//max 100 parameters
-
-  gMinuit->mncler();                     // Reset Minuit's list of paramters
-  gMinuit->SetPrintLevel(-1) ;           // No Printout
-  gMinuit->SetFCN(AliEMCALUnfolding::UnfoldingChiSquareV2) ;
-  // To set the address of the minimization function
-  TList * toMinuit = new TList();
-  toMinuit->AddAt(recPoint,0) ;
-  toMinuit->AddAt(fDigitsArr,1) ;
-  toMinuit->AddAt(fGeom,2) ;
-
-  gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
-
-  // filling initial values for fit parameters
-  AliEMCALDigit * digit ;
-
-  Int_t ierflg  = 0;
-  Int_t index   = 0 ;
-  Int_t nDigits = (Int_t) nPar / 3 ;
-
-  Int_t iDigit ;
-
-  Int_t iSupMod =  0 ;
-  Int_t iTower  =  0 ;
-  Int_t iIphi   =  0 ;
-  Int_t iIeta   =  0 ;
-  Int_t iphi    =  0 ;//x direction
-  Int_t ieta    =  0 ;//z direstion
-
-  for(iDigit = 0; iDigit < nDigits; iDigit++){
-    digit = maxAt[iDigit];
-    if(digit==0) AliError("energy of digit = 0!");
-    fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
-    fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
-                                      iIphi, iIeta,iphi,ieta);
-
-    Float_t energy = maxAtEnergy[iDigit] ;
-
-    //gMinuit->mnparm(index, "x",  iphi, 0.1, 0, 0, ierflg) ;//original
-    gMinuit->mnparm(index, "x",  iphi, 0.05, 0, 0, ierflg) ;
-    index++ ;
-    if(ierflg != 0){
-      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %d", iphi ) ;
-      toMinuit->Clear();
-      delete toMinuit ;
-      return kFALSE;
-    }
-    //gMinuit->mnparm(index, "z",  ieta, 0.1, 0, 0, ierflg) ;//original
-    gMinuit->mnparm(index, "z",  ieta, 0.05, 0, 0, ierflg) ;
-    index++ ;
-    if(ierflg != 0){
-      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %d", ieta) ;
-      toMinuit->Clear();
-      delete toMinuit ;
-      return kFALSE;
-    }
-    //gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;//original
-    gMinuit->mnparm(index, "Energy",  energy , 0.001*energy, 0., 5.*energy, ierflg) ;//was 0.05
-    index++ ;
-    if(ierflg != 0){
-      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;
-      toMinuit->Clear();
-      delete toMinuit ;
-      return kFALSE;
-    }
-  }
-
-  Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; 
-                      // The number of function call slightly depends on it.
-  //  Double_t p1 = 1.0 ;// par to gradient 
-  Double_t p2 = 0.0 ;
-  //  Double_t p3 = 3.0 ;
-  gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
-  //  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
-  gMinuit->SetMaxIterations(5);//was 5
-  gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
-  //gMinuit->mnexcm("SET PRI", &p3 , 3, ierflg) ;  // printouts
-
-  gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
-  //gMinuit->mnexcm("MINI", &p0, 0, ierflg) ;    // minimize
-  if(ierflg == 4){  // Minimum not found
-    AliDebug(1,"EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;
-    toMinuit->Clear();
-    delete toMinuit ;
-    return kFALSE ;
-  }
-  for(index = 0; index < nPar; index++){
-    Double_t err = 0. ;
-    Double_t val = 0. ;
-    gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
-    fitparameters[index] = val ;
-  }
-
-  toMinuit->Clear();
-  delete toMinuit ;
-  return kTRUE;
-
-}
-
-//____________________________________________________________________________
-Double_t  AliEMCALUnfolding::ShowerShapeV2(Double_t x, Double_t y)
-{ 
-  // extended to whole EMCAL 
-  // Shape of the shower
-  // If you change this function, change also the gradient evaluation in ChiSquare()
-
-  Double_t r = fgSSPars[7]*TMath::Sqrt(x*x+y*y);
-  Double_t rp1  = TMath::Power(r, fgSSPars[1]) ;
-  Double_t rp5  = TMath::Power(r, fgSSPars[5]) ;
-  Double_t shape = fgSSPars[0]*TMath::Exp( -rp1 * (1. / (fgSSPars[2] + fgSSPars[3] * rp1) + fgSSPars[4] / (1 + fgSSPars[6] * rp5) ) ) ;
-  return shape ;
-}
-
-//____________________________________________________________________________
-void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,
-                                            Double_t & fret,
-                                            Double_t * x, Int_t iflag)
-{
-  // Calculates the Chi square for the cluster unfolding minimization
-  // Number of parameters, Gradient, Chi squared, parameters, what to do
-  
-  TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
-  if(toMinuit){
-    AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) )  ;
-    TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
-    // A bit buggy way to get an access to the geometry
-    // To be revised!
-    AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));
-    
-    if(recPoint && digits && geom){
-      
-      Int_t * digitsList     = recPoint->GetDigitsList() ;
-      
-      Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;
-      
-      Float_t * energiesList = recPoint->GetEnergiesList() ;
-      
-      fret = 0. ;
-      Int_t iparam = 0 ;
-      
-      if(iflag == 2)
-        for(iparam = 0 ; iparam < nPar ; iparam++)
-          Grad[iparam] = 0 ; // Will evaluate gradient
-      
-      Double_t efit = 0. ;
-      
-      AliEMCALDigit * digit ;
-      Int_t iDigit ;
-      
-      Int_t iSupMod =  0 ;
-      Int_t iTower  =  0 ;
-      Int_t iIphi   =  0 ;
-      Int_t iIeta   =  0 ;
-      Int_t iphi    =  0 ;//x direction
-      Int_t ieta    =  0 ;//z direstion
-      
-      
-      for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
-        if(energiesList[iDigit]==0) continue;
-        
-        digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );
-        
-        if(digit){
-        geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
-        geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
-                                          iIphi, iIeta,iphi,ieta);
-        EvalParsPhiDependence(digit->GetId(),geom);
-        
-        if(iflag == 2){  // calculate gradient
-          Int_t iParam = 0 ;
-          efit = 0. ;
-          while(iParam < nPar ){
-            Double_t dx = ((Float_t)iphi - x[iParam]) ;
-            iParam++ ;
-            Double_t dz = ((Float_t)ieta - x[iParam]) ;
-            iParam++ ;
-            efit += x[iParam] * ShowerShapeV2(dx,dz) ;
-            iParam++ ;
-          }
-          
-          Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)
-          iParam = 0 ;
-          while(iParam < nPar ){
-            Double_t xpar = x[iParam] ;
-            Double_t zpar = x[iParam+1] ;
-            Double_t epar = x[iParam+2] ;
-            
-            Double_t dr = fgSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) );
-            Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
-            Double_t rp1  = TMath::Power(dr, fgSSPars[1]) ;
-            Double_t rp5  = TMath::Power(dr, fgSSPars[5]) ;
-            
-            Double_t deriv = -2 * TMath::Power(dr,fgSSPars[1]-2.) * fgSSPars[7] * fgSSPars[7] * 
-            (fgSSPars[1] * ( 1/(fgSSPars[2]+fgSSPars[3]*rp1) + fgSSPars[4]/(1+fgSSPars[6]*rp5) ) - 
-             (fgSSPars[1]*fgSSPars[3]*rp1/( (fgSSPars[2]+fgSSPars[3]*rp1)*(fgSSPars[2]+fgSSPars[3]*rp1) ) + 
-              fgSSPars[4]*fgSSPars[5]*fgSSPars[6]*rp5/( (1+fgSSPars[6]*rp5)*(1+fgSSPars[6]*rp5) ) ) );
-            
-            //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) )
-            //                                                   - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;
-            
-            Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ;  // Derivative over x
-            iParam++ ;
-            Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ;  // Derivative over z
-            iParam++ ;
-            Grad[iParam] += shape ;                                  // Derivative over energy
-            iParam++ ;
-          }
-        }
-        efit = 0;
-        iparam = 0 ;
-        
-        while(iparam < nPar ){
-          Double_t xpar = x[iparam] ;
-          Double_t zpar = x[iparam+1] ;
-          Double_t epar = x[iparam+2] ;
-          iparam += 3 ;
-          efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
-        }
-        
-        fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;
-        // Here we assume, that sigma = sqrt(E) 
-        } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!, nPar %d \n", nPar); // put nPar here to cheat coverity and rule checker
-      } // digit loop
-    } // recpoint, digits and geom not NULL
-  }// List is not NULL
-  
-}
-
-
-//____________________________________________________________________________
-void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){
-  for(UInt_t i=0;i<7;++i)
-    fgSSPars[i]=pars[i];
-  if(pars[2]==0. && pars[3]==0.) fgSSPars[2]=1.;//to avoid dividing by 0
-}
-
-//____________________________________________________________________________
-void AliEMCALUnfolding::SetPar5(Double_t *pars){
-  for(UInt_t i=0;i<3;++i)
-    fgPar5[i]=pars[i];
-}
-
-//____________________________________________________________________________
-void AliEMCALUnfolding::SetPar6(Double_t *pars){
-  for(UInt_t i=0;i<3;++i)
-    fgPar6[i]=pars[i];
-}
-
-//____________________________________________________________________________
-void AliEMCALUnfolding::EvalPar5(Double_t phi){
-  //
-  //Evaluate the 5th parameter of the shower shape function
-  //phi in degrees range (-10,10)
-  //
-  //fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936;
-  fgSSPars[5] = fgPar5[0] + phi * fgPar5[1] + phi*phi * fgPar5[2];
-}
-
-//____________________________________________________________________________
-void AliEMCALUnfolding::EvalPar6(Double_t phi){
-  //
-  //Evaluate the 6th parameter of the shower shape function
-  //phi in degrees range (-10,10)
-  //
-  //fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361;
-  fgSSPars[6] = fgPar6[0] + phi * fgPar6[1] + phi*phi * fgPar6[2];
-}
-
-//____________________________________________________________________________
-void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, const AliEMCALGeometry *geom){
-  //
-  // calculate params p5 and p6 depending on the phi angle in global coordinate
-  // for the cell with given absId index
-  //
-  Double_t etaGlob = 0.;//eta in global c.s. - unused
-  Double_t phiGlob = 0.;//phi in global c.s. in radians
-  geom->EtaPhiFromIndex(absId, etaGlob, phiGlob);
-  phiGlob*=180./TMath::Pi();
-  phiGlob-=90.;
-  phiGlob-= (Double_t)((Int_t)geom->GetSuperModuleNumber(absId)/2 * 20);
-
-  EvalPar5(phiGlob);
-  EvalPar6(phiGlob);
-}
-
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: The ALICE Off-line Project.                                    *\r
+ * Contributors are mentioned in the code where appropriate.              *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+//_________________________________________________________________________\r
+//  Base class for the cluster unfolding algorithm \r
+//*-- Author: Adam Matyja (SUBATECH)\r
+//  Based on unfolding in clusterizerv1 done by Cynthia Hadjidakis\r
+//-- Unfolding for eta~0: Cynthia Hadjidakis - still in AliEMCALCLusterizerv1\r
+//-- Unfolding extension for whole EMCAL: Adam Matyja (SUBATECH & INP PAN)\r
+//\r
+//  unfolds the clusters having several local maxima. \r
+//////////////////////////////////////////////////////////////////////////////\r
+\r
+// --- ROOT system ---\r
+#include "TClonesArray.h"\r
+#include <TMath.h> \r
+#include <TMinuit.h>\r
+\r
+// --- Standard library ---\r
+#include <cassert>\r
+\r
+// --- AliRoot header files ---\r
+#include "AliEMCALUnfolding.h"\r
+#include "AliEMCALGeometry.h"\r
+#include "AliRunLoader.h"\r
+#include "AliRun.h"\r
+#include "AliEMCAL.h"\r
+#include "AliEMCALRecParam.h"\r
+#include "AliEMCALRecPoint.h"\r
+#include "AliEMCALDigit.h"\r
+#include "AliEMCALReconstructor.h"\r
+\r
+#include "AliLog.h"\r
+#include "AliCDBManager.h"\r
+class AliCDBStorage;\r
+#include "AliCDBEntry.h"\r
+\r
+Double_t AliEMCALUnfolding::fgSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};\r
+Double_t AliEMCALUnfolding::fgPar5[3]={12.31,-0.007381,-0.06936};\r
+Double_t AliEMCALUnfolding::fgPar6[3]={0.05452,0.0001228,0.001361};\r
+\r
+ClassImp(AliEMCALUnfolding)\r
+  \r
+//____________________________________________________________________________\r
+AliEMCALUnfolding::AliEMCALUnfolding():\r
+  fNumberOfECAClusters(0),\r
+  fECALocMaxCut(0),\r
+  fThreshold(0.01),//10 MeV\r
+  fGeom(NULL),\r
+  fRecPoints(NULL),\r
+  fDigitsArr(NULL)\r
+{\r
+  // ctor with the indication of the file where header Tree and digits Tree are stored\r
\r
+  Init() ;\r
+}\r
+\r
+//____________________________________________________________________________\r
+AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry):\r
+  fNumberOfECAClusters(0),\r
+  fECALocMaxCut(0),\r
+  fThreshold(0.01),//10 MeV\r
+  fGeom(geometry),\r
+  fRecPoints(NULL),\r
+  fDigitsArr(NULL)\r
+{\r
+  // ctor with the indication of the file where header Tree and digits Tree are stored\r
+  // use this contructor to avoid usage of Init() which uses runloader\r
+  // change needed by HLT - MP\r
+  if (!fGeom)\r
+  {\r
+    AliFatal("AliEMCALUnfolding: Geometry not initialized.");\r
+  }\r
+  \r
+}\r
+\r
+//____________________________________________________________________________\r
+AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMaxCut,Double_t *SSPars,Double_t *Par5,Double_t *Par6):\r
+  fNumberOfECAClusters(0),\r
+  fECALocMaxCut(ECALocMaxCut),\r
+  fThreshold(0.01),//10 MeV\r
+  fGeom(geometry),\r
+  fRecPoints(NULL),\r
+  fDigitsArr(NULL)\r
+{\r
+  // ctor with the indication of the file where header Tree and digits Tree are stored\r
+  // use this contructor to avoid usage of Init() which uses runloader\r
+  // change needed by HLT - MP\r
+  if (!fGeom)\r
+  {\r
+    AliFatal("AliEMCALUnfolding: Geometry not initialized.");\r
+  }\r
+  Int_t i=0;\r
+  for (i = 0; i < 8; i++) fgSSPars[i] = SSPars[i];\r
+  for (i = 0; i < 3; i++) {\r
+    fgPar5[i] = Par5[i];\r
+    fgPar6[i] = Par6[i];\r
+  }\r
+\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::Init()\r
+{\r
+  // Make all memory allocations which can not be done in default constructor.\r
+  // Attach the Clusterizer task to the list of EMCAL tasks\r
+\r
+  AliRunLoader *rl = AliRunLoader::Instance();\r
+  if (rl && rl->GetAliRun()){\r
+    AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));\r
+    if(emcal)fGeom = emcal->GetGeometry();\r
+  }\r
+  \r
+  if(!fGeom)\r
+    fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());\r
+  \r
+  AliDebug(1,Form("geom %p",fGeom));\r
+  \r
+  if(!gMinuit) \r
+    gMinuit = new TMinuit(100) ;\r
+  \r
+}\r
+\r
+//____________________________________________________________________________\r
+  AliEMCALUnfolding::~AliEMCALUnfolding()\r
+{\r
+  // dtor\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,TClonesArray *digitsArr)\r
+{\r
+  //\r
+  //Set input for unfolding purposes\r
+  SetNumberOfECAClusters(numberOfECAClusters);\r
+  SetRecPoints(recPoints);\r
+  SetDigitsArr(digitsArr);\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::MakeUnfolding()\r
+{\r
+  // Unfolds clusters using the shape of an ElectroMagnetic shower\r
+  // Performs unfolding of all clusters\r
+  \r
+  if(fNumberOfECAClusters > 0){\r
+    if (fGeom==0)\r
+      AliFatal("Did not get geometry from EMCALLoader") ;\r
+    //Int_t nModulesToUnfold = fGeom->GetNCells();\r
+    \r
+    Int_t numberofNotUnfolded = fNumberOfECAClusters ;\r
+    Int_t index ;\r
+    for(index = 0 ; index < numberofNotUnfolded ; index++){\r
+      AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;\r
+      if(recPoint){\r
+        \r
+        Int_t nMultipl = recPoint->GetMultiplicity() ;\r
+        AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;\r
+        Float_t * maxAtEnergy = new Float_t[nMultipl] ;\r
+        Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;\r
+        \r
+        if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0\r
+          if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){\r
+            fRecPoints->Remove(recPoint);\r
+            fRecPoints->Compress() ;//is it really needed\r
+            index-- ;\r
+            fNumberOfECAClusters-- ;\r
+            numberofNotUnfolded-- ;\r
+          }\r
+        }\r
+        else{\r
+          recPoint->SetNExMax(1) ; //Only one local maximum\r
+        }\r
+        \r
+        delete[] maxAt ;\r
+        delete[] maxAtEnergy ;\r
+      } else AliError("RecPoint NULL");\r
+    } // rec point loop\r
+  }\r
+  // End of Unfolding of clusters\r
+}\r
+\r
+//____________________________________________________________________________\r
+Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, \r
+                                         Int_t nMax, \r
+                                         AliEMCALDigit ** maxAt, \r
+                                         Float_t * maxAtEnergy)\r
+{\r
+  // Extended to whole EMCAL \r
+  \r
+  //**************************** part 1 *******************************************\r
+  // Performs the unfolding of a cluster with nMax overlapping showers \r
+  \r
+  //printf("Original cluster E %f\n",iniTower->GetEnergy());\r
+  \r
+  Int_t nPar = 3 * nMax ;\r
+  Float_t * fitparameters = new Float_t[nPar] ;\r
+  \r
+  if (fGeom==0)\r
+    AliFatal("Did not get geometry from EMCALLoader") ;\r
+  \r
+  Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;\r
+  if( !rv ) \r
+  {\r
+    // Fit failed, return (and remove cluster? - why? I leave the cluster)\r
+    iniTower->SetNExMax(-1) ;\r
+    delete[] fitparameters ;\r
+    return kFALSE;\r
+  }\r
+  \r
+  //**************************** part 2 *******************************************\r
+  // create unfolded rec points and fill them with new energy lists\r
+  // First calculate energy deposited in each sell in accordance with\r
+  // fit (without fluctuations): efit[]\r
+  // and later correct this number in acordance with actual energy\r
+  // deposition\r
+  \r
+  Int_t nDigits = iniTower->GetMultiplicity() ;\r
+  Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells\r
+  Float_t xpar=0.,zpar=0.,epar=0.  ;//center of gravity in cell units\r
+  \r
+  AliEMCALDigit * digit = 0 ;\r
+  Int_t * digitsList = iniTower->GetDigitsList() ;\r
+  \r
+  Int_t iSupMod =  0 ;\r
+  Int_t iTower  =  0 ;\r
+  Int_t iIphi   =  0 ;\r
+  Int_t iIeta   =  0 ;\r
+  Int_t iphi    =  0 ;//x direction\r
+  Int_t ieta    =  0 ;//z direstion\r
+  \r
+  Int_t iparam = 0 ;\r
+  Int_t iDigit = 0 ;\r
+  \r
+  for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)\r
+  {\r
+    digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;\r
+    if(digit)\r
+    {\r
+      fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+      fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+                                         iIphi, iIeta,iphi,ieta);\r
+      EvalParsPhiDependence(digit->GetId(),fGeom);\r
+      \r
+      efit[iDigit] = 0.;\r
+      iparam = 0;\r
+      while(iparam < nPar )\r
+      {\r
+        xpar = fitparameters[iparam] ;\r
+        zpar = fitparameters[iparam+1] ;\r
+        epar = fitparameters[iparam+2] ;\r
+\r
+        efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
+\r
+        iparam += 3 ;\r
+      }\r
+\r
+    } else AliDebug(1,"Digit NULL part 2!");\r
+    \r
+  }//digit loop\r
+  \r
+  //**************************** part 3 *******************************************\r
+  // Now create new RecPoints and fill energy lists with efit corrected to fluctuations\r
+  // so that energy deposited in each cell is distributed between new clusters proportionally\r
+  // to its contribution to efit\r
+  \r
+  Float_t * energiesList = iniTower->GetEnergiesList() ;\r
+  Float_t ratio = 0 ;\r
+  Float_t eDigit = 0. ;\r
+  Int_t nSplittedClusters=(Int_t)nPar/3;\r
+  \r
+  Float_t * correctedEnergyList = new Float_t[nDigits*nSplittedClusters];\r
+  //above - temporary table with energies after unfolding.\r
+  //the order is following: \r
+  //first cluster <first cell - last cell>, \r
+  //second cluster <first cell - last cell>, etc.\r
+  \r
+  //**************************** sub-part 3.1 *************************************\r
+  //here we check if energy of the cell in the cluster after unfolding is above threshold. \r
+  //If not the energy from a given cell in the cluster is divided in correct proportions \r
+  //in accordance to the other clusters and added to them and set to 0.\r
+  \r
+  iparam = 0 ;\r
+  while(iparam < nPar )\r
+  {\r
+    xpar = fitparameters[iparam] ;\r
+    zpar = fitparameters[iparam+1] ;\r
+    epar = fitparameters[iparam+2] ;\r
+    \r
+    for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)\r
+    {\r
+      digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r
+      if(digit)\r
+      {\r
+        fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+        fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+                                           iIphi, iIeta,iphi,ieta);\r
+        \r
+        EvalParsPhiDependence(digit->GetId(),fGeom);\r
+       \r
\r
+        if(efit[iDigit]==0) \r
+        {//just for sure\r
+          correctedEnergyList[iparam/3*nDigits+iDigit] = 0;//correction here\r
+          continue;\r
+        }\r
+        \r
+        ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;\r
+        eDigit = energiesList[iDigit] * ratio ;\r
+        \r
+        //add energy to temporary matrix\r
+        correctedEnergyList[iparam/3*nDigits+iDigit] = eDigit;\r
+        \r
+      } else AliDebug(1,"NULL digit part 3");\r
+    }//digit loop \r
+    iparam += 3 ;\r
+  }//while\r
+  \r
+  //**************************** sub-part 3.2 *************************************\r
+  //here we correct energy for each cell and cluster\r
+  Float_t maximumEne=0;\r
+  Int_t maximumIndex=0;\r
+  Bool_t isAnyBelowThreshold=kFALSE;\r
+  //  Float_t Threshold=0.01;\r
+  Float_t * energyFraction = new Float_t[nSplittedClusters];\r
+  Int_t iparam2 = 0 ;\r
+  for(iDigit = 0 ; iDigit < nDigits ; iDigit++)\r
+  {\r
+    isAnyBelowThreshold=kFALSE;\r
+    maximumEne=0;\r
+    for(iparam = 0 ; iparam < nPar ; iparam+=3)\r
+    {\r
+\r
+      if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold ) isAnyBelowThreshold = kTRUE;\r
+      if(correctedEnergyList[iparam/3*nDigits+iDigit] > maximumEne) \r
+      {\r
+        maximumEne = correctedEnergyList[iparam/3*nDigits+iDigit];\r
+        maximumIndex = iparam;\r
+      }\r
+    }//end of loop over clusters after unfolding\r
+    \r
+    if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold \r
+    \r
+    //    printf("Correct E cell %f < %f, org Digit index %d, e = %f\n",correctedEnergyList[iparam/3*nDigits+iDigit],fThreshold,iDigit, energiesList[iDigit]);\r
+    //if( energiesList[iDigit] < correctedEnergyList[iparam/3*nDigits+iDigit]) printf("\t What? \n");\r
+\r
+    if(maximumEne < fThreshold) \r
+    {//add all cluster cells and put energy into max index, other set to 0\r
+      maximumEne=0.;\r
+      for(iparam = 0 ; iparam < nPar ; iparam+=3)\r
+      {\r
+        maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];\r
+        correctedEnergyList[iparam/3*nDigits+iDigit]=0;\r
+      }\r
+      correctedEnergyList[maximumIndex/3*nDigits+iDigit]=maximumEne;\r
+      continue;\r
+    }//end if\r
+    \r
+    //divide energy of cell below threshold in the correct proportion and add to other cells\r
+    maximumEne=0;//not used any more so use it for the energy sum\r
+    for(iparam = 0 ; iparam < nPar ; iparam+=3)\r
+    {//calculate energy sum\r
+      if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold) energyFraction[iparam/3]=0;\r
+      else \r
+      {\r
+        energyFraction[iparam/3]=1;\r
+        maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];\r
+      }\r
+    }//end of loop over clusters after unfolding\r
+    if(maximumEne>0)\r
+    {\r
+      for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction\r
+        energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3*nDigits+iDigit] / maximumEne;\r
+      }\r
+      \r
+      for(iparam = 0 ; iparam < nPar ; iparam+=3)\r
+      {//add energy from cells below threshold to others\r
+        if(energyFraction[iparam/3]>0) continue;\r
+        else\r
+        {\r
+          for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3)\r
+          {\r
+            correctedEnergyList[iparam2/3*nDigits+iDigit] += (energyFraction[iparam2/3] * \r
+                                                      correctedEnergyList[iparam/3*nDigits+iDigit]) ;\r
+          }//inner loop\r
+          correctedEnergyList[iparam/3*nDigits+iDigit] = 0;\r
+        }\r
+      }\r
+    }\r
+    else\r
+    {\r
+      //digit energy to be set to 0\r
+      for(iparam = 0 ; iparam < nPar ; iparam+=3)\r
+      {\r
+        correctedEnergyList[iparam/3*nDigits+iDigit] = 0;\r
+      }\r
+    }//new adam correction for is energy>0\r
+    \r
+  }//end of loop over digits\r
+  delete[] energyFraction;\r
+  \r
+  //**************************** sub-part 3.3 *************************************\r
+  //here we add digits to recpoints with corrected energy\r
+  iparam = 0 ;\r
+  while(iparam < nPar )\r
+  {\r
+    AliEMCALRecPoint * recPoint = 0 ;\r
+    \r
+    if(fNumberOfECAClusters >= fRecPoints->GetSize())\r
+      fRecPoints->Expand(2*fNumberOfECAClusters) ;\r
+    \r
+    //add recpoint\r
+    (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;\r
+    recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\r
+    \r
+    if(recPoint)\r
+    {\r
+      fNumberOfECAClusters++ ;\r
+      recPoint->SetNExMax(nSplittedClusters) ;\r
+      \r
+      for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)\r
+      {\r
+        digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r
+        \r
+        \r
+        if(digit && correctedEnergyList[iparam/3*nDigits+iDigit]>0. )\r
+        {\r
+          if(correctedEnergyList[iparam/3*nDigits+iDigit]<fThreshold) printf("Final E cell %f < %f\n",correctedEnergyList[iparam/3*nDigits+iDigit],fThreshold);\r
+          recPoint->AddDigit( *digit, correctedEnergyList[iparam/3*nDigits+iDigit], kFALSE ) ; //FIXME, need to study the shared case\r
+        } else \r
+        {\r
+          AliDebug(1,Form("NULL digit part3.3 or NULL energy=%f",correctedEnergyList[iparam/3*nDigits+iDigit]));\r
+        }\r
+      }//digit loop \r
+    } else AliError("NULL RecPoint");\r
+    \r
+    //protection from recpoint with no digits\r
+    if(recPoint->GetMultiplicity()==0)\r
+    {\r
+      delete (*fRecPoints)[fNumberOfECAClusters];\r
+      //cout<<"size fRecPoints before "<<fRecPoints->GetSize()<<endl;\r
+      fRecPoints->RemoveAt(fNumberOfECAClusters);\r
+      //cout<<"size fRecPoints after "<<fRecPoints->GetSize()<<endl;\r
+      fNumberOfECAClusters--;\r
+      nSplittedClusters--;\r
+      \r
+    }\r
+    \r
+    iparam += 3 ;\r
+  }//while\r
+  \r
+  delete[] fitparameters ;\r
+  delete[] efit ;\r
+  delete[] correctedEnergyList ;\r
+  \r
+  return kTRUE;\r
+}\r
+\r
+\r
+//____________________________________________________________________________\r
+Bool_t AliEMCALUnfolding::UnfoldClusterV2old(AliEMCALRecPoint * iniTower, \r
+                                         Int_t nMax, \r
+                                         AliEMCALDigit ** maxAt, \r
+                                         Float_t * maxAtEnergy)\r
+{\r
+  // Extended to whole EMCAL \r
+  // Performs the unfolding of a cluster with nMax overlapping showers \r
+  \r
+  Int_t nPar = 3 * nMax ;\r
+  Float_t * fitparameters = new Float_t[nPar] ;\r
+  \r
+  if (fGeom==0)\r
+    AliFatal("Did not get geometry from EMCALLoader") ;\r
+  \r
+  Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;\r
+  if( !rv ) {\r
+    // Fit failed, return (and remove cluster? - why? I leave the cluster)\r
+    iniTower->SetNExMax(-1) ;\r
+    delete[] fitparameters ;\r
+    return kFALSE;\r
+  }\r
+  \r
+  // create unfolded rec points and fill them with new energy lists\r
+  // First calculate energy deposited in each sell in accordance with\r
+  // fit (without fluctuations): efit[]\r
+  // and later correct this number in acordance with actual energy\r
+  // deposition\r
+  \r
+  Int_t nDigits = iniTower->GetMultiplicity() ;\r
+  Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells\r
+  Float_t xpar=0.,zpar=0.,epar=0.  ;//center of gravity in cell units\r
+  \r
+  AliEMCALDigit * digit = 0 ;\r
+  Int_t * digitsList = iniTower->GetDigitsList() ;\r
+  \r
+  Int_t iSupMod =  0 ;\r
+  Int_t iTower  =  0 ;\r
+  Int_t iIphi   =  0 ;\r
+  Int_t iIeta   =  0 ;\r
+  Int_t iphi    =  0 ;//x direction\r
+  Int_t ieta    =  0 ;//z direstion\r
+  \r
+  Int_t iparam = 0 ;\r
+  Int_t iDigit = 0 ;\r
+  \r
+  for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
+    digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;\r
+    if(digit){\r
+      fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+      fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+                                         iIphi, iIeta,iphi,ieta);\r
+      EvalParsPhiDependence(digit->GetId(),fGeom);\r
+      \r
+      efit[iDigit] = 0.;\r
+      iparam = 0;\r
+      while(iparam < nPar ){\r
+        xpar = fitparameters[iparam] ;\r
+        zpar = fitparameters[iparam+1] ;\r
+        epar = fitparameters[iparam+2] ;\r
+        iparam += 3 ;\r
+        \r
+        efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
+      }\r
+    } else AliError("Digit NULL!");\r
+    \r
+  }//digit loop\r
+  \r
+  // Now create new RecPoints and fill energy lists with efit corrected to fluctuations\r
+  // so that energy deposited in each cell is distributed between new clusters proportionally\r
+  // to its contribution to efit\r
+  \r
+  Float_t * energiesList = iniTower->GetEnergiesList() ;\r
+  Float_t ratio = 0 ;\r
+  \r
+  iparam = 0 ;\r
+  while(iparam < nPar ){\r
+    xpar = fitparameters[iparam] ;\r
+    zpar = fitparameters[iparam+1] ;\r
+    epar = fitparameters[iparam+2] ;\r
+    iparam += 3 ;\r
+    \r
+    AliEMCALRecPoint * recPoint = 0 ;\r
+    \r
+    if(fNumberOfECAClusters >= fRecPoints->GetSize())\r
+      fRecPoints->Expand(2*fNumberOfECAClusters) ;\r
+    \r
+    //add recpoint\r
+    (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;\r
+    recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\r
+    \r
+    if(recPoint){\r
+      \r
+      fNumberOfECAClusters++ ;\r
+      recPoint->SetNExMax((Int_t)nPar/3) ;\r
+      \r
+      Float_t eDigit = 0. ;\r
+      for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
+        digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r
+        if(digit){\r
+          fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+          fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+                                             iIphi, iIeta,iphi,ieta);\r
+          EvalParsPhiDependence(digit->GetId(),fGeom);\r
+          if(efit[iDigit]==0) continue;//just for sure\r
+          ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;\r
+          eDigit = energiesList[iDigit] * ratio ;\r
+          recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case\r
+        } else AliError("NULL digit");\r
+      }//digit loop \r
+    } else AliError("NULL RecPoint");\r
+  }//while\r
+  \r
+  delete[] fitparameters ;\r
+  delete[] efit ;\r
+  \r
+  return kTRUE;\r
+}\r
+\r
+\r
+//____________________________________________________________________________\r
+Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, \r
+                                       const Float_t* maxAtEnergy,\r
+                                       Int_t nPar, Float_t * fitparameters) const\r
+{\r
+  // Calls TMinuit to fit the energy distribution of a cluster with several maxima\r
+  // The initial values for fitting procedure are set equal to the\r
+  // positions of local maxima.       \r
+  // Cluster will be fitted as a superposition of nPar/3\r
+  // electromagnetic showers\r
+\r
+  if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");\r
+       \r
+  if(!gMinuit)\r
+    gMinuit = new TMinuit(100) ;//max 100 parameters\r
+\r
+  gMinuit->mncler();                     // Reset Minuit's list of paramters\r
+  gMinuit->SetPrintLevel(-1) ;           // No Printout\r
+  gMinuit->SetFCN(AliEMCALUnfolding::UnfoldingChiSquareV2) ;\r
+  // To set the address of the minimization function\r
+  TList * toMinuit = new TList();\r
+  toMinuit->AddAt(recPoint,0) ;\r
+  toMinuit->AddAt(fDigitsArr,1) ;\r
+  toMinuit->AddAt(fGeom,2) ;\r
+\r
+  gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare\r
+\r
+  // filling initial values for fit parameters\r
+  AliEMCALDigit * digit ;\r
+\r
+  Int_t ierflg  = 0;\r
+  Int_t index   = 0 ;\r
+  Int_t nDigits = (Int_t) nPar / 3 ;\r
+\r
+  Int_t iDigit ;\r
+\r
+  Int_t iSupMod =  0 ;\r
+  Int_t iTower  =  0 ;\r
+  Int_t iIphi   =  0 ;\r
+  Int_t iIeta   =  0 ;\r
+  Int_t iphi    =  0 ;//x direction\r
+  Int_t ieta    =  0 ;//z direstion\r
+\r
+  for(iDigit = 0; iDigit < nDigits; iDigit++){\r
+    digit = maxAt[iDigit];\r
+    if(digit==0) AliError("energy of digit = 0!");\r
+    fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+    fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+                                      iIphi, iIeta,iphi,ieta);\r
+\r
+    Float_t energy = maxAtEnergy[iDigit] ;\r
+\r
+    //gMinuit->mnparm(index, "x",  iphi, 0.1, 0, 0, ierflg) ;//original\r
+    gMinuit->mnparm(index, "x",  iphi, 0.05, 0, 0, ierflg) ;\r
+    index++ ;\r
+    if(ierflg != 0){\r
+      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %d", iphi ) ;\r
+      toMinuit->Clear();\r
+      delete toMinuit ;\r
+      return kFALSE;\r
+    }\r
+    //gMinuit->mnparm(index, "z",  ieta, 0.1, 0, 0, ierflg) ;//original\r
+    gMinuit->mnparm(index, "z",  ieta, 0.05, 0, 0, ierflg) ;\r
+    index++ ;\r
+    if(ierflg != 0){\r
+      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %d", ieta) ;\r
+      toMinuit->Clear();\r
+      delete toMinuit ;\r
+      return kFALSE;\r
+    }\r
+    //gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;//original\r
+    gMinuit->mnparm(index, "Energy",  energy , 0.001*energy, 0., 5.*energy, ierflg) ;//was 0.05\r
+    index++ ;\r
+    if(ierflg != 0){\r
+      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;\r
+      toMinuit->Clear();\r
+      delete toMinuit ;\r
+      return kFALSE;\r
+    }\r
+  }\r
+\r
+  Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; \r
+                      // The number of function call slightly depends on it.\r
+  //  Double_t p1 = 1.0 ;// par to gradient \r
+  Double_t p2 = 0.0 ;\r
+  //  Double_t p3 = 3.0 ;\r
+  gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls\r
+  //  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient\r
+  gMinuit->SetMaxIterations(5);//was 5\r
+  gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings\r
+  //gMinuit->mnexcm("SET PRI", &p3 , 3, ierflg) ;  // printouts\r
+\r
+  gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize\r
+  //gMinuit->mnexcm("MINI", &p0, 0, ierflg) ;    // minimize\r
+  if(ierflg == 4){  // Minimum not found\r
+    AliDebug(1,"EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;\r
+    toMinuit->Clear();\r
+    delete toMinuit ;\r
+    return kFALSE ;\r
+  }\r
+  for(index = 0; index < nPar; index++){\r
+    Double_t err = 0. ;\r
+    Double_t val = 0. ;\r
+    gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index\r
+    fitparameters[index] = val ;\r
+  }\r
+\r
+  toMinuit->Clear();\r
+  delete toMinuit ;\r
+  return kTRUE;\r
+\r
+}\r
+\r
+//____________________________________________________________________________\r
+Double_t  AliEMCALUnfolding::ShowerShapeV2(Double_t x, Double_t y)\r
+{ \r
+  // extended to whole EMCAL \r
+  // Shape of the shower\r
+  // If you change this function, change also the gradient evaluation in ChiSquare()\r
+\r
+  Double_t r = fgSSPars[7]*TMath::Sqrt(x*x+y*y);\r
+  Double_t rp1  = TMath::Power(r, fgSSPars[1]) ;\r
+  Double_t rp5  = TMath::Power(r, fgSSPars[5]) ;\r
+  Double_t shape = fgSSPars[0]*TMath::Exp( -rp1 * (1. / (fgSSPars[2] + fgSSPars[3] * rp1) + fgSSPars[4] / (1 + fgSSPars[6] * rp5) ) ) ;\r
+  return shape ;\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,\r
+                                            Double_t & fret,\r
+                                            Double_t * x, Int_t iflag)\r
+{\r
+  // Calculates the Chi square for the cluster unfolding minimization\r
+  // Number of parameters, Gradient, Chi squared, parameters, what to do\r
+  \r
+  TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;\r
+  if(toMinuit){\r
+    AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) )  ;\r
+    TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;\r
+    // A bit buggy way to get an access to the geometry\r
+    // To be revised!\r
+    AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));\r
+    \r
+    if(recPoint && digits && geom){\r
+      \r
+      Int_t * digitsList     = recPoint->GetDigitsList() ;\r
+      \r
+      Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;\r
+      \r
+      Float_t * energiesList = recPoint->GetEnergiesList() ;\r
+      \r
+      fret = 0. ;\r
+      Int_t iparam = 0 ;\r
+      \r
+      if(iflag == 2)\r
+        for(iparam = 0 ; iparam < nPar ; iparam++)\r
+          Grad[iparam] = 0 ; // Will evaluate gradient\r
+      \r
+      Double_t efit = 0. ;\r
+      \r
+      AliEMCALDigit * digit ;\r
+      Int_t iDigit ;\r
+      \r
+      Int_t iSupMod =  0 ;\r
+      Int_t iTower  =  0 ;\r
+      Int_t iIphi   =  0 ;\r
+      Int_t iIeta   =  0 ;\r
+      Int_t iphi    =  0 ;//x direction\r
+      Int_t ieta    =  0 ;//z direstion\r
+      \r
+      \r
+      for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {\r
+        if(energiesList[iDigit]==0) continue;\r
+        \r
+        digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );\r
+        \r
+        if(digit){\r
+        geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+        geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+                                          iIphi, iIeta,iphi,ieta);\r
+        EvalParsPhiDependence(digit->GetId(),geom);\r
+        \r
+        if(iflag == 2){  // calculate gradient\r
+          Int_t iParam = 0 ;\r
+          efit = 0. ;\r
+          while(iParam < nPar ){\r
+            Double_t dx = ((Float_t)iphi - x[iParam]) ;\r
+            iParam++ ;\r
+            Double_t dz = ((Float_t)ieta - x[iParam]) ;\r
+            iParam++ ;\r
+            efit += x[iParam] * ShowerShapeV2(dx,dz) ;\r
+            iParam++ ;\r
+          }\r
+          \r
+          Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)\r
+          iParam = 0 ;\r
+          while(iParam < nPar ){\r
+            Double_t xpar = x[iParam] ;\r
+            Double_t zpar = x[iParam+1] ;\r
+            Double_t epar = x[iParam+2] ;\r
+            \r
+            Double_t dr = fgSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) );\r
+            Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
+            Double_t rp1  = TMath::Power(dr, fgSSPars[1]) ;\r
+            Double_t rp5  = TMath::Power(dr, fgSSPars[5]) ;\r
+            \r
+            Double_t deriv = -2 * TMath::Power(dr,fgSSPars[1]-2.) * fgSSPars[7] * fgSSPars[7] * \r
+            (fgSSPars[1] * ( 1/(fgSSPars[2]+fgSSPars[3]*rp1) + fgSSPars[4]/(1+fgSSPars[6]*rp5) ) - \r
+             (fgSSPars[1]*fgSSPars[3]*rp1/( (fgSSPars[2]+fgSSPars[3]*rp1)*(fgSSPars[2]+fgSSPars[3]*rp1) ) + \r
+              fgSSPars[4]*fgSSPars[5]*fgSSPars[6]*rp5/( (1+fgSSPars[6]*rp5)*(1+fgSSPars[6]*rp5) ) ) );\r
+            \r
+            //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) )\r
+            //                                                   - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;\r
+            \r
+            Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ;  // Derivative over x\r
+            iParam++ ;\r
+            Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ;  // Derivative over z\r
+            iParam++ ;\r
+            Grad[iParam] += shape ;                                  // Derivative over energy\r
+            iParam++ ;\r
+          }\r
+        }\r
+        efit = 0;\r
+        iparam = 0 ;\r
+        \r
+        while(iparam < nPar ){\r
+          Double_t xpar = x[iparam] ;\r
+          Double_t zpar = x[iparam+1] ;\r
+          Double_t epar = x[iparam+2] ;\r
+          iparam += 3 ;\r
+          efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
+        }\r
+        \r
+        fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;\r
+        // Here we assume, that sigma = sqrt(E) \r
+        } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!, nPar %d \n", nPar); // put nPar here to cheat coverity and rule checker\r
+      } // digit loop\r
+    } // recpoint, digits and geom not NULL\r
+  }// List is not NULL\r
+  \r
+}\r
+\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){\r
+  for(UInt_t i=0;i<7;++i)\r
+    fgSSPars[i]=pars[i];\r
+  if(pars[2]==0. && pars[3]==0.) fgSSPars[2]=1.;//to avoid dividing by 0\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::SetPar5(Double_t *pars){\r
+  for(UInt_t i=0;i<3;++i)\r
+    fgPar5[i]=pars[i];\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::SetPar6(Double_t *pars){\r
+  for(UInt_t i=0;i<3;++i)\r
+    fgPar6[i]=pars[i];\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::EvalPar5(Double_t phi){\r
+  //\r
+  //Evaluate the 5th parameter of the shower shape function\r
+  //phi in degrees range (-10,10)\r
+  //\r
+  //fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936;\r
+  fgSSPars[5] = fgPar5[0] + phi * fgPar5[1] + phi*phi * fgPar5[2];\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::EvalPar6(Double_t phi){\r
+  //\r
+  //Evaluate the 6th parameter of the shower shape function\r
+  //phi in degrees range (-10,10)\r
+  //\r
+  //fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361;\r
+  fgSSPars[6] = fgPar6[0] + phi * fgPar6[1] + phi*phi * fgPar6[2];\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, const AliEMCALGeometry *geom){\r
+  //\r
+  // calculate params p5 and p6 depending on the phi angle in global coordinate\r
+  // for the cell with given absId index\r
+  //\r
+  Double_t etaGlob = 0.;//eta in global c.s. - unused\r
+  Double_t phiGlob = 0.;//phi in global c.s. in radians\r
+  geom->EtaPhiFromIndex(absId, etaGlob, phiGlob);\r
+  phiGlob*=180./TMath::Pi();\r
+  phiGlob-=90.;\r
+  phiGlob-= (Double_t)((Int_t)geom->GetSuperModuleNumber(absId)/2 * 20);\r
+\r
+  EvalPar5(phiGlob);\r
+  EvalPar6(phiGlob);\r
+}\r
+\r