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