]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSClusterizerv1.cxx
Correcting compiler warnings
[u/mrichter/AliRoot.git] / PHOS / AliPHOSClusterizerv1.cxx
index 534c9f86334067e126a6a568179afd51283f1136..744fe1b8b177bde235dda3f646770792f4c3d1e8 100644 (file)
 /* History of cvs commits:
  *
  * $Log$
+ * 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.
  *
 
 // --- 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 "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)
   
@@ -83,8 +163,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
   
@@ -93,6 +181,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()
 {
@@ -106,18 +208,56 @@ const TString AliPHOSClusterizerv1::BranchName() const
 }
  
 //____________________________________________________________________________
-Float_t  AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId) const
+Float_t  AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
 {  
-  //To be replaced later by the method, reading individual parameters from the database
-  if(fCalibrationDB){
-    return fCalibrationDB->Calibrate(amp,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) //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;
 }
 
 //____________________________________________________________________________
@@ -136,8 +276,13 @@ void AliPHOSClusterizerv1::Exec(Option_t *option)
     return ;
   }
 
-  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
   GetCalibrationParameters() ;
+
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  if (fRawReader == 0)
+    gime->SetRawDigits(kFALSE);
+  else
+    gime->SetRawDigits(kTRUE);
   
   if (fLastEvent == -1) 
     fLastEvent = gime->MaxEvent() - 1 ;
@@ -149,12 +294,19 @@ 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() ;
@@ -167,7 +319,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();
@@ -191,7 +343,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
@@ -277,21 +428,29 @@ 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();
-  if(gime->IsRawDigits()){
-    fCalibrationDB = gime->CalibrationDB();    
-  }
-  else{
-    if ( !gime->Digitizer() ) 
-      gime->LoadDigitizer();
-    AliPHOSDigitizer * dig = gime->Digitizer(); 
-    fADCchanelEmc   = dig->GetEMCchannel() ;
-    fADCpedestalEmc = dig->GetEMCpedestal();
+  // 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() ; 
-  }  
+      fADCchanelCpv   = dig->GetCPVchannel() ;
+      fADCpedestalCpv = dig->GetCPVpedestal() ; 
+    }
 }
 
 //____________________________________________________________________________
@@ -343,9 +502,11 @@ void AliPHOSClusterizerv1::InitParameters()
 
   fWrite                   = kTRUE ;
 
-  fCalibrationDB           = 0 ;
+  fCalibData               = 0 ;
 
   SetEventRange(0,-1) ;
+
+  fIsOldRCUFormat          = kFALSE;
 }
 
 //____________________________________________________________________________
@@ -392,11 +553,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() ;
@@ -405,6 +569,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
@@ -458,22 +629,24 @@ 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) ;
@@ -484,7 +657,7 @@ void AliPHOSClusterizerv1::WriteRecPoints()
     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) ;
@@ -533,9 +706,7 @@ void AliPHOSClusterizerv1::MakeClusters()
   //Remove digits below threshold
   CleanDigits(digits) ;
 
-
   TClonesArray * digitsC =  static_cast<TClonesArray*>( digits->Clone() ) ;
-  
  
   // Clusterization starts  
 
@@ -551,8 +722,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) ) {   
@@ -562,10 +735,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 { 
@@ -578,7 +751,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) ; 
@@ -611,7 +784,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) ;
@@ -726,13 +903,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 ;
 }
@@ -770,12 +950,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 ++){
@@ -790,9 +972,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) ;
     }
   }
   
@@ -810,6 +992,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 ;  
 
@@ -837,9 +1020,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 ) ;
     }        
@@ -860,8 +1042,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
 
@@ -872,6 +1053,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 ;
 
@@ -900,12 +1082,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) 
@@ -914,8 +1097,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) ) +
@@ -937,9 +1123,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] ; 
@@ -990,8 +1176,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)
 {