]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSClusterizerv1.cxx
Add digitization of amplitudes, baseline and noise
[u/mrichter/AliRoot.git] / PHOS / AliPHOSClusterizerv1.cxx
index 0a02d61a7ffbc8c9bb9d53c4b08719c5a8262eac..06ee60461440e7ff62550724c9121d1a1e0babdb 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:
-     AliPHOSClusterizerv1()
-     CPV clusterizing parameters added
-
-     MakeClusters()
-     After first PPSD digit remove EMC digits only once
-*/
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.99  2006/11/07 16:49:51  kharlov
+ * Corrections for next event switch in case of raw data (B.Polichtchouk)
+ *
+ * Revision 1.98  2006/10/27 17:14:27  kharlov
+ * Introduce AliDebug and AliLog (B.Polichtchouk)
+ *
+ * Revision 1.97  2006/08/29 11:41:19  kharlov
+ * Missing implementation of ctors and = operator are added
+ *
+ * Revision 1.96  2006/08/25 16:56:30  kharlov
+ * Compliance with Effective C++
+ *
+ * Revision 1.95  2006/08/11 12:36:26  cvetan
+ * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
+ *
+ * Revision 1.94  2006/08/07 12:27:49  hristov
+ * Removing obsolete code which affected the event numbering scheme
+ *
+ * Revision 1.93  2006/08/01 12:20:17  cvetan
+ * 1. Adding a possibility to read and reconstruct an old rcu formatted raw data. This is controlled by an option of AliReconstruction and AliPHOSReconstructor. 2. In case of raw data processing (without galice.root) create the default AliPHOSGeometry object. Most likely this should be moved to the CDB
+ *
+ * Revision 1.92  2006/04/29 20:26:46  hristov
+ * Separate EMC and CPV calibration (Yu.Kharlov)
+ *
+ * Revision 1.91  2006/04/22 10:30:17  hristov
+ * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
+ *
+ * Revision 1.90  2006/04/11 15:22:59  hristov
+ * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
+ *
+ * Revision 1.89  2006/03/13 14:05:42  kharlov
+ * Calibration objects for EMC and CPV
+ *
+ * Revision 1.88  2006/01/11 08:54:52  hristov
+ * Additional protection in case no calibration entry was found
+ *
+ * Revision 1.87  2005/11/22 08:46:43  kharlov
+ * Updated to new CDB (Boris Polichtchouk)
+ *
+ * Revision 1.86  2005/11/14 21:52:43  hristov
+ * Coding conventions
+ *
+ * Revision 1.85  2005/09/27 16:08:08  hristov
+ * New version of CDB storage framework (A.Colla)
+ *
+ * Revision 1.84  2005/09/21 10:02:47  kharlov
+ * Reading calibration from CDB (Boris Polichtchouk)
+ *
+ * Revision 1.82  2005/09/02 15:43:13  kharlov
+ * Add comments in GetCalibrationParameters and Calibrate
+ *
+ * Revision 1.81  2005/09/02 14:32:07  kharlov
+ * Calibration of raw data
+ *
+ * Revision 1.80  2005/08/24 15:31:36  kharlov
+ * Setting raw digits flag
+ *
+ * Revision 1.79  2005/07/25 15:53:53  kharlov
+ * Read raw data
+ *
+ * Revision 1.78  2005/05/28 14:19:04  schutz
+ * Compilation warnings fixed by T.P.
+ *
+ */
+
 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
 //////////////////////////////////////////////////////////////////////////////
 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and 
-//  unfolding of the clusters with several local maxima.  
-//  results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
+//  unfolds the clusters having several local maxima.  
+//  Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
 //  PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all 
-//  parameters including input digits branch file name thresholds etc.)
-//  This TTask normally called from Reconstructioner, but as well can be used it in 
-//  standalone mode:
-// root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root")  
-// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
-//               //reads gAlice from header file "..."                      
-// root [1] cl->ExecuteTask()  
+//  parameters including input digits branch title, thresholds etc.)
+//  This TTask is normally called from Reconstructioner, but can as well be used in 
+//  standalone mode.
+// Use Case:
+//  root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname")  
+//  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
+//               // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default")
+//               // and saves recpoints in branch named "recpointsname" (default = "digitsname")                       
+//  root [1] cl->ExecuteTask()  
 //               //finds RecPoints in all events stored in galice.root
-// root [2] cl->SetDigitsBranch("digits2") 
-//               //sets another input file
-// root [3] cl->SetRecPointsBranch("recp2")  
-//               //sets another aouput file
-// root [4] cl->SetEmcLocalMaxCut(0.03)  
+//  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->SetEmcLocalMaxCut(0.03)  
 //               //set clusterization parameters
-// root [5] cl->ExecuteTask("deb all time")  
+//  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 
 
 // --- 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"
 
 // --- Standard library ---
 
-#include <iostream.h>
-#include <iomanip.h>
-
 // --- AliRoot header files ---
-
+#include "AliLog.h"
+#include "AliRunLoader.h"
+#include "AliPHOSGetter.h"
+#include "AliPHOSGeometry.h" 
 #include "AliPHOSClusterizerv1.h"
+#include "AliPHOSEmcRecPoint.h"
 #include "AliPHOSCpvRecPoint.h"
 #include "AliPHOSDigit.h"
 #include "AliPHOSDigitizer.h"
-#include "AliPHOSEmcRecPoint.h"
-#include "AliPHOS.h"
-#include "AliPHOSPpsdRecPoint.h"
-#include "AliRun.h"
+#include "AliPHOSCalibrationDB.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
+#include "AliCDBEntry.h"
 
 ClassImp(AliPHOSClusterizerv1)
-
+  
 //____________________________________________________________________________
-  AliPHOSClusterizerv1::AliPHOSClusterizerv1():AliPHOSClusterizer()
+AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
+  AliPHOSClusterizer(),
+  fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
+  fWrite(0),                  fNumberOfEmcClusters(0),  fNumberOfCpvClusters(0),
+  fCalibData(0),              fADCchanelEmc(0),         fADCpedestalEmc(0),
+  fADCchanelCpv(0),           fADCpedestalCpv(0),       fEmcClusteringThreshold(0),
+  fCpvClusteringThreshold(0), fEmcMinE(0),              fCpvMinE(0),
+  fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
+  fW0CPV(0),                  fRecPointsInRun(0),       fEmcTimeGate(0),
+  fIsOldRCUFormat(0)
 {
-  // default ctor (to be used)
-  SetName("AliPHOSClusterizer");
-  SetTitle("Version 1") ;
-
-  fNumberOfCpvClusters     = 0 ; 
-  fNumberOfEmcClusters     = 0 ; 
-    
-  fCpvClusteringThreshold  = 0.0;
-  fEmcClusteringThreshold  = 0.2;   
-  fPpsdClusteringThreshold = 0.0000002 ;
-  
-  fEmcLocMaxCut            = 0.03 ;
-  fCpvLocMaxCut            = 0.03 ;
+  // default ctor (to be used mainly by Streamer)
   
-  fW0                      = 4.5 ;
-  fW0CPV                   = 4.0 ;
+  InitParameters() ; 
+  fDefaultInit = kTRUE ; 
+}
 
