]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSTrackSegmentMakerv1.cxx
Inserting TMath.h where required by the new version of ROOT
[u/mrichter/AliRoot.git] / PHOS / AliPHOSTrackSegmentMakerv1.cxx
index fc531e801bc70b52dbcb398b0cad8de6191260cd..203abda60f7c51fa5731b9b7d629c57afbcd8e3d 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.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.
+ *
+ */
 
 //_________________________________________________________________________
-// Algorithm class to construct track segments connection RecPoints in 
-// EMCA and Ppsd. Unfolds also the clusters in EMCA. 
-//*-- Author : D. Peressounko  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 "TObjectTable.h"
+#include "TTree.h"
+#include "TBenchmark.h"
 
 // --- Standard library ---
-
-#include <iostream>
-#include <cassert>
-
+#include "Riostream.h"
 // --- AliRoot header files ---
-
+#include "AliPHOSGeometry.h"
 #include "AliPHOSTrackSegmentMakerv1.h"
 #include "AliPHOSTrackSegment.h"
 #include "AliPHOSLink.h"
-#include "AliPHOSv0.h"
-#include "AliRun.h"
-
-extern void UnfoldingChiSquare(Int_t &NPar, Double_t *Grad, Double_t & fret, Double_t *x, Int_t iflag) ; 
+#include "AliPHOSGetter.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
 
 ClassImp( AliPHOSTrackSegmentMakerv1) 
 
 
 //____________________________________________________________________________
- AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() 
+AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() :
+  AliPHOSTrackSegmentMaker(),
+  fDefaultInit(kTRUE),
+  fWrite(kFALSE),
+  fNTrackSegments(0),
+  fRcpv(0.f),
+  fRtpc(0.f),
+  fLinkUpArray(0),
+  fEmcFirst(0),
+  fEmcLast(0),
+  fCpvFirst(0),
+  fCpvLast(0),
+  fModule(0),
+  fTrackSegmentsInRun(0)
+  
+{
+  // default ctor (to be used mainly by Streamer)
+  InitParameters() ; 
+}
+
+//____________________________________________________________________________
+AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const TString & alirunFileName, const TString & eventFolderName) :
+  AliPHOSTrackSegmentMaker(alirunFileName, eventFolderName),
+  fDefaultInit(kFALSE),
+  fWrite(kFALSE),
+  fNTrackSegments(0),
+  fRcpv(0.f),
+  fRtpc(0.f),
+  fLinkUpArray(0),
+  fEmcFirst(0),
+  fEmcLast(0),
+  fCpvFirst(0),
+  fCpvLast(0),
+  fModule(0),
+  fTrackSegmentsInRun(0)
 {
   // ctor
-  fR0 = 4. ;   
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-  //clusters are sorted in "rows" and "columns" of width geom->GetCrystalSize(0),
-  fDelta = fR0 + geom->GetCrystalSize(0) ;
-  fMinuit = new TMinuit(100) ; 
+  InitParameters() ; 
+  Init() ;
+  fESD = 0;
 }
 
+
+AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const AliPHOSTrackSegmentMakerv1 & tsm) :
+  AliPHOSTrackSegmentMaker(tsm),
+  fDefaultInit(kFALSE),
+  fWrite(kFALSE),
+  fNTrackSegments(0),
+  fRcpv(0.f),
+  fRtpc(0.f),
+  fLinkUpArray(0),
+  fEmcFirst(0),
+  fEmcLast(0),
+  fCpvFirst(0),
+  fCpvLast(0),
+  fModule(0),
+  fTrackSegmentsInRun(0)
+{
+  // cpy ctor: no implementation yet
+  // requested by the Coding Convention
+  Fatal("cpy ctor", "not implemented") ;
+}
+
+
 //____________________________________________________________________________
  AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
 { 
   // dtor
-  delete fMinuit ; 
+  // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
+  if (!fDefaultInit)  
+    delete fLinkUpArray ;
 }
