change definition of TRU fitter, remove unneeded class dependencies
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALUnfolding.cxx
index f59c609..c9d603f 100644 (file)
@@ -47,9 +47,9 @@
 class AliCDBStorage;
 #include "AliCDBEntry.h"
 
-Double_t AliEMCALUnfolding::fSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};
-Double_t AliEMCALUnfolding::fPar5[3]={12.31,-0.007381,-0.06936};
-Double_t AliEMCALUnfolding::fPar6[3]={0.05452,0.0001228,0.001361};
+Double_t AliEMCALUnfolding::fgSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};
+Double_t AliEMCALUnfolding::fgPar5[3]={12.31,-0.007381,-0.06936};
+Double_t AliEMCALUnfolding::fgPar6[3]={0.05452,0.0001228,0.001361};
 
 ClassImp(AliEMCALUnfolding)
   
@@ -58,12 +58,12 @@ AliEMCALUnfolding::AliEMCALUnfolding():
   fNumberOfECAClusters(0),
   fECALocMaxCut(0),
   fThreshold(0.01),//10 MeV
+  fRejectBelowThreshold(0),//split
   fGeom(NULL),
   fRecPoints(NULL),
   fDigitsArr(NULL)
 {
   // ctor with the indication of the file where header Tree and digits Tree are stored
   Init() ;
 }
 
@@ -72,6 +72,7 @@ AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry):
   fNumberOfECAClusters(0),
   fECALocMaxCut(0),
   fThreshold(0.01),//10 MeV
+  fRejectBelowThreshold(0),//split
   fGeom(geometry),
   fRecPoints(NULL),
   fDigitsArr(NULL)
@@ -83,7 +84,7 @@ AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry):
   {
     AliFatal("AliEMCALUnfolding: Geometry not initialized.");
   }
-  
+
 }
 
 //____________________________________________________________________________
@@ -91,6 +92,7 @@ AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMa
   fNumberOfECAClusters(0),
   fECALocMaxCut(ECALocMaxCut),
   fThreshold(0.01),//10 MeV
+  fRejectBelowThreshold(0),//split
   fGeom(geometry),
   fRecPoints(NULL),
   fDigitsArr(NULL)
@@ -103,10 +105,10 @@ AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMa
     AliFatal("AliEMCALUnfolding: Geometry not initialized.");
   }
   Int_t i=0;
-  for (i = 0; i < 8; i++) fSSPars[i] = SSPars[i];
+  for (i = 0; i < 8; i++) fgSSPars[i] = SSPars[i];
   for (i = 0; i < 3; i++) {
-    fPar5[i] = Par5[i];
-    fPar6[i] = Par6[i];
+    fgPar5[i] = Par5[i];
+    fgPar6[i] = Par6[i];
   }
 
 }
@@ -129,7 +131,8 @@ void AliEMCALUnfolding::Init()
   AliDebug(1,Form("geom %p",fGeom));
   
   if(!gMinuit) 
-    gMinuit = new TMinuit(100) ;
+    //    gMinuit = new TMinuit(100) ;//the same is in FindFitV2
+    gMinuit = new TMinuit(30) ;//the same is in FindFitV2
   
 }
 
@@ -144,6 +147,7 @@ void AliEMCALUnfolding::SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,
 {
   //
   //Set input for unfolding purposes
+  //
   SetNumberOfECAClusters(numberOfECAClusters);
   SetRecPoints(recPoints);
   SetDigitsArr(digitsArr);
@@ -155,54 +159,84 @@ void AliEMCALUnfolding::MakeUnfolding()
   // Unfolds clusters using the shape of an ElectroMagnetic shower
   // Performs unfolding of all clusters
   
+  AliDebug(4,Form(" V1: total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
   if(fNumberOfECAClusters > 0){
     if (fGeom==0)
       AliFatal("Did not get geometry from EMCALLoader") ;
     //Int_t nModulesToUnfold = fGeom->GetNCells();
     
-    Int_t numberofNotUnfolded = fNumberOfECAClusters ;
+    Int_t numberOfClustersToUnfold=fNumberOfECAClusters;
+    //we unfold only clusters present in the array untill now
+    //fNumberOfECAClusters may change due to unfilded clusters
+    //so 0 to numberOfClustersToUnfold-1: clusters before unfolding
+    //numberOfClustersToUnfold to the end: new clusters from unfolding
+    //of course numberOfClustersToUnfold also is decreased but we don't loop over clusters added in UF 
     Int_t index ;
-    for(index = 0 ; index < numberofNotUnfolded ; index++){
+    for(index = 0 ; index < numberOfClustersToUnfold ; index++){
       AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
       if(recPoint){
-        
         Int_t nMultipl = recPoint->GetMultiplicity() ;
         AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
         Float_t * maxAtEnergy = new Float_t[nMultipl] ;
         Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
-        
         if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
+         AliDebug(4,Form("  *** V1+UNFOLD *** Cluster index before UF %d",fNumberOfECAClusters));
           if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){
+           //if unfolding correct remove old recPoint
             fRecPoints->Remove(recPoint);
             fRecPoints->Compress() ;//is it really needed
             index-- ;
             fNumberOfECAClusters-- ;
-            numberofNotUnfolded-- ;
+           numberOfClustersToUnfold--;
           }
-        }
-        else{
+         AliDebug(4,Form("  Cluster index after UF %d",fNumberOfECAClusters));
+       } else{
           recPoint->SetNExMax(1) ; //Only one local maximum
         }
         
         delete[] maxAt ;
         delete[] maxAtEnergy ;
-      } else AliError("RecPoint NULL");
+      } else {
+       //AliError("RecPoint NULL"); //end of check if recPoint exist
+       Error("MakeUnfolding", "RecPoint NULL, index = %d, fNumberOfECAClusters = %d, numberOfClustersToUnfold = %d",index,fNumberOfECAClusters,numberOfClustersToUnfold) ;
+      }
     } // rec point loop
-  }
+  }//end of check fNumberOfECAClusters
   // End of Unfolding of clusters
