[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 "AliITSClusterFinder.h"
#include "AliITSclusterV2.h"
#include "AliITSRecPoint.h"
#include "AliITSdigitSPD.h"
#include "AliITSdigitSDD.h"
#include "AliITSdigitSSD.h"
#include "AliITSgeom.h"
#include "AliITSMap.h"
#include "AliRun.h"
#include "AliMC.h"
#include "AliITS.h"
#include "TParticle.h"

ClassImp(AliITSClusterFinder)

//----------------------------------------------------------------------
AliITSClusterFinder::AliITSClusterFinder():
TObject(),
fDebug(0),
fModule(0),
fDigits(0),
fNdigits(0),
fResponse(0),
fSegmentation(0),
fClusters(0),
fNRawClusters(0),
fMap(0),
fNperMax(0),
fDeclusterFlag(0),
fClusterSize(0),
fNPeaks(-1){
    // default cluster finder
    // Input:
    //   none.
    // Output:
    //   none.
    // Return:
    //   A default constructed AliITSCulsterFinder
  Init();
}
//----------------------------------------------------------------------
AliITSClusterFinder::AliITSClusterFinder(AliITSsegmentation *seg, 
                                         AliITSresponse *res):
TObject(),
fDebug(0),
fModule(0),
fDigits(0),
fNdigits(0),
fResponse(res),
fSegmentation(seg),
fClusters(0),
fNRawClusters(0),
fMap(0),
fNperMax(0),
fDeclusterFlag(0),
fClusterSize(0),
fNPeaks(-1){
    // 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

    SetNperMax();
    SetClusterSize();
    SetDeclusterFlag();
    Init();
}
//----------------------------------------------------------------------
AliITSClusterFinder::AliITSClusterFinder(AliITSsegmentation *seg, 
                                         AliITSresponse *response, 
                                         TClonesArray *digits):
TObject(),
fDebug(0),
fModule(0),
fDigits(digits),
fNdigits(0),
fResponse(response),
fSegmentation(seg),
fClusters(0),
fNRawClusters(0),
fMap(0),
fNperMax(0),
fDeclusterFlag(0),
fClusterSize(0),
fNPeaks(-1){
    // 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();
    SetNperMax();
    SetClusterSize();
    SetDeclusterFlag();
    Init();
}
//______________________________________________________________________
AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : TObject(source) {
  // Copy constructor
  // Copies are not allowed. The method is protected to avoid misuse.
  Fatal("AliITSClusterFinder","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.
    fSegmentation = 0;
    fResponse     = 0;
    fMap          = 0;
    fDigits       = 0;
    fNdigits      = 0;
    fNRawClusters = 0;
    fNperMax      = 0;
    fDeclusterFlag= 0;
    fClusterSize  = 0;
    fNPeaks       = 0;
    fITS          = 0;
}

//__________________________________________________________________________
void AliITSClusterFinder::Init(){

  //Initialisation of ITS geometry

  fITS = (AliITS*)gAlice->GetModule("ITS");
  AliITSgeom *geom=(AliITSgeom*)fITS->GetITSgeom();

  Int_t mmax=geom->GetIndexMax();
  if (mmax>2200) {
    Fatal("AliITSClusterFinder","Too many ITS subdetectors !"); 
  }
  Int_t m;
  for (m=0; m<mmax; m++) {
    Int_t lay,lad,det; geom->GetModuleId(m,lay,lad,det);
    Float_t x,y,z;     geom->GetTrans(lay,lad,det,x,y,z); 
    Double_t rot[9];   geom->GetRotMatrix(lay,lad,det,rot);
    Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
    Double_t ca=TMath::Cos(alpha), sa=TMath::Sin(alpha);
    fYshift[m] = x*ca + y*sa;
    fZshift[m] = (Double_t)z;
    fNdet[m] = (lad-1)*geom->GetNdetectors(lay) + (det-1);
    fNlayer[m] = lay-1;
  }

}


//----------------------------------------------------------------------
void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c){
    // Add a raw cluster copy to the list
    // Input:
    //   Int_t       branch  The branch to which the cluster is to be added to
    //   AliITSRawCluster *c The cluster to be added to the array of clusters
    // Output:
    //   none.
    // Return:
    //   none.

    fITS->AddCluster(branch,c); 
    fNRawClusters++;
}
//----------------------------------------------------------------------
void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c, 
				     AliITSRecPoint &rp){
    // Add a raw cluster copy to the list and the RecPoint
    // Input:
    //   Int_t       branch  The branch to which the cluster is to be added to
    //   AliITSRawCluster *c The cluster to be added to the array of clusters
    //   AliITSRecPoint  &rp The RecPoint to be added to the array of RecPoints
    // Output:
    //   none.
    // Return:
    //   none.

    fITS->AddCluster(branch,c); 
    fNRawClusters++;
    fITS->AddRecPoint(rp); 
}

