Loaders removed from the reconstruction code (C.Cheshkov)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSClusterizerv1.cxx
index 1c5d107..f801bf9 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.112  2007/08/22 09:20:50  hristov
+ * Updated QA classes (Yves)
+ *
+ * Revision 1.111  2007/08/08 12:11:28  kharlov
+ * Protection against uninitialized fQADM
+ *
+ * Revision 1.110  2007/08/07 14:16:00  kharlov
+ * Quality assurance added (Yves Schutz)
+ *
+ * Revision 1.109  2007/07/24 17:20:35  policheh
+ * Usage of RecoParam objects instead of hardcoded parameters in reconstruction.
+ * (See $ALICE_ROOT/PHOS/macros/BeamTest2006/RawReconstruction.C).
+ *
+ * Revision 1.108  2007/06/18 07:00:51  kharlov
+ * Bug fix for attempt to use AliPHOSEmcRecPoint after its deletion
+ *
+ * Revision 1.107  2007/05/25 14:12:26  policheh
+ * Local to tracking CS transformation added for CPV rec. points
+ *
+ * Revision 1.106  2007/05/24 13:01:22  policheh
+ * Local to tracking CS transformation invoked for each EMC rec.point
+ *
+ * Revision 1.105  2007/05/02 13:41:22  kharlov
+ * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
+ *
+ * Revision 1.104  2007/04/27 16:55:53  kharlov
+ * Calibration stops if PHOS CDB objects do not exist
+ *
+ * Revision 1.103  2007/04/11 11:55:45  policheh
+ * SetDistancesToBadChannels() added.
+ *
+ * Revision 1.102  2007/03/28 19:18:15  kharlov
+ * RecPoints recalculation in TSM removed
+ *
+ * Revision 1.101  2007/03/06 06:51:27  kharlov
+ * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
+ *
+ * Revision 1.100  2007/01/10 11:07:26  kharlov
+ * Raw digits writing to file (B.Polichtchouk)
+ *
+ * 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 
 //  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 [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1(<pointer_to_phos_geometry_onject>)  
+//  root [1] cl->Digits2Clusters(digitsTree,clusterTree)
+//               //finds RecPoints in the current event
 //  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")  
-//               //once more finds RecPoints options are 
-//               // deb - print number of found rec points
-//               // deb all - print number of found RecPoints and some their characteristics 
-//               // time - print benchmarking results
 
 // --- 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"
+#include "TClonesArray.h"
 
 // --- Standard library ---
 
 // --- AliRoot header files ---
+#include "AliLog.h"
+#include "AliConfig.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 "AliPHOSGetter.h"
-#include "AliRun.h"
-#include "AliPHOSCalibrManager.h"
-#include "AliPHOSCalibrationData.h"
+#include "AliPHOSCalibrationDB.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
+#include "AliCDBEntry.h"
+#include "AliPHOSRecoParam.h"
+#include "AliPHOSQualAssDataMaker.h" 
+#include "AliPHOSCalibData.h"
+#include "AliPHOSReconstructor.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),                  fEmcTimeGate(0),
+  fIsOldRCUFormat(0)
 {
   // default ctor (to be used mainly by Streamer)
   
@@ -92,15 +192,22 @@ ClassImp(AliPHOSClusterizerv1)
 }
 
 //____________________________________________________________________________
-AliPHOSClusterizerv1::AliPHOSClusterizerv1(const char* headerFile,const char* name, const Bool_t toSplit)
-:AliPHOSClusterizer(headerFile, name, toSplit)
+AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
+  AliPHOSClusterizer(geom),
+  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),                  fEmcTimeGate(0),
+  fIsOldRCUFormat(0)
 {
   // ctor with the indication of the file where header Tree and digits Tree are stored
   
   InitParameters() ;
   Init() ;
   fDefaultInit = kFALSE ; 
-
 }
 
 //____________________________________________________________________________
@@ -108,130 +215,109 @@ AliPHOSClusterizerv1::AliPHOSClusterizerv1(const char* headerFile,const char* na
 {
   // dtor
 
-  delete fPedestals ;
-  delete fGains ;
-
-  fSplitFile = 0 ; 
-
 }
 
 //____________________________________________________________________________
-const TString AliPHOSClusterizerv1::BranchName() const 
+Float_t  AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
 {  
-  // Returns branch name of the AliPHOSClusterizer
-
-  TString branchName(GetName() ) ;
-  branchName.Remove(branchName.Index(Version())-1) ;
-  return branchName ;
-}
-//____________________________________________________________________________
-Float_t  AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId) const
-{ 
-  // Returns energy value in the cell absId with the digit amplitude amp
-
-  if(fPedestals || fGains ){  //use calibration data
-    if(!fPedestals || !fGains ){
-      Error("Calibrate","Either Pedestals of Gains not set!") ;
-      return 0 ;
+  // 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];
+    fGeom->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 ;        
     }
-    Float_t en=(amp - fPedestals->Data(absId))*fGains->Data(absId) ;
-    if(en>0)
-      return en ;
-    else
-      return 0. ;
   }
   else{ //simulation
-    if(absId <= fEmcCrystals) //calibrate as EMC 
+    if(absId <= fEmcCrystals) // this is EMC 
       return fADCpedestalEmc + amp*fADCchanelEmc ;        
-    else //calibrate as CPV
+  }
+  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];
+    fGeom->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;
 }
 
 //____________________________________________________________________________
