]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALClusterizerv1.cxx
changes to cope with effC++
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerv1.cxx
index 6461b45967dc5fa6f2c09ee10b58c4451ab1b098..700aacc0857cdd5303d3247f992dad84bc3201de 100644 (file)
  **************************************************************************/
 
 /* $Id$ */
-/* $Log:
-   1 October 2000. Yuri Kharlov:
-     AreNeighbours()
-     PPSD upper layer is considered if number of layers>1
-   18 October 2000. Yuri Kharlov:
-     AliEMCALClusterizerv1()
-     CPV clusterizing parameters added
-     MakeClusters()
-     After first PPSD digit remove EMC digits only once
-*/
+
 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
 //  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
 //                           of new  IO (à la PHOS)
 
 // --- ROOT system ---
 
-#include "TROOT.h" 
-#include "TFile.h" 
-#include "TFolder.h" 
-#include "TMath.h" 
-#include "TMinuit.h"
-#include "TTree.h" 
-#include "TSystem.h" 
-#include "TBenchmark.h"
+class TROOT;
+#include <TH1.h>
+#include <TFile.h> 
+class TFolder;
+#include <TMath.h> 
+#include <TMinuit.h>
+#include <TTree.h> 
+class TSystem; 
+#include <TBenchmark.h>
+#include <TBrowser.h>
 
 // --- Standard library ---
 
 
 // --- AliRoot header files ---
-
+#include "AliRunLoader.h"
+#include "AliRun.h"
+#include "AliESD.h"
+#include "AliEMCALLoader.h"
 #include "AliEMCALClusterizerv1.h"
+#include "AliEMCALRecPoint.h"
 #include "AliEMCALDigit.h"
 #include "AliEMCALDigitizer.h"
-#include "AliEMCALTowerRecPoint.h"
 #include "AliEMCAL.h"
-#include "AliEMCALGetter.h"
 #include "AliEMCALGeometry.h"
-#include "AliRun.h"
+#include "AliEMCALHistoUtilities.h"
+#include "AliCDBManager.h"
+
+class AliCDBStorage;
+#include "AliCDBEntry.h"
 
 ClassImp(AliEMCALClusterizerv1)
-  
+
 //____________________________________________________________________________
-  AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
+AliEMCALClusterizerv1::AliEMCALClusterizerv1() 
+  : AliEMCALClusterizer(),
+    fHists(0),fPointE(0),fPointL1(0),fPointL2(0),
+    fPointDis(0),fPointMult(0),fDigitAmp(0),fMaxE(0),
+    fMaxL1(0),fMaxL2(0),fMaxDis(0),fGeom(0),
+    fDefaultInit(kTRUE),
+    fToUnfold(kFALSE),
+    fNumberOfECAClusters(0),fNTowerInGroup(0),fCalibData(0),
+    fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
+    fECAW0(0.),fRecPointsInRun(0),fTimeGate(0.),fMinECut(0.)
 {
   // default ctor (to be used mainly by Streamer)
   
   InitParameters() ; 
-  fDefaultInit = kTRUE ; 
+  fGeom = AliEMCALGeometry::GetInstance();
+  fGeom->GetTransformationForSM(); // Global <-> Local
 }
 
 //____________________________________________________________________________
-AliEMCALClusterizerv1::AliEMCALClusterizerv1(const char* headerFile, const char* name, const Bool_t toSplit)
-:AliEMCALClusterizer(headerFile, name, toSplit)
+AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName)
+  : AliEMCALClusterizer(alirunFileName, eventFolderName),
+    fHists(0),fPointE(0),fPointL1(0),fPointL2(0),
+    fPointDis(0),fPointMult(0),fDigitAmp(0),fMaxE(0),
+    fMaxL1(0),fMaxL2(0),fMaxDis(0),fGeom(0),
+    fDefaultInit(kFALSE),
+    fToUnfold(kFALSE),
+    fNumberOfECAClusters(0),fNTowerInGroup(0),fCalibData(0),
+    fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
+    fECAW0(0.),fRecPointsInRun(0),fTimeGate(0.),fMinECut(0.)
 {
   // ctor with the indication of the file where header Tree and digits Tree are stored
   
   InitParameters() ; 
   Init() ;
-  fDefaultInit = kFALSE ; 
+}
 
+//____________________________________________________________________________
+AliEMCALClusterizerv1::AliEMCALClusterizerv1(const AliEMCALClusterizerv1& clus)
+  : AliEMCALClusterizer(clus),
+    fHists(clus.fHists),
+    fPointE(clus.fPointE),
+    fPointL1(clus.fPointL1),
+    fPointL2(clus.fPointL2),
+    fPointDis(clus.fPointDis),
+    fPointMult(clus.fPointMult),
+    fDigitAmp(clus.fDigitAmp),
+    fMaxE(clus.fMaxE),
+    fMaxL1(clus.fMaxL1),
+    fMaxL2(clus.fMaxL2),
+    fMaxDis(clus.fMaxDis),
+    fGeom(clus.fGeom),
+    fDefaultInit(clus.fDefaultInit),
+    fToUnfold(clus.fToUnfold),
+    fNumberOfECAClusters(clus.fNumberOfECAClusters),
+    fNTowerInGroup(clus.fNTowerInGroup),
+    fCalibData(clus.fCalibData),
+    fADCchannelECA(clus.fADCchannelECA),
+    fADCpedestalECA(clus.fADCpedestalECA),
+    fECAClusteringThreshold(clus.fECAClusteringThreshold),
+    fECALocMaxCut(clus.fECALocMaxCut),
+    fECAW0(clus.fECAW0),
+    fRecPointsInRun(clus.fRecPointsInRun),
+    fTimeGate(clus.fTimeGate),
+    fMinECut(clus.fMinECut)
+{
+  //copy ctor
 }
 
 //____________________________________________________________________________
   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
 {
   // dtor
-  fSplitFile = 0 ; 
-  
 }
 
 //____________________________________________________________________________
 const TString AliEMCALClusterizerv1::BranchName() const 
 { 
-  TString branchName(GetName() ) ;
-  branchName.Remove(branchName.Index(Version())-1) ;
-  return branchName ;
+   return GetName();
 }
 
 //____________________________________________________________________________