-//____________________________________________________________________________
-Bool_t  AliPHOSTrackSegmentMakerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy,
-                                   Int_t NPar, Float_t * FitParameters)
-{ 
-  // gObjectTable->Print() ; 
-  // Calls TMinuit for fitting cluster with several maxima 
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-  assert( NPar < 100 ) ; 
 
-  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 ;
+//____________________________________________________________________________
+const TString AliPHOSTrackSegmentMakerv1::BranchName() const 
+{  
+  return GetName() ;
+}
 
-  Int_t iDigit ;
+//____________________________________________________________________________
+void  AliPHOSTrackSegmentMakerv1::FillOneModule()
+{
+  // Finds first and last indexes between which 
+  // clusters from one PHOS module are
 
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  
+  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
+  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
+  //First EMC clusters
+  Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
+  for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&  
+       ((dynamic_cast<AliPHOSRecPoint *>(emcRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule ); 
+      fEmcLast ++)  ;
+  
+  //Now CPV clusters
+  Int_t totalCpv = cpvRecPoints->GetEntriesFast() ;
 
-  for(iDigit = 0 ; iDigit < NDigits ; iDigit++){
-    digit = (AliPHOSDigit *) maxAt[iDigit]; 
+    for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) && 
+         ((dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule ); 
+       fCpvLast ++) ;
+      
+}
 
-    Int_t RelId[4] ;
-    Float_t x ;
-    Float_t z ;
-    geom->AbsToRelNumbering(digit->GetId(), RelId) ;
-    geom->RelPosInModule(RelId, x, z) ;
+//____________________________________________________________________________
+Float_t  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,AliPHOSCpvRecPoint * cpvClu, Int_t &trackindex) const
+{
+  // Calculates the distance between the EMC RecPoint and the CPV RecPoint
+  // Clusters are sorted in "rows" and "columns" of width 1 cm
 
-    Float_t Energy = maxAtEnergy[iDigit] ;
+  //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()
+  Float_t distance2Cpv   = fRcpv ;
+  Float_t distance2Track = fRtpc ; 
 
-    gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
-    index++ ;   
-    if(ierflg != 0){ 
-      cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : x = " << x << endl ;
-      return kFALSE;
-    }
-    gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
-    index++ ;   
-    if(ierflg != 0){
-      cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : z = " << z << endl ;
-      return kFALSE;
-    }
-    gMinuit->mnparm(index, "Energy",  Energy , 0.05*Energy, 0., 4.*Energy, ierflg) ;
-    index++ ;   
-    if(ierflg != 0){
-      cout << "PHOS Unfolding>  Unable to set initial value for fit procedure : Energy = " << Energy << endl ;      
-      return kFALSE;
-}
-  }
+  trackindex = -1 ; // closest track within fRCpv 
 
-  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.0 ;
-  Double_t p2 = 0.0 ;
-
-  gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TgMinuit to reduce function calls  
-  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient  
-  gMinuit->SetMaxIterations(5);
-  gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
-  gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize 
-  if(ierflg == 4){  // Minimum not found   
-    cout << "PHOS Unfolding>  Fit not converged, cluster abandoned "<< 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
-    FitParameters[index] = val ;
-   }
-  return kTRUE;
+  TVector3 vecEmc ;   // Local position of EMC recpoint
+  TVector3 vecCpv ;   // Local position of CPV recpoint propagated to EMC
+  TVector3 vecDist ;  // Distance between local positions of two points
+  
+  emcClu->GetLocalPosition(vecEmc) ;
+  cpvClu->GetLocalPosition(vecCpv) ;
 
-}
+  //toofar = kTRUE ;
+  if(emcClu->GetPHOSMod() == cpvClu->GetPHOSMod()){ 
 
-//____________________________________________________________________________
-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 arrays with clusters from one PHOS module
-  AliPHOSEmcRecPoint *  emcRecPoint  ; 
-  AliPHOSPpsdRecPoint * ppsdRecPoint ;
-  Int_t index ;
-  cout << "Fill 1" << endl ;
-  Int_t NemcUnfolded = emcIn->GetEntries() ;
-  for(index = emcStopedAt; index < NemcUnfolded; index++){
-    emcRecPoint = (AliPHOSEmcRecPoint *) (*emcIn)[index] ;
-    cout << "Fill 2" << endl ;
-   
-    if(emcRecPoint->GetPHOSMod() != PHOSMod )  
-       break ;
-    
+    // Find EMC-CPV distance
+    distance2Cpv = (vecCpv - vecEmc).Mag() ;
     
-    Int_t NMultipl = emcRecPoint->GetMultiplicity() ; 
-    int maxAt[NMultipl] ;
-    Float_t maxAtEnergy[NMultipl] ;
-    Int_t Nmax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy) ;
+    if (fESD != 0x0) {
+      AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+      const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
+
+      Double_t rPHOS = geom->GetIPtoCrystalSurface();
+
+      //PH Acceptance boundaries for each PHOS module
+      Int_t nModules = geom->GetNModules();
+      Double_t * thmin = new Double_t[nModules];// theta min
+      Double_t * thmax = new Double_t[nModules];// theta max
+      Double_t * phmin = new Double_t[nModules];// phi min
+      Double_t * phmax = new Double_t[nModules];// phi max
+      
+      for (Int_t imod=0; imod<nModules; imod++) {
+       geom->EmcModuleCoverage(imod,
+                               thmin[imod],thmax[imod],
+                               phmin[imod],phmax[imod]);
+      }
 
-    
+      // Extrapolate the global track direction if any to CPV and find the closest track
+      Int_t nTracks = fESD->GetNumberOfTracks();
+      Int_t iClosestTrack = -1;
+      Double_t minDistance = 1e6;
+      Double_t pxyz[3], xyz[3];
+      TVector3 inPHOS; //PH Used to calculate theta and phi
+
+      //PH Loop on tracks
+      AliESDtrack *track;
+      for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
+       track = fESD->GetTrack(iTrack);
+       if (!track->GetXYZAt(rPHOS, 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;
+       //PH Here one has to cut out the tracks which are not inside the PHOS
+       //PH acceptance
+       inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
+       Double_t inPhi = inPHOS.Phi();
+       Double_t inTheta = inPHOS.Theta();
+
+       Bool_t skip = kTRUE;
+       for (Int_t imod=0; imod<nModules; imod++) {
+         //PH Loop on modules to check if the track enters in the acceptance 
+         if (thmin[imod] < inTheta && thmax[imod] > inTheta && 
+             phmin[imod] < inPhi   && phmax[imod] > inPhi) {
+           skip = kFALSE;
+           break;
+         }
+       }
+       if (skip) continue; //PH Skip, if not in the PHOS acceptance
+
+       if (!track->GetPxPyPzAt(rPHOS, fESD->GetMagneticField(), pxyz))
+           continue; // track momentum ibid.
+       PropagateToPlane(vecDist,xyz,pxyz,"CPV",cpvClu->GetPHOSMod());
+       //      Info("GetDistanceInPHOSPlane","Track %d propagation to CPV = (%f,%f,%f)",
+       //     iTrack,vecDist.X(),vecDist.Y(),vecDist.Z());
+       vecDist -= vecCpv;
+       distance2Track = TMath::Sqrt(vecDist.X()*vecDist.X() + vecDist.Z()*vecDist.Z());
+       // Find the closest track to the EMC recpoint
+       if (distance2Track < minDistance) {
+         minDistance = distance2Track;
+         iClosestTrack = iTrack;
+       }
+      }
 
-    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); 
-      cout << "Fill 3" << endl ;
-      emcIn->Compress() ;
-      NemcUnfolded-- ;
-      index-- ;
+      delete [] thmin;
+      delete [] thmax;
+      delete [] phmin;
+      delete [] phmax;
+
+      if (iClosestTrack != -1) {
+       track = fESD->GetTrack(iClosestTrack);
+       if (track->GetPxPyPzAt(rPHOS, fESD->GetMagneticField(), pxyz)) { // track momentum ibid.
+       TVector3 vecCpvGlobal; // Global position of the CPV recpoint
+       geom->GetGlobal((AliRecPoint*)cpvClu,vecCpvGlobal);
+       for (Int_t ixyz=0; ixyz<3; ixyz++)
+         xyz[ixyz] = vecCpvGlobal[ixyz];
+       PropagateToPlane(vecDist,xyz,pxyz,"EMC",cpvClu->GetPHOSMod());
+//     Info("GetDistanceInPHOSPlane","Track %d propagation to EMC = (%f,%f,%f)",
+//          iClosestTrack,vecDist.X(),vecDist.Y(),vecDist.Z());
+       vecDist -= vecEmc;
+       distance2Track = TMath::Sqrt(vecDist.X()*vecDist.X() + vecDist.Z()*vecDist.Z());
+       }
+      }
+//     } else {
+//       // If no ESD exists, than simply find EMC-CPV distance
+//       distance = (vecCpv - vecEmc).Mag() ;
+    
+      //if(distance2Track < fRcpv + 2*delta )
+      if(distance2Track < fRtpc )
+       trackindex = iClosestTrack ; 
+      //      toofar = kFALSE ;
     }
+    //     Info("GetDistanceInPHOSPlane","cpv-emc distance is %f cm",
+    //          distance);
   }