-void AliPHOSClusterizerv1::Exec(Option_t * option)
+void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
 {
-  // Steering method
-  if( strcmp(GetName(), "")== 0 ) 
-    Init() ;
+  // Steering method to perform clusterization for one event
+  // The input is the tree with digits.
+  // The output is the tree with clusters.
 
   if(strstr(option,"tim"))
     gBenchmark->Start("PHOSClusterizer"); 
   
-  if(strstr(option,"print"))
-    Print("") ; 
-
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
-  if(gime->BranchExists("RecPoints"))
+  if(strstr(option,"print")) {
+    Print() ; 
     return ;
-
-  //test, if this is real data:
-  TBranch * eh = gAlice->TreeE()->GetBranch("AliPHOSBeamTestEvent") ;
-  if(eh){ //Read CalibrationData
-    AliPHOSCalibrManager * cmngr=AliPHOSCalibrManager::GetInstance() ;
-    if(!cmngr){ //not defined yet
-      cmngr=AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
-    } 
-    if(fPedestals)
-      delete fPedestals ;
-    fPedestals = new AliPHOSCalibrationData("Pedestals",fCalibrVersion) ;
-    cmngr->ReadFromRoot(*fPedestals,fCalibrRun) ;
-    if(fGains)
-      delete fGains ;
-    fGains = new AliPHOSCalibrationData("Gains",fCalibrVersion) ;
-    cmngr->ReadFromRoot(*fGains,fCalibrRun) ;
   }
 
-  Int_t nevents = gime->MaxEvent() ;
-  Int_t ievent ;
-  
-  for(ievent = 0; ievent < nevents; ievent++){
-    gime->Event(ievent,"D") ;
-    Int_t pattern = gime->EventPattern() ;
-    if(pattern){
-           printf("pattern %d \n",pattern) ;
-      if(!fPatterns){
-       Error("Exec","Patterns for reconstruction of events are not set, use SetPatterns() ") ;
-       return ;
-      }
-      Bool_t skip = kTRUE ;
-      Int_t npat = fPatterns->GetSize() ;
-      for(Int_t i = 0; i<npat;i++)
-       if(pattern==fPatterns->At(i)){
-         skip = kFALSE ;
-         break ;
-       }
-      if(skip)
-       continue ;
-    }
+  GetCalibrationParameters() ;
 
-    if(ievent == 0)
-      GetCalibrationParameters() ;
-    
-    fNumberOfEmcClusters  = fNumberOfCpvClusters  = 0 ;
-            
-    MakeClusters() ;
+  MakeClusters() ;
     
-    if(fToUnfold)             
-      MakeUnfolding() ;
+  AliDebug(2,Form(" ---- Printing clusters (%d)\n",
+                 fEMCRecPoints->GetEntries()));
+  if(AliLog::GetGlobalDebugLevel()>1)
+    fEMCRecPoints->Print();
+
+  if(fToUnfold)             
+    MakeUnfolding();
+
+    //makes the quality assurance data
+  if (GetQualAssDataMaker()) {
+    GetQualAssDataMaker()->SetData(fEMCRecPoints) ; 
+    GetQualAssDataMaker()->Exec(AliQualAss::kRECPOINTS) ; 
+    GetQualAssDataMaker()->SetData(fCPVRecPoints) ; 
+    GetQualAssDataMaker()->Exec(AliQualAss::kRECPOINTS) ; 
+  }
 
-    WriteRecPoints(ievent) ;
+  WriteRecPoints();
 
-    if(strstr(option,"deb"))  
-      PrintRecPoints(option) ;
+  if(strstr(option,"deb"))  
+    PrintRecPoints(option) ;
+
+  // PLEASE FIX BY MOVING IT TO ALIRECONSTRUCTION !!!
+  //Write the quality assurance data only after the last event 
+  //  if (GetQualAssDataMaker() && fEventCounter == gime->MaxEvent()) 
+  //    GetQualAssDataMaker()->Finish(AliQualAss::kRECPOINTS) ;
 
-    //increment the total number of digits per run 
-    fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;  
-    fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;  
- }
-  
   if(strstr(option,"tim")){
     gBenchmark->Stop("PHOSClusterizer");
-    Info("Exec", "  took %f seconds for Clusterizing %f seconds per event \n",
-        gBenchmark->GetCpuTime("PHOSClusterizer"), 
-        gBenchmark->GetCpuTime("PHOSClusterizer")/nevents ) ; 
-  }
+    AliInfo(Form("took %f seconds for Clusterizing\n",
+                gBenchmark->GetCpuTime("PHOSClusterizer"))); 
+  } 
 }
 
 //____________________________________________________________________________
