]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PHOS/AliPHOSTrackSegmentMakerv1.cxx
Interference with other detectors removed: read/fill only our branches
[u/mrichter/AliRoot.git] / PHOS / AliPHOSTrackSegmentMakerv1.cxx
index 07303afa218c6ef9da8bfcb7c3c0cfbd44350bca..a2cab4beccc3c68c8a7b55bc173d64e685f69f12 100644 (file)
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-
 /* $Id$ */
-
 //_________________________________________________________________________
 // Implementation version 1 of algorithm class to construct PHOS track segments
-// Associates EMC and PPSD clusters
-// Unfolds the EMC cluster   
-//                  
+// Track segment for PHOS is list of 
+//        EMC RecPoint + (possibly) CPV RecPoint + (possibly) PPSD RecPoint
+// To find TrackSegments we do the following: for each EMC RecPoints we look at
+// CPV/PPSD RecPoints in the radious fR0. If there is such a CPV RecPoint, 
+// we make "Link" it is just indexes of EMC and CPV/PPSD 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/PPSD RecPoint we make TrackSegment 
+// consisting from EMC along. There is no TrackSegments without EMC RecPoint.
+//
+// In principle this class should be called from AliPHOSReconstructioner, but 
+// one can use it as well in standalong mode.
+// User case:
+// root [0] AliPHOSTrackSegmentMakerv1 * t = new AliPHOSTrackSegmentMaker("galice.root")
+// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
+// root [1] t->ExecuteTask()
+// root [2] t->SetMaxEmcPpsdDistance(5)
+// root [3] t->SetTrackSegmentsBranch("max distance 5 cm")
+// root [4] t->ExecuteTask("deb all time") 
+//                 
 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH)
 //
 
 // --- ROOT system ---
-
-#include "TObjArray.h"
-#include "TClonesArray.h"
-#include "TObjectTable.h"
-
+#include "TROOT.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TSystem.h"
+#include "TBenchmark.h"
 // --- Standard library ---
 
-#include <iostream>
-#include <cassert>
+#include <iostream.h>
+#include <iomanip.h>
 
 // --- AliRoot header files ---
 
 #include "AliPHOSTrackSegmentMakerv1.h"
+#include "AliPHOSClusterizerv1.h"
 #include "AliPHOSTrackSegment.h"
+#include "AliPHOSCpvRecPoint.h"
+#include "AliPHOSPpsdRecPoint.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) ; 
-
 ClassImp( AliPHOSTrackSegmentMakerv1) 
 
 
@@ -51,553 +67,579 @@ ClassImp( AliPHOSTrackSegmentMakerv1)
  AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() : AliPHOSTrackSegmentMaker()
 {
   // 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) ;
-  fUnfoldFlag = kTRUE ; 
+  SetTitle("version 1") ;
+  SetName("AliPHOSTrackSegmentMaker") ;
+  fR0 = 10. ;   
+  fEmcFirst = 0 ;    
+  fEmcLast  = 0 ;   
+  fCpvFirst = 0 ;   
+  fCpvLast  = 0 ;   
+  fPpsdFirst= 0 ;   
+  fPpsdLast = 0 ;   
+  fLinkLowArray = 0 ;
+  fLinkUpArray  = 0 ;
+  fIsInitialized = kFALSE ;
 }
-
 //____________________________________________________________________________
- AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
-{ 
-  // dtor
-   delete fMinuit ; 
-}
+ AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const char* headerFile, const char* branchTitle): 
+AliPHOSTrackSegmentMaker()
+{
+  // ctor
+  SetTitle("version 1") ;
+  SetName("AliPHOSTrackSegmentMaker") ;
+  fR0 = 10. ;   
+  fEmcFirst = 0 ;    
+  fEmcLast  = 0 ;   
+  fCpvFirst = 0 ;   
+  fCpvLast  = 0 ;   
+  fPpsdFirst= 0 ;   
+  fPpsdLast = 0 ;   
+
+  fHeaderFileName = headerFile ;
+  fRecPointsBranchTitle = branchTitle ;
+    
+  TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ;
+  
+  if(file == 0){
+    file = new TFile(fHeaderFileName.Data(),"update") ;
+    gAlice = (AliRun *) file->Get("gAlice") ;
+  }
+  
+  AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ;    
+  fGeom  = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() );
+  
+  fEmcRecPoints = new TObjArray(200) ;
+  fCpvRecPoints = new TObjArray(200) ;
+  fClusterizer  = new AliPHOSClusterizerv1() ;
+  
+  fTrackSegments = new TClonesArray("AliPHOSTrackSegment",200) ;
+  
+  fLinkLowArray = new TClonesArray("AliPHOSLink", 1000);
+  fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000); 
+  
+  fIsInitialized = kTRUE ;
 
