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