]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSTrackSegmentMakerv1.cxx
new amplitude calibration up to 50 MIPs
[u/mrichter/AliRoot.git] / PHOS / AliPHOSTrackSegmentMakerv1.cxx
index f5e1f5204c5aa96a11405d1384a34d38aec2c939..8135da99385ab17e2f28a6e325d46a99a9d25d6c 100644 (file)
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.93  2007/10/10 09:05:10  schutz
+ * Changing name QualAss to QA
+ *
+ * Revision 1.92  2007/08/28 12:55:08  policheh
+ * Loaders removed from the reconstruction code (C.Cheshkov)
+ *
+ * Revision 1.91  2007/08/07 14:12:03  kharlov
+ * Quality assurance added (Yves Schutz)
+ *
+ * Revision 1.90  2007/07/11 13:43:30  hristov
+ * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
+ *
+ * Revision 1.89  2007/07/03 08:13:04  kharlov
+ * Bug fix in CPV local coordinates
+ *
+ * Revision 1.88  2007/06/27 09:11:07  kharlov
+ * Bug fix for CPV-EMC distance
+ *
+ * Revision 1.87  2007/05/04 14:49:29  policheh
+ * AliPHOSRecPoint inheritance from AliCluster
+ *
+ * Revision 1.86  2007/04/02 15:00:16  cvetan
+ * No more calls to gAlice in the reconstruction
+ *
+ * Revision 1.85  2007/03/28 19:18:15  kharlov
+ * RecPoints recalculation in TSM removed
+ *
+ * Revision 1.84  2007/03/07 07:01:21  hristov
+ * Fixing copy/paste erro. Additional protections
+ *
+ * Revision 1.83  2007/03/06 21:07:37  kharlov
+ * DP: xz CPV-EMC distance filled to TS
+ *
+ * Revision 1.82  2007/03/06 06:54:48  kharlov
+ * DP:Calculation of cluster properties dep. on vertex added
+ *
+ * Revision 1.81  2007/02/05 10:02:40  kharlov
+ * Module numbering is corrected
+ *
+ * Revision 1.80  2006/08/28 10:01:56  kharlov
+ * Effective C++ warnings fixed (Timur Pocheptsov)
+ *
+ * Revision 1.79  2006/04/25 12:41:15  hristov
+ * Moving non-persistent data to AliESDfriend (Yu.Belikov)
+ *
+ * Revision 1.78  2005/11/18 13:04:51  hristov
+ * Bug fix
+ *
+ * Revision 1.77  2005/11/17 23:34:36  hristov
+ * Corrected logics
+ *
+ * Revision 1.76  2005/11/17 22:29:12  hristov
+ * Faster version, no attempt to match tracks outside the PHOS acceptance
+ *
+ * Revision 1.75  2005/11/17 12:35:27  hristov
+ * Use references instead of objects. Avoid to create objects when they are not really needed
+ *
+ * Revision 1.74  2005/07/08 14:01:36  hristov
+ * Tracking in non-uniform nmagnetic field (Yu.Belikov)
+ *
+ * Revision 1.73  2005/05/28 14:19:05  schutz
+ * Compilation warnings fixed by T.P.
+ *
+ */
 
 //_________________________________________________________________________
-// A brief description of the class
-//*-- Author : Yves Schutz  SUBATECH 
-//////////////////////////////////////////////////////////////////////////////
+// Implementation version 1 of algorithm class to construct PHOS track segments
+// Track segment for PHOS is list of 
+//        EMC RecPoint + (possibly) CPV RecPoint
+// To find TrackSegments we do the following: 
+//  for each EMC RecPoints we look at
+//   CPV RecPoints in the radious fRcpv. 
+//  If there is such a CPV RecPoint, 
+//   we make "Link" it is just indexes of EMC and CPV RecPoint and distance
+//   between them in the PHOS plane. 
+//  Then we sort "Links" and starting from the 
+//   least "Link" pointing to the unassined EMC and CPV RecPoints assing them to 
+//   new TrackSegment. 
+// If there is no CPV RecPoint we make TrackSegment 
+// consisting from EMC alone. There is no TrackSegments without EMC RecPoint.
+//// In principle this class should be called from AliPHOSReconstructor, but 
+// one can use it as well in standalone mode.
+// Use  case:
+//  root [0] AliPHOSTrackSegmentMakerv1 * t = new AliPHOSTrackSegmentMaker("galice.root", "tracksegmentsname", "recpointsname")
+//  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
+//               // reads gAlice from header file "galice.root", uses recpoints stored in the branch names "recpointsname" (default = "Default")
+//               // and saves recpoints in branch named "tracksegmentsname" (default = "recpointsname")                       
+//  root [1] t->ExecuteTask()
+//  root [3] t->SetTrackSegmentsBranch("max distance 5 cm")
+//  root [4] t->ExecuteTask("deb all time") 
+//                 
+//*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH) 
+//
 
 // --- ROOT system ---
-
-#include "TObjArray.h"
-#include "TClonesArray.h"
-#include "TMinuit.h"
+#include "TVector3.h"
+#include "TTree.h"
+#include "TBenchmark.h"
 
 // --- Standard library ---
-
-#include "iostream.h"
-
+#include "Riostream.h"
 // --- AliRoot header files ---
-
+#include "AliPHOSGeometry.h"
 #include "AliPHOSTrackSegmentMakerv1.h"
 #include "AliPHOSTrackSegment.h"
 #include "AliPHOSLink.h"
-#include "AliPHOSv0.h"
-#include "AliRun.h"
+#include "AliESDEvent.h"
+#include "AliESDtrack.h"
+#include "AliPHOSEmcRecPoint.h"
+#include "AliPHOSCpvRecPoint.h"
+#include "AliLog.h"
+#include "AliMagF.h"
+#include "AliMagF.h"
+#include "AliTracker.h"
+#include "AliGeomManager.h"
+#include "AliCluster.h"
+#include "AliKalmanTrack.h"
+#include "AliGlobalQADataMaker.h"
+#include "AliVParticle.h"
+
 
 ClassImp( AliPHOSTrackSegmentMakerv1) 
 
 
 //____________________________________________________________________________
