[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSPDbari.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.                  *
 **************************************************************************/


#include "AliITSClusterFinderSPDbari.h"
#include "AliITS.h"
#include "AliITSgeom.h"
#include "AliITSdigit.h"
#include "AliITSRawCluster.h"
#include "AliITSRecPoint.h"
#include "AliITSsegmentation.h"
#include "AliITSresponse.h"
#include "AliRun.h"


ClassImp(AliITSClusterFinderSPDbari)

//----------------------------------------------------------
AliITSClusterFinderSPDbari::AliITSClusterFinderSPDbari
(AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp)   
{
  // constructor
    fSegmentation=seg;
    fDigits=digits;
    fClusters=recp;
    fNclusters= fClusters->GetEntriesFast();
    SetDx();
    SetDz();
}

//_____________________________________________________________________________
AliITSClusterFinderSPDbari::AliITSClusterFinderSPDbari()
{
  // constructor
  fSegmentation=0;
  fDigits=0;
  fClusters=0;
  fNclusters=0;
  SetDx();
  SetDz();
  
}

//__________________________________________________________________________
AliITSClusterFinderSPDbari::AliITSClusterFinderSPDbari(const
AliITSClusterFinderSPDbari &source){
  //     Copy Constructor 
  if(&source == this) return;
  this->fClusters = source.fClusters ;
  this->fNclusters = source.fNclusters ;
  this->fDz = source.fDz ;
  this->fDx = source.fDx ;
  return;
}

//_________________________________________________________________________
AliITSClusterFinderSPDbari& 
  AliITSClusterFinderSPDbari::operator=(const AliITSClusterFinderSPDbari &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 AliITSClusterFinderSPDbari::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;


    AliITSdigit *dig;
    AliITSdigitSPD *dig1;
    Int_t ndig;
    for(ndig=0; ndig<ndigits; ndig++) {
         dig= (AliITSdigit*)fDigits->UncheckedAt(ndig);
     
         digx[digitcount] = dig->fCoord2+1;  //starts at 1
         digz[digitcount] = dig->fCoord1+1;  //starts at 1

         dig1= (AliITSdigitSPD*)fDigits->UncheckedAt(ndig);

         digtr1[digitcount] = dig1->fTracks[0];
         digtr2[digitcount] = dig1->fTracks[1];
         digtr3[digitcount] = dig1->fTracks[2];
         digtr4[digitcount] = dig1->fSignal;

         digitcount++;
    }


    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 AliITSClusterFinderSPDbari::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[]) {
//
// 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, ndxmax, ndzmin, ndzmax;
  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;
  }

  ifpad[0]=0;
  for (i=0; i < ndigits-1 ; i++) {
    if ( ifpad[i] == -1 ) { 
	   	nclus++;
		ipos++;
     	ifpad[i]=nclus-1;
    }
    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 loop on digits   

  if ( ifpad[ndigits-1] == -1 ) {
	  nclus++;
	  ifpad[ndigits-1]=nclus-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++; 
        }
        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);
	         }	               
	    }

	    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);
	         }
	    }

        // 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)); 
             }   
	         errxcenter[i] = TMath::Sqrt(errxcenter[i]);
        }
	
	    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));
	        }
	        errzcenter[i] = TMath::Sqrt(errzcenter[i]);
	    }    
        // 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;
          }
        }
        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;
        }
     }
     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.);
     }

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