[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinder.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.                  *
 **************************************************************************/
////////////////////////////////////////////////////////////////////////////
//                                                                        //
// Base Class used to find                                                //
// the reconstructed points for ITS                                       //
// See also AliITSClusterFinderSPD, AliITSClusterFinderSDD,               //
// AliITSClusterFinderSDD  AliITSClusterFinderV2                          //
////////////////////////////////////////////////////////////////////////////

#include "AliRun.h"
#include "AliITSClusterFinder.h"
#include "AliITSRecPoint.h"
#include "AliITSdigit.h"
#include "AliITSDetTypeRec.h"
#include "AliITSMap.h"
#include "AliITSgeomTGeo.h"
#include <TParticle.h>
#include "AliMC.h"

ClassImp(AliITSClusterFinder)

extern AliRun *gAlice;

//----------------------------------------------------------------------
AliITSClusterFinder::AliITSClusterFinder():
TObject(),
fModule(0),
fDigits(0),
fNdigits(0),
fDetTypeRec(0),
fClusters(0),
fMap(0),
fNPeaks(-1),
fNModules(AliITSgeomTGeo::GetNModules()),
fEvent(0){
    // default cluster finder
    // Input:
    //   none.
    // Output:
    //   none.
    // Return:
    //   A default constructed AliITSCulsterFinder
}
//----------------------------------------------------------------------
AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
TObject(),
fModule(0),
fDigits(0),
fNdigits(0),
fDetTypeRec(dettyp),
fClusters(0),
fMap(0),
fNPeaks(-1),
fNModules(AliITSgeomTGeo::GetNModules()),
fEvent(0){
    // default cluster finder
    // Standard constructor for cluster finder
    // Input:
    //   AliITSsegmentation *seg  The segmentation class to be used
    //   AliITSresponse     *res  The response class to be used
    // Output:
    //   none.
    // Return:
    //   A Standard constructed AliITSCulsterFinder

}
//----------------------------------------------------------------------
AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
					 TClonesArray *digits):
TObject(),
fModule(0),
fDigits(digits),
fNdigits(0),
fDetTypeRec(dettyp),
fClusters(0),
fMap(0),
fNPeaks(-1),
fNModules(AliITSgeomTGeo::GetNModules()),
fEvent(0){
    // default cluster finder
    // Standard + cluster finder constructor
    // Input:
    //   AliITSsegmentation *seg  The segmentation class to be used
    //   AliITSresponse     *res  The response class to be used
    //   TClonesArray    *digits  Array of digits to be used
    // Output:
    //   none.
    // Return:
    //   A Standard constructed AliITSCulsterFinder

    fNdigits = fDigits->GetEntriesFast();
}

//______________________________________________________________________
AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : TObject(source),
fModule(source.fModule),
fDigits(),
fNdigits(source.fNdigits),
fDetTypeRec(),
fClusters(),
fMap(),
fNPeaks(source.fNPeaks),
fNModules(source.fNModules),
fEvent(source.fEvent) {
  // Copy constructor
  // Copies are not allowed. The method is protected to avoid misuse.
  AliError("Copy constructor not allowed\n");
}


//______________________________________________________________________
//AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
  // Assignment operator
  // Assignment is not allowed. The method is protected to avoid misuse.
//  Fatal("= operator","Assignment operator not allowed\n");
//  return *this;
//}

//----------------------------------------------------------------------
AliITSClusterFinder::~AliITSClusterFinder(){
    // destructor cluster finder
    // Input:
    //   none.
    // Output:
    //   none.
    // Return:
    //   none.

    if(fMap) {delete fMap;}
    // Zero local pointers. Other classes own these pointers.
    fMap          = 0;
    fDigits       = 0;
    fNdigits      = 0;
    fNPeaks       = 0;
    fDetTypeRec   = 0;

}
//__________________________________________________________________________
void AliITSClusterFinder::InitGeometry(){
 //
 // Initialisation of ITS geometry
 //
  Int_t mmax=AliITSgeomTGeo::GetNModules();
  for (Int_t m=0; m<mmax; m++) {
    Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
    fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
    fNlayer[m] = lay-1;
  }
}


