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