Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, 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 : // //
17 : // Base Class used to find //
18 : // the reconstructed points for ITS //
19 : // See also AliITSClusterFinderSPD, AliITSClusterFinderSDD, //
20 : // AliITSClusterFinderSDD AliITSClusterFinderV2 //
21 : ////////////////////////////////////////////////////////////////////////////
22 :
23 : #include "AliRun.h"
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>
31 : #include <TArrayI.h>
32 : #include "AliMC.h"
33 : #include "AliLog.h"
34 : #include <iostream>
35 :
36 : using std::endl;
37 :
38 118 : ClassImp(AliITSClusterFinder)
39 :
40 : extern AliRun *gAlice;
41 :
42 : //----------------------------------------------------------------------
43 : AliITSClusterFinder::AliITSClusterFinder():
44 0 : TObject(),
45 0 : fModule(0),
46 0 : fDigits(0),
47 0 : fNdigits(0),
48 0 : fDetTypeRec(0),
49 0 : fClusters(0),
50 0 : fMap(0),
51 0 : fNPeaks(-1),
52 0 : fNModules(AliITSgeomTGeo::GetNModules()),
53 0 : fEvent(0),
54 0 : fZmin(0),
55 0 : fZmax(0),
56 0 : fXmin(0),
57 0 : fXmax(0),
58 0 : fNClusters(0),
59 0 : fRawID2ClusID(0)
60 0 : {
61 : // default cluster finder
62 : // Input:
63 : // none.
64 : // Output:
65 : // none.
66 : // Return:
67 : // A default constructed AliITSCulsterFinder
68 0 : for(Int_t i=0; i<2200; i++){
69 0 : fNdet[i]=0;
70 0 : fNlayer[i]=0;
71 : }
72 0 : }
73 : //----------------------------------------------------------------------
74 : AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
75 6 : TObject(),
76 6 : fModule(0),
77 6 : fDigits(0),
78 6 : fNdigits(0),
79 6 : fDetTypeRec(dettyp),
80 6 : fClusters(0),
81 6 : fMap(0),
82 6 : fNPeaks(-1),
83 6 : fNModules(AliITSgeomTGeo::GetNModules()),
84 6 : fEvent(0),
85 6 : fZmin(0),
86 6 : fZmax(0),
87 6 : fXmin(0),
88 6 : fXmax(0),
89 6 : fNClusters(0),
90 6 : fRawID2ClusID(0)
91 18 : {
92 : // default cluster finder
93 : // Standard constructor for cluster finder
94 : // Input:
95 : // AliITSsegmentation *seg The segmentation class to be used
96 : // AliITSresponse *res The response class to be used
97 : // Output:
98 : // none.
99 : // Return:
100 : // A Standard constructed AliITSCulsterFinder
101 26412 : for(Int_t i=0; i<2200; i++){
102 13200 : fNdet[i]=0;
103 13200 : fNlayer[i]=0;
104 : }
105 6 : }
106 : //----------------------------------------------------------------------
107 : AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
108 : TClonesArray *digits):
109 0 : TObject(),
110 0 : fModule(0),
111 0 : fDigits(digits),
112 0 : fNdigits(0),
113 0 : fDetTypeRec(dettyp),
114 0 : fClusters(0),
115 0 : fMap(0),
116 0 : fNPeaks(-1),
117 0 : fNModules(AliITSgeomTGeo::GetNModules()),
118 0 : fEvent(0),
119 0 : fZmin(0),
120 0 : fZmax(0),
121 0 : fXmin(0),
122 0 : fXmax(0),
123 0 : fNClusters(0),
124 0 : fRawID2ClusID(0)
125 0 : {
126 : // default cluster finder
127 : // Standard + cluster finder constructor
128 : // Input:
129 : // AliITSsegmentation *seg The segmentation class to be used
130 : // AliITSresponse *res The response class to be used
131 : // TClonesArray *digits Array of digits to be used
132 : // Output:
133 : // none.
134 : // Return:
135 : // A Standard constructed AliITSCulsterFinder
136 :
137 0 : fNdigits = fDigits->GetEntriesFast();
138 0 : for(Int_t i=0; i<2200; i++){
139 0 : fNdet[i]=0;
140 0 : fNlayer[i]=0;
141 : }
142 0 : }
143 :
144 : //______________________________________________________________________
145 : AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) :
146 0 : TObject(source),
147 0 : fModule(source.fModule),
148 0 : fDigits(),
149 0 : fNdigits(source.fNdigits),
150 0 : fDetTypeRec(),
151 0 : fClusters(),
152 0 : fMap(),
153 0 : fNPeaks(source.fNPeaks),
154 0 : fNModules(source.fNModules),
155 0 : fEvent(source.fEvent),
156 0 : fZmin(source.fZmin),
157 0 : fZmax(source.fZmax),
158 0 : fXmin(source.fXmin),
159 0 : fXmax(source.fXmax),
160 0 : fNClusters(source.fNClusters),
161 0 : fRawID2ClusID(source.fRawID2ClusID)
162 0 : {
163 : // Copy constructor
164 : // Copies are not allowed. The method is protected to avoid misuse.
165 0 : AliError("Copy constructor not allowed\n");
166 0 : }
167 :
168 :
169 : //______________________________________________________________________
170 : //AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
171 : // Assignment operator
172 : // Assignment is not allowed. The method is protected to avoid misuse.
173 : // Fatal("= operator","Assignment operator not allowed\n");
174 : // return *this;
175 : //}
176 :
177 : //----------------------------------------------------------------------
178 12 : AliITSClusterFinder::~AliITSClusterFinder(){
179 : // destructor cluster finder
180 : // Input:
181 : // none.
182 : // Output:
183 : // none.
184 : // Return:
185 : // none.
186 :
187 6 : if(fMap) {delete fMap;}
188 : // Zero local pointers. Other classes own these pointers.
189 6 : fMap = 0;
190 6 : fDigits = 0;
191 6 : fNdigits = 0;
192 6 : fNPeaks = 0;
193 6 : fDetTypeRec = 0;
194 :
195 6 : }
196 : //__________________________________________________________________________
197 : void AliITSClusterFinder::InitGeometry(){
198 : //
199 : // Initialisation of ITS geometry
200 : //
201 12 : Int_t mmax=AliITSgeomTGeo::GetNModules();
202 26388 : for (Int_t m=0; m<mmax; m++) {
203 13188 : Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
204 13188 : fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
205 13188 : fNlayer[m] = lay-1;
206 13188 : }
207 6 : }
208 :
209 :
210 :
211 :
212 : //______________________________________________________________________
213 : Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
214 : // Locagical function which checks to see if digit i has a neighbor.
215 : // If so, then it returns kTRUE and its neighbor index j.
216 : // This routine checks if the digits are side by side or one before the
217 : // other. Requires that the array of digits be in proper order.
218 : // Returns kTRUE in the following cases.
219 : // ji 0j if kdiagonal j0 0i
220 : // 00 0i if kdiagonal 0i j0
221 : // Inputs:
222 : // TObjArray *digs Array to search for neighbors in
223 : // Int_t i Index of digit for which we are searching for
224 : // a neighbor of.
225 : // Output:
226 : // Int_t j[4] Index of one or more of the digits which is a
227 : // neighbor of digit a index i.
228 : // Return:
229 : // Bool_t kTRUE if a neighbor was found kFALSE otherwise.
230 : Int_t ix,jx,iz,jz,j;
231 : const Bool_t kdiagonal=kFALSE;
232 0 : Bool_t nei[4];
233 :
234 : // No neighbors found if array empty.
235 0 : if(digs->GetEntriesFast()<=0) return kFALSE;
236 : // can not be a digit with first element or elements out or range
237 0 : if(i<=0 || i>=digs->GetEntriesFast()) return kFALSE;
238 :
239 0 : for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
240 0 : ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
241 0 : iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
242 0 : for(j=0;j<i;j++){
243 0 : jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
244 0 : jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
245 0 : if(jx+1==ix && jz ==iz){n[0] = j;nei[0] = kTRUE;}
246 0 : if(jx ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
247 0 : if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
248 0 : if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
249 : } // end for k
250 0 : if(nei[0]||nei[1]) return kTRUE;
251 : if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
252 : // no Neighbors found.
253 0 : return kFALSE;
254 0 : }
255 :
256 : //______________________________________________________________________
257 : void AliITSClusterFinder::Print(ostream *os) const{
258 : //Standard output format for this class
259 : // Inputs:
260 : // ostream *os Output stream
261 : // Output:
262 : // ostream *os Output stream
263 : // Return:
264 : // none.
265 :
266 0 : *os << fModule<<",";
267 0 : *os << fNdigits<<",";
268 0 : *os << fNPeaks<<endl;
269 0 : }
270 : //______________________________________________________________________
271 : void AliITSClusterFinder::Read(istream *is) {
272 : //Standard input for this class
273 : // Inputs:
274 : // istream *is Input stream
275 : // Output:
276 : // istream *is Input stream
277 : // Return:
278 : // none.
279 :
280 0 : *is >> fModule;
281 0 : *is >> fNdigits;
282 0 : *is >> fNPeaks;
283 0 : }
284 : //______________________________________________________________________
285 : ostream &operator<<(ostream &os,AliITSClusterFinder &source){
286 : // Standard output streaming function.
287 : // Inputs:
288 : // ostream *os Output stream
289 : // AliITSClusterFinder &source Class to be printed
290 : // Output:
291 : // ostream *os Output stream
292 : // Return:
293 : // none.
294 :
295 0 : source.Print(&os);
296 0 : return os;
297 : }
298 : //______________________________________________________________________
299 : istream &operator>>(istream &is,AliITSClusterFinder &source){
300 : // Standard output streaming function.
301 : // Inputs:
302 : // istream *is Input stream
303 : // AliITSClusterFinder &source Class to be read in.
304 : // Output:
305 : // istream *is Input stream
306 : // Return:
307 : // none.
308 :
309 0 : source.Read(&is);
310 0 : return is;
311 : }
312 :
313 : //______________________________________________________________________
314 : void AliITSClusterFinder::CheckLabels2(Int_t lab[10])
315 : {
316 : //------------------------------------------------------------
317 : // Tries to find mother's labels
318 : //------------------------------------------------------------
319 1398 : AliRunLoader *rl = AliRunLoader::Instance();
320 699 : if(!rl) return;
321 699 : TTree *trK=(TTree*)rl->TreeK();
322 845 : if (!trK) return;
323 : //
324 553 : int labS[10];
325 : Int_t nlabels = 0;
326 553 : Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
327 18545 : for (Int_t i=0;i<10;i++) if (lab[i]>=0) labS[nlabels++] = lab[i];
328 557 : if (nlabels==0) return;
329 : //
330 549 : float mom[10];
331 2796 : for (Int_t i=0;i<nlabels;i++) {
332 849 : Int_t label = labS[i];
333 849 : mom[i] = 0;
334 849 : if (label>=ntracks) continue;
335 849 : TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
336 849 : mom[i] = part->P();
337 849 : if (part->P() < 0.02) { // reduce soft particles from the same cluster
338 162 : Int_t m=part->GetFirstMother();
339 162 : if (m<0) continue; // primary
340 : //
341 162 : if (part->GetStatusCode()>0) continue;
342 : //
343 : // if the parent is within the same cluster, reassign the label to it
344 1107 : for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break; }
345 162 : }
346 849 : }
347 : //
348 549 : if (nlabels>3) { // only 3 labels are stored in cluster, sort in decreasing momentum
349 10 : int ind[10],labSS[10];
350 10 : TMath::Sort(nlabels,mom,ind);
351 102 : for (int i=nlabels;i--;) labSS[i] = labS[i];
352 102 : for (int i=0;i<nlabels;i++) labS[i] = labSS[ind[i]];
353 10 : }
354 : //
355 : //compress labels -- if multi-times the same
356 12078 : for (Int_t i=0;i<10;i++) lab[i]=-2;
357 : int nlabFin=0,j=0;
358 2796 : for (int i=0;i<nlabels;i++) {
359 2084 : for (j=0;j<nlabFin;j++) if (labS[i]==lab[j]) break; // the label already there
360 1432 : if (j==nlabFin) lab[nlabFin++] = labS[i];
361 : }
362 : //
363 1801 : }
364 :
365 : //______________________________________________________________________
366 : void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
367 : //add label to the cluster
368 2552 : AliRunLoader *rl = AliRunLoader::Instance();
369 1276 : TTree *trK=(TTree*)rl->TreeK();
370 1276 : if(trK){
371 1276 : if(label<0) return; // In case of no label just exit
372 :
373 1276 : Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
374 1276 : if (label>ntracks) return;
375 3444 : for (Int_t i=0;i<10;i++){
376 : // if (label<0) break;
377 2484 : if (lab[i]==label) break;
378 960 : if (lab[i]<0) {
379 514 : lab[i]= label;
380 514 : break;
381 : }
382 : }
383 1276 : }
384 2552 : }
385 :
386 :
387 : //______________________________________________________________________
388 : void AliITSClusterFinder::
389 : FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
390 : //------------------------------------------------------------
391 : // returns an array of indices of digits belonging to the cluster
392 : // (needed when the segmentation is not regular)
393 : //------------------------------------------------------------
394 822 : if (n<200) idx[n++]=bins[k].GetIndex();
395 274 : bins[k].Use();
396 :
397 278 : if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
398 274 : if (bins[k-1 ].IsNotUsed()) FindCluster(k-1 ,maxz,bins,n,idx);
399 384 : if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
400 288 : if (bins[k+1 ].IsNotUsed()) FindCluster(k+1 ,maxz,bins,n,idx);
401 : /*
402 : if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
403 : if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
404 : if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
405 : if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
406 : */
407 274 : }
408 :
409 : //______________________________________________________________________
410 : Bool_t AliITSClusterFinder::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
411 : //------------------------------------------------------------
412 : //is this a local maximum ?
413 : //------------------------------------------------------------
414 0 : UShort_t q=bins[k].GetQ();
415 0 : if (q==1023) return kFALSE;
416 0 : if (bins[k-max].GetQ() > q) return kFALSE;
417 0 : if (bins[k-1 ].GetQ() > q) return kFALSE;
418 0 : if (bins[k+max].GetQ() > q) return kFALSE;
419 0 : if (bins[k+1 ].GetQ() > q) return kFALSE;
420 0 : if (bins[k-max-1].GetQ() > q) return kFALSE;
421 0 : if (bins[k+max-1].GetQ() > q) return kFALSE;
422 0 : if (bins[k+max+1].GetQ() > q) return kFALSE;
423 0 : if (bins[k-max+1].GetQ() > q) return kFALSE;
424 0 : return kTRUE;
425 0 : }
426 :
427 : //______________________________________________________________________
428 : void AliITSClusterFinder::
429 : FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
430 : //------------------------------------------------------------
431 : //find local maxima
432 : //------------------------------------------------------------
433 0 : if (n<31)
434 0 : if (IsMaximum(k,max,b)) {
435 0 : idx[n]=k; msk[n]=(2<<n);
436 0 : n++;
437 0 : }
438 0 : b[k].SetMask(0);
439 0 : if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
440 0 : if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
441 0 : if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
442 0 : if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
443 0 : }
444 :
445 : //______________________________________________________________________
446 : void AliITSClusterFinder::
447 : MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
448 : //------------------------------------------------------------
449 : //mark this peak
450 : //------------------------------------------------------------
451 0 : UShort_t q=bins[k].GetQ();
452 :
453 0 : bins[k].SetMask(bins[k].GetMask()|m);
454 :
455 0 : if (bins[k-max].GetQ() <= q)
456 0 : if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
457 0 : if (bins[k-1 ].GetQ() <= q)
458 0 : if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
459 0 : if (bins[k+max].GetQ() <= q)
460 0 : if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
461 0 : if (bins[k+1 ].GetQ() <= q)
462 0 : if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
463 0 : }
464 :
465 : //______________________________________________________________________
466 : void AliITSClusterFinder::
467 : MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSRecPoint &c) {
468 : //------------------------------------------------------------
469 : //make cluster using digits of this peak
470 : //------------------------------------------------------------
471 0 : Float_t q=(Float_t)bins[k].GetQ();
472 0 : Int_t i=k/max, j=k-i*max;
473 0 : if(c.GetQ()<0.01){ // first entry in cluster
474 0 : fXmin=i;
475 0 : fXmax=i;
476 0 : fZmin=j;
477 0 : fZmax=j;
478 0 : }else{ // check cluster extension
479 0 : if(i<fXmin) fXmin=i;
480 0 : if(i>fXmax) fXmax=i;
481 0 : if(j<fZmin) fZmin=j;
482 0 : if(j>fZmax) fZmax=j;
483 : }
484 0 : c.SetQ(c.GetQ()+q);
485 0 : c.SetY(c.GetY()+i*q);
486 0 : c.SetZ(c.GetZ()+j*q);
487 0 : c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
488 0 : c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
489 :
490 0 : bins[k].SetMask(0xFFFFFFFE);
491 0 : if (fRawID2ClusID) { // RS: Register cluster id in raw words list
492 0 : int rwid = bins[k].GetRawID();
493 0 : if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
494 0 : (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
495 0 : }
496 0 : if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
497 0 : if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
498 0 : if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
499 0 : if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);
500 0 : }
|
