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