Added option for different binning of DCAxy axis in THnSparse. Same width for all...
[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 ClassImp(AliITSClusterFinder)
36
37 extern AliRun *gAlice;
38
39 //----------------------------------------------------------------------
40 AliITSClusterFinder::AliITSClusterFinder():
41 TObject(),
42 fModule(0),
43 fDigits(0),
44 fNdigits(0),
45 fDetTypeRec(0),
46 fClusters(0),
47 fMap(0),
48 fNPeaks(-1),
49 fNModules(AliITSgeomTGeo::GetNModules()),
50 fEvent(0),
51 fZmin(0),
52 fZmax(0),
53 fXmin(0),
54 fXmax(0),
55 fNClusters(0),
56 fRawID2ClusID(0)
57 {
58     // default cluster finder
59     // Input:
60     //   none.
61     // Output:
62     //   none.
63     // Return:
64     //   A default constructed AliITSCulsterFinder
65   for(Int_t i=0; i<2200; i++){
66     fNdet[i]=0;
67     fNlayer[i]=0;
68   }
69 }
70 //----------------------------------------------------------------------
71 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
72 TObject(),
73 fModule(0),
74 fDigits(0),
75 fNdigits(0),
76 fDetTypeRec(dettyp),
77 fClusters(0),
78 fMap(0),
79 fNPeaks(-1),
80 fNModules(AliITSgeomTGeo::GetNModules()),
81 fEvent(0),
82 fZmin(0),
83 fZmax(0),
84 fXmin(0),
85 fXmax(0),
86 fNClusters(0),
87 fRawID2ClusID(0)
88 {
89     // default cluster finder
90     // Standard constructor for cluster finder
91     // Input:
92     //   AliITSsegmentation *seg  The segmentation class to be used
93     //   AliITSresponse     *res  The response class to be used
94     // Output:
95     //   none.
96     // Return:
97     //   A Standard constructed AliITSCulsterFinder
98   for(Int_t i=0; i<2200; i++){
99     fNdet[i]=0;
100     fNlayer[i]=0;
101   }
102 }
103 //----------------------------------------------------------------------
104 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
105                                          TClonesArray *digits):
106 TObject(),
107 fModule(0),
108 fDigits(digits),
109 fNdigits(0),
110 fDetTypeRec(dettyp),
111 fClusters(0),
112 fMap(0),
113 fNPeaks(-1),
114 fNModules(AliITSgeomTGeo::GetNModules()),
115 fEvent(0),
116 fZmin(0),
117 fZmax(0),
118 fXmin(0),
119 fXmax(0),
120 fNClusters(0),
121 fRawID2ClusID(0)
122 {
123     // default cluster finder
124     // Standard + cluster finder constructor
125     // Input:
126     //   AliITSsegmentation *seg  The segmentation class to be used
127     //   AliITSresponse     *res  The response class to be used
128     //   TClonesArray    *digits  Array of digits to be used
129     // Output:
130     //   none.
131     // Return:
132     //   A Standard constructed AliITSCulsterFinder
133
134   fNdigits = fDigits->GetEntriesFast();
135   for(Int_t i=0; i<2200; i++){
136     fNdet[i]=0;
137     fNlayer[i]=0;
138   }
139 }
140
141 //______________________________________________________________________
142 AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : 
143   TObject(source),
144   fModule(source.fModule),
145   fDigits(),
146   fNdigits(source.fNdigits),
147   fDetTypeRec(),
148   fClusters(),
149   fMap(),
150   fNPeaks(source.fNPeaks),
151   fNModules(source.fNModules),
152   fEvent(source.fEvent),
153   fZmin(source.fZmin),
154   fZmax(source.fZmax),
155   fXmin(source.fXmin),
156   fXmax(source.fXmax),
157   fNClusters(source.fNClusters),
158   fRawID2ClusID(source.fRawID2ClusID) 
159 {
160   // Copy constructor
161   // Copies are not allowed. The method is protected to avoid misuse.
162   AliError("Copy constructor not allowed\n");
163 }
164
165
166 //______________________________________________________________________
167 //AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
168   // Assignment operator
169   // Assignment is not allowed. The method is protected to avoid misuse.
170 //  Fatal("= operator","Assignment operator not allowed\n");
171 //  return *this;
172 //}
173
174 //----------------------------------------------------------------------
175 AliITSClusterFinder::~AliITSClusterFinder(){
176     // destructor cluster finder
177     // Input:
178     //   none.
179     // Output:
180     //   none.
181     // Return:
182     //   none.
183
184     if(fMap) {delete fMap;}
185     // Zero local pointers. Other classes own these pointers.
186     fMap          = 0;
187     fDigits       = 0;
188     fNdigits      = 0;
189     fNPeaks       = 0;
190     fDetTypeRec   = 0;
191
192 }
193 //__________________________________________________________________________
194 void AliITSClusterFinder::InitGeometry(){
195  //
196  // Initialisation of ITS geometry
197  //
198   Int_t mmax=AliITSgeomTGeo::GetNModules();
199   for (Int_t m=0; m<mmax; m++) {
200     Int_t lay,lad,det; AliITSgeomTGeo::GetModuleId(m,lay,lad,det);
201     fNdet[m] = (lad-1)*AliITSgeomTGeo::GetNDetectors(lay) + (det-1);
202     fNlayer[m] = lay-1;
203   }
204 }
205
206
207
208
209 //______________________________________________________________________
210 Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
211     // Locagical function which checks to see if digit i has a neighbor.
212     // If so, then it returns kTRUE and its neighbor index j.
213     // This routine checks if the digits are side by side or one before the
214     // other. Requires that the array of digits be in proper order.
215     // Returns kTRUE in the following cases.
216     //                 ji   0j   if kdiagonal  j0    0i
217     //                 00   0i   if kdiagonal  0i    j0
218     // Inputs:
219     //    TObjArray *digs   Array to search for neighbors in
220     //    Int_t      i      Index of digit for which we are searching for
221     //                      a neighbor of.
222     // Output:
223     //    Int_t      j[4]   Index of one or more of the digits which is a
224     //                      neighbor of digit a index i.
225     // Return:
226     //    Bool_t            kTRUE if a neighbor was found kFALSE otherwise.
227     Int_t ix,jx,iz,jz,j;
228     const Bool_t kdiagonal=kFALSE;
229     Bool_t nei[4];
230
231     // No neighbors found if array empty.
232     if(digs->GetEntriesFast()<=0) return kFALSE;
233     // can not be a digit with first element or elements out or range
234     if(i<=0 || i>=digs->GetEntriesFast()) return kFALSE;
235
236     for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
237     ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
238     iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
239     for(j=0;j<i;j++){
240         jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
241         jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
242         if(jx+1==ix && jz  ==iz){n[0] = j;nei[0] = kTRUE;}
243         if(jx  ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
244         if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
245         if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
246     } // end for k
247     if(nei[0]||nei[1]) return kTRUE;
248     if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
249     // no Neighbors found.
250     return kFALSE;
251 }
252
253 //______________________________________________________________________
254 void AliITSClusterFinder::Print(ostream *os) const{
255     //Standard output format for this class
256     // Inputs:
257     //    ostream *os   Output stream
258     // Output:
259     //    ostream *os   Output stream
260     // Return:
261     //    none.
262
263     *os << fModule<<",";
264     *os << fNdigits<<",";
265     *os << fNPeaks<<endl;
266 }
267 //______________________________________________________________________
268 void AliITSClusterFinder::Read(istream *is)  {
269     //Standard input for this class
270     // Inputs:
271     //    istream *is   Input stream
272     // Output:
273     //    istream *is   Input stream
274     // Return:
275     //    none.
276
277     *is >> fModule;
278     *is >> fNdigits;
279     *is >> fNPeaks;
280 }
281 //______________________________________________________________________
282 ostream &operator<<(ostream &os,AliITSClusterFinder &source){
283     // Standard output streaming function.
284     // Inputs:
285     //    ostream             *os     Output stream
286     //    AliITSClusterFinder &source Class to be printed
287     // Output:
288     //    ostream             *os     Output stream
289     // Return:
290     //    none.
291
292     source.Print(&os);
293     return os;
294 }
295 //______________________________________________________________________
296 istream &operator>>(istream &is,AliITSClusterFinder &source){
297     // Standard output streaming function.
298     // Inputs:
299     //    istream              *is      Input stream
300     //     AliITSClusterFinder &source  Class to be read in.
301     // Output:
302     //    istream              *is      Input stream
303     // Return:
304     //    none.
305
306     source.Read(&is);
307     return is;
308 }
309
310 //______________________________________________________________________
311 void AliITSClusterFinder::CheckLabels2(Int_t lab[10]) 
312 {
313   //------------------------------------------------------------
314   // Tries to find mother's labels
315   //------------------------------------------------------------
316   AliRunLoader *rl = AliRunLoader::Instance();
317   if(!rl) return;
318   TTree *trK=(TTree*)rl->TreeK();
319   if (!trK) return;
320   //
321   int labS[10];
322   Int_t nlabels = 0; 
323   Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
324   for (Int_t i=0;i<10;i++) if (lab[i]>=0) labS[nlabels++] = lab[i];
325   if (nlabels==0) return;
326   //
327   float mom[10];
328   for (Int_t i=0;i<nlabels;i++) {
329     Int_t label = labS[i];
330     mom[i] = 0;
331     if (label>=ntracks) continue;
332     TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
333     mom[i] = part->P();
334     if (part->P() < 0.02) {    // reduce soft particles from the same cluster
335       Int_t m=part->GetFirstMother();
336       if (m<0) continue; // primary
337       //
338       if (part->GetStatusCode()>0) continue;
339       //
340       // if the parent is within the same cluster, reassign the label to it
341       for (int j=0;j<nlabels;j++) if (labS[j]==m) { labS[i] = m; break; }
342     }
343   } 
344   //
345   if (nlabels>3) { // only 3 labels are stored in cluster, sort in decreasing momentum
346     int ind[10],labSS[10];
347     TMath::Sort(nlabels,mom,ind);
348     for (int i=nlabels;i--;) labSS[i] = labS[i];
349     for (int i=0;i<nlabels;i++) labS[i] = labSS[ind[i]]; 
350   }
351   //
352   //compress labels -- if multi-times the same
353   for (Int_t i=0;i<10;i++) lab[i]=-2;
354   int nlabFin=0,j=0;
355   for (int i=0;i<nlabels;i++) {
356     for (j=0;j<nlabFin;j++) if (labS[i]==lab[j]) break; // the label already there
357     if (j==nlabFin) lab[nlabFin++] = labS[i];
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   if (fRawID2ClusID) { // RS: Register cluster id in raw words list
489     int rwid = bins[k].GetRawID();
490     if (fRawID2ClusID->GetSize()<=rwid) fRawID2ClusID->Set( (rwid+10)<<1 );
491     (*fRawID2ClusID)[rwid] = fNClusters+1; // RS: store clID+1 as a reference to the cluster
492   }
493   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
494   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
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 }