-Float_t  AliEMCALClusterizerv1::Calibrate(Int_t amp, Bool_t inpresho) const
+Float_t  AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) 
 {
-  //To be replased later by the method, reading individual parameters from the database
-  if ( inpresho ) // calibrate as pre shower
-    return -fADCpedestalPreSho + amp * fADCchannelPreSho ; 
-  else //calibrate as tower 
-    return -fADCpedestalTower + amp * fADCchannelTower ;                
+  // Convert digitized amplitude into energy.
+  // Calibration parameters are taken from calibration data base for raw data,
+  // or from digitizer parameters for simulated data.
+
+  if(fCalibData){
+
+    //JLK 13-Mar-2006
+    //We now get geometry at a higher level
+    //
+    // Loader
+    //    AliRunLoader *rl = AliRunLoader::GetRunLoader();
+    
+    // Load EMCAL Geomtry
+    //    rl->LoadgAlice(); 
+    //AliRun   * gAlice = rl->GetAliRun(); 
+    //AliEMCAL * emcal  = (AliEMCAL*)gAlice->GetDetector("EMCAL");
+    //AliEMCALGeometry * geom = emcal->GetGeometry();
+    
+    if (fGeom==0)
+      AliFatal("Did not get geometry from EMCALLoader") ;
+    
+    Int_t iSupMod = -1;
+    Int_t nTower  = -1;
+    Int_t nIphi   = -1;
+    Int_t nIeta   = -1;
+    Int_t iphi    = -1;
+    Int_t ieta    = -1;
+    
+    Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nTower, nIphi, nIeta) ;
+    if(!bCell)
+      Error("DigitizeEnergy","Wrong cell id number") ;
+    fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nTower,nIphi, nIeta,iphi,ieta);
+
+    fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
+    fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
+    return -fADCpedestalECA + amp * fADCchannelECA ;        
+  }
+  else //Return energy with default parameters if calibration is not available
+    return -fADCpedestalECA + amp * fADCchannelECA ; 
+  
 }
 
 //____________________________________________________________________________
 void AliEMCALClusterizerv1::Exec(Option_t * option)
 {
-  // Steering method
-
-  if( strcmp(GetName(), "")== 0 ) 
-    Init() ;
+  // Steering method to perform clusterization for events
+  // in the range from fFirstEvent to fLastEvent.
+  // This range is optionally set by SetEventRange().
+  // if fLastEvent=-1 (by default), then process events until the end.
 
   if(strstr(option,"tim"))
     gBenchmark->Start("EMCALClusterizer"); 
@@ -141,55 +221,60 @@ void AliEMCALClusterizerv1::Exec(Option_t * option)
   if(strstr(option,"print"))
     Print("") ; 
 
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
-  if(gime->BranchExists("RecPoints"))
-    return ;
-  Int_t nevents = gime->MaxEvent() ;
-  Int_t ievent ;
+  AliRunLoader *rl = AliRunLoader::GetRunLoader();
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 
-  for(ievent = 0; ievent < nevents; ievent++){
+  //Get calibration parameters from file or digitizer default values.
+  GetCalibrationParameters() ;
 
-    gime->Event(ievent,"D") ;
+  if (fLastEvent == -1) 
+    fLastEvent = rl->GetNumberOfEvents() - 1;
+  Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 
-    if(ievent == 0)
-      GetCalibrationParameters() ;
+  Int_t ievent ;
+  rl->LoadDigits("EMCAL");
+  for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
+    rl->GetEvent(ievent);
+
+    fNumberOfECAClusters = 0;
+
+    if(strstr(option,"pseudo"))
+      MakeClusters("pseudo") ;  //both types
+    else
+      MakeClusters("") ;  //only the real clusters
 
-    fNumberOfTowerClusters = fNumberOfPreShoClusters = 0 ;
-           
-    MakeClusters() ;
-    
     if(fToUnfold)
       MakeUnfolding() ;
 
-    WriteRecPoints(ievent) ;
+    WriteRecPoints() ;
 
     if(strstr(option,"deb"))  
       PrintRecPoints(option) ;
 
-    //increment the total number of digits per run 
-    fRecPointsInRun += gime->TowerRecPoints()->GetEntriesFast() ;  
-    fRecPointsInRun += gime->PreShowerRecPoints()->GetEntriesFast() ;  
- }
+    //increment the total number of recpoints per run   
+    fRecPointsInRun += emcalLoader->RecPoints()->GetEntriesFast() ;  
+  }
   
+  Unload();
+
   if(strstr(option,"tim")){
     gBenchmark->Stop("EMCALClusterizer");
-    Info("Exec", "took %f seconds for Clusterizing %f seconds per event", 
-        gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nevents ) ;
+    printf("Exec took %f seconds for Clusterizing %f seconds per event", 
+        gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents );
   }
-  
+
 }
 
 //____________________________________________________________________________
-Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
+Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
                                    Int_t nPar, Float_t * fitparameters) const
 { 
   // Calls TMinuit to fit the energy distribution of a cluster with several maxima 
   // The initial values for fitting procedure are set equal to the positions of local maxima.
   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
-  TClonesArray * digits = gime->Digits() ; 
-  
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+  TClonesArray *digits = emcalLoader->Digits();
 
   gMinuit->mncler();                     // Reset Minuit's list of paramters
   gMinuit->SetPrintLevel(-1) ;           // No Printout
@@ -210,16 +295,15 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDig
 
   Int_t iDigit ;
 
-  AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
-
   for(iDigit = 0; iDigit < nDigits; iDigit++){
     digit = maxAt[iDigit]; 
 
-    Int_t relid[4] ;
     Float_t x = 0.;
     Float_t z = 0.;
-    geom->AbsToRelNumbering(digit->GetId(), relid) ;
-    geom->PosInAlice(relid, x, z) ;
+    //   have to be tune for TRD1; May 31,06
+    //   Int_t relid[2] ;
+    //   fGeom->AbsToRelNumbering(digit->GetId(), relid) ; // obsolete method
+    //   fGeom->PosInAlice(relid, x, z) ;                  // obsolete method
 
     Float_t energy = maxAtEnergy[iDigit] ;
 
@@ -274,14 +358,33 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDig
 //____________________________________________________________________________
 void AliEMCALClusterizerv1::GetCalibrationParameters() 
 {
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
-  const AliEMCALDigitizer * dig = gime->Digitizer(BranchName()) ;
-
-  fADCchannelTower   = dig->GetTowerchannel() ;
-  fADCpedestalTower  = dig->GetTowerpedestal();
-
-  fADCchannelPreSho  = dig->GetPreShochannel() ;
-  fADCpedestalPreSho = dig->GetPreShopedestal() ; 
+  // Set calibration parameters:
+  // if calibration database exists, they are read from database,
+  // otherwise, they are taken from digitizer.
+  //
+  // It is a user responsilibity to open CDB before reconstruction, 
+  // for example: 
+  // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
+
+  //Check if calibration is stored in data base
+   if(AliCDBManager::Instance()->IsDefaultStorageSet()){
+     AliCDBEntry *entry = (AliCDBEntry*) 
+     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
+     if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
+   }
+   if(!fCalibData)
+     {
+       //If calibration is not available use default parameters
+       //Loader
+       AliEMCALLoader *emcalLoader = 
+        dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+       if ( !emcalLoader->Digitizer() ) 
+        emcalLoader->LoadDigitizer();
+       AliEMCALDigitizer * dig = dynamic_cast<AliEMCALDigitizer*>(emcalLoader->Digitizer());
+       
+       fADCchannelECA   = dig->GetECAchannel() ;
+       fADCpedestalECA  = dig->GetECApedestal();
+     }
 }
 
 //____________________________________________________________________________
@@ -290,331 +393,346 @@ void AliEMCALClusterizerv1::Init()
   // Make all memory allocations which can not be done in default constructor.
   // Attach the Clusterizer task to the list of EMCAL tasks
   
-  if ( strcmp(GetTitle(), "") == 0 )
-    SetTitle("galice.root") ;
+  AliRunLoader *rl = AliRunLoader::GetRunLoader();
+  fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
+  fGeom->GetTransformationForSM(); // Global <-> Local
+  AliInfo(Form("geom 0x%x",fGeom));
 
-  TString branchname = GetName() ;
-  branchname.Remove(branchname.Index(Version())-1) ;
-
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), branchname.Data(), fToSplit ) ; 
-  if ( gime == 0 ) {
-    Error("Init", "Could not obtain the Getter object !" ) ; 
-    return ;
-  } 
-
-  fSplitFile = 0 ;
-  if(fToSplit){
-    // construct the name of the file as /path/EMCAL.SDigits.root
-    //First - extract full path if necessary
-    TString fileName(GetTitle()) ;
-    Ssiz_t islash = fileName.Last('/') ;
-    if(islash<fileName.Length())
-      fileName.Remove(islash+1,fileName.Length()) ;
-    else
-      fileName="" ;
-    // Next - append the file name 
-    fileName+="EMCAL.RecData." ;
-    if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
-      fileName+=branchname ;
-      fileName+="." ;
-    }
-    fileName+="root" ;
-    // Finally - check if the file already opened or open the file
-    fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));   
-    if(!fSplitFile)
-      fSplitFile =  TFile::Open(fileName.Data(),"update") ;
-  }
-
-  const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
-  fNTowers = geom->GetNZ() *  geom->GetNPhi() ;
   if(!gMinuit) 
     gMinuit = new TMinuit(100) ;
 
-  gime->PostClusterizer(this) ;
-  gime->PostRecPoints(branchname ) ;
+  fHists = BookHists();
 }
 
 //____________________________________________________________________________
 void AliEMCALClusterizerv1::InitParameters()
-{
-  fNumberOfPreShoClusters = fNumberOfTowerClusters = 0 ;   
-  fPreShoClusteringThreshold  = 0.0001;
-  fTowerClusteringThreshold   = 0.2;    
-  fTowerLocMaxCut  = 0.03 ;
-  fPreShoLocMaxCut = 0.03 ;
-  
-  fW0     = 4.5 ;
-  fW0CPV  = 4.0 ;
+{ 
+  // Initializes the parameters for the Clusterizer
+  fNumberOfECAClusters = 0;
+
+  fNTowerInGroup = 36;  //Produces maximum of 80 pseudoclusters per event
+
+  fECAClusteringThreshold   = 0.5;  // value obtained from Alexei
+  fECALocMaxCut = 0.03; // ??
 
+  fECAW0     = 4.5 ;
   fTimeGate = 1.e-8 ; 
-  
   fToUnfold = kFALSE ;
-  
-  TString clusterizerName( GetName()) ; 
-  if (clusterizerName.IsNull() ) 
-    clusterizerName = "Default" ; 
-  clusterizerName.Append(":") ; 
-  clusterizerName.Append(Version()) ; 
-  SetName(clusterizerName) ;
-  fRecPointsInRun          = 0 ; 
+  fRecPointsInRun  = 0 ;
+  fMinECut = 0.01; // have to be tune
 
+  fCalibData               = 0 ;
 }
 
 //____________________________________________________________________________
-Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
+Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2) const
 {
   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching 
   //                                       = 1 are neighbour
-  //                                       = 2 are not neighbour but do not continue searching
+  //                                       = 2 is in different SM, quit from searching
   // neighbours are defined as digits having at least a common vertex 
   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
   //                                      which is compared to a digit (d2)  not yet in a cluster  
 
-   AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ;
+  static Int_t rv; 
+  static Int_t nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
+  static Int_t nSupMod2=0, nTower2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
+  static Int_t rowdiff, coldiff;
+  rv = 0 ; 
 
-  Int_t rv = 0 ; 
+  fGeom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
+  fGeom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2);
+  if(nSupMod1 != nSupMod2) return 2; // different SM
 
