warning fixed
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV1.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 : PMDClusteringV1.cxx, Version 00      //
21 //                                                     //
22 //  Date   : September 26 2002                         //
23 //                                                     //
24 //  clustering code for alice pmd                      //
25 //                                                     //
26 //-----------------------------------------------------//
27
28 /* --------------------------------------------------------------------
29    Code developed by S. C. Phatak, Institute of Physics,
30    Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
31    ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
32    builds up superclusters and breaks them into clusters. The input is
33    in array edepcell[kNMX] and cluster information is in a
34    TObjarray. Integer clno gives total number of clusters in the
35    supermodule.
36
37    fClusters is the only global ( public ) variables.
38    Others are local ( private ) to the code.
39    At the moment, the data is read for whole detector ( all supermodules
40    and pmd as well as cpv. This will have to be modify later )
41    LAST UPDATE  :  October 23, 2002
42 -----------------------------------------------------------------------*/
43
44 #include <Riostream.h>
45 #include <TMath.h>
46 #include <TNtuple.h>
47 #include <TObjArray.h>
48 #include "TRandom.h"
49 #include <stdio.h>
50
51 #include "AliPMDcludata.h"
52 #include "AliPMDcluster.h"
53 #include "AliPMDClustering.h"
54 #include "AliPMDClusteringV1.h"
55 #include "AliLog.h"
56
57
58 ClassImp(AliPMDClusteringV1)
59
60 const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254;  // sqrt(3.)/2.
61
62 AliPMDClusteringV1::AliPMDClusteringV1():
63   fPMDclucont(new TObjArray()),
64   fCutoff(0.0)
65 {
66   for(Int_t i = 0; i < kNDIMX; i++)
67     {
68       for(Int_t j = 0; j < kNDIMY; j++)
69         {
70           fCoord[0][i][j] = i+j/2.;
71           fCoord[1][i][j] = fgkSqroot3by2*j;
72         }
73     }
74 }
75 // ------------------------------------------------------------------------ //
76 AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
77   AliPMDClustering(pmdclv1),
78   fPMDclucont(0),
79   fCutoff(0)
80 {
81   // copy constructor
82   AliError("Copy constructor not allowed ");
83   
84 }
85 // ------------------------------------------------------------------------ //
86 AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
87 {
88   // copy constructor
89   AliError("Assignment operator not allowed ");
90   return *this;
91 }
92 // ------------------------------------------------------------------------ //
93 AliPMDClusteringV1::~AliPMDClusteringV1()
94 {
95   delete fPMDclucont;
96 }
97 // ------------------------------------------------------------------------ //
98 void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, 
99                                  Double_t celladc[48][96], TObjArray *pmdcont)
100 {
101   // main function to call other necessary functions to do clustering
102   //
103
104   AliPMDcluster *pmdcl = 0;
105
106   const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
107   const Int_t kNmaxCell   = 19;     // # of cells surrounding a cluster center
108
109   Int_t    i,  j, nmx1, incr, id, jd;
110   Int_t    celldataX[kNmaxCell], celldataY[kNmaxCell];
111   Float_t  clusdata[6];
112   Double_t cutoff, ave;
113   Double_t edepcell[kNMX];
114   
115   
116   Double_t *cellenergy = new Double_t [11424];
117   
118
119   // ndimXr and ndimYr are different because of different module size
120
121   Int_t ndimXr = 0;
122   Int_t ndimYr = 0;
123
124   if (ismn < 12)
125     {
126       ndimXr = 96;
127       ndimYr = 48;
128     }
129   else if (ismn >= 12 && ismn <= 23)
130     {
131       ndimXr = 48;
132       ndimYr = 96;
133     }
134
135   for (i =0; i < 11424; i++)
136   {
137       cellenergy[i] = 0.;
138   }
139
140
141   Int_t kk = 0;
142   for (i = 0; i < kNDIMX; i++)
143     {
144       for (j = 0; j < kNDIMY; j++)
145         {
146           edepcell[kk] = 0.;
147           kk++;
148         }
149     }
150
151   for (id = 0; id < ndimXr; id++)
152     {
153       for (jd = 0; jd < ndimYr; jd++)
154         {
155           j = jd;
156           i = id+(ndimYr/2-1)-(jd/2);
157
158           Int_t ij = i + j*kNDIMX;
159           
160           if (ismn < 12)
161             {
162               cellenergy[ij]    = celladc[jd][id];
163             }
164           else if (ismn >= 12 && ismn <= 23)
165             {
166               cellenergy[ij]    = celladc[id][jd];
167             }
168         }
169     }
170   
171   for (i = 0; i < kNMX; i++)
172   {
173     edepcell[i] = cellenergy[i];
174   }
175
176   delete [] cellenergy;
177
178   Int_t iord1[kNMX];
179   TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
180   cutoff = fCutoff;                       // cutoff to discard cells
181   ave  = 0.;
182   nmx1 = -1;
183   for(i = 0;i < kNMX; i++)
184     {
185       if(edepcell[i] > 0.) 
186         {
187           ave += edepcell[i];
188         }
189       if(edepcell[i] > cutoff )
190         {
191           nmx1++;
192         }
193     }
194   
195   AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
196
197   if (nmx1 == 0) nmx1 = 1;
198   ave = ave/nmx1;
199   AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
200                   kNMX,ave));
201   
202   incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
203   RefClust(incr,edepcell);
204   Int_t nentries1 = fPMDclucont->GetEntries();
205   AliDebug(1,Form("Detector Plane = %d  Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
206   AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
207   
208   for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
209     {
210       AliPMDcludata *pmdcludata = 
211         (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
212       Float_t cluXC    = pmdcludata->GetClusX();
213       Float_t cluYC    = pmdcludata->GetClusY();
214       Float_t cluADC   = pmdcludata->GetClusADC();
215       Float_t cluCELLS = pmdcludata->GetClusCells();
216       Float_t cluSIGX  = pmdcludata->GetClusSigmaX();
217       Float_t cluSIGY  = pmdcludata->GetClusSigmaY();
218       
219       Float_t cluY0    = ktwobysqrt3*cluYC;
220       Float_t cluX0    = cluXC - cluY0/2.;
221       
222       // 
223       // Cluster X centroid is back transformed
224       //
225
226       if (ismn < 12)
227         {
228           clusdata[0] = cluX0 - (24-1) + cluY0/2.;
229         }
230       else if ( ismn >= 12 && ismn <= 23)
231         {
232           clusdata[0] = cluX0 - (48-1) + cluY0/2.;
233         }         
234       
235
236       clusdata[1]     = cluY0;
237       clusdata[2]     = cluADC;
238       clusdata[3]     = cluCELLS;
239       clusdata[4]     = cluSIGX;
240       clusdata[5]     = cluSIGY;
241       
242       //
243       // Cells associated with a cluster
244       //
245
246       for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
247         {
248           Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
249           Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
250
251           if (ismn < 12)
252             {
253               celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
254             }
255           else if (ismn >= 12 && ismn <= 23)
256             {
257               celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
258             }
259           
260           celldataY[ihit] = cellcol;
261         }
262
263
264       pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
265       pmdcont->Add(pmdcl);
266     }
267   
268   fPMDclucont->Delete();
269 }
270 // ------------------------------------------------------------------------ //
271 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
272                                   Int_t iord1[], Double_t edepcell[])
273 {
274   // Does crude clustering 
275   // Finds out only the big patch by just searching the
276   // connected cells
277   //
278   const Int_t kndim = 4609;
279   Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
280   Int_t jd1,jd2, icell, cellcount;
281   static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
282
283   AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
284   
285   for (j = 0; j < kNDIMX; j++)
286     {
287       for(k = 0; k < kNDIMY; k++)
288         {
289           fInfocl[0][j][k] = 0;
290           fInfocl[1][j][k] = 0;
291         }
292     }
293   for(i=0; i < kNMX; i++)
294     {
295       fInfcl[0][i] = -1;
296       
297       j  = iord1[i];
298       id2 = j/kNDIMX;
299       id1 = j-id2*kNDIMX;
300
301       if(edepcell[j] <= cutoff)
302         {
303           fInfocl[0][id1][id2] = -1;
304         }
305     }
306
307   // ---------------------------------------------------------------
308   // crude clustering begins. Start with cell having largest adc
309   // count and loop over the cells in descending order of adc count
310   // ---------------------------------------------------------------
311
312   icl       = -1;
313   cellcount = -1;
314
315   for(icell = 0; icell <= nmx1; icell++)
316     {
317       j  = iord1[icell];
318       id2 = j/kNDIMX;
319       id1 = j-id2*kNDIMX;
320
321       if(fInfocl[0][id1][id2] == 0 )
322         {
323           icl++;
324           numcell = 0;
325           cellcount++; 
326           fInfocl[0][id1][id2] = 1;
327           fInfocl[1][id1][id2] = icl;
328           fInfcl[0][cellcount] = icl;
329           fInfcl[1][cellcount] = id1;
330           fInfcl[2][cellcount] = id2;
331
332           clust[0][numcell] = id1;
333           clust[1][numcell] = id2;
334           
335           for(i = 1; i < kndim; i++)
336             {
337               clust[0][i]=0;
338             }
339           // ---------------------------------------------------------------
340           // check for adc count in neib. cells. If ne 0 put it in this clust
341           // ---------------------------------------------------------------
342           for(i = 0; i < 6; i++)
343             {
344               jd1 = id1 + neibx[i];
345               jd2 = id2 + neiby[i];
346               if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
347                   fInfocl[0][jd1][jd2] == 0)
348                 {
349                   numcell++;
350                   fInfocl[0][jd1][jd2] = 2;
351                   fInfocl[1][jd1][jd2] = icl;
352                   clust[0][numcell]    = jd1;
353                   clust[1][numcell]    = jd2;
354                   cellcount++;
355                   fInfcl[0][cellcount] = icl;
356                   fInfcl[1][cellcount] = jd1;
357                   fInfcl[2][cellcount] = jd2;
358                 }
359             }
360           // ---------------------------------------------------------------
361           // check adc count for neighbour's neighbours recursively and
362           // if nonzero, add these to the cluster.
363           // ---------------------------------------------------------------
364           for(i = 1; i < kndim;i++)
365             {
366               if(clust[0][i] != 0)
367                 {
368                   id1 = clust[0][i];
369                   id2 = clust[1][i];
370                   for(j = 0; j < 6 ; j++)
371                     {
372                       jd1 = id1 + neibx[j];
373                       jd2 = id2 + neiby[j];
374                       if( (jd1 >= 0 && jd1 < kNDIMX) && 
375                           (jd2 >= 0 && jd2 < kNDIMY) &&
376                           fInfocl[0][jd1][jd2] == 0 )
377                         {
378                           fInfocl[0][jd1][jd2] = 2;
379                           fInfocl[1][jd1][jd2] = icl;
380                           numcell++;
381                           clust[0][numcell]    = jd1;
382                           clust[1][numcell]    = jd2;
383                           cellcount++;
384                           fInfcl[0][cellcount] = icl;
385                           fInfcl[1][cellcount] = jd1;
386                           fInfcl[2][cellcount] = jd2;
387                         }
388                     }
389                 }
390             }
391         }
392     }
393   return cellcount;
394 }
395 // ------------------------------------------------------------------------ //
396 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
397 {
398   // Does the refining of clusters
399   // Takes the big patch and does gaussian fitting and
400   // finds out the more refined clusters
401   //
402   
403   AliPMDcludata *pmdcludata = 0;
404   
405   const Int_t kNmaxCell   = 19;    // # of cells surrounding a cluster center
406   
407   Int_t ndim = incr + 1;
408   
409   Int_t    *ncl  = 0x0;
410   Int_t    *clxy = 0x0;  
411   Int_t    i12, i22;
412   Int_t    i, j, k, i1, i2, id, icl,  itest,ihld, ig, nsupcl,clno, t1, t2;
413   Float_t  clusdata[6];
414   Double_t x1, y1, z1, x2, y2, z2, rr;
415   
416   ncl   = new Int_t [ndim];
417   clxy  = new Int_t [kNmaxCell];
418   
419   // Initialisation  
420   for(i = 0; i<ndim; i++)
421     {
422       ncl[i]  = -1; 
423       if (i < 6) clusdata[i] = 0.;
424       if (i < kNmaxCell) clxy[i]    = 0;
425     }
426
427   // clno counts the final clusters
428   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
429   // x, y and z store (x,y) coordinates of and energy deposited in a cell
430   // xc, yc store (x,y) coordinates of the cluster center
431   // zc stores the energy deposited in a cluster
432   // rc is cluster radius
433   
434   clno  = -1;
435   nsupcl = -1;
436
437   for(i = 0; i <= incr; i++)
438     {
439       if(fInfcl[0][i] != nsupcl)
440         {
441           nsupcl++;
442         }
443       if (nsupcl > ndim) 
444         {
445           AliWarning("RefClust: Too many superclusters!");
446           nsupcl = ndim;
447           break;
448         }
449       ncl[nsupcl]++;
450     }
451   
452   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
453                   incr+1,nsupcl+1));
454   id  = -1;
455   icl = -1;
456   
457   for(i = 0; i <= nsupcl; i++) 
458     {
459       if(ncl[i] == 0)
460         {
461           id++;
462           icl++;
463           if (clno >= 4608) 
464             {
465               AliWarning("RefClust: Too many clusters! more than 4608");
466               return;
467             }
468           clno++;
469           i1 = fInfcl[1][id];
470           i2 = fInfcl[2][id];
471           
472           i12 = i1 + i2*kNDIMX;
473           
474           clusdata[0] = fCoord[0][i1][i2];
475           clusdata[1] = fCoord[1][i1][i2];
476           clusdata[2] = edepcell[i12];
477           clusdata[3] = 1.;
478           clusdata[4] = 0.5;
479           clusdata[5] = 0.0;
480
481           clxy[0] = i1*10000 + i2;
482           
483           for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
484             {
485               clxy[icltr] = -1;
486             }
487           pmdcludata  = new AliPMDcludata(clusdata,clxy);
488           fPMDclucont->Add(pmdcludata);
489         }
490       else if(ncl[i] == 1) 
491         {
492           id++;
493           icl++;
494           if (clno >= 4608) 
495             {
496               AliWarning("RefClust: Too many clusters! more than 4608");
497               return;
498             }
499           clno++;
500           i1   = fInfcl[1][id];
501           i2   = fInfcl[2][id];
502           i12  = i1 + i2*kNDIMX;
503
504           x1   = fCoord[0][i1][i2];
505           y1   = fCoord[1][i1][i2];
506           z1   = edepcell[i12];
507
508           clxy[0] = i1*10000 + i2;
509           
510           id++;
511           i1   = fInfcl[1][id];
512           i2   = fInfcl[2][id];
513
514           i22  = i1 + i2*kNDIMX;
515           x2   = fCoord[0][i1][i2];
516           y2   = fCoord[1][i1][i2];
517           z2   = edepcell[i22];
518
519           clxy[1] = i1*10000 + i2;
520           
521
522           for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
523             {
524               clxy[icltr] = -1;
525             }
526           
527           clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
528           clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
529           clusdata[2] = z1+z2;
530           clusdata[3] = 2.;
531           clusdata[4] = 0.5;
532           clusdata[5] = 0.0;
533           pmdcludata  = new AliPMDcludata(clusdata,clxy);
534           fPMDclucont->Add(pmdcludata);
535         }
536       else
537         {
538           
539           Int_t    *iord, *tc, *t;
540           Double_t *x, *y, *z, *xc, *yc, *zc;
541
542           iord = new Int_t [ncl[i]+1];
543           tc   = new Int_t [ncl[i]+1];
544           t    = new Int_t [ncl[i]+1];
545           
546           x    = new Double_t [ncl[i]+1];
547           y    = new Double_t [ncl[i]+1];
548           z    = new Double_t [ncl[i]+1];
549           xc   = new Double_t [ncl[i]+1];
550           yc   = new Double_t [ncl[i]+1];
551           zc   = new Double_t [ncl[i]+1];
552           
553           for( k = 0; k < ncl[i]+1; k++)
554             {
555               iord[k] = -1;
556               t[k]    = -1;
557               tc[k]   = -1;
558               x[k]    = -1;
559               y[k]    = -1;
560               z[k]    = -1;
561               xc[k]   = -1;
562               yc[k]   = -1;
563               zc[k]   = -1;
564             }
565           id++;
566           // super-cluster of more than two cells - broken up into smaller
567           // clusters gaussian centers computed. (peaks separated by > 1 cell)
568           // Begin from cell having largest energy deposited This is first
569           // cluster center
570           i1      = fInfcl[1][id];
571           i2      = fInfcl[2][id];
572           i12     = i1 + i2*kNDIMX;
573           
574           x[0]    = fCoord[0][i1][i2];
575           y[0]    = fCoord[1][i1][i2];
576           z[0]    = edepcell[i12];
577           t[0]    = i1*10000 + i2;
578           
579
580           iord[0] = 0;
581           for(j = 1; j <= ncl[i]; j++)
582             {
583               id++;
584               i1      = fInfcl[1][id];
585               i2      = fInfcl[2][id];
586               i12     = i1 + i2*kNDIMX;
587
588               iord[j] = j;
589               x[j]    = fCoord[0][i1][i2];
590               y[j]    = fCoord[1][i1][i2];
591               z[j]    = edepcell[i12];
592               t[j]    = i1*10000 + i2;
593
594             }
595           
596           // arranging cells within supercluster in decreasing order
597           
598           for(j = 1;j <= ncl[i]; j++)
599             {
600               itest = 0;
601               ihld  = iord[j];
602               for(i1 = 0; i1 < j; i1++)
603                 {
604                   if(itest == 0 && z[iord[i1]] < z[ihld])
605                     {
606                       itest = 1;
607                       for(i2 = j-1; i2 >= i1; i2--)
608                         {
609                           iord[i2+1] = iord[i2];
610                         }
611                       iord[i1] = ihld;
612                     }
613                 }
614             }
615           /* MODIFICATION PART STARTS (Tapan July 2008)
616              iord[0] is the cell with highest ADC in the crude-cluster
617              ig is the number of local maxima in the crude-cluster
618              For the higest peak we make ig=0 which means first local maximum.
619              Next we go down in terms of the ADC sequence and find out if any
620              more of the cells form local maxima. The definition of local 
621              maxima is that all its neighbours are of less ADC compared to it. 
622           */
623           ig = 0;
624           xc[ig] = x[iord[0]];
625           yc[ig] = y[iord[0]];
626           zc[ig] = z[iord[0]];
627           tc[ig] = t[iord[0]];
628           Int_t ivalid = 0,  icount = 0;
629           
630           for(j=1;j<=ncl[i];j++)
631             {
632               x1  = x[iord[j]];
633               y1  = y[iord[j]]; 
634               z1  = z[iord[j]];
635               t1  = t[iord[j]];
636               rr=Distance(x1,y1,xc[ig],yc[ig]);
637               
638               // Check the cells which are outside the neighbours (rr>1.2)
639               if(rr>1.2 ) 
640                 {
641                   ivalid=0;
642                   icount=0;
643                   for(Int_t j1=1;j1<j;j1++)
644                     {
645                       icount++;
646                       Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);              
647                       if(rr1>1.2) ivalid++;
648                     }
649                   if(ivalid == icount && z1>0.5*zc[ig])
650                     {
651                       ig++;
652                       xc[ig]=x1; 
653                       yc[ig]=y1; 
654                       zc[ig]=z1;
655                       tc[ig]=t1;
656                     }
657                 }         
658             }
659           
660           icl=icl+ig+1;
661           
662           //  We use simple Gaussian weighting. (Tapan Jan 2005)
663           // compute the number of cells belonging to each cluster.
664           // cell can be shared between several clusters  
665           // in the ratio of cluster energy deposition
666           // To calculate: 
667           // (1) number of cells belonging to a cluster (ig) and 
668           // (2) total ADC of the cluster (ig) 
669           // (3) x and y positions of the cluster
670           
671           
672           Int_t *cellCount;
673           Int_t **cellXY;
674           
675           Int_t    *status;
676           Double_t *totaladc, *totaladc2, *ncell,*weight;
677           Double_t *xclust, *yclust, *sigxclust, *sigyclust;
678           Double_t *ax, *ay, *ax2, *ay2;
679           
680           
681           status    = new Int_t [ncl[i]+1];
682           cellXY    = new Int_t *[ncl[i]+1];
683           
684           cellCount = new Int_t [ig+1];
685           totaladc  = new Double_t [ig+1];
686           totaladc2 = new Double_t [ig+1];
687           ncell     = new Double_t [ig+1];
688           weight    = new Double_t [ig+1];
689           xclust    = new Double_t [ig+1];
690           yclust    = new Double_t [ig+1];
691           sigxclust = new Double_t [ig+1];
692           sigyclust = new Double_t [ig+1];
693           ax        = new Double_t [ig+1];
694           ay        = new Double_t [ig+1];
695           ax2       = new Double_t [ig+1];
696           ay2       = new Double_t [ig+1];
697           
698           for(j = 0; j < ncl[i]+1; j++) 
699             {
700               status[j]     = 0;
701               cellXY[j] = new Int_t[ig+1];
702             }
703           //initialization
704           for(Int_t kcl = 0; kcl < ig+1; kcl++)
705             {
706               cellCount[kcl] = 0;
707               totaladc[kcl]  = 0.;
708               totaladc2[kcl] = 0.;
709               ncell[kcl]     = 0.;
710               weight[kcl]    = 0.;      
711               xclust[kcl]    = 0.; 
712               yclust[kcl]    = 0.;
713               sigxclust[kcl] = 0.; 
714               sigyclust[kcl] = 0.;
715               ax[kcl]        = 0.;      
716               ay[kcl]        = 0.;      
717               ax2[kcl]       = 0.;      
718               ay2[kcl]       = 0.;    
719               for(j = 0; j < ncl[i]+1; j++)
720                 {
721                   cellXY[j][kcl] = 0;
722                 }
723             }
724           Double_t sumweight, gweight; 
725           
726           for(j = 0;j <= ncl[i]; j++)
727             {
728               x1 = x[iord[j]];
729               y1 = y[iord[j]];
730               z1 = z[iord[j]];
731               t1 = t[iord[j]];
732               
733               for(Int_t kcl=0; kcl<=ig; kcl++)
734                 {
735                   x2 = xc[kcl];
736                   y2 = yc[kcl];
737                   rr = Distance(x1,y1,x2,y2);
738                   t2 = tc[kcl];           
739                   
740                   if(rr==0)
741                     {
742                       ncell[kcl] = 1.;
743                       totaladc[kcl]  = z1;
744                       totaladc2[kcl]  = pow(z1,2);
745                       ax[kcl] =  x1 * z1;
746                       ay[kcl] =  y1 * z1;
747                       ax2[kcl]=  0.;
748                       ay2[kcl]=  0.;
749                       status[j] = 1;
750                     }
751                 }
752             }
753           
754           for(j = 0; j <= ncl[i]; j++)
755             {
756               Int_t   maxweight = 0;
757               Double_t     max  = 0.;
758               
759               if(status[j] == 0)
760                 { 
761                   x1 = x[iord[j]]; 
762                   y1 = y[iord[j]];
763                   z1 = z[iord[j]];
764                   t1 = t[iord[j]];
765                   sumweight = 0.;
766
767                   for(Int_t kcl = 0; kcl <= ig; kcl++)
768                     {
769                       x2 = xc[kcl]; 
770                       y2 = yc[kcl]; 
771                       rr = Distance(x1,y1,x2,y2);
772                       gweight     = exp(-(rr*rr)/(2*(1.2*1.2)));
773                       weight[kcl] = zc[kcl] * gweight;
774                       sumweight   = sumweight + weight[kcl];
775                       
776                       if(weight[kcl] > max)
777                         {
778                           max       =  weight[kcl];
779                           maxweight =  kcl;
780                         }
781                     }
782                   
783                   cellXY[cellCount[maxweight]][maxweight] = iord[j];
784                                   
785                   cellCount[maxweight]++;
786                   
787                   for(Int_t kcl = 0; kcl <= ig; kcl++)
788                     {
789                       x2 = xc[kcl];
790                       y2 = yc[kcl];
791                       rr = Distance(x1,y1,x2,y2);
792                       t2 = tc[kcl];
793                       // For calculating weighted mean:
794                       // Weighted_mean = (\sum w_i x_i) / (\sum w_i)
795                       
796                       if(sumweight>0.)
797                         {
798                           totaladc[kcl] = totaladc[kcl] + z1*(weight[kcl]/sumweight);
799                           ncell[kcl]    = ncell[kcl]    + (weight[kcl]/sumweight);        
800                           ax[kcl]       = ax[kcl]       + (x1 * z1 *weight[kcl]/sumweight);
801                           ay[kcl]       = ay[kcl]       + (y1 * z1 *weight[kcl]/sumweight);
802                           
803                           // For calculating weighted sigma:
804                           // Wieghted_sigma= ((\sum w_i)/((\sum w_i)^2 - \sum (w_i^2))) \sum w_i (x_i - \Hat\mu)^2
805                           // I assume here x1,y1 are cluster centers, otherwise life becomes too difficult
806                           totaladc2[kcl] = totaladc2[kcl] +  pow((z1*(weight[kcl]/sumweight)),2);
807                           ax2[kcl] = ax2[kcl] + (z1 *weight[kcl]/sumweight) * pow((x1-x2),2);
808                           ay2[kcl] = ay2[kcl] + (z1 *weight[kcl]/sumweight) * pow((y1-y2),2);
809                         }
810                     }
811                 }
812             }
813           
814           for(Int_t kcl = 0; kcl <= ig; kcl++)
815             {
816
817               if(totaladc[kcl]>0){
818                 if(totaladc[kcl]>0.)xclust[kcl] = (ax[kcl])/ totaladc[kcl];
819                 if(totaladc[kcl]>0.)yclust[kcl] = (ay[kcl])/ totaladc[kcl];
820                 
821                 if(totaladc[kcl]>0.)sigxclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ax2[kcl];
822                 if(totaladc[kcl]>0.)sigyclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ay2[kcl];
823               } 
824             
825               for(j = 0; j < cellCount[kcl]; j++) clno++; 
826               
827               if (clno >= 4608) 
828                 {
829                   AliWarning("RefClust: Too many clusters! more than 4608");
830                   return;
831                 }
832               clusdata[0] = xclust[kcl];
833               clusdata[1] = yclust[kcl];
834               clusdata[2] = totaladc[kcl];
835               clusdata[3] = ncell[kcl];
836               if(sigxclust[kcl] > sigyclust[kcl]) 
837                 {
838                   clusdata[4] = pow(sigxclust[kcl],0.5);
839                   clusdata[5] = pow(sigyclust[kcl],0.5);
840                 }
841               else
842                 {
843                   clusdata[4] = pow(sigyclust[kcl],0.5);
844                   clusdata[5] = pow(sigxclust[kcl],0.5);
845                 }
846
847               clxy[0] = tc[kcl];
848
849               Int_t Ncell=1;
850               for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
851                 {
852                   if(ii<18) 
853                     {   
854                       clxy[Ncell] = t[cellXY[ii][kcl]];
855                       Ncell++;
856                     }
857                 } 
858
859               pmdcludata = new AliPMDcludata(clusdata,clxy);
860               fPMDclucont->Add(pmdcludata);
861             }
862           
863           delete [] iord;
864           delete [] tc;   
865           delete [] t;
866           delete [] x;
867           delete [] y;
868           delete [] z;
869           delete [] xc;
870           delete [] yc;
871           delete [] zc;
872           
873           delete [] cellCount;
874           for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
875           
876           delete [] status;
877           delete [] totaladc;
878           delete [] totaladc2;
879           delete [] ncell;
880           delete [] xclust;
881           delete [] yclust;
882           delete [] sigxclust;
883           delete [] sigyclust;
884           delete [] ax;
885           delete [] ay;
886           delete [] ax2;
887           delete [] ay2;
888           delete [] weight;
889         }
890     }
891   delete [] ncl;
892   delete [] clxy;
893 }
894 // ------------------------------------------------------------------------ //
895 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1, 
896                                       Double_t x2, Double_t y2)
897 {
898   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
899 }
900 // ------------------------------------------------------------------------ //
901 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
902 {
903   fCutoff = decut;
904 }
905 // ------------------------------------------------------------------------ //