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

/*
$Log$
Revision 1.6  2002/10/22 14:45:33  alibrary
Introducing Riostream.h

Revision 1.5  2002/10/14 14:57:00  hristov
Merging the VirtualMC branch to the main development branch (HEAD)

Revision 1.2.8.1  2002/08/28 15:06:50  alibrary
Updating to v3-09-01

Revision 1.4  2002/08/22 13:18:46  nilsen
Added = operator and copy constructor in an attempt to fix complilation
problems found by Yves Schutz.

Revision 1.3  2002/08/21 22:04:27  nilsen
Added data member to SPD cluters and made related modifications to the
SPD Cluster finders. Generaly cleanded up some of the code.

Revision 1.2  2001/06/14 14:33:53  barbera
New version from B. Batyunya to get the Dubna model work with the present HEAD

Revision 1.15  2001/05/01 22:37:44  nilsen
Added a dummy argument to FindRawClusters. Argument used for SSD version.

Revision 1.14  2001/03/05 14:48:46  nilsen
Fixed a reoccuring bug. Array sizes must be declare const.

*/
//  Cluster Finder for
//  Silicon
//  Pixel
//  developed by the Dubna group
//
#include <Riostream.h>
#include "AliITSClusterFinderSPDdubna.h"
#include "AliITSMapA1.h"
#include "AliITS.h"
#include "AliITSdigit.h"
#include "AliITSRawCluster.h"
#include "AliITSRecPoint.h"
#include "AliITSsegmentation.h"
#include "AliRun.h"

//#define DEBUG

ClassImp(AliITSClusterFinderSPDdubna)

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

    fSegmentation = seg;
    fDigits       = digits;
    fClusters     = recp;
    fNclusters    = 0;
    fMap          = new AliITSMapA1(fSegmentation,fDigits);
    SetDx();
    SetDz();
    SetNCells();
}
//_______________________________________________________________________
AliITSClusterFinderSPDdubna::AliITSClusterFinderSPDdubna(){
    // constructor

    fSegmentation = 0;
    fDigits       = 0;
    fClusters     = 0;
    fNclusters    = 0;
    fMap          = 0;
    SetDx();
    SetDz();
    SetNCells();
}
//______________________________________________________________________
AliITSClusterFinderSPDdubna &AliITSClusterFinderSPDdubna::operator=(const
                                          AliITSClusterFinderSPDdubna &s){
    // The = operator for the class AliITSClusterFinderSPDdugna
    // Inputs:
    //    AliITSClusterFinderSPDdubna @s  The AliITSClusterFinderSPDdubna
    //                                    class to be copy constructed.
    // Outputs:
    //    none.
    // Returned:
    //    this, a copy of the class s.

    this->fNclusters = s.fNclusters;
    this->fDz        = s.fDz;
    this->fDx        = s.fDx;
    this->fMinNCells = s.fMinNCells;
    this->fClusters  = new TClonesArray(*(s.fClusters));
    return *this;
}
//______________________________________________________________________
AliITSClusterFinderSPDdubna::AliITSClusterFinderSPDdubna(const
                                         AliITSClusterFinderSPDdubna &s){
    // The Copy constructortor the class AliITSClusterFinderSPDdugna
    // It calles the = operator.
    // Inputs:
    //    AliITSClusterFinderSPDdubna @s  The AliITSClusterFinderSPDdubna
    //                                    class to be copy constructed.
    // Outputs:
    //    none.
    // Returned:
    //    none.

    *this = s;
    return;
}
//______________________________________________________________________
AliITSClusterFinderSPDdubna::~AliITSClusterFinderSPDdubna(){
    // destructor

    if (fMap) delete fMap;
}
//______________________________________________________________________
void AliITSClusterFinderSPDdubna::Find1DClusters(Int_t mod){
    // Find one dimensional clusters, i.e.
    // in r*phi(x) direction for each colunm in z direction

    static 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
	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;
		} // end if fMap
		if(lclx > 0 && fMap->TestHit(k,id) == kEmpty) {  
		    // end of cluster if a gap exists
		    xstop = id-1;
		    ilcl = 1;
		    break;
		} // end if lclx >0 ...
	    }   //  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;
	    } // end if id<
	    if(id >= fMaxNofSamples && lclx > 0) {  // the x row finished but
		xstop = fMaxNofSamples - 1;           // the end cluster exists
		ilcl = 1;
	    }  // end if id>=      
	    // ---  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);
		Int_t   dummy         = 0;
		Float_t clusterZ, dummyX; 
		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   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;
		} // end for its
		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(clusterZ,clusterX,clusterCharge,
					  clusterSizeZ,clusterSizeX,
					  xstart,xstop,zstart,zstop,k,mod);
		iTS->AddCluster(0,&clust);
	    }    // new cluster (ilcl=1)
	} // X direction loop (it)
    } // Z direction loop (k)
    //fMap->ClearMap();
    return;
}
//______________________________________________________________________
void  AliITSClusterFinderSPDdubna::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) {
#ifndef DEBUG
		cout << "clusters " << i << "," << j << " before grouping" << endl;
		clusterI->PrintInfo();
		clusterJ->PrintInfo();
#endif
		clusterI->Add(clusterJ);
#ifndef DEBUG
		cout << "remove cluster " << j << endl;
#endif
		label[j] = 1;
		fClusters->RemoveAt(j);
#ifndef DEBUG
		cout << "cluster  " << i << " after grouping" << endl;
		clusterI->PrintInfo();
#endif
	    }  // pair
	} // J clusters  
	label[i] = 1;
    } // I clusters
    fClusters->Compress();
