]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALClusterizerFixedWindow.cxx
Changes needed for evaporation and fragmentation
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerFixedWindow.cxx
index 289e397c58c8ebb69b641cca771858bf976b8b9f..6b0448004a172c47090e873fdb48170201bc6de8 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-// This class derives from AliEMCALClustrerizer but also keeps the API of AliEMCALClusterizerv1
-// Algorithm:
-// 1. peek the most energetic cell
-// 2. assign it as a center of the cluster and add cells surrounding it: 3x3, 5x5...
-// 3. remove the cells contributing to the cluster
-// 4. start from 1 for the remaining clusters
-// 5. cluster splitting (not implemented yet) - use the shape analysis to resolve the energy sharing
-// - for high energy clusters check the surrounding of the 3x3 clusters for extra energy 
-// (merge 3x3 clusters and resolve the internal energy sharing - case for 2 clusters merged)
-// Use Case:
-//  root [0] AliEMCALClusterizerFixedWindow * cl = new AliEMCALClusterizerFixedWindow("galice.root")  
-//  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
-//               //reads gAlice from header file "..."                      
-//  root [1] cl->ExecuteTask()  
-//               //finds RecPoints in all events stored in galice.root
-//  root [2] cl->SetDigitsBranch("digits2") 
-//               //sets another title for Digitis (input) branch
-//  root [3] cl->SetRecPointsBranch("recp2")  
-//               //sets another title four output branches
-//  root [4] cl->SetTowerLocalMaxCut(0.03)  
-//               //set clusterization parameters
-//  root [5] cl->ExecuteTask("deb all time")  
-//               //once more finds RecPoints options are 
-//               // deb - print number of found rec points
-//               // deb all - print number of found RecPoints and some their characteristics 
-//               // time - print benchmarking results
+// This class derives from AliEMCALClustrerizer
 
-#include "AliEMCALClusterizerFixedWindow.h"
-
-// --- ROOT system ---
+// --- Root ---
 #include <TMath.h> 
 #include <TMinuit.h>
 #include <TTree.h> 
 #include <TClonesArray.h>
 #include <TH1I.h>
 
-// --- Standard library ---
-#include <cassert>
-//#include <iostream>
-//#include <fstream>
-
-// --- AliRoot header files ---
+// --- AliRoot ---
 #include "AliLog.h"
 #include "AliEMCALRecPoint.h"
 #include "AliEMCALDigit.h"
 #include "AliEMCALCalibData.h"
 #include "AliESDCaloCluster.h"
 #include "AliEMCALUnfolding.h"
-#include "AliEMCALFixedWindowClusterInfo.h"
+
+#include "AliEMCALClusterizerFixedWindow.h"
 
 ClassImp(AliEMCALClusterizerFixedWindow)
 
-//____________________________________________________________________________
-AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow()
-: AliEMCALClusterizer(), 
-nPhi(4), 
-nEta(4), 
-shiftPhi(2), 
-shiftEta(2),
-fTRUshift(0), 
-clusters_array(0),
-fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo"))
+//__________________________________________________________________________________________
+  AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow()
+    : AliEMCALClusterizer(), 
+      fNphi(4), 
+      fNeta(4), 
+      fShiftPhi(2), 
+      fShiftEta(2),
+      fTRUshift(0),
+      fClustersArray(0)
 {
-       // ctor with the indication of the file where header Tree and digits Tree are stored
-
+  // Constructor
 }
 
-//____________________________________________________________________________
+//__________________________________________________________________________________________
 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry)
-: AliEMCALClusterizer(geometry), 
-nPhi(4), 
-nEta(4), 
-shiftPhi(2), 
-shiftEta(2), 
-fTRUshift(0), 
-clusters_array(0),
-fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo"))
+  : AliEMCALClusterizer(geometry), 
+    fNphi(4), 
+    fNeta(4), 
+    fShiftPhi(2), 
+    fShiftEta(2),
+    fTRUshift(0),
+    fClustersArray(0)
 {
-       // ctor with the indication of the file where header Tree and digits Tree are stored
-       // use this contructor to avoid usage of Init() which uses runloader
-       // change needed by HLT - MP
+  // Constructor
 }
 
-//____________________________________________________________________________
+//__________________________________________________________________________________________
 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
