]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSClusterFinderSPDbari.cxx
Fixed bug with HP and no unget in iostream.h. Now using putback instead.
[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;
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
0aedb77e 134 DigitToPoint(nclus,xcenterl,zcenterl,errxcenter,errzcenter,
a82b8484 135 tr1clus, tr2clus, tr3clus);
0aedb77e 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//-----------------------------------------------------------------
153void 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;
8a838ecf 176 Int_t i,k,ipos=0;
0aedb77e 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 }
8a838ecf 243 for (i=0 ; i < ndigits ; i++) iclus[ifpad[i]]++;
0aedb77e 244
245 min=0;
246 max=0;
247
248 // loop on found clusters
249
8a838ecf 250 for (i=0 ; i < nclus ; i++)
0aedb77e 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;
8a838ecf 260 for (k=min;k<min+nd;k++)
0aedb77e 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.;
8a838ecf 287 for (k=0;k<nd;k++) {
0aedb77e 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.;
8a838ecf 298 for (k=0;k<nd;k++) {
0aedb77e 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.;
8a838ecf 311 for (k=0;k<nd;k++){
0aedb77e 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.;
8a838ecf 324 for (k=0;k<nd;k++){
0aedb77e 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
8a838ecf 332 for (k=0;k<nd;k++){
0aedb77e 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//______________________________________________________
374void AliITSClusterFinderSPDbari::DigitToPoint(Int_t nclus,
375 Float_t *xcenter, Float_t *zcenter,
376 Float_t *errxcenter,Float_t *errzcenter,
a82b8484 377 Int_t *tr1clus, Int_t *tr2clus, Int_t *tr3clus){
0aedb77e 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
6ad50b01 384 Float_t l[3],xg,zg;
0aedb77e 385 const Float_t kconv = 1.0e-4; // micron -> cm
386
387 // get rec points
388 AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
389
0aedb77e 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];
a82b8484 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]);
0aedb77e 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}