@@ -242,10 +328,7 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit **
   // 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::GetInstance() ; 
-  TClonesArray * digits = gime->Digits() ; 
   
-
   gMinuit->mncler();                     // Reset Minuit's list of paramters
   gMinuit->SetPrintLevel(-1) ;           // No Printout
   gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;  
@@ -253,7 +336,8 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit **
 
   TList * toMinuit = new TList();
   toMinuit->AddAt(emcRP,0) ;
-  toMinuit->AddAt(digits,1) ;
+  toMinuit->AddAt(fDigitsArr,1) ;
+  toMinuit->AddAt(fGeom,2) ;
   
   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 
@@ -266,16 +350,14 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit **
 
   Int_t iDigit ;
 
-  const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
-
   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->RelPosInModule(relid, x, z) ;
+    fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+    fGeom->RelPosInModule(relid, x, z) ;
 
     Float_t energy = maxAtEnergy[iDigit] ;
 
@@ -330,24 +412,19 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit **
 //____________________________________________________________________________
 void AliPHOSClusterizerv1::GetCalibrationParameters() 
 {
-  // Get calibration parameters from AliPHOSDigitizer
-
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+  // Set calibration parameters:
+  // if calibration database exists, they are read from database,
+  // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
+  //
+  // It is a user responsilibity to open CDB before reconstruction, for example: 
+  // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
+
+  fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
+  if (fCalibData->GetCalibDataEmc() == 0)
+    AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
+  if (fCalibData->GetCalibDataCpv() == 0)
+    AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
 
-  const TTask * task = gime->Digitizer(BranchName()) ;
-  if(strcmp(task->IsA()->GetName(),"AliPHOSDigitizer")==0){
-    const AliPHOSDigitizer * dig = static_cast<const AliPHOSDigitizer *>(task) ;
-
-    fADCchanelEmc   = dig->GetEMCchannel() ;
-    fADCpedestalEmc = dig->GetEMCpedestal();
-    
-    fADCchanelCpv   = dig->GetCPVchannel() ;
-    fADCpedestalCpv = dig->GetCPVpedestal() ; 
-  }
-  else{
- ;
-//    fCalibrationDB = gime->CalibrationDB();
-  }
 }
 
 //____________________________________________________________________________
@@ -355,93 +432,48 @@ 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
-  
-  if ( strcmp(GetTitle(), "") == 0 )
-    SetTitle("galice.root") ;
-
-  TString branchname = GetName() ;
-  branchname.Remove(branchname.Index(Version())-1) ;
-
-  AliPHOSGetter * gime = AliPHOSGetter::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+="PHOS.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 AliPHOSGeometry * geom = gime->PHOSGeometry() ;
-  fEmcCrystals = geom->GetNModules() *  geom->GetNCristalsInModule() ;
+  fEmcCrystals = fGeom->GetNModules() *  fGeom->GetNCristalsInModule() ;
 
   if(!gMinuit) 
-    gMinuit = new TMinuit(100) ;
+    gMinuit = new TMinuit(100);
 
-  gime->PostClusterizer(this) ;
-  gime->PostRecPoints(branchname) ;
-  
 }
 
 //____________________________________________________________________________
 void AliPHOSClusterizerv1::InitParameters()
 {
-  // Set initial parameters for clusterization
 
   fNumberOfCpvClusters     = 0 ; 
   fNumberOfEmcClusters     = 0 ; 
+
+  const AliPHOSRecoParam* parEmc = AliPHOSReconstructor::GetRecoParamEmc();
+  if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!");
+
+  const AliPHOSRecoParam* parCpv = AliPHOSReconstructor::GetRecoParamCpv(); 
+  if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!");
+
+  fCpvClusteringThreshold  = parCpv->GetClusteringThreshold();
+  fEmcClusteringThreshold  = parEmc->GetClusteringThreshold();
   
-  fCpvClusteringThreshold  = 0.0;
-  fEmcClusteringThreshold  = 0.2;   
-  
-  fEmcLocMaxCut            = 0.03 ;
-  fCpvLocMaxCut            = 0.03 ;
-  
-  fW0                      = 4.5 ;
-  fW0CPV                   = 4.0 ;
+  fEmcLocMaxCut            = parEmc->GetLocalMaxCut();
+  fCpvLocMaxCut            = parCpv->GetLocalMaxCut();
+
+  fEmcMinE                 = parEmc->GetMinE();
+  fCpvMinE                 = parCpv->GetMinE();
+
+  fW0                      = parEmc->GetLogWeight();
+  fW0CPV                   = parCpv->GetLogWeight();
 
   fEmcTimeGate             = 1.e-8 ; 
   
   fToUnfold                = kTRUE ;
+    
+  fWrite                   = kTRUE ;
 
-  fPurifyThreshold         = 0.012 ;
-  
-  fPedestals = 0 ;
-  fGains     = 0 ;
-  fCalibrVersion= "v1" ;
-  fCalibrRun = 1 ;
-  fPatterns  = 0 ;
-
-  TString clusterizerName( GetName()) ;
-  if (clusterizerName.IsNull() ) 
-    clusterizerName = "Default" ; 
-  clusterizerName.Append(":") ; 
-  clusterizerName.Append(Version()) ; 
-  SetName(clusterizerName) ;
-  fRecPointsInRun          = 0 ;
-
+  fCalibData               = 0 ;
 
+  fIsOldRCUFormat          = kFALSE;
 }
 
 //____________________________________________________________________________
