V2 clusterer moved to the standard framework. V2 clusters and recpoints are still...
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderV2.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2003, 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 //            Implementation of the ITS clusterer V2 class                //
17 //                                                                        //
18 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch            //
19 //                                                                        //
20 ///////////////////////////////////////////////////////////////////////////
21
22
23 #include "AliRun.h"
24 #include "AliITSClusterFinderV2.h"
25 #include "AliITSRecPoint.h"
26 #include "AliITSclusterV2.h"
27 #include "AliITSsimulationFastPoints.h"
28 #include "AliITS.h"
29 #include "AliITSgeom.h"
30 #include <TParticle.h>
31 #include "AliMC.h"
32
33 ClassImp(AliITSClusterFinderV2)
34
35 extern AliRun *gAlice;
36
37 AliITSClusterFinderV2::AliITSClusterFinderV2():AliITSClusterFinder(){
38
39   //Default constructor
40   fEvent = 0;
41   fModule = 0;
42
43   AliITSgeom *geom=(AliITSgeom*)fITS->GetITSgeom();
44
45   fNModules = geom->GetIndexMax();
46 }
47
48 //______________________________________________________________________
49 AliITSClusterFinderV2::AliITSClusterFinderV2(const AliITSClusterFinderV2 &source) : AliITSClusterFinder(source) {
50   // Copy constructor
51   // Copies are not allowed. The method is protected to avoid misuse.
52   Fatal("AliITSClusterFinderV2","Copy constructor not allowed\n");
53 }
54
55 //______________________________________________________________________
56 AliITSClusterFinderV2& AliITSClusterFinderV2::operator=(const AliITSClusterFinderV2& /* source */){
57   // Assignment operator
58   // Assignment is not allowed. The method is protected to avoid misuse.
59   Fatal("= operator","Assignment operator not allowed\n");
60   return *this;
61 }
62
63
64 //______________________________________________________________________
65 void AliITSClusterFinderV2::CheckLabels2(Int_t lab[10]) {
66   //------------------------------------------------------------
67   // Tries to find mother's labels
68   //------------------------------------------------------------
69   Int_t nlabels =0; 
70   for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++;
71   if(nlabels == 0) return; // In case of no labels just exit
72
73
74   Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
75
76   for (Int_t i=0;i<10;i++){
77     Int_t label = lab[i];
78     if (label>=0 && label<ntracks) {
79       TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
80
81       if (part->P() < 0.02) {
82           Int_t m=part->GetFirstMother();
83           if (m<0) {    
84             continue;
85           }
86           if (part->GetStatusCode()>0) {
87             continue;
88           }
89           lab[i]=m;       
90       }
91       else
92         if (part->P() < 0.12 && nlabels>3) {
93           lab[i]=-2;
94           nlabels--;
95         } 
96     }
97     else{
98       if ( (label>ntracks||label <0) && nlabels>3) {
99         lab[i]=-2;
100         nlabels--;
101       } 
102     }
103   }  
104   if (nlabels>3){
105     for (Int_t i=0;i<10;i++){
106       if (nlabels>3){
107         Int_t label = lab[i];
108         if (label>=0 && label<ntracks) {
109           TParticle *part=(TParticle*)gAlice->GetMCApp()->Particle(label);
110           if (part->P() < 0.1) {
111             lab[i]=-2;
112             nlabels--;
113           }
114         }
115       }
116     }
117   }
118
119   //compress labels -- if multi-times the same
120   Int_t lab2[10];
121   for (Int_t i=0;i<10;i++) lab2[i]=-2;
122   for (Int_t i=0;i<10  ;i++){
123     if (lab[i]<0) continue;
124     for (Int_t j=0;j<10 &&lab2[j]!=lab[i];j++){
125       if (lab2[j]<0) {
126         lab2[j]= lab[i];
127         break;
128       }
129     }
130   }
131   for (Int_t j=0;j<10;j++) lab[j]=lab2[j];
132   
133 }
134
135 //______________________________________________________________________
136 void AliITSClusterFinderV2::AddLabel(Int_t lab[10], Int_t label) {
137
138   //add label to the cluster
139
140   if(label<0) return; // In case of no label just exit
141
142   Int_t ntracks = gAlice->GetMCApp()->GetNtrack();
143   if (label>ntracks) return;
144   for (Int_t i=0;i<10;i++){
145     //    if (label<0) break;
146     if (lab[i]==label) break;
147     if (lab[i]<0) {
148       lab[i]= label;
149       break;
150     }
151   }
152 }
153
154
155 //______________________________________________________________________
156 void AliITSClusterFinderV2:: 
157 FindCluster(Int_t k,Int_t maxz,AliBin *bins,Int_t &n,Int_t *idx) {
158   //------------------------------------------------------------
159   // returns an array of indices of digits belonging to the cluster
160   // (needed when the segmentation is not regular) 
161   //------------------------------------------------------------
162   if (n<200) idx[n++]=bins[k].GetIndex();
163   bins[k].Use();
164
165   if (bins[k-maxz].IsNotUsed()) FindCluster(k-maxz,maxz,bins,n,idx);
166   if (bins[k-1   ].IsNotUsed()) FindCluster(k-1   ,maxz,bins,n,idx);
167   if (bins[k+maxz].IsNotUsed()) FindCluster(k+maxz,maxz,bins,n,idx);
168   if (bins[k+1   ].IsNotUsed()) FindCluster(k+1   ,maxz,bins,n,idx);
169   /*
170   if (bins[k-maxz-1].IsNotUsed()) FindCluster(k-maxz-1,maxz,bins,n,idx);
171   if (bins[k-maxz+1].IsNotUsed()) FindCluster(k-maxz+1,maxz,bins,n,idx);
172   if (bins[k+maxz-1].IsNotUsed()) FindCluster(k+maxz-1,maxz,bins,n,idx);
173   if (bins[k+maxz+1].IsNotUsed()) FindCluster(k+maxz+1,maxz,bins,n,idx);
174   */
175 }
176
177 //______________________________________________________________________
178 Bool_t AliITSClusterFinderV2::IsMaximum(Int_t k,Int_t max,const AliBin *bins) {
179   //------------------------------------------------------------
180   //is this a local maximum ?
181   //------------------------------------------------------------
182   UShort_t q=bins[k].GetQ();
183   if (q==1023) return kFALSE;
184   if (bins[k-max].GetQ() > q) return kFALSE;
185   if (bins[k-1  ].GetQ() > q) return kFALSE; 
186   if (bins[k+max].GetQ() > q) return kFALSE; 
187   if (bins[k+1  ].GetQ() > q) return kFALSE; 
188   if (bins[k-max-1].GetQ() > q) return kFALSE;
189   if (bins[k+max-1].GetQ() > q) return kFALSE; 
190   if (bins[k+max+1].GetQ() > q) return kFALSE; 
191   if (bins[k-max+1].GetQ() > q) return kFALSE;
192   return kTRUE; 
193 }
194
195 //______________________________________________________________________
196 void AliITSClusterFinderV2::
197 FindPeaks(Int_t k,Int_t max,AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
198   //------------------------------------------------------------
199   //find local maxima
200   //------------------------------------------------------------
201   if (n<31)
202   if (IsMaximum(k,max,b)) {
203     idx[n]=k; msk[n]=(2<<n);
204     n++;
205   }
206   b[k].SetMask(0);
207   if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
208   if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
209   if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
210   if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
211 }
212
213 //______________________________________________________________________
214 void AliITSClusterFinderV2::
215 MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
216   //------------------------------------------------------------
217   //mark this peak
218   //------------------------------------------------------------
219   UShort_t q=bins[k].GetQ();
220
221   bins[k].SetMask(bins[k].GetMask()|m); 
222
223   if (bins[k-max].GetQ() <= q)
224      if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
225   if (bins[k-1  ].GetQ() <= q)
226      if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
227   if (bins[k+max].GetQ() <= q)
228      if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
229   if (bins[k+1  ].GetQ() <= q)
230      if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
231 }
232
233 //______________________________________________________________________
234 void AliITSClusterFinderV2::
235 MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,AliITSclusterV2 &c) {
236   //------------------------------------------------------------
237   //make cluster using digits of this peak
238   //------------------------------------------------------------
239   Float_t q=(Float_t)bins[k].GetQ();
240   Int_t i=k/max, j=k-i*max;
241
242   c.SetQ(c.GetQ()+q);
243   c.SetY(c.GetY()+i*q); 
244   c.SetZ(c.GetZ()+j*q); 
245   c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
246   c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);
247
248   bins[k].SetMask(0xFFFFFFFE);
249   
250   if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
251   if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
252   if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
253   if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
254 }
255
256 //______________________________________________________________________
257 Int_t AliITSClusterFinderV2::Hits2Clusters(TTree *hTree, TTree *cTree) {
258   //------------------------------------------------------------
259   // This function creates ITS clusters
260   //------------------------------------------------------------
261
262   AliITSgeom *geom=fITS->GetITSgeom();
263   Int_t mmax=geom->GetIndexMax();
264
265   fITS->InitModules(-1,mmax);
266   fITS->FillModules(hTree,0);
267
268   TClonesArray *clusters=new TClonesArray("AliITSclusterV2",1000);
269   TBranch *branch=cTree->GetBranch("Clusters");
270   if (!branch) cTree->Branch("Clusters",&clusters);
271   else branch->SetAddress(&clusters);
272
273   static TClonesArray *points=fITS->RecPoints();
274   AliITSsimulationFastPoints sim;
275   Int_t ncl=0;
276   for (Int_t m=0; m<mmax; m++) {
277     AliITSmodule *mod=fITS->GetModule(m);      
278     sim.CreateFastRecPoints(mod,m,gRandom);      
279
280     RecPoints2Clusters(points, m, clusters);
281     fITS->ResetRecPoints();
282
283     ncl+=clusters->GetEntriesFast();
284     cTree->Fill();
285     clusters->Clear();
286   }
287
288   Info("Hits2Clusters","Number of found fast clusters : %d",ncl);
289
290   //cTree->Write();
291
292   delete clusters;
293
294   return 0;
295 }
296