-  fGeom  = 0 ;
+//____________________________________________________________________________
+AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName) :
+  AliPHOSClusterizer(alirunFileName, eventFolderName),
+  fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
+  fWrite(0),                  fNumberOfEmcClusters(0),  fNumberOfCpvClusters(0),
+  fCalibData(0),              fADCchanelEmc(0),         fADCpedestalEmc(0),
+  fADCchanelCpv(0),           fADCpedestalCpv(0),       fEmcClusteringThreshold(0),
+  fCpvClusteringThreshold(0), fEmcMinE(0),              fCpvMinE(0),
+  fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
+  fW0CPV(0),                  fRecPointsInRun(0),       fEmcTimeGate(0),
+  fIsOldRCUFormat(0)
+{
+  // ctor with the indication of the file where header Tree and digits Tree are stored
   
-  fDigits = 0 ;
-  fDigitizer = 0 ;
-  fEmcRecPoints = 0 ;
-  fCpvRecPoints = 0 ;
+  InitParameters() ;
+  Init() ;
+  fDefaultInit = kFALSE ; 
+}
 
-  fIsInitialized = kFALSE ;
-  
+//____________________________________________________________________________
+AliPHOSClusterizerv1::AliPHOSClusterizerv1(const AliPHOSClusterizerv1 & obj) :
+  AliPHOSClusterizer(obj),
+  fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
+  fWrite(0),                  fNumberOfEmcClusters(0),  fNumberOfCpvClusters(0),
+  fCalibData(0),              fADCchanelEmc(0),         fADCpedestalEmc(0),
+  fADCchanelCpv(0),           fADCpedestalCpv(0),       fEmcClusteringThreshold(0),
+  fCpvClusteringThreshold(0), fEmcMinE(0),              fCpvMinE(0),
+  fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
+  fW0CPV(0),                  fRecPointsInRun(0),       fEmcTimeGate(0),
+  fIsOldRCUFormat(0)
+{
+  // Copy constructor
 }
 //____________________________________________________________________________
-  AliPHOSClusterizerv1::AliPHOSClusterizerv1(const char* headerFile,const char* digitsFile):AliPHOSClusterizer()
+  AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
 {
-  SetName("AliPHOSClusterizer");
-  SetTitle("Version 1") ;
-  
-  fNumberOfCpvClusters     = 0 ; 
-  fNumberOfEmcClusters     = 0 ; 
-    
-  fCpvClusteringThreshold  = 0.0;
-  fEmcClusteringThreshold  = 0.2;   
-  fPpsdClusteringThreshold = 0.0000002 ;
-  
-  fEmcLocMaxCut            = 0.03 ;
-  fCpvLocMaxCut            = 0.03 ;
-  
-  fW0                      = 4.5 ;
-  fW0CPV                   = 4.0 ;
-  
-  fToUnfold = kTRUE ;
-  
-  fHeaderFileName = headerFile  ;
-  fDigitsBranchTitle = digitsFile ;
-  
-  TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ;
-  
-  if(file == 0){
-    file = new TFile(fHeaderFileName.Data(),"update") ;
-    gAlice = (AliRun *) file->Get("gAlice") ;
-  }
-  
-  AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ;    
-  fGeom  = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() );
-  
-  fDigits = new TClonesArray("AliPHOSDigit",10) ;
-  fDigitizer = new AliPHOSDigitizer() ;
-  fEmcRecPoints = new TObjArray(200) ;
-  fCpvRecPoints = new TObjArray(200) ;
-  
-  if(!gMinuit) gMinuit = new TMinuit(100) ;
-  
-  // add Task to //root/Tasks folder
-  TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
-  roottasks->Add(this) ; 
-
-
-  fIsInitialized = kTRUE ;
+  // dtor
 
 }
 //____________________________________________________________________________
-void AliPHOSClusterizerv1::Exec(Option_t * option){
-  // Steerign function
+const TString AliPHOSClusterizerv1::BranchName() const 
+{  
+  return GetName();
+}
+//____________________________________________________________________________
+Float_t  AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
+{  
+  // Convert EMC measured amplitude into real energy.
+  // Calibration parameters are taken from calibration data base for raw data,
+  // or from digitizer parameters for simulated data.
+
+  if(fCalibData){
+    Int_t relId[4];
+    AliPHOSGetter *gime = AliPHOSGetter::Instance();
+    gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
+    Int_t   module = relId[0];
+    Int_t   column = relId[3];
+    Int_t   row    = relId[2];
+    if(absId <= fEmcCrystals) { // this is EMC 
+      fADCchanelEmc   = fCalibData->GetADCchannelEmc (module,column,row);
+      return amp*fADCchanelEmc ;        
+    }
+  }
+  else{ //simulation
+    if(absId <= fEmcCrystals) // this is EMC 
+      return fADCpedestalEmc + amp*fADCchanelEmc ;        
+  }
+  return 0;
+}
+
+//____________________________________________________________________________
+Float_t  AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
+{  
+  // Convert digitized CPV amplitude into charge.
+  // Calibration parameters are taken from calibration data base for raw data,
+  // or from digitizer parameters for simulated data.
+
+  if(fCalibData){
+    Int_t relId[4];
+    AliPHOSGetter *gime = AliPHOSGetter::Instance();
+    gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
+    Int_t   module = relId[0];
+    Int_t   column = relId[3];
+    Int_t   row    = relId[2];
+    if(absId > fEmcCrystals) { // this is CPV
+      fADCchanelCpv   = fCalibData->GetADCchannelCpv (module,column,row);
+      fADCpedestalCpv = fCalibData->GetADCpedestalCpv(module,column,row);
+      return fADCpedestalCpv + amp*fADCchanelCpv ;              
+    }     
+  }
+  else{ //simulation
+    if(absId > fEmcCrystals) // this is CPV
+      return fADCpedestalCpv+ amp*fADCchanelCpv ;       
+  }
+  return 0;
+}
 
-  if(!fIsInitialized) Init() ;
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::Exec(Option_t *option)
+{
+  // 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("PHOSClusterizer");  
+    gBenchmark->Start("PHOSClusterizer"); 
   
-  Int_t nEvents = (Int_t) gAlice->TreeE()->GetEntries() ;
+  if(strstr(option,"print")) {
+    Print() ; 
+    return ;
+  }
+
+  GetCalibrationParameters() ;
+
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  if (fRawReader == 0)
+    gime->SetRawDigits(kFALSE);
+  else
+    gime->SetRawDigits(kTRUE);
   
-  for(fEvent = 0;fEvent< nEvents; fEvent++){
-    if(!ReadDigits())  //reads digits for event fEvent
-      return;
+  if (fLastEvent == -1) 
+    fLastEvent = gime->MaxEvent() - 1 ;
+  else 
+    fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // one event at the time 
+  Int_t nEvents   = fLastEvent - fFirstEvent + 1;
+
+  Int_t ievent ;
+  for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
+    if (fRawReader == 0)
+      gime->Event(ievent    ,"D"); // Read digits from simulated data
+    else {
+      AliRunLoader * rl = AliRunLoader::GetRunLoader(gime->PhosLoader()->GetTitle());
+      rl->GetEvent(ievent);
+      gime->Event(fRawReader,"W",fIsOldRCUFormat); // Read digits from raw data
+    }
+    fNumberOfEmcClusters  = fNumberOfCpvClusters  = 0 ;
+    
     MakeClusters() ;
     
-    if(fToUnfold) MakeUnfolding() ;
-    WriteRecPoints() ;
-    if(strstr(option,"deb"))
+    AliDebug(2,Form(" ---- Printing clusters (%d) of event %d.\n",
+          gime->EmcRecPoints()->GetEntries(),ievent));
+    if(AliLog::GetGlobalDebugLevel()>1)
+      gime->EmcRecPoints()->Print();
+
+    if(fToUnfold)             
+      MakeUnfolding() ;
+
+    WriteRecPoints();
+
+    if(strstr(option,"deb"))  
       PrintRecPoints(option) ;
+
+    //increment the total number of recpoints per run 
+    fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;  
+    fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;  
   }
   
+  if(fWrite) //do not unload in "on flight" mode
+    Unload();
+  
   if(strstr(option,"tim")){
     gBenchmark->Stop("PHOSClusterizer");
-    cout << "AliPHOSClusterizer:" << endl ;
-    cout << "  took " << gBenchmark->GetCpuTime("PHOSClusterizer") << " seconds for Clusterizing " 
-        <<  gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents << " seconds per event " << endl ;
-    cout << endl ;
-  }
-  
-  
+    AliInfo(Form("took %f seconds for Clusterizing %f seconds per event \n",
+        gBenchmark->GetCpuTime("PHOSClusterizer"), 
+        gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents )) ; 
+  } 
 }
 
 //____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy,
