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