]>
Commit | Line | Data |
---|---|---|
1c53abe2 | 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 | ||
88cb7938 | 16 | /* $Id$ */ |
1c53abe2 | 17 | |
18 | //------------------------------------------------------- | |
19 | // Implementation of the TPC clusterer | |
20 | // | |
21 | // Origin: Marian Ivanov | |
22 | //------------------------------------------------------- | |
23 | ||
24 | #include "AliTPCclustererMI.h" | |
25 | #include "AliTPCclusterMI.h" | |
26 | #include <TObjArray.h> | |
27 | #include <TFile.h> | |
28 | #include "AliTPCClustersArray.h" | |
29 | #include "AliTPCClustersRow.h" | |
f8aae377 | 30 | #include "AliTPCRawStream.h" |
1c53abe2 | 31 | #include "AliDigits.h" |
32 | #include "AliSimDigits.h" | |
33 | #include "AliTPCParam.h" | |
f8aae377 | 34 | #include "AliRawReader.h" |
35 | #include "AliTPCRawStream.h" | |
36 | #include "AliRunLoader.h" | |
37 | #include "AliLoader.h" | |
cc5e9db0 | 38 | #include "Riostream.h" |
1c53abe2 | 39 | #include <TTree.h> |
40 | ||
41 | ClassImp(AliTPCclustererMI) | |
42 | ||
43 | ||
44 | ||
f8aae377 | 45 | AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par) |
1c53abe2 | 46 | { |
47 | fInput =0; | |
48 | fOutput=0; | |
f8aae377 | 49 | fParam = par; |
1c53abe2 | 50 | } |
51 | void AliTPCclustererMI::SetInput(TTree * tree) | |
52 | { | |
53 | // | |
54 | // set input tree with digits | |
55 | // | |
56 | fInput = tree; | |
57 | if (!fInput->GetBranch("Segment")){ | |
58 | cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n"; | |
59 | fInput=0; | |
60 | return; | |
61 | } | |
62 | } | |
63 | ||
64 | void AliTPCclustererMI::SetOutput(TTree * tree) | |
65 | { | |
66 | // | |
67 | // | |
68 | fOutput= tree; | |
69 | AliTPCClustersRow clrow; | |
70 | AliTPCClustersRow *pclrow=&clrow; | |
71 | clrow.SetClass("AliTPCclusterMI"); | |
72 | clrow.SetArray(1); // to make Clones array | |
73 | fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200); | |
74 | } | |
75 | ||
76 | ||
77 | Float_t AliTPCclustererMI::GetSigmaY2(Int_t iz){ | |
78 | // sigma y2 = in digits - we don't know the angle | |
79 | Float_t z = iz*fParam->GetZWidth(); | |
80 | Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/ | |
81 | (fPadWidth*fPadWidth); | |
82 | Float_t sres = 0.25; | |
83 | Float_t res = sd2+sres; | |
84 | return res; | |
85 | } | |
86 | ||
87 | ||
88 | Float_t AliTPCclustererMI::GetSigmaZ2(Int_t iz){ | |
89 | //sigma z2 = in digits - angle estimated supposing vertex constraint | |
90 | Float_t z = iz*fZWidth; | |
91 | Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth); | |
92 | Float_t angular = fPadLength*(fParam->GetZLength()-z)/(fRx*fZWidth); | |
93 | angular*=angular; | |
94 | angular/=12.; | |
95 | Float_t sres = fParam->GetZSigma()/fZWidth; | |
96 | sres *=sres; | |
97 | Float_t res = angular +sd2+sres; | |
98 | return res; | |
99 | } | |
100 | ||
176aff27 | 101 | void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Int_t *bins, UInt_t /*m*/, |
1c53abe2 | 102 | AliTPCclusterMI &c) |
103 | { | |
104 | Int_t i0=k/max; //central pad | |
105 | Int_t j0=k%max; //central time bin | |
106 | ||
107 | // set pointers to data | |
108 | //Int_t dummy[5] ={0,0,0,0,0}; | |
109 | Int_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2 | |
110 | Int_t * resmatrix[5]; | |
111 | for (Int_t di=-2;di<=2;di++){ | |
112 | matrix[di+2] = &bins[k+di*max]; | |
113 | resmatrix[di+2] = &fResBins[k+di*max]; | |
114 | } | |
115 | //build matrix with virtual charge | |
116 | Float_t sigmay2= GetSigmaY2(j0); | |
117 | Float_t sigmaz2= GetSigmaZ2(j0); | |
118 | ||
119 | Float_t vmatrix[5][5]; | |
120 | vmatrix[2][2] = matrix[2][0]; | |
121 | c.SetType(0); | |
122 | c.SetMax(Short_t(vmatrix[2][2])); // write maximal amplitude | |
123 | for (Int_t di =-1;di <=1;di++) | |
124 | for (Int_t dj =-1;dj <=1;dj++){ | |
125 | Float_t amp = matrix[di+2][dj]; | |
126 | if ( (amp<2) && (fLoop<2)){ | |
127 | // if under threshold - calculate virtual charge | |
128 | Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2); | |
129 | amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio; | |
130 | if (amp>2) amp = 2; | |
131 | vmatrix[2+di][2+dj]=amp; | |
132 | vmatrix[2+2*di][2+2*dj]=0; | |
133 | if ( (di*dj)!=0){ | |
134 | //DIAGONAL ELEMENTS | |
135 | vmatrix[2+2*di][2+dj] =0; | |
136 | vmatrix[2+di][2+2*dj] =0; | |
137 | } | |
138 | continue; | |
139 | } | |
140 | if (amp<4){ | |
141 | //if small amplitude - below 2 x threshold - don't consider other one | |
142 | vmatrix[2+di][2+dj]=amp; | |
143 | vmatrix[2+2*di][2+2*dj]=0; // don't take to the account next bin | |
144 | if ( (di*dj)!=0){ | |
145 | //DIAGONAL ELEMENTS | |
146 | vmatrix[2+2*di][2+dj] =0; | |
147 | vmatrix[2+di][2+2*dj] =0; | |
148 | } | |
149 | continue; | |
150 | } | |
151 | //if bigger then take everything | |
152 | vmatrix[2+di][2+dj]=amp; | |
153 | vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ; | |
154 | if ( (di*dj)!=0){ | |
155 | //DIAGONAL ELEMENTS | |
156 | vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj]; | |
157 | vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2]; | |
158 | } | |
159 | } | |
160 | ||
161 | ||
162 | ||
163 | Float_t sumw=0; | |
164 | Float_t sumiw=0; | |
165 | Float_t sumi2w=0; | |
166 | Float_t sumjw=0; | |
167 | Float_t sumj2w=0; | |
168 | // | |
169 | for (Int_t i=-2;i<=2;i++) | |
170 | for (Int_t j=-2;j<=2;j++){ | |
171 | Float_t amp = vmatrix[i+2][j+2]; | |
172 | ||
173 | sumw += amp; | |
174 | sumiw += i*amp; | |
175 | sumi2w += i*i*amp; | |
176 | sumjw += j*amp; | |
177 | sumj2w += j*j*amp; | |
178 | } | |
179 | // | |
180 | Float_t meani = sumiw/sumw; | |
181 | Float_t mi2 = sumi2w/sumw-meani*meani; | |
182 | Float_t meanj = sumjw/sumw; | |
183 | Float_t mj2 = sumj2w/sumw-meanj*meanj; | |
184 | // | |
185 | Float_t ry = mi2/sigmay2; | |
186 | Float_t rz = mj2/sigmaz2; | |
187 | ||
188 | // | |
189 | if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return; | |
190 | if ( (ry <1.2) && (rz<1.2) ) { | |
191 | //if cluster looks like expected | |
192 | //+1.2 deviation from expected sigma accepted | |
193 | // c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2)); | |
194 | ||
195 | meani +=i0; | |
196 | meanj +=j0; | |
197 | //set cluster parameters | |
198 | c.SetQ(sumw); | |
199 | c.SetY(meani*fPadWidth); | |
200 | c.SetZ(meanj*fZWidth); | |
201 | c.SetSigmaY2(mi2); | |
202 | c.SetSigmaZ2(mj2); | |
203 | AddCluster(c); | |
204 | //remove cluster data from data | |
205 | for (Int_t di=-2;di<=2;di++) | |
206 | for (Int_t dj=-2;dj<=2;dj++){ | |
207 | resmatrix[di+2][dj] -= Int_t(vmatrix[di+2][dj+2]); | |
208 | if (resmatrix[di+2][dj]<0) resmatrix[di+2][dj]=0; | |
209 | } | |
210 | resmatrix[2][0] =0; | |
211 | return; | |
212 | } | |
213 | // | |
214 | //unfolding when neccessary | |
215 | // | |
216 | ||
217 | Int_t * matrix2[7]; //7x7 matrix with digits - indexing i = 0 ..6 j = -3..3 | |
218 | Int_t dummy[7]={0,0,0,0,0,0}; | |
219 | for (Int_t di=-3;di<=3;di++){ | |
220 | matrix2[di+3] = &bins[k+di*max]; | |
221 | if ((k+di*max)<3) matrix2[di+3] = &dummy[3]; | |
222 | if ((k+di*max)>fMaxBin-3) matrix2[di+3] = &dummy[3]; | |
223 | } | |
224 | Float_t vmatrix2[5][5]; | |
225 | Float_t sumu; | |
226 | Float_t overlap; | |
227 | UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap); | |
228 | // | |
229 | // c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2)); | |
230 | meani +=i0; | |
231 | meanj +=j0; | |
232 | //set cluster parameters | |
233 | c.SetQ(sumu); | |
234 | c.SetY(meani*fPadWidth); | |
235 | c.SetZ(meanj*fZWidth); | |
236 | c.SetSigmaY2(mi2); | |
237 | c.SetSigmaZ2(mj2); | |
238 | c.SetType(Char_t(overlap)+1); | |
239 | AddCluster(c); | |
240 | ||
241 | //unfolding 2 | |
242 | meani-=i0; | |
243 | meanj-=j0; | |
244 | if (gDebug>4) | |
245 | printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]); | |
246 | } | |
247 | ||
248 | ||
249 | ||
250 | void AliTPCclustererMI::UnfoldCluster(Int_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj, | |
251 | Float_t & sumu, Float_t & overlap ) | |
252 | { | |
253 | // | |
254 | //unfold cluster from input matrix | |
255 | //data corresponding to cluster writen in recmatrix | |
256 | //output meani and meanj | |
257 | ||
258 | //take separatelly y and z | |
259 | ||
260 | Float_t sum3i[7] = {0,0,0,0,0,0,0}; | |
261 | Float_t sum3j[7] = {0,0,0,0,0,0,0}; | |
262 | ||
263 | for (Int_t k =0;k<7;k++) | |
264 | for (Int_t l = -1; l<=1;l++){ | |
265 | sum3i[k]+=matrix2[k][l]; | |
266 | sum3j[k]+=matrix2[l+3][k-3]; | |
267 | } | |
268 | Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}}; | |
269 | // | |
270 | //unfold y | |
271 | Float_t sum3wi = 0; //charge minus overlap | |
272 | Float_t sum3wio = 0; //full charge | |
273 | Float_t sum3iw = 0; //sum for mean value | |
274 | for (Int_t dk=-1;dk<=1;dk++){ | |
275 | sum3wio+=sum3i[dk+3]; | |
276 | if (dk==0){ | |
277 | sum3wi+=sum3i[dk+3]; | |
278 | } | |
279 | else{ | |
280 | Float_t ratio =1; | |
281 | if ( ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))|| | |
282 | sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){ | |
283 | Float_t xm2 = sum3i[-dk+3]; | |
284 | Float_t xm1 = sum3i[+3]; | |
285 | Float_t x1 = sum3i[2*dk+3]; | |
286 | Float_t x2 = sum3i[3*dk+3]; | |
cc5e9db0 | 287 | Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001); |
288 | Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.); | |
1c53abe2 | 289 | ratio = w11/(w11+w12); |
290 | for (Int_t dl=-1;dl<=1;dl++) | |
291 | mratio[dk+1][dl+1] *= ratio; | |
292 | } | |
293 | Float_t amp = sum3i[dk+3]*ratio; | |
294 | sum3wi+=amp; | |
295 | sum3iw+= dk*amp; | |
296 | } | |
297 | } | |
298 | meani = sum3iw/sum3wi; | |
299 | Float_t overlapi = (sum3wio-sum3wi)/sum3wio; | |
300 | ||
301 | ||
302 | ||
303 | //unfold z | |
304 | Float_t sum3wj = 0; //charge minus overlap | |
305 | Float_t sum3wjo = 0; //full charge | |
306 | Float_t sum3jw = 0; //sum for mean value | |
307 | for (Int_t dk=-1;dk<=1;dk++){ | |
308 | sum3wjo+=sum3j[dk+3]; | |
309 | if (dk==0){ | |
310 | sum3wj+=sum3j[dk+3]; | |
311 | } | |
312 | else{ | |
313 | Float_t ratio =1; | |
314 | if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) || | |
315 | (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){ | |
316 | Float_t xm2 = sum3j[-dk+3]; | |
317 | Float_t xm1 = sum3j[+3]; | |
318 | Float_t x1 = sum3j[2*dk+3]; | |
319 | Float_t x2 = sum3j[3*dk+3]; | |
cc5e9db0 | 320 | Float_t w11 = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001); |
321 | Float_t w12 = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.); | |
1c53abe2 | 322 | ratio = w11/(w11+w12); |
323 | for (Int_t dl=-1;dl<=1;dl++) | |
324 | mratio[dl+1][dk+1] *= ratio; | |
325 | } | |
326 | Float_t amp = sum3j[dk+3]*ratio; | |
327 | sum3wj+=amp; | |
328 | sum3jw+= dk*amp; | |
329 | } | |
330 | } | |
331 | meanj = sum3jw/sum3wj; | |
332 | Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo; | |
333 | overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3); | |
334 | sumu = (sum3wj+sum3wi)/2.; | |
335 | ||
336 | if (overlap ==3) { | |
337 | //if not overlap detected remove everything | |
338 | for (Int_t di =-2; di<=2;di++) | |
339 | for (Int_t dj =-2; dj<=2;dj++){ | |
340 | recmatrix[di+2][dj+2] = matrix2[3+di][dj]; | |
341 | } | |
342 | } | |
343 | else{ | |
344 | for (Int_t di =-1; di<=1;di++) | |
345 | for (Int_t dj =-1; dj<=1;dj++){ | |
346 | Float_t ratio =1; | |
347 | if (mratio[di+1][dj+1]==1){ | |
348 | recmatrix[di+2][dj+2] = matrix2[3+di][dj]; | |
349 | if (TMath::Abs(di)+TMath::Abs(dj)>1){ | |
350 | recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj]; | |
351 | recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj]; | |
352 | } | |
353 | recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj]; | |
354 | } | |
355 | else | |
356 | { | |
357 | //if we have overlap in direction | |
358 | recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj]; | |
359 | if (TMath::Abs(di)+TMath::Abs(dj)>1){ | |
cc5e9db0 | 360 | ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.); |
1c53abe2 | 361 | recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2]; |
362 | // | |
cc5e9db0 | 363 | ratio = TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.); |
1c53abe2 | 364 | recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2]; |
365 | } | |
366 | else{ | |
367 | ratio = recmatrix[di+2][dj+2]/matrix2[3][0]; | |
368 | recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2]; | |
369 | } | |
370 | } | |
371 | } | |
372 | } | |
373 | if (gDebug>4) | |
374 | printf("%f\n", recmatrix[2][2]); | |
375 | ||
376 | } | |
377 | ||
378 | Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz) | |
379 | { | |
380 | // | |
381 | // estimate max | |
382 | Float_t sumteor= 0; | |
383 | Float_t sumamp = 0; | |
384 | ||
385 | for (Int_t di = -1;di<=1;di++) | |
386 | for (Int_t dj = -1;dj<=1;dj++){ | |
387 | if (vmatrix[2+di][2+dj]>2){ | |
388 | Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2); | |
389 | sumteor += teor*vmatrix[2+di][2+dj]; | |
390 | sumamp += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj]; | |
391 | } | |
392 | } | |
393 | Float_t max = sumamp/sumteor; | |
394 | return max; | |
395 | } | |
396 | ||
397 | void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){ | |
398 | // | |
399 | // transform cluster to the global coordinata | |
400 | // add the cluster to the array | |
401 | // | |
402 | Float_t meani = c.GetY()/fPadWidth; | |
403 | Float_t meanj = c.GetZ()/fZWidth; | |
404 | ||
405 | Int_t ki = TMath::Nint(meani-3); | |
406 | if (ki<0) ki=0; | |
407 | if (ki>=fMaxPad) ki = fMaxPad-1; | |
408 | Int_t kj = TMath::Nint(meanj-3); | |
409 | if (kj<0) kj=0; | |
410 | if (kj>=fMaxTime-3) kj=fMaxTime-4; | |
411 | // ki and kj shifted to "real" coordinata | |
f8aae377 | 412 | if (fRowDig) { |
413 | c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0); | |
414 | c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1); | |
415 | c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2); | |
416 | } | |
1c53abe2 | 417 | |
418 | ||
419 | Float_t s2 = c.GetSigmaY2(); | |
420 | Float_t w=fParam->GetPadPitchWidth(fSector); | |
421 | ||
422 | c.SetSigmaY2(s2*w*w); | |
423 | s2 = c.GetSigmaZ2(); | |
424 | w=fZWidth; | |
425 | c.SetSigmaZ2(s2*w*w); | |
426 | c.SetY((meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector)); | |
427 | c.SetZ(fZWidth*(meanj-3)); | |
428 | c.SetZ(c.GetZ() - 3.*fParam->GetZSigma()); // PASA delay | |
429 | c.SetZ(fSign*(fParam->GetZLength() - c.GetZ())); | |
430 | ||
431 | if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) { | |
432 | //c.SetSigmaY2(c.GetSigmaY2()*25.); | |
433 | //c.SetSigmaZ2(c.GetSigmaZ2()*4.); | |
434 | c.SetType(-(c.GetType()+3)); //edge clusters | |
435 | } | |
436 | if (fLoop==2) c.SetType(100); | |
437 | ||
438 | TClonesArray * arr = fRowCl->GetArray(); | |
45bcf167 | 439 | // AliTPCclusterMI * cl = |
440 | new ((*arr)[fNcluster]) AliTPCclusterMI(c); | |
1c53abe2 | 441 | |
442 | fNcluster++; | |
443 | } | |
444 | ||
445 | ||
446 | //_____________________________________________________________________________ | |
f8aae377 | 447 | void AliTPCclustererMI::Digits2Clusters() |
1c53abe2 | 448 | { |
449 | //----------------------------------------------------------------- | |
450 | // This is a simple cluster finder. | |
451 | //----------------------------------------------------------------- | |
1c53abe2 | 452 | |
f8aae377 | 453 | if (!fInput) { |
454 | Error("Digits2Clusters", "input tree not initialised"); | |
1c53abe2 | 455 | return; |
456 | } | |
457 | ||
f8aae377 | 458 | if (!fOutput) { |
459 | Error("Digits2Clusters", "output tree not initialised"); | |
460 | return; | |
1c53abe2 | 461 | } |
462 | ||
463 | AliSimDigits digarr, *dummy=&digarr; | |
464 | fRowDig = dummy; | |
465 | fInput->GetBranch("Segment")->SetAddress(&dummy); | |
466 | Stat_t nentries = fInput->GetEntries(); | |
467 | ||
f8aae377 | 468 | fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3 after |
1c53abe2 | 469 | |
1c53abe2 | 470 | Int_t nclusters = 0; |
471 | ||
472 | for (Int_t n=0; n<nentries; n++) { | |
473 | fInput->GetEvent(n); | |
474 | Int_t row; | |
f8aae377 | 475 | if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,row)) { |
1c53abe2 | 476 | cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl; |
477 | continue; | |
478 | } | |
479 | ||
480 | AliTPCClustersRow *clrow= new AliTPCClustersRow(); | |
481 | fRowCl = clrow; | |
482 | clrow->SetClass("AliTPCclusterMI"); | |
483 | clrow->SetArray(1); | |
484 | ||
485 | clrow->SetID(digarr.GetID()); | |
486 | fOutput->GetBranch("Segment")->SetAddress(&clrow); | |
f8aae377 | 487 | fRx=fParam->GetPadRowRadii(fSector,row); |
1c53abe2 | 488 | |
489 | ||
f8aae377 | 490 | const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector(); |
1c53abe2 | 491 | fZWidth = fParam->GetZWidth(); |
492 | if (fSector < kNIS) { | |
f8aae377 | 493 | fMaxPad = fParam->GetNPadsLow(row); |
1c53abe2 | 494 | fSign = (fSector < kNIS/2) ? 1 : -1; |
f8aae377 | 495 | fPadLength = fParam->GetPadPitchLength(fSector,row); |
496 | fPadWidth = fParam->GetPadPitchWidth(); | |
1c53abe2 | 497 | } else { |
f8aae377 | 498 | fMaxPad = fParam->GetNPadsUp(row); |
1c53abe2 | 499 | fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; |
f8aae377 | 500 | fPadLength = fParam->GetPadPitchLength(fSector,row); |
501 | fPadWidth = fParam->GetPadPitchWidth(); | |
1c53abe2 | 502 | } |
503 | ||
504 | ||
505 | fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after | |
506 | fBins =new Int_t[fMaxBin]; | |
507 | fResBins =new Int_t[fMaxBin]; //fBins with residuals after 1 finder loop | |
508 | memset(fBins,0,sizeof(Int_t)*fMaxBin); | |
509 | ||
510 | if (digarr.First()) //MI change | |
511 | do { | |
512 | Short_t dig=digarr.CurrentDigit(); | |
f8aae377 | 513 | if (dig<=fParam->GetZeroSup()) continue; |
1c53abe2 | 514 | Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3; |
515 | fBins[i*fMaxTime+j]=dig; | |
516 | } while (digarr.Next()); | |
517 | digarr.ExpandTrackBuffer(); | |
518 | ||
f8aae377 | 519 | FindClusters(); |
8569a2b0 | 520 | |
521 | fOutput->Fill(); | |
88cb7938 | 522 | delete clrow; |
523 | nclusters+=fNcluster; | |
8569a2b0 | 524 | delete[] fBins; |
525 | delete[] fResBins; | |
88cb7938 | 526 | } |
f8aae377 | 527 | |
528 | Info("Digits2Clusters", "Number of found clusters : %d\n", nclusters); | |
529 | } | |
530 | ||
531 | void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) | |
532 | { | |
533 | //----------------------------------------------------------------- | |
534 | // This is a cluster finder for raw data. | |
535 | //----------------------------------------------------------------- | |
536 | ||
537 | if (!fOutput) { | |
538 | Error("Digits2Clusters", "output tree not initialised"); | |
539 | return; | |
540 | } | |
541 | ||
542 | rawReader->Reset(); | |
543 | AliTPCRawStream input(rawReader); | |
544 | ||
545 | fRowDig = NULL; | |
546 | ||
547 | Int_t nclusters = 0; | |
88cb7938 | 548 | |
f8aae377 | 549 | fMaxTime = fParam->GetMaxTBin() + 6; // add 3 virtual time bins before and 3 after |
550 | const Int_t kNIS = fParam->GetNInnerSector(); | |
551 | const Int_t kNOS = fParam->GetNOuterSector(); | |
552 | const Int_t kNS = kNIS + kNOS; | |
553 | fZWidth = fParam->GetZWidth(); | |
554 | Int_t zeroSup = fParam->GetZeroSup(); | |
555 | ||
556 | fBins = NULL; | |
557 | Int_t** splitRows = new Int_t* [kNS*2]; | |
558 | Int_t** splitRowsRes = new Int_t* [kNS*2]; | |
559 | for (Int_t iSector = 0; iSector < kNS*2; iSector++) | |
560 | splitRows[iSector] = NULL; | |
561 | Int_t iSplitRow = -1; | |
562 | ||
563 | Bool_t next = kTRUE; | |
564 | while (next) { | |
565 | next = input.Next(); | |
566 | ||
567 | // when the sector or row number has changed ... | |
568 | if (input.IsNewRow() || !next) { | |
569 | ||
570 | // ... find clusters in the previous pad row, and ... | |
571 | if (fBins) { | |
572 | if ((iSplitRow < 0) || splitRows[fSector + kNS*iSplitRow]) { | |
573 | fRowCl = new AliTPCClustersRow; | |
574 | fRowCl->SetClass("AliTPCclusterMI"); | |
575 | fRowCl->SetArray(1); | |
576 | fRowCl->SetID(fParam->GetIndex(fSector, input.GetPrevRow())); | |
577 | fOutput->GetBranch("Segment")->SetAddress(&fRowCl); | |
578 | ||
579 | FindClusters(); | |
580 | ||
581 | fOutput->Fill(); | |
582 | delete fRowCl; | |
583 | nclusters += fNcluster; | |
584 | delete[] fBins; | |
585 | delete[] fResBins; | |
586 | if (iSplitRow >= 0) splitRows[fSector + kNS*iSplitRow] = NULL; | |
587 | ||
588 | } else if (iSplitRow >= 0) { | |
589 | splitRows[fSector + kNS*iSplitRow] = fBins; | |
590 | splitRowsRes[fSector + kNS*iSplitRow] = fResBins; | |
591 | } | |
592 | } | |
593 | ||
594 | if (!next) break; | |
595 | ||
596 | // ... prepare for the next pad row | |
597 | fSector = input.GetSector(); | |
598 | Int_t iRow = input.GetRow(); | |
599 | fRx = fParam->GetPadRowRadii(fSector, iRow); | |
600 | ||
601 | iSplitRow = -1; | |
602 | if (fSector < kNIS) { | |
603 | fMaxPad = fParam->GetNPadsLow(iRow); | |
604 | fSign = (fSector < kNIS/2) ? 1 : -1; | |
605 | if (iRow == 30) iSplitRow = 0; | |
606 | } else { | |
607 | fMaxPad = fParam->GetNPadsUp(iRow); | |
608 | fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; | |
609 | if (iRow == 27) iSplitRow = 0; | |
610 | else if (iRow == 76) iSplitRow = 1; | |
611 | } | |
612 | fPadLength = fParam->GetPadPitchLength(fSector, iRow); | |
613 | fPadWidth = fParam->GetPadPitchWidth(); | |
614 | ||
615 | fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after | |
616 | if ((iSplitRow < 0) || !splitRows[fSector + kNS*iSplitRow]) { | |
617 | fBins = new Int_t[fMaxBin]; | |
618 | fResBins = new Int_t[fMaxBin]; //fBins with residuals after 1 finder loop | |
619 | memset(fBins, 0, sizeof(Int_t)*fMaxBin); | |
620 | } else { | |
621 | fBins = splitRows[fSector + kNS*iSplitRow]; | |
622 | fResBins = splitRowsRes[fSector + kNS*iSplitRow]; | |
623 | } | |
624 | } | |
625 | ||
626 | // fill fBins with digits data | |
627 | if (input.GetSignal() <= zeroSup) continue; | |
628 | Int_t i = input.GetPad() + 3; | |
629 | Int_t j = input.GetTime() + 3; | |
630 | fBins[i*fMaxTime+j] = input.GetSignal(); | |
631 | } | |
632 | ||
633 | // find clusters in split rows that were skipped until now. | |
634 | // this can happen if the rows were not splitted | |
635 | for (fSector = 0; fSector < kNS; fSector++) | |
636 | for (Int_t iSplit = 0; iSplit < 2; iSplit++) | |
637 | if (splitRows[fSector + kNS*iSplit]) { | |
638 | ||
639 | Int_t iRow = -1; | |
640 | if (fSector < kNIS) { | |
641 | iRow = 30; | |
642 | fMaxPad = fParam->GetNPadsLow(iRow); | |
643 | fSign = (fSector < kNIS/2) ? 1 : -1; | |
644 | } else { | |
645 | if (iSplit == 0) iRow = 27; else iRow = 76; | |
646 | fMaxPad = fParam->GetNPadsUp(iRow); | |
647 | fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; | |
648 | } | |
649 | fRx = fParam->GetPadRowRadii(fSector, iRow); | |
650 | fPadLength = fParam->GetPadPitchLength(fSector, iRow); | |
651 | fPadWidth = fParam->GetPadPitchWidth(); | |
652 | ||
653 | fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after | |
654 | fBins = splitRows[fSector + kNS*iSplit]; | |
655 | fResBins = splitRowsRes[fSector + kNS*iSplit]; | |
656 | ||
657 | fRowCl = new AliTPCClustersRow; | |
658 | fRowCl->SetClass("AliTPCclusterMI"); | |
659 | fRowCl->SetArray(1); | |
660 | fRowCl->SetID(fParam->GetIndex(fSector, iRow)); | |
661 | fOutput->GetBranch("Segment")->SetAddress(&fRowCl); | |
662 | ||
663 | FindClusters(); | |
664 | ||
665 | fOutput->Fill(); | |
666 | delete fRowCl; | |
667 | nclusters += fNcluster; | |
668 | delete[] fBins; | |
669 | delete[] fResBins; | |
670 | } | |
671 | ||
672 | delete[] splitRows; | |
673 | delete[] splitRowsRes; | |
674 | Info("Digits2Clusters", "Number of found clusters : %d\n", nclusters); | |
675 | } | |
676 | ||
677 | void AliTPCclustererMI::FindClusters() | |
678 | { | |
679 | //add virtual charge at the edge | |
680 | for (Int_t i=0; i<fMaxTime; i++){ | |
681 | Float_t amp1 = fBins[i+3*fMaxTime]; | |
682 | Float_t amp0 =0; | |
683 | if (amp1>0){ | |
684 | Float_t amp2 = fBins[i+4*fMaxTime]; | |
685 | if (amp2==0) amp2=0.5; | |
686 | Float_t sigma2 = GetSigmaY2(i); | |
687 | amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2); | |
688 | if (gDebug>4) printf("\n%f\n",amp0); | |
689 | } | |
690 | fBins[i+2*fMaxTime] = Int_t(amp0); | |
691 | amp0 = 0; | |
692 | amp1 = fBins[(fMaxPad+2)*fMaxTime+i]; | |
693 | if (amp1>0){ | |
694 | Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime]; | |
695 | if (amp2==0) amp2=0.5; | |
696 | Float_t sigma2 = GetSigmaY2(i); | |
697 | amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2); | |
698 | if (gDebug>4) printf("\n%f\n",amp0); | |
699 | } | |
700 | fBins[(fMaxPad+3)*fMaxTime+i] = Int_t(amp0); | |
701 | } | |
702 | ||
703 | // memcpy(fResBins,fBins, fMaxBin*2); | |
704 | memcpy(fResBins,fBins, fMaxBin); | |
705 | // | |
706 | fNcluster=0; | |
707 | //first loop - for "gold cluster" | |
708 | fLoop=1; | |
709 | Int_t *b=&fBins[-1]+2*fMaxTime; | |
710 | Int_t crtime = Int_t((fParam->GetZLength()-1.05*fRx)/fZWidth-5); | |
711 | ||
712 | for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) { | |
713 | b++; | |
714 | if (*b<8) continue; //threshold form maxima | |
715 | if (i%fMaxTime<crtime) { | |
716 | Int_t delta = -(i%fMaxTime)+crtime; | |
717 | b+=delta; | |
718 | i+=delta; | |
719 | continue; | |
720 | } | |
721 | ||
722 | if (!IsMaximum(*b,fMaxTime,b)) continue; | |
723 | AliTPCclusterMI c; | |
724 | Int_t dummy=0; | |
725 | MakeCluster(i, fMaxTime, fBins, dummy,c); | |
726 | //} | |
727 | } | |
728 | //memcpy(fBins,fResBins, fMaxBin*2); | |
729 | //second loop - for rest cluster | |
730 | /* | |
731 | fLoop=2; | |
732 | b=&fResBins[-1]+2*fMaxTime; | |
733 | for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) { | |
734 | b++; | |
735 | if (*b<25) continue; // bigger threshold for maxima | |
736 | if (!IsMaximum(*b,fMaxTime,b)) continue; | |
737 | AliTPCclusterMI c; | |
738 | Int_t dummy; | |
739 | MakeCluster(i, fMaxTime, fResBins, dummy,c); | |
740 | //} | |
741 | } | |
742 | */ | |
8569a2b0 | 743 | } |