Charged jets (pPb): Improved trackcut analysis
[u/mrichter/AliRoot.git] / PHOS / AliPHOSClusterizerv1.cxx
index 5532b13..840ffd0 100644 (file)
 
 /* History of cvs commits:
  *
- * $Log$
+ * $Log: AliPHOSClusterizerv1.cxx,v $
+ * Revision 1.118  2007/12/11 21:23:26  kharlov
+ * Added possibility to swith off unfolding
+ *
+ * Revision 1.117  2007/10/18 08:42:05  kharlov
+ * Bad channels cleaned before clusterization
+ *
+ * Revision 1.116  2007/10/01 20:24:08  kharlov
+ * Memory leaks fixed
+ *
+ * Revision 1.115  2007/09/26 14:22:17  cvetan
+ * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
+ *
+ * Revision 1.114  2007/09/06 16:06:44  kharlov
+ * Absence of sorting results in loose of all unfolded clusters
+ *
+ * Revision 1.113  2007/08/28 12:55:07  policheh
+ * Loaders removed from the reconstruction code (C.Cheshkov)
+ *
+ * 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
  *
 //  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 "TMinuit.h"
 #include "TTree.h" 
 #include "TBenchmark.h"
+#include "TClonesArray.h"
 
 // --- Standard library ---
 
 // --- AliRoot header files ---
 #include "AliLog.h"
-#include "AliPHOSGetter.h"
+#include "AliConfig.h"
 #include "AliPHOSGeometry.h" 
 #include "AliPHOSClusterizerv1.h"
 #include "AliPHOSEmcRecPoint.h"
 #include "AliPHOSCpvRecPoint.h"
 #include "AliPHOSDigit.h"
 #include "AliPHOSDigitizer.h"
-#include "AliPHOSCalibrationDB.h"
 #include "AliCDBManager.h"
 #include "AliCDBStorage.h"
 #include "AliCDBEntry.h"
+#include "AliPHOSRecoParam.h"
+#include "AliPHOSReconstructor.h"
+#include "AliPHOSCalibData.h"
 
 ClassImp(AliPHOSClusterizerv1)
   
 //____________________________________________________________________________
-  AliPHOSClusterizerv1::AliPHOSClusterizerv1() : AliPHOSClusterizer()
+AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
+  AliPHOSClusterizer(),
+  fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
+  fWrite(0),                  
+  fNumberOfEmcClusters(0),    fNumberOfCpvClusters(0),
+  fEmcClusteringThreshold(0), fCpvClusteringThreshold(0), 
+  fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
+  fW0CPV(0),                 
+  fTimeGateLowAmp(0.),        fTimeGateLow(0.),         fTimeGateHigh(0.),  
+  fEcoreRadius(0)
 {
   // default ctor (to be used mainly by Streamer)
   
-  InitParameters() ; 
-  fDefaultInit = kTRUE ; 
+  fDefaultInit = kTRUE ;
+  
+  for(Int_t i=0; i<53760; i++){
+    fDigitsUsed[i]=0 ;
+  }
 }
 
 //____________________________________________________________________________
-AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName)
-:AliPHOSClusterizer(alirunFileName, eventFolderName)
+AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
+  AliPHOSClusterizer(geom),
+  fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
+  fWrite(0),                
+  fNumberOfEmcClusters(0),    fNumberOfCpvClusters(0),
+  fEmcClusteringThreshold(0), fCpvClusteringThreshold(0), 
+  fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
+  fW0CPV(0),                  
+  fTimeGateLowAmp(0.),        fTimeGateLow(0.),         fTimeGateHigh(0.),  
+  fEcoreRadius(0) 
 {
   // ctor with the indication of the file where header Tree and digits Tree are stored
   
-  InitParameters() ;
+  for(Int_t i=0; i<53760; i++){
+    fDigitsUsed[i]=0 ;
+  }
+  
   Init() ;
   fDefaultInit = kFALSE ; 
 }
@@ -148,71 +239,11 @@ AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const T
 
 }
 //____________________________________________________________________________
-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;
-}
-
-//____________________________________________________________________________
-void AliPHOSClusterizerv1::Exec(Option_t *option)
+void AliPHOSClusterizerv1::Digits2Clusters(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.
+  // 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"); 
@@ -222,53 +253,28 @@ void AliPHOSClusterizerv1::Exec(Option_t *option)
     return ;
   }
 
-  GetCalibrationParameters() ;
-
-  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
-  if (fRawReader == 0)
-    gime->SetRawDigits(kFALSE);
-  else
-    gime->SetRawDigits(kTRUE);
-  
-  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 {
-      gime->Event(fRawReader,"W",fIsOldRCUFormat); // Read digits from raw data
-    }
-    fNumberOfEmcClusters  = fNumberOfCpvClusters  = 0 ;
+  MakeClusters() ;
     
-    MakeClusters() ;
+  AliDebug(2,Form(" ---- Printing clusters (%d)\n",
+                 fEMCRecPoints->GetEntries()));
+  if(AliLog::GetGlobalDebugLevel()>1)
+    fEMCRecPoints->Print();
 
-    if(fToUnfold)             
-      MakeUnfolding() ;
+  if(fToUnfold)             
+    MakeUnfolding();
 
-    WriteRecPoints();
+  WriteRecPoints();
 
-    if(strstr(option,"deb"))  
-      PrintRecPoints(option) ;
+  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");
-    AliInfo(Form("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"))); 
+  }
+  fEMCRecPoints->Delete();
+  fCPVRecPoints->Delete();
 }
 
 //____________________________________________________________________________
@@ -280,10 +286,8 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit **
   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 
   
-  AliPHOSGetter * gime = AliPHOSGetter::Instance();
-  TClonesArray * digits = gime->Digits(); 
-  
-
+  if(!gMinuit) //it was deleted by someone else
+    gMinuit = new TMinuit(100) ;
   gMinuit->mncler();                     // Reset Minuit's list of paramters
   gMinuit->SetPrintLevel(-1) ;           // No Printout
   gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;  
@@ -291,7 +295,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
 
@@ -304,16 +309,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] ;
 
@@ -365,33 +368,6 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit **
 
 }
 
-//____________________________________________________________________________
-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();
-    
-      fADCchanelCpv   = dig->GetCPVchannel() ;
-      fADCpedestalCpv = dig->GetCPVpedestal() ; 
-    }
-}
 
 //____________________________________________________________________________
 void AliPHOSClusterizerv1::Init()
@@ -399,20 +375,18 @@ 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());
-
-  AliPHOSGeometry * geom = gime->PHOSGeometry();
-
-  fEmcCrystals = geom->GetNModules() *  geom->GetNCristalsInModule() ;
+  fEmcCrystals = fGeom->GetNModules() *  fGeom->GetNCristalsInModule() ;
 
   if(!gMinuit) 
     gMinuit = new TMinuit(100);
 
-  if ( !gime->Clusterizer() ) {
-    gime->PostClusterizer(this);
-  }
+  if (!fgCalibData) 
+    fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
+  if (fgCalibData->GetCalibDataEmc() == 0)
+    AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
+  if (fgCalibData->GetCalibDataCpv() == 0)   
+    AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");   
+
 }
 
 //____________________________________________________________________________
@@ -421,32 +395,30 @@ 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  ;
+  const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
+  if(!recoParam) AliFatal("Reconstruction parameters are not set!");
 
-  fW0                      = 4.5 ;
-  fW0CPV                   = 4.0 ;
+  recoParam->Print();
 
-  fEmcTimeGate             = 1.e-8 ; 
+  fEmcClusteringThreshold  = recoParam->GetEMCClusteringThreshold();
+  fCpvClusteringThreshold  = recoParam->GetCPVClusteringThreshold();
   
-  fToUnfold                = kTRUE ;
-    
-  fRecPointsInRun          = 0 ;
+  fEmcLocMaxCut            = recoParam->GetEMCLocalMaxCut();
+  fCpvLocMaxCut            = recoParam->GetCPVLocalMaxCut();
 
-  fWrite                   = kTRUE ;
-
-  fCalibData               = 0 ;
+  fW0                      = recoParam->GetEMCLogWeight();
+  fW0CPV                   = recoParam->GetCPVLogWeight();
 
-  SetEventRange(0,-1) ;
+  fTimeGateLowAmp          = recoParam->GetTimeGateAmpThresh() ;
+  fTimeGateLow             = recoParam->GetTimeGateLow() ;
+  fTimeGateHigh            = recoParam->GetTimeGateHigh() ;
 
-  fIsOldRCUFormat          = kFALSE;
+  fEcoreRadius             = recoParam->GetEMCEcoreRadius();
+  
+  fToUnfold                = recoParam->EMCToUnfold() ;
+    
+  fWrite                   = kTRUE ;
 }
 
 //____________________________________________________________________________
@@ -455,58 +427,54 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)c
   // 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
+  //                                       =-1 are not neighbour, continue searching, but do not look before d2 next time
   // neighbours are defined as digits having at least a common vertex 
   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
   //                                      which is compared to a digit (d2)  not yet in a cluster  
 
-  AliPHOSGeometry * geom = AliPHOSGetter::Instance()->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] ) ;  
     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 ; 
+    if (( coldiff <= 1 )  && ( rowdiff <= 1 )){   //At least common vertex
+      //    if (( relid1[2]==relid2[2] && coldiff <= 1 )  || ( relid1[3]==relid2[3] &&  rowdiff <= 1 )){ //common side
+      if((relid1[1] != 0) || CheckTimeGate(d1->GetTime(),d1->GetEnergy(),d2->GetTime(),d2->GetEnergy())) 
+      return 1 ; 
     }
     else {
       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) 
-        rv = 2; //  Difference in row numbers is too large to look further 
+        return 2; //  Difference in row numbers is too large to look further 
     }
+    return 0 ;
 
   } 
   else {
+    if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
+      return -1 ;
+    if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
+      return -1 ;
     
-    if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )  
-      rv=2 ;
+    return 2 ;
 
   }
 
-  return rv ; 
+  return 0 ; 
 }
 //____________________________________________________________________________
-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) ;
+Bool_t AliPHOSClusterizerv1::CheckTimeGate(Float_t t1, Float_t amp1, Float_t t2, Float_t amp2)const{
+  //Check if two cells have reasonable time difference
+  //Note that at low amplitude time is defined up to 1 tick == 100 ns.
+  if(amp1<fTimeGateLowAmp || amp2<fTimeGateLowAmp){
+   return (TMath::Abs(t1 - t2 ) < fTimeGateLow) ;
   }
-  digits->Compress() ;
-  for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) { 
-    AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ; 
-    digit->SetIndexInList(i) ;     
+  else{ //Time should be measured with good accuracy
+   return (TMath::Abs(t1 - t2 ) < fTimeGateHigh) ;
   }
 
 }
@@ -516,9 +484,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 ; 
-  AliPHOSGeometry * geom = AliPHOSGetter::Instance()->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; 
 
@@ -532,9 +499,7 @@ Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
  
   Bool_t rv = kFALSE ; 
   
-  AliPHOSGeometry * geom = AliPHOSGetter::Instance()->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;
 
@@ -542,78 +507,59 @@ Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 }
 
 //____________________________________________________________________________
-void AliPHOSClusterizerv1::Unload() 
-{
-  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
-  gime->PhosLoader()->UnloadDigits() ; 
-  gime->PhosLoader()->UnloadRecPoints() ; 
-}
-
-//____________________________________________________________________________
 void AliPHOSClusterizerv1::WriteRecPoints()
 {
 
   // Creates new branches with given title
   // fills and writes into TreeR.
 
-  AliPHOSGetter * gime = AliPHOSGetter::Instance();
-
-  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();
+  Int_t nEmc = fEMCRecPoints->GetEntriesFast();
+  Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
+  TVector3 fakeVtx(0.,0.,0.) ;
   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) ;
+    AliPHOSEmcRecPoint * rp =
+      static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
+    rp->Purify(emcMinE,fDigitsArr) ;
+    if(rp->GetMultiplicity()==0){
+      fEMCRecPoints->RemoveAt(index) ;
       delete rp ;
+      continue;
     }
+
+    // No vertex is available now, calculate corrections in PID
+    rp->EvalAll(fDigitsArr) ;
+    rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
+    rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
+    rp->EvalLocal2TrackingCSTransform();
   }
-  emcRecPoints->Compress() ;
-  emcRecPoints->Sort() ;
-  //  emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
-  for(index = 0; index < emcRecPoints->GetEntries(); index++){
-    dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
+  fEMCRecPoints->Compress() ;
+  fEMCRecPoints->Sort() ; 
+  //  fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
+  for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
+    static_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++){
-    AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
-    rp->EvalAll(fW0CPV,digits) ;
+  for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
+    AliPHOSCpvRecPoint * rp = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
+    rp->EvalAll(fDigitsArr) ;
+    rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
+    rp->EvalLocal2TrackingCSTransform();
   }
-  cpvRecPoints->Sort() ;
+  fCPVRecPoints->Sort() ;
   
-  for(index = 0; index < cpvRecPoints->GetEntries(); index++)
-    dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
+  for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
+    static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
   
-  cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
+  fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
   
   if(fWrite){ //We write TreeR
-    TTree * treeR = gime->TreeR();
-    
-    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());
-    
-    emcBranch        ->Fill() ;
-    cpvBranch        ->Fill() ;
-    
-    gime->WriteRecPoints("OVERWRITE");
-    gime->WriteClusterizer("OVERWRITE");
+    fTreeR->Fill();
   }
 }
 
@@ -623,124 +569,101 @@ 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::Instance();
-
-  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
-  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
-
-  emcRecPoints->Delete() ;
-  cpvRecPoints->Delete() ;
-  
-  TClonesArray * digits = gime->Digits() ; 
-
-  //Remove digits below threshold
-  CleanDigits(digits) ;
+  fNumberOfCpvClusters     = 0 ;
+  fNumberOfEmcClusters     = 0 ;
 
-  TClonesArray * digitsC =  static_cast<TClonesArray*>( digits->Clone() ) ;
-  // Clusterization starts  
+  //Mark all digits as unused yet
+  const Int_t maxNDigits = 3584; // There is no clusters larger than PHOS module ;)
+  Int_t nDigits=fDigitsArr->GetEntriesFast() ;
 
-  TIter nextdigit(digitsC) ; 
+  for(Int_t i=0; i<nDigits; i++){
+    fDigitsUsed[i]=0 ;
+  }
+  Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
+                     //e.g. first digit in this module, first CPV digit etc.
   AliPHOSDigit * digit ; 
-  Bool_t notremoved = kTRUE ;
+  TArrayI clusterdigitslist(maxNDigits) ;   
+  AliPHOSRecPoint * clu = 0 ; 
+  for(Int_t i=0; i<nDigits; i++){
+    if(fDigitsUsed[i])
+      continue ;
 
-  while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
-    
+    digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
 
-    AliPHOSRecPoint * clu = 0 ; 
+    clu=0 ;
 
-    TArrayI clusterdigitslist(1500) ;   
     Int_t index ;
 
-    if (( IsInEmc (digit) &&
-         CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) || 
-        ( IsInCpv (digit) &&
-         CalibrateCPV(digit->GetAmp()   ,digit->GetId()) > fCpvClusteringThreshold ) ) {
+    //is this digit so energetic that start cluster?
+    if (( IsInEmc(digit) &&  Calibrate(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) || 
+        ( IsInCpv(digit) &&  Calibrate(digit->GetEnergy(),digit->GetId()) > fCpvClusteringThreshold ) ) {
       Int_t iDigitInCluster = 0 ; 
-      
       if  ( IsInEmc(digit) ) {   
         // start a new EMC RecPoint
-        if(fNumberOfEmcClusters >= emcRecPoints->GetSize()) 
-          emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
+        if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize()) 
+          fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
           
-        emcRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
-        clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ; 
+        fEMCRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
+        clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ; 
        fNumberOfEmcClusters++ ; 
-       clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
+       clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),CalibrateT(digit->GetTime(),digit->GetId(),digit->IsLG())) ;
         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
         iDigitInCluster++ ;
-        digitsC->Remove(digit) ; 
-
+        fDigitsUsed[i]=kTRUE ; 
       } else { 
-        
         // start a new CPV cluster
-        if(fNumberOfCpvClusters >= cpvRecPoints->GetSize()) 
-          cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
-
-        cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
+        if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize()) 
+          fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
 
-        clu =  dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;  
+        fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
+        clu =  static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;  
         fNumberOfCpvClusters++ ; 
-        clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;        
+        clu->AddDigit(*digit,  Calibrate(digit->GetEnergy(),digit->GetId()),0.) ; // no timing information in CPV
         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;        
         iDigitInCluster++ ; 
-        digitsC->Remove(digit) ; 
-        nextdigit.Reset() ;
-        
-        // Here we remove remaining EMC digits, which cannot make a cluster
-        
-        if( notremoved ) { 
-          while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
-            if( IsInEmc(digit) ) 
-              digitsC->Remove(digit) ;
-            else 
-              break ;
-          }
-          notremoved = kFALSE ;
-        }
-        
+        fDigitsUsed[i]=kTRUE ;
       } // else        
-      
-      nextdigit.Reset() ;
-      
+
+      //Now scan remaining digits in list to find neigbours of our seed
       AliPHOSDigit * digitN ; 
       index = 0 ;
       while (index < iDigitInCluster){ // scan over digits already in cluster 
-        digit =  dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) )  ;      
+        digit =  static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) )  ;      
         index++ ; 
-        while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits 
+        for(Int_t j=iFirst; j<nDigits; j++){
+         if (iDigitInCluster >= maxNDigits) {
+           AliError(Form("The number of digits in cluster is more than %d, skip the rest of event",
+                         maxNDigits));
+           return;
+         }
+          if(fDigitsUsed[j]) 
+            continue ;        //look through remaining digits
+          digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
           switch (ineb ) {
+          case -1:   //too early (e.g. previous module), do not look before j at subsequent passes
+            iFirst=j ;
+            break ;
           case 0 :   // not a neighbour
             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() ; 
+           clu->AddDigit(*digitN, Calibrate(digitN->GetEnergy(),digitN->GetId()),CalibrateT(digitN->GetTime(),digitN->GetId(),digit->IsLG())) ;
+            clusterdigitslist[iDigitInCluster] = j ; 
             iDigitInCluster++ ; 
-            digitsC->Remove(digitN) ;
+            fDigitsUsed[j]=kTRUE ;
             break ;
           case 2 :   // too far from each other
-            goto endofloop;   
+            goto endOfLoop;   
           } // switch
           
-        } // while digitN
+        }
         
-      endofloop: ;
-        nextdigit.Reset() ; 
+        endOfLoop: ; //scanned all possible neighbours for this digit
         
       } // loop over cluster     
-
     } // energy theshold  
-
-    
-  } // while digit
-
-  delete digitsC ;
+  }
 
 }
 
@@ -750,36 +673,29 @@ 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() ;
-
-  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
-  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
-  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 = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ;
+      AliPHOSEmcRecPoint * emcRecPoint = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
         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) ;
-        emcRecPoints->Remove(emcRecPoint); 
-        emcRecPoints->Compress() ;
+
+        fEMCRecPoints->Remove(emcRecPoint); 
+        fEMCRecPoints->Compress() ;
         index-- ;
         fNumberOfEmcClusters -- ;
         numberofNotUnfolded-- ;
@@ -798,28 +714,28 @@ 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 = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ;
+      AliPHOSRecPoint * recPoint = static_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
 
       if(recPoint->GetPHOSMod()> nModulesToUnfold)
         break ;
       
-      AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ; 
+      AliPHOSEmcRecPoint * emcRecPoint = static_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() ;
+        fCPVRecPoints->Remove(emcRecPoint); 
+        fCPVRecPoints->Compress() ;
         index-- ;
         numberofCpvNotUnfolded-- ;
         fNumberOfCpvClusters-- ;
@@ -834,13 +750,16 @@ 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 ;
 }
@@ -853,18 +772,11 @@ void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 { 
   // 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) ;
@@ -878,18 +790,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 = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ;   
-    geom->AbsToRelNumbering(digit->GetId(), relid) ;
-    geom->RelPosInModule(relid, xDigit, zDigit) ;
+    digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;   
+    fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+    fGeom->RelPosInModule(relid, xDigit, zDigit) ;
     efit[iDigit] = 0;
 
     iparam = 0 ;    
@@ -898,13 +812,12 @@ 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) ;
     }
   }
   
-
   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
   // so that energy deposited in each cell is distributed betwin new clusters proportionally
   // to its contribution to efit
@@ -918,38 +831,38 @@ 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 = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
+      (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
+      emcRP = static_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 = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
+      (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
+      emcRP = static_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
       fNumberOfCpvClusters++ ;
     }
     
     Float_t eDigit ;
     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
-      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] ; 
+      digit = static_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 ) ;
+      emcRP->AddDigit( *digit, eDigit,CalibrateT(digit->GetTime(),digit->GetId(),digit->IsLG()) ) ;
     }        
   }
  
@@ -964,14 +877,17 @@ 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 = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
-
-  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
-  TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
+  TList * toMinuit = static_cast<TList*>( gMinuit->GetObjectFit() ) ;
 
+  AliPHOSEmcRecPoint * emcRP = static_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
+  TClonesArray * digits = static_cast<TClonesArray*>( toMinuit->At(1) )  ;
+  // A bit buggy way to get an access to the geometry
+  // To be revised!
+  AliPHOSGeometry *geom = static_cast<AliPHOSGeometry *>(toMinuit->At(2));
 
+//  TVector3 * vtx = static_cast<TVector3*>(toMinuit->At(3)) ;  //Vertex position
   
-  //  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
+  //  AliPHOSEmcRecPoint * emcRP = static_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
 
   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 
@@ -979,7 +895,7 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou
 
   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
   
-  const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; 
+//  TVector3 vInc ;
   fret = 0. ;     
   Int_t iparam ;
 
@@ -994,7 +910,7 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou
 
   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 
-    digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ); 
+    digit = static_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ); 
 
     Int_t relid[4] ;
     Float_t xDigit ;
@@ -1008,12 +924,13 @@ 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]) ;
+         Double_t dx = (xDigit - x[iParam]) ;
          iParam++ ; 
-         distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ; 
-         distance = TMath::Sqrt( distance ) ; 
+         Double_t dz = (zDigit - x[iParam]) ; 
          iParam++ ;          
-         efit += x[iParam] * ShowerShape(distance) ;
+//         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 = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E) 
@@ -1022,8 +939,11 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou
          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(dr) ;
+//         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) ) +
@@ -1045,9 +965,9 @@ 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) ;
      }
 
      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ; 
@@ -1086,7 +1006,7 @@ void AliPHOSClusterizerv1::Print(const Option_t *)const
   else 
     message = " AliPHOSClusterizerv1 not initialized " ;
   
-  AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),  
+  AliInfo(Form("%s, %s %s %s %s %f %f %f %f %f %f", message.Data(),  
        taskName.Data(), 
        GetTitle(),
        taskName.Data(), 
@@ -1098,33 +1018,21 @@ void AliPHOSClusterizerv1::Print(const Option_t *)const
        fCpvLocMaxCut, 
        fW0CPV )) ; 
 }
-
-
 //____________________________________________________________________________
 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
 {
   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
 
-  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() ))  ;
+  AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints", 
+              fEMCRecPoints->GetEntriesFast(),
+              fCPVRecPoints->GetEntriesFast() ))  ;
  
-  fRecPointsInRun +=  emcRecPoints->GetEntriesFast() ; 
-  fRecPointsInRun +=  cpvRecPoints->GetEntriesFast() ; 
-  
-  
   if(strstr(option,"all")) {
     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]; 
@@ -1145,8 +1053,8 @@ void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
     //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) ; 
+    for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
+      AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ; 
       
       TVector3  locpos;  
       rp->GetLocalPosition(locpos);
@@ -1158,3 +1066,99 @@ void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
   }
 }
 
+
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::SetDistancesToBadChannels()
+{
+  //For each EMC rec. point set the distance to the nearest bad crystal.
+  //Author: Boris Polichtchouk 
+
+  if(!fgCalibData->GetNumOfEmcBadChannels()) return;
+
+  Int_t badIds[8000];
+  memset(badIds,0,8000*sizeof(Int_t));
+  fgCalibData->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;
+  Int_t relid[4]={0,0,0,0} ;
+  TVector3 lpos ;
+  for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
+    rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
+    //evaluate distance to border
+    relid[0]=rp->GetPHOSMod() ;
+    relid[2]=1 ;
+    relid[3]=1 ;
+    Float_t xcorner,zcorner;
+    fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
+    rp->GetLocalPosition(lpos) ;
+    minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
+    for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
+      fGeom->AbsToRelNumbering(badIds[iBad],relid)  ;
+      if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since 
+        continue ;                   //bad channels can be in the module which does not exist in simulations.
+      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); 
+  }
+
+}
+//==================================================================================
+Float_t AliPHOSClusterizerv1::Calibrate(Float_t amp, Int_t absId) const{
+  // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB
+
+  const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
+
+  //Determine rel.position of the cell absolute ID
+  Int_t relId[4];
+  geom->AbsToRelNumbering(absId,relId);
+  Int_t module=relId[0];
+  Int_t row   =relId[2];
+  Int_t column=relId[3];
+  if(relId[1]){ //CPV
+    Float_t calibration = fgCalibData->GetADCchannelCpv(module,column,row);
+    return amp*calibration ;
+  }   
+  else{ //EMC
+    Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row);
+    return amp*calibration ;
+  }
+}
+//==================================================================================
+Float_t AliPHOSClusterizerv1::CalibrateT(Float_t time, Int_t absId,Bool_t isLG)const{
+  // Calibrate time in EMC digit
+
+  const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
+
+  //Determine rel.position of the cell absolute ID
+  Int_t relId[4];
+  geom->AbsToRelNumbering(absId,relId);
+  Int_t module=relId[0];
+  Int_t row   =relId[2];
+  Int_t column=relId[3];
+  if(relId[1]){ //CPV
+    return 0. ;
+  }
+  else{ //EMC
+    if(isLG)
+      time += fgCalibData->GetLGTimeShiftEmc(module,column,row);
+    else
+      time += fgCalibData->GetTimeShiftEmc(module,column,row);
+    return time ;
+  }
+}
+