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