method Raw2SDigits added
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV2.cxx
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
49 ClassImp(AliPMDClusteringV2)
50
51 const Double_t AliPMDClusteringV2::fgkSqroot3by2=0.8660254;  // sqrt(3.)/2.
52
53 AliPMDClusteringV2::AliPMDClusteringV2():
54   fClno(0),
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 // ------------------------------------------------------------------------ //
68 AliPMDClusteringV2::~AliPMDClusteringV2()
69 {
70
71 }
72 // ------------------------------------------------------------------------ //
73 void 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 // ------------------------------------------------------------------------ //
197 void 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 // ------------------------------------------------------------------------ //
234 Int_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 // ------------------------------------------------------------------------ //
349 void 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
576 void 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 // ------------------------------------------------------------------------ //
826 Double_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 // ------------------------------------------------------------------------ //
832 void AliPMDClusteringV2::SetEdepCut(Float_t decut)
833 {
834   fCutoff = decut;
835 }
836 // ------------------------------------------------------------------------ //