+}
 //____________________________________________________________________________
-Bool_t  AliPHOSTrackSegmentMakerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy,
-                                   Int_t nPar, Float_t * fitparameters)
-{ 
-  // Calls TMinuit to fit the energy distribution of a cluster with several maxima 
-
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-
-  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 ;
+void  AliPHOSTrackSegmentMakerv1::Init(){
+  //Make all memory allokations not possible in default constructor
 
+  if(!fIsInitialized){
+    if(fHeaderFileName.IsNull())
+      fHeaderFileName = "galice.root" ;
+    
+    
+    TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ;
+    
+    if(file == 0){
+      file = new TFile(fHeaderFileName.Data(),"update") ;
+      gAlice = (AliRun *) file->Get("gAlice") ;
+    }
+    
+    AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ;    
+    fGeom  = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() );
 
-  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] ;
+    fEmcRecPoints = new TObjArray(200) ;
+    fCpvRecPoints = new TObjArray(200) ;
+    fClusterizer  = new AliPHOSClusterizerv1() ;
 
-    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;
-    }
-  }
+    
+    fTrackSegments = new TClonesArray("AliPHOSTrackSegment",200) ;
 
-  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 ;
+    fLinkLowArray = new TClonesArray("AliPHOSLink", 1000);
+    fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000); 
+    
+    fIsInitialized = kTRUE ;
    }
-  return kTRUE;
+}
 
+//____________________________________________________________________________
+ AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
+{ 
+  // dtor
+  if(fLinkLowArray) delete fLinkLowArray ;
+  if(fLinkUpArray)  delete fLinkUpArray  ;
 }
 
 //____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::FillOneModule(DigitsList * dl, RecPointsList * emcIn, TObjArray * emcOut, 
-                                       RecPointsList * ppsdIn, TObjArray * ppsdOutUp,
-                                       TObjArray * ppsdOutLow, Int_t & phosmod, Int_t & emcStopedAt, 
-                                       Int_t & ppsdStopedAt)
+void  AliPHOSTrackSegmentMakerv1::FillOneModule()
 {
-  // Unfold clusters and fill xxxOut arrays with clusters from one PHOS module
+  // Finds first and last indexes between which 
+  // clusters from one PHOS module are
  
-  AliPHOSEmcRecPoint *  emcRecPoint  ; 
-  AliPHOSPpsdRecPoint * ppsdRecPoint ;
-  Int_t index ;
+
+  //First EMC clusters
+  Int_t totalEmc = fEmcRecPoints->GetEntriesFast() ;
+  for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&  
+       (((AliPHOSRecPoint *) fEmcRecPoints->At(fEmcLast))->GetPHOSMod() == fModule ); 
+      fEmcLast ++)  ;
+  
   
-  Int_t nEmcUnfolded = emcIn->GetEntries() ;
-  for(index = emcStopedAt; index < nEmcUnfolded; index++){
-    emcRecPoint = (AliPHOSEmcRecPoint *) (*emcIn)[index] ;
+  //Now CPV clusters
+  Int_t totalCpv = fCpvRecPoints->GetEntriesFast() ;
 
-    if(emcRecPoint->GetPHOSMod() != phosmod )  
-       break ;
+  if(fModule <= fGeom->GetNCPVModules()){ // in CPV geometry
     
-    Int_t nMultipl = emcRecPoint->GetMultiplicity() ; 
-    Int_t * maxAt = new Int_t[nMultipl] ;
-    Float_t * maxAtEnergy = new Float_t[nMultipl] ;
-    Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy) ;
-
-    if(nMax <= 1 )     // if cluster is very flat (no pronounced maximum) then nMax = 0 
-      emcOut->Add(emcRecPoint) ;
-    else if (fUnfoldFlag) {
-      UnfoldClusters(dl, emcIn, emcRecPoint, nMax, maxAt, maxAtEnergy, emcOut) ;
-      emcIn->Remove(emcRecPoint); 
-      emcIn->Compress() ;
-      nEmcUnfolded-- ;
-      index-- ;
-    }
+    for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) && 
+         (((AliPHOSRecPoint *) fCpvRecPoints->At(fCpvLast))->GetPHOSMod() == fModule ); 
+       fCpvLast ++) ;
     
