V2 clusterer moved to the standard framework. V2 clusters and recpoints are still...
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSPD.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 //  Cluster finder                                                       ///
17 //  for Silicon pixels                                                    // 
18 //                                                                        //
19 ////////////////////////////////////////////////////////////////////////////        
20 #include "AliITSClusterFinderSPD.h"
21 #include "AliITS.h"
22 #include "AliITSdigitSPD.h"
23 #include "AliITSRawClusterSPD.h"
24 #include "AliITSRecPoint.h"
25 #include "AliITSresponseSPD.h"
26 #include "AliITSsegmentationSPD.h"
27 #include "AliRun.h"
28
29 //#define DEBUG
30
31 ClassImp(AliITSClusterFinderSPD)
32
33 //______________________________________________________________________
34 AliITSClusterFinderSPD::AliITSClusterFinderSPD():AliITSClusterFinder(),
35 fDz(0.0),
36 fDx(0.0),
37 fMinNCells(0){
38     // constructor
39 }
40 //----------------------------------------------------------
41 AliITSClusterFinderSPD::AliITSClusterFinderSPD(AliITSsegmentation *seg,
42                                                AliITSresponse *res):
43 AliITSClusterFinder(seg,res),
44 fDz(0.0),
45 fDx(0.0),
46 fMinNCells(0){
47     // constructor
48
49     SetDx();
50     SetDz();
51 }
52 //----------------------------------------------------------
53 AliITSClusterFinderSPD::AliITSClusterFinderSPD(AliITSsegmentation *seg,
54                                                TClonesArray *digits,
55                                                TClonesArray *recp):
56 AliITSClusterFinder(seg,0),
57 fDz(0.0),
58 fDx(0.0),
59 fMinNCells(0){
60     // constructor
61
62     SetDigits(digits);
63     SetClusters(recp);
64     SetDx();
65     SetDz();
66 }
67
68 //______________________________________________________________________
69 AliITSClusterFinderSPD::AliITSClusterFinderSPD(const AliITSClusterFinderSPD &source) : AliITSClusterFinder(source) {
70   // Copy constructor
71   // Copies are not allowed. The method is protected to avoid misuse.
72   Fatal("AliITSClusterFinderSPD","Copy constructor not allowed\n");
73 }
74
75 //______________________________________________________________________
76 AliITSClusterFinderSPD& AliITSClusterFinderSPD::operator=(const AliITSClusterFinderSPD& /* source */){
77   // Assignment operator
78   // Assignment is not allowed. The method is protected to avoid misuse.
79   Fatal("= operator","Assignment operator not allowed\n");
80   return *this;
81 }
82 //______________________________________________________________________
83 void AliITSClusterFinderSPD::FindRawClusters(Int_t module){   
84     // input of Cluster Finder
85     Int_t   digitcount  = 0;
86     Int_t   numberd     = 100000;
87     Int_t   *digx       = new Int_t[numberd];
88     Int_t   *digz       = new Int_t[numberd];
89     Int_t   *digtr1     = new Int_t[numberd];
90     Int_t   *digtr2     = new Int_t[numberd];
91     Int_t   *digtr3     = new Int_t[numberd];
92     Int_t   *digtr4     = new Int_t[numberd];
93     //  output of Cluster Finder    
94     Int_t   numberc     = 10000;
95     Double_t *xcenterl   = new Double_t[numberc];
96     Double_t *zcenterl   = new Double_t[numberc];
97     Double_t *errxcenter = new Double_t[numberc];
98     Double_t *errzcenter = new Double_t[numberc];
99     Int_t   *tr1clus    = new Int_t[numberc];
100     Int_t   *tr2clus    = new Int_t[numberc];
101     Int_t   *tr3clus    = new Int_t[numberc];
102     Int_t   nclus;
103
104     SetModule(module);
105     digitcount=0;
106     Int_t ndigits = Digits()->GetEntriesFast();  
107     if (!ndigits) return;
108
109     AliITSdigitSPD *dig;
110     Int_t ndig=0,i;
111     /*
112     if(GetDebug(4)){
113         cout << "SPD - FindRawclusters"<<endl;
114         scanf("%d",&ndig);
115     } // end if GetDebug
116     */
117     for(ndig=0; ndig<ndigits; ndig++) {
118         dig= (AliITSdigitSPD*)GetDigit(ndig);
119         digx[digitcount] = dig->GetCoord2()+1;  //starts at 1
120         digz[digitcount] = dig->GetCoord1()+1;  //starts at 1
121         digtr1[digitcount] = dig->GetTrack(0);
122         digtr2[digitcount] = -3;
123         digtr3[digitcount] = -3;
124         if(GetDebug(5)){
125             cout << "digtr1["<<digitcount <<"]="<<digtr1[digitcount];
126             cout << " fTracks["<<0<<"]="<<dig->GetTrack(0)<<": ";
127         } // end if GetDebug
128         i=1;
129         while(digtr1[digitcount]==dig->GetTrack(i) && i<dig->GetNTracks()) i++;
130         if(GetDebug(5)) cout << " fTracks["<<i<<"]="<<dig->GetTrack(i);
131         if(i<dig->GetNTracks()){
132             digtr2[digitcount] = dig->GetTrack(i);
133             if(GetDebug(5)) cout<<"digtr2["<<digitcount <<"]="
134                                 <<digtr2[digitcount]<<": ";
135             while((digtr1[digitcount]==dig->GetTrack(i) || 
136                    digtr2[digitcount]==dig->GetTrack(i))&&
137                   i<=dig->GetNTracks()) i++;
138             if(i<dig->GetNTracks()) digtr3[digitcount] = dig->GetTrack(i);
139             if(GetDebug(5)){
140                 cout << " fTracks["<<i<<"]=";
141                 if(i<dig->GetNTracks()) cout <<dig->GetTrack(i);
142                 cout << "digtr3["<<digitcount <<"]="<<digtr3[digitcount];
143             } // end if GetDebug
144         } // end if
145         if(GetDebug(4)) cout<<endl;
146         digtr4[digitcount] = dig->GetSignal();
147         digitcount++;
148     } // end for ndig
149     ClusterFinder(digitcount,digx,digz,digtr1,digtr2,digtr3,digtr4,
150                   nclus,xcenterl,zcenterl,errxcenter,errzcenter,
151                   tr1clus, tr2clus, tr3clus);
152     DigitToPoint(nclus,xcenterl,zcenterl,errxcenter,errzcenter,
153                  tr1clus, tr2clus, tr3clus);
154     delete[] digx;
155     delete[] digz;
156     delete[] digtr1;
157     delete[] digtr2;
158     delete[] digtr3;
159     delete[] digtr4;
160     delete[] xcenterl;
161     delete[] zcenterl;
162     delete[] errxcenter;
163     delete[] errzcenter;
164     delete[] tr1clus;
165     delete[] tr2clus;
166     delete[] tr3clus;
167 }
168 //----------------------------------------------------------------------
169 void AliITSClusterFinderSPD::ClusterFinder(Int_t ndigits,Int_t digx[],
170                                            Int_t digz[],Int_t digtr1[],
171                                            Int_t digtr2[],Int_t digtr3[],
172                                            Int_t digtr4[],Int_t &nclus,
173                                            Double_t xcenter[],Double_t zcenter[],
174                                            Double_t errxcenter[],
175                                            Double_t errzcenter[],
176                                            Int_t tr1clus[],Int_t tr2clus[],
177                                            Int_t tr3clus[]){
178     // Search for clusters of fired pixels (digits). Two digits are linked
179     // inside a cluster if they are countiguous both in row or column
180     // direction.  Diagonal digits are not linked.
181     // xcenter, ycenter, zcenter are the coordinates of the center
182     // of each found cluster, calculated from the averaging the corresponding
183     // coordinate of the center of the linked digits. The coordinates are
184     // given in the local reference sistem. 
185     // errxcenter, errycenter, errzcenter are the errors associated to
186     // the corresponding average.
187     Int_t   if1, min, max, nd;
188     Int_t   x1, z1, t1, t2, t3, t4;
189     Int_t   ndx, ndz, ndxmin=0, ndxmax=0, ndzmin=0, ndzmax=0;
190     Double_t dx, dz; 
191     Int_t   i,k,ipos=0;
192     Float_t xdum, zdum;      
193     Int_t   kmax, sigmax;
194     Double_t deltax, deltaz;
195     Double_t ndig;
196     Int_t   numberd = 10000;
197     Int_t   *ifpad  = new Int_t[numberd];
198     Int_t   *xpad   = new Int_t[numberd];
199     Int_t   *zpad   = new Int_t[numberd];
200     Int_t   *tr1pad = new Int_t[numberd];
201     Int_t   *tr2pad = new Int_t[numberd];
202     Int_t   *tr3pad = new Int_t[numberd];
203     Int_t   *tr4pad = new Int_t[numberd];
204     Int_t   *iclus  = new Int_t[numberd];
205
206     nclus=1;
207     for (i=0; i < ndigits ; i++){
208         ifpad[i] = -1;
209         iclus[i] = 0;
210     } // end for i
211
212     ifpad[0]=0;
213     for (i=0; i < ndigits-1 ; i++) {
214         if ( ifpad[i] == -1 ) { 
215             nclus++;
216             ipos++;
217             ifpad[i]=nclus-1;
218         } // end if ipad[i]
219         for (Int_t j=i+1 ; j < ndigits ; j++)  {  
220             if (ifpad[j]== -1 ) {
221                 dx = TMath::Abs(digx[i]-digx[j]);
222                 dz = TMath::Abs(digz[i]-digz[j]);
223                 // if ( ( dx+dz )==1 )  //clusters are not diagonal
224                 if(( dx+dz )==1 || (dx==1 && dz==1)){
225                     //diagonal clusters allowed
226                     ipos++;
227                     ifpad[j] = ifpad[i];
228                     
229                     x1         = digx[j];
230                     z1         = digz[j];
231                     digx[j]    = digx[ipos];
232                     digz[j]    = digz[ipos];
233                     digx[ipos] = x1;
234                     digz[ipos] = z1;
235                     
236                     t1 = digtr1[j];
237                     t2 = digtr2[j];
238                     t3 = digtr3[j];
239                     t4 = digtr4[j];
240                     digtr1[j] = digtr1[ipos];
241                     digtr2[j] = digtr2[ipos];
242                     digtr3[j] = digtr3[ipos];
243                     digtr4[j] = digtr4[ipos];
244                     digtr1[ipos] = t1;
245                     digtr2[ipos] = t2;
246                     digtr3[ipos] = t3;
247                     digtr4[ipos] = t4;
248                     
249                     if1 = ifpad[j];
250                     ifpad[j] = ifpad[ipos];
251                     ifpad[ipos] = if1;
252                 } // end dx+dx...
253             }// end if ifpad[j]== -1 
254         } // end for j
255     }//end loop on digits
256     if ( ifpad[ndigits-1] == -1 ) {
257         nclus++;
258         ifpad[ndigits-1]=nclus-1;
259     } // end if ifpad[ndigits-1] == -1
260     
261     for (i=0 ; i < ndigits ; i++) iclus[ifpad[i]]++;
262
263     min=0;
264     max=0;
265     // loop on found clusters 
266     for (i=0 ; i < nclus ; i++){
267         min = max;
268         max += iclus[i];
269         deltax = GetSeg()->Dpx(0);
270         if (iclus[i]!=1){
271             //cluster with more than one digit
272             nd=iclus[i];
273             ndig=(Double_t) nd;
274             Int_t count=0;
275             for (k=min;k<min+nd;k++){
276                 xpad[count] = digx[k];     
277                 zpad[count] = digz[k];
278
279                 tr1pad[count] = digtr1[k];         
280                 tr2pad[count] = digtr2[k];         
281                 tr3pad[count] = digtr3[k];         
282                 tr4pad[count] = digtr4[k];         
283                 
284                 count++; 
285             } // end for k
286             ndxmin = xpad[TMath::LocMin(nd,xpad)];
287             ndxmax = xpad[TMath::LocMax(nd,xpad)];
288             ndzmin = zpad[TMath::LocMin(nd,zpad)];
289             ndzmax = zpad[TMath::LocMax(nd,zpad)];
290             ndx = ndxmax - ndxmin+1;
291             ndz = ndzmax - ndzmin+1;
292             
293             // calculate x and z coordinates of the center of the cluster
294             GetSeg()->GetPadCxz(digx[min],digz[min]-1,xdum, zdum);
295             
296             if (ndx == 1) {
297                 xcenter[i] = xdum;
298             }else{ 
299                 xcenter[i] = 0.;
300                 for (k=0;k<nd;k++) {
301                     GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
302                     xcenter[i] += (xdum / nd);
303                 } // end for k               
304             } // end if ndx
305             
306             if (ndz == 1) {
307                 zcenter[i] = zdum;
308             } else {
309                 zcenter[i] = 0.;
310                 for (k=0;k<nd;k++) {          
311                     GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
312                     zcenter[i] += (zdum / nd);
313                 } // end for k
314             } // end if ndz
315             
316             // error on points in x and z directions
317             
318             if (ndx == 1) {
319                 errxcenter[i] = deltax / TMath::Sqrt(12.);
320             } else {
321                 errxcenter[i] = 0.;                     
322                 for (k=0;k<nd;k++){ 
323                     GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
324                     errxcenter[i] += ((xdum-xcenter[i])*(xdum-xcenter[i]))/
325                         (nd*(nd-1)); 
326                 } // end for k
327                 errxcenter[i] = TMath::Sqrt(errxcenter[i]);
328             } // end if ndx     
329             if (ndz == 1) {
330                 deltaz = GetSeg()->Dpz(digz[min]);                    
331                 errzcenter[i] = deltaz / TMath::Sqrt(12.);
332             } else {
333                 errzcenter[i] = 0.;
334                 for (k=0;k<nd;k++){ 
335                     GetSeg()->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
336                     errzcenter[i] += ((zdum-zcenter[i])*(zdum-zcenter[i]))/
337                         (nd*(nd-1));
338                 } // end for k
339                 errzcenter[i] = TMath::Sqrt(errzcenter[i]);
340             } // end if ndz
341             // take three track numbers for the cluster
342             // choose the track numbers of the digit with higher signal 
343             kmax = 0;
344             sigmax = 0;
345             for (k=0;k<nd;k++){
346                 if(tr4pad[k] > sigmax){
347                     sigmax = tr4pad[k];
348                     kmax   = k;
349                 } // end if tr4pad[k]
350             } // end for k
351             if(sigmax != 0) {
352                 tr1clus[i]= tr1pad[kmax];
353                 tr2clus[i]= tr2pad[kmax];
354                 tr3clus[i]= tr3pad[kmax];
355             } else {
356                 tr1clus[i]= -2;
357                 tr2clus[i]= -2;
358                 tr3clus[i]= -2;
359             } // end if sigmax
360         } else {
361             // cluster with single digit
362             ndig= 1.;
363             ndx = 1;
364             ndz = 1;
365             GetSeg()->GetPadCxz(digx[min],digz[min]-1,xdum,zdum);
366             xcenter[i] = xdum;
367             zcenter[i] = zdum;
368             tr1clus[i]=digtr1[min];
369             tr2clus[i]=digtr2[min];
370             tr3clus[i]=digtr3[min];
371             deltaz = GetSeg()->Dpz(digz[min]);
372             errxcenter[i] = deltax / TMath::Sqrt(12.);
373             errzcenter[i] = deltaz / TMath::Sqrt(12.);
374         } // end if iclus[i]
375
376         // store the cluster information to the AliITSRawCLusterSPD object
377    
378
379         //put the cluster center in local reference frame of the detector
380         // and in microns
381         xcenter[i] = xcenter[i] - GetSeg()->Dx()/2.; 
382         zcenter[i] = zcenter[i] - GetSeg()->Dz()/2.;
383
384         AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(zcenter[i], //f
385                                                              xcenter[i], //f
386                                                              ndig, //f
387                                                              ndz,ndx, //ii
388                                                              ndxmin,ndxmax,//ii
389                                                              (Double_t) ndzmin,
390                                                              (Double_t) ndzmax,
391                                                              0,GetModule());
392         fITS->AddCluster(0,clust);
393         delete clust;
394     }//end loop on clusters   
395     delete[] ifpad;
396     delete[] xpad ;
397     delete[] zpad ;
398     delete[] iclus;
399     delete[] tr1pad;
400     delete[] tr2pad;
401     delete[] tr3pad;
402     delete[] tr4pad;
403 }
404 //______________________________________________________----------------
405 void AliITSClusterFinderSPD::DigitToPoint(Int_t nclus,
406                                           Double_t *xcenter,Double_t *zcenter,
407                                           Double_t *errxcenter,
408                                           Double_t *errzcenter, 
409                                           Int_t *tr1clus, Int_t *tr2clus,
410                                           Int_t *tr3clus){
411     // A point is associated to each cluster of SPD digits. The points
412     // and their associated errors are stored in the file galiceSP.root.
413     Double_t l[3],xg,zg;
414     const Double_t kconv = 1.0e-4; // micron -> cm
415
416     // get rec points
417     for (Int_t i=0; i<nclus; i++){
418         l[0] = kconv*xcenter[i];
419         l[1] = kconv*GetSeg()->Dy()/2.;
420         l[2] = kconv*zcenter[i];
421
422         xg = l[0]; 
423         zg = l[2]; 
424
425         Double_t sigma2x = (kconv*errxcenter[i]) * (kconv*errxcenter[i]);
426         Double_t sigma2z = (kconv*errzcenter[i]) * (kconv*errzcenter[i]);
427         AliITSRecPoint rnew;
428         rnew.SetX(xg);
429         rnew.SetZ(zg);
430         rnew.SetQ(1.);
431         rnew.SetdEdX(0.);
432         rnew.SetSigmaX2(sigma2x);
433         rnew.SetSigmaZ2(sigma2z);
434         rnew.fTracks[0]=tr1clus[i];
435         rnew.fTracks[1]=tr2clus[i];
436         rnew.fTracks[2]=tr3clus[i];
437         fITS->AddRecPoint(rnew); 
438     } // end for i
439 }