[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 "AliITSsegmentation.h"
#include "AliRun.h"

//#define DEBUG

ClassImp(AliITSClusterFinderSPD)

//----------------------------------------------------------
AliITSClusterFinderSPD::AliITSClusterFinderSPD(AliITSsegmentation *seg,
					       TClonesArray *digits,
					       TClonesArray *recp){
    // constructor

    fSegmentation = seg;
    fDigits       = digits;
    fClusters     = recp;
    fNclusters    = fClusters->GetEntriesFast();
    SetDx();
    SetDz();
}
//______________________________________________________________________
AliITSClusterFinderSPD::AliITSClusterFinderSPD(){
    // constructor

    fSegmentation = 0;
    fDigits       = 0;
    fClusters     = 0;
    fNclusters    = 0;
    SetDx();
    SetDz();
}
//_____________________________________________________________________
AliITSClusterFinderSPD::AliITSClusterFinderSPD(
         const AliITSClusterFinderSPD &source): AliITSClusterFinder(source){
    //     Copy Constructor 

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

    if(&source == this) return *this;
    this->fClusters  = source.fClusters ;
    this->fNclusters = source.fNclusters ;
    this->fDz        = source.fDz ;
    this->fDx        = source.fDx ;
    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;
    Float_t *xcenterl   = new Float_t[numberc];
    Float_t *zcenterl   = new Float_t[numberc];
    Float_t *errxcenter = new Float_t[numberc];
    Float_t *errzcenter = new Float_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;

    digitcount=0;
    Int_t ndigits = fDigits->GetEntriesFast();  
    if (!ndigits) return;

    AliITSdigitSPD *dig;
    Int_t ndig,i;
//    cout << "FindRawcluters"<<endl;
//    scanf("%d",ndig);
    for(ndig=0; ndig<ndigits; ndig++) {
	dig= (AliITSdigitSPD*)fDigits->UncheckedAt(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;
	//cout << "digtr1["<<digitcount <<"]="<<digtr1[digitcount];
	//cout << " fTracks["<<0<<"]="<<dig->fTracks[0]<<": ";
	i=1;
	while(digtr1[digitcount]==dig->GetTrack(i) && i<dig->GetNTracks()) i++;
	//cout << " fTracks["<<i<<"]="<<dig->fTracks[i];
	if(i<dig->GetNTracks()){
	    digtr2[digitcount] = dig->GetTrack(i);
	    //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);
	    //cout << " fTracks["<<i<<"]=";
	    //if(i<dig->GetNTracks()) cout <<dig->fTracks[i];
	    //cout << "digtr3["<<digitcount <<"]="<<digtr3[digitcount];
	} // end if
	//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,module);
    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,
					   Float_t xcenter[],Float_t zcenter[],
					   Float_t errxcenter[],
					   Float_t errzcenter[],
					   Int_t tr1clus[],Int_t tr2clus[],
					   Int_t tr3clus[],Int_t module){
    // 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;
    Float_t dx, dz; 
    Int_t   i,k,ipos=0;
    Float_t xdum, zdum;      
    Int_t   kmax, sigmax;
    Float_t deltax, deltaz;
    Float_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 = fSegmentation->Dpx(0);
	if (iclus[i]!=1){
	    //cluster with more than one digit
	    nd=iclus[i];
	    ndig=(Float_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
	    fSegmentation->GetPadCxz(digx[min],digz[min]-1,xdum, zdum);

	    if (ndx == 1) {
		xcenter[i] = xdum;
	    }else{ 
		xcenter[i] = 0.;
		for (k=0;k<nd;k++) {
		    fSegmentation->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++) {	      
		    fSegmentation->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++){ 
		    fSegmentation->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 = fSegmentation->Dpz(digz[min]);	              
	        errzcenter[i] = deltaz / TMath::Sqrt(12.);
	    } else {
	        errzcenter[i] = 0.;
	        for (k=0;k<nd;k++){ 
		    fSegmentation->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;
	    fSegmentation->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 = fSegmentation->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] - fSegmentation->Dx()/2.; 
	zcenter[i] = zcenter[i] - fSegmentation->Dz()/2.;

	AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(zcenter[i], //f
							     xcenter[i], //f
							     ndig, //f
							     ndz,ndx, //ii
							     ndxmin,ndxmax,//ii
							     (Float_t) ndzmin,
							     (Float_t) ndzmax,
							     0,module); //ii
	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,
					  Float_t *xcenter,Float_t *zcenter,
					  Float_t *errxcenter,
					  Float_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.
    Float_t l[3],xg,zg;
    const Float_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*fSegmentation->Dy()/2.;
        l[2] = kconv*zcenter[i];

        xg = l[0]; 
        zg = l[2]; 

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