-    delete[] maxAt ; 
-    delete[] maxAtEnergy ; 
+    fPpsdFirst = fCpvLast ; //To avoid scanning RecPoints between fPpsdFirst and fPpsdLast
+    fPpsdLast  = fCpvLast ; //and to be ready to switch to mixed geometry 
   }
-  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) ;
+  else{  //in PPSD geometry    
+    fCpvLast = fPpsdLast ;
+    //Upper layer first
+    for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) &&  
+         (((AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fCpvLast))->GetPHOSMod() == fModule ) &&
+         (((AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fCpvLast))->GetUp()) ; 
+       fCpvLast ++)  ;
+    
+    fPpsdLast= fCpvLast ;
+    for(fPpsdFirst = fPpsdLast; (fPpsdLast < totalCpv)  &&
+         (((AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fPpsdLast))->GetPHOSMod() == fModule ) &&
+         (!((AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fPpsdLast))->GetUp()) ; 
+       fPpsdLast ++) ;
   }
-  ppsdStopedAt = index ;
-   
-  emcOut->Sort() ;
-  ppsdOutUp->Sort() ;
-  ppsdOutLow->Sort() ;   
+    
 }
 //____________________________________________________________________________
-Float_t  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcclu,AliPHOSPpsdRecPoint * PpsdClu, Bool_t &toofar)
+Float_t  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,AliPHOSRecPoint * cpvClu, Bool_t &toofar)
 {
   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
+  //clusters are sorted in "rows" and "columns" of width 1 cm
+
+  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 r = fR0 ;
  
   TVector3 vecEmc ;
-  TVector3 vecPpsd ;
-  
-  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
+  TVector3 vecCpv ;
+  
+  emcClu->GetLocalPosition(vecEmc) ;
+  cpvClu->GetLocalPosition(vecCpv)  ; 
+
+  if(emcClu->GetPHOSMod() == cpvClu->GetPHOSMod()){ 
+    if(vecCpv.X() <= vecEmc.X() + fR0 + 2*delta ){ 
+
+      vecCpv = vecCpv  - vecEmc ; 
+      r = vecCpv.Mag() ;
       toofar = kFALSE ;
-    } // if  xPpsd >= xEmc - fDelta
+
+    } // if  xPpsd >= xEmc + ...
     else 
       toofar = kTRUE ;
   } 
   else 
     toofar = kTRUE ;
+
+  //toofar = kFALSE ;
   
   return r ;
 }
 
 //____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::MakeLinks(TObjArray * emcRecPoints, TObjArray * ppsdRecPointsUp, 
-                                    TObjArray * ppsdRecPointsLow, TClonesArray * linklowArray, 
-                                    TClonesArray *linkupArray) 
+void  AliPHOSTrackSegmentMakerv1::MakeLinks()
 { 
-  // Finds distances (links) between all EMC and PPSD clusters, which are not further apart from each other than fR0 
-  
-  TIter nextEmc(emcRecPoints) ;
-  Int_t iEmcClu = 0 ; 
+  // Finds distances (links) between all EMC and PPSD clusters, 
+  // which are not further apart from each other than fR0 
+  // and sort them in accordance with this distance
   
-  AliPHOSPpsdRecPoint * ppsdlow ; 
-  AliPHOSPpsdRecPoint * ppsdup ;
+  fLinkUpArray->Clear() ;    
+  fLinkLowArray->Clear() ;
+
+  AliPHOSRecPoint * ppsd ; 
+  AliPHOSRecPoint * cpv ;
   AliPHOSEmcRecPoint * emcclu ;
-  
+
   Int_t iLinkLow = 0 ;
   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) ;
+  Int_t iEmcRP;
+  for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
+    emcclu = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(iEmcRP) ;
+
+    Bool_t toofar ;    
+    Int_t iPpsd ;
+    for(iPpsd = fPpsdFirst; iPpsd < fPpsdLast;iPpsd++ ) {
       
+      ppsd = (AliPHOSRecPoint *) fCpvRecPoints->At(iPpsd) ;
+      Float_t r = GetDistanceInPHOSPlane(emcclu, ppsd, toofar) ;
+
       if(toofar) 
        break ;  
-      if(r < fR0) 
-       new( (*linklowArray)[iLinkLow++]) AliPHOSLink(r, iEmcClu, iPpsdLow) ;
-      
-      iPpsdLow++ ;
-      
+      if(r < fR0)
+       new ((*fLinkLowArray)[iLinkLow++])  AliPHOSLink(r, iEmcRP, iPpsd) ;
     }
     