+
+  AliDebug(4,Form(" V1+UNFOLD: total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
+//  for(Int_t i=0;i<fNumberOfECAClusters;i++){
+//    AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(i));
+//    Int_t nMultipl = recPoint->GetMultiplicity() ;
+//    Double_t energy=recPoint->GetEnergy();
+//    Int_t absIdMaxDigit=recPoint->GetAbsIdMaxDigit();
+//    Int_t sm=recPoint->GetSuperModuleNumber();
+//    Double_t pointEne=recPoint->GetPointEnergy();
+//    Float_t maxEne=recPoint->GetMaximalEnergy();
+//    Int_t maxEneInd=recPoint->GetMaximalEnergyIndex();
+//    printf("  cluster %d,ncells %d,ene %f,absIdMaxCell %d,sm %d,pointEne %f,maxEne %f,maxEneInd %d\n",i,nMultipl,energy,absIdMaxDigit,sm,pointEne,maxEne,maxEneInd);
+//  }
+
 }
 
 //____________________________________________________________________________
-Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, 
+Int_t AliEMCALUnfolding::UnfoldOneCluster(AliEMCALRecPoint * iniTower, 
                                          Int_t nMax, 
                                          AliEMCALDigit ** maxAt, 
-                                         Float_t * maxAtEnergy)
+                                         Float_t * maxAtEnergy,
+                                         TObjArray *list)
 {
-  // Extended to whole EMCAL 
-
+  // Input one cluster
+  // Output list of clusters
+  // returns number of clusters
+  // if fit failed or unfolding is not applicable returns 0 and empty list
+  
   //**************************** part 1 *******************************************
   // Performs the unfolding of a cluster with nMax overlapping showers 
   
+  //cout<<"unfolding check here part 1"<<endl;
+  AliDebug(5,Form("  Original cluster E %f, nMax = %d",iniTower->GetEnergy(),nMax ));
+
   Int_t nPar = 3 * nMax ;
   Float_t * fitparameters = new Float_t[nPar] ;
   
@@ -210,13 +244,24 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
     AliFatal("Did not get geometry from EMCALLoader") ;
   
   Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
-  if( !rv ) {
+  if( !rv ) 
+  {
     // Fit failed, return (and remove cluster? - why? I leave the cluster)
     iniTower->SetNExMax(-1) ;
     delete[] fitparameters ;
-    return kFALSE;
+    return 0;//changed here
   }
   
+  //speed up solution for clusters with 2 maxima where one maximum is below threshold fThreshold
+  if(nMax==2){
+    if(fitparameters[2]<fThreshold || fitparameters[5]<fThreshold){
+      AliDebug(1,"One of fitted energy below threshold");
+      iniTower->SetNExMax(1) ;
+      delete[] fitparameters ;
+      return 0;//changed here
+    }
+  }
+
   //**************************** part 2 *******************************************
   // create unfolded rec points and fill them with new energy lists
   // First calculate energy deposited in each sell in accordance with
@@ -224,6 +269,7 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
   // and later correct this number in acordance with actual energy
   // deposition
   
+  //  cout<<"unfolding check here part 2"<<endl;
   Int_t nDigits = iniTower->GetMultiplicity() ;
   Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells
   Float_t xpar=0.,zpar=0.,epar=0.  ;//center of gravity in cell units
@@ -241,9 +287,11 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
   Int_t iparam = 0 ;
   Int_t iDigit = 0 ;
   
-  for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
+  for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
+  {
     digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
-    if(digit){
+    if(digit)
+    {
       fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
       fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
                                          iIphi, iIeta,iphi,ieta);
@@ -251,15 +299,17 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
       
       efit[iDigit] = 0.;
       iparam = 0;
-      while(iparam < nPar ){
+      while(iparam < nPar )
+      {
         xpar = fitparameters[iparam] ;
         zpar = fitparameters[iparam+1] ;
         epar = fitparameters[iparam+2] ;
-        iparam += 3 ;
-        
+
         efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
+        iparam += 3 ;
       }
-    } else AliError("Digit NULL part 2!");
+
+    } else AliDebug(1,"Digit NULL part 2!");
     
   }//digit loop
   
@@ -269,168 +319,294 @@ Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower,
   // to its contribution to efit
   
   Float_t * energiesList = iniTower->GetEnergiesList() ;
-  Float_t ratio = 0 ;
+  Float_t ratio = 0. ;
   Float_t eDigit = 0. ;
   Int_t nSplittedClusters=(Int_t)nPar/3;
   
   Float_t * correctedEnergyList = new Float_t[nDigits*nSplittedClusters];
   //above - temporary table with energies after unfolding.
-  //the orderis following: 
+  //the order is following: 
   //first cluster <first cell - last cell>, 
   //second cluster <first cell - last cell>, etc.
-
+  
   //**************************** sub-part 3.1 *************************************
-  //here we check if energy of the cell in the cluster after unfolding is above threshold. 
   //If not the energy from a given cell in the cluster is divided in correct proportions 
   //in accordance to the other clusters and added to them and set to 0.
+  
+  //  cout<<"unfolding check here part 3.1"<<endl;
 
   iparam = 0 ;
-  while(iparam < nPar ){
+  while(iparam < nPar )
+  {
     xpar = fitparameters[iparam] ;
     zpar = fitparameters[iparam+1] ;
     epar = fitparameters[iparam+2] ;
-
-
-    for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
+    
+    for(iDigit = 0 ; iDigit < nDigits ; iDigit ++)
+    {
       digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
-      if(digit){
-       fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
-       fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
-                                          iIphi, iIeta,iphi,ieta);
-       EvalParsPhiDependence(digit->GetId(),fGeom);
-       if(efit[iDigit]==0) {//just for sure
-         correctedEnergyList[iparam/3+iDigit] = 0;
-         continue;
-       }
-       ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;
-       eDigit = energiesList[iDigit] * ratio ;
-
-       //add energy to temporary matrix
-       correctedEnergyList[iparam/3+iDigit] = eDigit;
-
-      } else AliError("NULL digit part 3");
+      if(digit)
+      {
+        fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
+        fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
+                                           iIphi, iIeta,iphi,ieta);
+        
+        EvalParsPhiDependence(digit->GetId(),fGeom);
+       
+        if(efit[iDigit]==0) 
+        {//just for sure
+          correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;//correction here
+          continue;
+        }
+        
+        ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;
+        eDigit = energiesList[iDigit] * ratio ;
+        
+        //add energy to temporary matrix
+        correctedEnergyList[iparam/3*nDigits+iDigit] = eDigit;
+        
+      } else AliDebug(1,"NULL digit part 3");
     }//digit loop 
     iparam += 3 ;
   }//while
-
+  
   //**************************** sub-part 3.2 *************************************
+  //here we check if energy of the cell in the cluster after unfolding is above threshold. 
   //here we correct energy for each cell and cluster
-  Float_t maximumEne=0;
-  Int_t maximumIndex=0;
-  Bool_t isAnyBelowThreshold=kFALSE;
-  //  Float_t Threshold=0.01;
-  Float_t * energyFraction = new Float_t[nSplittedClusters];
-  Int_t iparam2 = 0 ;
-  for(iDigit = 0 ; iDigit < nDigits ; iDigit++){
-    isAnyBelowThreshold=kFALSE;
-    maximumEne=0;
-    for(iparam = 0 ; iparam < nPar ; iparam+=3){
-
-      if(correctedEnergyList[iparam/3+iDigit] < fThreshold ) isAnyBelowThreshold = kTRUE;
-      if(correctedEnergyList[iparam/3+iDigit] > maximumEne) {
-       maximumEne = correctedEnergyList[iparam/3+iDigit];
-       maximumIndex = iparam;
-      }
-    }//end of loop over clusters after unfolding
-
-    if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold 
-    if(maximumEne < fThreshold) {//add all cluster cells and put energy into max index, other set to 0
-      maximumEne=0.;
-      for(iparam = 0 ; iparam < nPar ; iparam+=3){
-       maximumEne+=correctedEnergyList[iparam/3+iDigit];
-       correctedEnergyList[iparam/3+iDigit]=0;
-      }
-      correctedEnergyList[maximumIndex/3+iDigit]=maximumEne;
-      continue;
-    }//end if
-
-    //divide energy of cell below threshold in the correct proportion and add to other cells
-    maximumEne=0;//not used any more so use it for the energy sum
-    for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate energy sum
-      if(correctedEnergyList[iparam/3+iDigit] < fThreshold) energyFraction[iparam/3]=0;
-      else {
-       energyFraction[iparam/3]=1;
-       maximumEne+=correctedEnergyList[iparam/3+iDigit];
-      }
-    }//end of loop over clusters after unfolding
-    if(maximumEne>0){
-      for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction
-       energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3+iDigit] / maximumEne;
-      }
-      
-      for(iparam = 0 ; iparam < nPar ; iparam+=3){//add energy from cells below threshold to others
-       if(energyFraction[iparam/3]>0) continue;
-       else{
-         for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3){
-           correctedEnergyList[iparam2/3+iDigit] += (energyFraction[iparam2/3] * 
-                                                     correctedEnergyList[iparam/3+iDigit]) ;
-         }//inner loop
-         correctedEnergyList[iparam/3+iDigit] = 0;
+  //  cout<<"unfolding check here part 3.2"<<endl;
+
+
+  //here we have 3 possibilities
+  //when after UF cell energy in cluster is below threshold:
+  //1 - keep it associated to cluster - equivalent of threshold=0
+  //2 - default - split (or add) energy of that cell into that cell in the other cluster(s)
+  //3 - reject that cell from cluster - fraction energy in cell=0 - breaks energy conservation
+  //Bool_t rejectBelowThreshold=kTRUE;//default option = 2 - split = kFALSE
+
+  if(fThreshold > 0){//option 2 or 3
+    if(fRejectBelowThreshold){//option 3
+      for(iDigit = 0 ; iDigit < nDigits ; iDigit++){//digit loop
+       for(iparam = 0 ; iparam < nPar ; iparam+=3){//param0 loop = energy loop
+         if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold ) correctedEnergyList[iparam/3*nDigits+iDigit]=0.;
        }
       }
-    }
-    else{
-      //digit energy to be set to 0
-      for(iparam = 0 ; iparam < nPar ; iparam+=3){
-       correctedEnergyList[iparam/3+iDigit] = 0;
-      }
-    }//new adam correction for is energy>0
-
-  }//end of loop over digits
-  delete[] energyFraction;
+    }else{//option 2
+      Float_t maximumEne=0.;
+      Int_t maximumIndex=0;
+      Bool_t isAnyBelowThreshold=kFALSE;
+      //  Float_t Threshold=0.01;
+      Float_t * energyFraction = new Float_t[nSplittedClusters];
+      Int_t iparam2 = 0 ;
+      for(iDigit = 0 ; iDigit < nDigits ; iDigit++){
+       isAnyBelowThreshold=kFALSE;
+       maximumEne=0.;
+       for(iparam = 0 ; iparam < nPar ; iparam+=3){
+         if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold ) isAnyBelowThreshold = kTRUE;
+         if(correctedEnergyList[iparam/3*nDigits+iDigit] > maximumEne) 
+           {
+             maximumEne = correctedEnergyList[iparam/3*nDigits+iDigit];
+             maximumIndex = iparam;
+           }
+       }//end of loop over clusters after unfolding
+       
+       if(!isAnyBelowThreshold) continue; //no cluster-cell below threshold 
+       
+       if(maximumEne < fThreshold) 
+         {//add all cluster cells and put energy into max index, other set to 0
+           maximumEne=0.;
+           for(iparam = 0 ; iparam < nPar ; iparam+=3)
+             {
+               maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];
+               correctedEnergyList[iparam/3*nDigits+iDigit]=0;
+             }
+           correctedEnergyList[maximumIndex/3*nDigits+iDigit]=maximumEne;
+           continue;
+         }//end if
+       
+       //divide energy of cell below threshold in the correct proportion and add to other cells
+       maximumEne=0.;//not used any more so use it for the energy sum 
+       for(iparam = 0 ; iparam < nPar ; iparam+=3)
+         {//calculate energy sum
+           if(correctedEnergyList[iparam/3*nDigits+iDigit] < fThreshold) energyFraction[iparam/3]=0;
+           else 
+             {
+               energyFraction[iparam/3]=1.;
+               maximumEne+=correctedEnergyList[iparam/3*nDigits+iDigit];
+             }
+         }//end of loop over clusters after unfolding
+       if(maximumEne>0.) {
+         for(iparam = 0 ; iparam < nPar ; iparam+=3){//calculate fraction
+           energyFraction[iparam/3] = energyFraction[iparam/3] * correctedEnergyList[iparam/3*nDigits+iDigit] / maximumEne;
+         }
+         
+         for(iparam = 0 ; iparam < nPar ; iparam+=3)
+           {//add energy from cells below threshold to others
+             if(energyFraction[iparam/3]>0.) continue;
+             else
+               {
+                 for(iparam2 = 0 ; iparam2 < nPar ; iparam2+=3)
+                   {
+                     correctedEnergyList[iparam2/3*nDigits+iDigit] += (energyFraction[iparam2/3] * 
+                                                                       correctedEnergyList[iparam/3*nDigits+iDigit]) ;
+                   }//inner loop
+                 correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;
+               }
+           }
+       } else {
+         //digit energy to be set to 0
+         for(iparam = 0 ; iparam < nPar ; iparam+=3)
+           {
+             correctedEnergyList[iparam/3*nDigits+iDigit] = 0.;
+           }
+       }//correction for: is energy>0
+       
+      }//end of loop over digits
+      delete[] energyFraction;
+      
+    }//end of option 2 or 3 
+  } else {//option 1
+    //do nothing
+  }
 
