[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSPD.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
//
//  Cluster finder
//  for Silicon pixels
//
//
#include "AliITSClusterFinderSPD.h"
#include "AliITS.h"
#include "AliITSdigitSPD.h"
#include "AliITSRawClusterSPD.h"
#include "AliITSRecPoint.h"
#include "AliITSresponseSPD.h"
#include "AliITSsegmentationSPD.h"
#include "AliRun.h"

//#define DEBUG

ClassImp(AliITSClusterFinderSPD)

//______________________________________________________________________
AliITSClusterFinderSPD::AliITSClusterFinderSPD():AliITSClusterFinder(),
fDz(0.0),
fDx(0.0),
fMinNCells(0){
    // constructor
}
//----------------------------------------------------------
AliITSClusterFinderSPD::AliITSClusterFinderSPD(AliITSsegmentation *seg,
                                               AliITSresponse *res):
AliITSClusterFinder(seg,res),
fDz(0.0),
fDx(0.0),
fMinNCells(0){
    // constructor

    SetDx();
    SetDz();
}
//----------------------------------------------------------
AliITSClusterFinderSPD::AliITSClusterFinderSPD(AliITSsegmentation *seg,
                                               TClonesArray *digits,
                                               TClonesArray *recp):
AliITSClusterFinder(seg,0),
fDz(0.0),
fDx(0.0),
fMinNCells(0){
    // constructor

    SetDigits(digits);
    SetClusters(recp);
    SetDx();
    SetDz();
}
//_____________________________________________________________________
AliITSClusterFinderSPD::AliITSClusterFinderSPD(
         const AliITSClusterFinderSPD &source): AliITSClusterFinder(source){
    //     Copy Constructor 

    if(&source == this) return;
    this->fDz        = source.fDz;
    this->fDx        = source.fDx;
    this->fMinNCells = source.fMinNCells;
    return;
}
//______________________________________________________________________
AliITSClusterFinderSPD& AliITSClusterFinderSPD::operator=(
                                       const AliITSClusterFinderSPD &source) {
    //    Assignment operator

    if(&source == this) return *this;
    this->fDz        = source.fDz;
    this->fDx        = source.fDx;
    this->fMinNCells = source.fMinNCells;
    return *this;
}
//______________________________________________________________________
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;
    /*
    if(GetDebug(4)){
        cout << "SPD - FindRawclusters"<<endl;
        scanf("%d",&ndig);
    } // end if GetDebug
    */
    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;
        if(GetDebug(5)){
            cout << "digtr1["<<digitcount <<"]="<<digtr1[digitcount];
            cout << " fTracks["<<0<<"]="<<dig->GetTrack(0)<<": ";
        } // end if GetDebug
        i=1;
        while(digtr1[digitcount]==dig->GetTrack(i) && i<dig->GetNTracks()) i++;
        if(GetDebug(5)) cout << " fTracks["<<i<<"]="<<dig->GetTrack(i);
        if(i<dig->GetNTracks()){
            digtr2[digitcount] = dig->GetTrack(i);
            if(GetDebug(5)) cout<<"digtr2["<<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);
            if(GetDebug(5)){
                cout << " fTracks["<<i<<"]=";
                if(i<dig->GetNTracks()) cout <<dig->GetTrack(i);
                cout << "digtr3["<<digitcount <<"]="<<digtr3[digitcount];
            } // end if GetDebug
        } // 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::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 (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
        static AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");

        //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());
        iTS->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::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
    static AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
    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];
        iTS->AddRecPoint(rnew); 
    } // end for i
}
Commit	Line	Data
c98c0281	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
bf3f2830	15	//
	16	// Cluster finder
	17	// for Silicon pixels
	18	//
	19	//
c98c0281	20	#include "AliITSClusterFinderSPD.h"
c98c0281	21	#include "AliITS.h"
e869281d	22	#include "AliITSdigitSPD.h"
41b19549	23	#include "AliITSRawClusterSPD.h"
c98c0281	24	#include "AliITSRecPoint.h"
aacedc3e	25	#include "AliITSresponseSPD.h"
aacedc3e	26	#include "AliITSsegmentationSPD.h"
c98c0281	27	#include "AliRun.h"
c98c0281	28
9e1e0cd7	29	//#define DEBUG
c98c0281	30
	31	ClassImp(AliITSClusterFinderSPD)
	32
