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