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