//______________________________________________________________________
void AliITSClusterFinder::CheckLabels(Int_t lab[3]) {
  //------------------------------------------------------------
  // Tries to find mother's labels
  //------------------------------------------------------------

  if(lab[0]<0 && lab[1]<0 && lab[2]<0) return; // In case of no labels just exit
  // Check if simulation
  AliMC* mc = gAlice->GetMCApp();
  if(!mc)return;

  Int_t ntracks = mc->GetNtrack();
  for (Int_t i=0;i<3;i++){
    Int_t label = lab[i];
    if (label>=0 && label<ntracks) {
      TParticle *part=(TParticle*)mc->Particle(label);
      if (part->P() < 0.005) {
	Int_t m=part->GetFirstMother();
	if (m<0) {	
	  continue;
	}
	if (part->GetStatusCode()>0) {
	  continue;
	}
	lab[i]=m;       
      }
    }    
  }
  
}

//______________________________________________________________________
void AliITSClusterFinder::FindRawClusters(Int_t module){
    // Default Cluster finder.
    // Input:
    //   Int_t module   Module number for which culster are to be found.
    // Output:
    //   none.
    // Return:
    //   none.
    const Int_t kelms = 10;
    Int_t ndigits = fDigits->GetEntriesFast();
    TObjArray *digs = new TObjArray(ndigits);
    TObjArray *clusts = new TObjArray(ndigits); // max # cluster
    TObjArray *clust0=0; // A spacific cluster of digits
    TObjArray *clust1=0; // A spacific cluster of digits
    AliITSdigit *dig=0; // locat pointer to a digit
    Int_t i=0,nc=0,j[4],k,k2=0;

    // Copy all digits for this module into a local TObjArray.
    for(i=0;i<ndigits;i++) digs->AddAt(new AliITSdigit(*(GetDigit(i))),i);
    digs->Sort();
    // First digit is a cluster.
    i  = 0;
    nc = 0;
    clusts->AddAt(new TObjArray(kelms),nc);
    clust0 = (TObjArray*)(clusts->At(nc));
    clust0->AddAtFree(digs->At(i)); // move owner ship from digs to clusts
    nc++;
    for(i=1;i<ndigits;i++){
        if(IsNeighbor(digs,i,j)){
            dig = (AliITSdigit*)(digs->At(j[0]));
            // Add to existing cluster. Find which cluster this digis 
            for(k=0;k<nc;k++){
                clust0 = ((TObjArray*)(clusts->At(k)));
                if(clust0->IndexOf(dig)>=0) break;
            } // end for k
            if(k>=nc){
                Fatal("FindRawClusters","Digit not found as expected");
            } // end if
            if(j[1]>=0){
                dig = (AliITSdigit*)(digs->At(j[1]));
                // Add to existing cluster. Find which cluster this digis 
                for(k2=0;k2<nc;k2++){
                    clust1 = ((TObjArray*)(clusts->At(k2)));
                    if(clust1->IndexOf(dig)>=0) break;
                } // end for k2
                if(k2>=nc){
                    Fatal("FindRawClusters","Digit not found as expected");
                } // end if
            } // end if j[1]>=0
            // Found cluster with neighboring digits add this one to it.
            if(clust0==clust1){ // same cluster
                clust0->AddAtFree(digs->At(i));
                clust0 = 0; // finished with cluster. zero for safty
                clust1 = 0; // finished wit hcluster. zero for safty
            }else{ // two different clusters which need to be merged.
                clust0->AddAtFree(digs->At(i)); // Add digit to this cluster.
                for(k=0;k<clust1->GetEntriesFast();k++){
                    // move clust1 into clust0
                    //move digit to this cluster
                    clust0->AddAtFree(clust1->At(k));
                    clust1->AddAt(0,k); // zero this one
                } // end for k
                delete clust1;
                clusts->AddAt(0,k2); // zero array of clusters element clust1
                clust0 = 0; // finished with cluster. zero for safty
                clust1 = 0; // finished wit hcluster. zero for safty
            } // end if clust0==clust1
        }else{// New cluster
            clusts->AddAt(new TObjArray(kelms),nc);
            clust0 = ((TObjArray*)(clusts->At(nc)));
            // move owner ship from digs to clusts
            clust0->AddAtFree(digs->At(i));
            clust0 = 0; // finished with cluster. zero for safty
            nc++;
        } // End if IsNeighbor
    } // end for i
    // There are now nc clusters in clusts. Each element of clust is an
    // array of digits which are clustered together.

    // For each cluster call detector specific CreateRecPoints
    for(i=0;i<nc;i++) CreateRecPoints((TObjArray*)(clusts->At(i)),module);

    // clean up at the end.
    for(i=0;i<nc;i++){ 
        clust0 =(TObjArray*)(clusts->At(i));
        // Digits deleted below, so zero this TObjArray
        for(k=0;k<clust0->GetEntriesFast();k++) clust0->AddAt(0,k);
        delete clust0; // Delete this TObjArray
        clusts->AddAt(0,i); // Contents deleted above, so zero it.
    } // end for i
    delete clusts; // Delete this TObjArray/
    // Delete the digits then the TObjArray which containted them.
    for(i=0;i<ndigits;i++) delete ((AliITSdigit*)(digs->At(i)));
    delete digs;
}
//______________________________________________________________________
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){
    //Standard output format for this class
    // Inputs:
    //    ostream *os   Output stream
    // Output:
    //    ostream *os   Output stream
    // Return:
    //    none.

    *os << fDebug<<",";
    *os << fModule<<",";
    *os << fNdigits<<",";
    *os << fNRawClusters<<",";
    *os << fNperMax<<",";
    *os << fDeclusterFlag<<",";
    *os << fClusterSize<<",";
    *os << fNPeaks<<endl;
}
//______________________________________________________________________
void AliITSClusterFinder::RecPoints2Clusters
(const TClonesArray *points, Int_t idx, TClonesArray *clusters) {
  //------------------------------------------------------------
  // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS 
  // subdetector indexed by idx 
  //------------------------------------------------------------
  AliITSgeom* geom = (AliITSgeom*)fITS->GetITSgeom();
  Int_t lastSPD1=geom->GetModuleIndex(2,1,1)-1;
  TClonesArray &cl=*clusters;
  Int_t ncl=points->GetEntriesFast();
  for (Int_t i=0; i<ncl; i++) {
    AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i);
    Float_t lp[5];
    lp[0]=-(-p->GetX()+fYshift[idx]); if (idx<=lastSPD1) lp[0]*=-1; //SPD1
    lp[1]=  -p->GetZ()+fZshift[idx];
    lp[2]=p->GetSigmaX2();
    lp[3]=p->GetSigmaZ2();
    lp[4]=p->GetQ()*36./23333.;  //electrons -> ADC
    Int_t lab[4]; 
    lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
    lab[3]=fNdet[idx];
    CheckLabels(lab);
    Int_t dummy[3]={0,0,0};
    new (cl[i]) AliITSclusterV2(lab,lp, dummy);
  }  
} 

