New version compatible with the current HEAD
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSPDbari.cxx
CommitLineData
0aedb77e 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
30ClassImp(AliITSClusterFinderSPDbari)
31
32//----------------------------------------------------------
33AliITSClusterFinderSPDbari::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//_____________________________________________________________________________
46AliITSClusterFinderSPDbari::AliITSClusterFinderSPDbari()
47{
48 // constructor
49 fSegmentation=0;
50 fDigits=0;
51 fClusters=0;
52 fNclusters=0;
53 SetDx();
54 SetDz();
55
56}
57
58//__________________________________________________________________________
59AliITSClusterFinderSPDbari::AliITSClusterFinderSPDbari(const
60AliITSClusterFinderSPDbari &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//_________________________________________________________________________
71AliITSClusterFinderSPDbari&
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//_____________________________________________________________________________
83void AliITSClusterFinderSPDbari::FindRawClusters(){
0aedb77e 84
85 // input of Cluster Finder
86 Int_t digitcount=0;
426edb57 87 Int_t numberd=100000;
0aedb77e 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];
426edb57 93 Int_t *digtr4 = new Int_t[numberd];
0aedb77e 94
95 // output of Cluster Finder
426edb57 96 Int_t numberc=10000;
0aedb77e 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];
426edb57 126 digtr4[digitcount] = dig1->fSignal;
0aedb77e 127
128 digitcount++;
129 }
130
131
426edb57 132 ClusterFinder(digitcount,digx,digz,digtr1,digtr2,digtr3,digtr4,
0aedb77e 133 nclus,xcenterl,zcenterl,errxcenter,errzcenter,
134 tr1clus, tr2clus, tr3clus);
135
426edb57 136 DigitToPoint(nclus,xcenterl,zcenterl,errxcenter,errzcenter,
a82b8484 137 tr1clus, tr2clus, tr3clus);
0aedb77e 138
139
426edb57 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
0aedb77e 154}
155//-----------------------------------------------------------------
156void AliITSClusterFinderSPDbari::ClusterFinder(Int_t ndigits,
157 Int_t digx[],Int_t digz[],
426edb57 158 Int_t digtr1[],Int_t digtr2[],Int_t digtr3[],Int_t digtr4[],
0aedb77e 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;
426edb57 176 Int_t x1, z1, t1, t2, t3, t4;
0aedb77e 177 Int_t ndx, ndz, ndxmin, ndxmax, ndzmin, ndzmax;
178 Float_t dx, dz;
8a838ecf 179 Int_t i,k,ipos=0;
0aedb77e 180 Float_t xdum, zdum;
426edb57 181 Int_t kmax, sigmax;
182 Float_t deltax, deltaz;
183 Float_t ndig;
0aedb77e 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];
426edb57 192 Int_t *tr4pad = new Int_t[numberd];
0aedb77e 193 Int_t *iclus = new Int_t[numberd];
194
0aedb77e 195 nclus=1;
426edb57 196 for (i=0; i < ndigits ; i++){
197 ifpad[i] = -1;
198 iclus[i] = 0;
199 }
0aedb77e 200
201 ifpad[0]=0;
426edb57 202 for (i=0; i < ndigits-1 ; i++) {
203 if ( ifpad[i] == -1 ) {
0aedb77e 204 nclus++;
205 ipos++;
206 ifpad[i]=nclus-1;
207 }
426edb57 208 for (Int_t j=i+1 ; j < ndigits ; j++) {
209 if (ifpad[j]== -1 ) {
0aedb77e 210 dx = TMath::Abs(digx[i]-digx[j]);
211 dz = TMath::Abs(digz[i]-digz[j]);
212
426edb57 213// if ( ( dx+dz )==1 ) //clusters are not diagonal
214 if ( ( dx+dz )==1 || (dx==1 && dz==1) ) { //diagonal clusters allowed
0aedb77e 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];
426edb57 228 t4=digtr4[j];
0aedb77e 229 digtr1[j]=digtr1[ipos];
230 digtr2[j]=digtr2[ipos];
231 digtr3[j]=digtr3[ipos];
426edb57 232 digtr4[j]=digtr4[ipos];
0aedb77e 233 digtr1[ipos]=t1;
234 digtr2[ipos]=t2;
235 digtr3[ipos]=t3;
426edb57 236 digtr4[ipos]=t4;
0aedb77e 237
238 if1=ifpad[j];
239 ifpad[j]=ifpad[ipos];
240 ifpad[ipos]=if1;
241 }
242 }
243 }
426edb57 244 }//end loop on digits
245
246 if ( ifpad[ndigits-1] == -1 ) {
0aedb77e 247 nclus++;
248 ifpad[ndigits-1]=nclus-1;
426edb57 249 }
0aedb77e 250
426edb57 251 for (i=0 ; i < ndigits ; i++) iclus[ifpad[i]]++;
0aedb77e 252
426edb57 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 {
0aedb77e 263 //cluster with more than one digit
264 nd=iclus[i];
426edb57 265 ndig=(Float_t) nd;
0aedb77e 266 Int_t count=0;
8a838ecf 267 for (k=min;k<min+nd;k++)
0aedb77e 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];
426edb57 275 tr4pad[count] = digtr4[k];
0aedb77e 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
426edb57 286
0aedb77e 287 // calculate x and z coordinates of the center of the cluster
288 fSegmentation->GetPadCxz(digx[min],digz[min]-1,xdum, zdum);
289
426edb57 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 }
0aedb77e 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.;
8a838ecf 319 for (k=0;k<nd;k++){
0aedb77e 320 fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
321 errxcenter[i] += ((xdum-xcenter[i])*(xdum-xcenter[i]))/(nd*(nd-1));
322 }
426edb57 323 errxcenter[i] = TMath::Sqrt(errxcenter[i]);
0aedb77e 324 }
325
426edb57 326 if (ndz == 1) {
327 deltaz = fSegmentation->Dpz(digz[min]);
328 errzcenter[i] = deltaz / TMath::Sqrt(12.);
0aedb77e 329 }
426edb57 330 else {
331 errzcenter[i] = 0.;
332 for (k=0;k<nd;k++){
0aedb77e 333 fSegmentation->GetPadCxz(xpad[k],zpad[k]-1,xdum,zdum);
426edb57 334 errzcenter[i] += ((zdum-zcenter[i])*(zdum-zcenter[i]))/(nd*(nd-1));
335 }
336 errzcenter[i] = TMath::Sqrt(errzcenter[i]);
337 }
0aedb77e 338 // take three track numbers for the cluster
426edb57 339 // choose the track numbers of the digit with higher signal
340 kmax = 0;
341 sigmax = 0;
8a838ecf 342 for (k=0;k<nd;k++){
426edb57 343 if(tr4pad[k] > sigmax){
344 sigmax = tr4pad[k];
345 kmax = k;
346 }
0aedb77e 347 }
426edb57 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 {
0aedb77e 360
361 // cluster with single digit
426edb57 362 ndig= 1.;
363 ndx = 1;
0aedb77e 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];
426edb57 371 deltaz = fSegmentation->Dpz(digz[min]);
0aedb77e 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
426edb57 384
385 AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(zcenter[i],xcenter[i],ndig,ndz,ndx,0.,0.,0.,0.,0.,0.,0.);
0aedb77e 386 iTS->AddCluster(0,clust);
a0d3be87 387 delete clust;
426edb57 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;
0aedb77e 397}
398//______________________________________________________
399void AliITSClusterFinderSPDbari::DigitToPoint(Int_t nclus,
400 Float_t *xcenter, Float_t *zcenter,
401 Float_t *errxcenter,Float_t *errzcenter,
a82b8484 402 Int_t *tr1clus, Int_t *tr2clus, Int_t *tr3clus){
0aedb77e 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
6ad50b01 409 Float_t l[3],xg,zg;
0aedb77e 410 const Float_t kconv = 1.0e-4; // micron -> cm
411
412 // get rec points
413 AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
414
0aedb77e 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];
a82b8484 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]);
0aedb77e 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 }
0aedb77e 438}