X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSClusterFinderSPD.cxx;h=a56bf3cd2ccfe839e6d47a2b0a9f8e03da97ebff;hb=c63c3c5d365bcab29e17aca25ea8c04edc250e2d;hp=f24fa7382d9961e439cf0ae92d8db550b444b5a7;hpb=ebf6f0eed7e18233c40951b9c787daf05a0dd956;p=u%2Fmrichter%2FAliRoot.git

diff --git a/ITS/AliITSClusterFinderSPD.cxx b/ITS/AliITSClusterFinderSPD.cxx
index f24fa7382d9..a56bf3cd2cc 100644
--- a/ITS/AliITSClusterFinderSPD.cxx
+++ b/ITS/AliITSClusterFinderSPD.cxx
@@ -12,424 +12,419 @@
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
+////////////////////////////////////////////////////////////////////////////
+//  Cluster finder                                                       ///
+//  for Silicon pixels                                                    // 
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////        
 
-/*
-$Log$
-Revision 1.14  2001/03/05 14:48:46  nilsen
-Fixed a reoccuring bug. Array sizes must be declare const.
-
-*/
-
-#include <iostream.h>
 #include "AliITSClusterFinderSPD.h"
-#include "AliITSMapA1.h"
-#include "AliITS.h"
-#include "AliITSdigit.h"
-#include "AliITSRawCluster.h"
+#include "AliITSDetTypeRec.h"
+#include "AliITSRawClusterSPD.h"
 #include "AliITSRecPoint.h"
-#include "AliITSsegmentation.h"
-#include "AliITSresponse.h"
-#include "AliRun.h"
-
+#include "AliITSdigitSPD.h"
+#include "AliITSsegmentationSPD.h"
+#include "AliLog.h"
 
+//#define DEBUG
 
 ClassImp(AliITSClusterFinderSPD)
 
+//______________________________________________________________________
+AliITSClusterFinderSPD::AliITSClusterFinderSPD():AliITSClusterFinder(),
+fDz(0.0),
+fDx(0.0),
+fMinNCells(0){
+    // constructor
+}
 //----------------------------------------------------------
