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