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>
35 ClassImp(AliITSClusterFinder)
37 extern AliRun *gAlice;
39 //----------------------------------------------------------------------
40 AliITSClusterFinder::AliITSClusterFinder():
49 fNModules(AliITSgeomTGeo::GetNModules()),
58 // default cluster finder
64 // A default constructed AliITSCulsterFinder
65 for(Int_t i=0; i<2200; i++){
70 //----------------------------------------------------------------------
71 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
80 fNModules(AliITSgeomTGeo::GetNModules()),
89 // default cluster finder
90 // Standard constructor for cluster finder
92 // AliITSsegmentation *seg The segmentation class to be used
93 // AliITSresponse *res The response class to be used
97 // A Standard constructed AliITSCulsterFinder
98 for(Int_t i=0; i<2200; i++){
103 //----------------------------------------------------------------------
104 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
105 TClonesArray *digits):
114 fNModules(AliITSgeomTGeo::GetNModules()),
123 // default cluster finder
124 // Standard + cluster finder constructor
126 // AliITSsegmentation *seg The segmentation class to be used
127 // AliITSresponse *res The response class to be used
128 // TClonesArray *digits Array of digits to be used
132 // A Standard constructed AliITSCulsterFinder
134 fNdigits = fDigits->GetEntriesFast();
135 for(Int_t i=0; i<2200; i++){
141 //______________________________________________________________________
142 AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) :
144 fModule(source.fModule),
146 fNdigits(source.fNdigits),
150 fNPeaks(source.fNPeaks),
151 fNModules(source.fNModules),
152 fEvent(source.fEvent),
157 fNClusters(source.fNClusters),
158 fRawID2ClusID(source.fRawID2ClusID)
161 // Copies are not allowed. The method is protected to avoid misuse.
162 AliError("Copy constructor not allowed\n");
166 //______________________________________________________________________
167 //AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
168 // Assignment operator
169 // Assignment is not allowed. The method is protected to avoid misuse.
170 // Fatal("= operator","Assignment operator not allowed\n");
174 //----------------------------------------------------------------------
175 AliITSClusterFinder::~AliITSClusterFinder(){
176 // destructor cluster finder
184 if(fMap) {delete fMap;}
185 // Zero local pointers. Other classes own these pointers.
193 //__________________________________________________________________________
194 void AliITSClusterFinder::InitGeometry(){
196 // Initialisation of ITS geometry
198 Int_t mmax=AliITSgeomTGeo::GetNModules();
199 for (Int_t m=0; m<mmax; m++) {
200 Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
201 fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
209 //______________________________________________________________________
210 Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
211 // Locagical function which checks to see if digit i has a neighbor.
212 // If so, then it returns kTRUE and its neighbor index j.
213 // This routine checks if the digits are side by side or one before the
214 // other. Requires that the array of digits be in proper order.
215 // Returns kTRUE in the following cases.
216 // ji 0j if kdiagonal j0 0i
217 // 00 0i if kdiagonal 0i j0
219 // TObjArray *digs Array to search for neighbors in
220 // Int_t i Index of digit for which we are searching for
223 // Int_t j[4] Index of one or more of the digits which is a
224 // neighbor of digit a index i.
226 // Bool_t kTRUE if a neighbor was found kFALSE otherwise.
228 const Bool_t kdiagonal=kFALSE;
231 // No neighbors found if array empty.
232 if(digs->GetEntriesFast()<=0) return kFALSE;
233 // can not be a digit with first element or elements out or range
234 if(i<=0 || i>=digs->GetEntriesFast()) return kFALSE;
236 for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
237 ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
238 iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
240 jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
241 jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
242 if(jx+1==ix && jz ==iz){n[0] = j;nei[0] = kTRUE;}
243 if(jx ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
244 if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
245 if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
247 if(nei[0]||nei[1]) return kTRUE;
248 if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
249 // no Neighbors found.
253 //______________________________________________________________________
254 void AliITSClusterFinder::Print(ostream *os) const{
255 //Standard output format for this class
257 // ostream *os Output stream
259 // ostream *os Output stream
264 *os << fNdigits<<",";
265 *os << fNPeaks<<endl;
267 //______________________________________________________________________
268 void AliITSClusterFinder::Read(istream *is) {
269 //Standard input for this class
271 // istream *is Input stream
273 // istream *is Input stream
281 //______________________________________________________________________
282 ostream &operator<<(ostream &os,AliITSClusterFinder &source){
283 // Standard output streaming function.
285 // ostream *os Output stream
286 // AliITSClusterFinder &source Class to be printed
288 // ostream *os Output stream
295 //______________________________________________________________________
296 istream &operator>>(istream &is,AliITSClusterFinder &source){
297 // Standard output streaming function.
299 // istream *is Input stream
300 // AliITSClusterFinder &source Class to be read in.
302 // istream *is Input stream
310 //______________________________________________________________________
311 void AliITSClusterFinder::CheckLabels2(Int_t lab[10])
313 //------------------------------------------------------------
314 // Tries to find mother's labels
315 //------------------------------------------------------------
316 AliRunLoader *rl = AliRunLoader::Instance();
318 TTree *trK=(TTree*)rl->TreeK();
323 Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
324 for (Int_t i=0;i<10;i++) if (lab[i]>=0) labS[nlabels++] = lab[i];
325 if (nlabels==0) return;
328 for (Int_t i=0;i<nlabels;i++) {
329 Int_t label = labS[i];
331 if (label>=ntracks) continue;
332 TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
334 if (part->P() < 0.02) { // reduce soft particles from the same cluster
335 Int_t m=part->GetFirstMother();
336 if (m<0) continue; // primary
338 if (part->GetStatusCode()>0) continue;
340 // if the parent is within the same cluster, reassign the label to it
341 for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break; }
345 if (nlabels>3) { // only 3 labels are stored in cluster, sort in decreasing momentum
346 int ind[10],labSS[10];
347 TMath::Sort(nlabels,mom,ind);
348 for (int i=nlabels;i--;) labSS[i] = labS[i];
349 for (int i=0;i<nlabels;i++) labS[i] = labSS[ind[i]];
352 //compress labels -- if multi-times the same
353 for (Int_t i=0;i<10;i++) lab[i]=-2;
355 for (int i=0;i<nlabels;i++) {
356 for (j=0;j<nlabFin;j++) if (labS[i]==lab[j]) break; // the label already there
357 if (j==nlabFin) lab[nlabFin++] = labS[i];
362 //______________________________________________________________________
363 void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
364 //add label to the cluster
365 AliRunLoader *rl = AliRunLoader::Instance();
366 TTree *trK=(TTree*)rl->TreeK();
368 if(label<0) return; // In case of no label just exit
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;
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();
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);
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);
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;
424 //______________________________________________________________________
425 void AliITSClusterFinder::
426 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
427 //------------------------------------------------------------
429 //------------------------------------------------------------
431 if (IsMaximum(k,max,b)) {
432 idx[n]=k; msk[n]=(2<<n);
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);
442 //______________________________________________________________________
443 void AliITSClusterFinder::
444 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
445 //------------------------------------------------------------
447 //------------------------------------------------------------
448 UShort_t q=bins[k].GetQ();
450 bins[k].SetMask(bins[k].GetMask()|m);
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);
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;
470 if(c.GetQ()<0.01){ // first entry in cluster
475 }else{ // check cluster extension
482 c.SetY(c.GetY()+i*q);
483 c.SetZ(c.GetZ()+j*q);
484 c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
485 c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
487 bins[k].SetMask(0xFFFFFFFE);
488 if (fRawID2ClusID) { // RS: Register cluster id in raw words list
489 int rwid = bins[k].GetRawID();
490 if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
491 (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
493 if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
494 if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
495 if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
496 if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);