-  Int_t relid1[4] ; 
-  geom->AbsToRelNumbering(d1->GetId(), relid1) ; 
+  fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, iphi1,ieta1);
+  fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nTower2,nIphi2,nIeta2, iphi2,ieta2);
 
-  Int_t relid2[4] ; 
-  geom->AbsToRelNumbering(d2->GetId(), relid2) ; 
+  rowdiff = TMath::Abs(iphi1 - iphi2);  
+  coldiff = TMath::Abs(ieta1 - ieta2) ;  
+  
+  if (( coldiff <= 1 )  && ( rowdiff <= 1 )) rv = 1;  // neighbours with at least commom vertex
  
-  if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same EMCAL Arm 
-    Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
-    Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
-    
-    if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
-      if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate))
-      rv = 1 ; 
-    }
-    else {
-      if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) 
-       rv = 2; //  Difference in row numbers is too large to look further 
-    }
-
-  } 
-  else {
-    
-    if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )  
-      rv=2 ;
-  }
+  if (gDebug == 2 && rv==1) 
+  printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n", 
+        rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);   
+  
   return rv ; 
 }
 
 //____________________________________________________________________________
-Bool_t AliEMCALClusterizerv1::IsInTower(AliEMCALDigit * digit) const
+Int_t AliEMCALClusterizerv1::AreInGroup(AliEMCALDigit * d1, AliEMCALDigit * d2) const
 {
-  // Tells if (true) or not (false) the digit is in a EMCAL-Tower 
-  Bool_t rv = kFALSE ; 
-  if (!digit->IsInPreShower()) 
-    rv = kTRUE; 
-  return rv ; 
+  // Tells whether two digits fall within the same supermodule and
+  // tower grouping.  The number of towers in a group is controlled by
+  // the parameter nTowersInGroup
+  //    = 0 are not in same group but continue searching 
+  //    = 1 same group
+  //    = 2 is in different SM, quit from searching
+  //    = 3 different tower group, quit from searching
+  //
+  // The order of d1 and d2 is important: first (d1) should be a digit 
+  // already in a cluster which is compared to a digit (d2)  not yet in a cluster  
+
+  //JLK Question: does the quit from searching assume that the digits
+  //are ordered, so that once you are in a different SM, you'll not
+  //find another in the list that will match?  How about my TowerGroup search?
+
+  static Int_t rv;
+  static Int_t nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
+  static Int_t nSupMod2=0, nTower2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
+  static Int_t towerGroup1 = -1, towerGroup2 = -1;
+  rv = 0 ;
+
+  fGeom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
+  fGeom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2);
+  if(nSupMod1 != nSupMod2) return 2; // different SM
+
+  static Int_t nTowerInSM = fGeom->GetNCellsInSupMod()/fGeom->GetNCellsInTower();
+
+  //figure out which tower grouping each digit belongs to
+  for(int it = 0; it < nTowerInSM/fNTowerInGroup; it++) {
+    if(nTower1 <= nTowerInSM - it*fNTowerInGroup) towerGroup1 = it;
+    if(nTower2 <= nTowerInSM - it*fNTowerInGroup) towerGroup2 = it;
+  }
+  if(towerGroup1 != towerGroup2) return 3; //different Towergroup
+
+  //same SM, same towergroup, we're happy
+  if(towerGroup1 == towerGroup2 && towerGroup2 >= 0)
+    rv = 1;
+
+  if (gDebug == 2 && rv==1)
+  printf("AreInGroup: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n",
+         rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);
+
+  return rv ;
 }
 
 //____________________________________________________________________________
-Bool_t AliEMCALClusterizerv1::IsInPreShower(AliEMCALDigit * digit) const
+void AliEMCALClusterizerv1::Unload() 
 {
-  // Tells if (true) or not (false) the digit is in a EMCAL-PreShower
-  Bool_t rv = kFALSE ; 
-  if (digit->IsInPreShower()) 
-    rv = kTRUE; 
-  return rv ; 
+  // Unloads the Digits and RecPoints
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+    
+  emcalLoader->UnloadDigits() ; 
+  emcalLoader->UnloadRecPoints() ; 
 }
-
 //____________________________________________________________________________
-void AliEMCALClusterizerv1::WriteRecPoints(Int_t event)
+void AliEMCALClusterizerv1::WriteRecPoints()
 {
 
   // Creates new branches with given title
   // fills and writes into TreeR.
 
-  AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; 
-  TObjArray * towerRecPoints = gime->TowerRecPoints() ; 
-  TObjArray * preshoRecPoints = gime->PreShowerRecPoints() ; 
-  TClonesArray * digits = gime->Digits() ; 
-  TTree * treeR ; 
-  
-  if(fToSplit){
-    if(!fSplitFile)
-      return ;
-    fSplitFile->cd() ;
-    TString name("TreeR") ;
-    name += event ; 
-    treeR = dynamic_cast<TTree*>(fSplitFile->Get(name)); 
-  }
-  else{
-    treeR = gAlice->TreeR();
-  }
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+
+  TObjArray * aECARecPoints = emcalLoader->RecPoints() ; 
 
-  if(!treeR){
-    gAlice->MakeTree("R", fSplitFile);
-    treeR = gAlice->TreeR() ;
+  TClonesArray * digits = emcalLoader->Digits() ; 
+  TTree * treeR = emcalLoader->TreeR();  
+  if ( treeR==0 ) {
+    emcalLoader->MakeRecPointsContainer();
+    treeR = emcalLoader->TreeR();  
   }
   Int_t index ;
-  //Evaluate position, dispersion and other RecPoint properties...
-  for(index = 0; index < towerRecPoints->GetEntries(); index++)
-    (dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(index)))->EvalAll(fW0,digits) ;
-  
-  towerRecPoints->Sort() ;
-
-  for(index = 0; index < towerRecPoints->GetEntries(); index++)
-    (dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(index)))->SetIndexInList(index) ;
 