@@ -454,15 +486,13 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)c
   // 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  
 
-  const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
-
   Int_t rv = 0 ; 
 
   Int_t relid1[4] ; 
-  geom->AbsToRelNumbering(d1->GetId(), relid1) ; 
+  fGeom->AbsToRelNumbering(d1->GetId(), relid1) ; 
 
   Int_t relid2[4] ; 
-  geom->AbsToRelNumbering(d2->GetId(), relid2) ; 
+  fGeom->AbsToRelNumbering(d2->GetId(), relid2) ; 
  
   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module 
     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
@@ -474,7 +504,7 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)c
     }
     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 
     }
 
   } 
@@ -487,7 +517,23 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)c
 
   return rv ; 
 }
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
+{
+  // Remove digits with amplitudes below threshold
 
+  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) ;     
+  }
+}
 
 //____________________________________________________________________________
 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
@@ -495,9 +541,8 @@ Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
  
   Bool_t rv = kFALSE ; 
-  const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
 
-  Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();  
+  Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();  
 
   if(digit->GetId() <= nEMC )   rv = kTRUE; 
 
@@ -510,9 +555,8 @@ 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 ; 
-  const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
-
-  Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
+  
+  Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 
   if(digit->GetId() > nEMC )   rv = kTRUE;
 
@@ -520,92 +564,57 @@ Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 }
 
 //____________________________________________________________________________
-void AliPHOSClusterizerv1::WriteRecPoints(Int_t event)
+void AliPHOSClusterizerv1::WriteRecPoints()
 {
 
   // Creates new branches with given title
   // fills and writes into TreeR.
-  
-  
-  AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; 
-  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
-  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
-  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();
-  }
-  
-  if(!treeR){
-    if(fSplitFile)
-      gAlice->MakeTree("R", fSplitFile);
-    else
-      gAlice->MakeTree("R", gROOT->GetFile(GetTitle()));
-    treeR = gAlice->TreeR() ;
-  }
 
   Int_t index ;
-  //Evaluate position, dispersion and other RecPoint properties...
-  for(index = 0; index < emcRecPoints->GetEntries(); index++){
-    AliPHOSEmcRecPoint * emcrp = static_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(index)) ;
-    emcrp->Purify(fPurifyThreshold) ; 
-    if(emcrp->GetMultiplicity())
-      emcrp->EvalAll(fW0,digits) ;
-    else
-      emcRecPoints->Remove(emcrp) ;
- }
+  //Evaluate position, dispersion and other RecPoint properties..
+  Int_t nEmc = fEMCRecPoints->GetEntriesFast();
+  for(index = 0; index < nEmc; index++){
+    AliPHOSEmcRecPoint * rp =
+      dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
+    rp->Purify(fEmcMinE) ;
+    if(rp->GetMultiplicity()==0){
+      fEMCRecPoints->RemoveAt(index) ;
+      delete rp ;
+      continue;
+    }
 
-  emcRecPoints->Compress() ;  
-  emcRecPoints->Sort() ;
-  for(index = 0; index < emcRecPoints->GetEntries(); index++)
-    ((AliPHOSEmcRecPoint *)emcRecPoints->At(index))->SetIndexInList(index) ;
-  
-  emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ; 
+    // No vertex is available now, calculate corrections in PID
+    rp->EvalAll(fW0,fDigitsArr) ;
+    TVector3 fakeVtx(0.,0.,0.) ;
+    rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
+    rp->EvalLocal2TrackingCSTransform();
+  }
+  fEMCRecPoints->Compress() ;
+  //  fEMCRecPoints->Sort() ; //Can not sort until position is calculated!
+  //  fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
+  for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
+    dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
+  }
   
+  //For each rec.point set the distance to the nearest bad crystal (BVP)
+  SetDistancesToBadChannels();
+
   //Now the same for CPV
-  for(index = 0; index < cpvRecPoints->GetEntries(); index++)
-    ((AliPHOSRecPoint *)cpvRecPoints->At(index))->EvalAll(fW0CPV,digits)  ;
-  
-  cpvRecPoints->Sort() ;
+  for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
+    AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
+    rp->EvalAll(fW0CPV,fDigitsArr) ;
+    rp->EvalLocal2TrackingCSTransform();
+  }
+//  fCPVRecPoints->Sort() ;
   