-  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) ;
-  }
-  ppsdStopedAt = index ;
-  
-  PHOSMod ++ ; 
-  
-  emcOut->Sort() ;
-  ppsdOutUp->Sort() ;
-  ppsdOutLow->Sort() ; 
   
+  return distance2Cpv ;
 }
+
 //____________________________________________________________________________
-Float_t  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * EmcClu,AliPHOSPpsdRecPoint * PpsdClu, Bool_t &TooFar)
+void AliPHOSTrackSegmentMakerv1::PropagateToPlane(TVector3& globalIntersection,
+                                                 Double_t *x,
+                                                 Double_t *p,
+                                                 const char *det,
+                                                 Int_t moduleNumber) const
 {
-  Float_t R = fR0 ;
-  TVector3 vecEmc ;
-  TVector3 vecPpsd ;
+  // Propagate a straight-line track from the origin point x
+  // along the direction p to the CPV or EMC module moduleNumber
+  // Returns a local position of such a propagation
+
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
+  TVector3 moduleCenter;
+  geom->GetModuleCenter(moduleCenter,det,moduleNumber);
+  TVector3 vertex; vertex.SetXYZ(x[0],x[1],x[2]);
+  TVector3 direction; direction.SetXYZ(p[0],p[1],p[2]);
+
+//   Info("PropagateToCPV","Center of the %s module %d is (%f,%f,%f)",
+//        det,moduleNumber,moduleCenter[0],moduleCenter[1],moduleCenter[2]);
+
+  Double_t time = (moduleCenter.Mag2() - vertex.Dot(moduleCenter)) /
+    (direction.Dot(moduleCenter));
+  vertex += direction*time;
+  geom->Global2Local(globalIntersection,vertex,moduleNumber);
+}
+
+//____________________________________________________________________________
+void  AliPHOSTrackSegmentMakerv1::Init()
+{
+  // Make all memory allocations that are not possible in default constructor
   
-  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 ;
+  AliPHOSGetter* gime = AliPHOSGetter::Instance();
+  if(!gime)
+    gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
   
-  return R ;
+  fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000); 
+  if ( !gime->TrackSegmentMaker() ) {
+    gime->PostTrackSegmentMaker(this);
+  }
 }
 
 //____________________________________________________________________________
