]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCClusterFinder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes suggested by Alessandra and Paolo to avoid overlapped
[u/mrichter/AliRoot.git] / TPC / AliTPCClusterFinder.cxx
CommitLineData
cc80f89e 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$Log$
67d37b32 18Revision 1.3 2000/07/10 20:57:39 hristov
19Update of TPC code and macros by M.Kowalski
20
37831078 21Revision 1.2 2000/06/30 12:07:49 kowal2
22Updated from the TPC-PreRelease branch
23
73042f01 24Revision 1.1.2.1 2000/06/25 08:52:51 kowal2
25replacing AliClusterFinder
cc80f89e 26
27*/
28
29//-----------------------------------------------------------------------------
30//
67d37b32 31// Implementation of class ALITPCCLUSTERFINDER
cc80f89e 32//
33//Class for cluster finding in two dimension.
34//In the present there are implemented two algorithm
35//primitive recursion algorithm. (FindPeaks)
36//Algorithm is not working in case of overlaping clusters
37//Maximum - minimum in direction algoritm (Find clusters)
38//In this algoritm we suppose that each cluster has local
39//maximum. From this local maximum I mus see each point
40//of cluster.
41//From maximum i can accept every point in radial
42//direction which is before border in direction
43//Border in direction occur if we have next in
44//direction nder threshold or response begin
45//to increase in given radial direction
46//-----------------------------------------------------------------------------
73042f01 47
cc80f89e 48#include "TMinuit.h"
49#include "AliArrayI.h"
50#include "TClonesArray.h"
51#include "AliTPC.h"
52#include "TRandom.h"
53#include "AliH2F.h"
54#include "TMarker.h"
67d37b32 55#include "AliComplexCluster.h"
73042f01 56#include "AliTPCClusterFinder.h"
57#include <fstream.h>
cc80f89e 58
59//direction constants possible direction in 8 different sectors
60//
61
62
63const Int_t kClStackSize =1000;
64
65
66
67
73042f01 68static AliTPCClusterFinder * gClusterFinder; //for fitting routine
cc80f89e 69
70void gauss(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
71{
72 AliArrayI * points = gClusterFinder->GetStack();
73 const Int_t nbins = gClusterFinder->GetStackIndex();
74 Int_t i;
75 //calculate chisquare
76 Double_t chisq = 0;
77 Double_t delta;
78 for (i=0;i<nbins; i++) {
79 Float_t x = points->At(i*3);
80 Float_t y = points->At(i*3+1);
81 Float_t z = points->At(i*3+2);
82 Float_t deltax2 = (x-par[1]);
83 deltax2*=deltax2;
84 deltax2*=par[3];
85 Float_t deltay2 = (y-par[2]);
86 deltay2*=deltay2;
87 deltay2*=par[4];
88
89 delta = z-par[0]*TMath::Exp(-deltax2-deltay2);
90 chisq += delta*delta;
91 }
92 f = chisq;
93}
94
95
73042f01 96ClassImp(AliTPCClusterFinder)
97 //ClassImp(AliCell)
cc80f89e 98
73042f01 99AliTPCClusterFinder::AliTPCClusterFinder()
cc80f89e 100{
101 fDigits =0;
102 fDimX = 0;
103 fDimY = 0;
cc80f89e 104 fNoiseTh = 3;
105 fMulSigma2 = 16; //4 sigma
106 fDirSigmaFac = 1.4;
107 fDirAmpFac =1.3;
108 fNType=8;
109 fThreshold = 2;
110 fStack = new AliArrayI;
111 fStack->Set(kClStackSize);
112 fClustersArray =0;
73042f01 113 SetSigmaX(1,0,0);
114 SetSigmaY(1,0,0);
115
cc80f89e 116
117 fDetectorParam = 0;
118 ResetStatus();
119 fBFit = kFALSE;
120 fMinuit= new TMinuit(5);
121 fMinuit->SetFCN(gauss);
122 gClusterFinder = this;
123
124}
125
126
73042f01 127AliTPCClusterFinder::~AliTPCClusterFinder()
cc80f89e 128{
129 if (fDigits != 0) delete fDigits;
130}
131
73042f01 132void AliTPCClusterFinder::SetSigmaX(Float_t s0, Float_t s1x, Float_t s1y)
133{
134 fSigmaX[0]=s0;
135 fSigmaX[1]=s1x;
136 fSigmaX[2]=s1y;
137
138}
139void AliTPCClusterFinder::SetSigmaY(Float_t s0, Float_t s1x, Float_t s1y)
140{
141 fSigmaY[0]=s0;
142 fSigmaY[1]=s1x;
143 fSigmaY[2]=s1y;
144}
145
146
147
148Bool_t AliTPCClusterFinder::SetSigma2(Int_t i, Int_t j, Float_t & sigmax2, Float_t &sigmay2)
149{
150 //
151 //set sigmax2 and sigma y2 accordig i and j position of cell
152 //
153
154 // Float_t x[3] = {ItoX(i),JtoY(j),0};
155 Float_t x= ItoX(i);
156 Float_t y= JtoY(j);
157
158 sigmax2= fSigmaX[0]+fSigmaX[1]*x+fSigmaX[2]*y;
159 sigmay2= fSigmaY[0]+fSigmaY[1]*x+fSigmaY[2]*y;
160 return kTRUE;
161}
162
163/*
164Bool_t AliTPCClusterFinder::SetSigma2(Int_t i, Int_t j, Float_t & sigmax2, Float_t &sigmay2)
165{
166 //
167 //set sigmax2 and sigma y2 accordig i and j position of cell
168 //
169 if (fDetectorParam==0) {
170 sigmax2=1;
171 sigmay2=1;
172 return kFALSE;
173 }
174 Float_t x[3] = {ItoX(i),JtoY(j),0};
175 Float_t sigma[2];
176 fDetectorParam->GetClusterSize(x,fDetectorIndex,0,0,sigma);
177 sigmax2=sigma[0]*(fX2-fX1)*(fX2-fX1)/(fDimX*fDimX);
178 sigmay2=sigma[1]*(fY2-fY1)*(fY2-fY1)/(fDimY*fDimY);
179 return kTRUE;
180}
181*/
182
183
184void AliTPCClusterFinder::GetHisto(TH2F * his2)
cc80f89e 185{
186
73042f01 187 UInt_t idim =his2->GetNbinsX();
188 UInt_t jdim =his2->GetNbinsY();
cc80f89e 189 fX1 = his2->GetXaxis()->GetXmin();
190 fX2 = his2->GetXaxis()->GetXmax();
191 fY1 = his2->GetYaxis()->GetXmin();
192 fY2 = his2->GetYaxis()->GetXmax();
193
194 if ( (idim>0) && (jdim>0))
195 {
196 rOK = kTRUE;
197 fDimX = idim;
198 fDimY = jdim;
73042f01 199 Int_t size =idim*jdim;
cc80f89e 200 if (fDigits !=0) delete fDigits;
73042f01 201 fDigits = (Int_t*) new Int_t[size];
202 fCells = (AliCell*) new AliCell[size];
203
cc80f89e 204 } else
205 rOK=kFALSE;
37831078 206 for (Int_t i = 0; i<(Int_t)idim;i++)
207 for (Int_t j = 0; j<(Int_t)jdim;j++)
cc80f89e 208 {
209 Int_t index = his2->GetBin(i+1,j+1);
73042f01 210 //AliCell * cell = GetCell(i,j);
211 //if (cell!=0) cell->SetSignal(his2->GetBinContent(index));
212 SetSignal(his2->GetBinContent(index),i,j);
cc80f89e 213 }
214
215}
216
217
218
219
73042f01 220void AliTPCClusterFinder::FindMaxima()
cc80f89e 221{
222 for (Int_t i=0; i<fDimX; i++)
223 for (Int_t j=0;j<fDimY; j++)
224 if (IsMaximum(i,j)) cout<<i<<" "<<j<<"\n";
225}
226
227
73042f01 228void AliTPCClusterFinder::Transform(AliDigitCluster * c)
cc80f89e 229{
230 //transform coordinata from bin coordinata to "normal coordinata"
231 //for example if we initialize finder with histogram
232 //it transform values from bin coordinata to the histogram coordinata
233 c->fX=ItoX(c->fX);
234 c->fY=JtoY(c->fY);
235 c->fMaxX=ItoX(c->fMaxX);
236 c->fMaxY=JtoY(c->fMaxY);
237
238 c->fSigmaX2=c->fSigmaX2*(fX2-fX1)*(fX2-fX1)/(fDimX*fDimX);
239 c->fSigmaY2=c->fSigmaY2*(fY2-fY1)*(fY2-fY1)/(fDimY*fDimY);
240 c->fArea =c->fArea*(fX2-fX1)*(fY2-fY1)/(fDimX*fDimY);
241}
73042f01 242void AliTPCClusterFinder::AddToStack(Int_t i, Int_t j, Int_t signal)
cc80f89e 243{
244 //
245 //add digit to stack
246 //
247 if ( ((fStackIndex+2)>=kClStackSize) || (fStackIndex<0) ) return;
248 fStack->AddAt(i,fStackIndex);
249 fStack->AddAt(j,fStackIndex+1);
250 fStack->AddAt(signal,fStackIndex+2);
251 fStackIndex+=3;
252}
253
73042f01 254void AliTPCClusterFinder::GetClusterStatistic(AliDigitCluster & cluster)
cc80f89e 255{
256 //
257 //calculate statistic of cluster
258 //
259 Double_t sumxw,sumyw,sumx2w,sumy2w,sumxyw,sumw;
73042f01 260 Int_t minx,maxx,miny,maxy;
cc80f89e 261 sumxw=sumyw=sumx2w=sumy2w=sumxyw=sumw=0;
262 minx=fDimX;
263 maxx=-fDimX;
264 miny=fDimY;
265 maxy=-fDimY;
266 Int_t x0=fStack->At(0);
267 Int_t y0=fStack->At(1);
73042f01 268 Int_t maxQx =x0;
269 Int_t maxQy =y0;
270 Int_t maxQ=fStack->At(2);
271
272
cc80f89e 273 for (Int_t i = 0; i<fStackIndex;i+=3){
274 Int_t x = fStack->At(i);
275 Int_t y = fStack->At(i+1);
276 Int_t dx=x-x0;
277 Int_t dy=y-y0;
278 Int_t w = fStack->At(i+2);
73042f01 279 if (w>maxQ){
280 maxQ = w;
281 maxQx = x;
282 maxQy=y;
283 }
cc80f89e 284 if (x<minx) minx=x;
285 if (y<miny) miny=y;
286 if (x>maxx) maxx=x;
287 if (y>maxy) maxy=y;
288 sumxw+=dx*w;
289 sumyw+=dy*w;
290 sumx2w+=dx*dx*w;
291 sumy2w+=dy*dy*w;
292 sumxyw+=dx*dy*w;
293 sumw+=w;
294 }
295 cluster.fQ = sumw;
296 if (sumw>0){
297 cluster.fX = sumxw/sumw;
298 cluster.fY = sumyw/sumw;
299 cluster.fQ = sumw;
300 cluster.fSigmaX2 = sumx2w/sumw-cluster.fX*cluster.fX;
301 cluster.fSigmaY2 = sumy2w/sumw-cluster.fY*cluster.fY;
302 cluster.fSigmaXY = sumxyw/sumw-cluster.fX*cluster.fY;
73042f01 303 cluster.fMaxX = maxQx;
304 cluster.fMaxY = maxQy;
cc80f89e 305 cluster.fMax = maxQ;
306 cluster.fArea = fStackIndex/3;
307 cluster.fNx = maxx-minx+1;
308 cluster.fNy = maxy-miny+1;
309 cluster.fX +=x0;
310 cluster.fY +=y0;
311 }
312}
73042f01 313void AliTPCClusterFinder::GetClusterFit(AliDigitCluster & cluster)
cc80f89e 314{
315 //
316 //calculate statistic of cluster
317 //
318 Double_t arglist[10];
319 Int_t ierflg = 0;
320
321 arglist[0] = 1;
322 fMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
323
324 //fistly find starting parameters
325 Int_t minx,maxx,miny,maxy,maxQ,maxQx,maxQy;
326 Int_t over =0;
327 Float_t sumxw,sumyw,sumw;
328 sumxw=sumyw=sumw=0;
329 minx=fDimX;
330 maxx=-fDimX;
331 miny=fDimY;
332 maxy=-fDimY;
333 maxQx=fStack->At(0);
334 maxQy=fStack->At(1);
335 maxQ=fStack->At(2);
336
337 for (Int_t i = 0; i<fStackIndex;i+=3){
338 Int_t x = fStack->At(i);
339 Int_t y = fStack->At(i+1);
340 Int_t w = fStack->At(i+2);
341 if (w>fThreshold) {
342 over++;
343 sumw+=w;
344 sumxw+=x*w;
345 sumyw+=y*w;
346 if (x<minx) minx=x;
347 if (y<miny) miny=y;
348 if (x>maxx) maxx=x;
349 if (y>maxy) maxy=y;
350 if (w>maxQ) {
351 maxQ=w;
352 maxQx=x;
353 maxQy=y;
354 }
355 }
356 }
357 Int_t nx = maxx-minx+1;
358 Int_t ny = maxy-miny+1;
359
360 SetSigma2(maxQx,maxQy,fCurrentSigmaX2,fCurrentSigmaY2);
361 Double_t vstart[5]={maxQ,sumxw/sumw,sumyw/sumw,1/(2*fCurrentSigmaX2),1/(2*fCurrentSigmaY2)};
362 Double_t step[5]={1.,0.01,0.01,0.01,0.01};
363 fMinuit->mnparm(0, "amp", vstart[0], step[0], 0,0,ierflg);
364 fMinuit->mnparm(1, "x0", vstart[1], step[1], 0,0,ierflg);
365 fMinuit->mnparm(2, "y0", vstart[2], step[2], 0,0,ierflg);
366 fMinuit->mnparm(3, "sx2", vstart[3], step[3], 0,0,ierflg);
367 fMinuit->mnparm(4, "sy2", vstart[4], step[4], 0,0,ierflg);
368 arglist[0] = 500;
369 arglist[1] = 1.;
370
371 fMinuit->mnfree(0); //set unfixed all parameters
372 //if we have area less then
373 if (over<=21) { //if we dont't have more then 7 points
374 fMinuit->FixParameter(3);
375 fMinuit->FixParameter(4);
376 }
377 else {
378 if (nx<3) fMinuit->FixParameter(3); //fix sigma x if no data in x direction
379 if (ny<3) fMinuit->FixParameter(4); //fix sigma y if no data in y direction
380 }
381 fMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
382
383 if (sumw>0){
384 Double_t x[5];
385 Double_t error[5];
386 fMinuit->GetParameter(0,x[0],error[0]);
387 fMinuit->GetParameter(1,x[1],error[1]);
388 fMinuit->GetParameter(2,x[2],error[2]);
389 fMinuit->GetParameter(3,x[3],error[3]);
390 fMinuit->GetParameter(4,x[4],error[4]);
391
392 cluster.fX = x[1];
393 cluster.fY = x[2];
394 cluster.fMaxX = maxQx;
395 cluster.fMaxY = maxQy;
396
397 cluster.fQ = sumw;
398 cluster.fSigmaX2 = 1/TMath::Sqrt(2*x[3]);
399 cluster.fSigmaY2 = 1/TMath::Sqrt(2*x[4]);
400 cluster.fSigmaXY = 0;
401 cluster.fMax = x[0];
402 cluster.fArea = over;
403 cluster.fNx = nx;
404 cluster.fNy = ny;
405 }
406}
407
73042f01 408Bool_t AliTPCClusterFinder::CheckIfDirBorder(Float_t x, Float_t y,
cc80f89e 409 Int_t i,Int_t j)
410{
411 //
412 //function which control if given cell with index i, j is the
413 //minimum in direction
414 // x and y are estimate of local maximum
415 //direction is given by the
416 Float_t virtualcell;
417 AliCell * cellor= GetCell(i,j);
73042f01 418 Int_t sigor = GetSignal(i,j);
419
cc80f89e 420 //control derivation in direction
421 //if function grows up in direction then there is border
422 Float_t dx = i-x;
423 Float_t dy = j-y;
424 Float_t dd = TMath::Sqrt(dx*dx+dy*dy);
425 Float_t ddx = TMath::Abs(dx);
426 ddx = (ddx>0.5) ? ddx-0.5: 0;
427 ddx*=ddx;
428 Float_t ddy = TMath::Abs(dy);
429 ddy = (ddy>0.5) ? ddy-0.5: 0;
430 ddy*=ddy;
431 Float_t d2 = ddx/(2*fDirSigmaFac*fCurrentSigmaX2)+ddy/(2*fDirSigmaFac*fCurrentSigmaY2); //safety factor
432 //I accept sigmax and sigma y bigge by factor sqrt(fDirsigmaFac)
433 Float_t amp = TMath::Exp(-d2)*fCurrentMaxAmp*fDirAmpFac; //safety factor fDirFac>1
434
73042f01 435 if (sigor>amp) return kTRUE;
cc80f89e 436 if (dd==0) return kFALSE;
437
438 dx/=dd;
439 dy/=dd;
440 virtualcell = GetVirtualSignal(i+dx,j+dy);
441 if (virtualcell <=fThreshold) return kFALSE;
73042f01 442 if (virtualcell>sigor)
443 if (virtualcell>(sigor+fNoiseTh))
444 {cellor->SetDirBorder(fIndex); return kTRUE;}
cc80f89e 445 else
446 {
447 virtualcell = GetVirtualSignal(i+2*dx,j+2*dy);
73042f01 448 if (virtualcell>sigor)
449 { cellor->SetDirBorder(fIndex); return kTRUE;}
cc80f89e 450 };
451 return kFALSE;
452}
453
454
455
73042f01 456
457
458Bool_t AliTPCClusterFinder::IsMaximum(Int_t i, Int_t j)
cc80f89e 459{
460 //there is maximum if given digits is 1 sigma over all adjacent
461 //in 8 neighborow
462 //or ther exist virual maximum
463 //is maximum on 24 points neighboring
464 // Bool_t res = kFALSE;
465 Int_t over =0;
466 Int_t overth=0;
467 Int_t oversigma =0;
468 AliCell * cell = GetCell(i,j);
73042f01 469 Int_t signal = GetSignal(i,j);
cc80f89e 470 if (cell == 0) return kFALSE;
471 for ( Int_t di=-1;di<=1;di++)
472 for ( Int_t dj=-1;dj<=1;dj++){
473 if ( (di!=0) || (dj!=0))
474 {
475 AliCell * cell2=GetCell(i+di,j+dj);
73042f01 476 Int_t signal2 = GetSignal(i+di,j+dj);
cc80f89e 477 if (cell2 == 0) {
478 over+=1;
479 oversigma+=1;
480 }
481 else
482 {
73042f01 483 if (signal2>signal) return kFALSE;
484 if (signal2>fThreshold) overth++;
485 if (signal2==signal) {
cc80f89e 486 if (di<0) return kFALSE;
487 if ( (di+dj)<0) return kFALSE;
488 }
73042f01 489 // if (signal>=signal2){
cc80f89e 490 over+=1;
73042f01 491 if (signal>fNoiseTh+signal2)
cc80f89e 492 oversigma+=1;
493 //}
494 }
495 }
496 }
497 //if I have only one neighborough over threshold
498 if (overth<2) return kFALSE;
499 if (over<8) return kFALSE;
500 if (oversigma==8) {
501 fCurrentMaxX = i;
502 fCurrentMaxY = j;
73042f01 503 fCurrentMaxAmp =signal;
504 SetMaximum(fIndex,i,j);
cc80f89e 505 return kTRUE;
506 }
507 //check if there exist virtual maximum
73042f01 508 for (Float_t ddi=0.;(ddi<1.);ddi+=0.5)
509 for (Float_t ddj=0.;(ddj<1.);ddj+=0.5)
510 if (IsVirtualMaximum(Float_t(i)+ddi,Float_t(j)+ddj)){
511 fCurrentMaxX = i+ddi;
512 fCurrentMaxY = j+ddj;
513 fCurrentMaxAmp =signal;
514 SetMaximum(fIndex,i,j);
cc80f89e 515 return kTRUE;
516 }
517 return kFALSE;
518}
519
73042f01 520Bool_t AliTPCClusterFinder::IsVirtualMaximum(Float_t x, Float_t y)
cc80f89e 521{
522 //there is maximum if given digits is 1 sigma over all adjacent
523 //in 8 neighborow or
524 //is maximum on 24 points neighboring
525 Bool_t res = kFALSE;
526 Int_t over =0;
527 Int_t overth=0;
528 Int_t oversigma =0;
529 Float_t virtualcell = GetVirtualSignal(x,y);
530 if (virtualcell < 0) return kFALSE;
531 for ( Int_t di=-1;di<=1;di++)
532 for ( Int_t dj=-1;dj<=1;dj++)
533 if ( (di!=0) || (dj!=0))
534 {
535 Float_t virtualcell2=GetVirtualSignal(x+di,y+dj);
536 if (virtualcell2 < 0) {
537 over+=1;
538 oversigma+=1;
539 }
540 else
541 {
542 if (virtualcell2>fThreshold) overth++;
543 if (virtualcell>=virtualcell2){
544 over+=1;
545 if (virtualcell>fNoiseTh+virtualcell2)
546 oversigma+=1;
547 }
548 }
549 }
550 if (overth<2) return kFALSE;
551 //if there exist only one or less neighboring above threshold
552 if (oversigma==8) res = kTRUE;
553 else if ((over==8)&&(GetNType()==8)) res=kTRUE;
554 else if (over ==8 )
73042f01 555 for ( Int_t di=-2;di<=2;di++)
cc80f89e 556 for ( Int_t dj=-2;dj<=2;dj++)
73042f01 557 if ( (di==2)||(di==-2) || (dj==2)|| (dj==-2) )
cc80f89e 558 {
73042f01 559 Float_t virtualcell2=GetVirtualSignal(x+di,y+dj);
cc80f89e 560 if (virtualcell2 < 0) {
561 over+=1;
562 oversigma+=1;
563 }
564 else
565 {
566 if (virtualcell>=virtualcell2) over+=1;
567 }
568 }
569 if (over == 24) res=kTRUE;
570 return res;
571
572}
573
574
73042f01 575void AliTPCClusterFinder::ResetSignal()
cc80f89e 576{
577 //reset dignals to 0
578 Int_t size = fDimX*fDimY;
73042f01 579 AliCell *dig=fCells;
cc80f89e 580 if (rOK==kTRUE) for (Int_t i=0 ; i<size;i++) dig[i] = 0;
581}
582
583
584
73042f01 585void AliTPCClusterFinder::ResetStatus()
cc80f89e 586{
587 //reset status of signals to not used
588 Int_t size = fDimX*fDimY;
73042f01 589 AliCell *dig=fCells;
cc80f89e 590 if (rOK==kTRUE) for (Int_t i=0 ; i<size;i++)
591 dig[i].SetStatus(0);
592}
593
594
73042f01 595AliCell * AliTPCClusterFinder::GetCell(Int_t i, Int_t j)
cc80f89e 596{
597 //return reference to the cell with index i,j
598 if (rOK == kTRUE)
599 if ( (i>=0) && (i<fDimX) && (j>=0) && (j<fDimY) )
73042f01 600 return &fCells[i+j*fDimX];
cc80f89e 601 return 0;
602}
603
73042f01 604Float_t AliTPCClusterFinder::GetVirtualSignal(Float_t ri, Float_t rj)
cc80f89e 605{
606 //it generate virtual cell as mean value from different cels
607 //after using it must be destructed !!!
608 Int_t i=(Int_t)ri;
609 Int_t j=(Int_t)rj;
610 Int_t ddi = (ri>i)? 1:0;
611 Int_t ddj = (rj>j)? 1:0;
612 Float_t sum = 0;
613 Float_t sumw= 0;
614 for (Int_t di=0;di<=ddi;di++)
615 for (Int_t dj=0;dj<=ddj;dj++)
616 {
617 Float_t w = (ri-i-di)*(ri-i-di)+(rj-j-dj)*(rj-j-dj);
618 if (w>0) w=1/TMath::Sqrt(w);
619 else w=9999999;
620 AliCell * cel2 =GetCell(i+di,j+dj);
73042f01 621 Int_t signal2 = GetSignal(i+di,j+dj);
cc80f89e 622 if (cel2!=0) {
623 sumw+=w;
73042f01 624 sum+= signal2*w;
cc80f89e 625 }
626 }
627 if (sumw>0) return (sum/sumw);
628 else
629 return -1;
630}
631
632
633
73042f01 634void AliTPCClusterFinder::Streamer(TBuffer & R__b)
cc80f89e 635{
636 if (R__b.IsReading()) {
637 // Version_t R__v = R__b.ReadVersion();
638 } else {
73042f01 639 R__b.WriteVersion(AliTPCClusterFinder::IsA());
cc80f89e 640 }
641}
642
643
cc80f89e 644
73042f01 645void AliTPCClusterFinder::SetBlockIndex(Int_t * index)
cc80f89e 646{
647 //
648 //calculate which indexes we must check for border
649 //
650 if (TMath::Abs(index[0])<2) index[2] = 0;
651 else {
652 index[2] = TMath::Abs(index[0])-1;
653 if (index[0]<0) index[2]*=-1; //first x block
654 }
655 if (TMath::Abs(index[1])<2) index[3] = 0;
656 else {
657 index[3] = TMath::Abs(index[1])-1;
658 if (index[1]<0) index[3]*=-1; //first y block
659 }
660 if (TMath::Abs(index[0])<TMath::Abs(index[1])){
661 index[4]=index[0];
662 index[5]=index[3];
663 }
664 else
665 if (index[0]==index[1]) {
666 index[4]=0;
667 index[5]=0;
668 }
669 else{
670 index[4]=index[2];
671 index[5]=index[1];
672 }
673 return;
674}
675
676//***********************************************************************
677//***********************************************************************
678
73042f01 679TClonesArray * AliTPCClusterFinder::FindPeaks1(TClonesArray *arr)
cc80f89e 680{
681 //find peaks and write it in form of AliTPCcluster to array
73042f01 682 if (arr==0){
683 fClustersArray=new TClonesArray("AliDigitCluster",300);
684 fIndex=1;
685 }
686 else {
687 fClustersArray = arr;
688 fIndex = fClustersArray->GetEntriesFast();
689 }
690
691 AliDigitCluster c;
cc80f89e 692 ResetStatus();
cc80f89e 693 for (Int_t i=0; i<fDimX; i++)
694 for (Int_t j=0;j<fDimY; j++)
695 {
696 fStackIndex=0;
697 fBDistType = kFALSE;
698 AliCell * cell = GetCell(i,j);
699 if (!(cell->IsChecked())) Adjacent(i,j);
700 //if there exists more then 2 digits cluster
701 if (fStackIndex >2 ){
702 if (fBFit==kFALSE) GetClusterStatistic(c);
703 else GetClusterFit(c);
704 //write some important chracteristic area of cluster
705 //
706 Transform(&c);
707 //write cluster information to array
73042f01 708 TClonesArray &lclusters = *fClustersArray;
709 new (lclusters[fIndex++]) AliDigitCluster(c);
cc80f89e 710 // cout<<"fx="<<c.fX<<" fy"<<c.fY<<"\n";
711 }
712 }
73042f01 713 return fClustersArray;
cc80f89e 714}
715
716
73042f01 717TClonesArray * AliTPCClusterFinder::FindPeaks2(TClonesArray *arr)
cc80f89e 718{
719 //find peaks and write it in form of AliTPCcluster to array
73042f01 720 if (arr==0){
721 fClustersArray=new TClonesArray("AliDigitCluster",300);
722 fIndex=1;
723 }
724 else {
725 fClustersArray = arr;
726 fIndex = fClustersArray->GetEntriesFast();
727 }
728
729 AliDigitCluster c;
cc80f89e 730 ResetStatus();
731
732 for (Int_t i=0; i<fDimX; i++)
733 for (Int_t j=0;j<fDimY; j++)
734 {
735 fStackIndex=0;
736 if (IsMaximum(i,j) == kTRUE){
737 SetSigma2(i,j,fCurrentSigmaX2,fCurrentSigmaY2);
738 fBDistType = kTRUE;
739 Adjacent(i,j);
740 //if there exists more then 2 digits cluster
741 if (fStackIndex >2 ){
742 if (fBFit==kFALSE) GetClusterStatistic(c);
743 else GetClusterFit(c);
744 //write some important chracteristic area of cluster
745 //
746 Transform(&c);
747 //write cluster information to array
73042f01 748 TClonesArray &lclusters = *fClustersArray;
cc80f89e 749 new(lclusters[fIndex++]) AliDigitCluster(c);
750 // cout<<"fx="<<c.fX<<" fy"<<c.fY<<"\n";
751 }
752 }
753 }
73042f01 754 return fClustersArray;
cc80f89e 755}
756
757
73042f01 758TClonesArray * AliTPCClusterFinder::FindPeaks3(TClonesArray *arr)
cc80f89e 759{
760 //find peaks and write it in form of AliTPCcluster to array
73042f01 761 if (arr==0){
762 fClustersArray=new TClonesArray("AliDigitCluster",300);
763 fIndex=1;
764 }
765 else {
766 fClustersArray = arr;
767 fIndex = fClustersArray->GetEntriesFast();
768 }
769
770 AliDigitCluster c;
cc80f89e 771 ResetStatus();
772
773 Int_t dmax=5;
774 Int_t naccepted =1;
775 for (Int_t i=0; i<fDimX; i++)
776 for (Int_t j=0;j<fDimY; j++)
777 {
778 fStackIndex=0;
779 if (IsMaximum(i,j) == kTRUE){
780 SetSigma2(i,j,fCurrentSigmaX2,fCurrentSigmaY2);
73042f01 781 AddToStack(i,j,GetSignal(i,j));
cc80f89e 782
783 //loop over different distance
784 naccepted =1;
785 for ( Int_t dd =1;((dd<=dmax) && (naccepted>0));dd++){
786 naccepted=0;
787 for (Int_t di = -dd;di<=dd;di++){
67d37b32 788 Int_t ddj = dd-TMath::Abs(di);
cc80f89e 789 Int_t sigstart = (ddj>0) ? -1 : 0;
790 for (Int_t sig = sigstart;sig<=1;sig+=2){
791 Int_t dj= sig*ddj;
792 AliCell *cell= GetCell(i+di,j+dj);
73042f01 793 Int_t signal = GetSignal(i+di,j+dj);
cc80f89e 794 if (cell==0) continue;
795 Int_t index[6];
796 index[0]=di;
797 index[1]=dj;
798 if (dd>2) {
799 SetBlockIndex(index); //adjust index to control
73042f01 800 if ( IsBorder(fIndex,i+index[2],j+index[3]) ||
801 IsBorder(fIndex,i+index[4],j+index[5])) {
802 cell->SetBorder(fIndex);
cc80f89e 803 continue;
804 }
805 }
73042f01 806 if ( signal<=fThreshold ){
cc80f89e 807 //if under threshold
73042f01 808 cell->SetThBorder(fIndex);
809 if (fBFit==kTRUE) AddToStack(i+di,j+dj,signal);
cc80f89e 810 continue;
811 }
812 naccepted++;
813 if (CheckIfDirBorder(fCurrentMaxX,fCurrentMaxY,i+di,j+dj) == kTRUE) {
73042f01 814 if (fBFit==kFALSE) AddToStack(i+di,j+dj,signal/2);
cc80f89e 815 continue;
816 }
73042f01 817 AddToStack(i+di,j+dj,signal);
cc80f89e 818
819 } //loop over sig dj
820 } //loop over di
821
822 }//loop over dd
823 } //if there is maximum
824 //if there exists more then 2 digits cluster
825 if (fStackIndex >2 ){
826 if (fBFit==kFALSE) GetClusterStatistic(c);
827 else GetClusterFit(c);
828 //write some important chracteristic area of cluster
829 //
830 Transform(&c);
831 //write cluster information to array
73042f01 832 TClonesArray &lclusters = *fClustersArray;
cc80f89e 833 new(lclusters[fIndex++]) AliDigitCluster(c);
834 // cout<<"fx="<<c.fX<<" fy"<<c.fY<<"\n";
835 }
836 } //lopp over all digits
837
73042f01 838 return fClustersArray;
cc80f89e 839}
840
841
842
843
844
845
73042f01 846void AliTPCClusterFinder::Adjacent(Int_t i,Int_t j)
cc80f89e 847{
848 //
849 //recursive agorithm program
850 //
851 if (fBDistType==kTRUE) {
852 Float_t delta = (i-fCurrentMaxX)*(i-fCurrentMaxX)/fCurrentSigmaX2;
853 delta+=(j-fCurrentMaxY)*(j-fCurrentMaxY)/fCurrentSigmaY2;
854 if (delta > fMulSigma2) {
73042f01 855 SetDirBorder(fIndex,i,j);
cc80f89e 856 return;
857 }
858 }
859 AliCell *cell = GetCell(i,j);
73042f01 860 Int_t signal = GetSignal(i,j);
861 Int_t q=signal;
862 cell->SetChecked(fIndex);
cc80f89e 863 if ( (q>fThreshold) || (fBFit==kTRUE)) AddToStack(i,j,q);
864 if ( q >fThreshold )
865 {
866
867 AliCell * newcel;
868 newcel = GetCell(i-1,j);
73042f01 869 if (newcel !=0) if (!newcel->IsChecked(fIndex) ) Adjacent(i-1,j);
cc80f89e 870 newcel = GetCell(i,j-1);
73042f01 871 if (newcel !=0) if (!newcel->IsChecked(fIndex) ) Adjacent(i,j-1);
cc80f89e 872 newcel = GetCell(i+1,j);
73042f01 873 if (newcel !=0) if (!newcel->IsChecked(fIndex) ) Adjacent(i+1,j);
cc80f89e 874 newcel = GetCell(i,j+1);
73042f01 875 if (newcel !=0) if (!newcel->IsChecked(fIndex) ) Adjacent(i,j+1);
cc80f89e 876 }
73042f01 877 else cell->SetThBorder(fIndex);
cc80f89e 878}
879
880
881
73042f01 882AliH2F * AliTPCClusterFinder::DrawHisto( const char *option=0,
cc80f89e 883 Float_t x1, Float_t x2, Float_t y1, Float_t y2)
884{
885 //
886 //draw digits in given array
887 //
888 //make digits histo
889 char ch[30];
890 sprintf(ch,"Cluster finder digits ");
891 if ( (fDimX<1)|| (fDimY<1)) {
892 return 0;
893 }
894 AliH2F * his = new AliH2F(ch,ch,fDimX,fX1,fX2,fDimY,fY1,fY2);
895 //set histogram values
896 for (Int_t i = 0; i<fDimX;i++)
897 for (Int_t j = 0; j<fDimY;j++){
898 Float_t x = ItoX(i);
899 Float_t y= JtoY(j);
900 his->Fill(x,y,GetSignal(i,j));
901 }
902 if (x1>=0) {
903 AliH2F *h2fsub = his->GetSubrange2d(x1,x2,y1,y2);
904 delete his;
905 his=h2fsub;
906 }
907 if (his==0) return 0;
908 if (option!=0) his->Draw(option);
909 else his->Draw();
910 return his;
911}
912
913
73042f01 914void AliTPCClusterFinder::DrawCluster(
cc80f89e 915 Int_t color, Int_t size, Int_t style)
916{
917
918 if (fClustersArray==0) return;
919 //draw marker for each of cluster
920 Int_t ncl=fClustersArray->GetEntriesFast();
921 for (Int_t i=0;i<ncl;i++){
67d37b32 922 AliComplexCluster *cl = (AliComplexCluster*)fClustersArray->UncheckedAt(i);
cc80f89e 923 TMarker * marker = new TMarker;
924 marker->SetX(cl->fX);
925 marker->SetY(cl->fY);
926 marker->SetMarkerSize(size);
927 marker->SetMarkerStyle(style);
928 marker->SetMarkerColor(color);
929 marker->Draw();
930 }
931}
932
933
934
73042f01 935AliH2F * AliTPCClusterFinder::DrawBorders( const char *option, AliH2F *h, Int_t type ,
cc80f89e 936 Float_t x1, Float_t x2, Float_t y1, Float_t y2)
937{
938 //
939 //draw digits in given array
940 //
941 //make digits histo
942 char ch[30];
943 sprintf(ch,"Cluster finder digits borders");
944 if ( (fDimX<1)|| (fDimY<1)) {
945 return 0;
946 }
947 AliH2F * his;
948 if (h!=0) his =h;
949 else his = new AliH2F(ch,ch,fDimX,fX1,fX2,fDimY,fY1,fY2);
950 //set histogram values
951 for (Int_t i = 0; i<fDimX;i++)
952 for (Int_t j = 0; j<fDimY;j++){
953 Float_t x = ItoX(i);
954 Float_t y= JtoY(j);
955 if (((type==1)||(type==0))&&IsMaximum(0,i,j)) his->Fill(x,y,16);
956 if (((type==3)||(type==0))&&(IsDirBorder(0,i,j))) his->Fill(x,y,8);
957 if (((type==4)||(type==0))&&(IsThBorder(0,i,j))) his->Fill(x,y,4);
958 if (((type==2)||(type==0))&&IsBorder(0,i,j)) his->Fill(x,y,1);
959
960 }
961
962 if (x1>=0) {
963 AliH2F *h2fsub = his->GetSubrange2d(x1,x2,y1,y2);
964 delete his;
965 his=h2fsub;
966 }
967 if (his==0) return 0;
968 if (option!=0) his->Draw(option);
969 else his->Draw();
970 return his;
971}