c47621d5791b070054630a26923129247475ad40
[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=10000;
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     
94     //  output of Cluster Finder    
95     Int_t numberc=1000;
96     Float_t *xcenterl   = new Float_t[numberc];
97     Float_t *zcenterl   = new Float_t[numberc];
98     Float_t *errxcenter = new Float_t[numberc];
99     Float_t *errzcenter = new Float_t[numberc];
100     Int_t   *tr1clus    = new Int_t[numberc];
101     Int_t   *tr2clus    = new Int_t[numberc];
102     Int_t   *tr3clus    = new Int_t[numberc];
103
104     Int_t nclus;
105
106     digitcount=0;
107     Int_t ndigits = fDigits->GetEntriesFast();  
108     if (!ndigits) return;
109
110
111     AliITSdigit *dig;
112     AliITSdigitSPD *dig1;
113     Int_t ndig;
114     for(ndig=0; ndig<ndigits; ndig++) {
115          dig= (AliITSdigit*)fDigits->UncheckedAt(ndig);
116      
117          digx[digitcount] = dig->fCoord2+1;  //starts at 1
118          digz[digitcount] = dig->fCoord1+1;  //starts at 1
119
120          dig1= (AliITSdigitSPD*)fDigits->UncheckedAt(ndig);
121
122          digtr1[digitcount] = dig1->fTracks[0];
123          digtr2[digitcount] = dig1->fTracks[1];
124          digtr3[digitcount] = dig1->fTracks[2];
125
126          digitcount++;
127     }
128
129
130         ClusterFinder(digitcount,digx,digz,digtr1,digtr2,digtr3,
131               nclus,xcenterl,zcenterl,errxcenter,errzcenter,
132               tr1clus, tr2clus, tr3clus);
133  
134         DigitToPoint(nclus,xcenterl,zcenterl,errxcenter,errzcenter,
135               tr1clus, tr2clus, tr3clus);
136
137
138   delete[] digx       ;
139   delete[] digz       ;
140   delete[] digtr1     ;
141   delete[] digtr2     ;
142   delete[] digtr3     ;
143   delete[] xcenterl   ;
144   delete[] zcenterl   ;
145   delete[] errxcenter ;
146   delete[] errzcenter ;
147   delete[] tr1clus    ;
148   delete[] tr2clus    ;
149   delete[] tr3clus    ;
150   
151 }
152 //-----------------------------------------------------------------
153 void AliITSClusterFinderSPDbari::ClusterFinder(Int_t ndigits,
154                     Int_t digx[],Int_t digz[],
155                     Int_t digtr1[],Int_t digtr2[],Int_t digtr3[],
156                     Int_t &nclus, Float_t xcenter[],Float_t zcenter[],
157                                     Float_t errxcenter[],Float_t errzcenter[],
158                     Int_t tr1clus[], Int_t tr2clus[], Int_t tr3clus[]) {
159 //
160 // Search for clusters of fired pixels (digits). Two digits are linked
161 // inside a cluster if they are countiguous both in row or column
162 // direction.  Diagonal digits are not linked.
163 // xcenter, ycenter, zcenter are the coordinates of the center
164 // of each found cluster, calculated from the averaging the corresponding
165 // coordinate of the center of the linked digits. The coordinates are
166 // given in the local reference sistem. 
167 // errxcenter, errycenter, errzcenter are the errors associated to
168 // the corresponding average.
169 //
170 //
171
172   Int_t if1, min, max, nd;
173   Int_t x1, z1, t1, t2, t3;
174   Int_t ndx, ndz, ndxmin, ndxmax, ndzmin, ndzmax;
175   Float_t dx, dz; 
176   Int_t i,k,ipos=0;
177   Float_t xdum, zdum;      
178
179   Int_t numberd=10000;
180   Int_t *ifpad = new Int_t[numberd];
181   Int_t *xpad  = new Int_t[numberd];
182   Int_t *zpad  = new Int_t[numberd];
183   Int_t *tr1pad  = new Int_t[numberd];
184   Int_t *tr2pad  = new Int_t[numberd];
185   Int_t *tr3pad  = new Int_t[numberd];
186   Int_t *iclus   = new Int_t[numberd];
187
188  
189   nclus=1;
190   for (i=0; i < ndigits ; i++) ifpad[i] = -1;
191
192   ifpad[0]=0;
193   for (i=0; i < ndigits-1 ; i++) 
194   {
195     if ( ifpad[i] == -1 ) 
196     { 
197                 nclus++;
198                 ipos++;
199         ifpad[i]=nclus-1;
200     }
201     for (Int_t j=i+1 ; j < ndigits ; j++)  
202     {  
203       if (ifpad[j]== -1 )
204       {
205              dx = TMath::Abs(digx[i]-digx[j]);
206              dz = TMath::Abs(digz[i]-digz[j]);
207
208 //           if ( ( dx+dz )==1 )  //  clusters are not diagonal
209              if ( ( dx+dz )==1 || (dx==1 && dz==1) )  //  diagonal clusters allowed
210              {
211                     ipos++;
212                     ifpad[j]=ifpad[i];
213
214                     x1=digx[j];
215                     z1=digz[j];
216                     digx[j]=digx[ipos];
217                     digz[j]=digz[ipos];
218                     digx[ipos]=x1;
219                     digz[ipos]=z1;
220
221                     t1=digtr1[j];
222                     t2=digtr2[j];
223                     t3=digtr3[j];
224                     digtr1[j]=digtr1[ipos];
225                     digtr2[j]=digtr2[ipos];
226                     digtr3[j]=digtr3[ipos];
227                     digtr1[ipos]=t1;
228                     digtr2[ipos]=t2;
229                     digtr3[ipos]=t3;
230
231                     if1=ifpad[j];
232                     ifpad[j]=ifpad[ipos];
233                     ifpad[ipos]=if1;
234              }
235       }
236     }
237    }   
238    if ( ifpad[ndigits-1] == -1 )
239    {
240           nclus++;
241           ifpad[ndigits-1]=nclus-1;
242    }
243    for (i=0 ; i < ndigits ; i++) iclus[ifpad[i]]++;
244
245    min=0;
246    max=0;
247
248    // loop on found clusters 
249
250    for (i=0 ; i < nclus ; i++)  
251    {
252       min = max;
253       max += iclus[i];
254       Float_t deltax = fSegmentation->Dpx(0);
255       if (iclus[i]!=1) 
256       {
257         //cluster with more than one digit
258         nd=iclus[i];
259             Int_t count=0;
260         for (k=min;k<min+nd;k++)
261         {
262                xpad[count] = digx[k];      
263                zpad[count] = digz[k];
264
265                tr1pad[count] = digtr1[k];          
266                tr2pad[count] = digtr2[k];          
267                tr3pad[count] = digtr3[k];          
268
269                count++; 
270         }
271         ndxmin = xpad[TMath::LocMin(nd,xpad)];
272         ndxmax = xpad[TMath::LocMax(nd,xpad)];
273         ndzmin = zpad[TMath::LocMin(nd,zpad)];
274         ndzmax = zpad[TMath::LocMax(nd,zpad)];
275         ndx = ndxmax - ndxmin+1;
276         ndz = ndzmax - ndzmin+1;
277
278         // calculate x and z coordinates of the center of the cluster
279         fSegmentation->GetPadCxz(digx[min],digz[min]-1,xdum, zdum);
280
281
282         if (ndx == 1) {     
283             xcenter[i] = xdum;
284         }    
285         else{ 
286            xcenter[i] = 0.;
287            for (k=0;k<nd;k++) {
288              fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
289              xcenter[i] += (xdum / nd);
290            }                   
291         }
292
293         if (ndz == 1) {
294             zcenter[i] = zdum;
295         }   
296         else {
297            zcenter[i] = 0.;
298            for (k=0;k<nd;k++) {       
299              fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
300              zcenter[i] += (zdum / nd);
301            }
302         }
303
304         // error on points in x and z directions
305
306         if (ndx == 1) {
307              errxcenter[i] = deltax / TMath::Sqrt(12.);
308         }
309         else {
310              errxcenter[i] = 0.;                        
311              for (k=0;k<nd;k++){ 
312                fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
313                errxcenter[i] += ((xdum-xcenter[i])*(xdum-xcenter[i]))/(nd*(nd-1)); 
314              }   
315              errxcenter[i] = TMath::Sqrt(errxcenter[i]);
316         }
317         
318         if (ndz == 1) {
319             Float_t deltaz = fSegmentation->Dpz(digz[min]);                   
320             errzcenter[i] = deltaz / TMath::Sqrt(12.);
321         }
322         else {
323              errzcenter[i] = 0.;
324              for (k=0;k<nd;k++){ 
325                fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
326                errzcenter[i] += ((zdum-zcenter[i])*(zdum-zcenter[i]))/(nd*(nd-1));
327              }
328              errzcenter[i] = TMath::Sqrt(errzcenter[i]);
329         }    
330
331         // take three track numbers for the cluster
332         for (k=0;k<nd;k++){
333           if(tr1pad[k] != -2) tr1clus[i]=tr1pad[k];
334           if(tr2pad[k] != -2) tr2clus[i]=tr2pad[k];
335           if(tr3pad[k] != -2) tr3clus[i]=tr3pad[k];
336         }
337         if(tr1clus[i] == 0) tr1clus[i]= -2;
338         if(tr2clus[i] == 0) tr2clus[i]= -2;
339         if(tr3clus[i] == 0) tr3clus[i]= -2;
340       }
341       else  {
342       
343         // cluster with single digit
344         ndx = 1;
345         ndz = 1;
346         fSegmentation->GetPadCxz(digx[min],digz[min]-1,xdum,zdum);
347         xcenter[i] = xdum;
348         zcenter[i] = zdum;
349         tr1clus[i]=digtr1[min];
350         tr2clus[i]=digtr2[min];
351         tr3clus[i]=digtr3[min];
352             Float_t deltaz = fSegmentation->Dpz(digz[min]);
353             errxcenter[i] = deltax / TMath::Sqrt(12.);
354             errzcenter[i] = deltaz / TMath::Sqrt(12.);
355      }
356
357      // store the cluster information to the AliITSRawCLusterSPD object
358      AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
359
360      //put the cluster center in local reference frame of the detector
361      // and in microns
362      xcenter[i] = xcenter[i] - fSegmentation->Dx()/2.; 
363      zcenter[i] = zcenter[i] - fSegmentation->Dz()/2.;
364
365      AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(zcenter[i],xcenter[i],1.,ndz,ndx,0.,0.,0.,0.,0.,0.,0.);
366      iTS->AddCluster(0,clust);
367    }      
368    delete[] ifpad;
369    delete[] xpad ;
370    delete[] zpad ;
371    delete[] iclus;
372 }
373 //______________________________________________________
374 void AliITSClusterFinderSPDbari::DigitToPoint(Int_t nclus,
375                    Float_t *xcenter,   Float_t *zcenter,
376                    Float_t *errxcenter,Float_t *errzcenter, 
377                    Int_t *tr1clus, Int_t *tr2clus, Int_t *tr3clus){
378  //
379  // A point is associated to each cluster of SPD digits. The points
380  // and their associated errors are stored in the file galiceSP.root.
381  //
382  //
383  
384      Float_t l[3],xg,zg;
385      const Float_t kconv = 1.0e-4; // micron -> cm
386
387      // get rec points
388      AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
389
390      for (Int_t i=0; i<nclus; i++)
391      {
392         l[0] = kconv*xcenter[i];
393         l[1] = kconv*fSegmentation->Dy()/2.;
394         l[2] = kconv*zcenter[i];
395
396         xg = l[0]; 
397         zg = l[2]; 
398
399             Float_t sigma2x = (kconv*errxcenter[i]) * (kconv*errxcenter[i]);
400             Float_t sigma2z = (kconv*errzcenter[i]) * (kconv*errzcenter[i]);
401         AliITSRecPoint rnew;
402         rnew.SetX(xg);
403         rnew.SetZ(zg);
404         rnew.SetQ(1.);
405         rnew.SetdEdX(0.);
406         rnew.SetSigmaX2(sigma2x);
407         rnew.SetSigmaZ2(sigma2z);
408         rnew.fTracks[0]=tr1clus[i];
409         rnew.fTracks[1]=tr2clus[i];
410         rnew.fTracks[2]=tr3clus[i];
411         iTS->AddRecPoint(rnew); 
412      }
413
414 }