- AliPHOSTrackSegmentMakerv1:: AliPHOSTrackSegmentMakerv1() // ctor
+AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() :
+  AliPHOSTrackSegmentMaker(),
+  fDefaultInit(kTRUE),
+  fWrite(kFALSE),
+  fNTrackSegments(0),
+  fRcpv(0.f),
+  fRtpc(0.f),
+  fVtx(0.f), 
+  fLinkUpArray(0),
+  fEmcFirst(0),
+  fEmcLast(0),
+  fCpvFirst(0),
+  fCpvLast(0),
+  fModule(0),
+  fTrackSegments(NULL)
 {
-  fR0 = 4. ;   
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-  //clusters are sorted in "rows" and "columns" of width geom->GetCrystalSize(0),
-  fDelta = fR0 + geom->GetCrystalSize(0) ;
+  // default ctor (to be used mainly by Streamer)
+  InitParameters() ; 
 }
 
 //____________________________________________________________________________
-Bool_t  AliPHOSTrackSegmentMakerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy,
-                                   Int_t NPar, Float_t * FitParametres)
-{ //Calls TMinuit for fitting cluster with several maxima 
-
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-  TMinuit *gMinuit = new TMinuit(NPar);  //initialize TMinuit with a maximum of 5 params
-  gMinuit->SetPrintLevel(-1) ;           //No PRIntout
-  gMinuit->SetFCN(UnfoldingChiSquare );  //To set the address of the minimization function 
-  gMinuit->SetObjectFit(emcRP) ;   //To tranfer pointer to UnfoldingChiSquare
-
-  //filling initial values for fit parameters
-  AliPHOSDigit * digit ;
-  Int_t ierflg = 0; 
-  Int_t index = 0 ;
-  Int_t NDigits = (Int_t) NPar / 3 ;
-  Int_t iDigit ;
-  for(iDigit = 0 ; iDigit < NDigits ; iDigit ++){
-    digit = (AliPHOSDigit *) maxAt[iDigit]; 
-
-    Int_t RelId[4] ;
-    Float_t x ;
-    Float_t z ;
-    geom->AbsToRelNumbering(digit->GetId(),RelId) ;
-    geom->RelPosInModule(RelId,x,z) ;
-
-    Float_t Energy = maxAtEnergy[iDigit] ;
-
-    gMinuit->mnparm(index++, " ",  x, 0.1, 0, 0, ierflg) ;  
-    if(ierflg != 0){ 
-      cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : x = " << x << endl ;
-      return kFALSE;
-    }
-    gMinuit->mnparm(index++, " ",  z, 0.1, 0, 0, ierflg) ;
-    if(ierflg != 0){
-      cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : z = " << z << endl ;
-      return kFALSE;
-    }
-    gMinuit->mnparm(index++, " ",  Energy , 0.05*Energy, 0., 4.*Energy, ierflg) ;
-    if(ierflg != 0){
-      cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : Energy = " << Energy << endl ;      
-      return kFALSE;
+AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(AliPHOSGeometry *geom) :
+  AliPHOSTrackSegmentMaker(geom),
+  fDefaultInit(kFALSE),
+  fWrite(kFALSE),
+  fNTrackSegments(0),
+  fRcpv(0.f),
+  fRtpc(0.f),
+  fVtx(0.f), 
+  fLinkUpArray(0),
+  fEmcFirst(0),
+  fEmcLast(0),
+  fCpvFirst(0),
+  fCpvLast(0),
+  fModule(0),
+  fTrackSegments(NULL)
+{
+  // ctor
+  InitParameters() ; 
+  Init() ;
+  fESD = 0;
 }
-  }
 
-  Double_t p0=0.1; //"Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
-                  //depends on it. 
-  Double_t p1 = 1.;
-  
-  gMinuit->mnexcm("SET STR", 0, 0, ierflg) ;   //force TMinuit to reduce function calls  
-  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; //force TMinuit to use my gradient  
-  gMinuit->SetMaxIterations(5);
-  gMinuit->mnexcm("SET NOW", 0 , 0, ierflg) ; //No Warnings
-  gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg);   // minimize 
-  if(ierflg == 4){  //Minimum not found   
-      cout << "PHOS Unfolding>  Fit not converged, cluster abondoned "<< endl ;      
-    return kFALSE ;
-  }            
-  for(index = 0; index < NPar; index++){
-    Double_t err ;
-    Double_t val ;
-    gMinuit->GetParameter(index, val, err) ;    //Returns value and error of parameter index
-    FitParametres[index] = val ;
-   }
-  gMinuit->Delete() ;
-  return kTRUE;
+
+AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const AliPHOSTrackSegmentMakerv1 & tsm) :
+  AliPHOSTrackSegmentMaker(tsm),
+  fDefaultInit(kFALSE),
+  fWrite(kFALSE),
+  fNTrackSegments(0),
+  fRcpv(0.f),
+  fRtpc(0.f),
+  fVtx(0.f), 
+  fLinkUpArray(0),
+  fEmcFirst(0),
+  fEmcLast(0),
+  fCpvFirst(0),
+  fCpvLast(0),
+  fModule(0),
+  fTrackSegments(NULL)
+{
+  // cpy ctor: no implementation yet
+  // requested by the Coding Convention
+  Fatal("cpy ctor", "not implemented") ;
 }
+
+
 //____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::FillOneModule(DigitsList * Dl, RecPointsList * emcIn, TObjArray * emcOut, 
