In case of unfolding, the splitted energy of the shared cells must not be lower than...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 9 Oct 2011 18:24:36 +0000 (18:24 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 9 Oct 2011 18:24:36 +0000 (18:24 +0000)
EMCAL/AliEMCALAfterBurnerUF.cxx
EMCAL/AliEMCALAfterBurnerUF.h
EMCAL/AliEMCALClusterizer.h
EMCAL/AliEMCALUnfolding.cxx
EMCAL/AliEMCALUnfolding.h

index c727fcb..d7f2b59 100644 (file)
@@ -80,6 +80,7 @@ AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF():
   fGeom(NULL),
   fLogWeight(4.5),      // correct?
   fECALocMaxCut(0.03),  // value suggested by Adam
+  fMinECut(0.01),
   fRecPoints(NULL),
   fDigitsArr(NULL),
   fClusterUnfolding(NULL)
@@ -90,10 +91,11 @@ AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF():
 }
 
 //------------------------------------------------------------------------
-AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ECALocMaxCut):
+AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ecaLocMaxCut, Float_t minECut):
   fGeom(NULL),
   fLogWeight(logWeight),
-  fECALocMaxCut(ECALocMaxCut),
+  fECALocMaxCut(ecaLocMaxCut),
+  fMinECut(minECut),
   fRecPoints(NULL),
   fDigitsArr(NULL),
   fClusterUnfolding(NULL)
@@ -135,7 +137,8 @@ void AliEMCALAfterBurnerUF::Init()
 
   fClusterUnfolding = new AliEMCALUnfolding(fGeom);
   fClusterUnfolding->SetECALocalMaxCut(fECALocMaxCut);
-
+  fClusterUnfolding->SetThreshold(fMinECut); 
+  
   // clusters --> recPoints, cells --> digits and back ;)
   fRecPoints = new TObjArray(100);
   fDigitsArr = new TClonesArray("AliEMCALDigit",1152);
index 9c456bf..ee0b938 100644 (file)
@@ -26,7 +26,7 @@ class AliEMCALAfterBurnerUF {
 
   public:
     AliEMCALAfterBurnerUF();
-    AliEMCALAfterBurnerUF(Float_t logWeight, Float_t locMaxCut);
+    AliEMCALAfterBurnerUF(Float_t logWeight, Float_t locMaxCut, Float_t minEcut);
     virtual ~AliEMCALAfterBurnerUF();
 
   private:
@@ -46,12 +46,13 @@ class AliEMCALAfterBurnerUF {
     AliEMCALGeometry  *fGeom;          // EMCAL geometry
     Float_t            fLogWeight;     // used in AliEMCALRecPoint::EvalGlobalPosition()
     Float_t            fECALocMaxCut;  // this amount of energy must distinguish a local maximum from its neighbours
+    Float_t            fMinECut;       // minimum energy of cell   
     TObjArray         *fRecPoints;     //! cluster <=> recPoint
     TClonesArray      *fDigitsArr;     //->   cell <=> digit
 
     AliEMCALUnfolding *fClusterUnfolding;  // unfolding class instance
 
-    ClassDef(AliEMCALAfterBurnerUF,1)
+    ClassDef(AliEMCALAfterBurnerUF,2)
 } ;
 
 #endif // AliEMCALAFTERBURNERUF_H
index 011b967..a772a18 100644 (file)
@@ -88,7 +88,8 @@ public:
   virtual void    SetPar5     (Int_t ipar, Double_t par)    { fPar5  [ipar] = par;             }
   virtual void    SetPar6     (Int_t ipar, Double_t par)    { fPar6  [ipar] = par;             }
   virtual void    InitClusterUnfolding()                    {
-    fClusterUnfolding=new AliEMCALUnfolding(fGeom,fECALocMaxCut,fSSPars,fPar5,fPar6);          }
+    fClusterUnfolding=new AliEMCALUnfolding(fGeom,fECALocMaxCut,fSSPars,fPar5,fPar6); 
+    fClusterUnfolding->SetThreshold(fMinECut);                                                 }
 
   //NxN (only used in NxN clusterizer)
   
