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