-                                   Int_t nPar, Float_t * fitparameters)
+Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** 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
+
+  
+  AliPHOSGetter * gime = AliPHOSGetter::Instance();
+  TClonesArray * digits = gime->Digits(); 
+  
 
   gMinuit->mncler();                     // Reset Minuit's list of paramters
   gMinuit->SetPrintLevel(-1) ;           // No Printout
@@ -207,7 +345,7 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Fl
 
   TList * toMinuit = new TList();
   toMinuit->AddAt(emcRP,0) ;
-  toMinuit->AddAt(fDigits,1) ;
+  toMinuit->AddAt(digits,1) ;
   
   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 
@@ -220,34 +358,35 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Fl
 
   Int_t iDigit ;
 
+  const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
 
   for(iDigit = 0; iDigit < nDigits; iDigit++){
-    digit = (AliPHOSDigit *) maxAt[iDigit]; 
+    digit = maxAt[iDigit]; 
 
     Int_t relid[4] ;
-    Float_t x ;
-    Float_t z ;
-    fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
-    fGeom->RelPosInModule(relid, x, z) ;
+    Float_t x = 0.;
+    Float_t z = 0.;
+    geom->AbsToRelNumbering(digit->GetId(), relid) ;
+    geom->RelPosInModule(relid, x, z) ;
 
     Float_t energy = maxAtEnergy[iDigit] ;
 
     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
     index++ ;   
     if(ierflg != 0){ 
-      cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : x = " << x << endl ;
+      Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
       return kFALSE;
     }
     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
     index++ ;   
     if(ierflg != 0){
-      cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : z = " << z << endl ;
+       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
       return kFALSE;
     }
     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
     index++ ;   
     if(ierflg != 0){
-      cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : energy = " << energy << endl ;      
+      Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;      
       return kFALSE;
     }
   }
@@ -265,7 +404,7 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Fl
   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize 
 
   if(ierflg == 4){  // Minimum not found   
-    cout << "PHOS Unfolding>  Fit not converged, cluster abandoned "<< endl ;      
+    Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );      
     return kFALSE ;
   }            
   for(index = 0; index < nPar; index++){
@@ -281,110 +420,168 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Fl
 }
 
 //____________________________________________________________________________
-void AliPHOSClusterizerv1::Init(){
-
-  if(!fIsInitialized){
-    if(fHeaderFileName.IsNull())
-      fHeaderFileName = "galice.root" ;
+void AliPHOSClusterizerv1::GetCalibrationParameters() 
+{
+  // 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");
+
+  AliPHOSGetter * gime = AliPHOSGetter::Instance();
+  // fCalibData = new AliPHOSCalibData(gAlice->GetRunNumber()); //original
+  fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
+  
+  if(!fCalibData)
+    {
+      if ( !gime->Digitizer() ) 
+       gime->LoadDigitizer();
+      AliPHOSDigitizer * dig = gime->Digitizer(); 
+      fADCchanelEmc   = dig->GetEMCchannel() ;
+      fADCpedestalEmc = dig->GetEMCpedestal();
     
-
-    TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ;
-
-    if(file == 0){
-      file = new TFile(fHeaderFileName.Data(),"update") ;
-      gAlice = (AliRun *) file->Get("gAlice") ;
+      fADCchanelCpv   = dig->GetCPVchannel() ;
+      fADCpedestalCpv = dig->GetCPVpedestal() ; 
     }
+}
+
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::Init()
+{
+  // Make all memory allocations which can not be done in default constructor.
+  // Attach the Clusterizer task to the list of PHOS tasks
+  AliPHOSGetter* gime = AliPHOSGetter::Instance() ;
+  if(!gime)
+    gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
 
-    AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ;    
-    fGeom  = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() );
+  AliPHOSGeometry * geom = gime->PHOSGeometry();
 
-    fDigits = new TClonesArray("AliPHOSDigit",10) ;
-    fDigitizer = new AliPHOSDigitizer() ;
-    fEmcRecPoints = new TObjArray(200) ;
-    fCpvRecPoints = new TObjArray(200) ;
-     
-    if(!gMinuit) gMinuit = new TMinuit(100) ;
+  fEmcCrystals = geom->GetNModules() *  geom->GetNCristalsInModule() ;
 
-    // add Task to //root/Tasks folder
-    TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
-    roottasks->Add(this) ; 
+  if(!gMinuit) 
+    gMinuit = new TMinuit(100);
 
-    fIsInitialized = kTRUE ;
+  if ( !gime->Clusterizer() ) {
+    gime->PostClusterizer(this);
   }
 }
+
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::InitParameters()
+{
+
+  fNumberOfCpvClusters     = 0 ; 
+  fNumberOfEmcClusters     = 0 ; 
+  
+  fCpvClusteringThreshold  = 0.0;
+  fEmcClusteringThreshold  = 0.2;   
+  
+  fEmcLocMaxCut            = 0.03 ;
+  fCpvLocMaxCut            = 0.03 ;
+
+  fEmcMinE                 = 0.01 ;
+  fCpvMinE                 = 0.0  ;
+
+  fW0                      = 4.5 ;
+  fW0CPV                   = 4.0 ;
+
+  fEmcTimeGate             = 1.e-8 ; 
+  
+  fToUnfold                = kTRUE ;
+    
+  fRecPointsInRun          = 0 ;
+
+  fWrite                   = kTRUE ;
+
+  fCalibData               = 0 ;
+
+  SetEventRange(0,-1) ;
+
+  fIsOldRCUFormat          = kFALSE;
+}
+
 //____________________________________________________________________________
 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * 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
-  // neighbours are defined as digits having at least common vertex 
+  // neighbours are defined as digits having at least 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  
 
+  AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
+
   Int_t rv = 0 ; 
 
   Int_t relid1[4] ; 
-  fGeom->AbsToRelNumbering(d1->GetId(), relid1) ; 
+  geom->AbsToRelNumbering(d1->GetId(), relid1) ; 
 
   Int_t relid2[4] ; 
-  fGeom->AbsToRelNumbering(d2->GetId(), relid2) ; 
+  geom->AbsToRelNumbering(d2->GetId(), relid2) ; 
  
-  if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module and the same PPSD Module 
+  if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module 
     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() ) < fEmcTimeGate))
       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 
+        rv = 2; //  Difference in row numbers is too large to look further 
     }
 
   } 
   else {
     
-    if( (relid1[0] < relid2[0]) || (relid1[1] < relid2[1]) )  
+    if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )  
       rv=2 ;
 
   }
 
-  //Do NOT clusterize upper PPSD  
-  if( IsInPpsd(d1) && IsInPpsd(d2) &&
-     relid1[1] > 0                 &&
-     relid1[1] < fGeom->GetNumberOfPadsPhi()*fGeom->GetNumberOfPadsPhi() ) rv = 2 ;
-
   return rv ; 
 }
-
-
 //____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
