new version clustering algo
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV2.cxx
CommitLineData
8c7250c5 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// //
18// Source File : PMDClusteringV2.cxx //
19// //
20// clustering code for alice pmd //
21// //
22//-----------------------------------------------------//
23
24/* --------------------------------------------------------------------
25 Code developed by S. C. Phatak, Institute of Physics,
26 Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
27 ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
28 builds up superclusters and breaks them into clusters. The input is
29 in array fEdepCell[kNDIMX][kNDIMY] and cluster information is in array
30 fClusters[5][5000]. integer fClno gives total number of clusters in the
31 supermodule.
32
33 fEdepCell, fClno and fClusters are the only global ( public ) variables.
34 Others are local ( private ) to the code.
35 At the moment, the data is read for whole detector ( all supermodules
36 and pmd as well as cpv. This will have to be modify later )
37 LAST UPDATE : October 23, 2002
38-----------------------------------------------------------------------*/
39
40#include "Riostream.h"
41#include <TObjArray.h>
42#include <stdio.h>
43
44#include "AliPMDcluster.h"
45#include "AliPMDClustering.h"
46#include "AliPMDClusteringV2.h"
47#include "AliLog.h"
48
49ClassImp(AliPMDClusteringV2)
50
51const Double_t AliPMDClusteringV2::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
52
53AliPMDClusteringV2::AliPMDClusteringV2():
54 fCutoff(0.0)
55{
56 for(int i = 0; i < kNDIMX; i++)
57 {
58 for(int j = 0; j < kNDIMY; j++)
59 {
60 fCoord[0][i][j] = i+j/2.;
61 fCoord[1][i][j] = fgkSqroot3by2*j;
62 fEdepCell[i][j] = 0;
63 }
64 }
65}
66// ------------------------------------------------------------------------ //
67AliPMDClusteringV2::~AliPMDClusteringV2()
68{
69
70}
71// ------------------------------------------------------------------------ //
72void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96], TObjArray *pmdcont)
73{
74 // main function to call other necessary functions to do clustering
75 //
76 AliPMDcluster *pmdcl = 0;
77
78 Int_t i, i1, i2, j, nmx1, incr, id, jd;
79 Int_t celldataX[15], celldataY[15];
80 Float_t clusdata[6];
81 Double_t cutoff, ave;
82
83 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
84
85 Int_t ndimXr =0;
86 Int_t ndimYr =0;
87
88 if (ismn < 12)
89 {
90 ndimXr = 96;
91 ndimYr = 48;
92 }
93 else if (ismn >= 12 && ismn <= 23)
94 {
95 ndimXr = 48;
96 ndimYr = 96;
97 }
98
99 for (Int_t i =0; i < kNDIMX; i++)
100 {
101 for (Int_t j =0; j < kNDIMY; j++)
102 {
103 fEdepCell[i][j] = 0;
104 }
105 }
106
107
108 for (id = 0; id < ndimXr; id++)
109 {
110 for (jd = 0; jd < ndimYr; jd++)
111 {
112 j=jd;
113 i=id+(ndimYr/2-1)-(jd/2);
114
115 if (ismn < 12)
116 {
117 fEdepCell[i][j] = celladc[jd][id];
118 }
119 else if (ismn >= 12 && ismn <= 23)
120 {
121 fEdepCell[i][j] = celladc[id][jd];
122 }
123
124 }
125 }
126
127 Order(); // order the data
128 cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
129 ave=0.;
130 nmx1=-1;
131
132 for(j=0;j<kNMX; j++)
133 {
134 i1 = fIord[0][j];
135 i2 = fIord[1][j];
136 if (fEdepCell[i1][i2] > 0.) {ave = ave + fEdepCell[i1][i2];}
137 if (fEdepCell[i1][i2] > cutoff ) nmx1 = nmx1 + 1;
138 }
139 // nmx1 --- number of cells having ener dep >= cutoff
140
141 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
142
143 if (nmx1 == 0) nmx1 = 1;
144 ave=ave/nmx1;
145
146 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
147 kNMX,ave));
148
149 incr = CrClust(ave, cutoff, nmx1);
150 RefClust(incr);
151
152 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, fClno));
153
154 for(i1=0; i1<=fClno; i1++)
155 {
156 Float_t cluXC = (Float_t) fClusters[0][i1];
157 Float_t cluYC = (Float_t) fClusters[1][i1];
158 Float_t cluADC = (Float_t) fClusters[2][i1];
159 Float_t cluCELLS = (Float_t) fClusters[3][i1];
160 Float_t sigmaX = (Float_t) fClusters[4][i1];
161 Float_t sigmaY = (Float_t) fClusters[5][i1];
162 Float_t cluY0 = ktwobysqrt3*cluYC;
163 Float_t cluX0 = cluXC - cluY0/2.;
164 //
165 // Cluster X centroid is back transformed
166 //
167 if (ismn < 12)
168 {
169 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
170 }
171 else if (ismn >= 12 && ismn <= 23)
172 {
173 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
174 }
175
176 clusdata[1] = cluY0;
177 clusdata[2] = cluADC;
178 clusdata[3] = cluCELLS;
179 clusdata[4] = sigmaX;
180 clusdata[5] = sigmaY;
181
182 //
183 // Cells associated with a cluster
184 //
185 for (Int_t ihit = 0; ihit < 15; ihit++)
186 {
187 celldataX[ihit] = 1; // dummy nos. -- will be changed
188 celldataY[ihit] = 1; // dummy nos. -- will be changed
189 }
190
191 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
192 pmdcont->Add(pmdcl);
193 }
194}
195// ------------------------------------------------------------------------ //
196void AliPMDClusteringV2::Order()
197{
198 // Sorting algorithm
199 // sorts the ADC values from higher to lower
200 //
201 double dd[kNMX];
202 // matrix fEdepCell converted into
203 // one dimensional array dd. adum a place holder for double
204 int i, j, i1, i2, iord1[kNMX];
205 // information of
206 // ordering is stored in iord1, original array not ordered
207 //
208 // define arrays dd and iord1
209 for(i1=0; i1 < kNDIMX; i1++)
210 {
211 for(i2=0; i2 < kNDIMY; i2++)
212 {
213 i = i1 + i2*kNDIMX;
214 iord1[i] = i;
215 dd[i] = fEdepCell[i1][i2];
216 }
217 }
218 // sort and store sorting information in iord1
219
220 TMath::Sort(kNMX,dd,iord1);
221
222 // store the sorted information in fIord for later use
223 for(i=0; i<kNMX; i++)
224 {
225 j = iord1[i];
226 i2 = j/kNDIMX;
227 i1 = j-i2*kNDIMX;
228 fIord[0][i]=i1;
229 fIord[1][i]=i2;
230 }
231}
232// ------------------------------------------------------------------------ //
233Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1)
234{
235 // Does crude clustering
236 // Finds out only the big patch by just searching the
237 // connected cells
238 //
239
240 int i,j,k,id1,id2,icl, numcell;
241 int jd1,jd2, icell, cellcount;
242 int clust[2][5000];
243 static int neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
244
245 // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
246 // cell. There are six neighbours.
247 // cellcount --- total number of cells having nonzero ener dep
248 // numcell --- number of cells in a given supercluster
249 // ofstream ofl0("cells_loc",ios::out);
250 // initialize fInfocl[2][kNDIMX][kNDIMY]
251
252 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
253
254 for (j=0; j < kNDIMX; j++){
255 for(k=0; k < kNDIMY; k++){
256 fInfocl[0][j][k] = 0;
257 fInfocl[1][j][k] = 0;
258 }
259 }
260 for(i=0; i < kNMX; i++){
261 fInfcl[0][i] = -1;
262 id1=fIord[0][i];
263 id2=fIord[1][i];
264 if(fEdepCell[id1][id2] <= cutoff){fInfocl[0][id1][id2]=-1;}
265 }
266 // ---------------------------------------------------------------
267 // crude clustering begins. Start with cell having largest adc
268 // count and loop over the cells in descending order of adc count
269 // ---------------------------------------------------------------
270 icl=-1;
271 cellcount=-1;
272 for(icell=0; icell <= nmx1; icell++){
273 id1=fIord[0][icell];
274 id2=fIord[1][icell];
275 if(fInfocl[0][id1][id2] == 0 ){
276 // ---------------------------------------------------------------
277 // icl -- cluster #, numcell -- # of cells in it, clust -- stores
278 // coordinates of the cells in a cluster, fInfocl[0][i1][i2] is 1 for
279 // primary and 2 for secondary cells,
280 // fInfocl[1][i1][i2] stores cluster #
281 // ---------------------------------------------------------------
282 icl=icl+1;
283 numcell=0;
284 cellcount = cellcount + 1;
285 fInfocl[0][id1][id2]=1;
286 fInfocl[1][id1][id2]=icl;
287 fInfcl[0][cellcount]=icl;
288 fInfcl[1][cellcount]=id1;
289 fInfcl[2][cellcount]=id2;
290
291 clust[0][numcell]=id1;
292 clust[1][numcell]=id2;
293 for(i=1; i<5000; i++)clust[0][i] = -1;
294 // ---------------------------------------------------------------
295 // check for adc count in neib. cells. If ne 0 put it in this clust
296 // ---------------------------------------------------------------
297 for(i=0; i<6; i++){
298 jd1=id1+neibx[i];
299 jd2=id2+neiby[i];
300 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
301 fInfocl[0][jd1][jd2] == 0){
302 numcell=numcell+1;
303 fInfocl[0][jd1][jd2]=2;
304 fInfocl[1][jd1][jd2]=icl;
305 clust[0][numcell]=jd1;
306 clust[1][numcell]=jd2;
307 cellcount=cellcount+1;
308 fInfcl[0][cellcount]=icl;
309 fInfcl[1][cellcount]=jd1;
310 fInfcl[2][cellcount]=jd2;
311 }
312 }
313 // ---------------------------------------------------------------
314 // check adc count for neighbour's neighbours recursively and
315 // if nonzero, add these to the cluster.
316 // ---------------------------------------------------------------
317 for(i=1;i < 5000;i++){
318 if(clust[0][i] != -1){
319 id1=clust[0][i];
320 id2=clust[1][i];
321 for(j=0; j<6 ; j++){
322 jd1=id1+neibx[j];
323 jd2=id2+neiby[j];
324 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
325 fInfocl[0][jd1][jd2] == 0 ){
326 fInfocl[0][jd1][jd2] = 2;
327 fInfocl[1][jd1][jd2] = icl;
328 numcell = numcell + 1;
329 clust[0][numcell] = jd1;
330 clust[1][numcell] = jd2;
331 cellcount = cellcount+1;
332 fInfcl[0][cellcount] = icl;
333 fInfcl[1][cellcount] = jd1;
334 fInfcl[2][cellcount] = jd2;
335 }
336 }
337 }
338 }
339 }
340 }
341 // for(icell=0; icell<=cellcount; icell++){
342 // ofl0 << fInfcl[0][icell] << " " << fInfcl[1][icell] << " " <<
343 // fInfcl[2][icell] << endl;
344 // }
345 return cellcount;
346}
347// ------------------------------------------------------------------------ //
348void AliPMDClusteringV2::RefClust(Int_t incr)
349{
350 // Does the refining of clusters
351 // Takes the big patch and does gaussian fitting and
352 // finds out the more refined clusters
353 //
354
355 const Int_t kndim = 4500;
356
357 int i, j, k, i1, i2, id, icl, itest;
358 int ihld;
359 int ig, nsupcl;
360 int ncl[kndim], iord[kndim];
361
362 double x1, y1, z1, x2, y2, z2;
363 double rr;
364
365 double x[kndim], y[kndim], z[kndim];
366 double xc[kndim], yc[kndim], zc[kndim], cells[kndim];
367 double rcl[kndim], rcs[kndim];
368
369 // fClno counts the final clusters
370 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
371 // x, y and z store (x,y) coordinates of and energy deposited in a cell
372 // xc, yc store (x,y) coordinates of the cluster center
373 // zc stores the energy deposited in a cluster
374 // rc is cluster radius
375 // finally the cluster information is put in 2-dimensional array clusters
376 // ofstream ofl1("checking.5",ios::app);
377
378 fClno = -1;
379 nsupcl = -1;
380 for(i=0; i<4500; i++){ncl[i]=-1;}
381 for(i=0; i<incr; i++){
382 if(fInfcl[0][i] != nsupcl){ nsupcl=nsupcl+1; }
383 if (nsupcl > 4500) {
384 AliWarning("RefClust: Too many superclusters!");
385 nsupcl = 4500;
386 break;
387 }
388 ncl[nsupcl]=ncl[nsupcl]+1;
389 }
390
391 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
392 incr+1,nsupcl+1));
393
394 id=-1;
395 icl=-1;
396 for(i=0; i<nsupcl; i++){
397 if(ncl[i] == 0){
398 id++;
399 icl++;
400 // one cell super-clusters --> single cluster
401 // cluster center at the centyer of the cell
402 // cluster radius = half cell dimension
403 if (fClno >= 5000) {
404 AliWarning("RefClust: Too many clusters! more than 5000");
405 return;
406 }
407 fClno++;
408 i1 = fInfcl[1][id];
409 i2 = fInfcl[2][id];
410 fClusters[0][fClno] = fCoord[0][i1][i2];
411 fClusters[1][fClno] = fCoord[1][i1][i2];
412 fClusters[2][fClno] = fEdepCell[i1][i2];
413 fClusters[3][fClno] = 1.;
414 fClusters[4][fClno] = 0.0;
415 fClusters[5][fClno] = 0.0;
416 //ofl1 << icl << " " << fCoord[0][i1][i2] << " " << fCoord[1][i1][i2] <<
417 //" " << fEdepCell[i1][i2] << " " << fClusters[3][fClno] <<endl;
418 }else if(ncl[i] == 1){
419 // two cell super-cluster --> single cluster
420 // cluster center is at ener. dep.-weighted mean of two cells
421 // cluster radius == half cell dimension
422 id++;
423 icl++;
424 if (fClno >= 5000) {
425 AliWarning("RefClust: Too many clusters! more than 5000");
426 return;
427 }
428 fClno++;
429 i1 = fInfcl[1][id];
430 i2 = fInfcl[2][id];
431 x1 = fCoord[0][i1][i2];
432 y1 = fCoord[1][i1][i2];
433 z1 = fEdepCell[i1][i2];
434
435 id++;
436 i1 = fInfcl[1][id];
437 i2 = fInfcl[2][id];
438 x2 = fCoord[0][i1][i2];
439 y2 = fCoord[1][i1][i2];
440 z2 = fEdepCell[i1][i2];
441
442 fClusters[0][fClno] = (x1*z1+x2*z2)/(z1+z2);
443 fClusters[1][fClno] = (y1*z1+y2*z2)/(z1+z2);
444 fClusters[2][fClno] = z1+z2;
445 fClusters[3][fClno] = 2.;
446 fClusters[4][fClno] = sqrt(z1*z2)/(z1+z2);
447 fClusters[5][fClno] = 0; // sigma large nonzero, sigma small zero
448
449 //ofl1 << icl << " " << fClusters[0][fClno] << " " << fClusters[1][fClno]
450 // << " " << fClusters[2][fClno] << " " <<fClusters[3][fClno] <<endl;
451 }
452 else{
453 id = id + 1;
454 iord[0] = 0;
455 // super-cluster of more than two cells - broken up into smaller
456 // clusters gaussian centers computed. (peaks separated by > 1 cell)
457 // Begin from cell having largest energy deposited This is first
458 // cluster center
459 // *****************************************************************
460 // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
461 // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
462 // SINCE WE EXPECT THE SUPERCLUSTER
463 // TO BE A SINGLE CLUSTER
464 //*******************************************************************
465
466 i1 = fInfcl[1][id];
467 i2 = fInfcl[2][id];
468 x[0] = fCoord[0][i1][i2];
469 y[0] = fCoord[1][i1][i2];
470 z[0] = fEdepCell[i1][i2];
471 iord[0] = 0;
472 for(j=1;j<=ncl[i];j++){
473
474 id = id + 1;
475 i1 = fInfcl[1][id];
476 i2 = fInfcl[2][id];
477 iord[j] = j;
478 x[j] = fCoord[0][i1][i2];
479 y[j] = fCoord[1][i1][i2];
480 z[j] = fEdepCell[i1][i2];
481 }
482 // arranging cells within supercluster in decreasing order
483 for(j=1;j<=ncl[i];j++)
484 {
485 itest = 0;
486 ihld = iord[j];
487 for(i1=0; i1<j; i1++)
488 {
489 if(itest == 0 && z[iord[i1]] < z[ihld])
490 {
491 itest = 1;
492 for(i2=j-1;i2>=i1;i2--)
493 {
494 iord[i2+1] = iord[i2];
495 }
496 iord[i1] = ihld;
497 }
498 }
499 }
500
501 // compute the number of clusters and their centers ( first
502 // guess )
503 // centers must be separated by cells having smaller ener. dep.
504 // neighbouring centers should be either strong or well-separated
505 ig = 0;
506 xc[ig] = x[iord[0]];
507 yc[ig] = y[iord[0]];
508 zc[ig] = z[iord[0]];
509 for(j=1;j<=ncl[i];j++){
510 itest = -1;
511 x1 = x[iord[j]];
512 y1 = y[iord[j]];
513 for(k=0;k<=ig;k++){
514 x2 = xc[k];
515 y2 = yc[k];
516 rr = Distance(x1,y1,x2,y2);
517 //***************************************************************
518 // finetuning cluster splitting
519 // the numbers zc/4 and zc/10 may need to be changed.
520 // Also one may need to add one more layer because our
521 // cells are smaller in absolute scale
522 //****************************************************************
523
524
525 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)
526 itest++;
527 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)
528 itest++;
529 if( rr >= 2.1)itest++;
530 }
531 if(itest == ig){
532 ig++;
533 xc[ig] = x1;
534 yc[ig] = y1;
535 zc[ig] = z[iord[j]];
536 }
537 }
538
539 ClustDetails(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0],
540 rcl[0], rcs[0], cells[0]);
541
542 icl = icl + ig + 1;
543
544 for(j=0; j<=ig; j++)
545 {
546 if (fClno >= 5000)
547 {
548 AliWarning("RefClust: Too many clusters! more than 5000");
549 return;
550 }
551 fClno++;
552 fClusters[0][fClno] = xc[j];
553 fClusters[1][fClno] = yc[j];
554 fClusters[2][fClno] = zc[j];
555 fClusters[4][fClno] = rcl[j];
556 fClusters[5][fClno] = rcs[j];
557 if(ig == 0)
558 {
559 fClusters[3][fClno] = ncl[i];
560 }
561 else
562 {
563 fClusters[3][fClno] = cells[j];
564 }
565 }
566
567
568 }
569 }
570}
571
572
573// ------------------------------------------------------------------------ //
574
575void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust,
576 Double_t &x, Double_t &y, Double_t &z,
577 Double_t &xc, Double_t &yc, Double_t &zc,
578 Double_t &rcl, Double_t &rcs,
579 Double_t &cells)
580{
581 // function begins
582 //
583
584 const Int_t kndim1 = 4500;
585 const Int_t kndim2 = 10;
586 const Int_t kndim3 = 100;
587
588 int i, j, k, i1, i2;
589 int cluster[kndim1][kndim2];
590
591 double x1, y1, x2, y2, rr;
592 double sumx, sumy, sumxy, sumxx;
593 double sum, sum1, sumyy;
594 double b, c, r1, r2;
595
596 double xx[kndim1], yy[kndim1], zz[kndim1];
597 double xxc[kndim1], yyc[kndim1];
598
599 double str[kndim1];
600
601 double str1[kndim1];
602 double xcl[kndim1], ycl[kndim1], cln[kndim1];
603 double clustcell[kndim1][kndim3];
604
605 for(i=0; i<=nclust; i++){
606 xxc[i]=*(&xc+i);
607 yyc[i]=*(&yc+i);
608 str[i]=0.;
609 str1[i]=0.;
610 }
611 for(i=0; i<=ncell; i++){
612 xx[i]=*(&x+i);
613 yy[i]=*(&y+i);
614 zz[i]=*(&z+i);
615 }
616 // INITIALIZE
617 for(i=0; i<4500; i++){
618 for(j=0; j<100; j++){
619 clustcell[i][j]=0.;
620 }
621 }
622
623 // INITIALIZE
624 for(i=0;i<4500;i++){
625 for(j=0;j<10;j++){
626 cluster[i][j]=0;
627 }
628 }
629
630
631 if(nclust > 0){
632 // more than one cluster
633 // checking cells shared between several clusters.
634 // First check if the cell is within
635 // one cell unit ( nearest neighbour). Else,
636 // if it is within 1.74 cell units ( next nearest )
637 // Else if it is upto 2 cell units etc.
638
639 for (i=0; i<=ncell; i++){
640 x1 = xx[i];
641 y1 = yy[i];
642 cluster[i][0] = 0;
643 // distance <= 1 cell unit
644 for(j=0; j<=nclust; j++)
645 {
646 x2 = xxc[j];
647 y2 = yyc[j];
648 rr = Distance(x1, y1, x2, y2);
649 if(rr <= 1.)
650 {
651 cluster[i][0]++;
652 i1 = cluster[i][0];
653 cluster[i][i1] = j;
654 }
655 }
656 // next nearest neighbour
657 if(cluster[i][0] == 0)
658 {
659 for(j=0; j<=nclust; j++)
660 {
661 x2 = xxc[j];
662 y2 = yyc[j];
663 rr = Distance(x1, y1, x2, y2);
664 if(rr <= sqrt(3.))
665 {
666 cluster[i][0]++;
667 i1 = cluster[i][0];
668 cluster[i][i1] = j;
669 }
670 }
671 }
672 // next-to-next nearest neighbour
673 if(cluster[i][0] == 0)
674 {
675 for(j=0; j<=nclust; j++)
676 {
677 x2 = xxc[j];
678 y2 = yyc[j];
679 rr = Distance(x1, y1, x2, y2);
680 if(rr <= 2.)
681 {
682 cluster[i][0]++;
683 i1 = cluster[i][0];
684 cluster[i][i1] = j;
685 }
686 }
687 }
688 // one more
689 if(cluster[i][0] == 0)
690 {
691 for(j=0; j<=nclust; j++)
692 {
693 x2 = xxc[j];
694 y2 = yyc[j];
695 rr = Distance(x1, y1, x2, y2);
696 if(rr <= 2.7)
697 {
698 cluster[i][0]++;
699 i1 = cluster[i][0];
700 cluster[i][i1] = j;
701 }
702 }
703 }
704 }
705
706
707 // computing cluster strength. Some cells are shared.
708 for(i=0; i<=ncell; i++){
709 if(cluster[i][0] != 0){
710 i1 = cluster[i][0];
711 for(j=1; j<=i1; j++){
712 i2 = cluster[i][j];
713 str[i2] = str[i2]+zz[i]/i1;
714 }
715 }
716 }
717
718 for(k=0; k<5; k++)
719 {
720 for(i=0; i<=ncell; i++)
721 {
722 if(cluster[i][0] != 0)
723 {
724 i1=cluster[i][0];
725 sum=0.;
726 for(j=1; j<=i1; j++)
727 {
728 sum=sum+str[cluster[i][j]];
729 }
730
731 for(j=1; j<=i1; j++)
732 {
733 i2 = cluster[i][j];
734 str1[i2] = str1[i2] + zz[i]*str[i2]/sum;
735 clustcell[i2][i] = zz[i]*str[i2]/sum;
736 }
737 }
738 }
739
740
741 for(j=0; j<=nclust; j++)
742 {
743 str[j]=str1[j];
744 str1[j]=0.;
745 }
746 }
747
748 for(i=0; i<=nclust; i++){
749 sumx = 0.;
750 sumy = 0.;
751 sum = 0.;
752 sum1 = 0.;
753 for(j=0; j<=ncell; j++){
754 if(clustcell[i][j] != 0){
755 sumx = sumx+clustcell[i][j]*xx[j];
756 sumy = sumy+clustcell[i][j]*yy[j];
757 sum = sum+clustcell[i][j];
758 sum1 = sum1+clustcell[i][j]/zz[j];
759 }
760 }
761 //***** xcl and ycl are cluster centroid positions ( center of gravity )
762
763 xcl[i] = sumx/sum;
764 ycl[i] = sumy/sum;
765 cln[i] = sum1;
766 sumxx = 0.;
767 sumyy = 0.;
768 sumxy = 0.;
769 for(j=0; j<=ncell; j++){
770 sumxx = sumxx+clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
771 sumyy = sumyy+clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
772 sumxy = sumxy+clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
773 }
774 b = sumxx+sumyy;
775 c = sumxx*sumyy-sumxy*sumxy;
776 // ******************r1 and r2 are major and minor axes ( r1 > r2 ).
777 r1 = b/2.+sqrt(b*b/4.-c);
778 r2 = b/2.-sqrt(b*b/4.-c);
779 // final assignments to proper external variables
780 *(&xc + i) = xcl[i];
781 *(&yc + i) = ycl[i];
782 *(&zc + i) = str[i];
783 *(&cells + i) = cln[i];
784 *(&rcl+i) = r1;
785 *(&rcs+i) = r2;
786 }
787 }else{
788 sumx = 0.;
789 sumy = 0.;
790 sum = 0.;
791 sum1 = 0.;
792 i = 0;
793 for(j=0; j<=ncell; j++){
794 sumx = sumx+zz[j]*xx[j];
795 sumy = sumy+zz[j]*yy[j];
796 sum = sum+zz[j];
797 sum1 = sum1+1.;
798 }
799 xcl[i] = sumx/sum;
800 ycl[i] = sumy/sum;
801 cln[i] = sum1;
802 sumxx = 0.;
803 sumyy = 0.;
804 sumxy = 0.;
805 for(j=0; j<=ncell; j++){
806 sumxx = sumxx+clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
807 sumyy = sumyy+clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
808 sumxy = sumxy+clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
809 }
810 b = sumxx+sumyy;
811 c = sumxx*sumyy-sumxy*sumxy;
812 r1 = b/2.+sqrt(b*b/4.-c);
813 r2 = b/2.-sqrt(b*b/4.-c);
814 // final assignments
815 *(&xc + i) = xcl[i];
816 *(&yc + i) = ycl[i];
817 *(&zc + i) = str[i];
818 *(&cells + i) = cln[i];
819 *(&rcl+i) = r1;
820 *(&rcs+i) = r2;
821 }
822}
823
824// ------------------------------------------------------------------------ //
825Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
826 Double_t x2, Double_t y2)
827{
828 return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
829}
830// ------------------------------------------------------------------------ //
831void AliPMDClusteringV2::SetEdepCut(Float_t decut)
832{
833 fCutoff = decut;
834}
835// ------------------------------------------------------------------------ //