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