16b332d6ef56222c4e62c1a47642764b2ddba92b
[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 /* $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 TObjarray  and cluster information is in TObjArray.
32    integer clno gives total number of clusters in the  supermodule.
33    fClusters is the  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 <TMath.h>
42 #include <TObjArray.h>
43 #include <TArrayI.h>
44
45 #include "AliPMDcludata.h"
46 #include "AliPMDcluster.h"
47 #include "AliPMDClustering.h"
48 #include "AliPMDClusteringV2.h"
49 #include "AliLog.h"
50
51 ClassImp(AliPMDClusteringV2)
52
53 const Double_t AliPMDClusteringV2::fgkSqroot3by2=0.8660254;  // sqrt(3.)/2.
54
55 AliPMDClusteringV2::AliPMDClusteringV2():
56   fPMDclucont(new TObjArray()),
57   fCutoff(0.0)
58 {
59   for(int i = 0; i < kNDIMX; i++)
60     {
61       for(int j = 0; j < kNDIMY; j++)
62         {
63           fCoord[0][i][j] = i+j/2.;
64           fCoord[1][i][j] = fgkSqroot3by2*j;
65         }
66     }
67 }
68 // ------------------------------------------------------------------------ //
69
70
71 AliPMDClusteringV2::AliPMDClusteringV2(const AliPMDClusteringV2& pmdclv2):
72   AliPMDClustering(pmdclv2),
73   fPMDclucont(0),
74   fCutoff(0)
75 {
76   // copy constructor
77   AliError("Copy constructor not allowed ");
78   
79 }
80 // ------------------------------------------------------------------------ //
81 AliPMDClusteringV2 &AliPMDClusteringV2::operator=(const AliPMDClusteringV2& /*pmdclv2*/)
82 {
83   // copy constructor
84   AliError("Assignment operator not allowed ");
85   return *this;
86 }
87 // ------------------------------------------------------------------------ //
88 AliPMDClusteringV2::~AliPMDClusteringV2()
89 {
90   delete fPMDclucont;
91 }
92 // ------------------------------------------------------------------------ //
93
94 void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn, 
95                                  Double_t celladc[48][96], TObjArray *pmdcont)
96 {
97   // main function to call other necessary functions to do clustering
98   //
99   AliPMDcluster *pmdcl = 0;
100
101   const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
102
103   Int_t    i, j, nmx1, incr, id, jd;
104   Int_t    ndimXr = 0;
105   Int_t    ndimYr = 0;
106   Int_t    celldataX[15], celldataY[15];
107   Float_t  clusdata[6];  
108   Double_t cutoff, ave;
109   Double_t edepcell[kNMX];
110
111
112   if (ismn < 12)
113     {
114       ndimXr = 96;
115       ndimYr = 48;
116     }
117   else if (ismn >= 12 && ismn <= 23)
118     {
119       ndimXr = 48;
120       ndimYr = 96;
121     }
122   
123   for (Int_t i =0; i < kNMX; i++)
124     {
125      edepcell[i] = 0.;
126     }
127     
128   for (id = 0; id < ndimXr; id++)
129     {
130       for (jd = 0; jd < ndimYr; jd++)
131         {
132           j = jd;
133           i = id + (ndimYr/2-1) - (jd/2);
134           Int_t ij = i + j*kNDIMX;
135           if (ismn < 12)
136             {
137               edepcell[ij]    = celladc[jd][id];
138             }
139           else if (ismn >= 12 && ismn <= 23)
140             {
141              edepcell[ij]    = celladc[id][jd];
142             }
143
144         }
145     }
146
147   Int_t iord1[kNMX];
148   TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
149   cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
150   ave  = 0.;
151   nmx1 = -1;
152
153   for(i = 0;i < kNMX; i++)
154     {
155       if(edepcell[i] > 0.) 
156         {
157           ave += edepcell[i];
158         }
159       if(edepcell[i] > cutoff )
160         {
161           nmx1++;
162         }
163     }
164   
165   AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
166   
167   if (nmx1 == 0) 
168     {
169       nmx1 = 1;
170     }
171   ave = ave/nmx1;
172   
173   AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
174                   kNMX,ave));
175   
176   incr = CrClust(ave, cutoff, nmx1,iord1, edepcell);
177   RefClust(incr,edepcell );
178   
179   Int_t nentries1 = fPMDclucont->GetEntries();
180   AliDebug(1,Form("Detector Plane = %d  Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
181   AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
182   for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
183     {
184       AliPMDcludata *pmdcludata = 
185         (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
186       Float_t cluXC    = pmdcludata->GetClusX();
187       Float_t cluYC    = pmdcludata->GetClusY();
188       Float_t cluADC   = pmdcludata->GetClusADC();
189       Float_t cluCELLS = pmdcludata->GetClusCells();
190       Float_t cluSIGX  = pmdcludata->GetClusSigmaX();
191       Float_t cluSIGY  = pmdcludata->GetClusSigmaY();
192       
193       Float_t cluY0    = ktwobysqrt3*cluYC;
194       Float_t cluX0    = cluXC - cluY0/2.;
195       
196       // 
197       // Cluster X centroid is back transformed
198       //
199       if (ismn < 12)
200         {
201           clusdata[0] = cluX0 - (24-1) + cluY0/2.;
202         }
203       else if (ismn  >= 12 && ismn <= 23)
204         {
205           clusdata[0] = cluX0 - (48-1) + cluY0/2.;
206         }         
207
208       clusdata[1]     = cluY0;
209       clusdata[2]     = cluADC;
210       clusdata[3]     = cluCELLS;
211       clusdata[4]     = cluSIGX;
212       clusdata[5]     = cluSIGY;
213       //
214       // Cells associated with a cluster
215       //
216       for (Int_t ihit = 0; ihit < 15; ihit++)
217         {
218           Int_t dummyXY = pmdcludata->GetCellXY(ihit);
219          
220           Int_t celldumY   = dummyXY%10000;
221           Int_t celldumX   = dummyXY/10000;
222           Float_t cellY    = (Float_t) (ktwobysqrt3*celldumY/10);
223           Float_t cellX    = (Float_t) (celldumX/10 - (celldumY/2.)/10);
224           
225           // 
226           // Cell X centroid is back transformed
227           //
228           if (ismn < 12)
229             {
230               celldataX[ihit] = (Int_t) (cellX - (24-1) + cellY/2.);
231             }
232           else if (ismn  >= 12 && ismn <= 23)
233             {
234               celldataX[ihit] = (Int_t) (cellX - (48-1) + cellY/2.);
235             }     
236           celldataY[ihit] = (Int_t) cellY;
237         }
238
239       pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
240       pmdcont->Add(pmdcl);
241     }
242   fPMDclucont->Clear();
243 }
244 // ------------------------------------------------------------------------ //
245 Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
246                                   Int_t iord1[], Double_t edepcell[])
247 {
248   // Does crude clustering
249   // Finds out only the big patch by just searching the
250   // connected cells
251   //
252
253   Int_t i,j,k,id1,id2,icl, numcell;
254   Int_t jd1,jd2, icell, cellcount;
255   Int_t clust[2][5000];
256   static Int_t neibx[6] = {1,0,-1,-1,0,1}, neiby[6] = {0,1,1,0,-1,-1};
257
258   // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
259   // cell. There are six neighbours.
260   // cellcount --- total number of cells having nonzero ener dep
261   // numcell --- number of cells in a given supercluster
262   
263   AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
264   
265   for (j=0; j < kNDIMX; j++)
266     {
267       for(k=0; k < kNDIMY; k++)
268         {
269           fInfocl[0][j][k] = 0;
270           fInfocl[1][j][k] = 0;
271         }
272     }
273  
274   for(i=0; i < kNMX; i++)
275     {
276       fInfcl[0][i] = -1;
277       
278       j  = iord1[i];
279       id2 = j/kNDIMX;
280       id1 = j-id2*kNDIMX;
281       
282       if(edepcell[j] <= cutoff)
283         {
284           fInfocl[0][id1][id2] = -1;
285         }
286     }
287   // ---------------------------------------------------------------
288   // crude clustering begins. Start with cell having largest adc
289   // count and loop over the cells in descending order of adc count
290   // ---------------------------------------------------------------
291   icl       = -1;
292   cellcount = -1;
293   for(icell=0; icell <= nmx1; icell++)
294     {
295       j  = iord1[icell];
296       id2 = j/kNDIMX;
297       id1 = j-id2*kNDIMX;
298       if(fInfocl[0][id1][id2] == 0 )
299         {
300           // ---------------------------------------------------------------
301           // icl -- cluster #, numcell -- # of cells in it, clust -- stores
302           // coordinates of the cells in a cluster, fInfocl[0][i1][i2] is 1 for
303           // primary and 2 for secondary cells,
304           // fInfocl[1][i1][i2] stores cluster #
305           // ---------------------------------------------------------------
306           icl++;
307           numcell = 0;
308           cellcount++;
309           fInfocl[0][id1][id2]  = 1;
310           fInfocl[1][id1][id2]  = icl;
311           fInfcl[0][cellcount]  = icl;
312           fInfcl[1][cellcount]  = id1;
313           fInfcl[2][cellcount]  = id2;
314
315           clust[0][numcell]     = id1;
316           clust[1][numcell]     = id2;
317           for(i = 1; i < 5000; i++)
318             {
319               clust[0][i] = -1;
320             }
321           // ---------------------------------------------------------------
322           // check for adc count in neib. cells. If ne 0 put it in this clust
323           // ---------------------------------------------------------------
324           for(i = 0; i < 6; i++)
325             {
326             jd1 = id1 + neibx[i];
327             jd2 = id2 + neiby[i];
328             if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
329                 fInfocl[0][jd1][jd2] == 0)
330               {
331                 numcell++;
332                 fInfocl[0][jd1][jd2] = 2;
333                 fInfocl[1][jd1][jd2] = icl;
334                 clust[0][numcell]    = jd1;
335                 clust[1][numcell]    = jd2;
336                 cellcount++;
337                 fInfcl[0][cellcount] = icl;
338                 fInfcl[1][cellcount] = jd1;
339                 fInfcl[2][cellcount] = jd2;
340               }
341             }
342           // ---------------------------------------------------------------
343           // check adc count for neighbour's neighbours recursively and
344           // if nonzero, add these to the cluster.
345           // ---------------------------------------------------------------
346           for(i = 1;i < 5000; i++)
347             { 
348               if(clust[0][i] != -1)
349                 {
350                   id1 = clust[0][i];
351                   id2 = clust[1][i];
352                   for(j = 0; j < 6 ; j++)
353                     {
354                       jd1 = id1 + neibx[j];
355                       jd2 = id2 + neiby[j];
356                       if( (jd1 >= 0 && jd1 < kNDIMX) && 
357                           (jd2 >= 0 && jd2 < kNDIMY) 
358                           && fInfocl[0][jd1][jd2] == 0 )
359                         {
360                           fInfocl[0][jd1][jd2] = 2;
361                           fInfocl[1][jd1][jd2] = icl;
362                           numcell++;
363                           clust[0][numcell]    = jd1;
364                           clust[1][numcell]    = jd2;
365                           cellcount++;
366                           fInfcl[0][cellcount] = icl;
367                           fInfcl[1][cellcount] = jd1;
368                           fInfcl[2][cellcount] = jd2;
369                         }
370                     }
371                 }
372             }
373         }
374     }
375   return cellcount;
376 }
377 // ------------------------------------------------------------------------ //
378   void AliPMDClusteringV2::RefClust(Int_t incr, Double_t edepcell[])
379 {
380   // Does the refining of clusters
381   // Takes the big patch and does gaussian fitting and
382   // finds out the more refined clusters
383
384   AliPMDcludata *pmdcludata = 0;
385   TArrayI testncl;
386   TArrayI testindex;
387   const Int_t kndim = 4500;
388   Int_t    i, j, k, i1, i2, id, icl, itest, ihld;
389   Int_t    ig, nsupcl, clno, clX,clY;
390   Int_t    clxy[15];
391   Int_t    ncl[kndim], iord[kndim];
392   Float_t  clusdata[6];
393   Double_t x1, y1, z1, x2, y2, z2, rr;
394   Double_t x[kndim], y[kndim], z[kndim];
395   Double_t xc[kndim], yc[kndim], zc[kndim], cells[kndim];
396   Double_t rcl[kndim], rcs[kndim];
397   
398   for(Int_t kk = 0; kk < 15; kk++)
399     {
400       if( kk < 6 )clusdata[kk] = 0.;
401     }
402    
403   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
404   // x, y and z store (x,y) coordinates of and energy deposited in a cell
405   // xc, yc store (x,y) coordinates of the cluster center
406   // zc stores the energy deposited in a cluster, rc is cluster radius
407
408   clno   = -1;
409   nsupcl = -1;
410   for(i = 0; i < 4500; i++)
411     {
412       ncl[i] = -1;
413     }
414   for(i = 0; i < incr; i++)
415     {
416       if(fInfcl[0][i] != nsupcl)
417         {
418           nsupcl++;
419         }
420       if (nsupcl > 4500) 
421         {
422           AliWarning("RefClust: Too many superclusters!");
423           nsupcl = 4500;
424           break;
425         }
426       ncl[nsupcl]++;
427     }
428   
429   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
430                   incr+1,nsupcl+1));
431   
432   id  = -1;
433   icl = -1;
434   for(i = 0; i < nsupcl; i++)
435     {
436       if(ncl[i] == 0)
437         {
438           id++;
439           icl++;
440           // one  cell super-clusters --> single cluster
441           // cluster center at the centyer of the cell
442           // cluster radius = half cell dimension
443           if (clno >= 5000) 
444             {
445               AliWarning("RefClust: Too many clusters! more than 5000");
446               return;
447             }
448           clno++;
449           i1          = fInfcl[1][id];
450           i2          = fInfcl[2][id];
451           Int_t i12 = i1 + i2*kNDIMX;
452           clusdata[0] = fCoord[0][i1][i2];
453           clusdata[1] = fCoord[1][i1][i2];
454           clusdata[2] = edepcell[i12];
455           clusdata[3] = 1.;
456           clusdata[4] = 0.0;
457           clusdata[5] = 0.0;
458           
459           //cell information
460           clX = (Int_t) fCoord[0][i1][i2]*10;
461           clY = (Int_t) fCoord[1][i1][i2]*10;
462           clxy[0] = clX*10000 + clY;
463           for(Int_t icltr = 1; icltr < 15; icltr++)
464             {
465               clxy[icltr] = -1;
466             }
467           pmdcludata  = new AliPMDcludata(clusdata,clxy);
468           fPMDclucont->Add(pmdcludata);
469           
470           
471         }
472       else if(ncl[i] == 1)
473         {
474           // two cell super-cluster --> single cluster
475           // cluster center is at ener. dep.-weighted mean of two cells
476           // cluster radius == half cell dimension
477           id++;
478           icl++;
479           if (clno >= 5000) 
480             {
481               AliWarning("RefClust: Too many clusters! more than 5000");
482               return;
483             }
484           clno++;
485           i1   = fInfcl[1][id];
486           i2   = fInfcl[2][id];
487           Int_t i12 = i1 + i2*kNDIMX;
488           
489           x1   = fCoord[0][i1][i2];
490           y1   = fCoord[1][i1][i2];
491           z1   = edepcell[i12];
492           
493           id++;
494           i1   = fInfcl[1][id];
495           i2   = fInfcl[2][id];
496           i12  = i1 + i2*kNDIMX;
497           
498           x2   = fCoord[0][i1][i2];
499           y2   = fCoord[1][i1][i2];
500           z2   = edepcell[i12];
501           
502           clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
503           clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
504           clusdata[2] = z1+z2;
505           clusdata[3] = 2.;
506           clusdata[4] = (TMath::Sqrt(z1*z2))/(z1+z2);
507           clusdata[5] = 0.0;
508
509           clX = (Int_t) x1*10;
510           clY = (Int_t) y1*10;
511           clxy[0] = clX*10000 + clY;
512           
513           clX = (Int_t) x2*10;
514           clY = (Int_t) y2*10;
515           clxy[1] = clX*10000 + clY;
516           
517           for(Int_t icltr = 2; icltr < 15; icltr++)
518             {
519               clxy[icltr] = -1;
520             }
521           pmdcludata  = new AliPMDcludata(clusdata, clxy);
522           fPMDclucont->Add(pmdcludata);
523         }
524       else{
525         id++;
526         iord[0] = 0;
527         // super-cluster of more than two cells - broken up into smaller
528         // clusters gaussian centers computed. (peaks separated by > 1 cell)
529         // Begin from cell having largest energy deposited This is first
530         // cluster center
531         // *****************************************************************
532         // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
533         // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
534         // SINCE WE EXPECT THE SUPERCLUSTER 
535         // TO BE A SINGLE CLUSTER
536         //*******************************************************************
537         
538         i1      = fInfcl[1][id];
539         i2      = fInfcl[2][id];
540         Int_t i12 = i1 + i2*kNDIMX;
541         
542         x[0]    = fCoord[0][i1][i2];
543         y[0]    = fCoord[1][i1][i2];
544         z[0]    = edepcell[i12];
545         
546         iord[0] = 0;
547         for(j = 1; j <= ncl[i]; j++)
548           {
549             
550             id++;
551             i1      = fInfcl[1][id];
552             i2      = fInfcl[2][id];
553             Int_t i12 = i1 + i2*kNDIMX;
554             iord[j] = j;
555             x[j]    = fCoord[0][i1][i2];
556             y[j]    = fCoord[1][i1][i2];
557             z[j]    = edepcell[i12];
558           }
559         
560         // arranging cells within supercluster in decreasing order
561         for(j = 1; j <= ncl[i];j++)
562           {
563             itest = 0;
564             ihld  = iord[j];
565             for(i1 = 0; i1 < j; i1++)
566               {
567                 if(itest == 0 && z[iord[i1]] < z[ihld])
568                   {
569                     itest = 1;
570                     for(i2 = j-1;i2 >= i1;i2--)
571                       {
572                         iord[i2+1] = iord[i2];
573                       }
574                     iord[i1] = ihld;
575                   }
576               }
577           }
578         
579         
580         // compute the number of clusters and their centers ( first
581         // guess )
582         // centers must be separated by cells having smaller ener. dep.
583         // neighbouring centers should be either strong or well-separated
584         ig     = 0;
585         xc[ig] = x[iord[0]];
586         yc[ig] = y[iord[0]];
587         zc[ig] = z[iord[0]];
588         for(j = 1; j <= ncl[i]; j++)
589           {
590             itest = -1;
591             x1    = x[iord[j]];
592             y1    = y[iord[j]];
593             for(k = 0; k <= ig; k++)
594               {
595                 x2 = xc[k];
596                 y2 = yc[k];
597                 rr = Distance(x1,y1,x2,y2);
598                 //************************************************************
599                 // finetuning cluster splitting
600                 // the numbers zc/4 and zc/10 may need to be changed. 
601                 // Also one may need to add one more layer because our 
602                 // cells are smaller in absolute scale
603                 //************************************************************
604                 
605                 
606                 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.) itest++;
607                 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.) itest++;
608                 if( rr >= 2.1)itest++;
609               }
610             
611             if(itest == ig)
612               {
613                 ig++;
614                 xc[ig] = x1;
615                 yc[ig] = y1;
616                 zc[ig] = z[iord[j]];
617               }
618           }
619         ClustDetails(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0],
620                      rcl[0], rcs[0], cells[0], testncl, testindex);
621         
622         Int_t pp = 0;
623         for(j = 0; j <= ig; j++)
624           { 
625             clno++;
626             if (clno >= 5000)
627               {
628                 AliWarning("RefClust: Too many clusters! more than 5000");
629                 return;
630               }
631             clusdata[0] = xc[j];
632             clusdata[1] = yc[j];
633             clusdata[2] = zc[j];
634             clusdata[4] = rcl[j];
635             clusdata[5] = rcs[j];
636             if(ig == 0)
637               {
638                 clusdata[3] = ncl[i];
639               }
640             else
641               {
642                 clusdata[3] = cells[j];
643               }
644             // cell information
645             Int_t ncellcls =  testncl[j];
646             if( ncellcls < 15 )
647               {
648                 for(Int_t kk = 1; kk <= ncellcls; kk++)
649                   {
650                     Int_t ll =  testindex[pp];
651                     clX = (Int_t) x[ll]*10;
652                     clY = (Int_t) y[ll]*10;
653                     clxy[kk-1] = (Int_t) clX*10000 + clY ;
654                     pp++;
655                   }
656                 for(Int_t icltr = ncellcls ; icltr < 15; icltr++)
657                   {
658                     clxy[icltr] = -1;
659                   }
660               }
661             pmdcludata = new AliPMDcludata(clusdata, clxy);
662             fPMDclucont->Add(pmdcludata);
663           }
664         testncl.Set(0);
665         testindex.Set(0);
666       }
667     }
668 }
669 // ------------------------------------------------------------------------ //
670 void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust, Double_t &x, 
671                                       Double_t &y, Double_t &z, Double_t &xc,
672                                       Double_t &yc, Double_t &zc,
673                                       Double_t &rcl, Double_t &rcs, 
674                                       Double_t &cells, TArrayI &testncl,
675                                       TArrayI &testindex)
676 {
677   // function begins
678   //
679
680   const Int_t kndim1 = 2000;
681   const Int_t kndim2 = 10;
682   const Int_t kndim3 = 400;
683   
684   Int_t    i, j, k, i1, i2;
685   Int_t    cluster[kndim1][kndim2], cell[kndim1][kndim3];
686   
687   Double_t x1, y1, x2, y2, rr, b, c, r1, r2;
688   Double_t sumx, sumy, sumxy, sumxx, sum, sum1, sumyy;
689   Double_t xx[kndim1], yy[kndim1], zz[kndim1];
690   Double_t xxc[kndim1], yyc[kndim1],clustcell[kndim3][kndim1];
691   Double_t str[kndim1], str1[kndim1],xcl[kndim1], ycl[kndim1], cln[kndim1];
692
693   for(i = 0; i <= nclust; i++)
694     {
695       xxc[i]  = *(&xc+i); 
696       yyc[i]  = *(&yc+i); 
697       str[i]  = 0.; 
698       str1[i] = 0.;
699     }
700   for(i = 0; i <= ncell; i++)
701     {
702       xx[i] = *(&x+i); 
703       yy[i] = *(&y+i); 
704       zz[i] = *(&z+i);
705     }
706
707   // INITIALIZE 
708
709   for(i = 0; i < kndim1; i++)
710     {
711       for(j = 0; j < kndim2; j++)
712         {
713           cluster[i][j] = 0;
714           
715         }
716       for(j = 0; j < kndim3; j++)
717         {
718           cell[i][j]      = 0;
719           clustcell[j][i] = 0.;
720         }
721     }
722
723
724   if(nclust > 0)
725     {
726       // more than one cluster
727       // checking cells shared between several  clusters.
728       // First check if the cell is within
729       // one cell unit ( nearest neighbour). Else, 
730       // if it is within 1.74 cell units ( next nearest )
731       // Else if it is upto 2 cell units etc.
732       
733       for (i = 0; i <= ncell; i++)
734         {
735           x1            = xx[i];
736           y1            = yy[i];
737           cluster[i][0] = 0;
738           // distance <= 1 cell unit
739           for(j = 0; j <= nclust; j++)
740             {
741               x2 = xxc[j];
742               y2 = yyc[j];
743               rr = Distance(x1, y1, x2, y2);
744               if(rr <= 1.)
745                 {
746                   cluster[i][0]++;
747                   i1             = cluster[i][0];
748                   cluster[i][i1] = j;
749                 }
750             }
751           // next nearest neighbour
752           if(cluster[i][0] == 0)
753             {
754               for(j=0; j<=nclust; j++)
755                 {
756                   x2 = xxc[j];
757                   y2 = yyc[j];
758                   rr = Distance(x1, y1, x2, y2);
759                   if(rr <= TMath::Sqrt(3.))
760                     {
761                       cluster[i][0]++;
762                       i1             = cluster[i][0];
763                       cluster[i][i1] = j;
764                     }
765                 }
766             }
767           // next-to-next nearest neighbour
768           if(cluster[i][0] == 0)
769             {
770               for(j=0; j<=nclust; j++)
771                 {
772                   x2 = xxc[j];
773                   y2 = yyc[j];
774                   rr = Distance(x1, y1, x2, y2);
775                   if(rr <= 2.)
776                     {
777                       cluster[i][0]++;
778                       i1             = cluster[i][0];
779                       cluster[i][i1] = j;
780                     }
781                 }
782             }
783           // one more
784           if(cluster[i][0] == 0)
785             {
786               for(j = 0; j <= nclust; j++)
787                 {
788                   x2 = xxc[j];
789                   y2 = yyc[j];
790                   rr = Distance(x1, y1, x2, y2);
791                   if(rr <= 2.7)
792                     {
793                       cluster[i][0]++;
794                       i1             = cluster[i][0];
795                       cluster[i][i1] = j;
796                     }
797                 }
798             }
799         }
800       
801       // computing cluster strength. Some cells are shared.
802       for(i = 0; i <= ncell; i++)
803         {
804           if(cluster[i][0] != 0)
805             {
806               i1 = cluster[i][0];
807               for(j = 1; j <= i1; j++)
808                 {
809                   i2       = cluster[i][j];
810                   str[i2] += zz[i]/i1;
811                 }
812             }
813         }
814       
815       for(k = 0; k < 5; k++)
816         {
817           for(i = 0; i <= ncell; i++)
818             {
819               if(cluster[i][0] != 0)
820                 {
821                   i1=cluster[i][0];
822                   sum=0.;
823                   for(j = 1; j <= i1; j++)
824                     {
825                       sum += str[cluster[i][j]];
826                     }
827                   
828                   for(j = 1; j <= i1; j++)
829                     {
830                       i2               = cluster[i][j]; 
831                       str1[i2]        +=  zz[i]*str[i2]/sum;
832                       clustcell[i2][i] = zz[i]*str[i2]/sum;
833                     }
834                 }
835             }
836           
837           
838           for(j = 0; j <= nclust; j++)
839             {
840               str[j]  = str1[j];
841               str1[j] = 0.;
842             }
843         }
844       
845       for(i = 0; i <= nclust; i++)
846         {
847           sumx = 0.;
848           sumy = 0.;
849           sum  = 0.;
850           sum1 = 0.;
851           for(j = 0; j <= ncell; j++)
852             {
853               if(clustcell[i][j] != 0)
854                 {
855                   sumx  +=  clustcell[i][j]*xx[j];
856                   sumy  +=  clustcell[i][j]*yy[j];
857                   sum   +=  clustcell[i][j];
858                   sum1  +=  clustcell[i][j]/zz[j];
859                 }
860             }
861           //** xcl and ycl are cluster centroid positions ( center of gravity )
862           
863           xcl[i] = sumx/sum;
864           ycl[i] = sumy/sum;
865           cln[i] = sum1;
866           sumxx = 0.;
867           sumyy = 0.;
868           sumxy = 0.;
869           for(j = 0; j <= ncell; j++)
870             {
871               sumxx += clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
872               sumyy += clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
873               sumxy += clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
874             }
875           b = sumxx+sumyy;
876           c = sumxx*sumyy-sumxy*sumxy;
877           // ******************r1 and r2 are major and minor axes ( r1 > r2 ). 
878           r1 = b/2.+TMath::Sqrt(b*b/4.-c);
879           r2 = b/2.-TMath::Sqrt(b*b/4.-c);
880           // final assignments to proper external variables
881           *(&xc + i) = xcl[i];
882           *(&yc + i) = ycl[i];
883           *(&zc + i) = str[i];
884           *(&cells + i) = cln[i];
885           *(&rcl+i) = r1;
886           *(&rcs+i) = r2;
887         }
888       
889       //To get the cell position in a cluster
890       
891       for(Int_t ii=0; ii<= ncell; ii++)
892         {
893           Int_t jj = cluster[ii][0]; 
894           for(Int_t kk=1; kk<= jj; kk++)
895             {
896               Int_t ll = cluster[ii][kk];
897               cell[ll][0]++;
898               cell[ll][cell[ll][0]] = ii;
899             }
900         }
901       
902       testncl.Set(nclust+1);
903       Int_t counter = 0;
904       
905       for(Int_t ii=0; ii <= nclust; ii++)
906         {
907           testncl[ii] =  cell[ii][0];
908           counter += testncl[ii];
909         }
910       testindex.Set(counter);
911       Int_t ll = 0;
912       for(Int_t ii=0; ii<= nclust; ii++)
913         {
914           for(Int_t jj = 1; jj<= testncl[ii]; jj++)
915             {
916               Int_t kk = cell[ii][jj];
917               testindex[ll] = kk;
918               ll++;
919             }
920         }
921       
922     }
923   else
924     {
925       sumx = 0.;
926       sumy = 0.;
927       sum  = 0.;
928       sum1 = 0.;
929       i    = 0;
930       for(j = 0; j <= ncell; j++)
931         {
932           sumx += zz[j]*xx[j];
933           sumy += zz[j]*yy[j];
934           sum  += zz[j];
935           sum1++;
936         }
937       xcl[i] = sumx/sum;
938       ycl[i] = sumy/sum;
939       cln[i] = sum1;
940       sumxx  = 0.;
941       sumyy  = 0.;
942       sumxy  = 0.;
943       for(j = 0; j <= ncell; j++)
944         {
945           sumxx += clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
946           sumyy += clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
947           sumxy += clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
948         }
949       b  = sumxx+sumyy;
950       c  = sumxx*sumyy-sumxy*sumxy;
951       r1 = b/2.+ TMath::Sqrt(b*b/4.-c);
952       r2 = b/2.- TMath::Sqrt(b*b/4.-c);
953       
954       // To get the cell position in a cluster
955       testncl.Set(nclust+1);
956       testindex.Set(ncell);
957       cell[0][0] = ncell;
958       testncl[0] = cell[0][0];
959       Int_t ll   = 0;
960       for(Int_t ii=1; ii<=ncell; ii++)
961         {
962           cell[0][ii]=ii;
963           //clustcell[0][ii]=1.;
964           Int_t kk = cell[0][ii];
965           testindex[ll] = kk;
966           ll++;
967         }
968       // final assignments
969       *(&xc + i)    = xcl[i];
970       *(&yc + i)    = ycl[i];
971       *(&zc + i)    = str[i];
972       *(&cells + i) = cln[i];
973       *(&rcl+i)     = r1;
974       *(&rcs+i)     = r2;
975     }
976 }
977
978 // ------------------------------------------------------------------------ //
979 Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
980                                       Double_t x2, Double_t y2)
981 {
982   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
983 }
984 // ------------------------------------------------------------------------ //
985 void AliPMDClusteringV2::SetEdepCut(Float_t decut)
986 {
987   fCutoff = decut;
988 }
989 // ------------------------------------------------------------------------ //