-    TIter nextPpsdUp(ppsdRecPointsUp ) ;
-    Int_t iPpsdUp = 0 ;
-    
-    while( (ppsdup = (AliPHOSPpsdRecPoint*)nextPpsdUp() ) ) { 
-      Float_t r = GetDistanceInPHOSPlane(emcclu, ppsdup, toofar) ;
+    Int_t iCpv = 0 ;    
+    for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) { 
+      
+      cpv = (AliPHOSRecPoint *) fCpvRecPoints->At(iCpv) ;
+      Float_t r = GetDistanceInPHOSPlane(emcclu, cpv, toofar) ;
       
       if(toofar)
        break ;  
-      if(r < fR0) 
-       new( (*linkupArray)[iLinkUp++]) AliPHOSLink(r, iEmcClu, iPpsdUp) ;
-      
-      iPpsdUp++ ;
-      
+      if(r < fR0) { 
+       new ((*fLinkUpArray)[iLinkUp++])  AliPHOSLink(r, iEmcRP, iCpv) ;
+      }      
     }
-    
-    iEmcClu++ ; 
-    
-  } // while nextEmC
+  } 
   
-  linklowArray->Sort() ; //first links with smallest distances
-  linkupArray->Sort() ;
+  fLinkLowArray->Sort() ; //first links with smallest distances
+  fLinkUpArray->Sort() ;
 }
-    
+
 //____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::MakePairs(TObjArray * emcRecPoints, TObjArray * ppsdRecPointsUp, 
-                                   TObjArray * ppsdRecPointsLow, TClonesArray * linklowArray, 
-                                   TClonesArray * linkupArray, TrackSegmentsList * trsl) 
+void  AliPHOSTrackSegmentMakerv1::MakePairs()
 { 
-
-  // Finds the smallest links and makes pairs of PPSD and EMC clusters with smallest distance 
+  // Using the previously made list of "links", we found the smallest link - i.e. 
+  // link with the least distance betwing 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 chousen
+  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 ;
+  
+  Bool_t * ppsdExist = 0;
+  if(fPpsdLast > fPpsdFirst)
+    ppsdExist = new Bool_t[fPpsdLast-fPpsdFirst] ;
+  for(index = 0; index <fPpsdLast-fPpsdFirst; index ++)
+    ppsdExist[index] = kTRUE ;
   
-  TIter nextLow(linklowArray) ;
-  TIter nextUp(linkupArray) ;
+  // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance 
+  TIter nextLow(fLinkLowArray) ;
+  TIter nextUp(fLinkUpArray) ;
   
   AliPHOSLink * linkLow ;
   AliPHOSLink * linkUp ;
 
-  AliPHOSEmcRecPoint * emc ;
-  AliPHOSPpsdRecPoint * ppsdLow ;
-  AliPHOSPpsdRecPoint * ppsdUp ;
 
   AliPHOSRecPoint * nullpointer = 0 ;
 
   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) ;  
-        new( (*trsl)[fNTrackSegments] ) AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
-        fNTrackSegments++ ;
-        emcRecPoints->AddAt(nullpointer,linkLow->GetEmc()) ;     
-        ppsdRecPointsLow->AddAt(nullpointer,linkLow->GetPpsd()) ;
-        
-        if(ppsdUp)  
-          ppsdRecPointsUp->AddAt(nullpointer,linkUp->GetPpsd()) ;
+  
+    if( (emcExist[linkLow->GetEmc()-fEmcFirst]> 0) && 
+       ppsdExist[linkLow->GetPpsd()-fPpsdFirst]  ){ // RecPoints not removed yet 
+      new ((*fTrackSegments)[fNTrackSegments]) AliPHOSTrackSegment((AliPHOSEmcRecPoint *) fEmcRecPoints->At(linkLow->GetEmc()), 
+                                                nullpointer, 
+                                               (AliPHOSRecPoint *)fCpvRecPoints->At(linkLow->GetPpsd()) ) ;
         
+      ((AliPHOSTrackSegment* )fTrackSegments->At(fNTrackSegments))->SetIndexInList(fNTrackSegments);    
+      //replace index of emc to negative and shifted index of TS      
+      emcExist[linkLow->GetEmc()-fEmcFirst] = -2 - fNTrackSegments ;  
+      //mark ppsd as used
+      ppsdExist[linkLow->GetPpsd()-fPpsdFirst] = kFALSE ; 
+      fNTrackSegments++ ;
     } 
   } 
