making the container name parameterized
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / AliCalorimeterUtils.cxx
index bf6417d..1322e8d 100755 (executable)
 
 
 // --- ROOT system ---
-#include "TGeoManager.h"
+#include <TGeoManager.h>
+#include <TH2F.h>
+#include <TCanvas.h>
+#include <TStyle.h>
+#include <TPad.h>
+#include <TFile.h>
 
 // --- ANALYSIS system ---
 #include "AliCalorimeterUtils.h"
@@ -34,6 +39,9 @@
 #include "AliVCluster.h"
 #include "AliVCaloCells.h"
 #include "AliMixedEvent.h"
+#include "AliAODCaloCluster.h"
+#include "AliOADBContainer.h"
+#include "AliAnalysisManager.h"
 
 // --- Detector ---
 #include "AliEMCALGeometry.h"
@@ -45,21 +53,24 @@ ClassImp(AliCalorimeterUtils)
 //____________________________________________
   AliCalorimeterUtils::AliCalorimeterUtils() : 
     TObject(), fDebug(0), 
-    fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
-    fPHOSGeoName ("PHOSgeo"), 
+    fEMCALGeoName(""),
+    fPHOSGeoName (""), 
     fEMCALGeo(0x0),                   fPHOSGeo(0x0), 
     fEMCALGeoMatrixSet(kFALSE),       fPHOSGeoMatrixSet(kFALSE), 
     fLoadEMCALMatrices(kFALSE),       fLoadPHOSMatrices(kFALSE),
     fRemoveBadChannels(kFALSE),       fPHOSBadChannelMap(0x0), 
     fNCellsFromPHOSBorder(0),
     fNMaskCellColumns(0),             fMaskCellColumns(0x0),
-    fRecalibration(kFALSE),           fPHOSRecalibrationFactors(),
+    fRecalibration(kFALSE),           fRunDependentCorrection(kFALSE), fPHOSRecalibrationFactors(),
     fEMCALRecoUtils(new AliEMCALRecoUtils),
     fRecalculatePosition(kFALSE),     fCorrectELinearity(kFALSE),
     fRecalculateMatching(kFALSE),
     fCutR(20),                        fCutZ(20),
     fCutEta(20),                      fCutPhi(20),
-    fLocMaxCutE(0),                   fLocMaxCutEDiff(0)
+    fLocMaxCutE(0),                   fLocMaxCutEDiff(0),
+    fPlotCluster(0),                  fOADBSet(kFALSE),
+    fOADBForEMCAL(kFALSE),            fOADBForPHOS(kFALSE),
+    fOADBFilePathEMCAL(""),           fOADBFilePathPHOS("")
 {
   //Ctor
   
@@ -93,6 +104,372 @@ AliCalorimeterUtils::~AliCalorimeterUtils()
   
 }
 