-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::InitParameters()
+{
+  //Initializes parameters
+  fRcpv      = 10. ;
+  fRtpc      = 4. ;
+  fEmcFirst  = 0 ;    
+  fEmcLast   = 0 ;   
+  fCpvFirst  = 0 ;   
+  fCpvLast   = 0 ;   
+  fLinkUpArray = 0 ;
+  fWrite                   = kTRUE ;
+  fTrackSegmentsInRun       = 0 ; 
+  SetEventRange(0,-1) ;
+}
 
-  TIter nextEmc(EmcRecPoints) ;
-  Int_t iEmcClu = 0 ; 
-  
-  AliPHOSPpsdRecPoint * PpsdLow ; 
-  AliPHOSPpsdRecPoint * PpsdUp ;
-  AliPHOSEmcRecPoint * EmcClu ;
+
+//____________________________________________________________________________
+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
   
-  Int_t iLinkLow = 0 ;
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
+  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
+
+  fLinkUpArray->Clear() ;    
+
+  AliPHOSCpvRecPoint * cpv ;
+  AliPHOSEmcRecPoint * emcclu ;
+
   Int_t iLinkUp  = 0 ;
   
-  while( (EmcClu = (AliPHOSEmcRecPoint*)nextEmc() ) ) {
-    Bool_t TooFar ;
-    TIter nextPpsdLow(PpsdRecPointsLow ) ;
-    Int_t iPpsdLow = 0 ;
-    
-    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++ ;
-      
-    }
-    
-    TIter nextPpsdUp(PpsdRecPointsUp ) ;
-    Int_t iPpsdUp = 0 ;
-    
-    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++ ;
+  Int_t iEmcRP;
+  for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
+    emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP)) ;
+
+    //Bool_t toofar ;        
+    Int_t iCpv = 0 ;    
+    for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) { 
       
+      cpv = dynamic_cast<AliPHOSCpvRecPoint *>(cpvRecPoints->At(iCpv)) ;
+      Int_t track = -1 ; 
+      Float_t r = GetDistanceInPHOSPlane(emcclu, cpv, track) ;     
+      //      if(toofar)
+      //       continue ;       
+      if(r < fRcpv) { 
+        new ((*fLinkUpArray)[iLinkUp++])  AliPHOSLink(r, iEmcRP, iCpv, track) ;
+      }      
     }