//______________________________________________________________________
void AliITSClusterFinder::Read(istream *is){
    //Standard input for this class
    // Inputs:
    //    istream *is   Input stream
    // Output:
    //    istream *is   Input stream
    // Return:
    //    none.

    *is >> fDebug;
    *is >> fModule;
    *is >> fNdigits;
    *is >> fNRawClusters;
    *is >> fNperMax;
    *is >> fDeclusterFlag;
    *is >> fClusterSize;
    *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::RecPoints2Clusters
(const TClonesArray *points, Int_t idx, TClonesArray *clusters) {
  //------------------------------------------------------------
  // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS 
  // subdetector indexed by idx 
  //------------------------------------------------------------
  TClonesArray &cl=*clusters;
  Int_t ncl=points->GetEntriesFast();
  for (Int_t i=0; i<ncl; i++) {
    AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i);
    Float_t lp[5];
    lp[0]=-(-p->GetX()+fYshift[idx]); if (idx<=fLastSPD1) lp[0]*=-1; //SPD1
    lp[1]=  -p->GetZ()+fZshift[idx];
    lp[2]=p->GetSigmaX2();
    lp[3]=p->GetSigmaZ2();
    lp[4]=p->GetQ()*36./23333.;  //electrons -> ADC
    Int_t lab[4]; 
    lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
    lab[3]=fNdet[idx];
    CheckLabels(lab);
    Int_t dummy[3]={0,0,0};
    new (cl[i]) AliITSclusterV2(lab,lp, dummy);
  }  
} 
*/
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
b0f5e3fc	23	#include "AliITSClusterFinder.h"
04366a57	24	#include "AliITSclusterV2.h"
04366a57	25	#include "AliITSRecPoint.h"
e869281d	26	#include "AliITSdigitSPD.h"
	27	#include "AliITSdigitSDD.h"
	28	#include "AliITSdigitSSD.h"