-   
-  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 ;
-      }
-      
-    }
-//     AliPHOSTrackSegment * subtr = new AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
-//     trsl->Add(subtr) ;   
-    new( (*trsl)[fNTrackSegments] ) AliPHOSTrackSegment(emc, ppsdUp, ppsdLow ) ;
-    fNTrackSegments++ ;
-    
+  while ( (linkUp =  (AliPHOSLink *)nextUp() ) ){  
+    if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){ //without ppsd Up yet 
 
-    if(ppsdUp)  
-      ppsdRecPointsUp->AddAt(nullpointer,linkUp->GetPpsd()) ;
+      if(cpvExist[linkUp->GetPpsd()-fCpvFirst]){ //CPV still exist
+       
+       if(emcExist[linkUp->GetEmc()-fEmcFirst] > 0){ //without ppsd Low => create new TS
+
+         new ((* fTrackSegments)[fNTrackSegments]) AliPHOSTrackSegment((AliPHOSEmcRecPoint *) fEmcRecPoints->At(linkUp->GetEmc()) , 
+                                                                     (AliPHOSRecPoint *)fCpvRecPoints->At(linkUp->GetPpsd()), 
+                                                                     nullpointer) ;
+         ((AliPHOSTrackSegment *) fTrackSegments->At(fNTrackSegments))->SetIndexInList(fNTrackSegments);
+         fNTrackSegments++ ;
+       }
+       else{ // append ppsd Up to existing TS
+         ((AliPHOSTrackSegment *)fTrackSegments->At(-2-emcExist[linkUp->GetEmc()-fEmcFirst]))->SetCpvRecPoint((AliPHOSCpvRecPoint *)fCpvRecPoints->At(linkUp->GetPpsd()));
+       }
+
+       emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc  that Cpv was found 
+       //mark CPV recpoint as already used 
+        cpvExist[linkUp->GetPpsd()-fCpvFirst] = kFALSE ;
+      } //if ppsdUp still exist
+    } 
+  }     
+
+  //look through emc recPoints left without CPV/PPSD
+  if(emcExist){ //if there is emc rec point
+    Int_t iEmcRP ;
+    for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst  ; iEmcRP++ ){
+      if(emcExist[iEmcRP] > 0 ){
+       new ((*fTrackSegments)[fNTrackSegments])  AliPHOSTrackSegment((AliPHOSEmcRecPoint *)fEmcRecPoints->At(iEmcRP+fEmcFirst), 
+                                                                   nullpointer, 
+                                                                   nullpointer ) ;
+       ((AliPHOSTrackSegment *) fTrackSegments->At(fNTrackSegments))->SetIndexInList(fNTrackSegments);
+       fNTrackSegments++;    
+      } 
+    }
   }
-     
+  
 }
 
 //____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::MakeTrackSegments(DigitsList * dl, RecPointsList * emcl, 
