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