04366a57	29	#include "AliITSgeom.h"
aacedc3e	30	#include "AliITSMap.h"
b0f5e3fc	31	#include "AliRun.h"
04366a57	32	#include "AliMC.h"
e8189707	33	#include "AliITS.h"
04366a57	34	#include "TParticle.h"
b0f5e3fc	35
b0f5e3fc	36	ClassImp(AliITSClusterFinder)
b0f5e3fc	37
9de0700b	38	//----------------------------------------------------------------------
aacedc3e	39	AliITSClusterFinder::AliITSClusterFinder():
	40	TObject(),
	41	fDebug(0),
	42	fModule(0),
	43	fDigits(0),
	44	fNdigits(0),
	45	fResponse(0),
	46	fSegmentation(0),
	47	fClusters(0),
	48	fNRawClusters(0),
	49	fMap(0),
	50	fNperMax(0),
	51	fDeclusterFlag(0),
	52	fClusterSize(0),
	53	fNPeaks(-1){
9de0700b	54	// default cluster finder
aacedc3e	55	// Input:
	56	// none.
	57	// Output:
	58	// none.
	59	// Return:
	60	// A default constructed AliITSCulsterFinder
04366a57	61	Init();
9de0700b	62	}
	63	//----------------------------------------------------------------------
	64	AliITSClusterFinder::AliITSClusterFinder(AliITSsegmentation *seg,
aacedc3e	65	AliITSresponse *res):
	66	TObject(),
	67	fDebug(0),
	68	fModule(0),
	69	fDigits(0),
	70	fNdigits(0),
	71	fResponse(res),
	72	fSegmentation(seg),
	73	fClusters(0),
	74	fNRawClusters(0),
	75	fMap(0),
	76	fNperMax(0),
	77	fDeclusterFlag(0),
	78	fClusterSize(0),
	79	fNPeaks(-1){
	80	// Standard constructor for cluster finder
	81	// Input:
	82	// AliITSsegmentation *seg The segmentation class to be used
	83	// AliITSresponse *res The response class to be used
	84	// Output:
	85	// none.
	86	// Return:
	87	// A Standard constructed AliITSCulsterFinder
b0f5e3fc	88
	89	SetNperMax();
	90	SetClusterSize();
	91	SetDeclusterFlag();
04366a57	92	Init();
aacedc3e	93	}
	94	//----------------------------------------------------------------------
	95	AliITSClusterFinder::AliITSClusterFinder(AliITSsegmentation *seg,
	96	AliITSresponse *response,
	97	TClonesArray *digits):
	98	TObject(),
	99	fDebug(0),
	100	fModule(0),
	101	fDigits(digits),
	102	fNdigits(0),
	103	fResponse(response),
	104	fSegmentation(seg),
	105	fClusters(0),
	106	fNRawClusters(0),
	107	fMap(0),
	108	fNperMax(0),
	109	fDeclusterFlag(0),
	110	fClusterSize(0),
	111	fNPeaks(-1){
	112	// Standard + cluster finder constructor
	113	// Input:
	114	// AliITSsegmentation *seg The segmentation class to be used
	115	// AliITSresponse *res The response class to be used
	116	// TClonesArray *digits Array of digits to be used
	117	// Output:
	118	// none.
	119	// Return:
	120	// A Standard constructed AliITSCulsterFinder
b0f5e3fc	121
aacedc3e	122	fNdigits = fDigits->GetEntriesFast();
	123	SetNperMax();
	124	SetClusterSize();
	125	SetDeclusterFlag();
04366a57	126	Init();
b0f5e3fc	127	}
04366a57	128	//______________________________________________________________________
	129	AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : TObject(source) {
	130	// Copy constructor
	131	// Copies are not allowed. The method is protected to avoid misuse.
	132	Fatal("AliITSClusterFinder","Copy constructor not allowed\n");
	133	}
	134
	135	//______________________________________________________________________
	136	AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
	137	// Assignment operator
	138	// Assignment is not allowed. The method is protected to avoid misuse.
	139	Fatal("= operator","Assignment operator not allowed\n");
	140	return *this;
	141	}
	142
