Added two missing includes to allow macro compilation (thanks to Laurent for remarkin...
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinder.cxx
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 "AliMC.h"
32
33 ClassImp(AliITSClusterFinder)
34
35 extern AliRun *gAlice;
36
37 //----------------------------------------------------------------------
38 AliITSClusterFinder::AliITSClusterFinder():
39 TObject(),
40 fModule(0),
41 fDigits(0),
42 fNdigits(0),
43 fDetTypeRec(0),
44 fClusters(0),
45 fMap(0),
46 fNPeaks(-1),
47 fNModules(AliITSgeomTGeo::GetNModules()),
48 fEvent(0){
49     // default cluster finder
50     // Input:
51     //   none.
52     // Output:
53     //   none.
54     // Return:
55     //   A default constructed AliITSCulsterFinder
56 }
57 //----------------------------------------------------------------------
58 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
59 TObject(),
60 fModule(0),
61 fDigits(0),
62 fNdigits(0),
63 fDetTypeRec(dettyp),
64 fClusters(0),
65 fMap(0),
66 fNPeaks(-1),
67 fNModules(AliITSgeomTGeo::GetNModules()),
68 fEvent(0){
69     // default cluster finder
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
78
79 }
80 //----------------------------------------------------------------------
81 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
82                                          TClonesArray *digits):
83 TObject(),
84 fModule(0),
85 fDigits(digits),
86 fNdigits(0),
87 fDetTypeRec(dettyp),
88 fClusters(0),
89 fMap(0),
90 fNPeaks(-1),
91 fNModules(AliITSgeomTGeo::GetNModules()),
92 fEvent(0){
93     // default cluster finder
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
103
104     fNdigits = fDigits->GetEntriesFast();
105 }
106
107 //______________________________________________________________________
108 AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : TObject(source),
109 fModule(source.fModule),
110 fDigits(),
111 fNdigits(source.fNdigits),
112 fDetTypeRec(),
113 fClusters(),
114 fMap(),
115 fNPeaks(source.fNPeaks),
116 fNModules(source.fNModules),
117 fEvent(source.fEvent) {
118   // Copy constructor
119   // Copies are not allowed. The method is protected to avoid misuse.
120   AliError("Copy constructor not allowed\n");
121 }
122
123
124 //______________________________________________________________________
125 //AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
126   // Assignment operator
127   // Assignment is not allowed. The method is protected to avoid misuse.
128 //  Fatal("= operator","Assignment operator not allowed\n");
129 //  return *this;
130 //}
131
132 //----------------------------------------------------------------------
133 AliITSClusterFinder::~AliITSClusterFinder(){
134     // destructor cluster finder
135     // Input:
136     //   none.
137     // Output:
138     //   none.
139     // Return:
140     //   none.
141
142     if(fMap) {delete fMap;}
143     // Zero local pointers. Other classes own these pointers.
144     fMap          = 0;
145     fDigits       = 0;
146     fNdigits      = 0;
147     fNPeaks       = 0;
148     fDetTypeRec   = 0;
149
150 }
151 //__________________________________________________________________________
152 void AliITSClusterFinder::InitGeometry(){
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);
159     fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
160     fNlayer[m] = lay-1;
161   }
162 }
163
164
165
166
167 //______________________________________________________________________
168 Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
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++){
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;}
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 }
210
211 //______________________________________________________________________
212 void AliITSClusterFinder::Print(ostream *os) const{
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
221     *os << fModule<<",";
222     *os << fNdigits<<",";
223     *os << fNPeaks<<endl;
224 }
225 //______________________________________________________________________
226 void AliITSClusterFinder::Read(istream *is)  {
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
235     *is >> fModule;
236     *is >> fNdigits;
237     *is >> fNPeaks;
238 }
239 //______________________________________________________________________
240 ostream &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 //______________________________________________________________________
254 istream &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 }
267 //______________________________________________________________________
268 void AliITSClusterFinder::CheckLabels2(Int_t lab[10]) {
269   //------------------------------------------------------------
270   // Tries to find mother's labels
271   //------------------------------------------------------------
272   AliRunLoader *rl = AliRunLoader::Instance();
273   if(!rl) return;
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 //______________________________________________________________________
345 void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
346   //add label to the cluster
347   AliRunLoader *rl = AliRunLoader::Instance();
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 //______________________________________________________________________
367 void AliITSClusterFinder:: 
368 FindCluster(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 //______________________________________________________________________
389 Bool_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 //______________________________________________________________________
407 void AliITSClusterFinder::
408 FindPeaks(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 //______________________________________________________________________
425 void AliITSClusterFinder::
426 MarkPeak(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 //______________________________________________________________________
445 void AliITSClusterFinder::
446 MakeCluster(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 }