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():
52 fNModules(AliITSgeomTGeo::GetNModules()),
54 // default cluster finder
60 // A default constructed AliITSCulsterFinder
62 //----------------------------------------------------------------------
63 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
77 fNModules(AliITSgeomTGeo::GetNModules()),
79 // default cluster finder
80 // Standard constructor for cluster finder
82 // AliITSsegmentation *seg The segmentation class to be used
83 // AliITSresponse *res The response class to be used
87 // A Standard constructed AliITSCulsterFinder
93 //----------------------------------------------------------------------
94 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
95 TClonesArray *digits):
109 fNModules(AliITSgeomTGeo::GetNModules()),
111 // default cluster finder
112 // Standard + cluster finder constructor
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
120 // A Standard constructed AliITSCulsterFinder
122 fNdigits = fDigits->GetEntriesFast();
128 //______________________________________________________________________
129 AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : TObject(source),
130 fDebug(source.fDebug),
131 fModule(source.fModule),
133 fNdigits(source.fNdigits),
136 fNRawClusters(source.fNRawClusters),
138 fNperMax(source.fNperMax),
139 fDeclusterFlag(source.fDeclusterFlag),
140 fClusterSize(source.fClusterSize),
141 fNPeaks(source.fNPeaks),
142 fNModules(source.fNModules),
143 fEvent(source.fEvent) {
145 // Copies are not allowed. The method is protected to avoid misuse.
146 AliError("Copy constructor not allowed\n");
150 //______________________________________________________________________
151 //AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
152 // Assignment operator
153 // Assignment is not allowed. The method is protected to avoid misuse.
154 // Fatal("= operator","Assignment operator not allowed\n");
158 //----------------------------------------------------------------------
159 AliITSClusterFinder::~AliITSClusterFinder(){
160 // destructor cluster finder
168 if(fMap) {delete fMap;}
169 // Zero local pointers. Other classes own these pointers.
181 //__________________________________________________________________________
182 void AliITSClusterFinder::InitGeometry(){
184 // Initialisation of ITS geometry
186 Int_t mmax=AliITSgeomTGeo::GetNModules();
187 for (Int_t m=0; m<mmax; m++) {
188 Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
189 fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
196 //----------------------------------------------------------------------
197 void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c){
198 // Add a raw cluster copy to the list
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
208 Error("AddCluster","fDetTypeRec is null!");
211 fDetTypeRec->AddCluster(branch,c);
214 //----------------------------------------------------------------------
215 void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c,
217 // Add a raw cluster copy to the list and the RecPoint
219 // Int_t branch The branch to which the cluster is to be added to
220 // AliITSRawCluster *c The cluster to be added to the array of clusters
221 // AliITSRecPoint &rp The RecPoint to be added to the array of RecPoints
227 Error("AddCluster","fDetTypeRec is null!");
231 fDetTypeRec->AddCluster(branch,c);
233 fDetTypeRec->AddRecPoint(rp);
237 //______________________________________________________________________
238 void AliITSClusterFinder::CheckLabels(Int_t lab[3]) {
239 //------------------------------------------------------------
240 // Tries to find mother's labels
241 //------------------------------------------------------------
243 if(lab[0]<0 && lab[1]<0 && lab[2]<0) return; // In case of no labels just exit
244 // Check if simulation
245 AliMC* mc = gAlice->GetMCApp();
248 Int_t ntracks = mc->GetNtrack();
249 for (Int_t i=0;i<3;i++){
250 Int_t label = lab[i];
251 if (label>=0 && label<ntracks) {
252 TParticle *part=(TParticle*)mc->Particle(label);
253 if (part->P() < 0.005) {
254 Int_t m=part->GetFirstMother();
258 if (part->GetStatusCode()>0) {
268 //______________________________________________________________________
269 void AliITSClusterFinder::FindRawClusters(Int_t module){
270 // Default Cluster finder.
272 // Int_t module Module number for which culster are to be found.
277 const Int_t kelms = 10;
278 Int_t ndigits = fDigits->GetEntriesFast();
279 TObjArray *digs = new TObjArray(ndigits);
280 TObjArray *clusts = new TObjArray(ndigits); // max # cluster
281 TObjArray *clust0=0; // A spacific cluster of digits
282 TObjArray *clust1=0; // A spacific cluster of digits
283 AliITSdigit *dig=0; // locat pointer to a digit
284 Int_t i=0,nc=0,j[4],k,k2=0;
286 // Copy all digits for this module into a local TObjArray.
287 for(i=0;i<ndigits;i++) digs->AddAt(new AliITSdigit(*(GetDigit(i))),i);
289 // First digit is a cluster.
292 clusts->AddAt(new TObjArray(kelms),nc);
293 clust0 = (TObjArray*)(clusts->At(nc));
294 clust0->AddAtFree(digs->At(i)); // move owner ship from digs to clusts
296 for(i=1;i<ndigits;i++){
297 if(IsNeighbor(digs,i,j)){
298 dig = (AliITSdigit*)(digs->At(j[0]));
299 // Add to existing cluster. Find which cluster this digis
301 clust0 = ((TObjArray*)(clusts->At(k)));
302 if(clust0->IndexOf(dig)>=0) break;
305 Fatal("FindRawClusters","Digit not found as expected");
308 dig = (AliITSdigit*)(digs->At(j[1]));
309 // Add to existing cluster. Find which cluster this digis
310 for(k2=0;k2<nc;k2++){
311 clust1 = ((TObjArray*)(clusts->At(k2)));
312 if(clust1->IndexOf(dig)>=0) break;
315 Fatal("FindRawClusters","Digit not found as expected");
318 // Found cluster with neighboring digits add this one to it.
319 if(clust0==clust1){ // same cluster
320 clust0->AddAtFree(digs->At(i));
321 clust0 = 0; // finished with cluster. zero for safty
322 clust1 = 0; // finished wit hcluster. zero for safty
323 }else{ // two different clusters which need to be merged.
324 clust0->AddAtFree(digs->At(i)); // Add digit to this cluster.
325 for(k=0;k<clust1->GetEntriesFast();k++){
326 // move clust1 into clust0
327 //move digit to this cluster
328 clust0->AddAtFree(clust1->At(k));
329 clust1->AddAt(0,k); // zero this one
332 clusts->AddAt(0,k2); // zero array of clusters element clust1
333 clust0 = 0; // finished with cluster. zero for safty
334 clust1 = 0; // finished wit hcluster. zero for safty
335 } // end if clust0==clust1
337 clusts->AddAt(new TObjArray(kelms),nc);
338 clust0 = ((TObjArray*)(clusts->At(nc)));
339 // move owner ship from digs to clusts
340 clust0->AddAtFree(digs->At(i));
341 clust0 = 0; // finished with cluster. zero for safty
343 } // End if IsNeighbor
345 // There are now nc clusters in clusts. Each element of clust is an
346 // array of digits which are clustered together.
348 // For each cluster call detector specific CreateRecPoints
349 for(i=0;i<nc;i++) CreateRecPoints((TObjArray*)(clusts->At(i)),module);
351 // clean up at the end.
353 clust0 =(TObjArray*)(clusts->At(i));
354 // Digits deleted below, so zero this TObjArray
355 for(k=0;k<clust0->GetEntriesFast();k++) clust0->AddAt(0,k);
356 delete clust0; // Delete this TObjArray
357 clusts->AddAt(0,i); // Contents deleted above, so zero it.
359 delete clusts; // Delete this TObjArray/
360 // Delete the digits then the TObjArray which containted them.
361 for(i=0;i<ndigits;i++) delete ((AliITSdigit*)(digs->At(i)));
364 //______________________________________________________________________
365 Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
366 // Locagical function which checks to see if digit i has a neighbor.
367 // If so, then it returns kTRUE and its neighbor index j.
368 // This routine checks if the digits are side by side or one before the
369 // other. Requires that the array of digits be in proper order.
370 // Returns kTRUE in the following cases.
371 // ji 0j if kdiagonal j0 0i
372 // 00 0i if kdiagonal 0i j0
374 // TObjArray *digs Array to search for neighbors in
375 // Int_t i Index of digit for which we are searching for
378 // Int_t j[4] Index of one or more of the digits which is a
379 // neighbor of digit a index i.
381 // Bool_t kTRUE if a neighbor was found kFALSE otherwise.
383 const Bool_t kdiagonal=kFALSE;
386 // No neighbors found if array empty.
387 if(digs->GetEntriesFast()<=0) return kFALSE;
388 // can not be a digit with first element or elements out or range
389 if(i<=0 || i>=digs->GetEntriesFast()) return kFALSE;
391 for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
392 ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
393 iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
395 jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
396 jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
397 if(jx+1==ix && jz ==iz){n[0] = j;nei[0] = kTRUE;}
398 if(jx ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
399 if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
400 if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
402 if(nei[0]||nei[1]) return kTRUE;
403 if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
404 // no Neighbors found.
408 //______________________________________________________________________
409 void AliITSClusterFinder::Print(ostream *os) const{
410 //Standard output format for this class
412 // ostream *os Output stream
414 // ostream *os Output stream
420 *os << fNdigits<<",";
421 *os << fNRawClusters<<",";
422 *os << fNperMax<<",";
423 *os << fDeclusterFlag<<",";
424 *os << fClusterSize<<",";
425 *os << fNPeaks<<endl;
427 //______________________________________________________________________
428 void AliITSClusterFinder::Read(istream *is) {
429 //Standard input for this class
431 // istream *is Input stream
433 // istream *is Input stream
440 *is >> fNRawClusters;
442 *is >> fDeclusterFlag;
446 //______________________________________________________________________
447 ostream &operator<<(ostream &os,AliITSClusterFinder &source){
448 // Standard output streaming function.
450 // ostream *os Output stream
451 // AliITSClusterFinder &source Class to be printed
453 // ostream *os Output stream
460 //______________________________________________________________________
461 istream &operator>>(istream &is,AliITSClusterFinder &source){
462 // Standard output streaming function.
464 // istream *is Input stream
465 // AliITSClusterFinder &source Class to be read in.
467 // istream *is Input stream
474 //______________________________________________________________________
475 void AliITSClusterFinder::CheckLabels2(Int_t lab[10]) {
476 //------------------------------------------------------------
477 // Tries to find mother's labels
478 //------------------------------------------------------------
479 AliRunLoader *rl = AliRunLoader::GetRunLoader();
480 TTree *trK=(TTree*)rl->TreeK();
484 for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
485 if(nlabels == 0) return; // In case of no labels just exit
488 Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
490 for (Int_t i=0;i<10;i++){
491 Int_t label = lab[i];
492 if (label>=0 && label<ntracks) {
493 TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
495 if (part->P() < 0.02) {
496 Int_t m=part->GetFirstMother();
500 if (part->GetStatusCode()>0) {
506 if (part->P() < 0.12 && nlabels>3) {
512 if ( (label>ntracks||label <0) && nlabels>3) {
519 for (Int_t i=0;i<10;i++){
521 Int_t label = lab[i];
522 if (label>=0 && label<ntracks) {
523 TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
524 if (part->P() < 0.1) {
533 //compress labels -- if multi-times the same
535 for (Int_t i=0;i<10;i++) lab2[i]=-2;
536 for (Int_t i=0;i<10 ;i++){
537 if (lab[i]<0) continue;
538 for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
545 for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
550 //______________________________________________________________________
551 void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
552 //add label to the cluster
553 AliRunLoader *rl = AliRunLoader::GetRunLoader();
554 TTree *trK=(TTree*)rl->TreeK();
556 if(label<0) return; // In case of no label just exit
558 Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
559 if (label>ntracks) return;
560 for (Int_t i=0;i<10;i++){
561 // if (label<0) break;
562 if (lab[i]==label) break;
572 //______________________________________________________________________
573 void AliITSClusterFinder::
574 FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
575 //------------------------------------------------------------
576 // returns an array of indices of digits belonging to the cluster
577 // (needed when the segmentation is not regular)
578 //------------------------------------------------------------
579 if (n<200) idx[n++]=bins[k].GetIndex();
582 if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
583 if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx);
584 if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
585 if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx);
587 if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
588 if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
589 if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
590 if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
594 //______________________________________________________________________
595 Bool_t AliITSClusterFinder::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
596 //------------------------------------------------------------
597 //is this a local maximum ?
598 //------------------------------------------------------------
599 UShort_t q=bins[k].GetQ();
600 if (q==1023) return kFALSE;
601 if (bins[k-max].GetQ() > q) return kFALSE;
602 if (bins[k-1 ].GetQ() > q) return kFALSE;
603 if (bins[k+max].GetQ() > q) return kFALSE;
604 if (bins[k+1 ].GetQ() > q) return kFALSE;
605 if (bins[k-max-1].GetQ() > q) return kFALSE;
606 if (bins[k+max-1].GetQ() > q) return kFALSE;
607 if (bins[k+max+1].GetQ() > q) return kFALSE;
608 if (bins[k-max+1].GetQ() > q) return kFALSE;
612 //______________________________________________________________________
613 void AliITSClusterFinder::
614 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
615 //------------------------------------------------------------
617 //------------------------------------------------------------
619 if (IsMaximum(k,max,b)) {
620 idx[n]=k; msk[n]=(2<<n);
624 if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
625 if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
626 if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
627 if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
630 //______________________________________________________________________
631 void AliITSClusterFinder::
632 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
633 //------------------------------------------------------------
635 //------------------------------------------------------------
636 UShort_t q=bins[k].GetQ();
638 bins[k].SetMask(bins[k].GetMask()|m);
640 if (bins[k-max].GetQ() <= q)
641 if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
642 if (bins[k-1 ].GetQ() <= q)
643 if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
644 if (bins[k+max].GetQ() <= q)
645 if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
646 if (bins[k+1 ].GetQ() <= q)
647 if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
650 //______________________________________________________________________
651 void AliITSClusterFinder::
652 MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSRecPoint &c) {
653 //------------------------------------------------------------
654 //make cluster using digits of this peak
655 //------------------------------------------------------------
656 Float_t q=(Float_t)bins[k].GetQ();
657 Int_t i=k/max, j=k-i*max;
660 c.SetY(c.GetY()+i*q);
661 c.SetZ(c.GetZ()+j*q);
662 c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
663 c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
665 bins[k].SetMask(0xFFFFFFFE);
667 if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
668 if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
669 if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
670 if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);