+void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
 {
-  // Tells if (true) or not (false) the digit is in a PHOS-EMC module
-  Bool_t rv = kFALSE ; 
-
-  Int_t relid[4] ; 
-  fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
+  // Remove digits with amplitudes below threshold
 
-  if ( relid[1] == 0  ) rv = kTRUE; 
+  for(Int_t i=0; i<digits->GetEntriesFast(); i++){
+    AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
+    if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
+        (IsInCpv(digit) && CalibrateCPV(digit->GetAmp()   ,digit->GetId()) < fCpvMinE) )
+      digits->RemoveAt(i) ;
+  }
+  digits->Compress() ;
+  for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) { 
+    AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ; 
+    digit->SetIndexInList(i) ;     
+  }
 
-  return rv ; 
+  //Overwrite digits tree
+  AliPHOSGetter* gime = AliPHOSGetter::Instance();
+  TTree * treeD = gime->TreeD();
+  treeD->Branch("PHOS", &digits);
+  treeD->Fill() ;
+  gime->WriteDigits("OVERWRITE");
+  gime->PhosLoader()->UnloadDigits() ;
 }
-
 //____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::IsInPpsd(AliPHOSDigit * digit) const
+Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
 {
-  // Tells if (true) or not (false) the digit is in a PHOS-PPSD module
+  // Tells if (true) or not (false) the digit is in a PHOS-EMC module
  
   Bool_t rv = kFALSE ; 
+  AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
 
-  Int_t relid[4] ; 
-  fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
+  Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();  
 
-  if ( relid[1] > 0 && relid[0] > fGeom->GetNCPVModules() ) rv = kTRUE; 
+  if(digit->GetId() <= nEMC )   rv = kTRUE; 
 
   return rv ; 
 }
@@ -395,206 +592,90 @@ Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
  
   Bool_t rv = kFALSE ; 
+  
+  AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
 
-  Int_t relid[4] ; 
-  fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
+  Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
 
-  if ( relid[1] > 0 && relid[0] <= fGeom->GetNCPVModules() ) rv = kTRUE; 
+  if(digit->GetId() > nEMC )   rv = kTRUE;
 
   return rv ; 
 }
-//____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::ReadDigits(){
 
-  fNumberOfEmcClusters  = 0 ;
-  fNumberOfCpvClusters  = 0 ;
-
-  // Get Digits Tree header from file
-  char treeName[20]; 
-  sprintf(treeName,"TreeD%d",fEvent);
-  gAlice->GetEvent(fEvent) ;
-  gAlice->SetEvent(fEvent) ;
-
-  TTree * treeD = gAlice->TreeD()  ; // (TTree*)file->Get(treeName);
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::Unload() 
+{
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  gime->PhosLoader()->UnloadDigits() ; 
+  gime->PhosLoader()->UnloadRecPoints() ; 
+}
 