-AliITSClusterFinderSPD::AliITSClusterFinderSPD
-(AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp)   
-{
-  // constructor
-    fSegmentation=seg;
-    fDigits=digits;
-    fClusters=recp;
-    fNclusters= fClusters->GetEntriesFast();
+AliITSClusterFinderSPD::AliITSClusterFinderSPD(AliITSDetTypeRec* dettyp):
+AliITSClusterFinder(dettyp),
+fDz(0.0),
+fDx(0.0),
+fMinNCells(0){
+    // constructor
+
     SetDx();
     SetDz();
-    SetMap();
-    SetNCells();
 }
-
-//_____________________________________________________________________________
-AliITSClusterFinderSPD::AliITSClusterFinderSPD()
-{
-  // constructor
-  fSegmentation=0;
-  fDigits=0;
-  fClusters=0;
-  fNclusters=0;
-  fMap=0;
-  SetDx();
-  SetDz();
-  SetNCells();
-  
+//----------------------------------------------------------
+AliITSClusterFinderSPD::AliITSClusterFinderSPD(AliITSDetTypeRec* dettyp,
+                                               TClonesArray *digits,
+                                               TClonesArray *recp):
+AliITSClusterFinder(dettyp,digits),
+fDz(0.0),
+fDx(0.0),
+fMinNCells(0){
+    // constructor
+
+    SetDigits(digits);
+    SetClusters(recp);
+    SetDx();
+    SetDz();
 }
 
-//_____________________________________________________________________________
-AliITSClusterFinderSPD::~AliITSClusterFinderSPD()
-{
-  // destructor
-  if (fMap) delete fMap;
-
-
-}
-//__________________________________________________________________________
-AliITSClusterFinderSPD::AliITSClusterFinderSPD(const AliITSClusterFinderSPD &source){
-  //     Copy Constructor 
-  if(&source == this) return;
-  this->fClusters = source.fClusters ;
-  this->fNclusters = source.fNclusters ;
-  this->fMap = source.fMap ;
-  this->fDz = source.fDz ;
-  this->fDx = source.fDx ;
-  this->fMinNCells = source.fMinNCells ;
-  return;
+//______________________________________________________________________
+AliITSClusterFinderSPD::AliITSClusterFinderSPD(const AliITSClusterFinderSPD &source) : AliITSClusterFinder(source) {
+  // Copy constructor
+  // Copies are not allowed. The method is protected to avoid misuse.
+  Fatal("AliITSClusterFinderSPD","Copy constructor not allowed\n");
 }
 
-//_________________________________________________________________________
-AliITSClusterFinderSPD& 
-  AliITSClusterFinderSPD::operator=(const AliITSClusterFinderSPD &source) {
-  //    Assignment operator
-  if(&source == this) return *this;
-  this->fClusters = source.fClusters ;
-  this->fNclusters = source.fNclusters ;
-  this->fMap = source.fMap ;
-  this->fDz = source.fDz ;
-  this->fDx = source.fDx ;
-  this->fMinNCells = source.fMinNCells ;
+//______________________________________________________________________
+AliITSClusterFinderSPD& AliITSClusterFinderSPD::operator=(const AliITSClusterFinderSPD& /* source */){
+  // Assignment operator
+  // Assignment is not allowed. The method is protected to avoid misuse.
+  Fatal("= operator","Assignment operator not allowed\n");
   return *this;
 }
-
-//_____________________________________________________________________________
-void AliITSClusterFinderSPD::SetMap()
-{
-  // set map
-
-  if(!fMap) fMap=new AliITSMapA1(fSegmentation,fDigits);
-
+//______________________________________________________________________
+void AliITSClusterFinderSPD::FindRawClusters(Int_t module){   
+    // input of Cluster Finder
+    Int_t   digitcount  = 0;
+    Int_t   numberd     = 100000;
+    Int_t   *digx       = new Int_t[numberd];
+    Int_t   *digz       = new Int_t[numberd];
+    Int_t   *digtr1     = new Int_t[numberd];
+    Int_t   *digtr2     = new Int_t[numberd];
+    Int_t   *digtr3     = new Int_t[numberd];
+    Int_t   *digtr4     = new Int_t[numberd];
+    //  output of Cluster Finder    
+    Int_t   numberc     = 10000;
+    Double_t *xcenterl   = new Double_t[numberc];
+    Double_t *zcenterl   = new Double_t[numberc];
+    Double_t *errxcenter = new Double_t[numberc];
+    Double_t *errzcenter = new Double_t[numberc];
+    Int_t   *tr1clus    = new Int_t[numberc];
+    Int_t   *tr2clus    = new Int_t[numberc];
+    Int_t   *tr3clus    = new Int_t[numberc];
+    Int_t   nclus;
+
+    SetModule(module);
+    digitcount=0;
+    Int_t ndigits = Digits()->GetEntriesFast();  
+    if (!ndigits) return;
+
+    AliITSdigitSPD *dig;
+    Int_t ndig=0,i;
+    /*
+    AliDebug(4," ");
+        scanf("%d",&ndig);
+    */
+    for(ndig=0; ndig<ndigits; ndig++) {
+        dig= (AliITSdigitSPD*)GetDigit(ndig);
+        digx[digitcount] = dig->GetCoord2()+1;  //starts at 1
+        digz[digitcount] = dig->GetCoord1()+1;  //starts at 1
+        digtr1[digitcount] = dig->GetTrack(0);
+        digtr2[digitcount] = -3;
+        digtr3[digitcount] = -3;
+        AliDebug(5,Form("digtr1[%d]=%d fTracks[%d]=%d: ",
+			digitcount,digtr1[digitcount],0,dig->GetTrack(0)));
+        i=1;
+        while(digtr1[digitcount]==dig->GetTrack(i) && i<dig->GetNTracks()) i++;
+        AliDebug(5,Form(" fTracks[%d]=%d",i,dig->GetTrack(i)));
+        if(i<dig->GetNTracks()){
+            digtr2[digitcount] = dig->GetTrack(i);
+            AliDebug(5,Form("digtr2[%d]=%d: ",digitcount,digtr2[digitcount]));
+            while((digtr1[digitcount]==dig->GetTrack(i) || 
+                   digtr2[digitcount]==dig->GetTrack(i))&&
+                  i<=dig->GetNTracks()) i++;
+            if(i<dig->GetNTracks()) digtr3[digitcount] = dig->GetTrack(i);
+	    AliDebug(5,Form(" fTracks[%d]=%d digtr3[%d]=%d",
+			    i,i<dig->GetNTracks()?dig->GetTrack(i):-1,digitcount,digtr3[digitcount]));
+        } // end if
+    //        if(GetDebug(4)) cout<<endl;
+        digtr4[digitcount] = dig->GetSignal();
+        digitcount++;
+    } // end for ndig
+    ClusterFinder(digitcount,digx,digz,digtr1,digtr2,digtr3,digtr4,
+                  nclus,xcenterl,zcenterl,errxcenter,errzcenter,
+                  tr1clus, tr2clus, tr3clus);
+    DigitToPoint(nclus,xcenterl,zcenterl,errxcenter,errzcenter,
+                 tr1clus, tr2clus, tr3clus);
+    delete[] digx;
+    delete[] digz;
+    delete[] digtr1;
+    delete[] digtr2;
+    delete[] digtr3;
+    delete[] digtr4;
+    delete[] xcenterl;
+    delete[] zcenterl;
+    delete[] errxcenter;
+    delete[] errzcenter;
+    delete[] tr1clus;
+    delete[] tr2clus;
+    delete[] tr3clus;
 }
-
-//_____________________________________________________________________________
-
-void AliITSClusterFinderSPD::Find1DClusters()
-{
-  // Find one dimensional clusters, i.e.
-  // in r*phi(x) direction for each colunm in z direction
-  
-  AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
-  
-  // retrieve the parameters 
-  Int_t fNofPixels = fSegmentation->Npz(); 
-  Int_t fMaxNofSamples = fSegmentation->Npx();
-  
-  // read in digits -> do not apply threshold 
-  // signal in fired pixels is always 1
-
-  fMap->FillMap();
-  
-  Int_t nofFoundClusters = 0;
-  
-  Int_t k,it,m;
-  for(k=0;k<fNofPixels;k++) {
-    
-    Int_t mmax = 10;  // a size of the window for the cluster finding
+//----------------------------------------------------------------------
+void AliITSClusterFinderSPD::ClusterFinder(Int_t ndigits,Int_t digx[],
+					   Int_t digz[],Int_t digtr1[],
+					   Int_t digtr2[],Int_t digtr3[],
+					   Int_t digtr4[],Int_t &nclus,
+					   Double_t xcenter[],Double_t zcenter[],
+					   Double_t errxcenter[],
+					   Double_t errzcenter[],
+					   Int_t tr1clus[],Int_t tr2clus[],
+					   Int_t tr3clus[]){
+    // Search for clusters of fired pixels (digits). Two digits are linked
+    // inside a cluster if they are countiguous both in row or column
+    // direction.  Diagonal digits are not linked.
+    // xcenter, ycenter, zcenter are the coordinates of the center
+    // of each found cluster, calculated from the averaging the corresponding
+    // coordinate of the center of the linked digits. The coordinates are
+    // given in the local reference sistem. 
+    // errxcenter, errycenter, errzcenter are the errors associated to
+    // the corresponding average.
+    Int_t   if1, min, max, nd;
+    Int_t   x1, z1, t1, t2, t3, t4;
+    Int_t   ndx, ndz, ndxmin=0, ndxmax=0, ndzmin=0, ndzmax=0;
+    Double_t dx, dz; 
+    Int_t   i,k,ipos=0;
+    Float_t xdum, zdum;      
+    Int_t   kmax, sigmax;
+    Double_t deltax, deltaz;
+    Double_t ndig;
+    Int_t   numberd = 10000;
+    Int_t   *ifpad  = new Int_t[numberd];
+    Int_t   *xpad   = new Int_t[numberd];
+    Int_t   *zpad   = new Int_t[numberd];
+    Int_t   *tr1pad = new Int_t[numberd];
+    Int_t   *tr2pad = new Int_t[numberd];
+    Int_t   *tr3pad = new Int_t[numberd];
+    Int_t   *tr4pad = new Int_t[numberd];
+    Int_t   *iclus  = new Int_t[numberd];
+
+    nclus=1;
+    for (i=0; i < ndigits ; i++){
+        ifpad[i] = -1;
+        iclus[i] = 0;
+    } // end for i
+
+    ifpad[0]=0;
+    for (i=0; i < ndigits-1 ; i++) {
+        if ( ifpad[i] == -1 ) { 
+            nclus++;
+            ipos++;
+            ifpad[i]=nclus-1;
+        } // end if ipad[i]
+        for (Int_t j=i+1 ; j < ndigits ; j++)  {  
+            if (ifpad[j]== -1 ) {
+                dx = TMath::Abs(digx[i]-digx[j]);
+                dz = TMath::Abs(digz[i]-digz[j]);
+                // if ( ( dx+dz )==1 )  //clusters are not diagonal
+                if(( dx+dz )==1 || (dx==1 && dz==1)){
+                    //diagonal clusters allowed
+                    ipos++;
+                    ifpad[j] = ifpad[i];
+                    
+                    x1         = digx[j];
+                    z1         = digz[j];
+                    digx[j]    = digx[ipos];
+                    digz[j]    = digz[ipos];
+                    digx[ipos] = x1;
+                    digz[ipos] = z1;
+                    
+                    t1 = digtr1[j];
+                    t2 = digtr2[j];
+                    t3 = digtr3[j];
+                    t4 = digtr4[j];
+                    digtr1[j] = digtr1[ipos];
+                    digtr2[j] = digtr2[ipos];
+                    digtr3[j] = digtr3[ipos];
+                    digtr4[j] = digtr4[ipos];
+                    digtr1[ipos] = t1;
+                    digtr2[ipos] = t2;
+                    digtr3[ipos] = t3;
+                    digtr4[ipos] = t4;
+                    
+                    if1 = ifpad[j];
+                    ifpad[j] = ifpad[ipos];
+                    ifpad[ipos] = if1;
+                } // end dx+dx...
+            }// end if ifpad[j]== -1 
+        } // end for j
+    }//end loop on digits
+    if ( ifpad[ndigits-1] == -1 ) {
+        nclus++;
+        ifpad[ndigits-1]=nclus-1;
+    } // end if ifpad[ndigits-1] == -1
     
-    for(it=0;it<fMaxNofSamples;it++) {
-      
-      Int_t lclx = 0;
-      Int_t xstart = 0;
-      Int_t xstop = 0;
-      Int_t id = 0;
-      Int_t ilcl =0;
-      
-      for(m=0;m<mmax;m++) {  // find the cluster inside the window
-	id = it+m;
-	if(id >= fMaxNofSamples) break;    // ! no possible for the fadc 
-	
-	if(fMap->TestHit(k,id) == kUnused) {   // start of the cluster
-	  lclx += 1;
-	  if(lclx == 1) xstart = id;
-	  
-	}
-	
-	if(lclx > 0 && fMap->TestHit(k,id) == kEmpty) {  
-	  // end of cluster if a gap exists
-	  xstop = id-1;
-	  ilcl = 1;
-	  break;
-	}            
-	
-      }   //  end of m-loop
-      
-      if(lclx == 0 && ilcl == 0) it = id; // no cluster in the window,
-      // continue the "it" loop
-      
-      if(id >= fMaxNofSamples && lclx == 0) break; // the x row finished
-      
-      if(id < fMaxNofSamples && ilcl == 0 && lclx > 0) {  
-	                           // cluster end is outside of the window,
-	mmax += 5;                 // increase mmax and repeat the cluster
-	                           // finding
-	it -= 1;
-      }
-      
-      if(id >= fMaxNofSamples && lclx > 0) {  // the x row finished but
-	xstop = fMaxNofSamples - 1;           // the end cluster exists
-	ilcl = 1;
-      } 
-      
-      // ---  Calculate z and x coordinates for one dimensional clusters
-      
-      if(ilcl == 1) {         // new cluster exists
-	it = id;
-	mmax = 10;
-	    nofFoundClusters++;
-	    Float_t clusterCharge = 0.;
-            Float_t zpitch = fSegmentation->Dpz(k+1); 
-	    Float_t clusterZ, dummyX; 
-            Int_t dummy=0;
-            fSegmentation->GetPadCxz(dummy,k,dummyX,clusterZ);
-            Float_t zstart = clusterZ - 0.5*zpitch;
-            Float_t zstop = clusterZ + 0.5*zpitch;
-	    Float_t clusterX = 0.;
-            Int_t xstartfull = xstart;
-            Int_t xstopfull = xstop;
-	    Int_t clusterSizeX = lclx;
-	    Int_t clusterSizeZ = 1;
-	    
-            Int_t its;
-	    for(its=xstart; its<=xstop; its++) {
-              Int_t firedpixel=0;
-              if (fMap->GetHitIndex(k,its)>=0) firedpixel=1; 
-	      clusterCharge += firedpixel;
-              clusterX +=its + 0.5;
-	    }
-            Float_t fRphiPitch = fSegmentation->Dpx(dummy);
-	    clusterX /= (clusterSizeX/fRphiPitch); // center of gravity for x 
-	    
-	    // Write the points (coordinates and some cluster information) to the
-	    // AliITSRawClusterSPD object
-	    
-	    AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(clusterZ,clusterX,clusterCharge,clusterSizeZ,clusterSizeX,xstart,xstop,xstartfull,xstopfull,zstart,zstop,k);
-
-	    iTS->AddCluster(0,clust);
-	    
-      }    // new cluster (ilcl=1)
-    } // X direction loop (it)
-  } // Z direction loop (k)
-
-  //fMap->ClearMap();
-  return;
-  
+    for (i=0 ; i < ndigits ; i++) iclus[ifpad[i]]++;
+
+    min=0;
+    max=0;
+    // loop on found clusters 
+    for (i=0 ; i < nclus ; i++){
+        min = max;
+        max += iclus[i];
+        deltax = GetSeg()->Dpx(0);
+        if (iclus[i]!=1){
+            //cluster with more than one digit
+            nd=iclus[i];
+            ndig=(Double_t) nd;
+            Int_t count=0;
+            for (k=min;k<min+nd;k++){
+                xpad[count] = digx[k];	   
+                zpad[count] = digz[k];
+
+                tr1pad[count] = digtr1[k];	   
+                tr2pad[count] = digtr2[k];	   
+                tr3pad[count] = digtr3[k];	   
+                tr4pad[count] = digtr4[k];	   
+                
+                count++; 
+            } // end for k
+            ndxmin = xpad[TMath::LocMin(nd,xpad)];
+            ndxmax = xpad[TMath::LocMax(nd,xpad)];
+            ndzmin = zpad[TMath::LocMin(nd,zpad)];
+            ndzmax = zpad[TMath::LocMax(nd,zpad)];
+            ndx = ndxmax - ndxmin+1;
+            ndz = ndzmax - ndzmin+1;
+            
+            // calculate x and z coordinates of the center of the cluster
+            GetSeg()->GetPadCxz(digx[min],digz[min]-1,xdum, zdum);
+            
+            if (ndx == 1) {
+                xcenter[i] = xdum;
+            }else{ 
+                xcenter[i] = 0.;
+                for (k=0;k<nd;k++) {
+                    GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
+                    xcenter[i] += (xdum / nd);
+                } // end for k               
+            } // end if ndx
+            
+            if (ndz == 1) {
+                zcenter[i] = zdum;
+            } else {
+                zcenter[i] = 0.;
+                for (k=0;k<nd;k++) {	      
+                    GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
+                    zcenter[i] += (zdum / nd);
+                } // end for k
+            } // end if ndz
+            
+            // error on points in x and z directions
+            
+            if (ndx == 1) {
+                errxcenter[i] = deltax / TMath::Sqrt(12.);
+            } else {
+                errxcenter[i] = 0.;	 		
+                for (k=0;k<nd;k++){ 
+                    GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
+                    errxcenter[i] += ((xdum-xcenter[i])*(xdum-xcenter[i]))/
+                        (nd*(nd-1)); 
+                } // end for k
+                errxcenter[i] = TMath::Sqrt(errxcenter[i]);
+            } // end if ndx	
+            if (ndz == 1) {
+                deltaz = GetSeg()->Dpz(digz[min]);	              
+                errzcenter[i] = deltaz / TMath::Sqrt(12.);
+            } else {
+                errzcenter[i] = 0.;
+                for (k=0;k<nd;k++){ 
+                    GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
+                    errzcenter[i] += ((zdum-zcenter[i])*(zdum-zcenter[i]))/
+                        (nd*(nd-1));
+                } // end for k
+                errzcenter[i] = TMath::Sqrt(errzcenter[i]);
+            } // end if ndz
+            // take three track numbers for the cluster
+            // choose the track numbers of the digit with higher signal 
+            kmax = 0;
+            sigmax = 0;
+            for (k=0;k<nd;k++){
+                if(tr4pad[k] > sigmax){
+                    sigmax = tr4pad[k];
+                    kmax   = k;
+                } // end if tr4pad[k]
+            } // end for k
+            if(sigmax != 0) {
+                tr1clus[i]= tr1pad[kmax];
+                tr2clus[i]= tr2pad[kmax];
+                tr3clus[i]= tr3pad[kmax];
+            } else {
+                tr1clus[i]= -2;
+                tr2clus[i]= -2;
+                tr3clus[i]= -2;
+            } // end if sigmax
+        } else {
+            // cluster with single digit
+            ndig= 1.;
+            ndx = 1;
+            ndz = 1;
+            GetSeg()->GetPadCxz(digx[min],digz[min]-1,xdum,zdum);
+            xcenter[i] = xdum;
+            zcenter[i] = zdum;
+            tr1clus[i]=digtr1[min];
+            tr2clus[i]=digtr2[min];
+            tr3clus[i]=digtr3[min];
+            deltaz = GetSeg()->Dpz(digz[min]);
+            errxcenter[i] = deltax / TMath::Sqrt(12.);
+            errzcenter[i] = deltaz / TMath::Sqrt(12.);
+        } // end if iclus[i]
+
+        // store the cluster information to the AliITSRawCLusterSPD object
+   
+
+        //put the cluster center in local reference frame of the detector
+        // and in microns
+        xcenter[i] = xcenter[i] - GetSeg()->Dx()/2.; 
+        zcenter[i] = zcenter[i] - GetSeg()->Dz()/2.;
+
+        AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(zcenter[i], //f
+                                                             xcenter[i], //f
+                                                             ndig, //f
+                                                             ndz,ndx, //ii
+                                                             ndxmin,ndxmax,//ii
+                                                             (Double_t) ndzmin,
+                                                             (Double_t) ndzmax,
+                                                             0,GetModule());
+	fDetTypeRec->AddCluster(0,clust);
+        delete clust;
+    }//end loop on clusters   
+    delete[] ifpad;
+    delete[] xpad ;
+    delete[] zpad ;
+    delete[] iclus;
+    delete[] tr1pad;
+    delete[] tr2pad;
+    delete[] tr3pad;
+    delete[] tr4pad;
 }
-
-//_____________________________________________________________________________
-void  AliITSClusterFinderSPD::GroupClusters()
-{
-  // Find two dimensional clusters, i.e. group one dimensional clusters
-  // into two dimensional ones (go both in x and z directions).
-  
-  // get number of clusters for this module
-  Int_t nofClusters = fClusters->GetEntriesFast();
-  nofClusters -= fNclusters;
-
-  AliITSRawClusterSPD *clusterI;
-  AliITSRawClusterSPD *clusterJ;
-  
-  Int_t *label=new Int_t[nofClusters];  
-  Int_t i,j;
-  for(i=0; i<nofClusters; i++) label[i] = 0;
-  for(i=0; i<nofClusters; i++) {
-    if(label[i] != 0) continue;
-    for(j=i+1; j<nofClusters; j++) { 
-      if(label[j] != 0) continue;
-      clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
-      clusterJ = (AliITSRawClusterSPD*) fClusters->At(j);
-      Bool_t pair = clusterI->Brother(clusterJ,fDz,fDx);
-      if(pair) {     
-	
-	//    if((clusterI->XStop() == clusterJ->XStart()-1)||(clusterI->XStart()==clusterJ->XStop()+1)) cout<<"!! Diagonal cluster"<<endl;
-
-	/*
-	      cout << "clusters " << i << "," << j << " before grouping" << endl;
-	      clusterI->PrintInfo();
-	      clusterJ->PrintInfo();
-	*/
-
-	clusterI->Add(clusterJ);
-	//	cout << "remove cluster " << j << endl;
-	label[j] = 1;
-	fClusters->RemoveAt(j);
-
-	/*
-	  cout << "cluster  " << i << " after grouping" << endl;
-	  clusterI->PrintInfo();
-	*/
-
-      }  // pair
-    } // J clusters  
-    label[i] = 1;
-  } // I clusters
-  fClusters->Compress();
-
-  /*
-    Int_t totalNofClusters = fClusters->GetEntriesFast();
-    cout << " Nomber of clusters at the group end ="<< totalNofClusters<<endl;
-  */
-
-  delete [] label;
-
-  return;
-  
-  
+//______________________________________________________----------------
+void AliITSClusterFinderSPD::DigitToPoint(Int_t nclus,
+					  Double_t *xcenter,Double_t *zcenter,
+					  Double_t *errxcenter,
+					  Double_t *errzcenter, 
+					  Int_t *tr1clus, Int_t *tr2clus,
+					  Int_t *tr3clus){
+    // A point is associated to each cluster of SPD digits. The points
+    // and their associated errors are stored in the file galiceSP.root.
+    Double_t l[3],xg,zg;
+    const Double_t kconv = 1.0e-4; // micron -> cm
+
+    // get rec points
+    for (Int_t i=0; i<nclus; i++){
+        l[0] = kconv*xcenter[i];
+        l[1] = kconv*GetSeg()->Dy()/2.;
+        l[2] = kconv*zcenter[i];
+
+        xg = l[0]; 
+        zg = l[2]; 
+
+        Double_t sigma2x = (kconv*errxcenter[i]) * (kconv*errxcenter[i]);
+        Double_t sigma2z = (kconv*errzcenter[i]) * (kconv*errzcenter[i]);
+        AliITSRecPoint rnew;
+        rnew.SetX(xg);
+        rnew.SetZ(zg);
+        rnew.SetQ(1.);
+        rnew.SetdEdX(0.);
+        rnew.SetSigmaX2(sigma2x);
+        rnew.SetSigmaZ2(sigma2z);
+        rnew.fTracks[0]=tr1clus[i];
+        rnew.fTracks[1]=tr2clus[i];
+        rnew.fTracks[2]=tr3clus[i];
+	fDetTypeRec->AddRecPoint(rnew); 
+    } // end for i
 }
-//_____________________________________________________________________________
-
-void AliITSClusterFinderSPD::TracksInCluster()
-{
-  
-  // Find tracks creating one cluster
-
-  // get number of clusters for this module
-  Int_t nofClusters = fClusters->GetEntriesFast();
-  nofClusters -= fNclusters;
-
-  Int_t i, ix, iz, jx, jz, xstart, xstop, zstart, zstop, nclx, nclz;
-  const Int_t trmax = 100;
-  Int_t cltracks[trmax], itr, tracki, ii, is, js, ie, ntr, tr0, tr1, tr2;
-
-  for(i=0; i<nofClusters; i++) { 
-    ii = 0;
-    memset(cltracks,-1,sizeof(int)*trmax);
-    tr0=tr1=tr2=-1;
-
-    AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
-
-    nclx = clusterI->NclX();
-    nclz = clusterI->NclZ();
-    xstart = clusterI->XStartf();
-    xstop = clusterI->XStopf();
-    zstart = clusterI->Zend()-nclz+1;
-    zstop = clusterI->Zend();
-    Int_t ind; 
-
-     for(iz=0; iz<nclz; iz++) { 
-         jz = zstart + iz;
-       for(ix=0; ix<nclx; ix++) { 
-	 jx = xstart + ix;
-	 ind = fMap->GetHitIndex(jz,jx);
-	 if(ind < 0) {
-          continue;
-         }
-	 if(ind == 0 && iz >= 0 && ix > 0) {
-          continue;
-         }
-	 if(ind == 0 && iz > 0 && ix >= 0) {
-          continue;
-         }
-	 if(ind == 0 && iz == 0 && ix == 0 && i > 0) {
-          continue;
-         }
-
-        AliITSdigitSPD *dig = (AliITSdigitSPD*)fMap->GetHit(jz,jx);
-
-	/*
-         signal=dig->fSignal;
-         track0=dig->fTracks[0];
-         track1=dig->fTracks[1];
-         track2=dig->fTracks[2];
-	*/
-
-          for(itr=0; itr<3; itr++) { 
-	    tracki = dig->fTracks[itr];
-            if(tracki >= 0) {
-	      ii += 1;
-	      if(ii > 90) { 
-	      }
-	      if(ii < 99) cltracks[ii-1] = tracki;
-            }
-	  }
-       } // ix pixel
-     }  // iz pixel
- 
-     for(is=0; is<trmax; is++) { 
-         if(cltracks[is]<0) continue;
-       for(js=is+1; js<trmax; js++) { 
-         if(cltracks[js]<0) continue;
-         if(cltracks[js]==cltracks[is]) cltracks[js]=-5;
-       }
-     }
-
-     ntr = 0;
-     for(ie=0; ie<trmax; ie++) { 
-       if(cltracks[ie] >= 0) {
-        ntr=ntr+1;
-        if(ntr==1) tr0=cltracks[ie];
-        if(ntr==2) tr1=cltracks[ie];
-        if(ntr==3) tr2=cltracks[ie];
-       }
-     }
-     // if delta ray only
-     if(ntr == 0) ntr = 1;
-
-     clusterI->SetNTracks(ntr);
-     clusterI->SetTracks(tr0,tr1,tr2);
-
-  } // I cluster
-
-}
-//_____________________________________________________________________________
-
-void AliITSClusterFinderSPD::GetRecPoints()
-{
-  // get rec points
-  AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
-  
-  // get number of clusters for this module
-  Int_t nofClusters = fClusters->GetEntriesFast();
-  nofClusters -= fNclusters;
-  const Float_t kconv = 1.0e-4;
-  const Float_t kRMSx = 12.0*kconv; // microns -> cm ITS TDR Table 1.3
-  const Float_t kRMSz = 70.0*kconv; // microns -> cm ITS TDR Table 1.3
-
-  Float_t spdLength = fSegmentation->Dz();
-  Float_t spdWidth = fSegmentation->Dx();
-
-  Int_t i;
-  Int_t track0, track1, track2;
-
-  for(i=0; i<nofClusters; i++) { 
-
-    AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
-    clusterI->GetTracks(track0, track1, track2); 
-    AliITSRecPoint rnew;
-
-    rnew.SetX((clusterI->X() - spdWidth/2)*kconv);
-    rnew.SetZ((clusterI->Z() - spdLength/2)*kconv);
-    rnew.SetQ(1.);
-    rnew.SetdEdX(0.);
-    rnew.SetSigmaX2(kRMSx*kRMSx);
-    rnew.SetSigmaZ2(kRMSz*kRMSz);
-    rnew.fTracks[0]=track0;
-    rnew.fTracks[1]=track1;
-    rnew.fTracks[2]=track2;
-    iTS->AddRecPoint(rnew);
-  } // I clusters
-
-  fMap->ClearMap();
-  
-}
-//_____________________________________________________________________________
-
-void AliITSClusterFinderSPD::FindRawClusters(Int_t mod)
-{
-  // find raw clusters
-  Find1DClusters();
-  GroupClusters();
-  TracksInCluster();
-  GetRecPoints();
-
-}
-
-
-