+//____________________________________________________
+void AliCalorimeterUtils::AccessOADB(AliVEvent* event)
+{
+  // Set the AODB calibration, bad channels etc. parameters at least once
+  // alignment matrices from OADB done in SetGeometryMatrices
+  
+  //Set it only once
+  if(fOADBSet) return ; 
+  
+  Int_t   runnumber = event->GetRunNumber() ;
+  TString pass      = GetPass();
+  
+  // EMCAL
+  if(fOADBForEMCAL)
+  {
+    printf("AliCalorimeterUtils::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePathEMCAL.Data(),runnumber,pass.Data());
+    
+    Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
+    
+    // Bad map
+    if(fRemoveBadChannels)
+    {
+      AliOADBContainer *contBC=new AliOADBContainer("");
+      contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePathEMCAL.Data()),"AliEMCALBadChannels"); 
+      
+      TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
+      
+      if(arrayBC)
+      {
+        SwitchOnDistToBadChannelRecalculation();
+        printf("AliCalorimeterUtils::SetOADBParameters() - Remove EMCAL bad cells \n");
+        
+        for (Int_t i=0; i<nSM; ++i) 
+        {
+          TH2I *hbm = GetEMCALChannelStatusMap(i);
+          
+          if (hbm)
+            delete hbm;
+          
+          hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
+          
+          if (!hbm) 
+          {
+            AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
+            continue;
+          }
+          
+          hbm->SetDirectory(0);
+          SetEMCALChannelStatusMap(i,hbm);
+          
+        } // loop
+      } else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT remove EMCAL bad channels\n"); // run array
+    }  // Remove bad
+    
+    // Energy Recalibration
+    if(fRecalibration)
+    {
+      AliOADBContainer *contRF=new AliOADBContainer("");
+      
+      contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRecalib");
+      
+      TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber); 
+      
+      if(recal)
+      {
+        TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
+        
+        if(recalpass)
+        {
+          TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
+          
+          if(recalib)
+          {
+            printf("AliCalorimeterUtils::SetOADBParameters() - Recalibrate EMCAL \n");
+            for (Int_t i=0; i<nSM; ++i) 
+            {
+              TH2F *h = GetEMCALChannelRecalibrationFactors(i);
+              
+              if (h)
+                delete h;
+              
+              h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
+              
+              if (!h) 
+              {
+                AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
+                continue;
+              }
+              
+              h->SetDirectory(0);
+              
+              SetEMCALChannelRecalibrationFactors(i,h);
+            } // SM loop
+          }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params object array\n"); // array ok
+        }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for pass\n"); // array pass ok
+      }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for run\n");  // run number array ok
+      
+      // once set, apply run dependent corrections if requested
+      //fEMCALRecoUtils->SetRunDependentCorrections(runnumber);
+            
+    } // Recalibration on
+    
+    // Energy Recalibration, apply on top of previous calibration factors
+    if(fRunDependentCorrection)
+    {
+      AliOADBContainer *contRFTD=new AliOADBContainer("");
+      
+      contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRunDepTempCalibCorrections");
+      
+      TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber); 
+      
+      if(htd)
+      {
+        printf("AliCalorimeterUtils::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
+        
+        for (Int_t ism=0; ism<nSM; ++ism) 
+        {        
+          for (Int_t icol=0; icol<48; ++icol) 
+          {        
+            for (Int_t irow=0; irow<24; ++irow) 
+            {
+              Float_t factor = GetEMCALChannelRecalibrationFactor(ism,icol,irow);
+              
+              Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
+              factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
+              //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID, 
+              //      GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
+              SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
+            } // columns
+          } // rows 
+        } // SM loop
+      }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n"); 
+    } // Run by Run T calibration    
+    
+    // Time Recalibration
+    if(fEMCALRecoUtils->IsTimeRecalibrationOn())
+    {
+      AliOADBContainer *contTRF=new AliOADBContainer("");
+      
+      contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALTimeCalib");
+      
+      TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber); 
+      
+      if(trecal)
+      {
+        TString passTmp = pass;
+        if(pass!="pass1" && pass!="pass2") passTmp = "pass2"; // TEMPORARY FIX FOR LHC11a analysis
+        
+        TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passTmp);
+        
+        if(trecalpass)
+        {
+          printf("AliCalorimeterUtils::SetOADBParameters() - Time Recalibrate EMCAL \n");
+          for (Int_t ibc = 0; ibc < 4; ++ibc) 
+          {
+            TH1F *h = GetEMCALChannelTimeRecalibrationFactors(ibc);
+            
+            if (h)
+              delete h;
+            
+            h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
+            
+            if (!h) 
+            {
+              AliError(Form("Could not load hAllTimeAvBC%d",ibc));
+              continue;
+            }
+            
+            h->SetDirectory(0);
+            
+            SetEMCALChannelTimeRecalibrationFactors(ibc,h);
+          } // bunch crossing loop
+        }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for pass\n"); // array pass ok
+      }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for run\n");  // run number array ok
+      
+    } // Recalibration on    
+    
+  }// EMCAL
+  
+  // PHOS
+  if(fOADBForPHOS)
+  {
+    printf("AliCalorimeterUtils::SetOADBParameters() - Get AODB parameters from PHOS in %s for run %d, and <%s> \n",fOADBFilePathPHOS.Data(),runnumber,pass.Data());
+    
+    // Bad map
+    if(fRemoveBadChannels)
+    {
+      AliOADBContainer badmapContainer(Form("phosBadMap"));
+      TString fileName="$ALICE_ROOT/OADB/PHOS/PHOSBadMaps.root";
+      badmapContainer.InitFromFile(Form("%s/PHOSBadMaps.root",fOADBFilePathPHOS.Data()),"phosBadMap");
+      
+      //Use a fixed run number from year 2010, this year not available yet.
+      TObjArray *maps = (TObjArray*)badmapContainer.GetObject(139000,"phosBadMap");
+      if(!maps)
+      {
+        printf("AliCalorimeterUtils::SetOADBParameters() - Can not read PHOS bad map for run %d.\n",runnumber) ;    
+      }
+      else
+      {
+        printf("AliCalorimeterUtils::SetOADBParameters() - Setting PHOS bad map with name %s \n",maps->GetName()) ;
+        for(Int_t mod=1; mod<5;mod++)
+        {
+          TH2I *hbmPH = GetPHOSChannelStatusMap(mod);
+          
+          if(hbmPH) 
+            delete hbmPH ;  
+          
+          hbmPH = (TH2I*)maps->At(mod);
+          
+          if(hbmPH) hbmPH->SetDirectory(0);
+          
+          SetPHOSChannelStatusMap(mod-1,hbmPH);
+          
+        } // modules loop  
+      } // maps exist
+    } // Remove bad channels
+  } // PHOS
+  
+  // Parameters already set once, so do not it again
+  fOADBSet = kTRUE;
+  
+}  
+
+//_____________________________________________________________
+void AliCalorimeterUtils::AccessGeometry(AliVEvent* inputEvent) 
+{
+  //Set the calorimeters transformation matrices and init geometry
+  
+  // First init the geometry, a priory not done before
+  Int_t runnumber = inputEvent->GetRunNumber() ;
+  InitPHOSGeometry (runnumber);
+  InitEMCALGeometry(runnumber);
+  
+  //Get the EMCAL transformation geometry matrices from ESD 
+  if(!fEMCALGeoMatrixSet && fEMCALGeo)
+  {
+    if(fLoadEMCALMatrices)
+    {
+      printf("AliCalorimeterUtils::AccessGeometry() - Load user defined EMCAL geometry matrices\n");
+      
+      // OADB if available
+      AliOADBContainer emcGeoMat("AliEMCALgeo");
+      emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePathEMCAL.Data()),"AliEMCALgeo");
+      TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
+      
+      for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
+      {
+        if (!fEMCALMatrix[mod]) // Get it from OADB
+        {
+          if(fDebug > 1 ) 
+            printf("AliCalorimeterUtils::AccessGeometry() - EMCAL matrices SM %d, %p\n",
+                   mod,((TGeoHMatrix*) matEMCAL->At(mod)));
+          //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
+          
+          fEMCALMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
+        }
+        
+        if(fEMCALMatrix[mod])
+        {
+          if(fDebug > 1) 
+            fEMCALMatrix[mod]->Print();
+          
+          fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;  
+        }
+      }//SM loop
+      
+      fEMCALGeoMatrixSet = kTRUE;//At least one, so good
+      
+    }//Load matrices
+    else if (!gGeoManager) 
+    { 
+      if(fDebug > 1) 
+        printf(" AliCalorimeterUtils::AccessGeometry() - Load EMCAL misalignment matrices. \n");
+      if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  
+      {
+        for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
+        { 
+          //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
+          if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) 
+          {
+            fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
+          }
+        }// loop over super modules    
+        
+        fEMCALGeoMatrixSet = kTRUE;//At least one, so good
+        
+      }//ESD as input
+      else 
+      {
+        if(fDebug > 1)
+          printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
+      }//AOD as input
+    }//Get matrix from data
+    else if(gGeoManager)
+    {
+      fEMCALGeoMatrixSet = kTRUE;
+    }
+  }//EMCAL geo && no geoManager
+  
+       //Get the PHOS transformation geometry matrices from ESD 
+  if(!fPHOSGeoMatrixSet && fPHOSGeo)
+  {
+    if(fLoadPHOSMatrices)
+    {
+      printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
+      
+      // OADB if available
+      AliOADBContainer geomContainer("phosGeo");
+      geomContainer.InitFromFile(Form("%s/PHOSGeometry.root",fOADBFilePathPHOS.Data()),"PHOSRotationMatrixes");
+      TObjArray *matPHOS = (TObjArray*)geomContainer.GetObject(139000,"PHOSRotationMatrixes");    
+      
+      for(Int_t mod = 0 ; mod < 5 ; mod++)
+      {
+        if (!fPHOSMatrix[mod]) // Get it from OADB
+        {
+          if(fDebug > 1 ) 
+            printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - PHOS matrices module %d, %p\n",
+                   mod,((TGeoHMatrix*) matPHOS->At(mod)));
+          //((TGeoHMatrix*) matPHOS->At(mod))->Print();
+          
+          fPHOSMatrix[mod] = (TGeoHMatrix*) matPHOS->At(mod) ;
+        }
+        
+        // Set it, if it exists
+        if(fPHOSMatrix[mod])
+        {
+          if(fDebug > 1 ) 
+            fPHOSMatrix[mod]->Print();
+          
+          fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;  
+        }      
+      }// SM loop
+      
+      fPHOSGeoMatrixSet = kTRUE;//At least one, so good
+      
+    }//Load matrices
+    else if (!gGeoManager) 
+    { 
+      if(fDebug > 1) 
+        printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
+                       if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  
+      {
+                               for(Int_t mod = 0; mod < 5; mod++)
+        { 
+                                       if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) 
+          {
+                                               //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
+                                               fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
+                                       }
+                               }// loop over modules
+        fPHOSGeoMatrixSet  = kTRUE; //At least one so good
+                       }//ESD as input
+                       else 
+      {
+                               if(fDebug > 1) 
+                                       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
+      }//AOD as input
+    }// get matrix from data
+    else if(gGeoManager)
+    {
+      fPHOSGeoMatrixSet = kTRUE;
+    }
+       }//PHOS geo     and  geoManager was not set
+  
+}
+
 //______________________________________________________________________________________
 Bool_t AliCalorimeterUtils::AreNeighbours(const TString calo, 
                                           const Int_t absId1, const Int_t absId2 ) const