-  if(treeD==0){
-    cout << "Error in AliPHOSClusterizerv1 : no "<<treeName << endl  ;
-    cout << "    Do nothing " << endl ;
-    return kFALSE ;
-  }
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::WriteRecPoints()
+{
 
-  TBranch * digitsBranch = 0;
-  TBranch * digitizerBranch = 0;
+  // Creates new branches with given title
+  // fills and writes into TreeR.
 
-  TObjArray * branches = treeD->GetListOfBranches() ;
-  Int_t ibranch;
-  Bool_t phosNotFound = kTRUE ;
-  Bool_t digitizerNotFound = kTRUE ;
-  
-  for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
+  AliPHOSGetter * gime = AliPHOSGetter::Instance();
 
-    if(phosNotFound){
-      digitsBranch=(TBranch *) branches->At(ibranch) ;
-      if( fDigitsBranchTitle.CompareTo(digitsBranch->GetTitle())==0 )
-       if( strcmp(digitsBranch->GetName(),"PHOS") == 0)
-         phosNotFound = kFALSE ;
-    }
-    
-    if(digitizerNotFound){
-      digitizerBranch = (TBranch *) branches->At(ibranch) ;
-      if( fDigitsBranchTitle.CompareTo(digitizerBranch->GetTitle()) == 0)
-       if( strcmp(digitizerBranch->GetName(),"AliPHOSDigitizer") == 0) 
-         digitizerNotFound = kFALSE ;
+  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
+  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
+  TClonesArray * digits = gime->Digits() ; 
+  Int_t index ;
+  //Evaluate position, dispersion and other RecPoint properties..
+  Int_t nEmc = emcRecPoints->GetEntriesFast();
+  for(index = 0; index < nEmc; index++){
+    AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
+    rp->Purify(fEmcMinE) ;
+    if(rp->GetMultiplicity()>0.) //If this RP is not empty
+      rp->EvalAll(fW0,digits) ;
+    else{
+      emcRecPoints->RemoveAt(index) ;
+      delete rp ;
     }
-    
   }
-  
-  if(digitizerNotFound || phosNotFound){
-    cout << "ERROR in AliPHOSClusterizerv1: " << endl ;
-    cout << "      Can't find Branch with digits or Digitizer "<< endl ; ;
-    cout << "      Do nothing" <<endl  ;
-    return kFALSE ;
+  emcRecPoints->Compress() ;
+  emcRecPoints->Sort() ;
+  //  emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
+  for(index = 0; index < emcRecPoints->GetEntries(); index++){
+    dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
   }
   
-  digitsBranch->SetAddress(&fDigits) ;
-  digitizerBranch->SetAddress(&fDigitizer) ;
-  
-  treeD->GetEvent(0) ;
-  
-  fPedestal = fDigitizer->GetPedestal() ;
-  fSlope    = fDigitizer->GetSlope() ;
-  return kTRUE ;
-}
-
-//____________________________________________________________________________
-void AliPHOSClusterizerv1::WriteRecPoints(){
-  
-  Int_t index ;
-  //Evaluate poisition, dispersion and other RecPoint properties...
-  for(index = 0; index < fEmcRecPoints->GetEntries(); index++)
-    ((AliPHOSEmcRecPoint *)fEmcRecPoints->At(index))->EvalAll(fW0,fDigits) ;
-
-  fEmcRecPoints->Sort() ;
-
-  for(index = 0; index < fEmcRecPoints->GetEntries(); index++)
-    ((AliPHOSEmcRecPoint *)fEmcRecPoints->At(index))->SetIndexInList(index) ;
-
   //Now the same for CPV
-  for(index = 0; index < fCpvRecPoints->GetEntries(); index++)
-    ((AliPHOSRecPoint *)fCpvRecPoints->At(index))->EvalAll(fW0CPV,fDigits)  ;
-
-  fCpvRecPoints->Sort() ;
-
-  for(index = 0; index < fCpvRecPoints->GetEntries(); index++)
-    ((AliPHOSRecPoint *)fCpvRecPoints->At(index))->SetIndexInList(index) ;
-
-  if(gAlice->TreeR()==0)
-    gAlice->MakeTree("R") ;
+  for(index = 0; index < cpvRecPoints->GetEntries(); index++){
+    AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
+    rp->EvalAll(fW0CPV,digits) ;
+  }
+  cpvRecPoints->Sort() ;
   
-
-  //Check, if branches already exist
-  TBranch * emcBranch = 0;
-  TBranch * cpvBranch = 0;
-  TBranch * clusterizerBranch = 0;
+  for(index = 0; index < cpvRecPoints->GetEntries(); index++)
+    dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
   
-  TObjArray * branches = gAlice->TreeR()->GetListOfBranches() ;
-  Int_t ibranch;
-  Bool_t emcNotFound = kTRUE ;
-  Bool_t cpvNotFound = kTRUE ;  
-  Bool_t clusterizerNotFound = kTRUE ;
+  cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
   
-  for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
+  if(fWrite){ //We write TreeR
+    TTree * treeR = gime->TreeR();
     
-    if(emcNotFound){
-      emcBranch=(TBranch *) branches->At(ibranch) ;
-      if( fRecPointsBranchTitle.CompareTo(emcBranch->GetTitle())==0 ){
-       if( strcmp(emcBranch->GetName(),"PHOSEmcRP") == 0) {
-         emcNotFound = kFALSE ;
-       }
-      }
-    }
-    if(cpvNotFound){
-      cpvBranch=(TBranch *) branches->At(ibranch) ;
-      if( fRecPointsBranchTitle.CompareTo(cpvBranch->GetTitle())==0 ){
-       if( strcmp(cpvBranch->GetName(),"PHOSCpvRP") == 0) {
-         cpvNotFound = kFALSE ;
-       }
-      }
-    }
-
-    if(clusterizerNotFound){
-      clusterizerBranch = (TBranch *) branches->At(ibranch) ;
-      if( fRecPointsBranchTitle.CompareTo(clusterizerBranch->GetTitle()) == 0){
-       if( strcmp(clusterizerBranch->GetName(),"AliPHOSClusterizer") == 0) {
-         clusterizerNotFound = kFALSE ;
-       }
-      }
-    }
+    Int_t bufferSize = 32000 ;    
+    Int_t splitlevel = 0 ;
     
-  }
-  
-  if(!(clusterizerNotFound && emcNotFound && cpvNotFound)){
-    cout << "AliPHOSClusterizer error" << endl;
-    cout << "       Branches PHOSEmcRP, PHOSCpvRP and AliPHOSClusterizer " << endl ;
-    cout << "       with title '" << fRecPointsBranchTitle.Data() <<"' already exist" << endl ;
-    cout << "       can not overwrite " << endl ;
-    return ;
-  }
-
-  //Make branches in TreeR for RecPoints and Clusterizer
-  char * filename = 0;
-  if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){   //generating file name
-    filename = new char[strlen(gAlice->GetBaseFile())+20] ;
-    sprintf(filename,"%s/PHOS.Reco.root",gAlice->GetBaseFile()) ;
-  }
-  
-  //Make new branches
-  TDirectory *cwd = gDirectory;
-  
-  //First EMC
-  Int_t bufferSize = 32000 ;    
-  Int_t splitlevel = 0 ;
-  emcBranch = gAlice->TreeR()->Branch("PHOSEmcRP","TObjArray",&fEmcRecPoints,bufferSize,splitlevel);
-  emcBranch->SetTitle(fRecPointsBranchTitle.Data());
-  if (filename) {
-    emcBranch->SetFile(filename);
-    TIter next( emcBranch->GetListOfBranches());
-    while ((emcBranch=(TBranch*)next())) {
-      emcBranch->SetFile(filename);
-    }   
-    cwd->cd();
-  }
+    //First EMC
+    TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
+    emcBranch->SetTitle(BranchName());
     
-  //Now CPV branch
-  cpvBranch = gAlice->TreeR()->Branch("PHOSCpvRP","TObjArray",&fCpvRecPoints,bufferSize,splitlevel);
-  cpvBranch->SetTitle(fRecPointsBranchTitle.Data());
-  if (filename) {
-    cpvBranch->SetFile(filename);
-    TIter next( cpvBranch->GetListOfBranches());
-    while ((cpvBranch=(TBranch*)next())) {
-      cpvBranch->SetFile(filename);
-    }   
-    cwd->cd();
-  } 
+    //Now CPV branch
+    TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
+    cpvBranch->SetTitle(BranchName());
     
-  //And Finally  clusterizer branch
-  AliPHOSClusterizerv1 * cl = this ;
-  clusterizerBranch = gAlice->TreeR()->Branch("AliPHOSClusterizer","AliPHOSClusterizerv1",
-                                             &cl,bufferSize,splitlevel);
-  clusterizerBranch->SetTitle(fRecPointsBranchTitle.Data());
-  if (filename) {
-    clusterizerBranch->SetFile(filename);
-    TIter next( clusterizerBranch->GetListOfBranches());
-    while ((clusterizerBranch=(TBranch*)next())) {
-      clusterizerBranch->SetFile(filename);
-    }   
-    cwd->cd();
+    emcBranch        ->Fill() ;
+    cpvBranch        ->Fill() ;
+    
+    gime->WriteRecPoints("OVERWRITE");
+    gime->WriteClusterizer("OVERWRITE");
   }
-  
-  gAlice->TreeR()->Fill() ;
-  
-  gAlice->TreeR()->Write(0,kOverwrite) ;  
-  
 }
 
 //____________________________________________________________________________
@@ -602,69 +683,84 @@ void AliPHOSClusterizerv1::MakeClusters()
 {
   // Steering method to construct the clusters stored in a list of Reconstructed Points
   // A cluster is defined as a list of neighbour digits
-  fEmcRecPoints->Clear() ;
-  fCpvRecPoints->Clear() ;
+
+  
+  AliPHOSGetter * gime = AliPHOSGetter::Instance();
+
+  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
+  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
+
+  emcRecPoints->Delete() ;
+  cpvRecPoints->Delete() ;
   
+  TClonesArray * digits = gime->Digits() ; 
 
+  //Remove digits below threshold
+  CleanDigits(digits) ;
+
+  TClonesArray * digitsC =  static_cast<TClonesArray*>( digits->Clone() ) ;
   // Clusterization starts  
-  TClonesArray * digits =  (TClonesArray*)fDigits->Clone() ;
 
-  TIter nextdigit(digits) ; 
+  TIter nextdigit(digitsC) ; 
   AliPHOSDigit * digit ; 
   Bool_t notremoved = kTRUE ;
 
-  while ( (digit = (AliPHOSDigit *)nextdigit()) ) { // scan over the list of digits
+  while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
+    
+
     AliPHOSRecPoint * clu = 0 ; 
 
-    TArrayI clusterdigitslist(1000) ;   
+    TArrayI clusterdigitslist(1500) ;   
     Int_t index ;
 
-    if (( IsInEmc (digit) && Calibrate(digit->GetAmp()) > fEmcClusteringThreshold  ) || 
-        ( IsInPpsd(digit) && Calibrate(digit->GetAmp()) > fPpsdClusteringThreshold ) ||
-        ( IsInCpv (digit) && Calibrate(digit->GetAmp()) > fCpvClusteringThreshold  ) ) {
-      
+    if (( IsInEmc (digit) &&
+         CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) || 
+        ( IsInCpv (digit) &&
+         CalibrateCPV(digit->GetAmp()   ,digit->GetId()) > fCpvClusteringThreshold ) ) {
       Int_t iDigitInCluster = 0 ; 
-
+      
       if  ( IsInEmc(digit) ) {   
-       // start a new EMC RecPoint
-       if(fNumberOfEmcClusters >= fEmcRecPoints->GetSize()) fEmcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
-       fEmcRecPoints->AddAt(new  AliPHOSEmcRecPoint(), fNumberOfEmcClusters) ;
-       clu = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(fNumberOfEmcClusters) ; 
+        // start a new EMC RecPoint
+        if(fNumberOfEmcClusters >= emcRecPoints->GetSize()) 
+          emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
+          
+        emcRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
+        clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ; 
        fNumberOfEmcClusters++ ; 
-       clu->AddDigit(*digit, Calibrate(digit->GetAmp())) ; 
-       clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;  
-       iDigitInCluster++ ; 
-       digits->Remove(digit) ; 
+       clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
+        clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
+        iDigitInCluster++ ;
+        digitsC->Remove(digit) ; 
 
       } else { 
-       
-       // start a new PPSD/CPV cluster
-       if(fNumberOfCpvClusters >= fCpvRecPoints->GetSize()) fCpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
-       if(IsInPpsd(digit)) 
-         fCpvRecPoints->AddAt(new AliPHOSPpsdRecPoint(),fNumberOfCpvClusters) ;
-       else
-         fCpvRecPoints->AddAt(new AliPHOSCpvRecPoint(), fNumberOfCpvClusters) ;
-       clu =  (AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fNumberOfCpvClusters)  ;  
-       fNumberOfCpvClusters++ ; 
-
-       clu->AddDigit(*digit, Calibrate(digit->GetAmp()) ) ;    
-       clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ; 
-       iDigitInCluster++ ; 
-       digits->Remove(digit) ; 
+        
+        // start a new CPV cluster
+        if(fNumberOfCpvClusters >= cpvRecPoints->GetSize()) 
+          cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
+
+        cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
+
+        clu =  dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;  
+        fNumberOfCpvClusters++ ; 
+        clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;        
+        clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;        
+        iDigitInCluster++ ; 
+        digitsC->Remove(digit) ; 
         nextdigit.Reset() ;
-       
-       // Here we remove resting EMC digits, which cannot make cluster
-       
+        
+        // Here we remove remaining EMC digits, which cannot make a cluster
+        
         if( notremoved ) { 
-         while( ( digit = (AliPHOSDigit *)nextdigit() ) ) {
+          while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
             if( IsInEmc(digit) ) 
-             digits->Remove(digit) ;
+              digitsC->Remove(digit) ;
             else 
-             break ;
-         }
-         notremoved = kFALSE ;
-       }
-       
+              break ;
+          }
+          notremoved = kFALSE ;
+        }
+        
       } // else        
       
       nextdigit.Reset() ;
@@ -672,28 +768,32 @@ void AliPHOSClusterizerv1::MakeClusters()
       AliPHOSDigit * digitN ; 
       index = 0 ;
       while (index < iDigitInCluster){ // scan over digits already in cluster 
-       digit =  (AliPHOSDigit*)fDigits->At(clusterdigitslist[index])  ;      
-       index++ ; 
-        while ( (digitN = (AliPHOSDigit *)nextdigit()) ) { // scan over the reduced list of digits 
-         Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
+        digit =  dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) )  ;      
+        index++ ; 
+        while ( (digitN = dynamic_cast<AliPHOSDigit *>( 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() ) ) ;
-           clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; 
-           iDigitInCluster++ ; 
-           digits->Remove(digitN) ;
-           break ;
+            break ;
+          case 1 :   // are neighbours 
+           if (IsInEmc (digitN))
+             clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
+           else
+             clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp()   , digitN->GetId() ) );
+
+            clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; 
+            iDigitInCluster++ ; 
+            digitsC->Remove(digitN) ;
+            break ;
           case 2 :   // too far from each other
-           goto endofloop;   
-         } // switch
-         
-       } // while digitN
-       
+            goto endofloop;   
+          } // switch
+          
+        } // while digitN
+        
       endofloop: ;