//______________________________________________________________________
Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
    // Locagical function which checks to see if digit i has a neighbor.
    // If so, then it returns kTRUE and its neighbor index j.
    // This routine checks if the digits are side by side or one before the
    // other. Requires that the array of digits be in proper order.
    // Returns kTRUE in the following cases.
    //                 ji   0j   if kdiagonal  j0    0i
    //                 00   0i   if kdiagonal  0i    j0
    // Inputs:
    //    TObjArray *digs   Array to search for neighbors in
    //    Int_t      i      Index of digit for which we are searching for
    //                      a neighbor of.
    // Output:
    //    Int_t      j[4]   Index of one or more of the digits which is a
    //                      neighbor of digit a index i.
    // Return:
    //    Bool_t            kTRUE if a neighbor was found kFALSE otherwise.
    Int_t ix,jx,iz,jz,j;
    const Bool_t kdiagonal=kFALSE;
    Bool_t nei[4];

    // No neighbors found if array empty.
    if(digs->GetEntriesFast()<=0) return kFALSE;
    // can not be a digit with first element or elements out or range
    if(i<=0 || i>=digs->GetEntriesFast()) return kFALSE;

    for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
    ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
    iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
    for(j=0;j<i;j++){
        jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
        jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
        if(jx+1==ix && jz  ==iz){n[0] = j;nei[0] = kTRUE;}
        if(jx  ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
        if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
        if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
    } // end for k
    if(nei[0]||nei[1]) return kTRUE;
    if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
    // no Neighbors found.
    return kFALSE;
}

//______________________________________________________________________
void AliITSClusterFinder::Print(ostream *os) const{
    //Standard output format for this class
    // Inputs:
    //    ostream *os   Output stream
    // Output:
    //    ostream *os   Output stream
    // Return:
    //    none.

    *os << fModule<<",";
    *os << fNdigits<<",";
    *os << fNPeaks<<endl;
}
//______________________________________________________________________
void AliITSClusterFinder::Read(istream *is)  {
    //Standard input for this class
    // Inputs:
    //    istream *is   Input stream
    // Output:
    //    istream *is   Input stream
    // Return:
    //    none.

    *is >> fModule;
    *is >> fNdigits;
    *is >> fNPeaks;
}
//______________________________________________________________________
ostream &operator<<(ostream &os,AliITSClusterFinder &source){
    // Standard output streaming function.
    // Inputs:
    //    ostream             *os     Output stream
    //    AliITSClusterFinder &source Class to be printed
    // Output:
    //    ostream             *os     Output stream
    // Return:
    //    none.

    source.Print(&os);
    return os;
}
//______________________________________________________________________
istream &operator>>(istream &is,AliITSClusterFinder &source){
    // Standard output streaming function.
    // Inputs:
    //    istream              *is      Input stream
    //     AliITSClusterFinder &source  Class to be read in.
    // Output:
    //    istream              *is      Input stream
    // Return:
    //    none.

    source.Read(&is);
    return is;
}
//______________________________________________________________________
void AliITSClusterFinder::CheckLabels2(Int_t lab[10]) {
  //------------------------------------------------------------
  // Tries to find mother's labels
  //------------------------------------------------------------
  AliRunLoader *rl = AliRunLoader::Instance();
  if(!rl) return;
  TTree *trK=(TTree*)rl->TreeK();

  if(trK){
    Int_t nlabels =0; 
    for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
    if(nlabels == 0) return; // In case of no labels just exit


    Int_t ntracks = gAlice->GetMCApp()->GetNtrack();

    for (Int_t i=0;i<10;i++){
      Int_t label = lab[i];
      if (label>=0 && label<ntracks) {
	TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);

	if (part->P() < 0.02) {
	  Int_t m=part->GetFirstMother();
	  if (m<0) {	
	    continue;
	  }
	  if (part->GetStatusCode()>0) {
	    continue;
	  }
	  lab[i]=m;       
	}
	else
	  if (part->P() < 0.12 && nlabels>3) {
	    lab[i]=-2;
	    nlabels--;
	  } 
      }
      else{
	if ( (label>ntracks||label <0) && nlabels>3) {
	  lab[i]=-2;
	  nlabels--;
	} 
      }
    }  
    if (nlabels>3){
      for (Int_t i=0;i<10;i++){
	if (nlabels>3){
	  Int_t label = lab[i];
	  if (label>=0 && label<ntracks) {
	    TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
	    if (part->P() < 0.1) {
	      lab[i]=-2;
	      nlabels--;
	    }
	  }
	}
      }
    }

    //compress labels -- if multi-times the same
    Int_t lab2[10];
    for (Int_t i=0;i<10;i++) lab2[i]=-2;
    for (Int_t i=0;i<10  ;i++){
      if (lab[i]<0) continue;
      for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
	if (lab2[j]<0) {
	  lab2[j]= lab[i];
	  break;
	}
      }
    }
    for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
  
  }
}