-    
-    iEmcClu++ ; 
-    
-  } // while nextEmC
+  } 
   
-  LinkLowArray->Sort() ; //first links with smallest distances
-  LinkUpArray->Sort() ;
+  fLinkUpArray->Sort() ;  //first links with smallest distances
 }
-    
+
 //____________________________________________________________________________
-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) ;
+void  AliPHOSTrackSegmentMakerv1::MakePairs()
+{ 
+  // 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". 
+
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+
+  TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
+  TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
+  TClonesArray * trackSegments = gime->TrackSegments();
+    
+  //Make arrays to mark clusters already chosen
+  Int_t * emcExist = 0;
+  if(fEmcLast > fEmcFirst)
+    emcExist = new Int_t[fEmcLast-fEmcFirst] ;
+  
+  Int_t index;
+  for(index = 0; index <fEmcLast-fEmcFirst; index ++)
+    emcExist[index] = 1 ;
+  
+  Bool_t * cpvExist = 0;
+  if(fCpvLast > fCpvFirst)
+    cpvExist = new Bool_t[fCpvLast-fCpvFirst] ;
+  for(index = 0; index <fCpvLast-fCpvFirst; index ++)
+    cpvExist[index] = kTRUE ;
+  
+  
+  // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance 
+  TIter nextUp(fLinkUpArray) ;
   
-  AliPHOSLink * linkLow ;
   AliPHOSLink * linkUp ;
+  
+  AliPHOSCpvRecPoint * nullpointer = 0 ;
+  
+  while ( (linkUp =  static_cast<AliPHOSLink *>(nextUp()) ) ){  
 
-  AliPHOSEmcRecPoint * emc ;
-  AliPHOSPpsdRecPoint * ppsdLow ;
-  AliPHOSPpsdRecPoint * ppsdUp ;
-
-  AliPHOSRecPoint * NullPointer = 0 ;
+    if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
 
-  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 ;
-          }
-       
-        } // while nextUp
-        
-        nextUp.Reset();
-         AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
-        trsl->Add(subtr) ;  
-        EmcRecPoints->AddAt(NullPointer,linkLow->GetEmc()) ;     
-        PpsdRecPointsLow->AddAt(NullPointer,linkLow->GetPpsd()) ;
+      if(cpvExist[linkUp->GetCpv()-fCpvFirst]){ //CPV still exist
+        new ((* trackSegments)[fNTrackSegments]) 
+          AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(linkUp->GetEmc())) , 
+                              dynamic_cast<AliPHOSCpvRecPoint *>(cpvRecPoints->At(linkUp->GetCpv())) , 
+                              linkUp->GetTrack()) ;
         