-  for(index = 0; index < cpvRecPoints->GetEntries(); index++)
-    ((AliPHOSRecPoint *)cpvRecPoints->At(index))->SetIndexInList(index) ;
+  for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
+    dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
   
-  cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
+  fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
   
-  Int_t bufferSize = 32000 ;    
-  Int_t splitlevel = 0 ;
-  //First EMC
-  TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
-  emcBranch->SetTitle(BranchName());
-    
-  //Now CPV branch
-  TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
-  cpvBranch->SetTitle(BranchName());
-    
-  //And Finally  clusterizer branch
-  AliPHOSClusterizerv1 * cl = this ;
-  TBranch * clusterizerBranch = treeR->Branch("AliPHOSClusterizer","AliPHOSClusterizerv1",
-                                             &cl,bufferSize,splitlevel);
-  clusterizerBranch->SetTitle(BranchName());
-
-  emcBranch        ->Fill() ;
-  cpvBranch        ->Fill() ;
-  clusterizerBranch->Fill() ;
-
-  treeR->AutoSave() ; 
-  if(gAlice->TreeR()!=treeR)
-    treeR->Delete();
+  if(fWrite){ //We write TreeR
+    fTreeR->Fill();
+  }
 }
 
 //____________________________________________________________________________
@@ -614,19 +623,10 @@ 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
 
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
+  //Remove digits below threshold
+  CleanDigits(fDigitsArr) ;
 
-  TObjArray * emcRecPoints = gime->EmcRecPoints(BranchName()) ; 
-  TObjArray * cpvRecPoints = gime->CpvRecPoints(BranchName()) ; 
-  emcRecPoints->Delete() ;
-  cpvRecPoints->Delete() ;
-  
-  TClonesArray * digits = gime->Digits() ; 
-   if ( !digits ) {
-   Fatal("MakeClusters", "Digits with name %s not found !", BranchName().Data() ) ;  
-  } 
-  TClonesArray * digitsC =  (TClonesArray*)digits->Clone() ;
-  
+  TClonesArray * digitsC =  static_cast<TClonesArray*>( fDigitsArr->Clone() ) ;
  
   // Clusterization starts  
 
@@ -634,58 +634,61 @@ void AliPHOSClusterizerv1::MakeClusters()
   AliPHOSDigit * digit ; 
   Bool_t notremoved = kTRUE ;
 
-  while ( (digit = (AliPHOSDigit *)nextdigit()) ) { // scan over the list of digitsC
+  while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
+    
+
     AliPHOSRecPoint * clu = 0 ; 
 
     TArrayI clusterdigitslist(1500) ;   
     Int_t index ;
-   // cout << "digit " << digit << " " << Calibrate(digit->GetAmp(),digit->GetId()) << endl ;
 
-    if (( IsInEmc (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fEmcClusteringThreshold  ) || 
-        ( IsInCpv (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > 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 >= emcRecPoints->GetSize()) 
-         emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
-       
-       emcRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
-       clu = (AliPHOSEmcRecPoint *) emcRecPoints->At(fNumberOfEmcClusters) ; 
-       fNumberOfEmcClusters++ ; 
-       clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId())) ; 
-       clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;  
-       iDigitInCluster++ ; 
-       digitsC->Remove(digit) ; 
+        // start a new EMC RecPoint
+        if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize()) 
+          fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
+          
+        fEMCRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
+        clu = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ; 
+       fNumberOfEmcClusters++ ; 
+       clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
+        clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
+        iDigitInCluster++ ;
+        digitsC->Remove(digit) ; 
 
       } else { 
-       
-       // start a new CPV cluster
-       if(fNumberOfCpvClusters >= cpvRecPoints->GetSize()) 
-         cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
-
-       cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
-
-       clu =  (AliPHOSCpvRecPoint *) cpvRecPoints->At(fNumberOfCpvClusters)  ;  
-       fNumberOfCpvClusters++ ; 
-       clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId()) ) ;     
-       clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ; 
-       iDigitInCluster++ ; 
-       digitsC->Remove(digit) ; 
+        
+        // start a new CPV cluster
+        if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize()) 
+          fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
+
+        fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
+
+        clu =  dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;  
+        fNumberOfCpvClusters++ ; 
+        clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;        
+        clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;        
+        iDigitInCluster++ ; 
+        digitsC->Remove(digit) ; 
         nextdigit.Reset() ;
-       
-       // Here we remove remaining EMC digits, which cannot make a 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) ) 
-             digitsC->Remove(digit) ;
+              digitsC->Remove(digit) ;
             else 
-             break ;
-         }
-         notremoved = kFALSE ;
-       }
-       
+              break ;
+          }
+          notremoved = kFALSE ;
+        }
+        
       } // else        
       
       nextdigit.Reset() ;