-  towerRecPoints->Expand(towerRecPoints->GetEntriesFast()) ; 
-  //Now the same for pre shower
-  for(index = 0; index < preshoRecPoints->GetEntries(); index++)
-    (dynamic_cast<AliEMCALRecPoint *>(preshoRecPoints->At(index)))->EvalAll(fW0CPV,digits)  ;
-  preshoRecPoints->Sort() ;
+  //Evaluate position, dispersion and other RecPoint properties for EC section
+  for(index = 0; index < aECARecPoints->GetEntries(); index++)
+    (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ;
+  
+  aECARecPoints->Sort() ;
 
-  for(index = 0; index < preshoRecPoints->GetEntries(); index++)
-    (dynamic_cast<AliEMCALRecPoint *>(preshoRecPoints->At(index)))->SetIndexInList(index) ;
+  for(index = 0; index < aECARecPoints->GetEntries(); index++) {
+    (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->SetIndexInList(index) ;
+    (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->Print();
+  }
 
-  preshoRecPoints->Expand(preshoRecPoints->GetEntriesFast()) ;
-  
   Int_t bufferSize = 32000 ;    
   Int_t splitlevel = 0 ; 
 
-  //First Tower branch
-  TBranch * towerBranch = treeR->Branch("EMCALTowerRP","TObjArray",&towerRecPoints,bufferSize,splitlevel);
-  towerBranch->SetTitle(BranchName());
+  //EC section branch
   
-  //Now Pre Shower branch 
-  TBranch * preshoBranch = treeR->Branch("EMCALPreShoRP","TObjArray",&preshoRecPoints,bufferSize,splitlevel);
-  preshoBranch->SetTitle(BranchName());
-    
-  //And Finally  clusterizer branch
-  AliEMCALClusterizerv1 * cl = (AliEMCALClusterizerv1*)gime->Clusterizer(BranchName()) ;
-  TBranch * clusterizerBranch = treeR->Branch("AliEMCALClusterizer","AliEMCALClusterizerv1",
-                                             &cl,bufferSize,splitlevel);
-  clusterizerBranch->SetTitle(BranchName());
-
-  towerBranch        ->Fill() ;
-  preshoBranch        ->Fill() ;
-  clusterizerBranch->Fill() ;
-
-  treeR->AutoSave() ; //Write(0,kOverwrite) ;  
-  if(gAlice->TreeR()!=treeR)
-    treeR->Delete(); 
+  TBranch * branchECA = 0;
+  if ((branchECA = treeR->GetBranch("EMCALECARP")))
+    branchECA->SetAddress(&aECARecPoints);
+  else
+    treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel);
+  treeR->Fill() ;
+
+  emcalLoader->WriteRecPoints("OVERWRITE");
+
 }
 
 //____________________________________________________________________________
-void AliEMCALClusterizerv1::MakeClusters()
+void AliEMCALClusterizerv1::MakeClusters(char* option)
 {
   // Steering method to construct the clusters stored in a list of Reconstructed Points
   // A cluster is defined as a list of neighbour digits
+
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+
+  TObjArray * aECARecPoints = emcalLoader->RecPoints() ; 
     
-  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
-  
-  TObjArray * towerRecPoints  = gime->TowerRecPoints(BranchName()) ; 
-  TObjArray * preshoRecPoints = gime->PreShowerRecPoints(BranchName()) ; 
-  towerRecPoints->Delete() ;
-  preshoRecPoints->Delete() ;
-  
-  TClonesArray * digits = gime->Digits() ; 
-  if ( !digits ) {
-    Fatal("MakeClusters -> Digits with name %s not found", GetName() ) ; 
-  } 
-  TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone()) ;
-  
-  
-  // Clusterization starts  
-  
-  TIter nextdigit(digitsC) ; 
-  AliEMCALDigit * digit ; 
-  Bool_t notremoved = kTRUE ;
-  
-  while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
-    AliEMCALRecPoint * clu = 0 ; 
-    
-    TArrayI clusterdigitslist(1500) ;   
-    Int_t index ;
-    if (( IsInTower (digit)  && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fTowerClusteringThreshold  ) || 
-        ( IsInPreShower (digit) && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fPreShoClusteringThreshold  ) ) {
-      
-      Int_t iDigitInCluster = 0 ; 
-      
-      if  ( IsInTower(digit) ) {   
-       // start a new Tower RecPoint
-       if(fNumberOfTowerClusters >= towerRecPoints->GetSize()) 
-         towerRecPoints->Expand(2*fNumberOfTowerClusters+1) ;
-       
-       towerRecPoints->AddAt(new  AliEMCALTowerRecPoint(""), fNumberOfTowerClusters) ;
-       clu = dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(fNumberOfTowerClusters)) ; 
-       fNumberOfTowerClusters++ ; 
-       clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower())) ; 
-       clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;  
-       iDigitInCluster++ ; 
-       digitsC->Remove(digit) ; 
-       
-      } else { 
-       
-       // start a new Pre Shower cluster
-       if(fNumberOfPreShoClusters >= preshoRecPoints->GetSize()) 
-         preshoRecPoints->Expand(2*fNumberOfPreShoClusters+1);
-       
-       preshoRecPoints->AddAt(new AliEMCALTowerRecPoint(""), fNumberOfPreShoClusters) ;
-       
-       clu =  dynamic_cast<AliEMCALTowerRecPoint *>(preshoRecPoints->At(fNumberOfPreShoClusters))  ;  
-       fNumberOfPreShoClusters++ ; 
-       clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower() ) );     
-       clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ; 
-       iDigitInCluster++ ; 
+  if (fGeom==0) 
+     AliFatal("Did not get geometry from EMCALLoader");
+
+  aECARecPoints->Clear();
+
+  //Start with pseudoclusters, if option
+  if(strstr(option,"pseudo")) {
+    TClonesArray * digits  = emcalLoader->Digits() ; 
+    TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone());
+
+    TIter nextdigit(digitsC) ;
+    AliEMCALDigit * digit;
+
+    AliEMCALRecPoint * recPoint = 0 ; 
+    int ineb=0;
+    nextdigit.Reset();
+
+    // PseudoClusterization starts    
+    while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
+      recPoint = 0 ; 
+      TArrayI clusterECAdigitslist(1000); // what is this   
+
+      if(fGeom->CheckAbsCellId(digit->GetId()) ) { //Is this an EMCAL digit? Just maing sure...
+       Int_t iDigitInECACluster = 0; // start a new Tower RecPoint
+       if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
+       AliEMCALRecPoint * rp = new  AliEMCALRecPoint("") ; 
+       aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
+       recPoint = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ; 
+       fNumberOfECAClusters++ ; 
+
+       recPoint->SetClusterType(AliESDCaloCluster::kPseudoCluster);
+
+       recPoint->AddDigit(*digit, digit->GetAmp()) ; 
+       clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;    
+       iDigitInECACluster++ ; 
        digitsC->Remove(digit) ; 
-       nextdigit.Reset() ;
-       
-       // Here we remove remaining Tower digits, which cannot make a cluster
-       
-       if( notremoved ) { 
-         while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
-           if( IsInTower(digit) )
-             digitsC->Remove(digit) ;
-           else 
-             break ; 
-         }
-         notremoved = kFALSE ;
-       }
-       
-      } // else        
+       AliDebug(1,Form("MakePseudoClusters: OK id = %d, adc = %f \n", digit->GetId(), digit->GetAmp()));
+       nextdigit.Reset(); // will start from beggining
       