-                                       RecPointsList * ppsdIn, TObjArray * ppsdOutUp,
-                                       TObjArray * ppsdOutLow, Int_t &PHOSMod, Int_t & emcStopedAt, 
-                                       Int_t & ppsdStopedAt)
-{// Unfold clusters and fill xxxOut arrais with clusters from ome PHOS modeule
-  AliPHOSEmcRecPoint * emcRecPoint  ; 
-  AliPHOSPpsdRecPoint * ppsdRecPoint ;
-  Int_t index ;
-
-  Int_t NemcUnfolded = emcIn->GetEntries() ;
-  for(index = emcStopedAt; index < NemcUnfolded; index++){
-
-    emcRecPoint = (AliPHOSEmcRecPoint *) (*emcIn)[index] ;
-    
-    if(emcRecPoint->GetPHOSMod() != PHOSMod )  
-      break ;
-    
-    Int_t NMultipl = emcRecPoint->GetMultiplicity() ; 
-    int maxAt[NMultipl] ;
-    Float_t maxAtEnergy[NMultipl] ;
-    Int_t Nmax = emcRecPoint->GetNumberOfLocalMax(maxAt,maxAtEnergy) ;
-
-    if(Nmax<=1)     // if cluster is very flat, so that no prononsed maximum, then Nmax = 0 
-      emcOut->Add(emcRecPoint) ;
-    else {
-      UnfoldClusters(Dl, emcIn, emcRecPoint, Nmax, maxAt, maxAtEnergy, emcOut) ;
-      emcIn->Remove(emcRecPoint); 
-      emcIn->Compress() ;
-      NemcUnfolded-- ;
-      index-- ;
-    }
-  }
-  emcStopedAt = index ;
-
-  for(index = ppsdStopedAt; index < ppsdIn->GetEntries(); index++){
-    ppsdRecPoint = (AliPHOSPpsdRecPoint *) (*ppsdIn)[index] ;
-    if(ppsdRecPoint->GetPHOSMod() != PHOSMod )   break ;
-    if(ppsdRecPoint->GetUp() ) 
-      ppsdOutUp->Add(ppsdRecPoint) ;
-    else  
-      ppsdOutLow->Add(ppsdRecPoint) ;
+ AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
+{ 
+  // dtor
+  // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
+  if (!fDefaultInit)  
+    delete fLinkUpArray ;
+  if (fTrackSegments) {
+    fTrackSegments->Delete();
+    delete fTrackSegments;
   }
-  ppsdStopedAt = index ;
-  
-  PHOSMod ++ ; 
-  
-  emcOut->Sort() ;
-  ppsdOutUp->Sort() ;
-  ppsdOutLow->Sort() ; 
-  
 }
+
 //____________________________________________________________________________
-Float_t  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * EmcClu,AliPHOSPpsdRecPoint * PpsdClu, Bool_t &TooFar)
+void  AliPHOSTrackSegmentMakerv1::FillOneModule()
 {
-  Float_t R = fR0 ;
-  TVector3 vecEmc ;
-  TVector3 vecPpsd ;
+  // Finds first and last indexes between which 
+  // clusters from one PHOS module are
+
+  //First EMC clusters
+  Int_t totalEmc = fEMCRecPoints->GetEntriesFast() ;
+  for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&  
+       ((dynamic_cast<AliPHOSRecPoint *>(fEMCRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule ); 
+      fEmcLast ++)  ;
   
-  EmcClu->GetLocalPosition(vecEmc) ;
-  PpsdClu->GetLocalPosition(vecPpsd)  ; 
-  if(EmcClu->GetPHOSMod() == PpsdClu->GetPHOSMod()){ 
-    if(vecPpsd.X() >= vecEmc.X() - fDelta ){ 
-      if(vecPpsd.Z() >= vecEmc.Z() - fDelta ){
-       AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-       //Correct to difference in CPV and EMC position due to different distance to center.
-       //we assume, that particle moves from center
-       Float_t DCPV = geom->GetIPtoOuterCoverDistance();
-       Float_t DEMC = geom->GetIPtoCrystalSurface() ;
-       DEMC         = DEMC / DCPV ;
-        vecPpsd = DEMC * vecPpsd  - vecEmc ; 
-        R = vecPpsd.Mag() ;
-      } // if  zPpsd >= zEmc - fDelta
-      TooFar = kFALSE ;
-    } // if  xPpsd >= xEmc - fDelta
-    else 
-      TooFar = kTRUE ;
-  } 
-  else 
-    TooFar = kTRUE ;
-  
-  return R ;
+  //Now CPV clusters
+  Int_t totalCpv = fCPVRecPoints->GetEntriesFast() ;
+
+    for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) && 
+         ((dynamic_cast<AliPHOSRecPoint *>(fCPVRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule ); 
+       fCpvLast ++) ;
+      
 }
 
 //____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::MakeLinks(TObjArray * EmcRecPoints, TObjArray * PpsdRecPointsUp, 
-                                    TObjArray * PpsdRecPointsLow, TClonesArray * LinkLowArray, 
-                                    TClonesArray *LinkUpArray) 
-{ 
-  //Finds distanses (links) between all EMC and PPSD clusters, which are not further from each other than fR0 
+void  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,
+                                                         AliPHOSCpvRecPoint * cpvClu, 
+                                                         Int_t &trackindex, 
+                                                         Float_t &dx, Float_t &dz) const
+{
+  // Calculates the distance between the EMC RecPoint and the CPV RecPoint
+  // If no CPV, calculates the distance between the EMC RecPoint and the track
+  // prolongation to the PHOS module plane.
+  // Clusters are sorted in "rows" and "columns" of width 1 cm
 
-  TIter nextEmc(EmcRecPoints) ;
-  Int_t iEmcClu = 0 ; 
-  
-  AliPHOSPpsdRecPoint * PpsdLow ; 
-  AliPHOSPpsdRecPoint * PpsdUp ;
-  AliPHOSEmcRecPoint * EmcClu ;
-  
-  Int_t iLinkLow = 0 ;
-  Int_t iLinkUp  = 0 ;
-  
-  while( (EmcClu = (AliPHOSEmcRecPoint*)nextEmc() ) ) {
-    Bool_t TooFar ;
-    TIter nextPpsdLow(PpsdRecPointsLow ) ;
-    Int_t iPpsdLow = 0 ;
+//  Float_t delta = 1 ;  // Width of the rows in sorting of RecPoints (in cm)
+//                       // if you change this value, change it as well in xxxRecPoint::Compare()
+
+  if(!cpvClu) {
     
-    while( (PpsdLow = (AliPHOSPpsdRecPoint*)nextPpsdLow() ) ) { 
-      Float_t R = GetDistanceInPHOSPlane(EmcClu, PpsdLow, TooFar) ;
-      
-      if(TooFar) 
-       break ;  
-      if(R < fR0) 
-       new( (*LinkLowArray)[iLinkLow++]) AliPHOSLink(R, iEmcClu, iPpsdLow) ;
-      
-      iPpsdLow++ ;
+    trackindex = -1;
+    dx=999.;
+    dz=999.;
+    
+    if(!emcClu) {
+      return;
+    }
+
+    // *** Start the matching
+    Int_t nt=fESD->GetNumberOfTracks();
+    Int_t iPHOSMod = emcClu->GetPHOSMod()  ;
+    //Calculate actual distance to PHOS module
+    TVector3 globaPos ;
+    fGeom->Local2Global(iPHOSMod, 0.,0., globaPos) ;
+    const Double_t rPHOS = globaPos.Pt() ; //Distance to PHOS module
+    const Double_t kYmax = 72.+10. ; //Size of the module (with some reserve) in phi direction
+    const Double_t kZmax = 64.+10. ; //Size of the module (with some reserve) in z direction
+    const Double_t kAlpha= 20./180.*TMath::Pi() ; //TPC sector angular size
+    Double_t minDistance = 1.e6;
+
+    TMatrixF gmat;
+    TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
+    emcClu->GetGlobalPosition(gposRecPoint,gmat);
+    Double_t gposTrack[3] ; 
+
+    Double_t bz = AliTracker::GetBz() ; //B-Field for approximate matching
+    Double_t b[3]; 
+    for (Int_t i=0; i<nt; i++) {
+      AliESDtrack *esdTrack=fESD->GetTrack(i);
+
+//     // Skip the tracks having "wrong" status (has to be checked/tuned)
+//     ULong_t status = esdTrack->GetStatus();
+//     if ((status & AliESDtrack::kTRDout)   == 0) continue;
+//     if ((status & AliESDtrack::kTRDrefit) == 1) continue;
+
+      AliExternalTrackParam t(*esdTrack);
+      Int_t isec=Int_t(t.GetAlpha()/kAlpha);
+      Int_t imod=-isec-2; // PHOS module
+
+      Double_t y;                       // Some tracks do not reach the PHOS
+      if (!t.GetYAt(rPHOS,bz,y)) continue; //    because of the bending
+
+      Double_t z; t.GetZAt(rPHOS,bz,z);
+      if (TMath::Abs(z) > kZmax) continue; // Some tracks miss the PHOS in Z
+
+      Bool_t ok=kTRUE;
+      while (TMath::Abs(y) > kYmax) {   // Find the matching module
+        Double_t alp=t.GetAlpha();
+        if (y > kYmax) {
+          if (!t.Rotate(alp+kAlpha)) {ok=kFALSE; break;}
+          imod--;
+        } else if (y < -kYmax) {
+          if (!t.Rotate(alp-kAlpha)) {ok=kFALSE; break;}
+          imod++;
+        }
+        if (!t.GetYAt(rPHOS,bz,y)) {ok=kFALSE; break;}
+      }
+      if (!ok) continue; // Track rotation failed
+
+
+      if(imod!= iPHOSMod-1) 
+       continue; //not even approximate coincidence
+
+      //t.CorrectForMaterial(...); // Correct for the TOF material, if needed
+      t.GetBxByBz(b) ;
+      t.PropagateToBxByBz(rPHOS,b);        // Propagate to the matching module
+      t.GetXYZ(gposTrack) ;
       
+      Double_t ddx = gposTrack[0] - gposRecPoint.X(), ddy = gposTrack[1] - gposRecPoint.Y(), ddz = gposTrack[2] - gposRecPoint.Z();
+      Double_t d2 = ddx*ddx + ddy*ddy + ddz*ddz;
+      if(d2 < minDistance) {
+       dx = TMath::Sign(TMath::Sqrt(ddx*ddx + ddy*ddy),ddx) ;
+       dz = ddz ;
+       trackindex=i;
+       minDistance=d2 ;
+      }
     }
+    return ;
+  }
+
+
+  TVector3 emcGlobal;
+  fGeom->GetGlobalPHOS((AliPHOSRecPoint*)emcClu,emcGlobal);
     
-    TIter nextPpsdUp(PpsdRecPointsUp ) ;
-    Int_t iPpsdUp = 0 ;
+  // Radius from IP to current point
+  Double_t rEMC = TMath::Abs(emcGlobal.Pt());
     
-    while( (PpsdUp = (AliPHOSPpsdRecPoint*)nextPpsdUp() ) ) { 
-      Float_t R = GetDistanceInPHOSPlane(EmcClu, PpsdUp, TooFar) ;
-      
-      if(TooFar)
-       break ;  
-      if(R < fR0) 
-       new( (*LinkUpArray)[iLinkUp++]) AliPHOSLink(R, iEmcClu, iPpsdUp) ;
-      
-      iPpsdUp++ ;
-      
-    }
+  // Extrapolate the global track direction to EMC
+  // and find the closest track
     
-    iEmcClu++ ; 
+  Int_t nTracks = fESD->GetNumberOfTracks();
     
-  } // while nextEmC
-  
-  LinkLowArray->Sort() ; //first links with smallest distances
-  LinkUpArray->Sort() ;
-}
+  AliESDtrack *track;
+  Double_t xyz[] = {-1,-1,-1};
+  Double_t pxyz[3];
+  Double_t zEMC,xEMC;
+  Int_t module;
+  TVector3 vecP;
+  TVector3 locClu;
+
+  Float_t minDistance = 1.e6;
+  Float_t dr;
+
+  for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
+    track = fESD->GetTrack(iTrack);
+    if (!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz)) continue;
     
-//____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::MakePairs(TObjArray * EmcRecPoints, TObjArray * PpsdRecPointsUp, 
-                                   TObjArray * PpsdRecPointsLow, TClonesArray * LinkLowArray, 
-                                   TClonesArray * LinkUpArray, TrackSegmentsList * trsl) 
-{ // Finds the smallest links and makes pairs of PPSD and EMC clusters with smallest distance 
-  TIter nextLow(LinkLowArray) ;
-  TIter nextUp(LinkUpArray) ;
-  
-  AliPHOSLink * linkLow ;
-  AliPHOSLink * linkUp ;
+    AliDebug(1,Form("Event %d, iTrack: %d, (%.3f,%.3f,%.3f)",
+                   fESD->GetEventNumberInFile(),iTrack,xyz[0],xyz[1],xyz[2]));
+      
+    if (track->GetPxPyPzAt(rEMC,fESD->GetMagneticField(),pxyz)) { 
 
-  AliPHOSEmcRecPoint * emc ;
-  AliPHOSPpsdRecPoint * ppsdLow ;
-  AliPHOSPpsdRecPoint * ppsdUp ;
+      vecP.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
+      fGeom->ImpactOnEmc(xyz,vecP.Theta(),vecP.Phi(),module,zEMC,xEMC) ;
 
-  while ( (linkLow =  (AliPHOSLink *)nextLow() ) ){
-    emc = (AliPHOSEmcRecPoint *) EmcRecPoints->At(linkLow->GetEmc()) ;
-    ppsdLow = (AliPHOSPpsdRecPoint *) PpsdRecPointsLow->At(linkLow->GetPpsd()) ;
-    if((emc)&&(ppsdLow)){ // RecPoints not removed yet 
-        ppsdUp = 0 ;
-        
-        while ( (linkUp =  (AliPHOSLink *)nextUp() ) ){  
-          if(linkLow->GetEmc() == linkUp->GetEmc() ){
-            ppsdUp = (AliPHOSPpsdRecPoint *) PpsdRecPointsUp->At(linkUp->GetPpsd()) ;
-            break ;
-          }
+      if(!module) continue;
+      AliDebug(1,Form("\t\tTrack hit PHOS! Module: %d, (x,z)=(%.3f,%.3f)",module,xEMC,zEMC));
        
-        } // while nextUp
-        
-        nextUp.Reset();
-         AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
-        trsl->Add(subtr) ;  
-        EmcRecPoints->RemoveAt(linkLow->GetEmc()) ;      
-        PpsdRecPointsLow->RemoveAt(linkLow->GetPpsd()) ;
-        
-        if(ppsdUp)  
-          PpsdRecPointsUp->RemoveAt(linkUp->GetPpsd()) ;
-        
-    } // if NLocMax
+      if(emcClu->GetPHOSMod() != module) continue;
+       
+      // match track to EMC cluster
+      emcClu->GetLocalPosition(locClu);
+       
+      Float_t delta_x = xEMC - locClu.X();
+      Float_t delta_z = zEMC - locClu.Z();
+      dr = TMath::Sqrt(delta_x*delta_x + delta_z*delta_z);
+      AliDebug(1,Form("\tMatch iTrack=%d: (dx,dz)=(%.3f,%.3f)",iTrack,delta_x,delta_z));
+       
+      if(dr<minDistance) {
+       trackindex = iTrack;
+       minDistance = dr;
+       dx = delta_x;
+       dz = delta_z;
+      }
+    }
+      
+  }
+    
+  if(trackindex>=0)
+    AliDebug(1,Form("\t\tBest match for (xClu,zClu,eClu)=(%.3f,%.3f,%.3f): iTrack=%d, dR=%.3f",
+                   locClu.X(),locClu.Z(),emcClu->GetEnergy(),
+                   trackindex,TMath::Sqrt(dx*dx+dz*dz)));        
+  return;
+  
+  Float_t distance2Track = fRtpc ; 
+  
+  trackindex = -1 ; // closest track within fRCpv 
+  
+  TVector3 vecEmc ;   // Local position of EMC recpoint
+  TVector3 vecPloc ;     // Momentum direction at CPV plain
+  
+  //toofar = kTRUE ;
+  if(emcClu->GetPHOSMod() != cpvClu->GetPHOSMod()){
+    dx=999. ;
+    dz=999. ;
+    return ;
+  }
+  
+  emcClu->GetLocalPosition(vecEmc) ;
+  
+  Double_t xCPV,zCPV ; //EMC-projected coordinates of CPV cluster 
+  TVector3 cpvGlobal; // Global position of the CPV recpoint
+  fGeom->GetGlobalPHOS((AliPHOSRecPoint*)cpvClu,cpvGlobal);
+  Double_t vtxCPV[3]={cpvGlobal.X(),cpvGlobal.Y(),cpvGlobal.Z()} ;
+  Int_t dummyMod ;
+
+  if (fESD == 0x0) {
+    //if no track information available, assume straight line from IP to emcal
+    fGeom->ImpactOnEmc(vtxCPV,cpvGlobal.Theta(),cpvGlobal.Phi(),dummyMod,zCPV,xCPV) ;
+    dx=xCPV - vecEmc.X() ;
+    dz=zCPV - vecEmc.Z() ;
+    return ;
   } 
-   
-  TIter nextEmc(EmcRecPoints) ;          
-  nextEmc.Reset() ;
+  
+  //if there is ESD try to correct distance using TPC information on particle direct in CPV 
+  if (fESD != 0x0) {
+
+    Double_t rCPV = cpvGlobal.Pt() ;// Radius from IP to current point 
+
+    // Extrapolate the global track direction if any to CPV and find the closest track
+    Int_t iClosestTrack = -1;
+    TVector3 inPHOS ; 
+
+    for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
+      track = fESD->GetTrack(iTrack);
+      if (!track->GetXYZAt(rCPV, fESD->GetMagneticField(), xyz))
+       continue; //track coord on the cylinder of PHOS radius
+      if ((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0)
+       continue;
+      //Check if this track hits PHOS
+      inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
+      distance2Track = inPHOS.Angle(cpvGlobal) ;
+      // Find the closest track to the CPV recpoint
+      if (distance2Track < minDistance) {
+       minDistance = distance2Track;
+       iClosestTrack = iTrack;
+      }
+    }
 
-  while( (emc = (AliPHOSEmcRecPoint*)nextEmc()) ){ //to create pairs if no PpsdLow
-    ppsdLow = NULL ; 
-    ppsdUp = NULL ;
-    
-    while ( (linkUp =  (AliPHOSLink *)nextUp() ) ){
-      
-      if(EmcRecPoints->IndexOf(emc) == linkUp->GetEmc() ){
-       ppsdUp = (AliPHOSPpsdRecPoint *) PpsdRecPointsUp->At(linkUp->GetPpsd()) ;
-       break ;
+    if (iClosestTrack != -1) {
+      track = fESD->GetTrack(iClosestTrack);
+      if (track->GetPxPyPzAt(rCPV, fESD->GetMagneticField(), pxyz)) { // track momentum ibid.
+        vecP.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
+       fGeom->ImpactOnEmc(vtxCPV,vecP.Theta(),vecP.Phi(),dummyMod,zCPV,xCPV) ;
       }
-      
     }
-    AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
-    trsl->Add(subtr) ;   
-    if(ppsdUp)  
-      PpsdRecPointsUp->RemoveAt(linkUp->GetPpsd()) ;
     
+    if(minDistance < fRtpc ){
+      trackindex = iClosestTrack ; 
+    }
+  }
+  if(trackindex!=-1){
+    // If the closest global track is found, calculate EMC-CPV distance from it
+    dx=xCPV - vecEmc.X() ;
+    dz=zCPV - vecEmc.Z() ;
   }
-     
+  else{
+    // If no global track was found, just take the nearest CPV point
+    fGeom->ImpactOnEmc(vtxCPV,cpvGlobal.Theta(),cpvGlobal.Phi(),dummyMod,zCPV,xCPV) ;
+    dx=xCPV - vecEmc.X() ;
+    dz=zCPV - vecEmc.Z() ;
+  }
+  return ;
 }
-
 //____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::MakeTrackSegments(DigitsList * DL, RecPointsList * emcl, 
-                                       RecPointsList * ppsdl, TrackSegmentsList * trsl)
+void  AliPHOSTrackSegmentMakerv1::Init()
 {
-  //main function, does the job
-
-  Int_t PHOSMod    = 1 ;
-  Int_t emcStopedAt = 0 ; 
-  Int_t ppsdStopedAt = 0 ; 
-  
-  TObjArray * EmcRecPoints = new TObjArray(100) ;     //these arrays keeps pointers 
-  TObjArray * PpsdRecPointsUp = new TObjArray(100) ;  //on RecPoints, which are 
-  TObjArray * PpsdRecPointsLow = new TObjArray(100) ; //kept in TClonesArray's emcl and ppsdl
-  
-  
-  TClonesArray * LinkLowArray = new TClonesArray("AliPHOSLink", 100);
-  TClonesArray * LinkUpArray  = new TClonesArray("AliPHOSLink", 100); 
+  // Make all memory allocations that are not possible in default constructor
   
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-  
-  while(PHOSMod <= geom->GetNModules() ){
-    
-    FillOneModule(DL, emcl, EmcRecPoints, ppsdl, PpsdRecPointsUp, PpsdRecPointsLow, PHOSMod , emcStopedAt, ppsdStopedAt) ;
-   
-    MakeLinks(EmcRecPoints,PpsdRecPointsUp, PpsdRecPointsLow, LinkLowArray, LinkUpArray) ; 
+  fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000); 
+  fTrackSegments = new TClonesArray("AliPHOSTrackSegment",100);
+  fTrackSegments->SetName("TRACKS");
+}
 
-    MakePairs(EmcRecPoints,PpsdRecPointsUp, PpsdRecPointsLow, LinkLowArray, LinkUpArray, trsl) ;
-    EmcRecPoints->Clear() ;
-    PpsdRecPointsUp->Clear() ;
-    PpsdRecPointsLow->Clear() ;
-    LinkUpArray->Delete();
-    LinkLowArray->Delete();
-  }
+//____________________________________________________________________________
+void  AliPHOSTrackSegmentMakerv1::InitParameters()
+{
+  //Initializes parameters
+  fRcpv      = 10. ;
+  fRtpc      = 4. ;
+  fEmcFirst  = 0 ;    
+  fEmcLast   = 0 ;   
+  fCpvFirst  = 0 ;   
+  fCpvLast   = 0 ;   
+  fLinkUpArray = 0 ;
+  fWrite                   = kTRUE ;
 }
 
+
 //____________________________________________________________________________
-Double_t  AliPHOSTrackSegmentMakerv1::ShowerShape(Double_t r){ 
-// If you change this function, change also gradiend evaluation  in ChiSquare()
-  Double_t r4 = r*r*r*r ;
-  Double_t r295 = TMath::Power(r, 2.95) ;
-  Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
-  return shape ;
+void  AliPHOSTrackSegmentMakerv1::MakeLinks()const
+{ 
+  // Finds distances (links) between all EMC and CPV clusters, 
+  // which are not further apart from each other than fRcpv 
+  // and sort them in accordance with this distance
+  
+  fLinkUpArray->Clear() ;    
+
+  AliPHOSCpvRecPoint * cpv ;
+  AliPHOSEmcRecPoint * emcclu ;
+
+  Int_t iLinkUp  = 0 ;
+  
+  Int_t iEmcRP;
+  for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
+    emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP)) ;
+
+    //Bool_t toofar ;        
+    Int_t iCpv = 0 ;    
+    for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) { 
+      
+      cpv = dynamic_cast<AliPHOSCpvRecPoint *>(fCPVRecPoints->At(iCpv)) ;
+      Int_t track = -1 ; 
+      Float_t dx,dz ;
+      GetDistanceInPHOSPlane(emcclu, cpv, track,dx,dz) ;     
+      if(TMath::Sqrt(dx*dx+dz*dz) < fRcpv ){ 
+        new ((*fLinkUpArray)[iLinkUp++])  AliPHOSLink(dx, dz, iEmcRP, iCpv, track) ;
+      }      
+    }
+  } 
+  
+  fLinkUpArray->Sort() ;  //first links with smallest distances
 }
 
 //____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::UnfoldClusters(DigitsList * DL, RecPointsList * emcIn,  AliPHOSEmcRecPoint * iniEmc, 