-: AliEMCALClusterizer(geometry, calib, caloped),
-nPhi(4), 
-nEta(4), 
-shiftPhi(2), 
-shiftEta(2), 
-fTRUshift(0), 
-clusters_array(0),
-fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo"))
+  : AliEMCALClusterizer(geometry, calib, caloped),
+    fNphi(4), 
+    fNeta(4), 
+    fShiftPhi(2), 
+    fShiftEta(2),
+    fTRUshift(0),
+    fClustersArray(0)
 {
-       // ctor, geometry and calibration are initialized elsewhere.
+  // Constructor
 }
 
-//____________________________________________________________________________
+//__________________________________________________________________________________________
 AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
 {
-       // dtor
-  delete fClustersInfo;
+  // Destructor
+  
+  delete fClustersArray;
+}
+
+//__________________________________________________________________________________________
+void AliEMCALClusterizerFixedWindow::SetNphi (Int_t n) 
+{
+  // Set fNphi; if clusterizer already initialized gives a warning and does nothing
+  
+  if (fClustersArray)
+    AliWarning("Clusterizer already initialized. Unable to change the parameters.");
+  else
+    fNphi = n;
+}
+
+//__________________________________________________________________________________________
+void AliEMCALClusterizerFixedWindow::SetNeta (Int_t n) 
+{
+  // Set fNeta; if clusterizer already initialized gives a warning and does nothing
+  
+  if (fClustersArray)
+    AliWarning("Clusterizer already initialized. Unable to change the parameters.");
+  else
+    fNeta = n;
 }
 
-//____________________________________________________________________________
+//__________________________________________________________________________________________
+void AliEMCALClusterizerFixedWindow::SetShiftPhi (Int_t s) 
+{
+  // Set fShiftPhi; if clusterizer already initialized gives a warning and does nothing
+  
+  if (fClustersArray)
+    AliWarning("Clusterizer already initialized. Unable to change the parameters.");
+  else
+    fShiftPhi = s;
+}
+
+//__________________________________________________________________________________________
+void AliEMCALClusterizerFixedWindow::SetShiftEta (Int_t s) 
+{
+  // Set fShiftEta; if clusterizer already initialized gives a warning and does nothing
+  
+  if (fClustersArray)
+    AliWarning("Clusterizer already initialized. Unable to change the parameters.");
+  else
+    fShiftEta = s;
+}
+
+//__________________________________________________________________________________________
+void AliEMCALClusterizerFixedWindow::SetTRUshift(Bool_t b) 
+{
+  // Set fTRUshift; if clusterizer already initialized gives a warning and does nothing
+  
+  if (fClustersArray)
+    AliWarning("Clusterizer already initialized. Unable to change the parameters.");
+  else
+    fTRUshift = b;
+}
+
+
+//__________________________________________________________________________________________
 void AliEMCALClusterizerFixedWindow::Digits2Clusters(Option_t * option)
 {
-       // Steering method to perform clusterization for the current event 
-       // in AliEMCALLoader
+  // Steering method to perform clusterization for the current event 
        
-       if (strstr(option,"tim"))
-               gBenchmark->Start("EMCALClusterizer"); 
+  if (strstr(option,"tim"))
+    gBenchmark->Start("EMCALClusterizer"); 
        
-       if (strstr(option,"print"))
-               Print(""); 
+  if (strstr(option,"print"))
+    Print(""); 
        
-       //Get calibration parameters from file or digitizer default values.
-       GetCalibrationParameters();
+  //Get calibration parameters from file or digitizer default values.
+  GetCalibrationParameters();
        
-       //Get dead channel map from file or digitizer default values.
-       GetCaloCalibPedestal();
+  //Get dead channel map from file or digitizer default values.
+  GetCaloCalibPedestal();
        
-       MakeClusters();  //only the real clusters
+  MakeClusters();  //only the real clusters
        
-       if (fToUnfold) {
-               fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
-               fClusterUnfolding->MakeUnfolding();
-       }
-       
-       //Evaluate position, dispersion and other RecPoint properties for EC section 
-       for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) { 
-               AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
-               if (rp) {
-                       rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
-                       AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
-                       //For each rec.point set the distance to the nearest bad crystal
-                       if (fCaloPed)
-                               rp->EvalDistanceToBadChannels(fCaloPed);
-               }
-       }
+  if (fToUnfold) {
+    fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
+    fClusterUnfolding->MakeUnfolding();
+  }
        