-        if(ppsdUp)  
-          PpsdRecPointsUp->AddAt(NullPointer,linkUp->GetPpsd()) ;
-        
-    } // if NLocMax
-  } 
-   
-  TIter nextEmc(EmcRecPoints) ;          
-  nextEmc.Reset() ;
-
-  while( (emc = (AliPHOSEmcRecPoint*)nextEmc()) ){ //to create pairs if no PpsdLow
-    ppsdLow = 0 ; 
-    ppsdUp  = 0 ;
-    
-    while ( (linkUp =  (AliPHOSLink *)nextUp() ) ){
-      
-      if(EmcRecPoints->IndexOf(emc) == linkUp->GetEmc() ){
-       ppsdUp = (AliPHOSPpsdRecPoint *) PpsdRecPointsUp->At(linkUp->GetPpsd()) ;
-       break ;
-      }
-      
+       (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->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 ){
+       new ((*trackSegments)[fNTrackSegments])  
+         AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP+fEmcFirst)), 
+                           nullpointer) ;
+       (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
+       fNTrackSegments++;    
+      } 
     }
-    AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
-    trsl->Add(subtr) ;   
-    if(ppsdUp)  
-      PpsdRecPointsUp->AddAt(NullPointer,linkUp->GetPpsd()) ;
   }
-     
+  delete [] emcExist ; 
+  delete [] cpvExist ; 
 }
 
 //____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::MakeTrackSegments(DigitsList * DL, RecPointsList * emcl, 