-      nextdigit.Reset() ;
+       AliEMCALDigit * digitN = 0; // digi neighbor
+       Int_t index = 0 ;
+
+       // Find the neighbours
+       while (index < iDigitInECACluster){ // scan over digits already in cluster 
+         digit =  (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]);
+         index++ ; 
+         while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits 
+           ineb = AreInGroup(digit, digitN);       // call (digit,digitN) in THAT oder !!!!! 
+           switch (ineb ) {
+              case 0 :   // not a neighbour
+              break ;
+             case 1 :   // are neighbours 
+              recPoint->AddDigit(*digitN, digitN->GetAmp() ) ;
+              clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; 
+              iDigitInECACluster++ ; 
+              digitsC->Remove(digitN) ;
+              break ;
+             case 2 :   // different SM
+              break ;
+              case 3 :   // different Tower group
+              break ;
+           } // switch
+         } // scan over the reduced list of digits
+       } // scan over digits already in cluster
+       nextdigit.Reset() ;  // will start from beggining
+      }
+    }
+    if(recPoint)
+      AliDebug(1,Form("MakeClusters: cl.e %f \n", recPoint->GetEnergy())); 
+    //if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl; 
+    delete digitsC ;
+  }
+
+  //Now do real clusters
+  TClonesArray * digits  = emcalLoader->Digits() ; 
+  TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone());
+
+  TIter nextdigit(digitsC) ;
+  AliEMCALDigit * digit;
+
+  AliEMCALRecPoint * recPoint = 0 ; 
+  int ineb=0;
+  nextdigit.Reset();
+
+  double e=0.0, ehs = 0.0;
+  while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // clean up digits
+    e = Calibrate(digit->GetAmp(), digit->GetId());
+    AliEMCALHistoUtilities::FillH1(fHists, 10, digit->GetAmp());
+    AliEMCALHistoUtilities::FillH1(fHists, 11, e);
+    if(e < fMinECut ) digitsC->Remove(digit);
+    else              ehs += e;
+  } 
+  AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %d\n",
+                 digits->GetEntries(),fMinECut,ehs));
+  //cout << " Number of digits " << digits->GetEntries() << " -> (e>" <<fMinECut <<")";
+  //cout << digitsC->GetEntries()<< " ehs "<<ehs<<endl; 
+
+  // Clusterization starts    
+  //  cout << "Outer Loop" << endl;
+  ineb=0;
+  nextdigit.Reset();
+  recPoint = 0 ; 
+  while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
+    recPoint = 0 ; 
+    TArrayI clusterECAdigitslist(1000); // what is this   
+
+    if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp(), digit->GetId()) > fECAClusteringThreshold  ) ){
+      Int_t iDigitInECACluster = 0; // start a new Tower RecPoint
+      if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
+      AliEMCALRecPoint * rp = new  AliEMCALRecPoint("") ; 
+      aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
+      recPoint = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ; 
+      fNumberOfECAClusters++ ; 
+
+      recPoint->SetClusterType(AliESDCaloCluster::kClusterv1);
+
+      recPoint->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())) ; 
+      clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;     
+      iDigitInECACluster++ ; 
+      digitsC->Remove(digit) ; 
+      AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(),
+      Calibrate(digit->GetAmp(),digit->GetId()), fECAClusteringThreshold));  
+      nextdigit.Reset(); // will start from beggining
       