-       //fRecPoints->Sort();
+  //Evaluate position, dispersion and other RecPoint properties for EC section 
+  for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) { 
+    AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
+    if (rp) {
+      rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
+      AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
+      //For each rec.point set the distance to the nearest bad crystal
+      if (fCaloPed)
+        rp->EvalDistanceToBadChannels(fCaloPed);
+    }
+  }
+  
+  fRecPoints->Sort();
        
-       for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
-               AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
-               if (rp) {
-                       rp->SetIndexInList(index);
-               }
-               else AliFatal("RecPoint NULL!!");
-       }
+  for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
+    AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
+    if (rp) {
+      rp->SetIndexInList(index);
+    }
+    else AliFatal("RecPoint NULL!!");
+  }
        
-       if (fTreeR)
-               fTreeR->Fill();
+  if (fTreeR)
+    fTreeR->Fill();
        
-       if (strstr(option,"deb") || strstr(option,"all"))  
-               PrintRecPoints(option);
+  if (strstr(option,"deb") || strstr(option,"all"))  
+    PrintRecPoints(option);
        
-       AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
+  AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
        
-       if (strstr(option,"tim")) {
-               gBenchmark->Stop("EMCALClusterizer");
-               printf("Exec took %f seconds for Clusterizing", 
-                          gBenchmark->GetCpuTime("EMCALClusterizer"));
-       }    
+  if (strstr(option,"tim")) {
+    gBenchmark->Stop("EMCALClusterizer");
+    printf("Exec took %f seconds for Clusterizing", 
+           gBenchmark->GetCpuTime("EMCALClusterizer"));
+  }    
 }
 