+  
   //**************************** sub-part 3.3 *************************************
   //here we add digits to recpoints with corrected energy
+  //  cout<<"unfolding check here part 3.3"<<endl;
+
+  Int_t newClusterIndex=0;
   iparam = 0 ;
-  while(iparam < nPar ){
+  while(iparam < nPar )
+  {
     AliEMCALRecPoint * recPoint = 0 ;
     
-    if(fNumberOfECAClusters >= fRecPoints->GetSize())
-      fRecPoints->Expand(2*fNumberOfECAClusters) ;
+    if(nSplittedClusters >= list->GetSize())
+      list->Expand(nSplittedClusters);
     
     //add recpoint
-    (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
-    recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
+    (*list)[newClusterIndex] = new AliEMCALRecPoint("") ;
+    recPoint = dynamic_cast<AliEMCALRecPoint *>( list->At(newClusterIndex) ) ;
     
-    if(recPoint){
-      
-      fNumberOfECAClusters++ ;
-      recPoint->SetNExMax(nSplittedClusters) ;
+    if(recPoint){//recPoint present -> good
+      recPoint->SetNExMax(nSplittedClusters) ;//can be wrong number, to be corrected in outer method
       
-      for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
+      for(iDigit = 0 ; iDigit < nDigits ; iDigit ++) {
         digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
-
-        if(digit && correctedEnergyList[iparam/3+iDigit]>0. ){
-          recPoint->AddDigit( *digit, correctedEnergyList[iparam/3+iDigit], kFALSE ) ; //FIXME, need to study the shared case
+        if(digit && correctedEnergyList[iparam/3*nDigits+iDigit]>0. ){
+          //if(correctedEnergyList[iparam/3*nDigits+iDigit]<fThreshold) printf("Final E cell %f < %f\n",correctedEnergyList[iparam/3*nDigits+iDigit],fThreshold);
+          recPoint->AddDigit( *digit, correctedEnergyList[iparam/3*nDigits+iDigit], kFALSE ) ; //FIXME, need to study the shared case
         } else {
-         AliError("NULL digit part3.3 or energy=0");
-         //cout<<"nDigits "<<nDigits<<" iParam/3 "<<iparam/3<< endl;
-       }
-      }//digit loop 
-    } else AliError("NULL RecPoint");
-    //protection from recpoint with no digits
-    //cout<<"multi rec "<<recPoint->GetMultiplicity()<<endl;
-    if(recPoint->GetMultiplicity()==0){
-      delete (*fRecPoints)[fNumberOfECAClusters];
-      //cout<<"size fRecPoints before "<<fRecPoints->GetSize()<<endl;
-      fRecPoints->RemoveAt(fNumberOfECAClusters);
-      //cout<<"size fRecPoints after "<<fRecPoints->GetSize()<<endl;
-      fNumberOfECAClusters--;
-      nSplittedClusters--;
+          AliDebug(1,Form("NULL digit part3.3 or NULL energy=%f",correctedEnergyList[iparam/3*nDigits+iDigit]));
+        }
+      }//digit loop
 
+      if(recPoint->GetMultiplicity()==0){//recpoint exists but no digits associated -> remove from list
+       delete (*list)[newClusterIndex];
+       list->RemoveAt(newClusterIndex);
+       nSplittedClusters--;
+       newClusterIndex--;//decrease cluster number
+      }else {//recPoint exists and has digits associated -> very good increase number of clusters 
+       AliDebug(5,Form("cluster %d, digit no %d, energy %f",iparam/3,(recPoint->GetDigitsList())[0],(recPoint->GetEnergiesList())[0]));
+      }
+      
+    } else {//recPoint empty -> remove from list
+      AliError("NULL RecPoint");
+      //protection from recpoint with no digits
+      delete (*list)[newClusterIndex];
+      list->RemoveAt(newClusterIndex);
+      nSplittedClusters--;
+      newClusterIndex--;//decrease cluster number
     }
 
     iparam += 3 ;
+    newClusterIndex++;
   }//while
   
   delete[] fitparameters ;
   delete[] efit ;
   delete[] correctedEnergyList ;
 
-  return kTRUE;
+//  print 
+  AliDebug(5,Form("  nSplittedClusters %d, fNumberOfECAClusters %d, newClusterIndex %d,list->Entries() %d\n",nSplittedClusters,fNumberOfECAClusters,newClusterIndex,list->GetEntriesFast() ));
+  
+  //  cout<<"end of unfolding check part 3.3"<<endl;
+  return nSplittedClusters;
+}
+
+//____________________________________________________________________________
+Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, 
+                                         Int_t nMax, 
+                                         AliEMCALDigit ** maxAt, 
+                                         Float_t * maxAtEnergy)
+{
+  // Extended to whole EMCAL 
+  // Performs the unfolding of a cluster with nMax overlapping showers 
+  // Returns true if success (1->several clusters), otherwise false (fit failed)
+
+  TObjArray *list =new TObjArray(2);//temporary object
+  Int_t nUnfoldedClusters=UnfoldOneCluster(iniTower,nMax,maxAt,maxAtEnergy,list);
+
+  // here we write new clusters from list to fRecPoints 
+  AliDebug(5,Form("Number of clusters after unfolding %d",list->GetEntriesFast()));
+  Int_t iDigit=0;
+  AliEMCALDigit * digit = 0 ;
+  for(Int_t i=0;i<list->GetEntriesFast();i++) {
+    AliEMCALRecPoint * recPoint = 0 ;
+    
+    if(fNumberOfECAClusters >= fRecPoints->GetSize())
+      fRecPoints->Expand(2*fNumberOfECAClusters) ;
+    
+    //add recpoint
+    (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;//fNumberOfECAClusters-1 is old cluster before unfolding
+    recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
+               AliEMCALRecPoint * rpUFOne = dynamic_cast<AliEMCALRecPoint *>(list->At(i)) ;
+               
+    if( recPoint && rpUFOne ){//recPoint present -> good
+      
+                       recPoint->SetNExMax(list->GetEntriesFast()) ;
+                       
+      Int_t   *digitsList = rpUFOne->GetDigitsList();
+      Float_t *energyList = rpUFOne->GetEnergiesList();
+
+      if(!digitsList || ! energyList)
+      {
+        AliDebug(-1,"No digits index or energy available");
+                               delete (*fRecPoints)[fNumberOfECAClusters];
+                               fRecPoints->RemoveAt(fNumberOfECAClusters);
+        continue;
+      }
+      
+      AliDebug(5,Form("cluster %d, digit no %d, energy %f\n",i,digitsList[0],energyList[0]));
+
+      for(iDigit = 0 ; iDigit < rpUFOne->GetMultiplicity(); iDigit ++) {
+        digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
+        if(digit) recPoint->AddDigit( *digit, energyList[iDigit], kFALSE ) ; //FIXME, need to study the shared case
+      }//digit loop
+      fNumberOfECAClusters++ ; 
+    } else {//recPoint empty -> remove from list
+      AliError("NULL RecPoint");
+      delete (*fRecPoints)[fNumberOfECAClusters];
+      fRecPoints->RemoveAt(fNumberOfECAClusters);
+    }
+
+  }//loop over unfolded clusters
+  
+  //print energy of new unfolded clusters
+  AliDebug(5,Form("  nUnfoldedClusters %d, fNumberOfECAClusters %d",nUnfoldedClusters,fNumberOfECAClusters ));
+  for(Int_t inewclus=0; inewclus<nUnfoldedClusters;inewclus++){
+               AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters-1-inewclus));
+    if(rp) AliDebug(5,Form("  Unfolded cluster %d E %f",inewclus, rp->GetEnergy() ));
+  }
+
+  //clear tables  
+  list->SetOwner(kTRUE);
+  list->Delete();
+  delete list;
+  if(nUnfoldedClusters>1) return kTRUE;
+  return kFALSE;
 }
 
 