index ebf6316..7be7bba 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 "TTree.h"\r
-//#include <TFile.h> \r
-//class TFolder;\r
-#include <TMath.h> \r
-#include <TMinuit.h>\r
-//#include <TTree.h> \r
-//class TSystem; \r
-//#include <TBenchmark.h>\r
-//#include <TBrowser.h>\r
-//#include <TROOT.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
-//#include "AliEMCALClusterizer.h"\r
-\r
-\r
-\r
-#include "AliLog.h"\r
-\r
-#include "AliCDBManager.h"\r
-//#include "AliCaloCalibPedestal.h"\r
-//#include "AliEMCALCalibData.h"\r
-class AliCDBStorage;\r
-#include "AliCDBEntry.h"\r
-\r
-Double_t AliEMCALUnfolding::fSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};\r
-Double_t AliEMCALUnfolding::fPar5[3]={12.31,-0.007381,-0.06936};\r
-Double_t AliEMCALUnfolding::fPar6[3]={0.05452,0.0001228,0.001361};\r
-\r
-ClassImp(AliEMCALUnfolding)\r
-  \r
-//____________________________________________________________________________\r
-AliEMCALUnfolding::AliEMCALUnfolding():\r
-  fNumberOfECAClusters(0),\r
-  fECALocMaxCut(0),\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
-  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
-  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++) fSSPars[i] = SSPars[i];\r
-  for (i = 0; i < 3; i++) {\r
-    fPar5[i] = Par5[i];\r
-    fPar6[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
-        //do we really need it?\r
-        //      TVector3 gpos;\r
-        //      Int_t absId = -1;\r
-        //      recPoint->GetGlobalPosition(gpos);\r
-        //      fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);\r
-        //      if(absId > nModulesToUnfold)\r
-        //        break ;\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
-  // 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
-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
-    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
-    Error("FindFit", "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 = fSSPars[7]*TMath::Sqrt(x*x+y*y);\r
-  Double_t rp1  = TMath::Power(r, fSSPars[1]) ;\r
-  Double_t rp5  = TMath::Power(r, fSSPars[5]) ;\r
-  Double_t shape = fSSPars[0]*TMath::Exp( -rp1 * (1. / (fSSPars[2] + fSSPars[3] * rp1) + fSSPars[4] / (1 + fSSPars[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 = fSSPars[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, fSSPars[1]) ;\r
-            Double_t rp5  = TMath::Power(dr, fSSPars[5]) ;\r
-            \r
-            Double_t deriv = -2 * TMath::Power(dr,fSSPars[1]-2.) * fSSPars[7] * fSSPars[7] * \r
-            (fSSPars[1] * ( 1/(fSSPars[2]+fSSPars[3]*rp1) + fSSPars[4]/(1+fSSPars[6]*rp5) ) - \r
-             (fSSPars[1]*fSSPars[3]*rp1/( (fSSPars[2]+fSSPars[3]*rp1)*(fSSPars[2]+fSSPars[3]*rp1) ) + \r
-              fSSPars[4]*fSSPars[5]*fSSPars[6]*rp5/( (1+fSSPars[6]*rp5)*(1+fSSPars[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!\n");\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
-    fSSPars[i]=pars[i];\r
-  if(pars[2]==0. && pars[3]==0.) fSSPars[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
-    fPar5[i]=pars[i];\r
-}\r
-\r
-//____________________________________________________________________________\r
-void AliEMCALUnfolding::SetPar6(Double_t *pars){\r
-  for(UInt_t i=0;i<3;++i)\r
-    fPar6[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
-  fSSPars[5] = fPar5[0] + phi * fPar5[1] + phi*phi * fPar5[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
-  fSSPars[6] = fPar6[0] + phi * fPar6[1] + phi*phi * fPar6[2];\r
-}\r
-\r
-//____________________________________________________________________________\r
-void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, 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 "TTree.h"\r\r
+//#include <TFile.h> \r\r
+//class TFolder;\r\r
+#include <TMath.h> \r\r
+#include <TMinuit.h>\r\r
+//#include <TTree.h> \r\r
+//class TSystem; \r\r
+//#include <TBenchmark.h>\r\r
+//#include <TBrowser.h>\r\r
+//#include <TROOT.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
+//#include "AliEMCALClusterizer.h"\r\r
+\r\r
+\r\r
+\r\r
+#include "AliLog.h"\r\r
+\r\r
+#include "AliCDBManager.h"\r\r
+//#include "AliCaloCalibPedestal.h"\r\r
+//#include "AliEMCALCalibData.h"\r\r
+class AliCDBStorage;\r\r
+#include "AliCDBEntry.h"\r\r
+\r\r
+Double_t AliEMCALUnfolding::fSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};\r\r
+Double_t AliEMCALUnfolding::fPar5[3]={12.31,-0.007381,-0.06936};\r\r
+Double_t AliEMCALUnfolding::fPar6[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
\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++) fSSPars[i] = SSPars[i];\r\r
+  for (i = 0; i < 3; i++) {\r\r
+    fPar5[i] = Par5[i];\r\r
+    fPar6[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) ;\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
+  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 numberofNotUnfolded = fNumberOfECAClusters ;\r\r
+    Int_t index ;\r\r
+    for(index = 0 ; index < numberofNotUnfolded ; index++){\r\r
+      AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;\r\r
+      if(recPoint){\r\r
+        //do we really need it?\r\r
+        //      TVector3 gpos;\r\r
+        //      Int_t absId = -1;\r\r
+        //      recPoint->GetGlobalPosition(gpos);\r\r
+        //      fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);\r\r
+        //      if(absId > nModulesToUnfold)\r\r
+        //        break ;\r\r
+        \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
+        \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
+            fRecPoints->Remove(recPoint);\r\r
+            fRecPoints->Compress() ;//is it really needed\r\r
+            index-- ;\r\r
+            fNumberOfECAClusters-- ;\r\r
+            numberofNotUnfolded-- ;\r\r
+          }\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 AliError("RecPoint NULL");\r\r
+    } // rec point loop\r\r
+  }\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
+  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
+  //**************************** 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
+  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
+  //**************************** 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 ;\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 orderis 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
+  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
+\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
+       if(efit[iDigit]==0) {//just for sure\r\r
+         correctedEnergyList[iparam/3+iDigit] = 0;\r\r
+         continue;\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+iDigit] = eDigit;\r\r
+\r\r
+      } else AliError("NULL digit");\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
+  Float_t maximumEne=0;\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
+    isAnyBelowThreshold=kFALSE;\r\r
+    maximumEne=0;\r\r
+    for(iparam = 0 ; iparam < nPar ; iparam+=3){\r\r
+\r\r
+      if(correctedEnergyList[iparam/3+iDigit] < fThreshold ) isAnyBelowThreshold = kTRUE;\r\r
+      if(correctedEnergyList[iparam/3+iDigit] > maximumEne) {\r\r
+       maximumEne = correctedEnergyList[iparam/3+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
+    if(maximumEne < fThreshold) {//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
+       maximumEne+=correctedEnergyList[iparam/3+iDigit];\r\r
+       correctedEnergyList[iparam/3+iDigit]=0;\r\r
+      }\r\r
+      correctedEnergyList[maximumIndex/3+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){//calculate energy sum\r\r
+      if(correctedEnergyList[iparam/3+iDigit] < fThreshold) energyFraction[iparam/3]=0;\r\r
+      else {\r\r
+       energyFraction[iparam/3]=1;\r\r
+       maximumEne+=correctedEnergyList[iparam/3+iDigit];\r\r
+      }\r\r
+    }//end of loop over clusters after unfolding\r\r
+    for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction\r\r
+      energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3+iDigit] / maximumEne;\r\r
+    }\r\r
+    for(iparam = 0 ; iparam < nPar ; iparam+=3){//add energy from cells below threshold to others\r\r
+      if(energyFraction[iparam/3]>0) continue;\r\r
+      else{\r\r
+       for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3){\r\r
+         correctedEnergyList[iparam2/3+iDigit] += (energyFraction[iparam2/3] * \r\r
+                                                   correctedEnergyList[iparam/3+iDigit]) ;\r\r
+       }//inner loop\r\r
+       correctedEnergyList[iparam/3+iDigit] = 0;\r\r
+      }\r\r
+    }\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
+  iparam = 0 ;\r\r
+  while(iparam < nPar ){\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(nSplittedClusters) ;\r\r
+      \r\r
+      for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r\r
+        digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r\r
+\r\r
+        if(digit && correctedEnergyList[iparam/3+iDigit]>0. ){\r\r
+          recPoint->AddDigit( *digit, correctedEnergyList[iparam/3+iDigit], 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
+    iparam += 3 ;\r\r
+  }//while\r\r
+  \r\r
+  delete[] fitparameters ;\r\r
+  delete[] efit ;\r\r
+  delete[] correctedEnergyList ;\r\r
+\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
+\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
+    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", iphi ) ;\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", ieta) ;\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", energy) ;\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
+    Error("FindFit", "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 of parameter index\r\r
+    fitparameters[index] = val ;\r\r
+  }\r\r
+\r\r
+  toMinuit->Clear();\r\r
+  delete toMinuit ;\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 = fSSPars[7]*TMath::Sqrt(x*x+y*y);\r\r
+  Double_t rp1  = TMath::Power(r, fSSPars[1]) ;\r\r
+  Double_t rp5  = TMath::Power(r, fSSPars[5]) ;\r\r
+  Double_t shape = fSSPars[0]*TMath::Exp( -rp1 * (1. / (fSSPars[2] + fSSPars[3] * rp1) + fSSPars[4] / (1 + fSSPars[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 = fSSPars[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, fSSPars[1]) ;\r\r
+            Double_t rp5  = TMath::Power(dr, fSSPars[5]) ;\r\r
+            \r\r
+            Double_t deriv = -2 * TMath::Power(dr,fSSPars[1]-2.) * fSSPars[7] * fSSPars[7] * \r\r
+            (fSSPars[1] * ( 1/(fSSPars[2]+fSSPars[3]*rp1) + fSSPars[4]/(1+fSSPars[6]*rp5) ) - \r\r
+             (fSSPars[1]*fSSPars[3]*rp1/( (fSSPars[2]+fSSPars[3]*rp1)*(fSSPars[2]+fSSPars[3]*rp1) ) + \r\r
+              fSSPars[4]*fSSPars[5]*fSSPars[6]*rp5/( (1+fSSPars[6]*rp5)*(1+fSSPars[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!\n");\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
+    fSSPars[i]=pars[i];\r\r
+  if(pars[2]==0. && pars[3]==0.) fSSPars[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
+    fPar5[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
+    fPar6[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
+  fSSPars[5] = fPar5[0] + phi * fPar5[1] + phi*phi * fPar5[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
+  fSSPars[6] = fPar6[0] + phi * fPar6[1] + phi*phi * fPar6[2];\r\r
+}\r\r
+\r\r
+//____________________________________________________________________________\r\r
+void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, 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
index aada310..8f69142 100644 (file)
@@ -1,91 +1,97 @@
-#ifndef ALIEMCALUNFOLDING_H\r
-#define ALIEMCALUNFOLDING_H\r
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- * See cxx source for full Copyright notice                               */\r
-     \r
-//_________________________________________________________________________\r
-//  Base class for the cluster unfolding algorithm \r
-//*-- Author: Adam Matyja (SUBATECH)\r
-\r
-// --- ROOT system ---\r
-#include "AliLog.h"\r
-#include "TObject.h" \r
-//class TTree;\r
-\r
-// --- Standard library ---\r
-\r
-// --- AliRoot header files ---\r
-class AliEMCALGeometry ;\r
-//class AliEMCALCalibData ;\r
-//class AliCaloCalibPedestal ;\r
-class AliEMCALRecPoint ; \r
-class AliEMCALDigit ;\r
-\r
-\r
-class AliEMCALUnfolding : public TObject {\r
-\r
-public:\r
-\r
-  AliEMCALUnfolding() ;        // default ctor\r
-  virtual ~AliEMCALUnfolding() ; // dtorEM\r
-  AliEMCALUnfolding(AliEMCALGeometry* geometry);// constructor\r
-  AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMaxCut,Double_t *SSPars,Double_t *Par5,Double_t *Par6);// constructor\r
-\r
-  virtual void Init() ;\r
-  virtual void SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,TClonesArray *digitsArr);\r
-\r
-  //setters and getters\r
-  virtual void SetNumberOfECAClusters(Int_t n) { fNumberOfECAClusters = n; }\r
-  virtual Int_t GetNumberOfECAClusters() const { return fNumberOfECAClusters; }\r
-  virtual void SetRecPoints(TObjArray *rec) { fRecPoints = rec; }\r
-  virtual TObjArray * GetRecPoints() const { return fRecPoints; }\r
-  virtual void SetDigitsArr(TClonesArray *digit) { fDigitsArr = digit; }\r
-  virtual TClonesArray * GetDigitsArr() const { return fDigitsArr; }\r
-  virtual void SetECALocalMaxCut(Float_t cut) { fECALocMaxCut = cut ; }\r
-  virtual Float_t GetECALocalMaxCut() const { return fECALocMaxCut; }\r
-\r
-  //unfolding main methods\r
-  virtual void   MakeUnfolding();\r
-  static Double_t ShowerShapeV2(Double_t x, Double_t y) ; // Shape of EM shower used in unfolding; \r
-                                              //class member function (not object member function)\r
-  static void UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)  ;\r
-                                            // Chi^2 of the fit. Should be static to be passes to MINUIT\r
-  virtual void SetShowerShapeParams(Double_t *pars) ;\r
-  virtual Double_t* GetShowerShapeParams() const { return fSSPars ; }\r
-  virtual void SetPar5(Double_t *pars) ;\r
-  virtual Double_t* GetPar5() const { return fPar5 ; }\r
-  virtual void SetPar6(Double_t *pars) ;\r
-  virtual Double_t* GetPar6() const { return fPar6 ; }\r
-\r
-protected:\r
-  Int_t   fNumberOfECAClusters ;     // number of clusters found in EC section\r
-  Float_t fECALocMaxCut ;            // minimum energy difference to distinguish local maxima in a cluster\r
-  AliEMCALGeometry     * fGeom;       //! pointer to geometry for utilities\r
-  TObjArray    *fRecPoints; // Array with EMCAL clusters\r
-  TClonesArray *fDigitsArr; // Array with EMCAL digits\r
-\r
-private:\r
-  AliEMCALUnfolding(const AliEMCALUnfolding &); //copy ctor\r
-  AliEMCALUnfolding & operator = (const AliEMCALUnfolding &);\r
-  \r
-  Bool_t         UnfoldClusterV2(AliEMCALRecPoint * iniEmc, Int_t Nmax, \r
-                                AliEMCALDigit ** maxAt,\r
-                                Float_t * maxAtEnergy ); //Unfolds cluster using TMinuit package\r
-  Bool_t  FindFitV2(AliEMCALRecPoint * emcRP, AliEMCALDigit ** MaxAt, const Float_t * maxAtEnergy, \r
-                   Int_t NPar, Float_t * FitParametres) const; //Used in UnfoldClusters, calls TMinuit\r
-\r
-  static Double_t fSSPars[8];//! Unfolding shower shape parameters\r
-  // function:\r
-  // f(r)=exp(-(p0*r)^p1 * (1/(p2+p3*(p0*r)^p1)+p4/(1+p6*(p0*r)^p5) ) )\r
-  // p0,p1,p2,p3,p4 are fixed\r
-  // params p5 and p6 are phi-dependent and set in ShowerShapeV2\r
-  static Double_t fPar5[3];//! UF SSPar nr 5 = p0 + phi*p1 + phi^2 *p2\r
-  static Double_t fPar6[3];//! UF SSPar nr 6 = p0 + phi*p1 + phi^2 *p2\r
-  static void EvalPar5(Double_t phi);\r
-  static void EvalPar6(Double_t phi);\r
-  static void EvalParsPhiDependence(Int_t absId, AliEMCALGeometry *geom);\r
-\r
-  ClassDef(AliEMCALUnfolding,1)  // Unfolding algorithm class \r
-} ;\r
-\r
-#endif // AliEMCALUNFOLDING_H\r
+#ifndef ALIEMCALUNFOLDING_H\r\r
+#define ALIEMCALUNFOLDING_H\r\r
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r\r
+ * See cxx source for full Copyright notice                               */\r\r
+     \r\r
+//_________________________________________________________________________\r\r
+//  Base class for the cluster unfolding algorithm \r\r
+//*-- Author: Adam Matyja (SUBATECH)\r\r
+\r\r
+// --- ROOT system ---\r\r
+#include "AliLog.h"\r\r
+#include "TObject.h" \r\r
+//class TTree;\r\r
+\r\r
+// --- Standard library ---\r\r
+\r\r
+// --- AliRoot header files ---\r\r
+class AliEMCALGeometry ;\r\r
+//class AliEMCALCalibData ;\r\r
+//class AliCaloCalibPedestal ;\r\r
+class AliEMCALRecPoint ; \r\r
+class AliEMCALDigit ;\r\r
+\r\r
+\r\r
+class AliEMCALUnfolding : public TObject {\r\r
+\r\r
+public:\r\r
+\r\r
+  AliEMCALUnfolding() ;        // default ctor\r\r
+  virtual ~AliEMCALUnfolding() ; // dtorEM\r\r
+  AliEMCALUnfolding(AliEMCALGeometry* geometry);// constructor\r\r
+  AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMaxCut,Double_t *SSPars,Double_t *Par5,Double_t *Par6);// constructor\r\r
+\r\r
+  virtual void Init() ;\r\r
+  virtual void SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,TClonesArray *digitsArr);\r\r
+\r\r
+  //setters and getters\r\r
+  virtual void SetNumberOfECAClusters(Int_t n) { fNumberOfECAClusters = n; }\r\r
+  virtual Int_t GetNumberOfECAClusters() const { return fNumberOfECAClusters; }\r\r
+  virtual void SetRecPoints(TObjArray *rec) { fRecPoints = rec; }\r\r
+  virtual TObjArray * GetRecPoints() const { return fRecPoints; }\r\r
+  virtual void SetDigitsArr(TClonesArray *digit) { fDigitsArr = digit; }\r\r
+  virtual TClonesArray * GetDigitsArr() const { return fDigitsArr; }\r\r
+  virtual void SetECALocalMaxCut(Float_t cut) { fECALocMaxCut = cut ; }\r\r
+  virtual Float_t GetECALocalMaxCut() const { return fECALocMaxCut; }\r\r
+  virtual void SetThreshold(Float_t energy) { fThreshold = energy; }\r\r
+  virtual Float_t GetThreshold() const { return fThreshold; }\r\r
+\r\r
+  //unfolding main methods\r\r
+  virtual void   MakeUnfolding();\r\r
+  static Double_t ShowerShapeV2(Double_t x, Double_t y) ; // Shape of EM shower used in unfolding; \r\r
+                                              //class member function (not object member function)\r\r
+  static void UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)  ;\r\r
+                                            // Chi^2 of the fit. Should be static to be passes to MINUIT\r\r
+  virtual void SetShowerShapeParams(Double_t *pars) ;\r\r
+  virtual Double_t* GetShowerShapeParams() const { return fSSPars ; }\r\r
+  virtual void SetPar5(Double_t *pars) ;\r\r
+  virtual Double_t* GetPar5() const { return fPar5 ; }\r\r
+  virtual void SetPar6(Double_t *pars) ;\r\r
+  virtual Double_t* GetPar6() const { return fPar6 ; }\r\r
+\r\r
+protected:\r\r
+  Int_t   fNumberOfECAClusters ;     // number of clusters found in EC section\r\r
+  Float_t fECALocMaxCut ;            // minimum energy difference to distinguish local maxima in a cluster\r\r
+  Float_t fThreshold ; //minimum energy for cell to be joined to a cluster\r\r
+  AliEMCALGeometry     * fGeom;       //! pointer to geometry for utilities\r\r
+  TObjArray    *fRecPoints; // Array with EMCAL clusters\r\r
+  TClonesArray *fDigitsArr; // Array with EMCAL digits\r\r
+\r\r
+private:\r\r
+  AliEMCALUnfolding(const AliEMCALUnfolding &); //copy ctor\r\r
+  AliEMCALUnfolding & operator = (const AliEMCALUnfolding &);\r\r
+  \r\r
+  Bool_t         UnfoldClusterV2(AliEMCALRecPoint * iniEmc, Int_t Nmax, \r\r
+                                AliEMCALDigit ** maxAt,\r\r
+                                Float_t * maxAtEnergy ); //Unfolds cluster using TMinuit package\r\r
+  Bool_t         UnfoldClusterV2old(AliEMCALRecPoint * iniEmc, Int_t Nmax, \r\r
+                                   AliEMCALDigit ** maxAt,\r\r
+                                   Float_t * maxAtEnergy ); //Unfolds cluster using TMinuit package\r\r
+  Bool_t  FindFitV2(AliEMCALRecPoint * emcRP, AliEMCALDigit ** MaxAt, const Float_t * maxAtEnergy, \r\r
+                   Int_t NPar, Float_t * FitParametres) const; //Used in UnfoldClusters, calls TMinuit\r\r
+\r\r
+  static Double_t fSSPars[8];//! Unfolding shower shape parameters\r\r
+  // function:\r\r
+  // f(r)=exp(-(p0*r)^p1 * (1/(p2+p3*(p0*r)^p1)+p4/(1+p6*(p0*r)^p5) ) )\r\r
+  // p0,p1,p2,p3,p4 are fixed\r\r
+  // params p5 and p6 are phi-dependent and set in ShowerShapeV2\r\r
+  static Double_t fPar5[3];//! UF SSPar nr 5 = p0 + phi*p1 + phi^2 *p2\r\r
+  static Double_t fPar6[3];//! UF SSPar nr 6 = p0 + phi*p1 + phi^2 *p2\r\r
+  static void EvalPar5(Double_t phi);\r\r
+  static void EvalPar6(Double_t phi);\r\r
+  static void EvalParsPhiDependence(Int_t absId, AliEMCALGeometry *geom);\r\r
+\r\r
+  ClassDef(AliEMCALUnfolding,2)  // Unfolding algorithm class \r\r
+} ;\r\r
+\r\r
+#endif // AliEMCALUNFOLDING_H\r\r