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