Update for Ds
[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 #include "AliLog.h"
33
34 ClassImp(AliITSClusterFinder)
35
36 extern AliRun *gAlice;
37
38 //----------------------------------------------------------------------
39 AliITSClusterFinder::AliITSClusterFinder():
40 TObject(),
41 fModule(0),
42 fDigits(0),
43 fNdigits(0),
44 fDetTypeRec(0),
45 fClusters(0),
46 fMap(0),
47 fNPeaks(-1),
48 fNModules(AliITSgeomTGeo::GetNModules()),
49 fEvent(0),
50 fZmin(0),
51 fZmax(0),
52 fXmin(0),
53 fXmax(0){
54     // default cluster finder
55     // Input:
56     //   none.
57     // Output:
58     //   none.
59     // Return:
60     //   A default constructed AliITSCulsterFinder
61   for(Int_t i=0; i<2200; i++){
62     fNdet[i]=0;
63     fNlayer[i]=0;
64   }
65 }
66 //----------------------------------------------------------------------
67 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
68 TObject(),
69 fModule(0),
70 fDigits(0),
71 fNdigits(0),
72 fDetTypeRec(dettyp),
73 fClusters(0),
74 fMap(0),
75 fNPeaks(-1),
76 fNModules(AliITSgeomTGeo::GetNModules()),
77 fEvent(0),
78 fZmin(0),
79 fZmax(0),
80 fXmin(0),
81 fXmax(0){
82     // default cluster finder
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
91   for(Int_t i=0; i<2200; i++){
92     fNdet[i]=0;
93     fNlayer[i]=0;
94   }
95 }
96 //----------------------------------------------------------------------
97 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
98                                          TClonesArray *digits):
99 TObject(),
100 fModule(0),
101 fDigits(digits),
102 fNdigits(0),
103 fDetTypeRec(dettyp),
104 fClusters(0),
105 fMap(0),
106 fNPeaks(-1),
107 fNModules(AliITSgeomTGeo::GetNModules()),
108 fEvent(0),
109 fZmin(0),
110 fZmax(0),
111 fXmin(0),
112 fXmax(0){
113     // default cluster finder
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
123
124   fNdigits = fDigits->GetEntriesFast();
125   for(Int_t i=0; i<2200; i++){
126     fNdet[i]=0;
127     fNlayer[i]=0;
128   }
129 }
130
131 //______________________________________________________________________
132 AliITSClusterFinder::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 {
148   // Copy constructor
149   // Copies are not allowed. The method is protected to avoid misuse.
150   AliError("Copy constructor not allowed\n");
151 }
152
153
154 //______________________________________________________________________
155 //AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
156   // Assignment operator
157   // Assignment is not allowed. The method is protected to avoid misuse.
158 //  Fatal("= operator","Assignment operator not allowed\n");
159 //  return *this;
160 //}
161
162 //----------------------------------------------------------------------
163 AliITSClusterFinder::~AliITSClusterFinder(){
164     // destructor cluster finder
165     // Input:
166     //   none.
167     // Output:
168     //   none.
169     // Return:
170     //   none.
171
172     if(fMap) {delete fMap;}
173     // Zero local pointers. Other classes own these pointers.
174     fMap          = 0;
175     fDigits       = 0;
176     fNdigits      = 0;
177     fNPeaks       = 0;
178     fDetTypeRec   = 0;
179
180 }
181 //__________________________________________________________________________
182 void AliITSClusterFinder::InitGeometry(){
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);
189     fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
190     fNlayer[m] = lay-1;
191   }
192 }
193
194
195
196
197 //______________________________________________________________________
198 Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
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++){
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;}
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 }
240
241 //______________________________________________________________________
242 void AliITSClusterFinder::Print(ostream *os) const{
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
251     *os << fModule<<",";
252     *os << fNdigits<<",";
253     *os << fNPeaks<<endl;
254 }
255 //______________________________________________________________________
256 void AliITSClusterFinder::Read(istream *is)  {
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
265     *is >> fModule;
266     *is >> fNdigits;
267     *is >> fNPeaks;
268 }
269 //______________________________________________________________________
270 ostream &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 //______________________________________________________________________
284 istream &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 }
297
298 //______________________________________________________________________
299 void AliITSClusterFinder::CheckLabels2(Int_t lab[10]) 
300 {
301   //------------------------------------------------------------
302   // Tries to find mother's labels
303   //------------------------------------------------------------
304   AliRunLoader *rl = AliRunLoader::Instance();
305   if(!rl) return;
306   TTree *trK=(TTree*)rl->TreeK();
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; }
330     }
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];
346   }
347   //
348 }
349
350 //______________________________________________________________________
351 void AliITSClusterFinder::AddLabel(Int_t lab[10], Int_t label) {
352   //add label to the cluster
353   AliRunLoader *rl = AliRunLoader::Instance();
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 //______________________________________________________________________
373 void AliITSClusterFinder:: 
374 FindCluster(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 //______________________________________________________________________
395 Bool_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 //______________________________________________________________________
413 void AliITSClusterFinder::
414 FindPeaks(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 //______________________________________________________________________
431 void AliITSClusterFinder::
432 MarkPeak(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 //______________________________________________________________________
451 void AliITSClusterFinder::
452 MakeCluster(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;
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   }
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 }