-                                       RecPointsList * ppsdl, TrackSegmentsList * trsl)
+void  AliPHOSTrackSegmentMakerv1::Exec(Option_t * option)
 {
-  // Makes the track segments out of the list of EMC and PPSD Recpoints and stores them in a list
+  //STEERing method
 
-  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); 
-  
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-  
-  while(phosmod <= geom->GetNModules() ){
-    
-    FillOneModule(dl, emcl, emcRecPoints, ppsdl, ppsdRecPointsUp, ppsdRecPointsLow, phosmod, emcStopedAt, ppsdStopedAt) ;
+  if(! fIsInitialized) Init() ;
 
-    MakeLinks(emcRecPoints, ppsdRecPointsUp, ppsdRecPointsLow, linklowArray, linkupArray) ; 
+  if(strstr(option,"tim"))
+    gBenchmark->Start("PHOSTSMaker");  
 
-    MakePairs(emcRecPoints, ppsdRecPointsUp, ppsdRecPointsLow, linklowArray, linkupArray, trsl) ;
-    emcRecPoints->Clear() ;
-    ppsdRecPointsUp->Clear() ;
+  Int_t nEvents = (Int_t) gAlice->TreeE()->GetEntries() ;
   
-    ppsdRecPointsLow->Clear() ;
+  for(fEvent = 0;fEvent< nEvents; fEvent++){
+    if(!ReadRecPoints())  //reads RecPoints for event fEvent
+      return;
+    
+    for(fModule = 1; fModule <= fGeom->GetNModules() ; fModule++ ){
+      
+      FillOneModule() ; 
+      
+      MakeLinks() ;
+      
+      MakePairs() ;
+      
+    }
 
-    linkupArray->Clear() ;
-   
-    linklowArray->Clear() ;
-   
-    phosmod++ ; 
+    WriteTrackSegments() ;
+    if(strstr(option,"deb"))
+      PrintTrackSegments(option) ;
   }
-  delete emcRecPoints ; 
-  emcRecPoints = 0 ; 
-
-  delete ppsdRecPointsUp ; 
-  ppsdRecPointsUp = 0 ; 
 
-  delete ppsdRecPointsLow ; 
-  ppsdRecPointsLow = 0 ; 
+  if(strstr(option,"tim")){
+    gBenchmark->Stop("PHOSTSMaker");
+    cout << "AliPHOSTSMaker:" << endl ;
+    cout << "  took " << gBenchmark->GetCpuTime("PHOSTSMaker") << " seconds for making TS " 
+        <<  gBenchmark->GetCpuTime("PHOSTSMaker")/nEvents << " seconds per event " << endl ;
+    cout << endl ;
+  }
 
-  delete linkupArray ; 
-  linkupArray = 0  ; 
 
-  delete linklowArray ; 
-  linklowArray = 0 ; 
 }
-
 //____________________________________________________________________________
-Double_t  AliPHOSTrackSegmentMakerv1::ShowerShape(Double_t r)
-{ 
-  // Shape of the shower (see PHOS TDR)
-  // If you change this function, change also the gradien 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::Print(Option_t * option)const {
+  if(fIsInitialized){
+    cout <<  "======== AliPHOSTrackSegmentMakerv1 ========" << endl ;
+    cout <<  "Making Track segments "<< endl ;
+    cout <<  "    Headers file:                   " << fHeaderFileName.Data() << endl ;
+    cout <<  "    RecPoints branch file name:     " << fRecPointsBranchTitle.Data() << endl ;
+    cout <<  "    TrackSegments Branch file name: " << fTSBranchTitle.Data() << endl ;
+    cout <<  "with parameters: " << endl ;
+    cout <<  "    Maximal EMC - CPV (PPSD) distance (cm)" << fR0 << endl ;
+    cout <<  "============================================" << endl ;
+  }
+  else
+    cout << "AliPHOSTrackSegmentMakerv1 not initialized " << endl ;
 }
-
 //____________________________________________________________________________