-      AliEMCALDigit * digitN ; 
-      index = 0 ;
-      while (index < iDigitInCluster){ // scan over digits already in cluster 
-       digit =  (AliEMCALDigit*)digits->At(clusterdigitslist[index])  ;      
+      AliEMCALDigit * digitN = 0; // digi neighbor
+      Int_t index = 0 ;
+
+      // Find the neighbours
+      while (index < iDigitInECACluster){ // scan over digits already in cluster 
+       digit =  (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]);
        index++ ; 
-        while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits 
-         Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
-         switch (ineb ) {
-          case 0 :   // not a neighbour
-           break ;
-         case 1 :   // are neighbours 
-           clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->IsInPreShower() ) ) ;
-           clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; 
-           iDigitInCluster++ ; 
-           digitsC->Remove(digitN) ;
-           break ;
-          case 2 :   // too far from each other
-           goto endofloop;   
+        while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits 
+         ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!! 
+          switch (ineb ) {
+            case 0 :   // not a neighbour
+             break ;
+           case 1 :   // are neighbours 
+             recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmp(),digitN->GetId()) ) ;
+             clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; 
+             iDigitInECACluster++ ; 
+             digitsC->Remove(digitN) ;
+              break ;
+             case 2 :   // different SM
+              break ;
          } // switch
-         
-       } // while digitN
-       
-      endofloop: ;
-       nextdigit.Reset() ; 
-      } // loop over cluster     
-    } // energy theshold     
-  } // while digit  
+        } // scan over the reduced list of digits
+      } // scan over digits already in cluster
+      nextdigit.Reset() ;  // will start from beggining
+    }
+  } // while digit 
+  if(recPoint)
+    AliDebug(1,Form("MakeClusters: cl.e %f \n", recPoint->GetEnergy())); 
+  //if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl; 
   delete digitsC ;
+
+  AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,digits->GetEntriesFast())); 
 }
 
 //____________________________________________________________________________
-void AliEMCALClusterizerv1::MakeUnfolding()
+void AliEMCALClusterizerv1::MakeUnfolding() const
 {
   Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
  
@@ -633,10 +751,10 @@ Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t r)
 }
 
 //____________________________________________________________________________
-void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower
-                                                Int_t nMax
-                                                AliEMCALDigit ** maxAt
-                                                Float_t * maxAtEnergy)
+void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/
+                                          Int_t /*nMax*/
+                                          AliEMCALDigit ** /*maxAt*/
+                                          Float_t * /*maxAtEnergy*/) const
 {
   // Performs the unfolding of a cluster with nMax overlapping showers 
   
@@ -645,7 +763,9 @@ void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower,
 }
 
 //_____________________________________________________________________________
-void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
+void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * /*Grad*/,
+                                              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
@@ -653,7 +773,7 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Do
   ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
 }
 //____________________________________________________________________________
-void AliEMCALClusterizerv1::Print(Option_t * option)const
+void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const
 {
   // Print clusterizer parameters
 
@@ -663,122 +783,128 @@ void AliEMCALClusterizerv1::Print(Option_t * option)const
     
     // Print parameters
  
-    TString taskName(GetName()) ; 
+    TString taskName(GetName()) ;
     taskName.ReplaceAll(Version(), "") ;
     
-    message += "--------------- " ; 
-    message += taskName.Data() ; 
-    message += " " ; 
-    message += GetTitle() ; 
-    message += "-----------\n" ;  
-    message += "Clusterizing digits from the file: " ; 
-    message += taskName.Data() ;  
-    message += "\n                           Branch: " ; 
-    message += GetName() ;  
-    message += "\n                       EMC Clustering threshold = " ; 
-    message += fTowerClusteringThreshold ; 
-    message += "\n                       EMC Local Maximum cut    = " ;
-    message += fTowerLocMaxCut ; 
-    message += "\n                       EMC Logarothmic weight   = " ;
-    message += fW0 ;
-    message += "\n                       Pre Shower Clustering threshold = " ; 
-    message += fPreShoClusteringThreshold ;
-    message += "\n                       Pre Shower  Local Maximum cut    = " ;
-    message += fPreShoLocMaxCut ;
-    message += "\n                       Pre Shower Logarothmic weight   = " ; 
-    message += fW0CPV ;
+    printf("--------------- "); 
+    printf(taskName.Data()) ; 
+    printf(" "); 
+    printf(GetTitle()) ; 
+    printf("-----------\n");  
+    printf("Clusterizing digits from the file: "); 
+    printf(taskName.Data());  
+    printf("\n                           Branch: "); 
+    printf(GetName()); 
+    printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut); 
+    printf("\n                       ECA Logarithmic weight   = %f", fECAW0); 
     if(fToUnfold)
-      message +="\nUnfolding on\n" ;
+      printf("\nUnfolding on\n");
     else
-      message += "\nUnfolding off\n";
+      printf("\nUnfolding off\n");
     
-    message += "------------------------------------------------------------------" 
+    printf("------------------------------------------------------------------")
   }
   else