9de0700b	143	//----------------------------------------------------------------------
	144	AliITSClusterFinder::~AliITSClusterFinder(){
	145	// destructor cluster finder
aacedc3e	146	// Input:
	147	// none.
	148	// Output:
	149	// none.
	150	// Return:
	151	// none.
b0f5e3fc	152
aacedc3e	153	if(fMap) {delete fMap;}
9de0700b	154	// Zero local pointers. Other classes own these pointers.
55cd883f	155	fSegmentation = 0;
9de0700b	156	fResponse = 0;
	157	fMap = 0;
	158	fDigits = 0;
	159	fNdigits = 0;
	160	fNRawClusters = 0;
	161	fNperMax = 0;
	162	fDeclusterFlag= 0;
	163	fClusterSize = 0;
	164	fNPeaks = 0;
04366a57	165	fITS = 0;
b0f5e3fc	166	}
04366a57	167
b0f5e3fc	168	//__________________________________________________________________________
04366a57	169	void AliITSClusterFinder::Init(){
	170
	171	//Initialisation of ITS geometry
	172
	173	fITS = (AliITS*)gAlice->GetModule("ITS");
	174	AliITSgeom geom=(AliITSgeom)fITS->GetITSgeom();
	175
	176	Int_t mmax=geom->GetIndexMax();
	177	if (mmax>2200) {
	178	Fatal("AliITSClusterFinder","Too many ITS subdetectors !");
	179	}
	180	Int_t m;
	181	for (m=0; m<mmax; m++) {
	182	Int_t lay,lad,det; geom->GetModuleId(m,lay,lad,det);
	183	Float_t x,y,z; geom->GetTrans(lay,lad,det,x,y,z);
	184	Double_t rot[9]; geom->GetRotMatrix(lay,lad,det,rot);
	185	Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
	186	Double_t ca=TMath::Cos(alpha), sa=TMath::Sin(alpha);
	187	fYshift[m] = xca + ysa;
	188	fZshift[m] = (Double_t)z;
	189	fNdet[m] = (lad-1)*geom->GetNdetectors(lay) + (det-1);
	190	fNlayer[m] = lay-1;
	191	}
aacedc3e	192
b0f5e3fc	193	}
04366a57	194
04366a57	195
9de0700b	196	//----------------------------------------------------------------------
	197	void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c){
	198	// Add a raw cluster copy to the list
aacedc3e	199	// Input:
	200	// Int_t branch The branch to which the cluster is to be added to
	201	// AliITSRawCluster *c The cluster to be added to the array of clusters
	202	// Output:
	203	// none.
	204	// Return:
	205	// none.
b0f5e3fc	206
04366a57	207	fITS->AddCluster(branch,c);
b0f5e3fc	208	fNRawClusters++;
9355b256	209	}
9de0700b	210	//----------------------------------------------------------------------
	211	void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c,
	212	AliITSRecPoint &rp){
aacedc3e	213	// Add a raw cluster copy to the list and the RecPoint
	214	// Input:
	215	// Int_t branch The branch to which the cluster is to be added to
	216	// AliITSRawCluster *c The cluster to be added to the array of clusters
	217	// AliITSRecPoint &rp The RecPoint to be added to the array of RecPoints
	218	// Output:
	219	// none.
	220	// Return:
	221	// none.
9355b256	222
04366a57	223	fITS->AddCluster(branch,c);
9355b256	224	fNRawClusters++;
04366a57	225	fITS->AddRecPoint(rp);
b0f5e3fc	226	}
04366a57	227
	228	//______________________________________________________________________
	229	void AliITSClusterFinder::CheckLabels(Int_t lab[3]) {
	230	//------------------------------------------------------------
	231	// Tries to find mother's labels
	232	//------------------------------------------------------------
	233
	234	if(lab[0]<0 && lab[1]<0 && lab[2]<0) return; // In case of no labels just exit
	235	// Check if simulation
	236	AliMC* mc = gAlice->GetMCApp();
	237	if(!mc)return;
	238
	239	Int_t ntracks = mc->GetNtrack();
	240	for (Int_t i=0;i<3;i++){
	241	Int_t label = lab[i];
	242	if (label>=0 && label<ntracks) {
	243	TParticle part=(TParticle)mc->Particle(label);
	244	if (part->P() < 0.005) {
	245	Int_t m=part->GetFirstMother();
	246	if (m<0) {
	247	continue;
	248	}
	249	if (part->GetStatusCode()>0) {
	250	continue;
	251	}
	252	lab[i]=m;
	253	}
	254	}
	255	}
	256
	257	}
	258
