Added two missing includes to allow macro compilation (thanks to Laurent for remarkin...
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinder.cxx
CommitLineData
b0f5e3fc 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 **************************************************************************/
04366a57 15////////////////////////////////////////////////////////////////////////////
16// //
17// Base Class used to find //
18// the reconstructed points for ITS //
19// See also AliITSClusterFinderSPD, AliITSClusterFinderSDD, //
20// AliITSClusterFinderSDD AliITSClusterFinderV2 //
21////////////////////////////////////////////////////////////////////////////
88cb7938 22
5d2c2f86 23#include "AliRun.h"
b0f5e3fc 24#include "AliITSClusterFinder.h"
04366a57 25#include "AliITSRecPoint.h"
7d62fb64 26#include "AliITSdigit.h"
27#include "AliITSDetTypeRec.h"
aacedc3e 28#include "AliITSMap.h"
1f3e997f 29#include "AliITSgeomTGeo.h"
5d2c2f86 30#include <TParticle.h>
31#include "AliMC.h"
b0f5e3fc 32
b0f5e3fc 33ClassImp(AliITSClusterFinder)
34
5d2c2f86 35extern AliRun *gAlice;
36
9de0700b 37//----------------------------------------------------------------------
aacedc3e 38AliITSClusterFinder::AliITSClusterFinder():
39TObject(),
aacedc3e 40fModule(0),
41fDigits(0),
42fNdigits(0),
8ba39da9 43fDetTypeRec(0),
aacedc3e 44fClusters(0),
aacedc3e 45fMap(0),
5d2c2f86 46fNPeaks(-1),
47fNModules(AliITSgeomTGeo::GetNModules()),
48fEvent(0){
9de0700b 49 // default cluster finder
aacedc3e 50 // Input:
51 // none.
52 // Output:
53 // none.
54 // Return:
55 // A default constructed AliITSCulsterFinder
9de0700b 56}
57//----------------------------------------------------------------------
8ba39da9 58AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
aacedc3e 59TObject(),
aacedc3e 60fModule(0),
61fDigits(0),
62fNdigits(0),
e56160b8 63fDetTypeRec(dettyp),
aacedc3e 64fClusters(0),
aacedc3e 65fMap(0),
5d2c2f86 66fNPeaks(-1),
67fNModules(AliITSgeomTGeo::GetNModules()),
68fEvent(0){
69 // default cluster finder
aacedc3e 70 // Standard constructor for cluster finder
71 // Input:
72 // AliITSsegmentation *seg The segmentation class to be used
73 // AliITSresponse *res The response class to be used
74 // Output:
75 // none.
76 // Return:
77 // A Standard constructed AliITSCulsterFinder
b0f5e3fc 78
aacedc3e 79}
80//----------------------------------------------------------------------
8ba39da9 81AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
82 TClonesArray *digits):
aacedc3e 83TObject(),
aacedc3e 84fModule(0),
85fDigits(digits),
86fNdigits(0),
8ba39da9 87fDetTypeRec(dettyp),
aacedc3e 88fClusters(0),
aacedc3e 89fMap(0),
5d2c2f86 90fNPeaks(-1),
91fNModules(AliITSgeomTGeo::GetNModules()),
92fEvent(0){
93 // default cluster finder
aacedc3e 94 // Standard + cluster finder constructor
95 // Input:
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
99 // Output:
100 // none.
101 // Return:
102 // A Standard constructed AliITSCulsterFinder
b0f5e3fc 103
aacedc3e 104 fNdigits = fDigits->GetEntriesFast();
b0f5e3fc 105}
a971ed8e 106
04366a57 107//______________________________________________________________________
a971ed8e 108AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : TObject(source),
a971ed8e 109fModule(source.fModule),
110fDigits(),
111fNdigits(source.fNdigits),
112fDetTypeRec(),
113fClusters(),
a971ed8e 114fMap(),
5d2c2f86 115fNPeaks(source.fNPeaks),
116fNModules(source.fNModules),
117fEvent(source.fEvent) {
04366a57 118 // Copy constructor
119 // Copies are not allowed. The method is protected to avoid misuse.
a971ed8e 120 AliError("Copy constructor not allowed\n");
04366a57 121}
a971ed8e 122
04366a57 123
124//______________________________________________________________________
e56160b8 125//AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
04366a57 126 // Assignment operator
127 // Assignment is not allowed. The method is protected to avoid misuse.
e56160b8 128// Fatal("= operator","Assignment operator not allowed\n");
129// return *this;
130//}
04366a57 131
9de0700b 132//----------------------------------------------------------------------
133AliITSClusterFinder::~AliITSClusterFinder(){
134 // destructor cluster finder
aacedc3e 135 // Input:
136 // none.
137 // Output:
138 // none.
139 // Return:
140 // none.
b0f5e3fc 141
aacedc3e 142 if(fMap) {delete fMap;}
9de0700b 143 // Zero local pointers. Other classes own these pointers.
9de0700b 144 fMap = 0;
145 fDigits = 0;
146 fNdigits = 0;
9de0700b 147 fNPeaks = 0;
7d62fb64 148 fDetTypeRec = 0;
04366a57 149
7d62fb64 150}
b0f5e3fc 151//__________________________________________________________________________
7d62fb64 152void AliITSClusterFinder::InitGeometry(){
1f3e997f 153 //
154 // Initialisation of ITS geometry
155 //
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);
1f3e997f 159 fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
04366a57 160 fNlayer[m] = lay-1;
161 }
b0f5e3fc 162}
04366a57 163
164
7d62fb64 165
04366a57 166
f8d9a5b8 167//______________________________________________________________________
aacedc3e 168Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
f8d9a5b8 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
176 // Inputs:
177 // TObjArray *digs Array to search for neighbors in
178 // Int_t i Index of digit for which we are searching for
179 // a neighbor of.
180 // Output:
181 // Int_t j[4] Index of one or more of the digits which is a
182 // neighbor of digit a index i.
183 // Return:
184 // Bool_t kTRUE if a neighbor was found kFALSE otherwise.
185 Int_t ix,jx,iz,jz,j;
186 const Bool_t kdiagonal=kFALSE;
187 Bool_t nei[4];
188
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;
193
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();
197 for(j=0;j<i;j++){
aacedc3e 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;}
f8d9a5b8 204 } // end for k
205 if(nei[0]||nei[1]) return kTRUE;
206 if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
207 // no Neighbors found.
208 return kFALSE;
209}
aacedc3e 210
211//______________________________________________________________________
7d62fb64 212void AliITSClusterFinder::Print(ostream *os) const{
aacedc3e 213 //Standard output format for this class
214 // Inputs:
215 // ostream *os Output stream
216 // Output:
217 // ostream *os Output stream
218 // Return:
219 // none.
220
aacedc3e 221 *os << fModule<<",";
222 *os << fNdigits<<",";
aacedc3e 223 *os << fNPeaks<<endl;
224}
04366a57 225//______________________________________________________________________
7d62fb64 226void AliITSClusterFinder::Read(istream *is) {
aacedc3e 227 //Standard input for this class
228 // Inputs:
229 // istream *is Input stream
230 // Output:
231 // istream *is Input stream
232 // Return:
233 // none.
234
aacedc3e 235 *is >> fModule;
236 *is >> fNdigits;
aacedc3e 237 *is >> fNPeaks;
238}
239//______________________________________________________________________
240ostream &operator<<(ostream &os,AliITSClusterFinder &source){
241 // Standard output streaming function.
242 // Inputs:
243 // ostream *os Output stream
244 // AliITSClusterFinder &source Class to be printed
245 // Output:
246 // ostream *os Output stream
247 // Return:
248 // none.
249
250 source.Print(&os);
251 return os;
252}
253//______________________________________________________________________
254istream &operator>>(istream &is,AliITSClusterFinder &source){
255 // Standard output streaming function.
256 // Inputs:
257 // istream *is Input stream
258 // AliITSClusterFinder &source Class to be read in.
259 // Output:
260 // istream *is Input stream
261 // Return:
262 // none.
263
264 source.Read(&is);
265 return is;
266}
5d2c2f86 267//______________________________________________________________________
268void AliITSClusterFinder::CheckLabels2(Int_t lab[10]) {
269 //------------------------------------------------------------
270 // Tries to find mother's labels
271 //------------------------------------------------------------
33c3c91a 272 AliRunLoader *rl = AliRunLoader::Instance();
c605e5f7 273 if(!rl) return;
5d2c2f86 274 TTree *trK=(TTree*)rl->TreeK();
275
276 if(trK){
277 Int_t nlabels =0;
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
280
281
282 Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
283
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);
288
289 if (part->P() < 0.02) {
290 Int_t m=part->GetFirstMother();
291 if (m<0) {
292 continue;
293 }
294 if (part->GetStatusCode()>0) {
295 continue;
296 }
297 lab[i]=m;
298 }
299 else
300 if (part->P() < 0.12 && nlabels>3) {
301 lab[i]=-2;
302 nlabels--;
303 }
304 }
305 else{
306 if ( (label>ntracks||label <0) && nlabels>3) {
307 lab[i]=-2;
308 nlabels--;
309 }
310 }
311 }
312 if (nlabels>3){
313 for (Int_t i=0;i<10;i++){
314 if (nlabels>3){
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) {
319 lab[i]=-2;
320 nlabels--;
321 }
322 }
323 }
324 }
325 }
326
327 //compress labels -- if multi-times the same
328 Int_t lab2[10];
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++){
333 if (lab2[j]<0) {
334 lab2[j]= lab[i];
335 break;
336 }
337 }
338 }
339 for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
340
341 }
342}
343
344//______________________________________________________________________
345void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
346 //add label to the cluster
33c3c91a 347 AliRunLoader *rl = AliRunLoader::Instance();
5d2c2f86 348 TTree *trK=(TTree*)rl->TreeK();
349 if(trK){
350 if(label<0) return; // In case of no label just exit
351
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;
357 if (lab[i]<0) {
358 lab[i]= label;
359 break;
360 }
361 }
362 }
363}
364
365
366//______________________________________________________________________
367void AliITSClusterFinder::
368FindCluster(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();
374 bins[k].Use();
375
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);
380 /*
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);
385 */
386}
387
388//______________________________________________________________________
389Bool_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;
403 return kTRUE;
404}
405
406//______________________________________________________________________
407void AliITSClusterFinder::
408FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
409 //------------------------------------------------------------
410 //find local maxima
411 //------------------------------------------------------------
412 if (n<31)
413 if (IsMaximum(k,max,b)) {
414 idx[n]=k; msk[n]=(2<<n);
415 n++;
416 }
417 b[k].SetMask(0);
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);
422}
423
424//______________________________________________________________________
425void AliITSClusterFinder::
426MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
427 //------------------------------------------------------------
428 //mark this peak
429 //------------------------------------------------------------
430 UShort_t q=bins[k].GetQ();
431
432 bins[k].SetMask(bins[k].GetMask()|m);
433
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);
442}
443
444//______________________________________________________________________
445void AliITSClusterFinder::
446MakeCluster(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;
452
453 c.SetQ(c.GetQ()+q);
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);
458
459 bins[k].SetMask(0xFFFFFFFE);
460
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);
465}