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