@@ -693,28 +696,32 @@ void AliPHOSClusterizerv1::MakeClusters()
       AliPHOSDigit * digitN ; 
       index = 0 ;
       while (index < iDigitInCluster){ // scan over digits already in cluster 
-       digit =  (AliPHOSDigit*)digits->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*>( fDigitsArr->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(), digitN->GetId() ) ) ;
-           clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; 
-           iDigitInCluster++ ; 
-           digitsC->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  
@@ -732,44 +739,35 @@ void AliPHOSClusterizerv1::MakeUnfolding()
   // Unfolds clusters using the shape of an ElectroMagnetic shower
   // Performs unfolding of all EMC/CPV clusters
 
-  AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
-  
-
-  const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
-  TObjArray * emcRecPoints = gime->EmcRecPoints(BranchName()) ; 
-  TObjArray * cpvRecPoints = gime->CpvRecPoints(BranchName()) ; 
-  TClonesArray * digits = gime->Digits() ; 
-  
   // Unfold first EMC clusters 
   if(fNumberOfEmcClusters > 0){
 
-    Int_t nModulesToUnfold = geom->GetNModules() ; 
+    Int_t nModulesToUnfold = fGeom->GetNModules() ; 
 
     Int_t numberofNotUnfolded = fNumberOfEmcClusters ; 
     Int_t index ;   
     for(index = 0 ; index < numberofNotUnfolded ; index++){
       
-      AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint *) emcRecPoints->At(index) ;
+      AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
-       break ;
+        break ;
       
       Int_t nMultipl = emcRecPoint->GetMultiplicity() ; 
       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
-      Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
+      Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
       
       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0       
-       UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
-       if(emcRecPoint->GetNExMax()!=-1){ //Do not remove clusters which failed to unfold
-         emcRecPoints->Remove(emcRecPoint); 
-         emcRecPoints->Compress() ;
-         index-- ;
-         fNumberOfEmcClusters -- ;
-         numberofNotUnfolded-- ;
-       }
+        UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
+        fEMCRecPoints->Remove(emcRecPoint); 
+        fEMCRecPoints->Compress() ;
+        index-- ;
+        fNumberOfEmcClusters -- ;
+        numberofNotUnfolded-- ;
+      }
+      else{
+        emcRecPoint->SetNExMax(1) ; //Only one local maximum
       }
-      else
-       emcRecPoint->SetNExMax(1) ; //Only one local maximum
       
       delete[] maxAt ; 
       delete[] maxAtEnergy ; 
@@ -781,31 +779,31 @@ void AliPHOSClusterizerv1::MakeUnfolding()
   // Unfold now CPV clusters
   if(fNumberOfCpvClusters > 0){
     
-    Int_t nModulesToUnfold = geom->GetNModules() ;
+    Int_t nModulesToUnfold = fGeom->GetNModules() ;
 
     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;     
     Int_t index ;   
     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
       
-      AliPHOSRecPoint * recPoint = (AliPHOSRecPoint *) cpvRecPoints->At(index) ;
+      AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
 
       if(recPoint->GetPHOSMod()> nModulesToUnfold)
-       break ;
+        break ;
       
-      AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint*) recPoint ; 
+      AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ; 
       
       Int_t nMultipl = emcRecPoint->GetMultiplicity() ; 
       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
-      Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
+      Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
       
       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0       
-       UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
-       cpvRecPoints->Remove(emcRecPoint); 
-       cpvRecPoints->Compress() ;
-       index-- ;
-       numberofCpvNotUnfolded-- ;
-       fNumberOfCpvClusters-- ;
+        UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
+        fCPVRecPoints->Remove(emcRecPoint); 
+        fCPVRecPoints->Compress() ;
+        index-- ;
+        numberofCpvNotUnfolded-- ;
+        fNumberOfCpvClusters-- ;
       }
       
       delete[] maxAt ; 
@@ -817,38 +815,34 @@ void AliPHOSClusterizerv1::MakeUnfolding()
 }
 
 //____________________________________________________________________________
-Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t r)
+Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
 { 
   // Shape of the shower (see PHOS TDR)
   // 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) ;
+  //for the moment we neglect dependence on the incident angle.  
+
+  Double_t r2    = x*x + z*z ;
+  Double_t r4    = r2*r2 ;
+  Double_t r295  = TMath::Power(r2, 2.95/2.) ;
   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
   return shape ;
 }
 
 //____________________________________________________________________________
 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, 
-                                                Int_t nMax, 
-                                                AliPHOSDigit ** 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::GetInstance() ; 
-
-  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 
+    // Fit failed, return and remove cluster
     iniEmc->SetNExMax(-1) ;
     delete[] fitparameters ; 
     return ;
@@ -860,18 +854,20 @@ void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 
   Int_t nDigits = iniEmc->GetMultiplicity() ;  
   Float_t * efit = new Float_t[nDigits] ;
-  Float_t xDigit=0.,zDigit=0.,distance=0. ;
+  Float_t xDigit=0.,zDigit=0. ;
   Float_t xpar=0.,zpar=0.,epar=0.  ;
   Int_t relid[4] ;
   AliPHOSDigit * digit = 0 ;
   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 
+  TVector3 vIncid ;
+
   Int_t iparam ;
   Int_t iDigit ;
   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
