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>
34 ClassImp(AliITSClusterFinder)
36 extern AliRun *gAlice;
38 //----------------------------------------------------------------------
39 AliITSClusterFinder::AliITSClusterFinder():
48 fNModules(AliITSgeomTGeo::GetNModules()),
54 // default cluster finder
60 // A default constructed AliITSCulsterFinder
61 for(Int_t i=0; i<2200; i++){
66 //----------------------------------------------------------------------
67 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
76 fNModules(AliITSgeomTGeo::GetNModules()),
82 // default cluster finder
83 // Standard constructor for cluster finder
85 // AliITSsegmentation *seg The segmentation class to be used
86 // AliITSresponse *res The response class to be used
90 // A Standard constructed AliITSCulsterFinder
91 for(Int_t i=0; i<2200; i++){
96 //----------------------------------------------------------------------
97 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
98 TClonesArray *digits):
107 fNModules(AliITSgeomTGeo::GetNModules()),
113 // default cluster finder
114 // Standard + cluster finder constructor
116 // AliITSsegmentation *seg The segmentation class to be used
117 // AliITSresponse *res The response class to be used
118 // TClonesArray *digits Array of digits to be used
122 // A Standard constructed AliITSCulsterFinder
124 fNdigits = fDigits->GetEntriesFast();
125 for(Int_t i=0; i<2200; i++){
131 //______________________________________________________________________
132 AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) :
134 fModule(source.fModule),
136 fNdigits(source.fNdigits),
140 fNPeaks(source.fNPeaks),
141 fNModules(source.fNModules),
142 fEvent(source.fEvent),
149 // Copies are not allowed. The method is protected to avoid misuse.
150 AliError("Copy constructor not allowed\n");
154 //______________________________________________________________________
155 //AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
156 // Assignment operator
157 // Assignment is not allowed. The method is protected to avoid misuse.
158 // Fatal("= operator","Assignment operator not allowed\n");
162 //----------------------------------------------------------------------
163 AliITSClusterFinder::~AliITSClusterFinder(){
164 // destructor cluster finder
172 if(fMap) {delete fMap;}
173 // 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);
197 //______________________________________________________________________
198 Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
199 // Locagical function which checks to see if digit i has a neighbor.
200 // If so, then it returns kTRUE and its neighbor index j.
201 // This routine checks if the digits are side by side or one before the
202 // other. Requires that the array of digits be in proper order.
203 // Returns kTRUE in the following cases.
204 // ji 0j if kdiagonal j0 0i
205 // 00 0i if kdiagonal 0i j0
207 // TObjArray *digs Array to search for neighbors in
208 // Int_t i Index of digit for which we are searching for
211 // Int_t j[4] Index of one or more of the digits which is a
212 // neighbor of digit a index i.
214 // Bool_t kTRUE if a neighbor was found kFALSE otherwise.
216 const Bool_t kdiagonal=kFALSE;
219 // No neighbors found if array empty.
220 if(digs->GetEntriesFast()<=0) return kFALSE;
221 // can not be a digit with first element or elements out or range
222 if(i<=0 || i>=digs->GetEntriesFast()) return kFALSE;
224 for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
225 ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
226 iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
228 jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
229 jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
230 if(jx+1==ix && jz ==iz){n[0] = j;nei[0] = kTRUE;}
231 if(jx ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
232 if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
233 if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
235 if(nei[0]||nei[1]) return kTRUE;
236 if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
237 // no Neighbors found.
241 //______________________________________________________________________
242 void AliITSClusterFinder::Print(ostream *os) const{
243 //Standard output format for this class
245 // ostream *os Output stream
247 // ostream *os Output stream
252 *os << fNdigits<<",";
253 *os << fNPeaks<<endl;
255 //______________________________________________________________________
256 void AliITSClusterFinder::Read(istream *is) {
257 //Standard input for this class
259 // istream *is Input stream
261 // istream *is Input stream
269 //______________________________________________________________________
270 ostream &operator<<(ostream &os,AliITSClusterFinder &source){
271 // Standard output streaming function.
273 // ostream *os Output stream
274 // AliITSClusterFinder &source Class to be printed
276 // ostream *os Output stream
283 //______________________________________________________________________
284 istream &operator>>(istream &is,AliITSClusterFinder &source){
285 // Standard output streaming function.
287 // istream *is Input stream
288 // AliITSClusterFinder &source Class to be read in.
290 // istream *is Input stream
298 //______________________________________________________________________
299 void AliITSClusterFinder::CheckLabels2(Int_t lab[10])
301 //------------------------------------------------------------
302 // Tries to find mother's labels
303 //------------------------------------------------------------
304 AliRunLoader *rl = AliRunLoader::Instance();
306 TTree *trK=(TTree*)rl->TreeK();
311 Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
312 for (Int_t i=0;i<10;i++) if (lab[i]>=0) labS[nlabels++] = lab[i];
313 if (nlabels==0) return;
316 for (Int_t i=0;i<nlabels;i++) {
317 Int_t label = labS[i];
319 if (label>=ntracks) continue;
320 TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
322 if (part->P() < 0.02) { // reduce soft particles from the same cluster
323 Int_t m=part->GetFirstMother();
324 if (m<0) continue; // primary
326 if (part->GetStatusCode()>0) continue;
328 // if the parent is within the same cluster, reassign the label to it
329 for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break; }
333 if (nlabels>3) { // only 3 labels are stored in cluster, sort in decreasing momentum
334 int ind[10],labSS[10];
335 TMath::Sort(nlabels,mom,ind);
336 for (int i=nlabels;i--;) labSS[i] = labS[i];
337 for (int i=0;i<nlabels;i++) labS[i] = labSS[ind[i]];
340 //compress labels -- if multi-times the same
341 for (Int_t i=0;i<10;i++) lab[i]=-2;
343 for (int i=0;i<nlabels;i++) {
344 for (j=0;j<nlabFin;j++) if (labS[i]==lab[j]) break; // the label already there
345 if (j==nlabFin) lab[nlabFin++] = labS[i];
350 //______________________________________________________________________
351 void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
352 //add label to the cluster
353 AliRunLoader *rl = AliRunLoader::Instance();
354 TTree *trK=(TTree*)rl->TreeK();
356 if(label<0) return; // In case of no label just exit
358 Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
359 if (label>ntracks) return;
360 for (Int_t i=0;i<10;i++){
361 // if (label<0) break;
362 if (lab[i]==label) break;
372 //______________________________________________________________________
373 void AliITSClusterFinder::
374 FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
375 //------------------------------------------------------------
376 // returns an array of indices of digits belonging to the cluster
377 // (needed when the segmentation is not regular)
378 //------------------------------------------------------------
379 if (n<200) idx[n++]=bins[k].GetIndex();
382 if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
383 if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx);
384 if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
385 if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx);
387 if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
388 if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
389 if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
390 if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
394 //______________________________________________________________________
395 Bool_t AliITSClusterFinder::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
396 //------------------------------------------------------------
397 //is this a local maximum ?
398 //------------------------------------------------------------
399 UShort_t q=bins[k].GetQ();
400 if (q==1023) return kFALSE;
401 if (bins[k-max].GetQ() > q) return kFALSE;
402 if (bins[k-1 ].GetQ() > q) return kFALSE;
403 if (bins[k+max].GetQ() > q) return kFALSE;
404 if (bins[k+1 ].GetQ() > q) return kFALSE;
405 if (bins[k-max-1].GetQ() > q) return kFALSE;
406 if (bins[k+max-1].GetQ() > q) return kFALSE;
407 if (bins[k+max+1].GetQ() > q) return kFALSE;
408 if (bins[k-max+1].GetQ() > q) return kFALSE;
412 //______________________________________________________________________
413 void AliITSClusterFinder::
414 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
415 //------------------------------------------------------------
417 //------------------------------------------------------------
419 if (IsMaximum(k,max,b)) {
420 idx[n]=k; msk[n]=(2<<n);
424 if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
425 if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
426 if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
427 if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
430 //______________________________________________________________________
431 void AliITSClusterFinder::
432 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
433 //------------------------------------------------------------
435 //------------------------------------------------------------
436 UShort_t q=bins[k].GetQ();
438 bins[k].SetMask(bins[k].GetMask()|m);
440 if (bins[k-max].GetQ() <= q)
441 if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
442 if (bins[k-1 ].GetQ() <= q)
443 if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
444 if (bins[k+max].GetQ() <= q)
445 if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
446 if (bins[k+1 ].GetQ() <= q)
447 if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
450 //______________________________________________________________________
451 void AliITSClusterFinder::
452 MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSRecPoint &c) {
453 //------------------------------------------------------------
454 //make cluster using digits of this peak
455 //------------------------------------------------------------
456 Float_t q=(Float_t)bins[k].GetQ();
457 Int_t i=k/max, j=k-i*max;
458 if(c.GetQ()<0.01){ // first entry in cluster
463 }else{ // check cluster extension
470 c.SetY(c.GetY()+i*q);
471 c.SetZ(c.GetZ()+j*q);
472 c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
473 c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
475 bins[k].SetMask(0xFFFFFFFE);
477 if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
478 if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
479 if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
480 if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);