f8d9a5b8	259	//______________________________________________________________________
	260	void AliITSClusterFinder::FindRawClusters(Int_t module){
	261	// Default Cluster finder.
	262	// Input:
	263	// Int_t module Module number for which culster are to be found.
	264	// Output:
	265	// none.
	266	// Return:
	267	// none.
	268	const Int_t kelms = 10;
	269	Int_t ndigits = fDigits->GetEntriesFast();
	270	TObjArray *digs = new TObjArray(ndigits);
	271	TObjArray *clusts = new TObjArray(ndigits); // max # cluster
	272	TObjArray *clust0=0; // A spacific cluster of digits
	273	TObjArray *clust1=0; // A spacific cluster of digits
	274	AliITSdigit *dig=0; // locat pointer to a digit
f824982a	275	Int_t i=0,nc=0,j[4],k,k2=0;
f8d9a5b8	276
f8d9a5b8	277	// Copy all digits for this module into a local TObjArray.
aacedc3e	278	for(i=0;i<ndigits;i++) digs->AddAt(new AliITSdigit(*(GetDigit(i))),i);
f8d9a5b8	279	digs->Sort();
	280	// First digit is a cluster.
	281	i = 0;
	282	nc = 0;
	283	clusts->AddAt(new TObjArray(kelms),nc);
	284	clust0 = (TObjArray*)(clusts->At(nc));
	285	clust0->AddAtFree(digs->At(i)); // move owner ship from digs to clusts
	286	nc++;
	287	for(i=1;i<ndigits;i++){
aacedc3e	288	if(IsNeighbor(digs,i,j)){
	289	dig = (AliITSdigit*)(digs->At(j[0]));
	290	// Add to existing cluster. Find which cluster this digis
	291	for(k=0;k<nc;k++){
	292	clust0 = ((TObjArray*)(clusts->At(k)));
	293	if(clust0->IndexOf(dig)>=0) break;
	294	} // end for k
	295	if(k>=nc){
	296	Fatal("FindRawClusters","Digit not found as expected");
	297	} // end if
	298	if(j[1]>=0){
	299	dig = (AliITSdigit*)(digs->At(j[1]));
	300	// Add to existing cluster. Find which cluster this digis
	301	for(k2=0;k2<nc;k2++){
	302	clust1 = ((TObjArray*)(clusts->At(k2)));
	303	if(clust1->IndexOf(dig)>=0) break;
	304	} // end for k2
	305	if(k2>=nc){
	306	Fatal("FindRawClusters","Digit not found as expected");
	307	} // end if
	308	} // end if j[1]>=0
	309	// Found cluster with neighboring digits add this one to it.
	310	if(clust0==clust1){ // same cluster
	311	clust0->AddAtFree(digs->At(i));
	312	clust0 = 0; // finished with cluster. zero for safty
	313	clust1 = 0; // finished wit hcluster. zero for safty
	314	}else{ // two different clusters which need to be merged.
	315	clust0->AddAtFree(digs->At(i)); // Add digit to this cluster.
	316	for(k=0;k<clust1->GetEntriesFast();k++){
	317	// move clust1 into clust0
	318	//move digit to this cluster
	319	clust0->AddAtFree(clust1->At(k));
	320	clust1->AddAt(0,k); // zero this one
	321	} // end for k
	322	delete clust1;
	323	clusts->AddAt(0,k2); // zero array of clusters element clust1
	324	clust0 = 0; // finished with cluster. zero for safty
	325	clust1 = 0; // finished wit hcluster. zero for safty
	326	} // end if clust0==clust1
	327	}else{// New cluster
	328	clusts->AddAt(new TObjArray(kelms),nc);
	329	clust0 = ((TObjArray*)(clusts->At(nc)));
	330	// move owner ship from digs to clusts
	331	clust0->AddAtFree(digs->At(i));
	332	clust0 = 0; // finished with cluster. zero for safty
	333	nc++;
	334	} // End if IsNeighbor
f8d9a5b8	335	} // end for i
	336	// There are now nc clusters in clusts. Each element of clust is an
	337	// array of digits which are clustered together.
	338
	339	// For each cluster call detector specific CreateRecPoints
	340	for(i=0;i<nc;i++) CreateRecPoints((TObjArray*)(clusts->At(i)),module);
	341
	342	// clean up at the end.
	343	for(i=0;i<nc;i++){
aacedc3e	344	clust0 =(TObjArray*)(clusts->At(i));
	345	// Digits deleted below, so zero this TObjArray
	346	for(k=0;k<clust0->GetEntriesFast();k++) clust0->AddAt(0,k);
	347	delete clust0; // Delete this TObjArray
	348	clusts->AddAt(0,i); // Contents deleted above, so zero it.
f8d9a5b8	349	} // end for i
	350	delete clusts; // Delete this TObjArray/
	351	// Delete the digits then the TObjArray which containted them.
	352	for(i=0;i<ndigits;i++) delete ((AliITSdigit*)(digs->At(i)));
	353	delete digs;
	354	}
	355	//______________________________________________________________________