-    digit = (AliPHOSDigit*) digits->At(emcDigits[iDigit] ) ;   
-    geom->AbsToRelNumbering(digit->GetId(), relid) ;
-    geom->RelPosInModule(relid, xDigit, zDigit) ;
+    digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;   
+    fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+    fGeom->RelPosInModule(relid, xDigit, zDigit) ;
     efit[iDigit] = 0;
 
     iparam = 0 ;    
@@ -880,9 +876,9 @@ void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
       zpar = fitparameters[iparam+1] ;
       epar = fitparameters[iparam+2] ;
       iparam += 3 ;
-      distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
-      distance =  TMath::Sqrt(distance) ;
-      efit[iDigit] += epar * ShowerShape(distance) ;
+//      fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
+//      efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
+      efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
     }
   }
   
@@ -900,39 +896,39 @@ void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
     zpar = fitparameters[iparam+1] ;
     epar = fitparameters[iparam+2] ;
     iparam += 3 ;    
+//    fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
     
     AliPHOSEmcRecPoint * emcRP = 0 ;  
 
-    if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
+    if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
       
-      if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
-       emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
+      if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
+        fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
       
-      (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
-      emcRP = (AliPHOSEmcRecPoint *) emcRecPoints->At(fNumberOfEmcClusters);
+      (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
+      emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
       fNumberOfEmcClusters++ ;
       emcRP->SetNExMax((Int_t)nPar/3) ;
     }
-    else{//create new entries in fCpvRecPoints
-      if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
-       cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
+    else{//create new entries in fCPVRecPoints
+      if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
+        fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
       
-      (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
-      emcRP = (AliPHOSEmcRecPoint *) cpvRecPoints->At(fNumberOfCpvClusters);
+      (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
+      emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
       fNumberOfCpvClusters++ ;
     }
     
     Float_t eDigit ;
     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
-      digit = (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] ; 
+      digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ; 
+      fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+      fGeom->RelPosInModule(relid, xDigit, zDigit) ;
+//      ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ; 
+      ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ; 
       eDigit = emcEnergies[iDigit] * ratio ;
       emcRP->AddDigit( *digit, eDigit ) ;
-    }  
+    }        
   }
  
   delete[] fitparameters ; 
@@ -946,22 +942,25 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou
   // Calculates the Chi square for the cluster unfolding minimization
   // Number of parameters, Gradient, Chi squared, parameters, what to do
 
-  TList * toMinuit = (TList*) gMinuit->GetObjectFit() ;
-
-  AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint*) toMinuit->At(0)  ;
-  TClonesArray * digits = (TClonesArray*)toMinuit->At(1)  ;
+  TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
 
+  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
+  TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
+  // A bit buggy way to get an access to the geometry
+  // To be revised!
+  AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
 
+//  TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ;  //Vertex position
   
-  //  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 nOdigits = emcRP->GetDigitsMultiplicity() ; 
 
   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
-
-  const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ; 
+  
+//  TVector3 vInc ;
   fret = 0. ;     
   Int_t iparam ;
 
@@ -976,7 +975,7 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou
 
   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 ;
@@ -990,36 +989,37 @@ 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 dx = (xDigit - x[iParam]) ;
+         iParam++ ; 
+         Double_t dz = (zDigit - x[iParam]) ; 
+         iParam++ ;          
+//         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
+//         efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
+         efit += x[iParam] * ShowerShape(dx,dz) ;
+         iParam++ ;
        }
-       Double_t sum=0 ;        
-       if(emcEnergies[iDigit] > 0)
-        sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; 
-       // Here we assume, that sigma = sqrt(E) 
+       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] ;
+//         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
+         Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
+//         Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
+         Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
+//DP: No incident angle dependence in gradient yet!!!!!!
+         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;
@@ -1030,21 +1030,19 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou
        Double_t zpar = x[iparam+1] ;
        Double_t epar = x[iparam+2] ;
        iparam += 3 ;
-       Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
-       distance =  TMath::Sqrt(distance) ;
-       efit += epar * ShowerShape(distance) ;
+//       fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
+//       efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
+       efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
      }
 
-     if(emcEnergies[iDigit] > 0) //for the case if rounding error
-       fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/
-               emcEnergies[iDigit] ; 
+     fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ; 
      // Here we assume, that sigma = sqrt(E)
   }
 
 }
 
 //____________________________________________________________________________
-void AliPHOSClusterizerv1::Print(Option_t * option)const
+void AliPHOSClusterizerv1::Print(const Option_t *)const
 {
   // Print clusterizer parameters
 
@@ -1073,7 +1071,7 @@ void AliPHOSClusterizerv1::Print(Option_t * option)const
   else 
     message = " AliPHOSClusterizerv1 not initialized " ;
   
-  Info("Print", message.Data(),  
+  AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),  
        taskName.Data(), 
        GetTitle(),
        taskName.Data(), 
@@ -1083,32 +1081,44 @@ void AliPHOSClusterizerv1::Print(Option_t * option)const
        fW0, 
        fCpvClusteringThreshold, 
        fCpvLocMaxCut, 
-       fW0CPV ) ; 
+       fW0CPV )) ; 
 }