-                                       RecPointsList * ppsdl, TrackSegmentsList * trsl)
+void  AliPHOSTrackSegmentMakerv1::Exec(Option_t *option)
 {
-  // main function, does the job
-
-  Int_t PHOSMod      = 1 ;
-  Int_t emcStopedAt  = 0 ; 
-  Int_t ppsdStopedAt = 0 ; 
-  
-  TObjArray * EmcRecPoints     = new TObjArray(100) ;  // these arrays keep pointers 
-  TObjArray * PpsdRecPointsUp  = new TObjArray(100) ;  // to 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); 
+  // Steering method to perform track segment construction for events
+  // in the range from fFirstEvent to fLastEvent.
+  // This range is optionally set by SetEventRange().
+  // if fLastEvent=-1 (by default), then process events until the end.
   
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-  
-  while(PHOSMod <= geom->GetNModules() ){
-
-    cout << PHOSMod << " Track1 " << endl ;
-    
-    FillOneModule(DL, emcl, EmcRecPoints, ppsdl, PpsdRecPointsUp, PpsdRecPointsLow, PHOSMod , emcStopedAt, ppsdStopedAt) ;
-    cout << PHOSMod << " Track2 " << endl ;
-    MakeLinks(EmcRecPoints, PpsdRecPointsUp, PpsdRecPointsLow, LinkLowArray, LinkUpArray) ; 
-
-    cout << PHOSMod << " Track3 " << endl ;
-    MakePairs(EmcRecPoints, PpsdRecPointsUp, PpsdRecPointsLow, LinkLowArray, LinkUpArray, trsl) ;
+  if(strstr(option,"tim"))
+    gBenchmark->Start("PHOSTSMaker");
  
-    EmcRecPoints->Clear() ;
-    PpsdRecPointsUp->Clear() ;
-  
-    PpsdRecPointsLow->Clear() ;
-
-    LinkUpArray->Clear();
-   
-    LinkLowArray->Clear();
-   
+  if(strstr(option,"print")) {
+    Print() ; 
+    return ; 
   }
-  delete EmcRecPoints ; 
-  EmcRecPoints = 0 ; 
+  
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ;  
+  const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
 
-  delete PpsdRecPointsUp ; 
-  PpsdRecPointsUp = 0 ; 
+  if (fLastEvent == -1) 
+    fLastEvent = gime->MaxEvent() - 1 ;
+  else 
+    fLastEvent = TMath::Min(fFirstEvent,gime->MaxEvent());
+  Int_t nEvents   = fLastEvent - fFirstEvent + 1;
+
+  Int_t ievent ; 
+  for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
+    gime->Event(ievent,"R") ;
+   //Make some initializations 
+    fNTrackSegments = 0 ;
+    fEmcFirst = 0 ;    
+    fEmcLast  = 0 ;   
+    fCpvFirst = 0 ;   
+    fCpvLast  = 0 ;   
+    
+    gime->TrackSegments()->Clear();
 
-  delete PpsdRecPointsLow ; 
-  PpsdRecPointsLow = 0 ; 
+    //    if(!ReadRecPoints(ievent))   continue; //reads RecPoints for event ievent
+    
+    for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ) {
+      FillOneModule() ; 
+      MakeLinks() ;
+      MakePairs() ;
+    }
 
-  delete LinkUpArray ; 
-  LinkUpArray = 0  ; 
+    WriteTrackSegments() ;
 
-  delete LinkLowArray ; 
-  LinkLowArray = 0 ; 
+    if(strstr(option,"deb"))
+      PrintTrackSegments(option);
+    
+    //increment the total number of track segments per run 
+    fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ; 
+  }
+  
+  if(strstr(option,"tim")){
+    gBenchmark->Stop("PHOSTSMaker");
+    Info("Exec", "took %f seconds for making TS %f seconds per event", 
+          gBenchmark->GetCpuTime("PHOSTSMaker"), 
+          gBenchmark->GetCpuTime("PHOSTSMaker")/nEvents) ;
+   }
+  if(fWrite) //do not unload in "on flight" mode
+    Unload();
 }
 
 //____________________________________________________________________________
-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::Unload() 
+{
+  // Unloads the task from the folder
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ;  
+  gime->PhosLoader()->UnloadRecPoints() ;
+  gime->PhosLoader()->UnloadTracks() ;
 }
 
 //____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::UnfoldClusters(DigitsList * DL, RecPointsList * emcIn,  AliPHOSEmcRecPoint * iniEmc, 
-                                        Int_t Nmax, int * maxAt, Float_t * maxAtEnergy, TObjArray * emcList)
-{ 
-  // fits cluster with Nmax overlapping showers 
-  
-  Int_t NPar = 3 * Nmax ;
-  Float_t FitParameters[NPar] ;
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-
-  Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, NPar, FitParameters) ;
-  if( !rv )  // Fit failed, return and remove cluster
-    return ;
-  
-  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) ;
-    }
-
+void AliPHOSTrackSegmentMakerv1::Print(const Option_t *)const
+{
+  //  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 ") ;
+}
 