+
 //____________________________________________________________________________
 Bool_t AliEMCALUnfolding::UnfoldClusterV2old(AliEMCALRecPoint * iniTower, 
                                          Int_t nMax, 
@@ -564,8 +740,17 @@ Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit *
 
   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
        
-  if(!gMinuit)
-    gMinuit = new TMinuit(100) ;//max 100 parameters
+  if(!gMinuit){
+    //    gMinuit = new TMinuit(100) ;//max 100 parameters
+    if(nPar<30) gMinuit = new TMinuit(30);
+    else gMinuit = new TMinuit(nPar) ;//max nPar parameters
+    //
+  } else {
+    if(gMinuit->fMaxpar < nPar) {
+      delete gMinuit;
+      gMinuit = new TMinuit(nPar);
+    }
+  }
 
   gMinuit->mncler();                     // Reset Minuit's list of paramters
   gMinuit->SetPrintLevel(-1) ;           // No Printout
@@ -594,10 +779,16 @@ Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit *
   Int_t iphi    =  0 ;//x direction
   Int_t ieta    =  0 ;//z direstion
 
-  for(iDigit = 0; iDigit < nDigits; iDigit++){
+  for(iDigit = 0; iDigit < nDigits; iDigit++)
+  {
     digit = maxAt[iDigit];
-    if(digit==0) AliError("energy of digit = 0!");
-    fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
+    if(!digit)
+    {
+      AliError("Null digit pointer");
+      continue;
+    }
+    
+    fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta);
     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
                                       iIphi, iIeta,iphi,ieta);
 
@@ -607,7 +798,7 @@ Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit *
     gMinuit->mnparm(index, "x",  iphi, 0.05, 0, 0, ierflg) ;
     index++ ;
     if(ierflg != 0){
-      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %d", iphi ) ;
+      Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: x=%d, param.id=%d, nMaxima=%d",iphi,index-1,nPar/3 ) ;
       toMinuit->Clear();
       delete toMinuit ;
       return kFALSE;
@@ -616,7 +807,7 @@ Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit *
     gMinuit->mnparm(index, "z",  ieta, 0.05, 0, 0, ierflg) ;
     index++ ;
     if(ierflg != 0){
-      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %d", ieta) ;
+      Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: z=%d, param.id=%d, nMaxima=%d", ieta, index-1,nPar/3) ;
       toMinuit->Clear();
       delete toMinuit ;
       return kFALSE;
@@ -625,7 +816,7 @@ Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit *
     gMinuit->mnparm(index, "Energy",  energy , 0.001*energy, 0., 5.*energy, ierflg) ;//was 0.05
     index++ ;
     if(ierflg != 0){
-      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;
+      Error("FindFit", "EMCAL Unfolding unable to set initial value for fit procedure: energy = %f, param.id=%d, nMaxima=%d", energy, index-1, nPar/3) ;
       toMinuit->Clear();
       delete toMinuit ;
       return kFALSE;
@@ -646,7 +837,7 @@ Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit *
   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
   //gMinuit->mnexcm("MINI", &p0, 0, ierflg) ;    // minimize
   if(ierflg == 4){  // Minimum not found
-    Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;
+    AliDebug(1,"EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;
     toMinuit->Clear();
     delete toMinuit ;
     return kFALSE ;
@@ -654,12 +845,15 @@ Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit *
   for(index = 0; index < nPar; index++){
     Double_t err = 0. ;
     Double_t val = 0. ;
-    gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
+    gMinuit->GetParameter(index, val, err) ;    // Returns value and error ofOA parameter index
     fitparameters[index] = val ;
   }
 
   toMinuit->Clear();
   delete toMinuit ;
+
+  if(gMinuit->fMaxpar>30) delete gMinuit;
+
   return kTRUE;
 
 }
@@ -671,17 +865,17 @@ Double_t  AliEMCALUnfolding::ShowerShapeV2(Double_t x, Double_t y)
   // Shape of the shower
   // If you change this function, change also the gradient evaluation in ChiSquare()
 
-  Double_t r = fSSPars[7]*TMath::Sqrt(x*x+y*y);
-  Double_t rp1  = TMath::Power(r, fSSPars[1]) ;
-  Double_t rp5  = TMath::Power(r, fSSPars[5]) ;
-  Double_t shape = fSSPars[0]*TMath::Exp( -rp1 * (1. / (fSSPars[2] + fSSPars[3] * rp1) + fSSPars[4] / (1 + fSSPars[6] * rp5) ) ) ;
+  Double_t r = fgSSPars[7]*TMath::Sqrt(x*x+y*y);
+  Double_t rp1  = TMath::Power(r, fgSSPars[1]) ;
+  Double_t rp5  = TMath::Power(r, fgSSPars[5]) ;
+  Double_t shape = fgSSPars[0]*TMath::Exp( -rp1 * (1. / (fgSSPars[2] + fgSSPars[3] * rp1) + fgSSPars[4] / (1 + fgSSPars[6] * rp5) ) ) ;
   return shape ;
 }
 
 //____________________________________________________________________________
 void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,
-                                                Double_t & fret,
-                                                Double_t * x, Int_t iflag)
+                                            Double_t & fret,
+                                            Double_t * x, Int_t iflag)
 {
   // Calculates the Chi square for the cluster unfolding minimization
   // Number of parameters, Gradient, Chi squared, parameters, what to do
@@ -728,65 +922,65 @@ void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,
         digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );
         
         if(digit){
-        geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
-        geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
-                                          iIphi, iIeta,iphi,ieta);
-        EvalParsPhiDependence(digit->GetId(),geom);
-        
-        if(iflag == 2){  // calculate gradient
-          Int_t iParam = 0 ;
-          efit = 0. ;
-          while(iParam < nPar ){
-            Double_t dx = ((Float_t)iphi - x[iParam]) ;
-            iParam++ ;
-            Double_t dz = ((Float_t)ieta - x[iParam]) ;
-            iParam++ ;
-            efit += x[iParam] * ShowerShapeV2(dx,dz) ;
-            iParam++ ;
-          }
-          
-          Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)
-          iParam = 0 ;
-          while(iParam < nPar ){
-            Double_t xpar = x[iParam] ;
-            Double_t zpar = x[iParam+1] ;
-            Double_t epar = x[iParam+2] ;
-            
-            Double_t dr = fSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) );
-            Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
-            Double_t rp1  = TMath::Power(dr, fSSPars[1]) ;
-            Double_t rp5  = TMath::Power(dr, fSSPars[5]) ;
-            
-            Double_t deriv = -2 * TMath::Power(dr,fSSPars[1]-2.) * fSSPars[7] * fSSPars[7] * 
-            (fSSPars[1] * ( 1/(fSSPars[2]+fSSPars[3]*rp1) + fSSPars[4]/(1+fSSPars[6]*rp5) ) - 
-             (fSSPars[1]*fSSPars[3]*rp1/( (fSSPars[2]+fSSPars[3]*rp1)*(fSSPars[2]+fSSPars[3]*rp1) ) + 
-              fSSPars[4]*fSSPars[5]*fSSPars[6]*rp5/( (1+fSSPars[6]*rp5)*(1+fSSPars[6]*rp5) ) ) );
-            
-            //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) )
-            //                                                   - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;
-            
-            Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ;  // Derivative over x
-            iParam++ ;
-            Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ;  // Derivative over z
-            iParam++ ;
-            Grad[iParam] += shape ;                                  // Derivative over energy
-            iParam++ ;
-          }
-        }
-        efit = 0;
-        iparam = 0 ;
-        
-        while(iparam < nPar ){
-          Double_t xpar = x[iparam] ;
-          Double_t zpar = x[iparam+1] ;
-          Double_t epar = x[iparam+2] ;
-          iparam += 3 ;
-          efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
-        }
-        
-        fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;
-        // Here we assume, that sigma = sqrt(E) 
-        } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!\n");
+         geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); 
+         geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
+                                           iIphi, iIeta,iphi,ieta);
+         EvalParsPhiDependence(digit->GetId(),geom);
+         
+         if(iflag == 2){  // calculate gradient
+           Int_t iParam = 0 ;
+           efit = 0. ;
+           while(iParam < nPar ){
+             Double_t dx = ((Float_t)iphi - x[iParam]) ;
+             iParam++ ;
+             Double_t dz = ((Float_t)ieta - x[iParam]) ;
+             iParam++ ;
+             efit += x[iParam] * ShowerShapeV2(dx,dz) ;
+             iParam++ ;
+           }
+           
+           Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)
+           iParam = 0 ;
+           while(iParam < nPar ){
+             Double_t xpar = x[iParam] ;
+             Double_t zpar = x[iParam+1] ;
+             Double_t epar = x[iParam+2] ;
+             
+             Double_t dr = fgSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) );
+             Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
+             Double_t rp1  = TMath::Power(dr, fgSSPars[1]) ;
+             Double_t rp5  = TMath::Power(dr, fgSSPars[5]) ;
+             
+             Double_t deriv = -2 * TMath::Power(dr,fgSSPars[1]-2.) * fgSSPars[7] * fgSSPars[7] * 
+               (fgSSPars[1] * ( 1/(fgSSPars[2]+fgSSPars[3]*rp1) + fgSSPars[4]/(1+fgSSPars[6]*rp5) ) - 
+                (fgSSPars[1]*fgSSPars[3]*rp1/( (fgSSPars[2]+fgSSPars[3]*rp1)*(fgSSPars[2]+fgSSPars[3]*rp1) ) + 
+                 fgSSPars[4]*fgSSPars[5]*fgSSPars[6]*rp5/( (1+fgSSPars[6]*rp5)*(1+fgSSPars[6]*rp5) ) ) );
+             
+             //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) )
+             //                                                   - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;
+             
+             Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ;  // Derivative over x
+             iParam++ ;
+             Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ;  // Derivative over z
+             iParam++ ;
+             Grad[iParam] += shape ;                                  // Derivative over energy
+             iParam++ ;
+           }
+         }
+         efit = 0;
+         iparam = 0 ;
+         
+         while(iparam < nPar ){
+           Double_t xpar = x[iparam] ;
+           Double_t zpar = x[iparam+1] ;
+           Double_t epar = x[iparam+2] ;
+           iparam += 3 ;
+           efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;
+         }
+         
+         fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;
+         // Here we assume, that sigma = sqrt(E) 
+        } else printf("AliEMCALUnfoding::UnfoldingChiSquareV2 - NULL digit!, nPar %d \n", nPar); // put nPar here to cheat coverity and rule checker
       } // digit loop
     } // recpoint, digits and geom not NULL
   }// List is not NULL