-    message += "AliEMCALClusterizerv1 not initialized " ;
-  
-  Info("Print", message.Data() ) ; 
+    printf("AliEMCALClusterizerv1 not initialized ") ;
 }
 
 //____________________________________________________________________________
 void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
 {
   // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
+  TObjArray * aECARecPoints = emcalLoader->RecPoints() ; 
+    
+  printf("PrintRecPoints: Clusterization result:") ; 
+  
+  printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ;
+  printf("           Found %d ECA Rec Points\n ", 
+        aECARecPoints->GetEntriesFast()) ; 
 
-  TObjArray * towerRecPoints = AliEMCALGetter::GetInstance()->TowerRecPoints() ; 
-  TObjArray * preshoRecPoints = AliEMCALGetter::GetInstance()->PreShowerRecPoints() ; 
-
-  TString message("\n")  ;
-
-  message += "event " ; 
-  message += gAlice->GetEvNumber() ;
-  message += "\n       Found " ; 
-  message += towerRecPoints->GetEntriesFast() ; 
-  message += " TOWER Rec Points and " ;
-  message += preshoRecPoints->GetEntriesFast() ; 
-  message += " PRE SHOWER RecPoints\n" ;
-
-  fRecPointsInRun +=  towerRecPoints->GetEntriesFast() ; 
-  fRecPointsInRun +=  preshoRecPoints->GetEntriesFast() ; 
-
-  char * tempo = new char[8192]; 
+  fRecPointsInRun +=  aECARecPoints->GetEntriesFast() ; 
   
   if(strstr(option,"all")) {
-
-    message += "Tower clusters\n" ;
-    message += "Index    Ene(MeV) Multi Module     phi     r   theta  Lambda 1 Lambda 2  # of prim  Primaries list\n" ;      
-    
-    Int_t index ;
-    for (index = 0 ; index < towerRecPoints->GetEntries() ; index++) {
-      AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(towerRecPoints->At(index)) ; 
+    Int_t index =0;
+    printf("\n-----------------------------------------------------------------------\n") ;
+    printf("Clusters in ECAL section\n") ;
+    printf("Index    Ene(GeV) Multi Module     phi     r   theta    X    Y      Z   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;      
+   Float_t maxE=0; 
+   Float_t maxL1=0; 
+   Float_t maxL2=0; 
+   Float_t maxDis=0; 
+
+    for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) {
+      AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(aECARecPoints->At(index)) ; 
       TVector3  globalpos;  
-      rp->GetGlobalPosition(globalpos);
+      //rp->GetGlobalPosition(globalpos);
+      TVector3  localpos;  
+      rp->GetLocalPosition(localpos);
       Float_t lambda[2]; 
       rp->GetElipsAxis(lambda);
       Int_t * primaries; 
       Int_t nprimaries;
       primaries = rp->GetPrimaries(nprimaries);
-      sprintf(tempo, "\n%6d  %8.2f  %3d     %2d     %4.1f    %4.1f %4.1f %4f  %4f    %2d     : ", 
-             rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(), 
-             globalpos.X(), globalpos.Y(), globalpos.Z(), lambda[0], lambda[1], nprimaries) ; 
-      message += tempo ; 
-     
+      printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ", 
+            rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
+            globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), 
+            rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; 
+  /////////////
+      if(rp->GetEnergy()>maxE){
+             maxE=rp->GetEnergy();
+             maxL1=lambda[0];
+             maxL2=lambda[1];
+             maxDis=rp->GetDispersion();
+      }
+      fPointE->Fill(rp->GetEnergy());
+      fPointL1->Fill(lambda[0]);
+      fPointL2->Fill(lambda[1]);
+      fPointDis->Fill(rp->GetDispersion());
+      fPointMult->Fill(rp->GetMultiplicity());
+      ///////////// 
       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
-       sprintf(tempo, "%d ", primaries[iprimary] ) ; 
-       message += tempo ;
+       printf("%d ", primaries[iprimary] ) ; 
       } 
     }
 
-    //Now plot Pre shower recPoints
-
-    message += "\n-----------------------------------------------------------------------\n" ;
+      fMaxE->Fill(maxE);
+      fMaxL1->Fill(maxL1);
+      fMaxL2->Fill(maxL2);
+      fMaxDis->Fill(maxDis);
 
-    message += "PreShower clusters\n" ;
-    message += "Index    Ene(MeV) Multi Module     phi     r   theta  Lambda 1 Lambda 2  # of prim  Primaries list\n" ;      
-    
-    for (index = 0 ; index < preshoRecPoints->GetEntries() ; index++) {
-      AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint *>(preshoRecPoints->At(index)) ; 
-      TVector3  globalpos;  
-      rp->GetGlobalPosition(globalpos);
-      Float_t lambda[2]; 
-      rp->GetElipsAxis(lambda);
-      Int_t * primaries;
-      Int_t nprimaries;
-      primaries = rp->GetPrimaries(nprimaries);
-      sprintf(tempo, "\n%6d  %8.2f  %3d     %2d     %4.1f    %4.1f %4.1f %4f  %4f    %2d     : ", 
-             rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(), 
-             globalpos.X(), globalpos.Y(), globalpos.Z(), lambda[0], lambda[1], nprimaries) ; 
-      message += tempo ; 
-  
-      for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
-       sprintf(tempo, "%d ", primaries[iprimary] ) ; 
-       message += tempo ;
-      }         
-    }
 
-    message += "\n-----------------------------------------------------------------------" ;
+    printf("\n-----------------------------------------------------------------------\n");
   }
-  delete tempo ; 
-  Info("PrintRecPoints", message.Data() ) ; 
-  
+}
+TList* AliEMCALClusterizerv1::BookHists()
+{
+  //set up histograms for monitoring clusterizer performance
+
+  gROOT->cd();
+
+       fPointE = new TH1F("pointE","point energy", 2000, 0.0, 150.);
+       fPointL1 = new TH1F("pointL1","point L1", 1000, 0.0, 3.);
+       fPointL2 = new TH1F("pointL2","point L2", 1000, 0.0, 3.);
+       fPointDis = new TH1F("pointDis","point Dis", 1000, 0.0, 3.);
+       fPointMult = new TH1F("pointMult","point Mult", 100, 0.0, 100.);
+       fDigitAmp = new TH1F("digitAmp","Digit Amplitude", 2000, 0.0, 5000.);
+       fMaxE = new TH1F("maxE","Max point energy", 2000, 0.0, 150.);
+       fMaxL1 = new TH1F("maxL1","Max point L1", 1000, 0.0, 3.);
+       fMaxL2 = new TH1F("maxL2","Max point L2", 1000, 0.0, 3.);
+       fMaxDis = new TH1F("maxDis","Max point Dis", 1000, 0.0, 3.); // 9
+       //
+        new TH1F("adcOfDigits","adc of digits(threshold control)", 1001, -0.5, 1000.5);   // 10
+        new TH1F("energyOfDigits","energy of digits(threshold control)", 1000, 0.0, 1.);  // 11
+
+  return AliEMCALHistoUtilities::MoveHistsToList("EmcalClusterizerv1ControlHists", kFALSE);
+}
+
+void AliEMCALClusterizerv1::SaveHists(const char *fn)
+{
+  AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE);
+}
+
+void AliEMCALClusterizerv1::Browse(TBrowser* b)
+{
+  if(fHists) b->Add(fHists);
+  TTask::Browse(b);
 }