- fixing warnings/coverity
[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"
45a7397a 32#include "AliLog.h"
b0f5e3fc 33
b0f5e3fc 34ClassImp(AliITSClusterFinder)
35
5d2c2f86 36extern AliRun *gAlice;
37
9de0700b 38//----------------------------------------------------------------------
aacedc3e 39AliITSClusterFinder::AliITSClusterFinder():
40TObject(),
aacedc3e 41fModule(0),
42fDigits(0),
43fNdigits(0),
8ba39da9 44fDetTypeRec(0),
aacedc3e 45fClusters(0),
aacedc3e 46fMap(0),
5d2c2f86 47fNPeaks(-1),
48fNModules(AliITSgeomTGeo::GetNModules()),
d59a647b 49fEvent(0),
50fZmin(0),
51fZmax(0),
52fXmin(0),
53fXmax(0){
9de0700b 54 // default cluster finder
aacedc3e 55 // Input:
56 // none.
57 // Output:
58 // none.
59 // Return:
60 // A default constructed AliITSCulsterFinder
40f5092a 61 for(Int_t i=0; i<2200; i++){
62 fNdet[i]=0;
63 fNlayer[i]=0;
64 }
9de0700b 65}
66//----------------------------------------------------------------------
8ba39da9 67AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
aacedc3e 68TObject(),
aacedc3e 69fModule(0),
70fDigits(0),
71fNdigits(0),
e56160b8 72fDetTypeRec(dettyp),
aacedc3e 73fClusters(0),
aacedc3e 74fMap(0),
5d2c2f86 75fNPeaks(-1),
76fNModules(AliITSgeomTGeo::GetNModules()),
d59a647b 77fEvent(0),
78fZmin(0),
79fZmax(0),
80fXmin(0),
81fXmax(0){
5d2c2f86 82 // default cluster finder
aacedc3e 83 // Standard constructor for cluster finder
84 // Input:
85 // AliITSsegmentation *seg The segmentation class to be used
86 // AliITSresponse *res The response class to be used
87 // Output:
88 // none.
89 // Return:
90 // A Standard constructed AliITSCulsterFinder
40f5092a 91 for(Int_t i=0; i<2200; i++){
92 fNdet[i]=0;
93 fNlayer[i]=0;
94 }
aacedc3e 95}
96//----------------------------------------------------------------------
8ba39da9 97AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
98 TClonesArray *digits):
aacedc3e 99TObject(),
aacedc3e 100fModule(0),
101fDigits(digits),
102fNdigits(0),
8ba39da9 103fDetTypeRec(dettyp),
aacedc3e 104fClusters(0),
aacedc3e 105fMap(0),
5d2c2f86 106fNPeaks(-1),
107fNModules(AliITSgeomTGeo::GetNModules()),
d59a647b 108fEvent(0),
109fZmin(0),
110fZmax(0),
111fXmin(0),
112fXmax(0){
5d2c2f86 113 // default cluster finder
aacedc3e 114 // Standard + cluster finder constructor
115 // Input:
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
119 // Output:
120 // none.
121 // Return:
122 // A Standard constructed AliITSCulsterFinder
b0f5e3fc 123
40f5092a 124 fNdigits = fDigits->GetEntriesFast();
125 for(Int_t i=0; i<2200; i++){
126 fNdet[i]=0;
127 fNlayer[i]=0;
128 }
b0f5e3fc 129}
a971ed8e 130
04366a57 131//______________________________________________________________________
d59a647b 132AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) :
133 TObject(source),
134 fModule(source.fModule),
135 fDigits(),
136 fNdigits(source.fNdigits),
137 fDetTypeRec(),
138 fClusters(),
139 fMap(),
140 fNPeaks(source.fNPeaks),
141 fNModules(source.fNModules),
142 fEvent(source.fEvent),
143 fZmin(source.fZmin),
144 fZmax(source.fZmax),
145 fXmin(source.fXmin),
146 fXmax(source.fXmax)
147{
04366a57 148 // Copy constructor
149 // Copies are not allowed. The method is protected to avoid misuse.
a971ed8e 150 AliError("Copy constructor not allowed\n");
04366a57 151}
a971ed8e 152
04366a57 153
154//______________________________________________________________________
e56160b8 155//AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
04366a57 156 // Assignment operator
157 // Assignment is not allowed. The method is protected to avoid misuse.
e56160b8 158// Fatal("= operator","Assignment operator not allowed\n");
159// return *this;
160//}
04366a57 161
9de0700b 162//----------------------------------------------------------------------
163AliITSClusterFinder::~AliITSClusterFinder(){
164 // destructor cluster finder
aacedc3e 165 // Input:
166 // none.
167 // Output:
168 // none.
169 // Return:
170 // none.
b0f5e3fc 171
aacedc3e 172 if(fMap) {delete fMap;}
9de0700b 173 // Zero local pointers. Other classes own these pointers.
9de0700b 174 fMap = 0;
175 fDigits = 0;
176 fNdigits = 0;
9de0700b 177 fNPeaks = 0;
7d62fb64 178 fDetTypeRec = 0;
04366a57 179
7d62fb64 180}
b0f5e3fc 181//__________________________________________________________________________
7d62fb64 182void AliITSClusterFinder::InitGeometry(){
1f3e997f 183 //
184 // Initialisation of ITS geometry
185 //
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);
1f3e997f 189 fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
04366a57 190 fNlayer[m] = lay-1;
191 }
b0f5e3fc 192}
04366a57 193
194
7d62fb64 195
04366a57 196
f8d9a5b8 197//______________________________________________________________________
aacedc3e 198Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
f8d9a5b8 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
206 // Inputs:
207 // TObjArray *digs Array to search for neighbors in
208 // Int_t i Index of digit for which we are searching for
209 // a neighbor of.
210 // Output:
211 // Int_t j[4] Index of one or more of the digits which is a
212 // neighbor of digit a index i.
213 // Return:
214 // Bool_t kTRUE if a neighbor was found kFALSE otherwise.
215 Int_t ix,jx,iz,jz,j;
216 const Bool_t kdiagonal=kFALSE;
217 Bool_t nei[4];
218
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;
223
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();
227 for(j=0;j<i;j++){
aacedc3e 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;}
f8d9a5b8 234 } // end for k
235 if(nei[0]||nei[1]) return kTRUE;
236 if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
237 // no Neighbors found.
238 return kFALSE;
239}
aacedc3e 240
241//______________________________________________________________________
7d62fb64 242void AliITSClusterFinder::Print(ostream *os) const{
aacedc3e 243 //Standard output format for this class
244 // Inputs:
245 // ostream *os Output stream
246 // Output:
247 // ostream *os Output stream
248 // Return:
249 // none.
250
aacedc3e 251 *os << fModule<<",";
252 *os << fNdigits<<",";
aacedc3e 253 *os << fNPeaks<<endl;
254}
04366a57 255//______________________________________________________________________
7d62fb64 256void AliITSClusterFinder::Read(istream *is) {
aacedc3e 257 //Standard input for this class
258 // Inputs:
259 // istream *is Input stream
260 // Output:
261 // istream *is Input stream
262 // Return:
263 // none.
264
aacedc3e 265 *is >> fModule;
266 *is >> fNdigits;
aacedc3e 267 *is >> fNPeaks;
268}
269//______________________________________________________________________
270ostream &operator<<(ostream &os,AliITSClusterFinder &source){
271 // Standard output streaming function.
272 // Inputs:
273 // ostream *os Output stream
274 // AliITSClusterFinder &source Class to be printed
275 // Output:
276 // ostream *os Output stream
277 // Return:
278 // none.
279
280 source.Print(&os);
281 return os;
282}
283//______________________________________________________________________
284istream &operator>>(istream &is,AliITSClusterFinder &source){
285 // Standard output streaming function.
286 // Inputs:
287 // istream *is Input stream
288 // AliITSClusterFinder &source Class to be read in.
289 // Output:
290 // istream *is Input stream
291 // Return:
292 // none.
293
294 source.Read(&is);
295 return is;
296}
45a7397a 297
5d2c2f86 298//______________________________________________________________________
45a7397a 299void AliITSClusterFinder::CheckLabels2(Int_t lab[10])
300{
5d2c2f86 301 //------------------------------------------------------------
302 // Tries to find mother's labels
303 //------------------------------------------------------------
33c3c91a 304 AliRunLoader *rl = AliRunLoader::Instance();
c605e5f7 305 if(!rl) return;
5d2c2f86 306 TTree *trK=(TTree*)rl->TreeK();
45a7397a 307 if (!trK) return;
308 //
309 int labS[10];
310 Int_t nlabels = 0;
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;
314 //
315 float mom[10];
316 for (Int_t i=0;i<nlabels;i++) {
317 Int_t label = labS[i];
318 mom[i] = 0;
319 if (label>=ntracks) continue;
320 TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
321 mom[i] = part->P();
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
325 //
326 if (part->GetStatusCode()>0) continue;
327 //
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; }
5d2c2f86 330 }
45a7397a 331 }
332 //
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]];
338 }
339 //
340 //compress labels -- if multi-times the same
341 for (Int_t i=0;i<10;i++) lab[i]=-2;
342 int nlabFin=0,j=0;
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];
5d2c2f86 346 }
45a7397a 347 //
5d2c2f86 348}
349
350//______________________________________________________________________
351void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
352 //add label to the cluster
33c3c91a 353 AliRunLoader *rl = AliRunLoader::Instance();
5d2c2f86 354 TTree *trK=(TTree*)rl->TreeK();
355 if(trK){
356 if(label<0) return; // In case of no label just exit
357
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;
363 if (lab[i]<0) {
364 lab[i]= label;
365 break;
366 }
367 }
368 }
369}
370
371
372//______________________________________________________________________
373void AliITSClusterFinder::
374FindCluster(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();
380 bins[k].Use();
381
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);
386 /*
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);
391 */
392}
393
394//______________________________________________________________________
395Bool_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;
409 return kTRUE;
410}
411
412//______________________________________________________________________
413void AliITSClusterFinder::
414FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
415 //------------------------------------------------------------
416 //find local maxima
417 //------------------------------------------------------------
418 if (n<31)
419 if (IsMaximum(k,max,b)) {
420 idx[n]=k; msk[n]=(2<<n);
421 n++;
422 }
423 b[k].SetMask(0);
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);
428}
429
430//______________________________________________________________________
431void AliITSClusterFinder::
432MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
433 //------------------------------------------------------------
434 //mark this peak
435 //------------------------------------------------------------
436 UShort_t q=bins[k].GetQ();
437
438 bins[k].SetMask(bins[k].GetMask()|m);
439
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);
448}
449
450//______________________________________________________________________
451void AliITSClusterFinder::
452MakeCluster(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;
d59a647b 458 if(c.GetQ()<0.01){ // first entry in cluster
459 fXmin=i;
460 fXmax=i;
461 fZmin=j;
462 fZmax=j;
463 }else{ // check cluster extension
464 if(i<fXmin) fXmin=i;
465 if(i>fXmax) fXmax=i;
466 if(j<fZmin) fZmin=j;
467 if(j>fZmax) fZmax=j;
468 }
5d2c2f86 469 c.SetQ(c.GetQ()+q);
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);
474
475 bins[k].SetMask(0xFFFFFFFE);
476
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);
481}