[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderV2.cxx

/**************************************************************************
 * Copyright(c) 1998-2003, 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.                  *
 **************************************************************************/
////////////////////////////////////////////////////////////////////////////
//            Implementation of the ITS clusterer V2 class                //
//                                                                        //
//          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch            //
//                                                                        //
///////////////////////////////////////////////////////////////////////////


#include "AliRun.h"
#include "AliITSClusterFinderV2.h"
#include "AliITSRecPoint.h"
#include "AliITSclusterV2.h"
#include "AliITSsimulationFastPoints.h"
#include "AliITS.h"
#include "AliITSgeom.h"
#include <TParticle.h>
#include "AliMC.h"

ClassImp(AliITSClusterFinderV2)

extern AliRun *gAlice;

AliITSClusterFinderV2::AliITSClusterFinderV2():AliITSClusterFinder(){

  //Default constructor
  fEvent = 0;
  fModule = 0;

  AliITSgeom *geom=(AliITSgeom*)fITS->GetITSgeom();

  fNModules = geom->GetIndexMax();
}

//______________________________________________________________________
AliITSClusterFinderV2::AliITSClusterFinderV2(const AliITSClusterFinderV2 &source) : AliITSClusterFinder(source) {
  // Copy constructor
  // Copies are not allowed. The method is protected to avoid misuse.
  Fatal("AliITSClusterFinderV2","Copy constructor not allowed\n");
}

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


//______________________________________________________________________
void AliITSClusterFinderV2::CheckLabels2(Int_t lab[10]) {
  //------------------------------------------------------------
  // Tries to find mother's labels
  //------------------------------------------------------------
  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 AliITSClusterFinderV2::AddLabel(Int_t lab[10], Int_t label) {

  //add label to the cluster

  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 AliITSClusterFinderV2:: 
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 AliITSClusterFinderV2::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 AliITSClusterFinderV2::
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 AliITSClusterFinderV2::
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 AliITSClusterFinderV2::
MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &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);
}

//______________________________________________________________________
Int_t AliITSClusterFinderV2::Hits2Clusters(TTree *hTree, TTree *cTree) {
  //------------------------------------------------------------
  // This function creates ITS clusters
  //------------------------------------------------------------

  AliITSgeom *geom=fITS->GetITSgeom();
  Int_t mmax=geom->GetIndexMax();

  fITS->InitModules(-1,mmax);
  fITS->FillModules(hTree,0);

  TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000);
  TBranch *branch=cTree->GetBranch("Clusters");
  if (!branch) cTree->Branch("Clusters",&clusters);
  else branch->SetAddress(&clusters);

  static TClonesArray *points=fITS->RecPoints();
  AliITSsimulationFastPoints sim;
  Int_t ncl=0;
  for (Int_t m=0; m<mmax; m++) {
    AliITSmodule *mod=fITS->GetModule(m);      
    sim.CreateFastRecPoints(mod,m,gRandom);      

    RecPoints2Clusters(points, m, clusters);
    fITS->ResetRecPoints();

    ncl+=clusters->GetEntriesFast();
    cTree->Fill();
    clusters->Clear();
  }

  Info("Hits2Clusters","Number of found fast clusters : %d",ncl);

  //cTree->Write();

  delete clusters;

  return 0;
}
Commit	Line	Data
04366a57	1	/**************************************************************************
	2	* Copyright(c) 1998-2003, 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	// Implementation of the ITS clusterer V2 class //
	17	// //
	18	// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch //
	19	// //
	20	///////////////////////////////////////////////////////////////////////////
	21
	22
	23	#include "AliRun.h"
	24	#include "AliITSClusterFinderV2.h"
	25	#include "AliITSRecPoint.h"
	26	#include "AliITSclusterV2.h"
	27	#include "AliITSsimulationFastPoints.h"
	28	#include "AliITS.h"
	29	#include "AliITSgeom.h"
	30	#include <TParticle.h>
	31	#include "AliMC.h"
	32
	33	ClassImp(AliITSClusterFinderV2)
	34
	35	extern AliRun *gAlice;
	36
	37	AliITSClusterFinderV2::AliITSClusterFinderV2():AliITSClusterFinder(){
	38
	39	//Default constructor
	40	fEvent = 0;
	41	fModule = 0;
	42
	43	AliITSgeom geom=(AliITSgeom)fITS->GetITSgeom();
	44
	45	fNModules = geom->GetIndexMax();
	46	}
	47
	48	//______________________________________________________________________
	49	AliITSClusterFinderV2::AliITSClusterFinderV2(const AliITSClusterFinderV2 &source) : AliITSClusterFinder(source) {
	50	// Copy constructor
	51	// Copies are not allowed. The method is protected to avoid misuse.
	52	Fatal("AliITSClusterFinderV2","Copy constructor not allowed\n");
	53	}
	54
	55	//______________________________________________________________________
	56	AliITSClusterFinderV2& AliITSClusterFinderV2::operator=(const AliITSClusterFinderV2& /* source */){
	57	// Assignment operator
	58	// Assignment is not allowed. The method is protected to avoid misuse.
	59	Fatal("= operator","Assignment operator not allowed\n");
	60	return *this;
	61	}
	62
	63
	64	//______________________________________________________________________
