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