-void  AliPHOSTrackSegmentMakerv1::UnfoldClusters(DigitsList * dl, RecPointsList * emcIn,  AliPHOSEmcRecPoint * iniEmc, 
-                                        Int_t nMax, int * maxAt, Float_t * maxAtEnergy, TObjArray * emcList)
-{ 
-  // Performs the unfolding of a cluster with nMax overlapping showers 
-  // This is time consuming (use the (Un)SetUnfolFlag()  )
+Bool_t AliPHOSTrackSegmentMakerv1::ReadRecPoints(){
+  // Reads Emc and CPV recPoints with given title (fRecPointsBranchTitle) 
+  // made priveously with Clusterizer.
+
+
+  //Make some initializations 
+  fEmcRecPoints->Clear() ;
+  fCpvRecPoints->Clear() ;
+  fTrackSegments->Clear() ;
+  fNTrackSegments = 0 ;
+  fEmcFirst = 0 ;    
+  fEmcLast  = 0 ;   
+  fCpvFirst = 0 ;   
+  fCpvLast  = 0 ;   
+  fPpsdFirst= 0 ;   
+  fPpsdLast = 0 ;   
+
+
+  gAlice->GetEvent(fEvent) ;
+
+  // Get TreeR header from file
+  if(gAlice->TreeR()==0){
+    char treeName[20]; 
+    sprintf(treeName,"TreeR%d",fEvent);
+    cout << "Error in AliPHOSTrackSegmentMakerv1 : no "<<treeName << endl  ;
+    cout << "   Do nothing " << endl ;
+    return kFALSE ;
+  }
 
-  Int_t nPar = 3 * nMax ;
-  Float_t * fitparameters = new Float_t[nPar] ;
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
 
-  Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
-  if( !rv ) {
-    // Fit failed, return and remove cluster
-    delete[] fitparameters ; 
-    return ;
-  }
+  //Find RecPoints with title fRecPointsBranchTitle
+  TBranch * emcBranch = 0;
+  TBranch * cpvBranch = 0;
+  TBranch * clusterizerBranch = 0;
+
+  TObjArray * branches = gAlice->TreeR()->GetListOfBranches() ;
+  Int_t ibranch;
+  Bool_t emcNotFound = kTRUE ;
+  Bool_t cpvNotFound = kTRUE ;  
+  Bool_t clusterizerNotFound = 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 = new Float_t[nDigits] ;
-  Int_t iparam ;
-  Int_t iDigit ;
-  
-  AliPHOSDigit * digit ;
-  AliPHOSEmcRecPoint * emcRP ;  
-  Int_t * 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 ;
+  for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
+
+    if(emcNotFound){
+      emcBranch=(TBranch *) branches->At(ibranch) ;
+      if( fRecPointsBranchTitle.CompareTo(emcBranch->GetTitle())==0 )
+       if( strcmp(emcBranch->GetName(),"PHOSEmcRP") == 0) {
+         emcNotFound = kFALSE ;
+       }
+    }
+    
+    if(cpvNotFound){
+      cpvBranch=(TBranch *) branches->At(ibranch) ;
+      if( fRecPointsBranchTitle.CompareTo(cpvBranch->GetTitle())==0 )
+       if( strcmp(cpvBranch->GetName(),"PHOSCpvRP") == 0) 
+         cpvNotFound = kFALSE ;
+    }
     
-    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 * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) ;
+    if(clusterizerNotFound){
+      clusterizerBranch = (TBranch *) branches->At(ibranch) ;
+      if( fRecPointsBranchTitle.CompareTo(clusterizerBranch->GetTitle()) == 0)
+       if( strcmp(clusterizerBranch->GetName(),"AliPHOSClusterizer") == 0) 
+         clusterizerNotFound = kFALSE ;
     }
+    
   }
 
-  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 * AliPHOSTrackSegmentMakerv1::ShowerShape(distance) / efit[iDigit] ; 
-      eDigit = emcEnergies[iDigit] * ratio ;
-      emcRP->AddDigit( *digit, eDigit ) ;
+  if(clusterizerNotFound || emcNotFound || cpvNotFound){
+    cout << "AliPHOSTrackSegmentMakerv1: " << endl ;
+    cout << "    Can't find Branch with RecPoints or Clusterizer " ;
+    cout << "    Do nothing" <<endl  ;
+    return kFALSE ;
+  }
+  
+  emcBranch->SetAddress(&fEmcRecPoints) ;
+  cpvBranch->SetAddress(&fCpvRecPoints) ;
+  clusterizerBranch->SetAddress(&fClusterizer) ;
+
+  emcBranch->GetEntry(0) ;
+  cpvBranch->GetEntry(0) ;
+  clusterizerBranch->GetEntry(0) ;
+  
+  return kTRUE ;
+  
+}
+//____________________________________________________________________________
+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 chesk, if branches with the same title alredy exist.
+  // If yes - exits without writing.
+  
+  //First, check, if branches already exist
+  TBranch * tsMakerBranch = 0;
+  TBranch * tsBranch = 0;
+  
+  TObjArray * branches = gAlice->TreeR()->GetListOfBranches() ;
+  Int_t ibranch;
+  Bool_t tsMakerNotFound = kTRUE ;
+  Bool_t tsNotFound = kTRUE ;
+  
+  for(ibranch = 0;(ibranch <branches->GetEntries())&&(tsMakerNotFound||tsNotFound);ibranch++){
+    if(tsMakerNotFound){
+      tsMakerBranch=(TBranch *) branches->At(ibranch) ;
+      if( (strcmp(tsMakerBranch->GetName(),"AliPHOSTrackSegmentMaker") == 0) &&
+         (fTSBranchTitle.CompareTo( tsMakerBranch->GetTitle())==0 ))
+       tsMakerNotFound = kFALSE ;
+    }
+    if(tsNotFound){
+      tsBranch=(TBranch *) branches->At(ibranch) ;
+      if( (strcmp(tsBranch->GetName(),"PHOSTS") == 0)  &&
+         (fTSBranchTitle.CompareTo( tsBranch->GetTitle())==0 ))
+       tsNotFound = kFALSE ;
     }
+  }
 