-                                        Int_t Nmax, int * maxAt, Float_t * maxAtEnergy, TObjArray * emcList)
+void  AliPHOSTrackSegmentMakerv1::MakePairs()
 { 
-  //fits cluster with Nmax overlapping showers 
+  // Using the previously made list of "links", we found the smallest link - i.e. 
+  // link with the least distance between EMC and CPV and pointing to still 
+  // unassigned RecParticles. We assign these RecPoints to TrackSegment and 
+  // remove them from the list of "unassigned". 
+
+  //Make arrays to mark clusters already chosen
+  Int_t * emcExist = 0;
+  if(fEmcLast > fEmcFirst)
+    emcExist = new Int_t[fEmcLast-fEmcFirst] ;
   
-  Int_t NPar = 3 * Nmax ;
-  Float_t FitParameters[NPar] ;
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
+  Int_t index;
+  for(index = 0; index <fEmcLast-fEmcFirst; index ++)
+    emcExist[index] = 1 ;
   
-  if( !FindFit(iniEmc, maxAt, maxAtEnergy, NPar, FitParameters) )  //Fit failed, return and remove cluster
-    return ;
+  Bool_t * cpvExist = 0;
+  if(fCpvLast > fCpvFirst)
+    cpvExist = new Bool_t[fCpvLast-fCpvFirst] ;
+  for(index = 0; index <fCpvLast-fCpvFirst; index ++)
+    cpvExist[index] = kTRUE ;
   
-  Float_t xDigit ;
-  Float_t zDigit ;
-  Int_t RelId[4] ;
-
-  Int_t Ndigits = iniEmc->GetMultiplicity() ;  
-  Float_t xpar  ;
-  Float_t zpar  ;
-  Float_t Epar  ;
-  Float_t Distance ;
-  Float_t Ratio ;
-  Float_t Efit[Ndigits] ;
-  Int_t iparam ;
-  Int_t iDigit ;
   
-  AliPHOSDigit * digit ;
-  AliPHOSEmcRecPoint * emcRP ;  
-  int * emcDigits = iniEmc->GetDigitsList() ;
-  Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
-
-  Int_t iRecPoint = emcIn->GetEntries() ;
-
-  for(iDigit = 0 ; iDigit < Ndigits ; iDigit ++){
-    digit = (AliPHOSDigit *) emcDigits[iDigit];
-    geom->AbsToRelNumbering(digit->GetId(), RelId) ;
-    geom->RelPosInModule(RelId, xDigit, zDigit) ;
-    Efit[iDigit] = 0;
-    iparam = 0 ;
-    
-    while(iparam < NPar ){
-      xpar = FitParameters[iparam] ;
-      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) ;
+  // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance 
+  TIter nextUp(fLinkUpArray) ;
+  
+  AliPHOSLink * linkUp ;
+  
+  AliPHOSCpvRecPoint * nullpointer = 0 ;
+  
+  while ( (linkUp =  static_cast<AliPHOSLink *>(nextUp()) ) ){  
+
+    if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
+
+      if(cpvExist[linkUp->GetCpv()-fCpvFirst]){ //CPV still exist
+         Float_t dx,dz ;
+         linkUp->GetXZ(dx,dz) ;
+        new ((* fTrackSegments)[fNTrackSegments]) 
+          AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(linkUp->GetEmc())) , 
+                              dynamic_cast<AliPHOSCpvRecPoint *>(fCPVRecPoints->At(linkUp->GetCpv())) , 
+                              linkUp->GetTrack(),dx,dz) ;
+        
+       (dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
+       fNTrackSegments++ ;
+       emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc  that Cpv was found 
+       //mark CPV recpoint as already used 
+       cpvExist[linkUp->GetCpv()-fCpvFirst] = kFALSE ;
+      } //if CpvUp still exist
+    } 
+  }        
+
+  //look through emc recPoints left without CPV
+  if(emcExist){ //if there is emc rec point
+    Int_t iEmcRP ;
+    for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst  ; iEmcRP++ ){
+      if(emcExist[iEmcRP] > 0 ){
+        Int_t track = -1 ;
+        Float_t dx,dz ;
+        AliPHOSEmcRecPoint *emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP+fEmcFirst));
+        GetDistanceInPHOSPlane(emcclu, 0, track,dx,dz);
+        if(track<0)
+         new ((*fTrackSegments)[fNTrackSegments]) AliPHOSTrackSegment(emcclu,nullpointer) ;
+        else
+          new ((*fTrackSegments)[fNTrackSegments]) AliPHOSTrackSegment(emcclu,0,track,dx,dz);
+       (dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
+       fNTrackSegments++;    
+      } 
     }
