Fix double removal of rec points, existence of null pointers in rec points list;...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 23 Jul 2012 15:56:05 +0000 (15:56 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 23 Jul 2012 15:56:05 +0000 (15:56 +0000)
EMCAL/AliEMCALUnfolding.cxx

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