//______________________________________________________________________
void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
  //add label to the cluster
  AliRunLoader *rl = AliRunLoader::Instance();
  TTree *trK=(TTree*)rl->TreeK();
  if(trK){
    if(label<0) return; // In case of no label just exit

    Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
    if (label>ntracks) return;
    for (Int_t i=0;i<10;i++){
      //    if (label<0) break;
      if (lab[i]==label) break;
      if (lab[i]<0) {
	lab[i]= label;
	break;
      }
    }
  }
}


//______________________________________________________________________
void AliITSClusterFinder:: 
FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
  //------------------------------------------------------------
  // returns an array of indices of digits belonging to the cluster
  // (needed when the segmentation is not regular) 
  //------------------------------------------------------------
  if (n<200) idx[n++]=bins[k].GetIndex();
  bins[k].Use();

  if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
  if (bins[k-1   ].IsNotUsed()) FindCluster(k-1   ,maxz,bins,n,idx);
  if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
  if (bins[k+1   ].IsNotUsed()) FindCluster(k+1   ,maxz,bins,n,idx);
  /*
  if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
  if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
  if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
  if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
  */
}

//______________________________________________________________________
Bool_t AliITSClusterFinder::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
  //------------------------------------------------------------
  //is this a local maximum ?
  //------------------------------------------------------------
  UShort_t q=bins[k].GetQ();
  if (q==1023) return kFALSE;
  if (bins[k-max].GetQ() > q) return kFALSE;
  if (bins[k-1  ].GetQ() > q) return kFALSE; 
  if (bins[k+max].GetQ() > q) return kFALSE; 
  if (bins[k+1  ].GetQ() > q) return kFALSE; 
  if (bins[k-max-1].GetQ() > q) return kFALSE;
  if (bins[k+max-1].GetQ() > q) return kFALSE; 
  if (bins[k+max+1].GetQ() > q) return kFALSE; 
  if (bins[k-max+1].GetQ() > q) return kFALSE;
  return kTRUE; 
}

//______________________________________________________________________
void AliITSClusterFinder::
FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
  //------------------------------------------------------------
  //find local maxima
  //------------------------------------------------------------
  if (n<31)
  if (IsMaximum(k,max,b)) {
    idx[n]=k; msk[n]=(2<<n);
    n++;
  }
  b[k].SetMask(0);
  if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
  if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
  if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
  if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
}

//______________________________________________________________________
void AliITSClusterFinder::
MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
  //------------------------------------------------------------
  //mark this peak
  //------------------------------------------------------------
  UShort_t q=bins[k].GetQ();

  bins[k].SetMask(bins[k].GetMask()|m); 

  if (bins[k-max].GetQ() <= q)
     if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
  if (bins[k-1  ].GetQ() <= q)
     if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
  if (bins[k+max].GetQ() <= q)
     if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
  if (bins[k+1  ].GetQ() <= q)
     if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
}