aacedc3e	356	Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
f8d9a5b8	357	// Locagical function which checks to see if digit i has a neighbor.
	358	// If so, then it returns kTRUE and its neighbor index j.
	359	// This routine checks if the digits are side by side or one before the
	360	// other. Requires that the array of digits be in proper order.
	361	// Returns kTRUE in the following cases.
	362	// ji 0j if kdiagonal j0 0i
	363	// 00 0i if kdiagonal 0i j0
	364	// Inputs:
	365	// TObjArray *digs Array to search for neighbors in
	366	// Int_t i Index of digit for which we are searching for
	367	// a neighbor of.
	368	// Output:
	369	// Int_t j[4] Index of one or more of the digits which is a
	370	// neighbor of digit a index i.
	371	// Return:
	372	// Bool_t kTRUE if a neighbor was found kFALSE otherwise.
	373	Int_t ix,jx,iz,jz,j;
	374	const Bool_t kdiagonal=kFALSE;
	375	Bool_t nei[4];
	376
	377	// No neighbors found if array empty.
	378	if(digs->GetEntriesFast()<=0) return kFALSE;
	379	// can not be a digit with first element or elements out or range
	380	if(i<=0 \|\| i>=digs->GetEntriesFast()) return kFALSE;
	381
	382	for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
	383	ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
	384	iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
	385	for(j=0;j<i;j++){
aacedc3e	386	jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
	387	jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
	388	if(jx+1==ix && jz ==iz){n[0] = j;nei[0] = kTRUE;}
	389	if(jx ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
	390	if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
	391	if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
f8d9a5b8	392	} // end for k
	393	if(nei[0]\|\|nei[1]) return kTRUE;
	394	if(kdiagonal&&(nei[2]\|\|nei[3])) return kTRUE;
	395	// no Neighbors found.
	396	return kFALSE;
	397	}
aacedc3e	398
	399	//______________________________________________________________________
	400	void AliITSClusterFinder::Print(ostream *os){
	401	//Standard output format for this class
	402	// Inputs:
	403	// ostream *os Output stream
	404	// Output:
	405	// ostream *os Output stream
	406	// Return:
	407	// none.
	408
	409	*os << fDebug<<",";
	410	*os << fModule<<",";
	411	*os << fNdigits<<",";
	412	*os << fNRawClusters<<",";
	413	*os << fNperMax<<",";
	414	*os << fDeclusterFlag<<",";
	415	*os << fClusterSize<<",";
	416	*os << fNPeaks<<endl;
	417	}