65	void AliITSClusterFinderV2::CheckLabels2(Int_t lab[10]) {
66	//------------------------------------------------------------
67	// Tries to find mother's labels
68	//------------------------------------------------------------
69	Int_t nlabels =0;
70	for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
71	if(nlabels == 0) return; // In case of no labels just exit
72
73
74	Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
75
76	for (Int_t i=0;i<10;i++){
77	Int_t label = lab[i];
78	if (label>=0 && label<ntracks) {
79	TParticle part=(TParticle)gAlice->GetMCApp()->Particle(label);
80
81	if (part->P() < 0.02) {
82	Int_t m=part->GetFirstMother();
83	if (m<0) {
84	continue;
85	}
86	if (part->GetStatusCode()>0) {
87	continue;
88	}
89	lab[i]=m;
90	}
91	else
92	if (part->P() < 0.12 && nlabels>3) {
93	lab[i]=-2;
94	nlabels--;
95	}
96	}
97	else{
98	if ( (label>ntracks\|\|label <0) && nlabels>3) {
99	lab[i]=-2;
100	nlabels--;
101	}
102	}
103	}
104	if (nlabels>3){
105	for (Int_t i=0;i<10;i++){
106	if (nlabels>3){
107	Int_t label = lab[i];
108	if (label>=0 && label<ntracks) {
109	TParticle part=(TParticle)gAlice->GetMCApp()->Particle(label);
110	if (part->P() < 0.1) {
111	lab[i]=-2;
112	nlabels--;
113	}
114	}
115	}
116	}
117	}
118
119	//compress labels -- if multi-times the same
120	Int_t lab2[10];
121	for (Int_t i=0;i<10;i++) lab2[i]=-2;
122	for (Int_t i=0;i<10 ;i++){
123	if (lab[i]<0) continue;
124	for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
125	if (lab2[j]<0) {
126	lab2[j]= lab[i];
127	break;
128	}
129	}
130	}
131	for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
132
133	}
134
135	//______________________________________________________________________
136	void AliITSClusterFinderV2::AddLabel(Int_t lab[10], Int_t label) {
137
138	//add label to the cluster
139
140	if(label<0) return; // In case of no label just exit
141
142	Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
143	if (label>ntracks) return;
144	for (Int_t i=0;i<10;i++){
145	// if (label<0) break;
146	if (lab[i]==label) break;
147	if (lab[i]<0) {
148	lab[i]= label;
149	break;
150	}
151	}
152	}
153
154
155	//______________________________________________________________________
156	void AliITSClusterFinderV2::
157	FindCluster(Int_t k,Int_t maxz,AliBin bins,Int_t &n,Int_t idx) {
158	//------------------------------------------------------------
159	// returns an array of indices of digits belonging to the cluster
160	// (needed when the segmentation is not regular)
161	//------------------------------------------------------------
162	if (n<200) idx[n++]=bins[k].GetIndex();
163	bins[k].Use();
164
165	if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
166	if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx);
167	if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
168	if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx);
169	/*
170	if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
171	if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
172	if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
173	if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
174	*/
175	}
176
177	//______________________________________________________________________
178	Bool_t AliITSClusterFinderV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
179	//------------------------------------------------------------
180	//is this a local maximum ?
181	//------------------------------------------------------------
182	UShort_t q=bins[k].GetQ();
183	if (q==1023) return kFALSE;
184	if (bins[k-max].GetQ() > q) return kFALSE;
185	if (bins[k-1 ].GetQ() > q) return kFALSE;
186	if (bins[k+max].GetQ() > q) return kFALSE;
187	if (bins[k+1 ].GetQ() > q) return kFALSE;
188	if (bins[k-max-1].GetQ() > q) return kFALSE;
189	if (bins[k+max-1].GetQ() > q) return kFALSE;
190	if (bins[k+max+1].GetQ() > q) return kFALSE;
191	if (bins[k-max+1].GetQ() > q) return kFALSE;
192	return kTRUE;
193	}
194
195	//______________________________________________________________________
196	void AliITSClusterFinderV2::
197	FindPeaks(Int_t k,Int_t max,AliBin b,Int_t idx,UInt_t *msk,Int_t& n) {
198	//------------------------------------------------------------
199	//find local maxima
200	//------------------------------------------------------------
201	if (n<31)
202	if (IsMaximum(k,max,b)) {
203	idx[n]=k; msk[n]=(2<<n);
204	n++;
205	}
206	b[k].SetMask(0);
207	if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
208	if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
209	if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
210	if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
211	}
212
213	//______________________________________________________________________
214	void AliITSClusterFinderV2::
215	MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
216	//------------------------------------------------------------
217	//mark this peak
218	//------------------------------------------------------------
219	UShort_t q=bins[k].GetQ();
220
221	bins[k].SetMask(bins[k].GetMask()\|m);
222
223	if (bins[k-max].GetQ() <= q)
224	if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
225	if (bins[k-1 ].GetQ() <= q)
226	if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
227	if (bins[k+max].GetQ() <= q)
228	if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
229	if (bins[k+1 ].GetQ() <= q)
230	if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
231	}
232
233	//______________________________________________________________________
234	void AliITSClusterFinderV2::
235	MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) {
236	//------------------------------------------------------------
237	//make cluster using digits of this peak
238	//------------------------------------------------------------
239	Float_t q=(Float_t)bins[k].GetQ();
240	Int_t i=k/max, j=k-i*max;
241
242	c.SetQ(c.GetQ()+q);
243	c.SetY(c.GetY()+i*q);
244	c.SetZ(c.GetZ()+j*q);
245	c.SetSigmaY2(c.GetSigmaY2()+iiq);
246	c.SetSigmaZ2(c.GetSigmaZ2()+jjq);
247
248	bins[k].SetMask(0xFFFFFFFE);
249
250	if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
251	if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
252	if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
253	if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);
254	}
255
256	//______________________________________________________________________
257	Int_t AliITSClusterFinderV2::Hits2Clusters(TTree hTree, TTree cTree) {
258	//------------------------------------------------------------
259	// This function creates ITS clusters
260	//------------------------------------------------------------
261
262	AliITSgeom *geom=fITS->GetITSgeom();
263	Int_t mmax=geom->GetIndexMax();
264
265	fITS->InitModules(-1,mmax);
266	fITS->FillModules(hTree,0);
267
268	TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000);
269	TBranch *branch=cTree->GetBranch("Clusters");
270	if (!branch) cTree->Branch("Clusters",&clusters);
271	else branch->SetAddress(&clusters);
272
273	static TClonesArray *points=fITS->RecPoints();
274	AliITSsimulationFastPoints sim;
275	Int_t ncl=0;
276	for (Int_t m=0; m<mmax; m++) {
277	AliITSmodule *mod=fITS->GetModule(m);
278	sim.CreateFastRecPoints(mod,m,gRandom);
279
280	RecPoints2Clusters(points, m, clusters);
281	fITS->ResetRecPoints();
282
283	ncl+=clusters->GetEntriesFast();
284	cTree->Fill();
285	clusters->Clear();
286	}
287
288	Info("Hits2Clusters","Number of found fast clusters : %d",ncl);
289
290	//cTree->Write();
291
292	delete clusters;
293
294	return 0;
295	}
296