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