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