@@ -797,20 +991,20 @@ void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,
 //____________________________________________________________________________
 void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){
   for(UInt_t i=0;i<7;++i)
-    fSSPars[i]=pars[i];
-  if(pars[2]==0. && pars[3]==0.) fSSPars[2]=1.;//to avoid dividing by 0
+    fgSSPars[i]=pars[i];
+  if(pars[2]==0. && pars[3]==0.) fgSSPars[2]=1.;//to avoid dividing by 0
 }
 
 //____________________________________________________________________________
 void AliEMCALUnfolding::SetPar5(Double_t *pars){
   for(UInt_t i=0;i<3;++i)
-    fPar5[i]=pars[i];
+    fgPar5[i]=pars[i];
 }
 
 //____________________________________________________________________________
 void AliEMCALUnfolding::SetPar6(Double_t *pars){
   for(UInt_t i=0;i<3;++i)
-    fPar6[i]=pars[i];
+    fgPar6[i]=pars[i];
 }
 
 //____________________________________________________________________________
@@ -820,7 +1014,7 @@ void AliEMCALUnfolding::EvalPar5(Double_t phi){
   //phi in degrees range (-10,10)
   //
   //fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936;
-  fSSPars[5] = fPar5[0] + phi * fPar5[1] + phi*phi * fPar5[2];
+  fgSSPars[5] = fgPar5[0] + phi * fgPar5[1] + phi*phi * fgPar5[2];
 }
 
 //____________________________________________________________________________
@@ -830,11 +1024,11 @@ void AliEMCALUnfolding::EvalPar6(Double_t phi){
   //phi in degrees range (-10,10)
   //
   //fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361;
-  fSSPars[6] = fPar6[0] + phi * fPar6[1] + phi*phi * fPar6[2];
+  fgSSPars[6] = fgPar6[0] + phi * fgPar6[1] + phi*phi * fgPar6[2];
 }
 
 //____________________________________________________________________________
-void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, AliEMCALGeometry *geom){
+void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, const AliEMCALGeometry *geom){
   //
   // calculate params p5 and p6 depending on the phi angle in global coordinate
   // for the cell with given absId index