Adding some further mother volumes to speed-up the overlap checking and particle...
[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   TTree *trK=(TTree*)rl->TreeK();
274
275   if(trK){
276     Int_t nlabels =0; 
277     for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
278     if(nlabels == 0) return; // In case of no labels just exit
279
280
281     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
282
283     for (Int_t i=0;i<10;i++){
284       Int_t label = lab[i];
285       if (label>=0 && label<ntracks) {
286         TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
287
288         if (part->P() < 0.02) {
289           Int_t m=part->GetFirstMother();
290           if (m<0) {    
291             continue;
292           }
293           if (part->GetStatusCode()>0) {
294             continue;
295           }
296           lab[i]=m;       
297         }
298         else
299           if (part->P() < 0.12 && nlabels>3) {
300             lab[i]=-2;
301             nlabels--;
302           } 
303       }
304       else{
305         if ( (label>ntracks||label <0) && nlabels>3) {
306           lab[i]=-2;
307           nlabels--;
308         } 
309       }
310     }  
311     if (nlabels>3){
312       for (Int_t i=0;i<10;i++){
313         if (nlabels>3){
314           Int_t label = lab[i];
315           if (label>=0 && label<ntracks) {
316             TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
317             if (part->P() < 0.1) {
318               lab[i]=-2;
319               nlabels--;
320             }
321           }
322         }
323       }
324     }
325
326     //compress labels -- if multi-times the same
327     Int_t lab2[10];
328     for (Int_t i=0;i<10;i++) lab2[i]=-2;
329     for (Int_t i=0;i<10  ;i++){
330       if (lab[i]<0) continue;
331       for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
332         if (lab2[j]<0) {
333           lab2[j]= lab[i];
334           break;
335         }
336       }
337     }
338     for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
339   
340   }
341 }
342
343 //______________________________________________________________________
344 void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
345   //add label to the cluster
346   AliRunLoader *rl = AliRunLoader::Instance();
347   TTree *trK=(TTree*)rl->TreeK();
348   if(trK){
349     if(label<0) return; // In case of no label just exit
350
351     Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
352     if (label>ntracks) return;
353     for (Int_t i=0;i<10;i++){
354       //    if (label<0) break;
355       if (lab[i]==label) break;
356       if (lab[i]<0) {
357         lab[i]= label;
358         break;
359       }
360     }
361   }
362 }
363
364
365 //______________________________________________________________________
366 void AliITSClusterFinder:: 
367 FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
368   //------------------------------------------------------------
369   // returns an array of indices of digits belonging to the cluster
370   // (needed when the segmentation is not regular) 
371   //------------------------------------------------------------
372   if (n<200) idx[n++]=bins[k].GetIndex();
373   bins[k].Use();
374
375   if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
376   if (bins[k-1   ].IsNotUsed()) FindCluster(k-1   ,maxz,bins,n,idx);
377   if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
378   if (bins[k+1   ].IsNotUsed()) FindCluster(k+1   ,maxz,bins,n,idx);
379   /*
380   if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
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   */
385 }
386
387 //______________________________________________________________________
388 Bool_t AliITSClusterFinder::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
389   //------------------------------------------------------------
390   //is this a local maximum ?
391   //------------------------------------------------------------
392   UShort_t q=bins[k].GetQ();
393   if (q==1023) return kFALSE;
394   if (bins[k-max].GetQ() > q) return kFALSE;
395   if (bins[k-1  ].GetQ() > q) return kFALSE; 
396   if (bins[k+max].GetQ() > q) return kFALSE; 
397   if (bins[k+1  ].GetQ() > q) return kFALSE; 
398   if (bins[k-max-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   return kTRUE; 
403 }
404
405 //______________________________________________________________________
406 void AliITSClusterFinder::
407 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
408   //------------------------------------------------------------
409   //find local maxima
410   //------------------------------------------------------------
411   if (n<31)
412   if (IsMaximum(k,max,b)) {
413     idx[n]=k; msk[n]=(2<<n);
414     n++;
415   }
416   b[k].SetMask(0);
417   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
418   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
419   if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
420   if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
421 }
422
423 //______________________________________________________________________
424 void AliITSClusterFinder::
425 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
426   //------------------------------------------------------------
427   //mark this peak
428   //------------------------------------------------------------
429   UShort_t q=bins[k].GetQ();
430
431   bins[k].SetMask(bins[k].GetMask()|m); 
432
433   if (bins[k-max].GetQ() <= q)
434      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
435   if (bins[k-1  ].GetQ() <= q)
436      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
437   if (bins[k+max].GetQ() <= q)
438      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
439   if (bins[k+1  ].GetQ() <= q)
440      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
441 }
442
443 //______________________________________________________________________
444 void AliITSClusterFinder::
445 MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSRecPoint &c) {
446   //------------------------------------------------------------
447   //make cluster using digits of this peak
448   //------------------------------------------------------------
449   Float_t q=(Float_t)bins[k].GetQ();
450   Int_t i=k/max, j=k-i*max;
451
452   c.SetQ(c.GetQ()+q);
453   c.SetY(c.GetY()+i*q);
454   c.SetZ(c.GetZ()+j*q);
455   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
456   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
457
458   bins[k].SetMask(0xFFFFFFFE);
459   
460   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
461   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
462   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
463   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
464 }