e399286285352179c7756333f7eec9e9c693a17e
[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           /* MODIFICATION PART STARTS (Tapan July 2008)
655              iord[0] is the cell with highest ADC in the crude-cluster
656              ig is the number of local maxima in the crude-cluster
657              For the higest peak we make ig=0 which means first local maximum.
658              Next we go down in terms of the ADC sequence and find out if any
659              more of the cells form local maxima. The definition of local 
660              maxima is that all its neighbours are of less ADC compared to it. 
661           */
662           ig = 0;
663           xc[ig] = x[iord[0]];
664           yc[ig] = y[iord[0]];
665           zc[ig] = z[iord[0]];
666           tc[ig] = t[iord[0]];
667           Int_t ivalid = 0,  icount = 0;
668           
669           for(j=1;j<=ncl[i];j++)
670             {
671               x1  = x[iord[j]];
672               y1  = y[iord[j]]; 
673               z1  = z[iord[j]];
674               t1  = t[iord[j]];
675               rr=Distance(x1,y1,xc[ig],yc[ig]);
676               
677               // Check the cells which are outside the neighbours (rr>1.2)
678               if(rr>1.2 ) 
679                 {
680                   ivalid=0;
681                   icount=0;
682                   for(Int_t j1=1;j1<j;j1++)
683                     {
684                       icount++;
685                       Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);              
686                       if(rr1>1.2) ivalid++;
687                     }
688                   if(ivalid == icount && z1>0.5*zc[ig])
689                     {
690                       ig++;
691                       xc[ig]=x1; 
692                       yc[ig]=y1; 
693                       zc[ig]=z1;
694                       tc[ig]=t1;
695                     }
696                 }         
697             }
698           
699           icl=icl+ig+1;
700           
701           //  We use simple Gaussian weighting. (Tapan Jan 2005)
702           // compute the number of cells belonging to each cluster.
703           // cell can be shared between several clusters  
704           // in the ratio of cluster energy deposition
705           // To calculate: 
706           // (1) number of cells belonging to a cluster (ig) and 
707           // (2) total ADC of the cluster (ig) 
708           // (3) x and y positions of the cluster
709           
710           
711           Int_t *cellCount;
712           Int_t **cellXY;
713           
714           Int_t    *status;
715           Double_t *totaladc, *totaladc2, *ncell,*weight;
716           Double_t *xclust, *yclust, *sigxclust, *sigyclust;
717           Double_t *ax, *ay, *ax2, *ay2;
718           
719           
720           status    = new Int_t [ncl[i]+1];
721           cellXY    = new Int_t *[ncl[i]+1];
722           
723           cellCount = new Int_t [ig+1];
724           totaladc  = new Double_t [ig+1];
725           totaladc2 = new Double_t [ig+1];
726           ncell     = new Double_t [ig+1];
727           weight    = new Double_t [ig+1];
728           xclust    = new Double_t [ig+1];
729           yclust    = new Double_t [ig+1];
730           sigxclust = new Double_t [ig+1];
731           sigyclust = new Double_t [ig+1];
732           ax        = new Double_t [ig+1];
733           ay        = new Double_t [ig+1];
734           ax2       = new Double_t [ig+1];
735           ay2       = new Double_t [ig+1];
736           
737           for(j = 0; j < ncl[i]+1; j++) 
738             {
739               status[j]     = 0;
740               cellXY[j] = new Int_t[ig+1];
741             }
742           //initialization
743           for(Int_t kcl = 0; kcl < ig+1; kcl++)
744             {
745               cellCount[kcl] = 0;
746               totaladc[kcl]  = 0.;
747               totaladc2[kcl] = 0.;
748               ncell[kcl]     = 0.;
749               weight[kcl]    = 0.;      
750               xclust[kcl]    = 0.; 
751               yclust[kcl]    = 0.;
752               sigxclust[kcl] = 0.; 
753               sigyclust[kcl] = 0.;
754               ax[kcl]        = 0.;      
755               ay[kcl]        = 0.;      
756               ax2[kcl]       = 0.;      
757               ay2[kcl]       = 0.;    
758               for(j = 0; j < ncl[i]+1; j++)
759                 {
760                   cellXY[j][kcl] = 0;
761                 }
762             }
763           Double_t sumweight, gweight; 
764           
765           for(j = 0;j <= ncl[i]; j++)
766             {
767               x1 = x[iord[j]];
768               y1 = y[iord[j]];
769               z1 = z[iord[j]];
770               t1 = t[iord[j]];
771               
772               for(Int_t kcl=0; kcl<=ig; kcl++)
773                 {
774                   x2 = xc[kcl];
775                   y2 = yc[kcl];
776                   rr = Distance(x1,y1,x2,y2);
777                   t2 = tc[kcl];           
778                   
779                   if(rr==0)
780                     {
781                       ncell[kcl]     = 1.;
782                       totaladc[kcl]  = z1;
783                       totaladc2[kcl] = z1*z1;
784                       ax[kcl]        = x1 * z1;
785                       ay[kcl]        = y1 * z1;
786                       ax2[kcl]       = 0.;
787                       ay2[kcl]       = 0.;
788                       status[j]      = 1;
789                     }
790                 }
791             }
792           
793           for(j = 0; j <= ncl[i]; j++)
794             {
795               Int_t   maxweight = 0;
796               Double_t     max  = 0.;
797               
798               if(status[j] == 0)
799                 { 
800                   x1 = x[iord[j]]; 
801                   y1 = y[iord[j]];
802                   z1 = z[iord[j]];
803                   t1 = t[iord[j]];
804                   sumweight = 0.;
805
806                   for(Int_t kcl = 0; kcl <= ig; kcl++)
807                     {
808                       x2 = xc[kcl]; 
809                       y2 = yc[kcl]; 
810                       rr = Distance(x1,y1,x2,y2);
811                       gweight     = exp(-(rr*rr)/(2*(1.2*1.2)));
812                       weight[kcl] = zc[kcl] * gweight;
813                       sumweight   = sumweight + weight[kcl];
814                       
815                       if(weight[kcl] > max)
816                         {
817                           max       =  weight[kcl];
818                           maxweight =  kcl;
819                         }
820                     }
821                   
822                   cellXY[cellCount[maxweight]][maxweight] = iord[j];
823                                   
824                   cellCount[maxweight]++;
825                   
826                   x2 = xc[maxweight];
827                   y2 = yc[maxweight];
828                   totaladc[maxweight]  +=  z1;
829                   ax[maxweight]        +=  x1*z1;
830                   ay[maxweight]        +=  y1*z1;
831                   totaladc2[maxweight] +=  z1*z1;
832                   ax2[maxweight]       +=  z1*(x1-x2)*(x1-x2);
833                   ay2[maxweight]       +=  z1*(y1-y2)*(y1-y2);
834                   ncell[maxweight]++;
835
836                 }
837             }
838           
839           for(Int_t kcl = 0; kcl <= ig; kcl++)
840             {
841               
842               if(totaladc[kcl] > 0.)
843                 {
844                   xclust[kcl] = (ax[kcl])/ totaladc[kcl];
845                   yclust[kcl] = (ay[kcl])/ totaladc[kcl];
846                   
847                   //natasha
848                   Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
849                   if(totaladc2[kcl] >= sqtotadc)
850                     {
851                       sigxclust[kcl] = 0.25;
852                       sigyclust[kcl] = 0.25;
853                     }
854                   else
855                     {
856                       sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
857                       sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
858                     }   
859                 }
860               
861               for(j = 0; j < cellCount[kcl]; j++) clno++; 
862               
863               if (clno >= 4608) 
864                 {
865                   AliWarning("RefClust: Too many clusters! more than 4608");
866                   return;
867                 }
868               clusdata[0] = xclust[kcl];
869               clusdata[1] = yclust[kcl];
870               clusdata[2] = totaladc[kcl];
871               clusdata[3] = ncell[kcl];
872               
873               
874               if(sigxclust[kcl] > sigyclust[kcl]) 
875                 {
876                   clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
877                   clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
878                 }
879               else
880                 {
881                   clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
882                   clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
883                 }
884               
885               clxy[0] = tc[kcl];
886               
887               Int_t Ncell=1;
888               for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
889                 {
890                   if(ii<18) 
891                     {   
892                       clxy[Ncell] = t[cellXY[ii][kcl]];
893                       Ncell++;
894                     }
895                 } 
896               
897               pmdcludata = new AliPMDcludata(clusdata,clxy);
898               fPMDclucont->Add(pmdcludata);
899             }
900           
901           delete [] iord;
902           delete [] tc;   
903           delete [] t;
904           delete [] x;
905           delete [] y;
906           delete [] z;
907           delete [] xc;
908           delete [] yc;
909           delete [] zc;
910           
911           delete [] cellCount;
912           for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
913           
914           delete [] status;
915           delete [] totaladc;
916           delete [] totaladc2;
917           delete [] ncell;
918           delete [] xclust;
919           delete [] yclust;
920           delete [] sigxclust;
921           delete [] sigyclust;
922           delete [] ax;
923           delete [] ay;
924           delete [] ax2;
925           delete [] ay2;
926           delete [] weight;
927         }
928     }
929   delete [] ncl;
930   delete [] clxy;
931 }
932 // ------------------------------------------------------------------------ //
933 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1, 
934                                       Double_t x2, Double_t y2)
935 {
936   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
937 }
938 // ------------------------------------------------------------------------ //
939 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
940 {
941   fCutoff = decut;
942 }
943 // ------------------------------------------------------------------------ //
944 void AliPMDClusteringV1::SetClusteringParam(Int_t cluspar)
945 {
946   fClusParam = cluspar;
947 }
948 // ------------------------------------------------------------------------ //