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