#ifndef DEBUG
    Int_t totalNofClusters = fClusters->GetEntriesFast();
    cout << " Nomber of clusters at the group end ="<< totalNofClusters<<endl;
#endif
    delete [] label;
    return;
}
//______________________________________________________________________
void AliITSClusterFinderSPDdubna::TracksInCluster(){
    // Find tracks creating one cluster
    // get number of clusters for this module
    Int_t nofClusters = fClusters->GetEntriesFast();
    Int_t i, ix, iz, jx, jz, xstart, xstop, zstart, zstop, nclx, nclz;
    const Int_t kTrmax = 100;
    Int_t cltracks[kTrmax], itr, tracki, ii, is, js, ie, ntr, tr0, tr1, tr2;

    nofClusters -= fNclusters;
    for(i=0; i<nofClusters; i++) { 
	ii = 0;
	memset(cltracks,-1,sizeof(int)*kTrmax);
	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);
		for(itr=0; itr<3; itr++) { 
		    tracki = dig->fTracks[itr];
		    if(tracki >= 0) {
			ii += 1;
			if(ii < 99) cltracks[ii-1] = tracki;
		    } // end if tracki >-0
		} // end for itr
	    } // ix pixel
	}  // iz pixel
	for(is=0; is<kTrmax; is++) { 
	    if(cltracks[is]<0) continue;
	    for(js=is+1; js<kTrmax; js++) { 
		if(cltracks[js]<0) continue;
		if(cltracks[js]==cltracks[is]) cltracks[js]=-5;
	    } // end for js
	} // end for is
	ntr = 0;
	for(ie=0; ie<kTrmax; 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];
	    } // end if cltracks[ie] >-=
	} // end for ie
	// if delta ray only
	if(ntr == 0) ntr = 1;
	clusterI->SetNTracks(ntr);
	clusterI->SetTracks(tr0,tr1,tr2);
    } // I cluster
}
//______________________________________________________________________
void AliITSClusterFinderSPDdubna::GetRecPoints(){
    // get rec points
    static AliITS *iTS = (AliITS*)gAlice->GetModule("ITS");
    // get number of clusters for this module
    Int_t nofClusters   = fClusters->GetEntriesFast();
    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   dummy       = 0;
    Float_t xpitch      = fSegmentation->Dpx(dummy);
    Int_t i;
    Int_t track0, track1, track2;

    nofClusters -= fNclusters;
    for(i=0; i<nofClusters; i++) {
	AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
	Int_t clustersizex = clusterI->NclX();
	Int_t clustersizez = clusterI->NclZ();
	Int_t xstart = clusterI->XStartf();
	Float_t clusterx = 0;
	Float_t clusterz = clusterI->Z();
	for(Int_t ii=0; ii<clustersizex; ii++) clusterx += (xstart+0.5+ii)*xpitch;
	clusterx /= clustersizex;
	clusterz /= clustersizez;   
	clusterI->GetTracks(track0, track1, track2); 
	AliITSRecPoint rnew;
	rnew.SetX((clusterx - spdWidth/2)*kconv);
	rnew.SetZ((clusterz - 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 AliITSClusterFinderSPDdubna::FindRawClusters(Int_t mod){
    // find raw clusters

#ifndef DEBUG
    cout<<"SPDdubna: module ="<<mod<<endl;
#endif
    Find1DClusters(mod);
    GroupClusters();
    TracksInCluster();
    GetRecPoints();
}
Commit	Line	Data
409f8c84	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	$Log$
bf3f2830	18	Revision 1.6 2002/10/22 14:45:33 alibrary
	19	Introducing Riostream.h
	20
4ae5bbc4	21	Revision 1.5 2002/10/14 14:57:00 hristov
	22	Merging the VirtualMC branch to the main development branch (HEAD)
	23
b9d0a01d	24	Revision 1.2.8.1 2002/08/28 15:06:50 alibrary
	25	Updating to v3-09-01
	26
	27	Revision 1.4 2002/08/22 13:18:46 nilsen
	28	Added = operator and copy constructor in an attempt to fix complilation
	29	problems found by Yves Schutz.
	30
579fdb1c	31	Revision 1.3 2002/08/21 22:04:27 nilsen
	32	Added data member to SPD cluters and made related modifications to the
	33	SPD Cluster finders. Generaly cleanded up some of the code.
	34
9e1e0cd7	35	Revision 1.2 2001/06/14 14:33:53 barbera
	36	New version from B. Batyunya to get the Dubna model work with the present HEAD
	37
409f8c84	38	Revision 1.15 2001/05/01 22:37:44 nilsen
	39	Added a dummy argument to FindRawClusters. Argument used for SSD version.
	40
	41	Revision 1.14 2001/03/05 14:48:46 nilsen
	42	Fixed a reoccuring bug. Array sizes must be declare const.
	43
	44	*/
bf3f2830	45	// Cluster Finder for
	46	// Silicon
	47	// Pixel
	48	// developed by the Dubna group
	49	//
4ae5bbc4	50	#include <Riostream.h>
409f8c84	51	#include "AliITSClusterFinderSPDdubna.h"
	52	#include "AliITSMapA1.h"
	53	#include "AliITS.h"
	54	#include "AliITSdigit.h"
	55	#include "AliITSRawCluster.h"
	56	#include "AliITSRecPoint.h"
	57	#include "AliITSsegmentation.h"
409f8c84	58	#include "AliRun.h"
409f8c84	59
9e1e0cd7	60	//#define DEBUG
409f8c84	61
	62	ClassImp(AliITSClusterFinderSPDdubna)
	63
9e1e0cd7	64	//----------------------------------------------------------------------
	65	AliITSClusterFinderSPDdubna::AliITSClusterFinderSPDdubna(
	66	AliITSsegmentation seg, TClonesArray digits, TClonesArray *recp){
	67	// constructor
	68
	69	fSegmentation = seg;
	70	fDigits = digits;
	71	fClusters = recp;
	72	fNclusters = 0;
	73	fMap = new AliITSMapA1(fSegmentation,fDigits);
409f8c84	74	SetDx();
	75	SetDz();
	76	SetNCells();
	77	}
9e1e0cd7	78	//_______________________________________________________________________
	79	AliITSClusterFinderSPDdubna::AliITSClusterFinderSPDdubna(){
	80	// constructor
	81
	82	fSegmentation = 0;
	83	fDigits = 0;
	84	fClusters = 0;
	85	fNclusters = 0;
	86	fMap = 0;
	87	SetDx();
	88	SetDz();
	89	SetNCells();
409f8c84	90	}
9e1e0cd7	91	//______________________________________________________________________
579fdb1c	92	AliITSClusterFinderSPDdubna &AliITSClusterFinderSPDdubna::operator=(const
	93	AliITSClusterFinderSPDdubna &s){
	94	// The = operator for the class AliITSClusterFinderSPDdugna
	95	// Inputs:
	96	// AliITSClusterFinderSPDdubna @s The AliITSClusterFinderSPDdubna
	97	// class to be copy constructed.
	98	// Outputs:
	99	// none.
	100	// Returned:
	101	// this, a copy of the class s.
	102
	103	this->fNclusters = s.fNclusters;
	104	this->fDz = s.fDz;
	105	this->fDx = s.fDx;
	106	this->fMinNCells = s.fMinNCells;
	107	this->fClusters = new TClonesArray(*(s.fClusters));
	108	return *this;
	109	}
	110	//______________________________________________________________________
	111	AliITSClusterFinderSPDdubna::AliITSClusterFinderSPDdubna(const
	112	AliITSClusterFinderSPDdubna &s){
	113	// The Copy constructortor the class AliITSClusterFinderSPDdugna
	114	// It calles the = operator.
	115	// Inputs:
	116	// AliITSClusterFinderSPDdubna @s The AliITSClusterFinderSPDdubna
	117	// class to be copy constructed.
	118	// Outputs:
	119	// none.
	120	// Returned:
	121	// none.
	122
	123	*this = s;
	124	return;
	125	}
	126	//______________________________________________________________________
9e1e0cd7	127	AliITSClusterFinderSPDdubna::~AliITSClusterFinderSPDdubna(){
9e1e0cd7	128	// destructor
409f8c84	129
9e1e0cd7	130	if (fMap) delete fMap;
409f8c84	131	}
9e1e0cd7	132	//______________________________________________________________________
	133	void AliITSClusterFinderSPDdubna::Find1DClusters(Int_t mod){
	134	// Find one dimensional clusters, i.e.
	135	// in r*phi(x) direction for each colunm in z direction
338e4f06	136
9e1e0cd7	137	static AliITS iTS=(AliITS)gAlice->GetModule("ITS");
409f8c84	138
9e1e0cd7	139	// retrieve the parameters
	140	Int_t fNofPixels = fSegmentation->Npz();
	141	Int_t fMaxNofSamples = fSegmentation->Npx();
	142
	143	// read in digits -> do not apply threshold
	144	// signal in fired pixels is always 1
	145	fMap->FillMap();
	146
	147	Int_t nofFoundClusters = 0;
	148
	149	Int_t k,it,m;
	150	for(k=0;k<fNofPixels;k++) {
	151	Int_t mmax = 10; // a size of the window for the cluster finding
	152	for(it=0;it<fMaxNofSamples;it++) {
	153	Int_t lclx = 0;
	154	Int_t xstart = 0;
	155	Int_t xstop = 0;
	156	Int_t id = 0;
	157	Int_t ilcl = 0;
	158	for(m=0;m<mmax;m++) { // find the cluster inside the window
	159	id = it+m;
	160	if(id >= fMaxNofSamples) break; // ! no possible for the fadc
	161	if(fMap->TestHit(k,id) == kUnused) { // start of the cluster
	162	lclx += 1;
	163	if(lclx == 1) xstart = id;
	164	} // end if fMap
	165	if(lclx > 0 && fMap->TestHit(k,id) == kEmpty) {
	166	// end of cluster if a gap exists
	167	xstop = id-1;
	168	ilcl = 1;
	169	break;
	170	} // end if lclx >0 ...
	171	} // end of m-loop
	172	if(lclx == 0 && ilcl == 0) it = id; // no cluster in the window,
	173	// continue the "it" loop
	174	if(id >= fMaxNofSamples && lclx == 0) break; // the x row finished
	175	if(id < fMaxNofSamples && ilcl == 0 && lclx > 0) {
	176	// cluster end is outside of the window,
	177	mmax += 5; // increase mmax and repeat the cluster finding
	178	it -= 1;
	179	} // end if id<
	180	if(id >= fMaxNofSamples && lclx > 0) { // the x row finished but
	181	xstop = fMaxNofSamples - 1; // the end cluster exists
	182	ilcl = 1;
	183	} // end if id>=
	184	// --- Calculate z and x coordinates for one dimensional clusters
	185	if(ilcl == 1) { // new cluster exists
	186	it = id;
	187	mmax = 10;
	188	nofFoundClusters++;
	189	Float_t clusterCharge = 0.;
	190	Float_t zpitch = fSegmentation->Dpz(k+1);
	191	Int_t dummy = 0;
	192	Float_t clusterZ, dummyX;
	193	fSegmentation->GetPadCxz(dummy,k,dummyX,clusterZ);
	194	Float_t zstart = clusterZ - 0.5*zpitch;
	195	Float_t zstop = clusterZ + 0.5*zpitch;
	196	Float_t clusterX = 0.;
	197	Int_t clusterSizeX = lclx;
	198	Int_t clusterSizeZ = 1;
	199	Int_t its;
	200	for(its=xstart; its<=xstop; its++) {
	201	Int_t firedpixel=0;
	202	if (fMap->GetHitIndex(k,its)>=0) firedpixel=1;
203	clusterCharge += firedpixel;
204	clusterX += its + 0.5;
205	} // end for its
206	Float_t fRphiPitch = fSegmentation->Dpx(dummy);
207	clusterX /= (clusterSizeX/fRphiPitch); // center of gravity for x
208	// Write the points (coordinates and some cluster information)
209	// to the AliITSRawClusterSPD object
210
211	AliITSRawClusterSPD clust(clusterZ,clusterX,clusterCharge,
212	clusterSizeZ,clusterSizeX,
213	xstart,xstop,zstart,zstop,k,mod);
214	iTS->AddCluster(0,&clust);
215	} // new cluster (ilcl=1)
216	} // X direction loop (it)
217	} // Z direction loop (k)
218	//fMap->ClearMap();
219	return;
409f8c84	220	}
9e1e0cd7	221	//______________________________________________________________________
	222	void AliITSClusterFinderSPDdubna::GroupClusters(){
	223	// Find two dimensional clusters, i.e. group one dimensional clusters
	224	// into two dimensional ones (go both in x and z directions).
	225	// get number of clusters for this module
	226
	227	Int_t nofClusters = fClusters->GetEntriesFast();
	228	nofClusters -= fNclusters;
	229	AliITSRawClusterSPD *clusterI;
	230	AliITSRawClusterSPD *clusterJ;
	231	Int_t *label=new Int_t[nofClusters];
	232	Int_t i,j;
	233	for(i=0; i<nofClusters; i++) label[i] = 0;
	234	for(i=0; i<nofClusters; i++) {
	235	if(label[i] != 0) continue;
	236	for(j=i+1; j<nofClusters; j++) {
	237	if(label[j] != 0) continue;
	238	clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
	239	clusterJ = (AliITSRawClusterSPD*) fClusters->At(j);
	240	Bool_t pair = clusterI->Brother(clusterJ,fDz,fDx);
	241	if(pair) {
	242	#ifndef DEBUG
	243	cout << "clusters " << i << "," << j << " before grouping" << endl;
	244	clusterI->PrintInfo();
	245	clusterJ->PrintInfo();
	246	#endif
	247	clusterI->Add(clusterJ);
	248	#ifndef DEBUG
	249	cout << "remove cluster " << j << endl;
	250	#endif
	251	label[j] = 1;
	252	fClusters->RemoveAt(j);
	253	#ifndef DEBUG
	254	cout << "cluster " << i << " after grouping" << endl;
	255	clusterI->PrintInfo();
	256	#endif
	257	} // pair
	258	} // J clusters
	259	label[i] = 1;
	260	} // I clusters
	261	fClusters->Compress();
	262	#ifndef DEBUG
409f8c84	263	Int_t totalNofClusters = fClusters->GetEntriesFast();
409f8c84	264	cout << " Nomber of clusters at the group end ="<< totalNofClusters<<endl;
9e1e0cd7	265	#endif
	266	delete [] label;
	267	return;
409f8c84	268	}
9e1e0cd7	269	//______________________________________________________________________
	270	void AliITSClusterFinderSPDdubna::TracksInCluster(){
	271	// Find tracks creating one cluster
	272	// get number of clusters for this module
	273	Int_t nofClusters = fClusters->GetEntriesFast();
	274	Int_t i, ix, iz, jx, jz, xstart, xstop, zstart, zstop, nclx, nclz;
bf3f2830	275	const Int_t kTrmax = 100;
bf3f2830	276	Int_t cltracks[kTrmax], itr, tracki, ii, is, js, ie, ntr, tr0, tr1, tr2;
9e1e0cd7	277
	278	nofClusters -= fNclusters;
	279	for(i=0; i<nofClusters; i++) {
	280	ii = 0;
bf3f2830	281	memset(cltracks,-1,sizeof(int)*kTrmax);
9e1e0cd7	282	tr0=tr1=tr2=-1;
	283	AliITSRawClusterSPD clusterI = (AliITSRawClusterSPD) fClusters->At(i);
	284	nclx = clusterI->NclX();
	285	nclz = clusterI->NclZ();
	286	xstart = clusterI->XStartf();
	287	xstop = clusterI->XStopf();
	288	zstart = clusterI->Zend()-nclz+1;
	289	zstop = clusterI->Zend();
	290	Int_t ind;
	291	for(iz=0; iz<nclz; iz++) {
	292	jz = zstart + iz;
	293	for(ix=0; ix<nclx; ix++) {
	294	jx = xstart + ix;
	295	ind = fMap->GetHitIndex(jz,jx);
	296	if(ind < 0) continue;
	297	if(ind == 0 && iz >= 0 && ix > 0) continue;
	298	if(ind == 0 && iz > 0 && ix >= 0) continue;
	299	if(ind == 0 && iz == 0 && ix == 0 && i > 0) continue;
	300	AliITSdigitSPD dig = (AliITSdigitSPD)fMap->GetHit(jz,jx);
	301	for(itr=0; itr<3; itr++) {
	302	tracki = dig->fTracks[itr];
	303	if(tracki >= 0) {
	304	ii += 1;
	305	if(ii < 99) cltracks[ii-1] = tracki;
	306	} // end if tracki >-0
	307	} // end for itr
	308	} // ix pixel
	309	} // iz pixel
bf3f2830	310	for(is=0; is<kTrmax; is++) {
9e1e0cd7	311	if(cltracks[is]<0) continue;
bf3f2830	312	for(js=is+1; js<kTrmax; js++) {
9e1e0cd7	313	if(cltracks[js]<0) continue;
	314	if(cltracks[js]==cltracks[is]) cltracks[js]=-5;
	315	} // end for js
	316	} // end for is
	317	ntr = 0;
bf3f2830	318	for(ie=0; ie<kTrmax; ie++) {
9e1e0cd7	319	if(cltracks[ie] >= 0) {
	320	ntr=ntr+1;
	321	if(ntr==1) tr0=cltracks[ie];
	322	if(ntr==2) tr1=cltracks[ie];
	323	if(ntr==3) tr2=cltracks[ie];
	324	} // end if cltracks[ie] >-=
	325	} // end for ie
	326	// if delta ray only
	327	if(ntr == 0) ntr = 1;
	328	clusterI->SetNTracks(ntr);
	329	clusterI->SetTracks(tr0,tr1,tr2);
	330	} // I cluster
409f8c84	331	}
9e1e0cd7	332	//______________________________________________________________________
	333	void AliITSClusterFinderSPDdubna::GetRecPoints(){
	334	// get rec points
	335	static AliITS iTS = (AliITS)gAlice->GetModule("ITS");
	336	// get number of clusters for this module
	337	Int_t nofClusters = fClusters->GetEntriesFast();
	338	const Float_t kconv = 1.0e-4;
	339	const Float_t kRMSx = 12.0*kconv; // microns -> cm ITS TDR Table 1.3
	340	const Float_t kRMSz = 70.0*kconv; // microns -> cm ITS TDR Table 1.3
	341	Float_t spdLength = fSegmentation->Dz();
	342	Float_t spdWidth = fSegmentation->Dx();
	343	Int_t dummy = 0;
	344	Float_t xpitch = fSegmentation->Dpx(dummy);
	345	Int_t i;
	346	Int_t track0, track1, track2;
	347
	348	nofClusters -= fNclusters;
	349	for(i=0; i<nofClusters; i++) {
	350	AliITSRawClusterSPD clusterI = (AliITSRawClusterSPD) fClusters->At(i);
	351	Int_t clustersizex = clusterI->NclX();
	352	Int_t clustersizez = clusterI->NclZ();
	353	Int_t xstart = clusterI->XStartf();
	354	Float_t clusterx = 0;
	355	Float_t clusterz = clusterI->Z();
	356	for(Int_t ii=0; ii<clustersizex; ii++) clusterx += (xstart+0.5+ii)*xpitch;
	357	clusterx /= clustersizex;
	358	clusterz /= clustersizez;
	359	clusterI->GetTracks(track0, track1, track2);
	360	AliITSRecPoint rnew;
	361	rnew.SetX((clusterx - spdWidth/2)*kconv);
	362	rnew.SetZ((clusterz - spdLength/2)*kconv);
	363	rnew.SetQ(1.);
	364	rnew.SetdEdX(0.);
	365	rnew.SetSigmaX2(kRMSx*kRMSx);
	366	rnew.SetSigmaZ2(kRMSz*kRMSz);
	367	rnew.fTracks[0]=track0;
	368	rnew.fTracks[1]=track1;
	369	rnew.fTracks[2]=track2;
	370	iTS->AddRecPoint(rnew);
	371	} // I clusters
	372	fMap->ClearMap();
409f8c84	373	}
9e1e0cd7	374	//______________________________________________________________________
	375	void AliITSClusterFinderSPDdubna::FindRawClusters(Int_t mod){
	376	// find raw clusters
	377
	378	#ifndef DEBUG
	379	cout<<"SPDdubna: module ="<<mod<<endl;
	380	#endif
	381	Find1DClusters(mod);
	382	GroupClusters();
	383	TracksInCluster();
	384	GetRecPoints();
409f8c84	385	}