aacedc3e	33	//______________________________________________________________________
	34	AliITSClusterFinderSPD::AliITSClusterFinderSPD():AliITSClusterFinder(),
	35	fDz(0.0),
	36	fDx(0.0),
	37	fMinNCells(0){
	38	// constructor
	39	}
c98c0281	40	//----------------------------------------------------------
9e1e0cd7	41	AliITSClusterFinderSPD::AliITSClusterFinderSPD(AliITSsegmentation *seg,
aacedc3e	42	AliITSresponse *res):
	43	AliITSClusterFinder(seg,res),
	44	fDz(0.0),
	45	fDx(0.0),
	46	fMinNCells(0){
9e1e0cd7	47	// constructor
9e1e0cd7	48
c98c0281	49	SetDx();
	50	SetDz();
	51	}
aacedc3e	52	//----------------------------------------------------------
	53	AliITSClusterFinderSPD::AliITSClusterFinderSPD(AliITSsegmentation *seg,
	54	TClonesArray *digits,
	55	TClonesArray *recp):
	56	AliITSClusterFinder(seg,0),
	57	fDz(0.0),
	58	fDx(0.0),
	59	fMinNCells(0){
9e1e0cd7	60	// constructor
9e1e0cd7	61
aacedc3e	62	SetDigits(digits);
aacedc3e	63	SetClusters(recp);
9e1e0cd7	64	SetDx();
9e1e0cd7	65	SetDz();
c98c0281	66	}
9e1e0cd7	67	//_____________________________________________________________________
ac74f489	68	AliITSClusterFinderSPD::AliITSClusterFinderSPD(
ac74f489	69	const AliITSClusterFinderSPD &source): AliITSClusterFinder(source){
9e1e0cd7	70	// Copy Constructor
	71
	72	if(&source == this) return;
aacedc3e	73	this->fDz = source.fDz;
	74	this->fDx = source.fDx;
	75	this->fMinNCells = source.fMinNCells;
9e1e0cd7	76	return;
c98c0281	77	}
9e1e0cd7	78	//______________________________________________________________________
	79	AliITSClusterFinderSPD& AliITSClusterFinderSPD::operator=(
	80	const AliITSClusterFinderSPD &source) {
	81	// Assignment operator
	82
	83	if(&source == this) return *this;
aacedc3e	84	this->fDz = source.fDz;
	85	this->fDx = source.fDx;
	86	this->fMinNCells = source.fMinNCells;
9e1e0cd7	87	return *this;
c98c0281	88	}
9e1e0cd7	89	//______________________________________________________________________
	90	void AliITSClusterFinderSPD::FindRawClusters(Int_t module){
	91	// input of Cluster Finder
	92	Int_t digitcount = 0;
	93	Int_t numberd = 100000;
c98c0281	94	Int_t *digx = new Int_t[numberd];
	95	Int_t *digz = new Int_t[numberd];
	96	Int_t *digtr1 = new Int_t[numberd];
	97	Int_t *digtr2 = new Int_t[numberd];
	98	Int_t *digtr3 = new Int_t[numberd];
	99	Int_t *digtr4 = new Int_t[numberd];
c98c0281	100	// output of Cluster Finder
9e1e0cd7	101	Int_t numberc = 10000;
aacedc3e	102	Double_t *xcenterl = new Double_t[numberc];
	103	Double_t *zcenterl = new Double_t[numberc];
	104	Double_t *errxcenter = new Double_t[numberc];
	105	Double_t *errzcenter = new Double_t[numberc];
c98c0281	106	Int_t *tr1clus = new Int_t[numberc];
	107	Int_t *tr2clus = new Int_t[numberc];
	108	Int_t *tr3clus = new Int_t[numberc];
9e1e0cd7	109	Int_t nclus;
c98c0281	110
aacedc3e	111	SetModule(module);
c98c0281	112	digitcount=0;
aacedc3e	113	Int_t ndigits = Digits()->GetEntriesFast();
c98c0281	114	if (!ndigits) return;
c98c0281	115
ee86d557	116	AliITSdigitSPD *dig;
aacedc3e	117	Int_t ndig=0,i;
5d766c75	118	/*
aacedc3e	119	if(GetDebug(4)){
5d766c75	120	cout << "SPD - FindRawclusters"<<endl;
fc9f1fd2	121	scanf("%d",&ndig);
aacedc3e	122	} // end if GetDebug
5d766c75	123	*/
c98c0281	124	for(ndig=0; ndig<ndigits; ndig++) {
aacedc3e	125	dig= (AliITSdigitSPD*)GetDigit(ndig);
	126	digx[digitcount] = dig->GetCoord2()+1; //starts at 1
	127	digz[digitcount] = dig->GetCoord1()+1; //starts at 1
	128	digtr1[digitcount] = dig->GetTrack(0);
	129	digtr2[digitcount] = -3;
	130	digtr3[digitcount] = -3;
	131	if(GetDebug(5)){
	132	cout << "digtr1["<<digitcount <<"]="<<digtr1[digitcount];
	133	cout << " fTracks["<<0<<"]="<<dig->GetTrack(0)<<": ";
	134	} // end if GetDebug
	135	i=1;
	136	while(digtr1[digitcount]==dig->GetTrack(i) && i<dig->GetNTracks()) i++;
	137	if(GetDebug(5)) cout << " fTracks["<<i<<"]="<<dig->GetTrack(i);
	138	if(i<dig->GetNTracks()){
	139	digtr2[digitcount] = dig->GetTrack(i);
	140	if(GetDebug(5)) cout<<"digtr2["<<digitcount <<"]="
	141	<<digtr2[digitcount]<<": ";
	142	while((digtr1[digitcount]==dig->GetTrack(i) \|\|
	143	digtr2[digitcount]==dig->GetTrack(i))&&
	144	i<=dig->GetNTracks()) i++;
	145	if(i<dig->GetNTracks()) digtr3[digitcount] = dig->GetTrack(i);
	146	if(GetDebug(5)){
	147	cout << " fTracks["<<i<<"]=";
	148	if(i<dig->GetNTracks()) cout <<dig->GetTrack(i);
	149	cout << "digtr3["<<digitcount <<"]="<<digtr3[digitcount];
	150	} // end if GetDebug
	151	} // end if
	152	if(GetDebug(4)) cout<<endl;
	153	digtr4[digitcount] = dig->GetSignal();
	154	digitcount++;
9e1e0cd7	155	} // end for ndig
c98c0281	156	ClusterFinder(digitcount,digx,digz,digtr1,digtr2,digtr3,digtr4,
aacedc3e	157	nclus,xcenterl,zcenterl,errxcenter,errzcenter,
aacedc3e	158	tr1clus, tr2clus, tr3clus);
c98c0281	159	DigitToPoint(nclus,xcenterl,zcenterl,errxcenter,errzcenter,
aacedc3e	160	tr1clus, tr2clus, tr3clus);
9e1e0cd7	161	delete[] digx;
	162	delete[] digz;
	163	delete[] digtr1;
	164	delete[] digtr2;
	165	delete[] digtr3;
	166	delete[] digtr4;
	167	delete[] xcenterl;
	168	delete[] zcenterl;
	169	delete[] errxcenter;
	170	delete[] errzcenter;
	171	delete[] tr1clus;
	172	delete[] tr2clus;
	173	delete[] tr3clus;
c98c0281	174	}
9e1e0cd7	175	//----------------------------------------------------------------------
	176	void AliITSClusterFinderSPD::ClusterFinder(Int_t ndigits,Int_t digx[],
	177	Int_t digz[],Int_t digtr1[],
	178	Int_t digtr2[],Int_t digtr3[],
	179	Int_t digtr4[],Int_t &nclus,
aacedc3e	180	Double_t xcenter[],Double_t zcenter[],
	181	Double_t errxcenter[],
	182	Double_t errzcenter[],
9e1e0cd7	183	Int_t tr1clus[],Int_t tr2clus[],
aacedc3e	184	Int_t tr3clus[]){
9e1e0cd7	185	// Search for clusters of fired pixels (digits). Two digits are linked
	186	// inside a cluster if they are countiguous both in row or column
	187	// direction. Diagonal digits are not linked.
	188	// xcenter, ycenter, zcenter are the coordinates of the center
	189	// of each found cluster, calculated from the averaging the corresponding
	190	// coordinate of the center of the linked digits. The coordinates are
	191	// given in the local reference sistem.
	192	// errxcenter, errycenter, errzcenter are the errors associated to
	193	// the corresponding average.
	194	Int_t if1, min, max, nd;
	195	Int_t x1, z1, t1, t2, t3, t4;
	196	Int_t ndx, ndz, ndxmin=0, ndxmax=0, ndzmin=0, ndzmax=0;
aacedc3e	197	Double_t dx, dz;
9e1e0cd7	198	Int_t i,k,ipos=0;
	199	Float_t xdum, zdum;
	200	Int_t kmax, sigmax;
aacedc3e	201	Double_t deltax, deltaz;
aacedc3e	202	Double_t ndig;
9e1e0cd7	203	Int_t numberd = 10000;
	204	Int_t *ifpad = new Int_t[numberd];
	205	Int_t *xpad = new Int_t[numberd];
	206	Int_t *zpad = new Int_t[numberd];
	207	Int_t *tr1pad = new Int_t[numberd];
	208	Int_t *tr2pad = new Int_t[numberd];
	209	Int_t *tr3pad = new Int_t[numberd];
	210	Int_t *tr4pad = new Int_t[numberd];
	211	Int_t *iclus = new Int_t[numberd];
	212
	213	nclus=1;
	214	for (i=0; i < ndigits ; i++){
aacedc3e	215	ifpad[i] = -1;
aacedc3e	216	iclus[i] = 0;
9e1e0cd7	217	} // end for i
	218
	219	ifpad[0]=0;
	220	for (i=0; i < ndigits-1 ; i++) {
aacedc3e	221	if ( ifpad[i] == -1 ) {
	222	nclus++;
	223	ipos++;
	224	ifpad[i]=nclus-1;
	225	} // end if ipad[i]
	226	for (Int_t j=i+1 ; j < ndigits ; j++) {
	227	if (ifpad[j]== -1 ) {
	228	dx = TMath::Abs(digx[i]-digx[j]);
	229	dz = TMath::Abs(digz[i]-digz[j]);
	230	// if ( ( dx+dz )==1 ) //clusters are not diagonal
	231	if(( dx+dz )==1 \|\| (dx==1 && dz==1)){
	232	//diagonal clusters allowed
	233	ipos++;
	234	ifpad[j] = ifpad[i];
	235
	236	x1 = digx[j];
	237	z1 = digz[j];
	238	digx[j] = digx[ipos];
	239	digz[j] = digz[ipos];
	240	digx[ipos] = x1;
	241	digz[ipos] = z1;
	242
	243	t1 = digtr1[j];
	244	t2 = digtr2[j];
	245	t3 = digtr3[j];
	246	t4 = digtr4[j];
	247	digtr1[j] = digtr1[ipos];
	248	digtr2[j] = digtr2[ipos];
	249	digtr3[j] = digtr3[ipos];
	250	digtr4[j] = digtr4[ipos];
	251	digtr1[ipos] = t1;
	252	digtr2[ipos] = t2;
	253	digtr3[ipos] = t3;
	254	digtr4[ipos] = t4;
	255
	256	if1 = ifpad[j];
	257	ifpad[j] = ifpad[ipos];
	258	ifpad[ipos] = if1;
	259	} // end dx+dx...
	260	}// end if ifpad[j]== -1
	261	} // end for j
9e1e0cd7	262	}//end loop on digits
9e1e0cd7	263	if ( ifpad[ndigits-1] == -1 ) {
aacedc3e	264	nclus++;
aacedc3e	265	ifpad[ndigits-1]=nclus-1;
9e1e0cd7	266	} // end if ifpad[ndigits-1] == -1
aacedc3e	267
9e1e0cd7	268	for (i=0 ; i < ndigits ; i++) iclus[ifpad[i]]++;
	269
	270	min=0;
	271	max=0;
	272	// loop on found clusters
	273	for (i=0 ; i < nclus ; i++){
aacedc3e	274	min = max;
	275	max += iclus[i];
	276	deltax = GetSeg()->Dpx(0);
	277	if (iclus[i]!=1){
	278	//cluster with more than one digit
	279	nd=iclus[i];
	280	ndig=(Double_t) nd;
	281	Int_t count=0;
	282	for (k=min;k<min+nd;k++){
	283	xpad[count] = digx[k];
	284	zpad[count] = digz[k];
9e1e0cd7	285
aacedc3e	286	tr1pad[count] = digtr1[k];
	287	tr2pad[count] = digtr2[k];
	288	tr3pad[count] = digtr3[k];
	289	tr4pad[count] = digtr4[k];
	290
	291	count++;
	292	} // end for k
	293	ndxmin = xpad[TMath::LocMin(nd,xpad)];
	294	ndxmax = xpad[TMath::LocMax(nd,xpad)];
	295	ndzmin = zpad[TMath::LocMin(nd,zpad)];
	296	ndzmax = zpad[TMath::LocMax(nd,zpad)];
	297	ndx = ndxmax - ndxmin+1;
	298	ndz = ndzmax - ndzmin+1;
	299
	300	// calculate x and z coordinates of the center of the cluster
	301	GetSeg()->GetPadCxz(digx[min],digz[min]-1,xdum, zdum);
	302
	303	if (ndx == 1) {
	304	xcenter[i] = xdum;
	305	}else{
	306	xcenter[i] = 0.;
	307	for (k=0;k<nd;k++) {
	308	GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
	309	xcenter[i] += (xdum / nd);
	310	} // end for k
	311	} // end if ndx
	312
	313	if (ndz == 1) {
	314	zcenter[i] = zdum;
	315	} else {
	316	zcenter[i] = 0.;
	317	for (k=0;k<nd;k++) {
	318	GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
	319	zcenter[i] += (zdum / nd);
	320	} // end for k
	321	} // end if ndz
	322
	323	// error on points in x and z directions
	324
	325	if (ndx == 1) {
	326	errxcenter[i] = deltax / TMath::Sqrt(12.);
	327	} else {
	328	errxcenter[i] = 0.;
	329	for (k=0;k<nd;k++){
	330	GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
	331	errxcenter[i] += ((xdum-xcenter[i])*(xdum-xcenter[i]))/
	332	(nd*(nd-1));
	333	} // end for k
	334	errxcenter[i] = TMath::Sqrt(errxcenter[i]);
	335	} // end if ndx
	336	if (ndz == 1) {
	337	deltaz = GetSeg()->Dpz(digz[min]);
	338	errzcenter[i] = deltaz / TMath::Sqrt(12.);
	339	} else {
	340	errzcenter[i] = 0.;
	341	for (k=0;k<nd;k++){
	342	GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
	343	errzcenter[i] += ((zdum-zcenter[i])*(zdum-zcenter[i]))/
	344	(nd*(nd-1));
	345	} // end for k
	346	errzcenter[i] = TMath::Sqrt(errzcenter[i]);
	347	} // end if ndz
	348	// take three track numbers for the cluster
	349	// choose the track numbers of the digit with higher signal
350	kmax = 0;
351	sigmax = 0;
352	for (k=0;k<nd;k++){
353	if(tr4pad[k] > sigmax){
354	sigmax = tr4pad[k];
355	kmax = k;
356	} // end if tr4pad[k]
357	} // end for k
358	if(sigmax != 0) {
359	tr1clus[i]= tr1pad[kmax];
360	tr2clus[i]= tr2pad[kmax];
361	tr3clus[i]= tr3pad[kmax];
362	} else {
363	tr1clus[i]= -2;
364	tr2clus[i]= -2;
365	tr3clus[i]= -2;
366	} // end if sigmax
367	} else {
368	// cluster with single digit
369	ndig= 1.;
370	ndx = 1;
371	ndz = 1;
372	GetSeg()->GetPadCxz(digx[min],digz[min]-1,xdum,zdum);
373	xcenter[i] = xdum;
374	zcenter[i] = zdum;
375	tr1clus[i]=digtr1[min];
376	tr2clus[i]=digtr2[min];
377	tr3clus[i]=digtr3[min];
378	deltaz = GetSeg()->Dpz(digz[min]);
379	errxcenter[i] = deltax / TMath::Sqrt(12.);
380	errzcenter[i] = deltaz / TMath::Sqrt(12.);
381	} // end if iclus[i]
9e1e0cd7	382
aacedc3e	383	// store the cluster information to the AliITSRawCLusterSPD object
aacedc3e	384	static AliITS iTS=(AliITS)gAlice->GetModule("ITS");
9e1e0cd7	385
aacedc3e	386	//put the cluster center in local reference frame of the detector
	387	// and in microns
	388	xcenter[i] = xcenter[i] - GetSeg()->Dx()/2.;
	389	zcenter[i] = zcenter[i] - GetSeg()->Dz()/2.;
9e1e0cd7	390
aacedc3e	391	AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(zcenter[i], //f
	392	xcenter[i], //f
	393	ndig, //f
	394	ndz,ndx, //ii
	395	ndxmin,ndxmax,//ii
	396	(Double_t) ndzmin,
	397	(Double_t) ndzmax,
	398	0,GetModule());
	399	iTS->AddCluster(0,clust);
	400	delete clust;
9e1e0cd7	401	}//end loop on clusters
	402	delete[] ifpad;
	403	delete[] xpad ;
	404	delete[] zpad ;
	405	delete[] iclus;
	406	delete[] tr1pad;
	407	delete[] tr2pad;
	408	delete[] tr3pad;
	409	delete[] tr4pad;
c98c0281	410	}
9e1e0cd7	411	//______________________________________________________----------------
c98c0281	412	void AliITSClusterFinderSPD::DigitToPoint(Int_t nclus,
aacedc3e	413	Double_t xcenter,Double_t zcenter,
	414	Double_t *errxcenter,
	415	Double_t *errzcenter,
9e1e0cd7	416	Int_t tr1clus, Int_t tr2clus,
	417	Int_t *tr3clus){
	418	// A point is associated to each cluster of SPD digits. The points
	419	// and their associated errors are stored in the file galiceSP.root.
aacedc3e	420	Double_t l[3],xg,zg;
aacedc3e	421	const Double_t kconv = 1.0e-4; // micron -> cm
9e1e0cd7	422
	423	// get rec points
	424	static AliITS iTS=(AliITS)gAlice->GetModule("ITS");
	425	for (Int_t i=0; i<nclus; i++){
c98c0281	426	l[0] = kconv*xcenter[i];
aacedc3e	427	l[1] = kconv*GetSeg()->Dy()/2.;
c98c0281	428	l[2] = kconv*zcenter[i];
	429
	430	xg = l[0];
	431	zg = l[2];
	432
aacedc3e	433	Double_t sigma2x = (kconverrxcenter[i]) (kconv*errxcenter[i]);
aacedc3e	434	Double_t sigma2z = (kconverrzcenter[i]) (kconv*errzcenter[i]);
c98c0281	435	AliITSRecPoint rnew;
	436	rnew.SetX(xg);
	437	rnew.SetZ(zg);
	438	rnew.SetQ(1.);
	439	rnew.SetdEdX(0.);
	440	rnew.SetSigmaX2(sigma2x);
	441	rnew.SetSigmaZ2(sigma2z);
	442	rnew.fTracks[0]=tr1clus[i];
	443	rnew.fTracks[1]=tr2clus[i];
	444	rnew.fTracks[2]=tr3clus[i];
	445	iTS->AddRecPoint(rnew);
9e1e0cd7	446	} // end for i
c98c0281	447	}