-
-
+//____________________________________________________________________________
+//void AliPHOSClusterizerv1::GetVertex(void)
+//{ //Extracts vertex posisition
+//  
+  //ESD
+//DP - todo  if(){
+//
+//  }
+
+//  //MC Generator
+//  if(gAlice && gAlice->GetMCApp() && gAlice->Generator()){
+//    Float_t x,y,z ;
+//    gAlice->Generator()->GetOrigin(x,y,z) ;
+//    fVtx.SetXYZ(x,y,z) ;
+//    return ; 
+//  }
+//
+//  //No any source
+//  fVtx[0]=fVtx[1]=fVtx[2]=0. ;
+//
+//}
 //____________________________________________________________________________
 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
 {
   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
 
-  TObjArray * emcRecPoints = AliPHOSGetter::GetInstance()->EmcRecPoints(BranchName()) ; 
-  TObjArray * cpvRecPoints = AliPHOSGetter::GetInstance()->CpvRecPoints(BranchName()) ; 
-
-  if(strstr(option,"deb")){
-    printf("event %d, found %d EmcRecPoints and %d CPV RecPoints \n",
-           gAlice->GetEvNumber(), 
-          emcRecPoints->GetEntriesFast(),
-          cpvRecPoints->GetEntriesFast()) ; 
-  }
-    
-  
+  AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints", 
+              fEMCRecPoints->GetEntriesFast(),
+              fCPVRecPoints->GetEntriesFast() ))  ;
   if(strstr(option,"all")) {
-    printf("\nEMC clusters\n" );
-    printf("Index    Ene(MeV) Multi Module     X     Y   Z  Lambda 1 Lambda 2  # of prim  Primaries list\n" );      
+    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 < emcRecPoints->GetEntries() ; index++) {
-      AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ; 
+    for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
+      AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ; 
       TVector3  locpos;  
       rp->GetLocalPosition(locpos);
       Float_t lambda[2]; 
@@ -1116,34 +1126,70 @@ void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
       Int_t * primaries; 
       Int_t nprimaries;
       primaries = rp->GetPrimaries(nprimaries);
-      printf("\n%6d  %8.2f  %3d     %2d     %4.1f    %4.1f %4.1f %4f  %4f    %2d     : ", 
+      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] ) ; 
       }
-    }      
+      printf("\n") ;
+    }
+
     //Now plot CPV recPoints
-    printf("\n CPV RecPoints \n") ;
-    printf("Index    Ene(MeV) Module     X     Y   Z  # of prim  Primaries list\n") ;      
-    for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
-        AliPHOSRecPoint * rp = (AliPHOSRecPoint * )cpvRecPoints->At(index) ; 
-       
-       TVector3  locpos;  
-       rp->GetLocalPosition(locpos);
-       
-       Int_t * primaries; 
-       Int_t nprimaries ; 
-       primaries = rp->GetPrimaries(nprimaries);
-       printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f %2d     : ", 
-               rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(), 
-               locpos.X(), locpos.Y(), locpos.Z(), nprimaries) ; 
-       primaries = rp->GetPrimaries(nprimaries);
-       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
-         printf("%d ", primaries[iprimary] ) ; 
-       }
-      }
+    printf("\n CPV clusters \n") ;
+    printf("Index    Ene(MeV) Module     X     Y    Z  \n") ;      
+    for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
+      AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ; 
+      
+      TVector3  locpos;  
+      rp->GetLocalPosition(locpos);
+      
+      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()) ; 
+    }
   }
 }
 
+
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::SetDistancesToBadChannels()
+{
+  //For each EMC rec. point set the distance to the nearest bad crystal.
+  //Author: Boris Polichtchouk 
+
+  if(!fCalibData->GetNumOfEmcBadChannels()) return;
+  AliInfo(Form("%d bad channel(s) found.\n",fCalibData->GetNumOfEmcBadChannels()));
+
+  Int_t badIds[8000];
+  fCalibData->EmcBadChannelIds(badIds);
+
+  AliPHOSEmcRecPoint* rp;
+
+  TMatrixF gmat;
+  TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
+  TVector3 gposBadChannel; // global position of bad crystal
+  TVector3 dR;
+
+  Float_t dist,minDist;
+
+  for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
+    rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
+    minDist = 1.e+07;
+
+    for(Int_t iBad=0; iBad<fCalibData->GetNumOfEmcBadChannels(); iBad++) {
+      rp->GetGlobalPosition(gposRecPoint,gmat);
+      fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
+      AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
+                     gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
+                     gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
+      dR = gposBadChannel-gposRecPoint;
+      dist = dR.Mag();
+      if(dist<minDist) minDist = dist;
+    }
+
+    rp->SetDistanceToBadCrystal(minDist); 
+  }
+
+}