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