04366a57	418	//______________________________________________________________________
	419	void AliITSClusterFinder::RecPoints2Clusters
	420	(const TClonesArray points, Int_t idx, TClonesArray clusters) {
	421	//------------------------------------------------------------
	422	// Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS
	423	// subdetector indexed by idx
	424	//------------------------------------------------------------
	425	AliITSgeom* geom = (AliITSgeom*)fITS->GetITSgeom();
	426	Int_t lastSPD1=geom->GetModuleIndex(2,1,1)-1;
	427	TClonesArray &cl=*clusters;
	428	Int_t ncl=points->GetEntriesFast();
	429	for (Int_t i=0; i<ncl; i++) {
	430	AliITSRecPoint p = (AliITSRecPoint )points->UncheckedAt(i);
	431	Float_t lp[5];
	432	lp[0]=-(-p->GetX()+fYshift[idx]); if (idx<=lastSPD1) lp[0]*=-1; //SPD1
	433	lp[1]= -p->GetZ()+fZshift[idx];
	434	lp[2]=p->GetSigmaX2();
	435	lp[3]=p->GetSigmaZ2();
	436	lp[4]=p->GetQ()*36./23333.; //electrons -> ADC
	437	Int_t lab[4];
	438	lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
	439	lab[3]=fNdet[idx];
	440	CheckLabels(lab);
	441	Int_t dummy[3]={0,0,0};
	442	new (cl[i]) AliITSclusterV2(lab,lp, dummy);
	443	}
	444	}
	445
aacedc3e	446	//______________________________________________________________________
	447	void AliITSClusterFinder::Read(istream *is){
	448	//Standard input for this class
	449	// Inputs:
	450	// istream *is Input stream
	451	// Output:
	452	// istream *is Input stream
	453	// Return:
	454	// none.
	455
	456	*is >> fDebug;
	457	*is >> fModule;
	458	*is >> fNdigits;
	459	*is >> fNRawClusters;
	460	*is >> fNperMax;
	461	*is >> fDeclusterFlag;
	462	*is >> fClusterSize;
	463	*is >> fNPeaks;
	464	}
	465	//______________________________________________________________________
	466	ostream &operator<<(ostream &os,AliITSClusterFinder &source){
	467	// Standard output streaming function.
	468	// Inputs:
	469	// ostream *os Output stream
	470	// AliITSClusterFinder &source Class to be printed
	471	// Output:
	472	// ostream *os Output stream
	473	// Return:
	474	// none.
	475
	476	source.Print(&os);
	477	return os;
	478	}
	479	//______________________________________________________________________
	480	istream &operator>>(istream &is,AliITSClusterFinder &source){
	481	// Standard output streaming function.
	482	// Inputs:
	483	// istream *is Input stream
	484	// AliITSClusterFinder &source Class to be read in.
	485	// Output:
	486	// istream *is Input stream
	487	// Return:
	488	// none.
	489
	490	source.Read(&is);
	491	return is;
	492	}
04366a57	493	/*
	494	void AliITSClusterFinder::RecPoints2Clusters
	495	(const TClonesArray points, Int_t idx, TClonesArray clusters) {
	496	//------------------------------------------------------------
	497	// Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS
	498	// subdetector indexed by idx
	499	//------------------------------------------------------------
	500	TClonesArray &cl=*clusters;
	501	Int_t ncl=points->GetEntriesFast();
	502	for (Int_t i=0; i<ncl; i++) {
	503	AliITSRecPoint p = (AliITSRecPoint )points->UncheckedAt(i);
	504	Float_t lp[5];
	505	lp[0]=-(-p->GetX()+fYshift[idx]); if (idx<=fLastSPD1) lp[0]*=-1; //SPD1
	506	lp[1]= -p->GetZ()+fZshift[idx];
	507	lp[2]=p->GetSigmaX2();
	508	lp[3]=p->GetSigmaZ2();
	509	lp[4]=p->GetQ()*36./23333.; //electrons -> ADC
	510	Int_t lab[4];
	511	lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
	512	lab[3]=fNdet[idx];
	513	CheckLabels(lab);
	514	Int_t dummy[3]={0,0,0};
	515	new (cl[i]) AliITSclusterV2(lab,lp, dummy);
	516	}
	517	}
	518	*/