]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSClusterizerv1.cxx
Coding convensions satisfied (T.Pocheptsov)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSClusterizerv1.cxx
index d47ad3855f02189cf804e4caedc1bf246dad3203..73a91a6f70aa29ac7f38ab896b04342389e547b5 100644 (file)
 /* History of cvs commits:
  *
  * $Log$
+ * 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
  *
 
 // --- AliRoot header files ---
 #include "AliLog.h"
+#include "AliRunLoader.h"
+#include "AliGenerator.h"
 #include "AliPHOSGetter.h"
 #include "AliPHOSGeometry.h" 
 #include "AliPHOSClusterizerv1.h"
 #include "AliPHOSDigit.h"
 #include "AliPHOSDigitizer.h"
 #include "AliPHOSCalibrationDB.h"
+#include "AliCDBManager.h"
 #include "AliCDBStorage.h"
-#include "AliCDBLocal.h"
+#include "AliCDBEntry.h"
 
 ClassImp(AliPHOSClusterizerv1)
   
 //____________________________________________________________________________
-  AliPHOSClusterizerv1::AliPHOSClusterizerv1() : AliPHOSClusterizer()
+AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
+  AliPHOSClusterizer(),
+  fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
+  fWrite(0),                  fNumberOfEmcClusters(0),  fNumberOfCpvClusters(0),
+  fCalibData(0),              fADCchanelEmc(0),         fADCpedestalEmc(0),
+  fADCchanelCpv(0),           fADCpedestalCpv(0),       fEmcClusteringThreshold(0),
+  fCpvClusteringThreshold(0), fEmcMinE(0),              fCpvMinE(0),
+  fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
+  fW0CPV(0),                  fRecPointsInRun(0),       fEmcTimeGate(0),
+  fIsOldRCUFormat(0)
 {
   // default ctor (to be used mainly by Streamer)
   
@@ -97,8 +172,16 @@ ClassImp(AliPHOSClusterizerv1)
 }
 
 //____________________________________________________________________________
-AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName)
-:AliPHOSClusterizer(alirunFileName, eventFolderName)
+AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName) :
+  AliPHOSClusterizer(alirunFileName, eventFolderName),
+  fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
+  fWrite(0),                  fNumberOfEmcClusters(0),  fNumberOfCpvClusters(0),
+  fCalibData(0),              fADCchanelEmc(0),         fADCpedestalEmc(0),
+  fADCchanelCpv(0),           fADCpedestalCpv(0),       fEmcClusteringThreshold(0),
+  fCpvClusteringThreshold(0), fEmcMinE(0),              fCpvMinE(0),
+  fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
+  fW0CPV(0),                  fRecPointsInRun(0),       fEmcTimeGate(0),
+  fIsOldRCUFormat(0)
 {
   // ctor with the indication of the file where header Tree and digits Tree are stored
   
@@ -107,6 +190,20 @@ AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const T
   fDefaultInit = kFALSE ; 
 }
 
+//____________________________________________________________________________
+AliPHOSClusterizerv1::AliPHOSClusterizerv1(const AliPHOSClusterizerv1 & obj) :
+  AliPHOSClusterizer(obj),
+  fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
+  fWrite(0),                  fNumberOfEmcClusters(0),  fNumberOfCpvClusters(0),
+  fCalibData(0),              fADCchanelEmc(0),         fADCpedestalEmc(0),
+  fADCchanelCpv(0),           fADCpedestalCpv(0),       fEmcClusteringThreshold(0),
+  fCpvClusteringThreshold(0), fEmcMinE(0),              fCpvMinE(0),
+  fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
+  fW0CPV(0),                  fRecPointsInRun(0),       fEmcTimeGate(0),
+  fIsOldRCUFormat(0)
+{
+  // Copy constructor
+}
 //____________________________________________________________________________
   AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
 {
@@ -120,9 +217,9 @@ const TString AliPHOSClusterizerv1::BranchName() const
 }
  
 //____________________________________________________________________________
-Float_t  AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId)
+Float_t  AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
 {  
-  // Convert digitized amplitude into energy.
+  // 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.
 
@@ -133,20 +230,43 @@ Float_t  AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId)
     Int_t   module = relId[0];
     Int_t   column = relId[3];
     Int_t   row    = relId[2];
-    if(absId <= fEmcCrystals) { //calibrate as EMC 
+    if(absId <= fEmcCrystals) { // this is EMC 
       fADCchanelEmc   = fCalibData->GetADCchannelEmc (module,column,row);
-      fADCpedestalEmc = fCalibData->GetADCpedestalEmc(module,column,row);
-      return fADCpedestalEmc + amp*fADCchanelEmc ;        
+      return amp*fADCchanelEmc ;        
     }
-    else //calibrate as CPV, not implemented yet
-      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];
+    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;
 }
 
 //____________________________________________________________________________
@@ -183,16 +303,23 @@ void AliPHOSClusterizerv1::Exec(Option_t *option)
   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
     if (fRawReader == 0)
       gime->Event(ievent    ,"D"); // Read digits from simulated data
-    else
-      gime->Event(fRawReader,"W"); // Read digits from raw data
-    
+    else {
+      AliRunLoader * rl = AliRunLoader::GetRunLoader(gime->PhosLoader()->GetTitle());
+      rl->GetEvent(ievent);
+      gime->Event(fRawReader,"W",fIsOldRCUFormat); // Read digits from raw data
+    }
     fNumberOfEmcClusters  = fNumberOfCpvClusters  = 0 ;
     
     MakeClusters() ;
+    
+    AliDebug(2,Form(" ---- Printing clusters (%d) of event %d.\n",
+          gime->EmcRecPoints()->GetEntries(),ievent));
+    if(AliLog::GetGlobalDebugLevel()>1)
+      gime->EmcRecPoints()->Print();
 
     if(fToUnfold)             
       MakeUnfolding() ;
-
+    
     WriteRecPoints();
 
     if(strstr(option,"deb"))  
@@ -201,7 +328,7 @@ void AliPHOSClusterizerv1::Exec(Option_t *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();
@@ -225,7 +352,6 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit **
   
   AliPHOSGetter * gime = AliPHOSGetter::Instance();
   TClonesArray * digits = gime->Digits(); 
-  
 
   gMinuit->mncler();                     // Reset Minuit's list of paramters
   gMinuit->SetPrintLevel(-1) ;           // No Printout
@@ -313,28 +439,17 @@ void AliPHOSClusterizerv1::GetCalibrationParameters()
 {
   // Set calibration parameters:
   // if calibration database exists, they are read from database,
-  // otherwise, they are taken from digitizer.
+  // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
   //
-  // It is a user responsilibity to open CDB before reconstruction:
-  // AliCDBLocal *loc = new AliCDBLocal("CalibDB");
+  // 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(-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");
 
-  if(AliCDBStorage::Instance())
-    fCalibData = (AliPHOSCalibData*)AliCDBStorage::Instance()
-      ->Get("PHOS/Calib/GainFactors_and_Pedestals",gAlice->GetRunNumber());
-
-  if(!fCalibData)
-    {
-      if ( !gime->Digitizer() ) 
-       gime->LoadDigitizer();
-      AliPHOSDigitizer * dig = gime->Digitizer(); 
-      fADCchanelEmc   = dig->GetEMCchannel() ;
-      fADCpedestalEmc = dig->GetEMCpedestal();
-    
-      fADCchanelCpv   = dig->GetCPVchannel() ;
-      fADCpedestalCpv = dig->GetCPVpedestal() ; 
-    }
 }
 
 //____________________________________________________________________________
@@ -389,6 +504,8 @@ void AliPHOSClusterizerv1::InitParameters()
   fCalibData               = 0 ;
 
   SetEventRange(0,-1) ;
+
+  fIsOldRCUFormat          = kFALSE;
 }
 
 //____________________________________________________________________________
@@ -435,11 +552,14 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)c
   return rv ; 
 }
 //____________________________________________________________________________
-void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits){
+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)) ;
-    Float_t cut = IsInEmc(digit) ? fEmcMinE : fCpvMinE ;
-    if(Calibrate(digit->GetAmp(),digit->GetId()) < cut)
+    if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
+        (IsInCpv(digit) && CalibrateCPV(digit->GetAmp()   ,digit->GetId()) < fCpvMinE) )
       digits->RemoveAt(i) ;
   }
   digits->Compress() ;
@@ -448,6 +568,13 @@ void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits){
     digit->SetIndexInList(i) ;     
   }
 
+  //Overwrite digits tree
+  AliPHOSGetter* gime = AliPHOSGetter::Instance();
+  TTree * treeD = gime->TreeD();
+  treeD->Branch("PHOS", &digits);
+  treeD->Fill() ;
+  gime->WriteDigits("OVERWRITE");
+  gime->PhosLoader()->UnloadDigits() ;
 }
 //____________________________________________________________________________
 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
@@ -501,33 +628,38 @@ void AliPHOSClusterizerv1::WriteRecPoints()
   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
   TClonesArray * digits = gime->Digits() ; 
  
   Int_t index ;
   //Evaluate position, dispersion and other RecPoint properties..
   Int_t nEmc = emcRecPoints->GetEntriesFast();
   for(index = 0; index < nEmc; index++){
     AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
     rp->Purify(fEmcMinE) ;
-    if(rp->GetMultiplicity()>0.) //If this RP is not empty
-      rp->EvalAll(fW0,digits) ;
-    else{
+    if(rp->GetMultiplicity()==0){
       emcRecPoints->RemoveAt(index) ;
       delete rp ;
     }
+
+// No vertex is available now, calculate cirrections in PID
+      rp->EvalAll(fW0,digits) ;
+      TVector3 fakeVtx(0.,0.,0.) ;
+      rp->EvalAll(fW0,fakeVtx,digits) ;
   }
   emcRecPoints->Compress() ;
-  emcRecPoints->Sort() ;
+//  emcRecPoints->Sort() ; //Can not sort until position is calculated!
   //  emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
   for(index = 0; index < emcRecPoints->GetEntries(); index++){
     dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->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) ;
   }
-  cpvRecPoints->Sort() ;
+//  cpvRecPoints->Sort() ;
   
   for(index = 0; index < cpvRecPoints->GetEntries(); index++)
     dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
@@ -576,9 +708,7 @@ void AliPHOSClusterizerv1::MakeClusters()
   //Remove digits below threshold
   CleanDigits(digits) ;
 
-
   TClonesArray * digitsC =  static_cast<TClonesArray*>( digits->Clone() ) ;
-  
  
   // Clusterization starts  
 
@@ -594,8 +724,10 @@ void AliPHOSClusterizerv1::MakeClusters()
     TArrayI clusterdigitslist(1500) ;   
     Int_t index ;
 
-    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) ) {   
@@ -605,10 +737,10 @@ void AliPHOSClusterizerv1::MakeClusters()
           
         emcRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
         clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ; 
-          fNumberOfEmcClusters++ ; 
-        clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId())) ; 
-        clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;        
-        iDigitInCluster++ ; 
+       fNumberOfEmcClusters++ ; 
+       clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
+        clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
+        iDigitInCluster++ ;
         digitsC->Remove(digit) ; 
 
       } else { 
@@ -621,7 +753,7 @@ void AliPHOSClusterizerv1::MakeClusters()
 
         clu =  dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;  
         fNumberOfCpvClusters++ ; 
-        clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId()) ) ;        
+        clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;        
         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;        
         iDigitInCluster++ ; 
         digitsC->Remove(digit) ; 
@@ -654,7 +786,11 @@ void AliPHOSClusterizerv1::MakeClusters()
           case 0 :   // not a neighbour
             break ;
           case 1 :   // are neighbours 
-            clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->GetId() ) ) ;
+           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) ;
@@ -769,13 +905,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 ;
 }
@@ -813,12 +952,14 @@ 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 ++){
@@ -833,9 +974,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) ;
+//      geom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
+//      efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
+      efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
     }
   }
   
@@ -853,6 +994,7 @@ void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
     zpar = fitparameters[iparam+1] ;
     epar = fitparameters[iparam+2] ;
     iparam += 3 ;    
+//    geom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
     
     AliPHOSEmcRecPoint * emcRP = 0 ;  
 
@@ -880,9 +1022,8 @@ void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
       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] ; 
+//      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 ) ;
     }        
@@ -903,8 +1044,7 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou
 
   AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
   TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
-
-
+//  TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(2)) ;  //Vertex position
   
   //  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
 
@@ -915,6 +1055,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 ;
 
@@ -943,12 +1084,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) ;
+//         geom->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) 
@@ -957,8 +1099,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] ;
+//         geom->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) ) +
@@ -980,9 +1125,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) ;
+//       geom->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] ; 
@@ -1033,8 +1178,27 @@ void AliPHOSClusterizerv1::Print(const Option_t *)const
        fCpvLocMaxCut, 
        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)
 {
@@ -1093,3 +1257,49 @@ 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(!fCalibData->GetNumOfEmcBadChannels()) return;
+  AliInfo(Form("%d bad channel(s) found.\n",fCalibData->GetNumOfEmcBadChannels()));
+
+  AliPHOSGetter* gime = AliPHOSGetter::Instance();
+  AliPHOSGeometry* geom = gime->PHOSGeometry();
+  TObjArray * emcRecPoints = gime->EmcRecPoints() ;
+  
+  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<emcRecPoints->GetEntries(); iRP++){
+    rp = (AliPHOSEmcRecPoint*)emcRecPoints->At(iRP);
+    minDist = 1.e+07;
+
+    for(Int_t iBad=0; iBad<fCalibData->GetNumOfEmcBadChannels(); iBad++) {
+      rp->GetGlobalPosition(gposRecPoint,gmat);
+      geom->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); 
+  }
+
+}