-//____________________________________________________________________________
+//__________________________________________________________________________________________
 void AliEMCALClusterizerFixedWindow::MakeClusters()
 {
-       // Make clusters
+  // Make clusters
        
-       if (fGeom == 0) 
-               AliFatal("Did not get geometry from EMCALLoader");
+  if (fGeom == 0) 
+    AliFatal("Did not get geometry from EMCALLoader");
        
-       fNumberOfECAClusters = 0;
-       fRecPoints->Delete();
-  
-  if (fClustersInfo->GetLastElementId() > 0)
-    fClustersInfo->Clear();
+  fNumberOfECAClusters = 0;
+  fRecPoints->Delete();
        
-       Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0, iphi=0, ieta=0;
+  Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0, iphi=0, ieta=0;
        
-       // Defining geometry and clusterization parameter
-       Int_t nEtaDigitsSupMod = fGeom->GetNEta() * fGeom->GetNETAdiv(); // always 48?;
-       Int_t nPhiDigitsSupMod = fGeom->GetNPhi() * fGeom->GetNPHIdiv(); // always 24?;
+  // Defining geometry and clusterization parameter
+  Int_t nEtaDigitsSupMod = fGeom->GetNEta() * fGeom->GetNETAdiv(); // always 48?;
+  Int_t nPhiDigitsSupMod = fGeom->GetNPhi() * fGeom->GetNPHIdiv(); // always 24?;
   
   Int_t nTRUPhi = 1;
   Int_t nTRUEta = 1;
@@ -210,8 +226,7 @@ void AliEMCALClusterizerFixedWindow::MakeClusters()
   Int_t nEtaDigits = nEtaDigitsSupMod * fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
   Int_t nPhiDigits = nPhiDigitsSupMod * fGeom->GetNPhiSuperModule();    
   
-  if (fTRUshift)
-  {
+  if (fTRUshift){
     nTRUPhi = fGeom->GetNPhiSuperModule() * 3;
     nTRUEta = fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
     nEtaDigits /= nTRUEta;
@@ -219,184 +234,181 @@ void AliEMCALClusterizerFixedWindow::MakeClusters()
   }
 
   // Check if clusterizer parameter are compatible with calorimeter geometry
-  if (nEtaDigits < nEta)
-       {
-               AliFatal(Form("Error: nEta = %d is greater than nEtaDigits = %d.",nEta,nEtaDigits));
-               return;
-       }
-       if (nPhiDigits < nPhi)
-       {
-               AliFatal(Form("Error: nPhi = %d is greater than nPhiDigits = %d.",nPhi,nPhiDigits));
-               return;
-       }
-       if (nEtaDigits % shiftEta != 0)
-       {
-               AliFatal(Form("Error: shiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",shiftEta,nEtaDigits));
-               return;
-       }
-       if (nPhiDigits % shiftPhi != 0)
-       {
-               AliFatal(Form("Error: shiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",shiftPhi,nPhiDigits));
-               return;
-       }
-       if (nEta % shiftEta != 0)
-       {
-               AliFatal(Form("Error: shiftEta = %d is not divisor of nEta = %d.",shiftEta,nEta));
-               return;
-       }
-       if (nPhi % shiftPhi != 0)
-       {
-               AliFatal(Form("Error: shiftPhi = %d is not divisor of nPhi = %d).",shiftPhi,nPhi));
-               return;
-       }
+  if (nEtaDigits < fNeta){
+    AliFatal(Form("Error: fNeta = %d is greater than nEtaDigits = %d.",fNeta,nEtaDigits));
+    return;
+  }
+  if (nPhiDigits < fNphi){
+    AliFatal(Form("Error: fNphi = %d is greater than nPhiDigits = %d.",fNphi,nPhiDigits));
+    return;
+  }
+  if (nEtaDigits % fShiftEta != 0){
+    AliFatal(Form("Error: fShiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",fShiftEta,nEtaDigits));
+    return;
+  }
+  if (nPhiDigits % fShiftPhi != 0){
+    AliFatal(Form("Error: fShiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",fShiftPhi,nPhiDigits));
+    return;
+  }
+  if (fNeta % fShiftEta != 0){
+    AliFatal(Form("Error: fShiftEta = %d is not divisor of fNeta = %d.",fShiftEta,fNeta));
+    return;
+  }
+  if (fNphi % fShiftPhi != 0){
+    AliFatal(Form("Error: fShiftPhi = %d is not divisor of fNphi = %d).",fShiftPhi,fNphi));
+    return;
+  }
   
-  Int_t maxiShiftPhi = nPhi / shiftPhi;
-  Int_t maxiShiftEta = nEta / shiftEta;
+  Int_t maxiShiftPhi = fNphi / fShiftPhi;
+  Int_t maxiShiftEta = fNeta / fShiftEta;
        
-       Int_t nDigitsCluster = nPhi * nEta;
+  Int_t nDigitsCluster = fNphi * fNeta;
   
-  Int_t nClusEtaNoShift = nEtaDigits / nEta;
-  Int_t nClusPhiNoShift = nPhiDigits / nPhi;
+  Int_t nClusEtaNoShift = nEtaDigits / fNeta;
+  Int_t nClusPhiNoShift = nPhiDigits / fNphi;
   
   Int_t nClusters =  nClusEtaNoShift * nClusPhiNoShift * nTRUEta * nTRUPhi;
   
   Int_t nTotalClus = nClusters * maxiShiftEta * maxiShiftPhi;
   
-  if (!clusters_array)
-  {
-    clusters_array = new AliEMCALDigit**[nTotalClus];
+  if (!fClustersArray) {
+    fClustersArray = new AliEMCALDigit**[nTotalClus];
     for (Int_t i = 0; i < nTotalClus; i++)
     {
-      clusters_array[i] = NULL;
+      fClustersArray[i] = NULL;
     }
   }
   
-  AliEMCALDigit *digit = 0;
+  // Set up TObjArray with pointers to digits to work on calibrated digits 
+  TObjArray *digitsC = new TObjArray();
+  AliEMCALDigit *digit;
+  Float_t dEnergyCalibrated = 0.0, ehs = 0.0, time = 0.0;
+  TIter nextdigit(fDigitsArr);
+  while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigit()))) { // calibrate and clean up digits
+    dEnergyCalibrated =  digit->GetAmplitude();
+    time              =  digit->GetTime();
+    Calibrate(dEnergyCalibrated, time, digit->GetId());
+    digit->SetCalibAmp(dEnergyCalibrated);
+    digit->SetTime(time);
+    if (dEnergyCalibrated < fMinECut) {
+      continue;
+    }
+    else if (!fGeom->CheckAbsCellId(digit->GetId())) {
+      continue;
+    }
+    else {
+      ehs += dEnergyCalibrated;
+      digitsC->AddLast(digit);
+    }
+  } 
   