//______________________________________________________________________
void AliITSClusterFinder::
MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSRecPoint &c) {
  //------------------------------------------------------------
  //make cluster using digits of this peak
  //------------------------------------------------------------
  Float_t q=(Float_t)bins[k].GetQ();
  Int_t i=k/max, j=k-i*max;

  c.SetQ(c.GetQ()+q);
  c.SetY(c.GetY()+i*q);
  c.SetZ(c.GetZ()+j*q);
  c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
  c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);

  bins[k].SetMask(0xFFFFFFFE);
  
  if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
  if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
  if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
  if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
}
Commit	Line	Data
	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	// Base Class used to find //
	18	// the reconstructed points for ITS //
	19	// See also AliITSClusterFinderSPD, AliITSClusterFinderSDD, //
	20	// AliITSClusterFinderSDD AliITSClusterFinderV2 //
	21	////////////////////////////////////////////////////////////////////////////
	22
	23	#include "AliRun.h"
	24	#include "AliITSClusterFinder.h"
	25	#include "AliITSRecPoint.h"
	26	#include "AliITSdigit.h"
	27	#include "AliITSDetTypeRec.h"
	28	#include "AliITSMap.h"
	29	#include "AliITSgeomTGeo.h"
	30	#include <TParticle.h>
	31	#include "AliMC.h"
	32
	33	ClassImp(AliITSClusterFinder)
	34
	35	extern AliRun *gAlice;
	36
	37	//----------------------------------------------------------------------
	38	AliITSClusterFinder::AliITSClusterFinder():
	39	TObject(),
	40	fModule(0),
	41	fDigits(0),
	42	fNdigits(0),
	43	fDetTypeRec(0),
	44	fClusters(0),
	45	fMap(0),
	46	fNPeaks(-1),
	47	fNModules(AliITSgeomTGeo::GetNModules()),
	48	fEvent(0){
	49	// default cluster finder
	50	// Input:
	51	// none.
	52	// Output:
	53	// none.
	54	// Return:
	55	// A default constructed AliITSCulsterFinder
	56	}
	57	//----------------------------------------------------------------------
	58	AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
	59	TObject(),
	60	fModule(0),
	61	fDigits(0),
	62	fNdigits(0),
	63	fDetTypeRec(dettyp),
	64	fClusters(0),
	65	fMap(0),
	66	fNPeaks(-1),
	67	fNModules(AliITSgeomTGeo::GetNModules()),
	68	fEvent(0){
	69	// default cluster finder
	70	// Standard constructor for cluster finder
	71	// Input:
	72	// AliITSsegmentation *seg The segmentation class to be used
	73	// AliITSresponse *res The response class to be used
	74	// Output:
	75	// none.
	76	// Return:
	77	// A Standard constructed AliITSCulsterFinder
	78
	79	}
	80	//----------------------------------------------------------------------
	81	AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
	82	TClonesArray *digits):
	83	TObject(),
	84	fModule(0),
	85	fDigits(digits),
	86	fNdigits(0),
	87	fDetTypeRec(dettyp),
	88	fClusters(0),
	89	fMap(0),
	90	fNPeaks(-1),
	91	fNModules(AliITSgeomTGeo::GetNModules()),
	92	fEvent(0){
	93	// default cluster finder
	94	// Standard + cluster finder constructor
	95	// Input:
	96	// AliITSsegmentation *seg The segmentation class to be used
	97	// AliITSresponse *res The response class to be used
	98	// TClonesArray *digits Array of digits to be used
	99	// Output:
	100	// none.
	101	// Return:
	102	// A Standard constructed AliITSCulsterFinder
	103
	104	fNdigits = fDigits->GetEntriesFast();
	105	}
	106
	107	//______________________________________________________________________
	108	AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : TObject(source),
	109	fModule(source.fModule),
	110	fDigits(),
	111	fNdigits(source.fNdigits),
	112	fDetTypeRec(),
	113	fClusters(),
	114	fMap(),
	115	fNPeaks(source.fNPeaks),
	116	fNModules(source.fNModules),
	117	fEvent(source.fEvent) {
	118	// Copy constructor
	119	// Copies are not allowed. The method is protected to avoid misuse.
	120	AliError("Copy constructor not allowed\n");
	121	}
	122
	123
	124	//______________________________________________________________________
	125	//AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
	126	// Assignment operator
	127	// Assignment is not allowed. The method is protected to avoid misuse.
	128	// Fatal("= operator","Assignment operator not allowed\n");
	129	// return *this;
	130	//}
	131
	132	//----------------------------------------------------------------------
	133	AliITSClusterFinder::~AliITSClusterFinder(){
	134	// destructor cluster finder
	135	// Input:
	136	// none.
	137	// Output:
	138	// none.
	139	// Return:
	140	// none.
	141
	142	if(fMap) {delete fMap;}
	143	// Zero local pointers. Other classes own these pointers.
	144	fMap = 0;
	145	fDigits = 0;
	146	fNdigits = 0;
	147	fNPeaks = 0;
	148	fDetTypeRec = 0;
	149
	150	}
	151	//__________________________________________________________________________
	152	void AliITSClusterFinder::InitGeometry(){
	153	//
	154	// Initialisation of ITS geometry
	155	//
	156	Int_t mmax=AliITSgeomTGeo::GetNModules();
	157	for (Int_t m=0; m<mmax; m++) {
	158	Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
	159	fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
	160	fNlayer[m] = lay-1;
	161	}
	162	}
	163
	164
	165
	166
	167	//______________________________________________________________________
	168	Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
	169	// Locagical function which checks to see if digit i has a neighbor.
	170	// If so, then it returns kTRUE and its neighbor index j.
	171	// This routine checks if the digits are side by side or one before the
	172	// other. Requires that the array of digits be in proper order.
	173	// Returns kTRUE in the following cases.
	174	// ji 0j if kdiagonal j0 0i
	175	// 00 0i if kdiagonal 0i j0
	176	// Inputs:
	177	// TObjArray *digs Array to search for neighbors in
	178	// Int_t i Index of digit for which we are searching for
	179	// a neighbor of.
	180	// Output:
	181	// Int_t j[4] Index of one or more of the digits which is a
	182	// neighbor of digit a index i.
	183	// Return:
	184	// Bool_t kTRUE if a neighbor was found kFALSE otherwise.
	185	Int_t ix,jx,iz,jz,j;
	186	const Bool_t kdiagonal=kFALSE;
	187	Bool_t nei[4];
	188
	189	// No neighbors found if array empty.
	190	if(digs->GetEntriesFast()<=0) return kFALSE;
	191	// can not be a digit with first element or elements out or range
	192	if(i<=0 \|\| i>=digs->GetEntriesFast()) return kFALSE;
	193
	194	for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
	195	ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
	196	iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
	197	for(j=0;j<i;j++){
	198	jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
	199	jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
	200	if(jx+1==ix && jz ==iz){n[0] = j;nei[0] = kTRUE;}
	201	if(jx ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
	202	if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
	203	if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
	204	} // end for k
	205	if(nei[0]\|\|nei[1]) return kTRUE;
	206	if(kdiagonal&&(nei[2]\|\|nei[3])) return kTRUE;
	207	// no Neighbors found.
	208	return kFALSE;
	209	}
	210
	211	//______________________________________________________________________
	212	void AliITSClusterFinder::Print(ostream *os) const{
	213	//Standard output format for this class
	214	// Inputs:
	215	// ostream *os Output stream
	216	// Output:
	217	// ostream *os Output stream
	218	// Return:
	219	// none.
	220
	221	*os << fModule<<",";
	222	*os << fNdigits<<",";
	223	*os << fNPeaks<<endl;
	224	}
	225	//______________________________________________________________________
	226	void AliITSClusterFinder::Read(istream *is) {
	227	//Standard input for this class
	228	// Inputs:
	229	// istream *is Input stream
	230	// Output:
	231	// istream *is Input stream
	232	// Return:
	233	// none.
	234
	235	*is >> fModule;
	236	*is >> fNdigits;
	237	*is >> fNPeaks;
	238	}
	239	//______________________________________________________________________
	240	ostream &operator<<(ostream &os,AliITSClusterFinder &source){
	241	// Standard output streaming function.
	242	// Inputs:
	243	// ostream *os Output stream
	244	// AliITSClusterFinder &source Class to be printed
	245	// Output:
	246	// ostream *os Output stream
	247	// Return:
	248	// none.
	249
	250	source.Print(&os);
	251	return os;
	252	}
	253	//______________________________________________________________________
	254	istream &operator>>(istream &is,AliITSClusterFinder &source){
	255	// Standard output streaming function.
	256	// Inputs:
	257	// istream *is Input stream
	258	// AliITSClusterFinder &source Class to be read in.
	259	// Output:
	260	// istream *is Input stream
	261	// Return:
	262	// none.
	263
	264	source.Read(&is);
	265	return is;
	266	}
	267	//______________________________________________________________________
	268	void AliITSClusterFinder::CheckLabels2(Int_t lab[10]) {
	269	//------------------------------------------------------------
	270	// Tries to find mother's labels
	271	//------------------------------------------------------------
	272	AliRunLoader *rl = AliRunLoader::Instance();
	273	if(!rl) return;
	274	TTree trK=(TTree)rl->TreeK();
	275
	276	if(trK){
	277	Int_t nlabels =0;
	278	for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
	279	if(nlabels == 0) return; // In case of no labels just exit
	280
	281
	282	Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
	283
	284	for (Int_t i=0;i<10;i++){
	285	Int_t label = lab[i];
	286	if (label>=0 && label<ntracks) {
	287	TParticle part=(TParticle)gAlice->GetMCApp()->Particle(label);
	288
	289	if (part->P() < 0.02) {
	290	Int_t m=part->GetFirstMother();
	291	if (m<0) {
	292	continue;
	293	}
	294	if (part->GetStatusCode()>0) {
	295	continue;
	296	}
	297	lab[i]=m;
	298	}
	299	else
	300	if (part->P() < 0.12 && nlabels>3) {
	301	lab[i]=-2;
	302	nlabels--;
	303	}
	304	}
	305	else{
	306	if ( (label>ntracks\|\|label <0) && nlabels>3) {
	307	lab[i]=-2;
	308	nlabels--;
	309	}
	310	}
	311	}
	312	if (nlabels>3){
	313	for (Int_t i=0;i<10;i++){
	314	if (nlabels>3){
	315	Int_t label = lab[i];
	316	if (label>=0 && label<ntracks) {
	317	TParticle part=(TParticle)gAlice->GetMCApp()->Particle(label);
	318	if (part->P() < 0.1) {
	319	lab[i]=-2;
	320	nlabels--;
	321	}
	322	}
	323	}
	324	}
	325	}
	326
	327	//compress labels -- if multi-times the same
	328	Int_t lab2[10];
	329	for (Int_t i=0;i<10;i++) lab2[i]=-2;
	330	for (Int_t i=0;i<10 ;i++){
	331	if (lab[i]<0) continue;
	332	for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
	333	if (lab2[j]<0) {
	334	lab2[j]= lab[i];
	335	break;
	336	}
	337	}
	338	}
	339	for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
	340
	341	}
	342	}
	343
	344	//______________________________________________________________________
	345	void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
	346	//add label to the cluster
	347	AliRunLoader *rl = AliRunLoader::Instance();
	348	TTree trK=(TTree)rl->TreeK();
	349	if(trK){
	350	if(label<0) return; // In case of no label just exit
	351
	352	Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
	353	if (label>ntracks) return;
	354	for (Int_t i=0;i<10;i++){
	355	// if (label<0) break;
	356	if (lab[i]==label) break;
	357	if (lab[i]<0) {
	358	lab[i]= label;
	359	break;
	360	}
	361	}
	362	}
	363	}
	364
	365
	366	//______________________________________________________________________
	367	void AliITSClusterFinder::
	368	FindCluster(Int_t k,Int_t maxz,AliBin bins,Int_t &n,Int_t idx) {
	369	//------------------------------------------------------------
	370	// returns an array of indices of digits belonging to the cluster
	371	// (needed when the segmentation is not regular)
	372	//------------------------------------------------------------
	373	if (n<200) idx[n++]=bins[k].GetIndex();
	374	bins[k].Use();
	375
	376	if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
	377	if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx);
	378	if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
	379	if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx);
	380	/*
	381	if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
	382	if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
	383	if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
	384	if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
	385	*/
	386	}
	387
	388	//______________________________________________________________________
	389	Bool_t AliITSClusterFinder::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
	390	//------------------------------------------------------------
	391	//is this a local maximum ?
	392	//------------------------------------------------------------
	393	UShort_t q=bins[k].GetQ();
	394	if (q==1023) return kFALSE;
	395	if (bins[k-max].GetQ() > q) return kFALSE;
	396	if (bins[k-1 ].GetQ() > q) return kFALSE;
	397	if (bins[k+max].GetQ() > q) return kFALSE;
	398	if (bins[k+1 ].GetQ() > q) return kFALSE;
	399	if (bins[k-max-1].GetQ() > q) return kFALSE;
	400	if (bins[k+max-1].GetQ() > q) return kFALSE;
	401	if (bins[k+max+1].GetQ() > q) return kFALSE;
	402	if (bins[k-max+1].GetQ() > q) return kFALSE;
	403	return kTRUE;
	404	}
	405
	406	//______________________________________________________________________
	407	void AliITSClusterFinder::
	408	FindPeaks(Int_t k,Int_t max,AliBin b,Int_t idx,UInt_t *msk,Int_t& n) {
	409	//------------------------------------------------------------
	410	//find local maxima
	411	//------------------------------------------------------------
	412	if (n<31)
	413	if (IsMaximum(k,max,b)) {
	414	idx[n]=k; msk[n]=(2<<n);
	415	n++;
	416	}
	417	b[k].SetMask(0);
	418	if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
	419	if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
	420	if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
	421	if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
	422	}
	423
	424	//______________________________________________________________________
	425	void AliITSClusterFinder::
	426	MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
	427	//------------------------------------------------------------
	428	//mark this peak
	429	//------------------------------------------------------------
	430	UShort_t q=bins[k].GetQ();
	431
	432	bins[k].SetMask(bins[k].GetMask()\|m);
	433
	434	if (bins[k-max].GetQ() <= q)
	435	if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
	436	if (bins[k-1 ].GetQ() <= q)
	437	if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
	438	if (bins[k+max].GetQ() <= q)
	439	if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
	440	if (bins[k+1 ].GetQ() <= q)
	441	if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
	442	}
	443
	444	//______________________________________________________________________
	445	void AliITSClusterFinder::
	446	MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSRecPoint &c) {
	447	//------------------------------------------------------------
	448	//make cluster using digits of this peak
	449	//------------------------------------------------------------
	450	Float_t q=(Float_t)bins[k].GetQ();
	451	Int_t i=k/max, j=k-i*max;
	452
	453	c.SetQ(c.GetQ()+q);
	454	c.SetY(c.GetY()+i*q);
	455	c.SetZ(c.GetZ()+j*q);
	456	c.SetSigmaY2(c.GetSigmaY2()+iiq);
	457	c.SetSigmaZ2(c.GetSigmaZ2()+jjq);
	458
	459	bins[k].SetMask(0xFFFFFFFE);
	460
	461	if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
	462	if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
	463	if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
	464	if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);
	465	}