-       nextdigit.Reset() ; 
-       
+        nextdigit.Reset() ; 
+        
       } // loop over cluster     
 
     } // energy theshold  
@@ -701,77 +801,89 @@ void AliPHOSClusterizerv1::MakeClusters()
     
   } // while digit
 
-  delete digits ;
+  delete digitsC ;
 
 }
 
 //____________________________________________________________________________
-void AliPHOSClusterizerv1::MakeUnfolding(){
-  //Unfolds clusters using the shape of ElectroMagnetic shower
-  // Performs unfolding of all EMC/CPV but NOT ppsd clusters
+void AliPHOSClusterizerv1::MakeUnfolding()
+{
+  // Unfolds clusters using the shape of an ElectroMagnetic shower
+  // Performs unfolding of all EMC/CPV clusters
+
+  AliPHOSGetter * gime = AliPHOSGetter::Instance();
+
+  const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 
-  //Unfold first EMC clusters 
+  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
+  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
+  TClonesArray * digits = gime->Digits() ; 
+  
+  // Unfold first EMC clusters 
   if(fNumberOfEmcClusters > 0){
 
-    Int_t nModulesToUnfold = fGeom->GetNModules() ; 
+    Int_t nModulesToUnfold = geom->GetNModules() ; 
 
     Int_t numberofNotUnfolded = fNumberOfEmcClusters ; 
     Int_t index ;   
     for(index = 0 ; index < numberofNotUnfolded ; index++){
       
-      AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(index) ;
+      AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ;
       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
-       break ;
+        break ;
       
       Int_t nMultipl = emcRecPoint->GetMultiplicity() ; 
-      Int_t * maxAt = new Int_t[nMultipl] ;
+      AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
-      Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigits) ;
+      Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
       
       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0       
-       UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
-       fEmcRecPoints->Remove(emcRecPoint); 
-       fEmcRecPoints->Compress() ;
-       index-- ;
-       fNumberOfEmcClusters -- ;
-       numberofNotUnfolded-- ;
+        UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
+        emcRecPoints->Remove(emcRecPoint); 
+        emcRecPoints->Compress() ;
+        index-- ;
+        fNumberOfEmcClusters -- ;
+        numberofNotUnfolded-- ;
+      }
+      else{
+        emcRecPoint->SetNExMax(1) ; //Only one local maximum
       }
       
       delete[] maxAt ; 
       delete[] maxAtEnergy ; 
     }
   } 
-  //Unfolding of EMC clusters finished
+  // Unfolding of EMC clusters finished
 
 
-  //Unfold now CPV clusters
+  // Unfold now CPV clusters
   if(fNumberOfCpvClusters > 0){
     
-    Int_t nModulesToUnfold = fGeom->GetNCPVModules() ;
+    Int_t nModulesToUnfold = geom->GetNModules() ;
 
     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;     
     Int_t index ;   
     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
       
-      AliPHOSRecPoint * recPoint = (AliPHOSRecPoint *) fCpvRecPoints->At(index) ;
+      AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ;
 
       if(recPoint->GetPHOSMod()> nModulesToUnfold)
-       break ;
+        break ;
       
-      AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint*) recPoint ; 
+      AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ; 
       
       Int_t nMultipl = emcRecPoint->GetMultiplicity() ; 
-      Int_t * maxAt = new Int_t[nMultipl] ;
+      AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
-      Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigits) ;
+      Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
       
       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0       
-       UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
-       fCpvRecPoints->Remove(emcRecPoint); 
-       fCpvRecPoints->Compress() ;
-       index-- ;
-       numberofCpvNotUnfolded-- ;
-       fNumberOfCpvClusters-- ;
+        UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
+        cpvRecPoints->Remove(emcRecPoint); 
+        cpvRecPoints->Compress() ;
+        index-- ;
+        numberofCpvNotUnfolded-- ;
+        fNumberOfCpvClusters-- ;
       }
       
       delete[] maxAt ; 
@@ -782,24 +894,11 @@ void AliPHOSClusterizerv1::MakeUnfolding(){
   
 }
 
-//____________________________________________________________________________
-void AliPHOSClusterizerv1::SetDigitsBranch(const char * title){
-  
-    fDigitsBranchTitle = title  ; 
-
-}
-//____________________________________________________________________________
-void AliPHOSClusterizerv1::SetRecPointsBranch(const char * title){
-  
-    fRecPointsBranchTitle = title;
-
-}
-
 //____________________________________________________________________________
 Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t r)
 { 
   // Shape of the shower (see PHOS TDR)
-  // If you change this function, change also the gradien evaluation  in ChiSquare()
+  // If you change this function, change also the gradient evaluation in ChiSquare()
 
   Double_t r4    = r*r*r*r ;
   Double_t r295  = TMath::Power(r, 2.95) ;
@@ -809,18 +908,27 @@ Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t r)
 
 //____________________________________________________________________________
 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, 