@@ -169,10 +546,9 @@ Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
       for(Int_t i = 0; i < cluster->GetNCells() ; i++){
         Int_t absId = cluster->GetCellAbsId(i) ;
         for (Int_t j = startCell; j < endCell ;  j++) {
-          Short_t cellNumber; 
-          Double_t amp ; 
-          Double_t time; 
-          cells->GetCell(j, cellNumber, amp, time) ; 
+          Short_t cellNumber, mclabel; 
+          Double_t amp, time, efrac; 
+          cells->GetCell(j, cellNumber, amp, time,mclabel,efrac) ; 
           if (absId == cellNumber) {
             if(amp > ampMax){
               ampMax   = amp;
@@ -218,25 +594,39 @@ Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
     
                //Check rows/phi
     Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
-               if(iSM < 10){
+               if(iSM < 10)
+    {
                        if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE; 
     }
-               else{
-                       if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE; 
+               else
+    {
+      if(fEMCALGeoName.Contains("12SM")) // 1/3 SM
+      {
+        if(iphi >= nborder && iphi < 8-nborder) okrow =kTRUE; 
+      }
+      else // 1/2 SM
+      {
+        if(iphi >= nborder && iphi <12-nborder) okrow =kTRUE; 
+      }
                }
                
                //Check columns/eta
-               if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
+               if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0())
+    {
                        if(ieta  > nborder && ieta < 48-nborder) okcol =kTRUE; 
                }
-               else{
-                       if(iSM%2==0){
+               else
+    {
+                       if(iSM%2==0)
+      {
                                if(ieta >= nborder)     okcol = kTRUE;  
                        }
-                       else {
+                       else 
+      {
                                if(ieta <  48-nborder)  okcol = kTRUE;  
                        }
                }//eta 0 not checked
+    
                if(fDebug > 1)
                {
                        printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
@@ -437,7 +827,8 @@ Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEv
              particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
                return fEMCALGeo->GetSuperModuleNumber(absId) ;
        }//EMCAL
-       else if(particle->GetDetector()=="PHOS"){
+       else if(particle->GetDetector()=="PHOS")
+  {
     // In case we use the MC reader, the input are TParticles, 
     // in this case use the corresponing method in PHOS Geometry to get the particle.
     if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
@@ -480,10 +871,12 @@ Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
        //Get the EMCAL/PHOS module number that corresponds to this cluster
        TLorentzVector lv;
        Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
-  if(!cluster){
+  if(!cluster)
+  {
     if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
     return -1;
   }
+  
        cluster->GetMomentum(lv,v);
        Float_t phi = lv.Phi();
        if(phi < 0) phi+=TMath::TwoPi();        
@@ -495,9 +888,11 @@ Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
              lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
                return fEMCALGeo->GetSuperModuleNumber(absId) ;
        }//EMCAL
-       else if(cluster->IsPHOS()) {
+       else if(cluster->IsPHOS()) 
+  {
                Int_t    relId[4];
-               if ( cluster->GetNCells() > 0) {
+               if ( cluster->GetNCells() > 0) 
+    {
                        absId = cluster->GetCellAbsId(0);
                        if(fDebug > 2) 
                                printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
@@ -505,7 +900,8 @@ Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
                }
                else return -1;
                
-               if ( absId >= 0) {
+               if ( absId >= 0) 
+    {
                        fPHOSGeo->AbsToRelNumbering(absId,relId);
                        if(fDebug > 2) 
                          printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
@@ -523,32 +919,47 @@ Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const T
 {
        //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
        Int_t imod = -1;
-       if ( absId >= 0) {
-               if(calo=="EMCAL"){
+       if ( absId >= 0) 
+  {
+               if(calo=="EMCAL")
+    {
                        Int_t iTower = -1, iIphi = -1, iIeta = -1; 
                        fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
                        fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
-      if(imod < 0 || irow < 0 || icol < 0 ) {
+      if(imod < 0 || irow < 0 || icol < 0 ) 
+      {
         Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
       }
       
                        //RCU0
-                       if (0<=irow&&irow<8) iRCU=0; // first cable row
-                       else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half; 
-                       //second cable row
-                       //RCU1
-                       else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half; 
-                       //second cable row
-                       else if(16<=irow&&irow<24) iRCU=1; // third cable row
-                       
-                       if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
-                       if (iRCU<0) {
+      if(imod < 10 )
+      {
+        if      (0<=irow&&irow<8)                       iRCU=0; // first cable row
+        else if (8<=irow&&irow<16 &&  0<=icol&&icol<24) iRCU=0; // first half; 
+        //second cable row
+        //RCU1
+        else if (8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half; 
+        //second cable row
+        else if (16<=irow&&irow<24)                     iRCU=1; // third cable row
+        
+        if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
+      }
+      else 
+      {
+        // Last 2 SM have one single SRU, just assign RCU 0
+        iRCU = 0 ;
+      }
+
+                       if (iRCU<0) 
+      {
                                Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
                        }                       
                        
                        return imod ;
+      
                }//EMCAL
-               else{//PHOS
+               else //PHOS
+    {
                        Int_t    relId[4];
                        fPHOSGeo->AbsToRelNumbering(absId,relId);
                        irow = relId[2];
@@ -556,7 +967,8 @@ Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const T
                        imod = relId[0]-1;
                        iRCU= (Int_t)(relId[2]-1)/16 ;
                        //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
-                       if (iRCU >= 4) {
+                       if (iRCU >= 4) 
+      {
                                Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
                        }                       
                        return imod;
@@ -572,14 +984,11 @@ Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCalo
   // Find local maxima in cluster
   
   const Int_t   nc = cluster->GetNCells();
-  Int_t   *absIdList = new Int_t  [nc]; 
-  Float_t *maxEList  = new Float_t[nc]; 
+  Int_t   absIdList[nc]; 
+  Float_t maxEList[nc]; 
   
   Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
   
-  delete [] absIdList;
-  delete [] maxEList;
-  
   return nMax;
   
 }
@@ -589,7 +998,7 @@ Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCalo
                                                   Int_t *absIdList,     Float_t *maxEList) 
 {
   // Find local maxima in cluster
-
+  
   Int_t iDigitN = 0 ;
   Int_t iDigit  = 0 ;
   Int_t absId1 = -1 ;
@@ -600,14 +1009,22 @@ Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCalo
   if(!cluster->IsEMCAL()) calorimeter = "PHOS";
   
   //printf("cluster : ncells %d \n",nCells);
+  
+  Float_t emax  = 0;
+  Int_t   idmax =-1;
   for(iDigit = 0; iDigit < nCells ; iDigit++)
   {
     absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit]  ; 
-     //Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
-     //RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);  
-     //Int_t icol = -1, irow = -1, iRCU = -1;
-     //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
-     //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
+    Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
+    RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);  
+    if( en > emax )
+    {
+      emax  = en ;
+      idmax = absIdList[iDigit] ;
+    }
+    //Int_t icol = -1, irow = -1, iRCU = -1;
+    //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
+    //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
   }
   
   for(iDigit = 0 ; iDigit < nCells; iDigit++) 
@@ -618,7 +1035,7 @@ Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCalo
       
       Float_t en1 = cells->GetCellAmplitude(absId1);
       RecalibrateCellAmplitude(en1,calorimeter,absId1);  
-              
+      
       //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
       
       for(iDigitN = 0; iDigitN < nCells; iDigitN++) 
@@ -631,7 +1048,7 @@ Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCalo
         
         Float_t en2 = cells->GetCellAmplitude(absId2);
         RecalibrateCellAmplitude(en2,calorimeter,absId2);
-                
+        
         //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
         
         if ( AreNeighbours(calorimeter, absId1, absId2) ) 
@@ -679,20 +1096,68 @@ Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCalo
     }
   }
   
-  //printf("**********N maxima %d \n",iDigitN);
-  //for(Int_t imax = 0; imax < iDigitN; imax++) printf("imax %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
+  if(iDigitN == 0)
+  {
+    if(fDebug > 0) 
+      printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f\n",
+             idmax,emax,cluster->E());
+    iDigitN      = 1     ;
+    maxEList[0]  = emax  ;
+    absIdList[0] = idmax ; 
+  }
+  
+  if(fDebug > 0) 
+  {    
+    printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - In cluster E %2.2f, M02 %2.2f, M20 %2.2f, N maxima %d \n", 
+           cluster->E(),cluster->GetM02(),cluster->GetM20(), iDigitN);
+  
+    if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++) 
+    {
+      printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
+    }
+  }
   
   return iDigitN ;
   
 }
 
+//____________________________________
+TString AliCalorimeterUtils::GetPass()
+{
+  // Get passx from filename.
+    
+  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()) 
+  {
+    AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null\n");
+    return TString("");
+  }
+  
+  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()) 
+  {
+    AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null\n");
+    return TString("");
+  }
+  
+  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
+  if      (pass.Contains("ass1")) return TString("pass1");
+  else if (pass.Contains("ass2")) return TString("pass2");
+  else if (pass.Contains("ass3")) return TString("pass3");
+  else if (pass.Contains("ass4")) return TString("pass4");
+  else if (pass.Contains("ass5")) return TString("pass5");
+
+  // No condition fullfilled, give a default value
+  printf("AliCalorimeterUtils::GetPass() - Pass number string not found \n");
+  return TString("");            
+  
+}
+
 //________________________________________
 void AliCalorimeterUtils::InitParameters()
 {
   //Initialize the parameters of the analysis.
   
-  fEMCALGeoName         = "EMCAL_COMPLETEV1";
-  fPHOSGeoName          = "PHOSgeo";
+  fEMCALGeoName         = "";
+  fPHOSGeoName          = "";
   
   fEMCALGeoMatrixSet    = kFALSE;
   fPHOSGeoMatrixSet     = kFALSE;
@@ -711,6 +1176,13 @@ void AliCalorimeterUtils::InitParameters()
   //  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols; 
   //  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols; 
   
+  fOADBSet      = kFALSE;
+  fOADBForEMCAL = kTRUE ;            
+  fOADBForPHOS  = kFALSE;
+  
+  fOADBFilePathEMCAL = "$ALICE_ROOT/OADB/EMCAL" ;          
+  fOADBFilePathPHOS  = "$ALICE_ROOT/OADB/PHOS"  ;     
+  
 }
 
 
@@ -760,34 +1232,76 @@ void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
 }
 
 
-//___________________________________________
-void AliCalorimeterUtils::InitEMCALGeometry()
+//__________________________________________________________
+void AliCalorimeterUtils::InitEMCALGeometry(Int_t runnumber)
 {
        //Initialize EMCAL geometry if it did not exist previously
-       if (!fEMCALGeo){
+  
+       if (!fEMCALGeo)
+  {
+    if(fEMCALGeoName=="")
+    {
+      if     (runnumber <  140000 && 
+              runnumber >= 100000)   fEMCALGeoName = "EMCAL_FIRSTYEARV1";
+      else if(runnumber >= 140000 &&
+              runnumber <  171000)   fEMCALGeoName = "EMCAL_COMPLETEV1";
+      else                           fEMCALGeoName = "EMCAL_COMPLETE12SMV1";  
+      printf("AliCalorimeterUtils::InitEMCALGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fEMCALGeoName.Data(),runnumber);
+    }
+    
                fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
-               if(fDebug > 0){
-                       printf("AliCalorimeterUtils::InitEMCALGeometry()");
+    
+               if(fDebug > 0)
+    {
+                       printf("AliCalorimeterUtils::InitEMCALGeometry(run=%d)",runnumber);
                        if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
                        printf("\n");
                }
        }
 }
 
-//__________________________________________
-void AliCalorimeterUtils::InitPHOSGeometry()
+//_________________________________________________________
+void AliCalorimeterUtils::InitPHOSGeometry(Int_t runnumber)
 {
        //Initialize PHOS geometry if it did not exist previously
-       if (!fPHOSGeo){
+  
+       if (!fPHOSGeo)
+  {
+    if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
+      
                fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName); 
-               if(fDebug > 0){
-                       printf("AliCalorimeterUtils::InitPHOSGeometry()");
+    
+               if(fDebug > 0)
+    {
+                       printf("AliCalorimeterUtils::InitPHOSGeometry(run=%d)",runnumber);
                        if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
                        printf("\n");
                }       
        }       
 }
 
+//__________________________________________________________________
+Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,  
+                                             const Int_t ieta) const 
+{
+  //Check if cell is in one of the regions where we have significant amount 
+  //of material in front. Only EMCAL
+  
+  Int_t icol = ieta;
+  if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
+  
+  if (fNMaskCellColumns && fMaskCellColumns) 
+  {
+    for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) 
+    {
+      if(icol==fMaskCellColumns[imask]) return kTRUE;
+    }
+  }
+  
+  return kFALSE;
+  
+}
+
 //_________________________________________________________
 void AliCalorimeterUtils::Print(const Option_t * opt) const
 {
@@ -801,7 +1315,7 @@ void AliCalorimeterUtils::Print(const Option_t * opt) const
   printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
          fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
   if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
-  printf("Recalibrate Clusters? %d\n",fRecalibration);
+  printf("Recalibrate Clusters? %d, run by run  %d\n",fRecalibration,fRunDependentCorrection);
   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
@@ -812,26 +1326,6 @@ void AliCalorimeterUtils::Print(const Option_t * opt) const
   printf("    \n") ;
 } 
 
-//__________________________________________________________________
-Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,  
-                                             const Int_t ieta) const 
-{
-  //Check if cell is in one of the regions where we have significant amount 
-  //of material in front. Only EMCAL
-  
-  Int_t icol = ieta;
-  if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
-  
-  if (fNMaskCellColumns && fMaskCellColumns) {
-    for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
-      if(icol==fMaskCellColumns[imask]) return kTRUE;
-    }
-  }
-  
-  return kFALSE;
-  
-}
-
 //__________________________________________________________________________________________
 void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp,
                                                    const TString calo, const Int_t id) const
@@ -869,7 +1363,6 @@ void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time,
   
 }
 
-
 //__________________________________________________________________________
 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster, 
                                                       AliVCaloCells * cells)
@@ -922,89 +1415,6 @@ Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
        return energy;
 }
 
-//________________________________________________________________________________
-void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent) 
-{
-  //Set the calorimeters transformation matrices
-  
-  //Get the EMCAL transformation geometry matrices from ESD 
-  if(!fEMCALGeoMatrixSet && fEMCALGeo){
-    if(fLoadEMCALMatrices){
-      printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined EMCAL geometry matrices\n");
-      for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
-        if(fEMCALMatrix[mod]){
-          if(fDebug > 1) 
-            fEMCALMatrix[mod]->Print();
-          fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;  
-        }
-      }//SM loop
-      fEMCALGeoMatrixSet = kTRUE;//At least one, so good
-      
-    }//Load matrices
-    else if (!gGeoManager) { 
-      
-      if(fDebug > 1) 
-        printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
-      if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
-        for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ 
-          //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
-          if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) {
-            fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
-          }
-        }// loop over super modules    
-        
-        fEMCALGeoMatrixSet = kTRUE;//At least one, so good
-        
-      }//ESD as input
-      else {
-        if(fDebug > 1)
-          printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
-      }//AOD as input
-    }//Get matrix from data
-    else if(gGeoManager){
-      fEMCALGeoMatrixSet = kTRUE;
-    }
-  }//EMCAL geo && no geoManager
-  
-       //Get the PHOS transformation geometry matrices from ESD 
-  if(!fPHOSGeoMatrixSet && fPHOSGeo){
-
-    if(fLoadPHOSMatrices){
-      printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
-      for(Int_t mod = 0 ; mod < 5 ; mod++){
-        if(fPHOSMatrix[mod]){
-          if(fDebug > 1 ) 
-            fPHOSMatrix[mod]->Print();
-          fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;  
-        }
-      }//SM loop
-      fPHOSGeoMatrixSet = kTRUE;//At least one, so good
-
-    }//Load matrices
-    else if (!gGeoManager) { 
-      if(fDebug > 1) 
-        printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
-                       if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
-                               for(Int_t mod = 0; mod < 5; mod++){ 
-                                       if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) {
-                                               //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
-                                               fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
-                                       }
-                               }// loop over modules
-        fPHOSGeoMatrixSet  = kTRUE; //At least one so good
-                       }//ESD as input
-                       else {
-                               if(fDebug > 1) 
-                                       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
-      }//AOD as input
-    }// get matrix from data
-    else if(gGeoManager){
-      fPHOSGeoMatrixSet = kTRUE;
-    }
-       }//PHOS geo     and  geoManager was not set
-  
-}
-
 //__________________________________________________________________________________________
 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
 {
@@ -1021,7 +1431,8 @@ void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
 { 
   //Recalculate track matching
   
-  if (fRecalculateMatching) {
+  if (fRecalculateMatching) 
+  {
     fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo)   ; 
     //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
     //if(esdevent){
@@ -1030,3 +1441,263 @@ void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
     //}
   }
 }
+
+//___________________________________________________________________________
+void AliCalorimeterUtils::SplitEnergy(const Int_t absId1, const Int_t absId2,
+                                      AliVCluster* cluster, 
+                                      AliVCaloCells* cells,
+                                      //Float_t & e1, Float_t & e2,
+                                      AliAODCaloCluster* cluster1,
+                                      AliAODCaloCluster* cluster2,
+                                      const Int_t nMax,
+                                      const Int_t eventNumber)
+{
+  
+  // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2 
+  // maxima are too close and have common cells, split the energy between the 2
+  
+  TH2F* hClusterMap    = 0 ;
+  TH2F* hClusterLocMax = 0 ;
+  TH2F* hCluster1      = 0 ;
+  TH2F* hCluster2      = 0 ;
+  
+  if(fPlotCluster)
+  {
+    hClusterMap    = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
+    hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
+    hCluster1      = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
+    hCluster2      = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
+    hClusterMap    ->SetXTitle("column");
+    hClusterMap    ->SetYTitle("row");
+    hClusterLocMax ->SetXTitle("column");
+    hClusterLocMax ->SetYTitle("row");
+    hCluster1      ->SetXTitle("column");
+    hCluster1      ->SetYTitle("row");
+    hCluster2      ->SetXTitle("column");
+    hCluster2      ->SetYTitle("row");
+  }
+  
+  TString calorimeter = "EMCAL";
+  if(cluster->IsPHOS())
+  {
+    calorimeter="PHOS";
+    printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
+    return;
+  }
+  
+  const Int_t ncells  = cluster->GetNCells();  
+  Int_t absIdList[ncells]; 
+  
+  Float_t e1 = 0,  e2   = 0 ;
+  Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;  
+  Float_t eCluster = 0;
+  Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1; 
+  for(Int_t iDigit  = 0; iDigit < ncells; iDigit++ ) 
+  {
+    absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
+    
+    Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
+    RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
+    eCluster+=ec;
+    
+    if(fPlotCluster) 
+    {
+      //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
+      sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
+      if(icol > maxCol) maxCol = icol;
+      if(icol < minCol) minCol = icol;
+      if(irow > maxRow) maxRow = irow;
+      if(irow < minRow) minRow = irow;
+      hClusterMap->Fill(icol,irow,ec);
+    }
+    
+  }
+  
+  // Init counters and variables
+  Int_t ncells1 = 1 ;
+  UShort_t absIdList1[9] ;  
+  Double_t fracList1 [9] ;  
+  absIdList1[0] = absId1 ;
+  fracList1 [0] = 1. ;
+  
+  Float_t ecell1 = cells->GetCellAmplitude(absId1);
+  RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
+  e1 =  ecell1;  
+  
+  Int_t ncells2 = 1 ;
+  UShort_t absIdList2[9] ;  
+  Double_t fracList2 [9] ; 
+  absIdList2[0] = absId2 ;
+  fracList2 [0] = 1. ;
+  
+  Float_t ecell2 = cells->GetCellAmplitude(absId2);
+  RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
+  e2 =  ecell2;  
+  
+  if(fPlotCluster)
+  {
+    Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
+    sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
+    hClusterLocMax->Fill(icol1,irow1,ecell1);
+    sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
+    hClusterLocMax->Fill(icol2,irow2,ecell2);
+  }
+  
+  // Very rough way to share the cluster energy
+  Float_t eRemain = (eCluster-ecell1-ecell2)/2;
+  Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
+  Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
+  
+  for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
+  {
+    Int_t absId = absIdList[iDigit];
+    
+    if(absId==absId1 || absId==absId2 || absId < 0) continue;
+    
+    Float_t ecell = cells->GetCellAmplitude(absId);
+    RecalibrateCellAmplitude(ecell, calorimeter, absId);
+    
+    if(AreNeighbours(calorimeter, absId1,absId ))
+    { 
+      absIdList1[ncells1]= absId;
+      
+      if(AreNeighbours(calorimeter, absId2,absId ))
+      { 
+        fracList1[ncells1] = shareFraction1; 
+        e1 += ecell*shareFraction1;
+      }
+      else 
+      {
+        fracList1[ncells1] = 1.; 
+        e1 += ecell;
+      }
+      
+      ncells1++;
+      
+    } // neigbour to cell1
+    
+    if(AreNeighbours(calorimeter, absId2,absId ))
+    { 
+      absIdList2[ncells2]= absId;
+      
+      if(AreNeighbours(calorimeter, absId1,absId ))
+      { 
+        fracList2[ncells2] = shareFraction2; 
+        e2 += ecell*shareFraction2;
+      }
+      else
+      { 
+        fracList2[ncells2] = 1.; 
+        e2 += ecell;
+      }
+      
+      ncells2++;
+      
+    } // neigbour to cell2
+    
+  }
+  
+  if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::SplitEnergy() - n Local Max %d, Cluster energy  = %f, Ecell1 = %f, Ecell2 = %f, Enew1 = %f, Enew2 = %f, Remain %f, \n ncells %d, ncells1 %d, ncells2 %d, f1 %f, f2  %f, sum f12 = %f \n",
+                            nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
+  
+  cluster1->SetE(e1);
+  cluster2->SetE(e2);  
+  
+  cluster1->SetNCells(ncells1);
+  cluster2->SetNCells(ncells2);  
+  
+  cluster1->SetCellsAbsId(absIdList1);
+  cluster2->SetCellsAbsId(absIdList2);
+  
+  cluster1->SetCellsAmplitudeFraction(fracList1);
+  cluster2->SetCellsAmplitudeFraction(fracList2);
+  
+  if(calorimeter=="EMCAL")
+  {
+    GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
+    GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
+  }
+  
+  if(fPlotCluster)
+  {
+    //printf("Cells of cluster1: ");
+    for(Int_t iDigit  = 0; iDigit < ncells1; iDigit++ ) 
+    {
+      //printf(" %d ",absIdList1[iDigit]);
+      
+      sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
+      
+      if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) )
+        hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit])*shareFraction1);
+      else 
+        hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
+    }
+    
+    //printf(" \n ");
+    //printf("Cells of cluster2: ");
+    
+    for(Int_t iDigit  = 0; iDigit < ncells2; iDigit++ ) 
+    {
+      //printf(" %d ",absIdList2[iDigit]);
+      
+      sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
+      if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) )
+        hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit])*shareFraction2);
+      else
+        hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
+      
+    }
+    //printf(" \n ");
+    
+    gStyle->SetPadRightMargin(0.1);
+    gStyle->SetPadLeftMargin(0.1);
+    gStyle->SetOptStat(0);
+    gStyle->SetOptFit(000000);
+    
+    if(maxCol-minCol > maxRow-minRow)
+    {
+      maxRow+= maxCol-minCol;
+    }
+    else 
+    {
+      maxCol+= maxRow-minRow;
+    }
+    
+    TCanvas  * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
+    c->Divide(2,2);  
+    c->cd(1);
+    gPad->SetGridy();
+    gPad->SetGridx();
+    hClusterMap    ->SetAxisRange(minCol, maxCol,"X");
+    hClusterMap    ->SetAxisRange(minRow, maxRow,"Y");
+    hClusterMap    ->Draw("colz");
+    c->cd(2);
+    gPad->SetGridy();
+    gPad->SetGridx();
+    hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
+    hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
+    hClusterLocMax ->Draw("colz");
+    c->cd(3);
+    gPad->SetGridy();
+    gPad->SetGridx();
+    hCluster1      ->SetAxisRange(minCol, maxCol,"X");
+    hCluster1      ->SetAxisRange(minRow, maxRow,"Y");
+    hCluster1      ->Draw("colz");
+    c->cd(4);
+    gPad->SetGridy();
+    gPad->SetGridx();
+    hCluster2      ->SetAxisRange(minCol, maxCol,"X");
+    hCluster2      ->SetAxisRange(minRow, maxRow,"Y");
+    hCluster2      ->Draw("colz");
+    
+    if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
+                                   eventNumber,cluster->E(),nMax,ncells1,ncells2));
+    
+    delete c;
+    delete hClusterMap;
+    delete hClusterLocMax;
+    delete hCluster1;
+    delete hCluster2;
+  }
+}
+