1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 ////////////////////////////////////////////////////////////////////////////
17 // Base Class used to find //
18 // the reconstructed points for ITS //
19 // See also AliITSClusterFinderSPD, AliITSClusterFinderSDD, //
20 // AliITSClusterFinderSDD AliITSClusterFinderV2 //
21 ////////////////////////////////////////////////////////////////////////////
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>
33 ClassImp(AliITSClusterFinder)
35 extern AliRun *gAlice;
37 //----------------------------------------------------------------------
38 AliITSClusterFinder::AliITSClusterFinder():
47 fNModules(AliITSgeomTGeo::GetNModules()),
49 // default cluster finder
55 // A default constructed AliITSCulsterFinder
57 //----------------------------------------------------------------------
58 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
67 fNModules(AliITSgeomTGeo::GetNModules()),
69 // default cluster finder
70 // Standard constructor for cluster finder
72 // AliITSsegmentation *seg The segmentation class to be used
73 // AliITSresponse *res The response class to be used
77 // A Standard constructed AliITSCulsterFinder
80 //----------------------------------------------------------------------
81 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
82 TClonesArray *digits):
91 fNModules(AliITSgeomTGeo::GetNModules()),
93 // default cluster finder
94 // Standard + cluster finder constructor
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
102 // A Standard constructed AliITSCulsterFinder
104 fNdigits = fDigits->GetEntriesFast();
107 //______________________________________________________________________
108 AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : TObject(source),
109 fModule(source.fModule),
111 fNdigits(source.fNdigits),
115 fNPeaks(source.fNPeaks),
116 fNModules(source.fNModules),
117 fEvent(source.fEvent) {
119 // Copies are not allowed. The method is protected to avoid misuse.
120 AliError("Copy constructor not allowed\n");
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");
132 //----------------------------------------------------------------------
133 AliITSClusterFinder::~AliITSClusterFinder(){
134 // destructor cluster finder
142 if(fMap) {delete fMap;}
143 // Zero local pointers. Other classes own these pointers.
151 //__________________________________________________________________________
152 void AliITSClusterFinder::InitGeometry(){
154 // Initialisation of ITS geometry
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);
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
177 // TObjArray *digs Array to search for neighbors in
178 // Int_t i Index of digit for which we are searching for
181 // Int_t j[4] Index of one or more of the digits which is a
182 // neighbor of digit a index i.
184 // Bool_t kTRUE if a neighbor was found kFALSE otherwise.
186 const Bool_t kdiagonal=kFALSE;
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;
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();
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;}
205 if(nei[0]||nei[1]) return kTRUE;
206 if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
207 // no Neighbors found.
211 //______________________________________________________________________
212 void AliITSClusterFinder::Print(ostream *os) const{
213 //Standard output format for this class
215 // ostream *os Output stream
217 // ostream *os Output stream
222 *os << fNdigits<<",";
223 *os << fNPeaks<<endl;
225 //______________________________________________________________________
226 void AliITSClusterFinder::Read(istream *is) {
227 //Standard input for this class
229 // istream *is Input stream
231 // istream *is Input stream
239 //______________________________________________________________________
240 ostream &operator<<(ostream &os,AliITSClusterFinder &source){
241 // Standard output streaming function.
243 // ostream *os Output stream
244 // AliITSClusterFinder &source Class to be printed
246 // ostream *os Output stream
253 //______________________________________________________________________
254 istream &operator>>(istream &is,AliITSClusterFinder &source){
255 // Standard output streaming function.
257 // istream *is Input stream
258 // AliITSClusterFinder &source Class to be read in.
260 // istream *is Input stream
267 //______________________________________________________________________
268 void AliITSClusterFinder::CheckLabels2(Int_t lab[10]) {
269 //------------------------------------------------------------
270 // Tries to find mother's labels
271 //------------------------------------------------------------
272 AliRunLoader *rl = AliRunLoader::Instance();
274 TTree *trK=(TTree*)rl->TreeK();
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
282 Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
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);
289 if (part->P() < 0.02) {
290 Int_t m=part->GetFirstMother();
294 if (part->GetStatusCode()>0) {
300 if (part->P() < 0.12 && nlabels>3) {
306 if ( (label>ntracks||label <0) && nlabels>3) {
313 for (Int_t i=0;i<10;i++){
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) {
327 //compress labels -- if multi-times the same
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++){
339 for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
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();
350 if(label<0) return; // In case of no label just exit
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;
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();
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);
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);
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;
406 //______________________________________________________________________
407 void AliITSClusterFinder::
408 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
409 //------------------------------------------------------------
411 //------------------------------------------------------------
413 if (IsMaximum(k,max,b)) {
414 idx[n]=k; msk[n]=(2<<n);
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);
424 //______________________________________________________________________
425 void AliITSClusterFinder::
426 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
427 //------------------------------------------------------------
429 //------------------------------------------------------------
430 UShort_t q=bins[k].GetQ();
432 bins[k].SetMask(bins[k].GetMask()|m);
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);
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;
454 c.SetY(c.GetY()+i*q);
455 c.SetZ(c.GetZ()+j*q);
456 c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
457 c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
459 bins[k].SetMask(0xFFFFFFFE);
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);