-                                                Int_t nMax, 
-                                                int * maxAt, 
-                                                Float_t * maxAtEnergy)
+                                          Int_t nMax, 
+                                          AliPHOSDigit ** maxAt, 
+                                          Float_t * maxAtEnergy)
 { 
   // Performs the unfolding of a cluster with nMax overlapping showers 
 
+  AliPHOSGetter * gime = AliPHOSGetter::Instance();
+
+  const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
+
+  const TClonesArray * digits = gime->Digits() ; 
+  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
+  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
+
   Int_t nPar = 3 * nMax ;
   Float_t * fitparameters = new Float_t[nPar] ;
 
   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
   if( !rv ) {
     // Fit failed, return and remove cluster
+    iniEmc->SetNExMax(-1) ;
     delete[] fitparameters ; 
     return ;
   }
@@ -831,18 +939,18 @@ void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 
   Int_t nDigits = iniEmc->GetMultiplicity() ;  
   Float_t * efit = new Float_t[nDigits] ;
-  Float_t xDigit,zDigit,distance ;
-  Float_t xpar,zpar,epar  ;
+  Float_t xDigit=0.,zDigit=0.,distance=0. ;
+  Float_t xpar=0.,zpar=0.,epar=0.  ;
   Int_t relid[4] ;
-  AliPHOSDigit * digit ;
+  AliPHOSDigit * digit = 0 ;
   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 
   Int_t iparam ;
   Int_t iDigit ;
   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
-    digit = (AliPHOSDigit*) fDigits->At(emcDigits[iDigit] ) ;   
-    fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
-    fGeom->RelPosInModule(relid, xDigit, zDigit) ;
+    digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ;   
+    geom->AbsToRelNumbering(digit->GetId(), relid) ;
+    geom->RelPosInModule(relid, xDigit, zDigit) ;
     efit[iDigit] = 0;
 
     iparam = 0 ;    
@@ -872,37 +980,38 @@ void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
     epar = fitparameters[iparam+2] ;
     iparam += 3 ;    
     
-    AliPHOSEmcRecPoint * emcRP ;  
+    AliPHOSEmcRecPoint * emcRP = 0 ;  
 
     if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
       
-      if(fNumberOfEmcClusters >= fEmcRecPoints->GetSize())
-       fEmcRecPoints->Expand(2*fNumberOfEmcClusters) ;
+      if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
+        emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
       
-      (*fEmcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint() ;
-      emcRP = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(fNumberOfEmcClusters);
+      (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
+      emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
       fNumberOfEmcClusters++ ;
+      emcRP->SetNExMax((Int_t)nPar/3) ;
     }
     else{//create new entries in fCpvRecPoints
-      if(fNumberOfCpvClusters >= fCpvRecPoints->GetSize())
-       fCpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
+      if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
+        cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
       
-      (*fCpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint() ;
-      emcRP = (AliPHOSEmcRecPoint *) fCpvRecPoints->At(fNumberOfCpvClusters);
+      (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
+      emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
       fNumberOfCpvClusters++ ;
     }
     
     Float_t eDigit ;
     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
-      digit = (AliPHOSDigit*) fDigits->At( emcDigits[iDigit] ) ; 
-      fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
-      fGeom->RelPosInModule(relid, xDigit, zDigit) ;
+      digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ; 
+      geom->AbsToRelNumbering(digit->GetId(), relid) ;
+      geom->RelPosInModule(relid, xDigit, zDigit) ;
       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
       distance =  TMath::Sqrt(distance) ;
       ratio = epar * ShowerShape(distance) / efit[iDigit] ; 
       eDigit = emcEnergies[iDigit] * ratio ;
       emcRP->AddDigit( *digit, eDigit ) ;
-    }  
+    }        
   }
  
   delete[] fitparameters ; 
@@ -913,27 +1022,25 @@ void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 //_____________________________________________________________________________
 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 {
-  // Calculates th Chi square for the cluster unfolding minimization
+  // Calculates the Chi square for the cluster unfolding minimization
   // Number of parameters, Gradient, Chi squared, parameters, what to do
 
+  TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
 
-  TList * toMinuit = (TList*) gMinuit->GetObjectFit() ;
-
-  AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint*) toMinuit->At(0)  ;
-  TClonesArray * digits = (TClonesArray*)toMinuit->At(1)  ;
+  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
+  TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
 
 
   
-  //  AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit
+  //  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
 
   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 
-  Int_t nOfDigits = emcRP->GetDigitsMultiplicity() ; 
+  Int_t nOdigits = emcRP->GetDigitsMultiplicity() ; 
 
   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
-
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-
+  
+  const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; 
   fret = 0. ;     
   Int_t iparam ;
 
@@ -946,9 +1053,9 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou
   AliPHOSDigit * digit ;
   Int_t iDigit ;
 
-  for( iDigit = 0 ; iDigit < nOfDigits ; iDigit++) {
+  for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 
-    digit = (AliPHOSDigit*) digits->At( emcDigits[iDigit] ) 
+    digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) )
 
     Int_t relid[4] ;
     Float_t xDigit ;
@@ -962,33 +1069,33 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou
        Int_t iParam = 0 ;
        efit = 0 ;
        while(iParam < nPar ){
-        Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
-        iParam++ ; 
-        distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ; 
-        distance = TMath::Sqrt( distance ) ; 
-        iParam++ ;      
-        efit += x[iParam] * ShowerShape(distance) ;
-        iParam++ ;
+         Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
+         iParam++ ; 
+         distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ; 
+         distance = TMath::Sqrt( distance ) ; 
+         iParam++ ;          
+         efit += x[iParam] * ShowerShape(distance) ;
+         iParam++ ;
        }
        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E) 
        iParam = 0 ;
        while(iParam < nPar ){
-        Double_t xpar = x[iParam] ;
-        Double_t zpar = x[iParam+1] ;
-        Double_t epar = x[iParam+2] ;
-        Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
-        Double_t shape = sum * ShowerShape(dr) ;
-        Double_t r4 = dr*dr*dr*dr ;
-        Double_t r295 = TMath::Power(dr,2.95) ;
-        Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
-                                        0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
-        
-        Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x    
-        iParam++ ; 
-        Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z         
-        iParam++ ; 
-        Grad[iParam] += shape ;                                  // Derivative over energy             
-        iParam++ ; 
+         Double_t xpar = x[iParam] ;
+         Double_t zpar = x[iParam+1] ;
+         Double_t epar = x[iParam+2] ;
+         Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
+         Double_t shape = sum * ShowerShape(dr) ;
+         Double_t r4 = dr*dr*dr*dr ;
+         Double_t r295 = TMath::Power(dr,2.95) ;
+         Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
+                                         0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
+         
+         Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x    
+         iParam++ ; 
+         Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z         
+         iParam++ ; 
+         Grad[iParam] += shape ;                                  // Derivative over energy             
+         iParam++ ; 
        }
      }
      efit = 0;
@@ -1011,135 +1118,104 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou
 }
 
 //____________________________________________________________________________