-  for (Int_t ishiftPhi = 0; ishiftPhi < maxiShiftPhi; ishiftPhi++)
-  {
-    Int_t nClusPhi = (nPhiDigits - shiftPhi * ishiftPhi) / nPhi;
+  AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %f\n",
+                  fDigitsArr->GetEntries(),fMinECut,ehs));
+   
+  for (Int_t ishiftPhi = 0; ishiftPhi < maxiShiftPhi; ishiftPhi++){
+    Int_t nClusPhi = (nPhiDigits - fShiftPhi * ishiftPhi) / fNphi;
     
-    for (Int_t ishiftEta = 0; ishiftEta < maxiShiftEta; ishiftEta++)
-    {
+    for (Int_t ishiftEta = 0; ishiftEta < maxiShiftEta; ishiftEta++) {
       
-      Int_t nClusEta = (nEtaDigits - shiftEta * ishiftEta) / nEta; 
+      Int_t nClusEta = (nEtaDigits - fShiftEta * ishiftEta) / fNeta; 
       
       Int_t iTotalClus = nClusters * (ishiftPhi * maxiShiftEta + ishiftEta);
       
-      TIter nextdigit(fDigitsArr);
-      
-      nextdigit.Reset();
-      
-      while (digit = static_cast<AliEMCALDigit*>(nextdigit()))
-      {
+      TIter nextdigitC(digitsC);
+      while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigitC()))) { // scan over the list of digitsC
+        
         fGeom->GetCellIndex (digit->GetId(), nSupMod, nModule, nIphi, nIeta);
         fGeom->GetCellPhiEtaIndexInSModule (nSupMod, nModule, nIphi, nIeta, iphi, ieta);
         
-        Int_t iphi_eff = iphi - shiftPhi * ishiftPhi + nPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
+        Int_t iphi_eff = iphi - fShiftPhi * ishiftPhi + nPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
         
         Int_t iTRUphi = iphi_eff / nPhiDigits;
         
         iphi_eff -= iTRUphi * nPhiDigits;
         
-        Int_t iClusPhi = iphi_eff / nPhi; 
+        Int_t iClusPhi = iphi_eff / fNphi; 
         
         if (iphi_eff < 0 || iClusPhi >= nClusPhi) 
           continue;
         
-        Int_t ieta_eff = ieta - shiftEta * ishiftEta + nEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
+        Int_t ieta_eff = ieta - fShiftEta * ishiftEta + nEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
         
         Int_t iTRUeta = ieta_eff / nEtaDigits;
         
         ieta_eff -= iTRUeta * nEtaDigits;
         
-        Int_t iClusEta = ieta_eff / nEta; 
+        Int_t iClusEta = ieta_eff / fNeta; 
         
         if (ieta_eff < 0 || iClusEta >= nClusEta) 
           continue;
         
         iphi_eff += iTRUphi * nPhiDigits;
-        iClusPhi = iphi_eff / nPhi; 
+        iClusPhi = iphi_eff / fNphi; 
         
         ieta_eff += iTRUeta * nEtaDigits;
-        iClusEta = ieta_eff / nEta; 
+        iClusEta = ieta_eff / fNeta; 
         
         Int_t iCluster = iClusPhi + iClusEta * nClusPhiNoShift * nTRUPhi; 
-        Int_t iDigit = iphi_eff % nPhi + (ieta_eff % nEta) * nPhi;
+        Int_t iDigit = iphi_eff % fNphi + (ieta_eff % fNeta) * fNphi;
 
-        
-        if (iCluster >= nClusters)
-        {
-          AliFatal(Form("ERROR: iCluster out of range! iCluster = %d, nClusters = %d", iCluster, nClusters));
+        if (iCluster >= nClusters){
+          AliWarning(Form("iCluster out of range! iCluster = %d, nClusters = %d", iCluster, nClusters));
           return;
         }
         
         iCluster += iTotalClus;
         
-        if (clusters_array[iCluster] == NULL)
-        {
+        if (fClustersArray[iCluster] == NULL){
           fNumberOfECAClusters++;
-          clusters_array[iCluster] = new AliEMCALDigit*[nDigitsCluster];
-          for (Int_t i = 0; i < nDigitsCluster; i++)
-          {
-            clusters_array[iCluster][i] = NULL;
+          fClustersArray[iCluster] = new AliEMCALDigit*[nDigitsCluster];
+          for (Int_t i = 0; i < nDigitsCluster; i++){
+            fClustersArray[iCluster][i] = NULL;
           }
-          
-          fClustersInfo->Add(iCluster, -1, iClusEta, iClusPhi);
         }
         
-        if (clusters_array[iCluster][iDigit] != NULL)
-        {
-          AliFatal("ERROR: digit already added!");
+        if (fClustersArray[iCluster][iDigit] != NULL){
+          AliWarning("Digit already added!");
           return;
         }
         
-        clusters_array[iCluster][iDigit] = digit;
+        fClustersArray[iCluster][iDigit] = digit;
         
       } // loop on digit
       
     } // loop on eta shift
     
-       } // loop on phi shift
+  } // loop on phi shift
        
+  if(fNumberOfECAClusters >= fRecPoints->GetSize()) 
+    fRecPoints->Expand(fNumberOfECAClusters+1);
+  
   Int_t iRecPoint = 0;
-       for (Int_t iCluster = 0; iCluster < nTotalClus; iCluster++)
-       {
+  for (Int_t iCluster = 0; iCluster < nTotalClus; iCluster++) {
     
-               if (clusters_array[iCluster] == NULL) continue;
-               
-               (*fRecPoints)[iRecPoint] = new AliEMCALRecPoint("");
-               AliEMCALRecPoint *recPoint = dynamic_cast<AliEMCALRecPoint*> (fRecPoints->At(iRecPoint));
+    if (fClustersArray[iCluster] == NULL) 
+      continue;
+
+    (*fRecPoints)[iRecPoint] = new AliEMCALRecPoint("");
+    AliEMCALRecPoint *recPoint = dynamic_cast<AliEMCALRecPoint*> (fRecPoints->At(iRecPoint));
                
-               if (recPoint) 
-               {
-      if (fClustersInfo->ContainsIndex(iRecPoint))
-        AliFatal(Form("ERROR: index present already, %d", iRecPoint));
-      
-      fClustersInfo->SetIndexFromId(iCluster, iRecPoint);
-      
-                       iRecPoint++;       
-                       recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
-                       // note: this way the sharing info is lost!
-                       for (Int_t iDigit = 0; iDigit < nDigitsCluster; iDigit++)
-                       {
-                               if (clusters_array[iCluster][iDigit] == NULL) continue;
-                               digit = clusters_array[iCluster][iDigit];
-        Float_t dEnergyCalibrated = digit->GetAmplitude();
-        Float_t time              = digit->GetTime();
-        Calibrate(dEnergyCalibrated,time,digit->GetId());
-                               digit->SetCalibAmp(dEnergyCalibrated);
-                               recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE); //Time or TimeR?
-        clusters_array[iCluster][iDigit] = NULL;
-                       }
-               }
+    if (recPoint) {
+      iRecPoint++;       
+      recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
+      recPoint->SetUniqueID(iCluster);
+
+      for (Int_t iDigit = 0; iDigit < nDigitsCluster; iDigit++){
+        if (fClustersArray[iCluster][iDigit] == NULL) continue;
+        digit = fClustersArray[iCluster][iDigit];
+        recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
+        fClustersArray[iCluster][iDigit] = NULL;
+      }
+    }
     
-    delete[] clusters_array[iCluster];
-    clusters_array[iCluster] = NULL;
-       }
-  
-       AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f)\n",
-                                       fDigitsArr->GetEntries(),fMinECut));
-       
-       AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); 
+    delete[] fClustersArray[iCluster];
+    fClustersArray[iCluster] = NULL;
+  }
+
+  delete digitsC;
+  AliDebug(1, Form("MakeClusters: Number of digits %d  -> (e %f)\n", fDigitsArr->GetEntries(),fMinECut));
+  AliDebug(1, Form("total no of clusters %d from %d digits", fNumberOfECAClusters, fDigitsArr->GetEntriesFast())); 
 }