-  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 ) ;
-    }
+//____________________________________________________________________________
+void AliPHOSTrackSegmentMakerv1::WriteTrackSegments()
+{
+  // Writes found TrackSegments to TreeR. Creates branches 
+  // "PHOSTS" and "AliPHOSTrackSegmentMaker" with the same title.
+  // In the former branch found TrackSegments are stored, while 
+  // in the latter all parameters, with which TS were made. 
+  // ROOT does not allow overwriting existing branches, therefore
+  // first we check, if branches with the same title already exist.
+  // If yes - exits without writing.
 
-    emcList->Add(emcRP) ;
+  AliPHOSGetter *gime = AliPHOSGetter::Instance() ; 
 
+  TClonesArray * trackSegments = gime->TrackSegments() ; 
+  trackSegments->Expand(trackSegments->GetEntriesFast()) ;
+
+  if(fWrite){ //We write TreeT
+    TTree * treeT = gime->TreeT();
+    
+    //First TS
+    Int_t bufferSize = 32000 ; 
+    TBranch * tsBranch = treeT->Branch("PHOSTS",&trackSegments,bufferSize);
+    tsBranch->Fill() ;  
+    
+    gime->WriteTracks("OVERWRITE");
+    gime->WriteTrackSegmentMaker("OVERWRITE");
   }
 }
-//______________________________________________________________________________
-void UnfoldingChiSquare(Int_t &NPar, Double_t *Grad, Double_t & fret, Double_t *x, Int_t iflag)
-{
 
-// Number of paramters, Gradient , Chi squared, parameters, what to do
 
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
+//____________________________________________________________________________
+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
 
-  AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit
-  int * emcDigits = emcRP->GetDigitsList() ;
-  Float_t * emcEnergies = emcRP->GetEnergiesList() ;
-  fret = 0. ;     
-  Int_t iparam ;
+  TClonesArray * trackSegments = AliPHOSGetter::Instance()->TrackSegments() ; 
 
-  if(iflag == 2)
-    for(iparam = 0 ; iparam < NPar ; iparam++)    
-      Grad[iparam] = 0 ; // Will evaluate gradient
-  
-  Double_t Efit ;  
+  Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ; 
+  printf("nevent: %d\n", gAlice->GetEvNumber()) ; 
+  printf("        Found %d TrackSegments\n", trackSegments->GetEntriesFast() ); 
   
-  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 = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
-        iParam++ ; 
-        Distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ; 
-        Distance = TMath::Sqrt( Distance ) ; 
-        iParam++ ;      
-        Efit += x[iParam] * AliPHOSTrackSegmentMakerv1::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 * AliPHOSTrackSegmentMakerv1::ShowerShape(dr) ;
-        Double_t r4 = dr*dr*dr*dr ;
-        Double_t r295 = TMath::Power(dr,2.95) ;
-        Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
-                                        0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
-        
-        Grad[iParam] += Epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x    
-        iParam++ ; 
-        Grad[iParam] += Epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z         
-        iParam++ ; 
-        Grad[iParam] += shape ;                                  // Derivative over energy             
-        iParam++ ; 
-       }
-     }
-     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 * AliPHOSTrackSegmentMakerv1::ShowerShape(Distance) ;
-     }
-     fret += (Efit-emcEnergies[iDigit])*(Efit-emcEnergies[iDigit])/emcEnergies[iDigit] ; 
-     // Here we assume, that sigma = sqrt(E)
-     iDigit++ ;
+  if(strstr(option,"all")) {  // printing found TS
+    printf("TrackSegment #  EMC RP#  CPV RP#\n") ; 
+    Int_t index;
+    for (index = 0 ; index <trackSegments->GetEntriesFast() ; index++) {
+      AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )trackSegments->At(index) ; 
+      printf("   %d           %d        %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetCpvIndex() ) ; 
+    }  
   }
 }