-void AliPHOSClusterizerv1::Print(Option_t * option)const
+void AliPHOSClusterizerv1::Print(const Option_t *)const
 {
-  if(fIsInitialized){
-    
+  // Print clusterizer parameters
+
+  TString message ; 
+  TString taskName(GetName()) ; 
+  taskName.ReplaceAll(Version(), "") ;
+  
+  if( strcmp(GetName(), "") !=0 ) {  
     // Print parameters
-    
-    cout << "---------------"<< GetName() << " " << GetTitle()<< "-----------" << endl 
-        << "Clusterizing digits from the file: " << fHeaderFileName.Data() << endl 
-        << "                           Branch: " << fDigitsBranchTitle.Data() << endl 
-        << endl 
-        << "                       EMC Clustering threshold = " << fEmcClusteringThreshold << endl
-        << "                       EMC Local Maximum cut    = " << fEmcLocMaxCut << endl
-        << "                       EMC Logarothmic weight   = " << fW0 << endl
-        << endl
-        << "                       CPV Clustering threshold = " << fCpvClusteringThreshold << endl
-        << "                       CPV Local Maximum cut    = " << fCpvLocMaxCut << endl
-       << "                       CPV Logarothmic weight   = " << fW0CPV << endl
-        << endl
-        << "                      PPSD  Clustering threshold = " << fPpsdClusteringThreshold << endl;
+    message  = "\n--------------- %s %s -----------\n" ; 
+    message += "Clusterizing digits from the file: %s\n" ;
+    message += "                           Branch: %s\n" ; 
+    message += "                       EMC Clustering threshold = %f\n" ; 
+    message += "                       EMC Local Maximum cut    = %f\n" ; 
+    message += "                       EMC Logarothmic weight   = %f\n" ;
+    message += "                       CPV Clustering threshold = %f\n" ;
+    message += "                       CPV Local Maximum cut    = %f\n" ;
+    message += "                       CPV Logarothmic weight   = %f\n" ;
     if(fToUnfold)
-      cout << " Unfolding on " << endl ;
+      message += " Unfolding on\n" ;
     else
-      cout << " Unfolding off " << endl ;
+      message += " Unfolding off\n" ;
     
-    cout << "------------------------------------------------------------------" <<endl ;
+    message += "------------------------------------------------------------------" ;
   }
-  else
-    cout << " AliPHOSClusterizerv1 not initialized " << endl ;
+  else 
+    message = " AliPHOSClusterizerv1 not initialized " ;
+  
+  AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),  
+       taskName.Data(), 
+       GetTitle(),
+       taskName.Data(), 
+       GetName(), 
+       fEmcClusteringThreshold, 
+       fEmcLocMaxCut, 
+       fW0, 
+       fCpvClusteringThreshold, 
+       fCpvLocMaxCut, 
+       fW0CPV )) ; 
 }
+
+
 //____________________________________________________________________________
-void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option){
-  //Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
+void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
+{
+  // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
 
-  cout << "AliPHOSClusterizerv1: " << endl ;
-  cout << "       Found "<< fEmcRecPoints->GetEntriesFast() << " EMC Rec Points and " 
-          << fCpvRecPoints->GetEntriesFast() << " CPV RecPoints" << endl ;
+  AliPHOSGetter * gime = AliPHOSGetter::Instance();
+  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
+  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
 
+  AliInfo(Form("\nevent %d \n    Found %d EMC RecPoints and %d CPV RecPoints", 
+              gAlice->GetEvNumber(),
+              emcRecPoints->GetEntriesFast(),
+              cpvRecPoints->GetEntriesFast() ))  ;
+  fRecPointsInRun +=  emcRecPoints->GetEntriesFast() ; 
+  fRecPointsInRun +=  cpvRecPoints->GetEntriesFast() ; 
+  
+  
   if(strstr(option,"all")) {
-    cout << "EMC clusters " << endl ;
-    cout << "  Index    " 
-        << "  Ene(MeV) " 
-        << "  Multi    "
-        << "  Module   "  
-        << "    X      "
-        << "    Y      "
-        << "    Z      "
-        << " Lambda 1  "
-        << " Lambda 2  "
-        << " MaxEnergy "
-        << " # of prim "
-        << " Primaries list "      <<  endl;      
-    
+    printf("\n EMC clusters \n") ;
+    printf("Index    Ene(MeV) Multi Module     X    Y   Z    Lambdas_1  Lambda_2  # of prim  Primaries list\n") ;      
     Int_t index ;
-    for (index = 0 ; index < fEmcRecPoints->GetEntries() ; index++) {
-      AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEmcRecPoints->At(index) ; 
-
-      cout << setw(6) << rp->GetIndexInList() << "     ";
-      cout << setw(6) << rp->GetEnergy()      << "     ";
-      cout << setw(6) << rp->GetMultiplicity()<< "     ";
-      cout << setw(6) << rp->GetPHOSMod()     << "     ";
-
+    for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
+      AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ; 
       TVector3  locpos;  
       rp->GetLocalPosition(locpos);
-      cout << setw(8) <<  locpos.X()          << "     ";
-      cout << setw(8) <<  locpos.Y()          << "     ";
-      cout << setw(8) <<  locpos.Z()          << "     ";
-
       Float_t lambda[2]; 
       rp->GetElipsAxis(lambda);
-      cout << setw(10)<<  lambda[0]           << "     ";
-      cout << setw(10)<<  lambda[1]           << "     ";
-      
-      
       Int_t * primaries; 
       Int_t nprimaries;
       primaries = rp->GetPrimaries(nprimaries);
-      cout << setw(8) <<    primaries         << "     ";
-
-      for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
-       cout << setw(4)  <<  primaries[iprimary] << " ";
-      cout << endl;     
-    }
-    cout << endl ;
-
-    //Now plot PCV/PPSD recPoints
-    cout << "EMC clusters " << endl ;
-    cout << "  Index    " 
-        << "  Multi    "
-        << "  Module   "  
-        << "  Layer    "  
-        << "    X      "
-        << "    Y      "
-        << "    Z      "
-        << " # of prim "
-        << " Primaries list "      <<  endl;      
-    
-    for (index = 0 ; index < fEmcRecPoints->GetEntries() ; index++) {
-      AliPHOSRecPoint * rp = (AliPHOSRecPoint * )fCpvRecPoints->At(index) ; 
-      cout << setw(6) << rp->GetIndexInList() << "     ";
-      cout << setw(6) << rp->GetPHOSMod()     << "     ";
-
-      if( (strcmp(rp->ClassName() , "AliPHOSPpsdRecPoint" )) == 0){
-       AliPHOSPpsdRecPoint * ppsd = (AliPHOSPpsdRecPoint*) rp ;
-       if(ppsd->GetUp())
-         cout <<"        CPV     ";
-       else
-         cout <<"       PPSD     ";
+      printf("\n%6d  %8.2f  %3d     %2d     %4.1f   %4.1f   %4.1f   %4f  %4f    %2d     : ", 
+             rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(), 
+             locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ; 
+      
+      for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
+       printf("%d ", primaries[iprimary] ) ; 
       }
-      else
-       cout <<"        CPV     ";
+      printf("\n") ;
+    }
+
+    //Now plot CPV recPoints
+    printf("\n CPV clusters \n") ;
+    printf("Index    Ene(MeV) Module     X     Y    Z  \n") ;      
+    for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
+      AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )cpvRecPoints->At(index) ; 
       
       TVector3  locpos;  
       rp->GetLocalPosition(locpos);
-      cout << setw(8) <<  locpos.X()          << "     ";
-      cout << setw(8) <<  locpos.Y()          << "     ";
-      cout << setw(8) <<  locpos.Z()          << "     ";
       
-      Int_t * primaries; 
-      Int_t nprimaries;
-      primaries = rp->GetPrimaries(nprimaries);
-      cout << setw(8) <<    primaries         << "     ";
-
-      for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
-       cout << setw(4)  <<  primaries[iprimary] << " ";
-      cout << endl;     
+      printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f \n", 
+            rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(), 
+            locpos.X(), locpos.Y(), locpos.Z()) ; 
     }
-
-
-    cout << "-------------------------------------------------"<<endl ;
   }
 }