+  }
+  delete [] emcExist ; 
+  delete [] cpvExist ; 
+}
 
+//____________________________________________________________________________
+void AliPHOSTrackSegmentMakerv1::Clusters2TrackSegments(Option_t *option)
+{
+  // Steering method to perform track segment construction for the current event
+  // Returns an array with the found track-segments.
+  
+  if(strstr(option,"tim"))
+    gBenchmark->Start("PHOSTSMaker");
+  if(strstr(option,"print")) {
+    Print() ; 
+    return ; 
   }
+  
+  //Make some initializations 
+  fNTrackSegments = 0 ;
+  fEmcFirst = 0 ;    
+  fEmcLast  = 0 ;   
+  fCpvFirst = 0 ;   
+  fCpvLast  = 0 ;   
 
-  iparam = 0 ;
-  Float_t eDigit ;
-
-  while(iparam < NPar ){
-    xpar = FitParameters[iparam] ;
-    zpar = FitParameters[iparam+1] ;
-    Epar = FitParameters[iparam+2] ;
-    iparam += 3 ;
-    new ((*emcIn)[iRecPoint]) AliPHOSEmcRecPoint( iniEmc->GetLogWeightCut(), iniEmc->GetLocMaxCut() ) ;
-    emcRP = (AliPHOSEmcRecPoint *) emcIn->At(iRecPoint++);
-
-    for(iDigit = 0 ; iDigit < Ndigits ; iDigit ++){
-      digit = (AliPHOSDigit *) emcDigits[iDigit];
-      geom->AbsToRelNumbering(digit->GetId(), RelId) ;
-      geom->RelPosInModule(RelId, xDigit, zDigit) ;
-      Distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
-      Distance =  TMath::Sqrt(Distance) ;
-      Ratio = Epar * ShowerShape(Distance) / Efit[iDigit] ; 
-      eDigit = emcEnergies[iDigit] * Ratio ;
-      emcRP->AddDigit( *digit,eDigit ) ;
-    }
+  fTrackSegments->Clear();
 
-    emcList->Add(emcRP) ;
+  //   if(!ReadRecPoints(ievent))   continue; //reads RecPoints for event ievent
 
+  for(fModule = 1; fModule <= fGeom->GetNModules() ; fModule++ ) {
+    FillOneModule() ; 
+    MakeLinks() ;
+    MakePairs() ;
+  }
+    
+  if(strstr(option,"deb"))
+    PrintTrackSegments(option);
+
+  if(strstr(option,"tim")){
+    gBenchmark->Stop("PHOSTSMaker");
+    Info("Exec", "took %f seconds for making TS", 
+        gBenchmark->GetCpuTime("PHOSTSMaker")); 
   }
 }
-//______________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::UnfoldingChiSquare(Int_t &NPar, Double_t *Grad, Double_t & fret, Double_t *x, Int_t iflag)
+
+//____________________________________________________________________________
+void AliPHOSTrackSegmentMakerv1::Print(const Option_t *)const
 {
-// NUmber of paramters, Gradient , Chi squared, parameters, what to do
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-  AliPHOSTrackSegmentMakerv1 TrS ;
-
-  AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; //EmcRecPoint to fit
-  int * emcDigits = emcRP->GetDigitsList() ;
-  Float_t * emcEnergies = emcRP->GetEnergiesList() ;
-  fret = 0. ;     
-  Int_t iparam ;
-
-  if(iflag==2) 
-    for(iparam = 0 ; iparam < NPar ; iparam ++)    
-      Grad[iparam] = 0 ; //Will evaluate gradient
-  Double_t Efit ;  
-
-  AliPHOSDigit * digit ;
-  Int_t iDigit = 0 ;
-  while ( (digit = (AliPHOSDigit *)emcDigits[iDigit] )){
-    Int_t RelId[4] ;
-    Float_t xDigit ;
-    Float_t zDigit ;
-    geom->AbsToRelNumbering(digit->GetId(), RelId) ;
-    geom->RelPosInModule(RelId, xDigit, zDigit) ;
-    
-    if(iflag == 2){  //calculate gradient
-      Int_t iParam = 0 ;
-      Efit = 0 ;
-      while(iParam < NPar ){
-       Double_t Distance = TMath::Sqrt( (xDigit - x[iParam]) * (xDigit - x[iParam]) + 
-                                       (zDigit - x[++iParam]) * (zDigit - x[iParam]) ) ;
-       Efit += x[++iParam] * TrS.ShowerShape(Distance) ;
-       iParam++ ;
-      }
-      Double_t sum = 2. * (Efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; //Here we assume, that sigma = sqrt(E) 
-      iParam = 0 ;
-      while(iParam < NPar ){
-       Double_t xpar = x[iParam] ;
-       Double_t zpar = x[iParam+1] ;
-       Double_t Epar = x[iParam+2] ;
-       Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
-       Double_t shape = sum * TrS.ShowerShape(dr) ;
-       Double_t r4 = dr*dr*dr*dr ;
-       Double_t r295 = TMath::Power(dr,2.95) ;
-       Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
-                                  0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
-       
-       Grad[iParam++] += Epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x    
-       Grad[iParam++] += Epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z         
-       Grad[iParam++] += shape ;                           // Derivative over energy           
-      }
-    }
-    Efit = 0;
-    iparam = 0 ;
-    while(iparam < NPar ){
-      Double_t xpar = x[iparam] ;
-      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 * TrS.ShowerShape(Distance) ;
-    }
-    fret += (Efit-emcEnergies[iDigit])*(Efit-emcEnergies[iDigit])/emcEnergies[iDigit] ; 
-    //Here we assume, that sigma = sqrt(E)
-    iDigit++ ;
+  //  Print TrackSegmentMaker parameters
+
+  TString message("") ;
+  if( strcmp(GetName(), "") != 0 ) {
+    message = "\n======== AliPHOSTrackSegmentMakerv1 ========\n" ; 
+    message += "Making Track segments\n" ;
+    message += "with parameters:\n" ; 
+    message += "     Maximal EMC - CPV distance (cm) %f\n" ;
+    message += "============================================\n" ;
+    Info("Print", message.Data(),fRcpv) ;
   }
+  else
+    Info("Print", "AliPHOSTrackSegmentMakerv1 not initialized ") ;
 }
+
+//____________________________________________________________________________
+void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
+{
+  // option deb - prints # of found TrackSegments
+  // option deb all - prints as well indexed of found RecParticles assigned to the TS
+
+  Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ; 
+  printf("        Found %d TrackSegments\n", fTrackSegments->GetEntriesFast() ); 
+  
+  if(strstr(option,"all")) {  // printing found TS
+    printf("TrackSegment #  EMC RP#  CPV RP#\n") ; 
+    Int_t index;
+    for (index = 0 ; index <fTrackSegments->GetEntriesFast() ; index++) {
+      AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )fTrackSegments->At(index) ; 
+      printf("   %d           %d        %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetCpvIndex() ) ; 
+    }  
+  }
+}
+