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