-    emcList->Add(emcRP) ;
+  if(!(tsMakerNotFound && tsNotFound )){ 
+    cout << "AliPHOSTrackSegmentMakerv1 error:"<< endl ;
+    cout << "       Branches PHOSTS and AliPHOSTrackSegementMaker " << endl ;
+    cout << "       with title '"<<fTSBranchTitle.Data() << "' already exist " << endl ;
+    cout << "       can not overwrite " << endl ;
+    return ;
+  }
 
+  //Make branch in TreeR for TrackSegments 
+  char * filename = 0;
+  if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){   //generating file name
+    filename = new char[strlen(gAlice->GetBaseFile())+20] ;
+    sprintf(filename,"%s/PHOS.Reco.root",gAlice->GetBaseFile()) ; 
   }
-  
-  delete[] fitparameters ; 
-  delete[] efit ; 
 
+  TDirectory *cwd = gDirectory;
+  
+  //First TS
+  Int_t bufferSize = 32000 ;    
+  tsBranch = gAlice->TreeR()->Branch("PHOSTS",&fTrackSegments,bufferSize);
+  tsBranch->SetTitle(fTSBranchTitle.Data());
+  if (filename) {
+    tsBranch->SetFile(filename);
+    TIter next( tsBranch->GetListOfBranches());
+    TBranch * sb ;
+    while ((sb=(TBranch*)next())) {
+      sb->SetFile(filename);
+    }   
+    cwd->cd();
+  } 
+  
+  //Second -TSMaker
+  Int_t splitlevel = 0 ;
+  AliPHOSTrackSegmentMakerv1 * ts = this ;
+  tsMakerBranch = gAlice->TreeR()->Branch("AliPHOSTrackSegmentMaker","AliPHOSTrackSegmentMakerv1",
+                                         &ts,bufferSize,splitlevel);
+  tsMakerBranch->SetTitle(fTSBranchTitle.Data());
+  if (filename) {
+    tsMakerBranch->SetFile(filename);
+    TIter next( tsMakerBranch->GetListOfBranches());
+    TBranch * sb;
+    while ((sb=(TBranch*)next())) {
+      sb->SetFile(filename);
+    }   
+    cwd->cd();
+  } 
+  
+  tsBranch->Fill() ;  
+  tsMakerBranch->Fill() ;
+  gAlice->TreeR()->Write(0,kOverwrite) ;  
+  
 }
 
-//______________________________________________________________________________
-void UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
-{
-  // Calculates th Chi square for the cluster unfolding minimization
-  // Number of parameters, Gradient, Chi squared, parameters, what to do
-
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
 
-  AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit
-  Int_t * emcDigits     = emcRP->GetDigitsList() ;
-  Float_t * emcEnergies = emcRP->GetEnergiesList() ;
-  fret = 0. ;     
-  Int_t iparam ;
+//____________________________________________________________________________
+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
 
-  if(iflag == 2)
-    for(iparam = 0 ; iparam < nPar ; iparam++)    
-      Grad[iparam] = 0 ; // Will evaluate gradient
   
-  Double_t efit ;  
+  cout << "AliPHOSTrackSegmentMakerv1: " << endl ;
+  cout << "       Found " << fTrackSegments->GetEntriesFast() << "  trackSegments " << endl ;
   
-  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(strstr(option,"all")) {  // printing found TS
+    cout << "TrackSegment # " << "    EMC RP#    " << "    CPV RP#    " << "     PPSD RP#" << endl ; 
     
-     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++ ;
+    Int_t index;
+    for (index = 0 ; index <fTrackSegments->GetEntriesFast() ; index++) {
+      AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )fTrackSegments->At(index) ; 
+      cout<<"   "<< setw(4) << ts->GetIndexInList() << "            " 
+         <<setw(4) << ts->GetEmcIndex()<< "            " 
+         <<setw(4) << ts->GetCpvIndex()<< "            " 
+         <<setw(4) << ts->GetPpsdIndex()<< endl ;
+    }  
+    
+    cout << "-------------------------------------------------------"<< endl ;
   }
 }