defects fixed
[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 = 0,  j = 0, nmx1 = 0;
115   Int_t    incr = 0, id = 0, jd = 0;
116   Int_t    celldataX[kNmaxCell], celldataY[kNmaxCell];
117   Int_t    celldataTr[kNmaxCell], celldataPid[kNmaxCell];
118   Float_t  celldataAdc[kNmaxCell];
119   Float_t  clusdata[6] = {0.,0.,0.,0.,0.,0.};
120   Double_t cutoff, ave;
121   Double_t edepcell[kNMX];
122   Double_t cellenergy[kNMX];
123   
124   // ndimXr and ndimYr are different because of different module size
125
126   Int_t ndimXr = 0;
127   Int_t ndimYr = 0;
128
129   if (ismn < 12)
130     {
131       ndimXr = 96;
132       ndimYr = 48;
133     }
134   else if (ismn >= 12 && ismn <= 23)
135     {
136       ndimXr = 48;
137       ndimYr = 96;
138     }
139
140   for (i = 0; i < kNMX; i++)
141   {
142       cellenergy[i] = 0.;
143   }
144
145   Int_t kk = 0;
146   for (i = 0; i < kNDIMX; i++)
147     {
148       for (j = 0; j < kNDIMY; j++)
149         {
150           edepcell[kk] = 0.;
151           kk++;
152         }
153     }
154
155   for (id = 0; id < ndimXr; id++)
156     {
157       for (jd = 0; jd < ndimYr; jd++)
158         {
159           j = jd;
160           i = id+(ndimYr/2-1)-(jd/2);
161
162           Int_t ij = i + j*kNDIMX;
163           
164           if (ismn < 12)
165             {
166               cellenergy[ij]    = celladc[jd][id];
167             }
168           else if (ismn >= 12 && ismn <= 23)
169             {
170               cellenergy[ij]    = celladc[id][jd];
171             }
172         }
173     }
174   
175   for (i = 0; i < kNMX; i++)
176     {
177       edepcell[i] = cellenergy[i];
178     }
179
180   Bool_t jsort = true;
181   Int_t iord1[kNMX];
182   TMath::Sort((Int_t)kNMX,edepcell,iord1,jsort);// 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=0,j=0,k=0,id1=0,id2=0,icl=0, numcell=0;
318   Int_t jd1=0,jd2=0, icell=0, cellcount=0;
319   Int_t clust[2][kndim];
320   static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
321
322   AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
323   
324   for (j = 0; j < kNDIMX; j++)
325     {
326       for(k = 0; k < kNDIMY; k++)
327         {
328           fInfocl[0][j][k] = 0;
329           fInfocl[1][j][k] = 0;
330         }
331     }
332   for(i=0; i < kNMX; i++)
333     {
334       fInfcl[0][i] = -1;
335       
336       j  = iord1[i];
337       id2 = j/kNDIMX;
338       id1 = j-id2*kNDIMX;
339
340       if(edepcell[j] <= cutoff)
341         {
342           fInfocl[0][id1][id2] = -1;
343         }
344     }
345
346   // ---------------------------------------------------------------
347   // crude clustering begins. Start with cell having largest adc
348   // count and loop over the cells in descending order of adc count
349   // ---------------------------------------------------------------
350
351   icl       = -1;
352   cellcount = -1;
353
354   for(icell = 0; icell <= nmx1; icell++)
355     {
356       j  = iord1[icell];
357       id2 = j/kNDIMX;
358       id1 = j-id2*kNDIMX;
359
360       if(fInfocl[0][id1][id2] == 0 )
361         {
362           icl++;
363           numcell = 0;
364           cellcount++; 
365           fInfocl[0][id1][id2] = 1;
366           fInfocl[1][id1][id2] = icl;
367           fInfcl[0][cellcount] = icl;
368           fInfcl[1][cellcount] = id1;
369           fInfcl[2][cellcount] = id2;
370
371           clust[0][numcell] = id1;
372           clust[1][numcell] = id2;
373           
374           for(i = 1; i < kndim; i++)
375             {
376               clust[0][i]=0;
377             }
378           // ---------------------------------------------------------------
379           // check for adc count in neib. cells. If ne 0 put it in this clust
380           // ---------------------------------------------------------------
381           for(i = 0; i < 6; i++)
382             {
383               jd1 = id1 + neibx[i];
384               jd2 = id2 + neiby[i];
385               if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
386                   fInfocl[0][jd1][jd2] == 0)
387                 {
388                   numcell++;
389                   fInfocl[0][jd1][jd2] = 2;
390                   fInfocl[1][jd1][jd2] = icl;
391                   clust[0][numcell]    = jd1;
392                   clust[1][numcell]    = jd2;
393                   cellcount++;
394                   fInfcl[0][cellcount] = icl;
395                   fInfcl[1][cellcount] = jd1;
396                   fInfcl[2][cellcount] = jd2;
397                 }
398             }
399           // ---------------------------------------------------------------
400           // check adc count for neighbour's neighbours recursively and
401           // if nonzero, add these to the cluster.
402           // ---------------------------------------------------------------
403           for(i = 1; i < kndim;i++)
404             {
405               if(clust[0][i] != 0)
406                 {
407                   id1 = clust[0][i];
408                   id2 = clust[1][i];
409                   for(j = 0; j < 6 ; j++)
410                     {
411                       jd1 = id1 + neibx[j];
412                       jd2 = id2 + neiby[j];
413                       if( (jd1 >= 0 && jd1 < kNDIMX) && 
414                           (jd2 >= 0 && jd2 < kNDIMY) &&
415                           fInfocl[0][jd1][jd2] == 0 )
416                         {
417                           fInfocl[0][jd1][jd2] = 2;
418                           fInfocl[1][jd1][jd2] = icl;
419                           numcell++;
420                           clust[0][numcell]    = jd1;
421                           clust[1][numcell]    = jd2;
422                           cellcount++;
423                           fInfcl[0][cellcount] = icl;
424                           fInfcl[1][cellcount] = jd1;
425                           fInfcl[2][cellcount] = jd2;
426                         }
427                     }
428                 }
429             }
430         }
431     }
432   return cellcount;
433 }
434 // ------------------------------------------------------------------------ //
435 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
436 {
437   // Does the refining of clusters
438   // Takes the big patch and does gaussian fitting and
439   // finds out the more refined clusters
440   //
441   
442   AliPMDcludata *pmdcludata = 0;
443   
444   const Int_t kNmaxCell   = 19;    // # of cells surrounding a cluster center
445   
446   Int_t ndim = incr + 1;
447   
448   Int_t    *ncl  = 0x0;
449   Int_t    *clxy = 0x0;  
450   Int_t    i12 = 0, i22 = 0;
451   Int_t    i = 0, j = 0, k = 0;
452   Int_t    i1 = 0, i2 = 0, id = 0, icl = 0;
453   Int_t    itest = 0, ihld = 0, ig = 0;
454   Int_t    nsupcl = 0, clno = 0, t1 = 0, t2 = 0;
455   Float_t  clusdata[6];
456   Double_t x1 = 0, y1 = 0, z1 = 0, x2 = 0, y2 = 0, z2 = 0, rr = 0;
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     }
466   for(i = 0; i<6; i++)
467     {
468       clusdata[i] = 0.;
469     }
470   for(i = 0; i<19; i++)
471     {
472       clxy[i] = 0;
473     }
474
475   // clno counts the final clusters
476   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
477   // x, y and z store (x,y) coordinates of and energy deposited in a cell
478   // xc, yc store (x,y) coordinates of the cluster center
479   // zc stores the energy deposited in a cluster
480   // rc is cluster radius
481   
482   clno  = -1;
483   nsupcl = -1;
484
485   for(i = 0; i <= incr; i++)
486     {
487       if(fInfcl[0][i] != nsupcl)
488         {
489           nsupcl++;
490         }
491       if (nsupcl > ndim) 
492         {
493           AliWarning("RefClust: Too many superclusters!");
494           nsupcl = ndim;
495           break;
496         }
497       ncl[nsupcl]++;
498     }
499   
500   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
501                   incr+1,nsupcl+1));
502   id  = -1;
503   icl = -1;
504   
505   for(i = 0; i <= nsupcl; i++) 
506     {
507       if(ncl[i] == 0)
508         {
509           id++;
510           icl++;
511           if (clno >= 4608) 
512             {
513               AliWarning("RefClust: Too many clusters! more than 4608");
514               return;
515             }
516           clno++;
517           i1 = fInfcl[1][id];
518           i2 = fInfcl[2][id];
519           
520           i12 = i1 + i2*kNDIMX;
521           
522           clusdata[0] = fCoord[0][i1][i2];
523           clusdata[1] = fCoord[1][i1][i2];
524           clusdata[2] = edepcell[i12];
525           clusdata[3] = 1.;
526           clusdata[4] = 999.5;
527           clusdata[5] = 999.5;
528
529           clxy[0] = i1*10000 + i2;
530           
531           for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
532             {
533               clxy[icltr] = -1;
534             }
535
536           pmdcludata  = new AliPMDcludata(clusdata,clxy);
537           fPMDclucont->Add(pmdcludata);
538         }
539       else if(ncl[i] == 1) 
540         {
541           id++;
542           icl++;
543           if (clno >= 4608) 
544             {
545               AliWarning("RefClust: Too many clusters! more than 4608");
546               return;
547             }
548           clno++;
549           i1   = fInfcl[1][id];
550           i2   = fInfcl[2][id];
551           i12  = i1 + i2*kNDIMX;
552
553           x1   = fCoord[0][i1][i2];
554           y1   = fCoord[1][i1][i2];
555           z1   = edepcell[i12];
556
557           clxy[0] = i1*10000 + i2;
558           
559           id++;
560           i1   = fInfcl[1][id];
561           i2   = fInfcl[2][id];
562
563           i22  = i1 + i2*kNDIMX;
564           x2   = fCoord[0][i1][i2];
565           y2   = fCoord[1][i1][i2];
566           z2   = edepcell[i22];
567
568           clxy[1] = i1*10000 + i2;
569           
570
571           for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
572             {
573               clxy[icltr] = -1;
574             }
575           
576           clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
577           clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
578           clusdata[2] = z1+z2;
579           clusdata[3] = 2.;
580           clusdata[4] = 0.5;
581           clusdata[5] = 0.0;
582           pmdcludata  = new AliPMDcludata(clusdata,clxy);
583           fPMDclucont->Add(pmdcludata);
584         }
585       else
586         {
587           
588           Int_t    *iord, *tc, *t;
589           Double_t *x, *y, *z, *xc, *yc, *zc;
590
591           iord = new Int_t [ncl[i]+1];
592           tc   = new Int_t [ncl[i]+1];
593           t    = new Int_t [ncl[i]+1];
594           
595           x    = new Double_t [ncl[i]+1];
596           y    = new Double_t [ncl[i]+1];
597           z    = new Double_t [ncl[i]+1];
598           xc   = new Double_t [ncl[i]+1];
599           yc   = new Double_t [ncl[i]+1];
600           zc   = new Double_t [ncl[i]+1];
601           
602           for( k = 0; k < ncl[i]+1; k++)
603             {
604               iord[k] = -1;
605               t[k]    = -1;
606               tc[k]   = -1;
607               x[k]    = -1;
608               y[k]    = -1;
609               z[k]    = -1;
610               xc[k]   = -1;
611               yc[k]   = -1;
612               zc[k]   = -1;
613             }
614           id++;
615           // super-cluster of more than two cells - broken up into smaller
616           // clusters gaussian centers computed. (peaks separated by > 1 cell)
617           // Begin from cell having largest energy deposited This is first
618           // cluster center
619           i1      = fInfcl[1][id];
620           i2      = fInfcl[2][id];
621           i12     = i1 + i2*kNDIMX;
622           
623           x[0]    = fCoord[0][i1][i2];
624           y[0]    = fCoord[1][i1][i2];
625           z[0]    = edepcell[i12];
626           t[0]    = i1*10000 + i2;
627           
628
629           iord[0] = 0;
630           for(j = 1; j <= ncl[i]; j++)
631             {
632               id++;
633               i1      = fInfcl[1][id];
634               i2      = fInfcl[2][id];
635               i12     = i1 + i2*kNDIMX;
636
637               iord[j] = j;
638               x[j]    = fCoord[0][i1][i2];
639               y[j]    = fCoord[1][i1][i2];
640               z[j]    = edepcell[i12];
641               t[j]    = i1*10000 + i2;
642
643             }
644           
645           // arranging cells within supercluster in decreasing order
646           
647           for(j = 1;j <= ncl[i]; j++)
648             {
649               itest = 0;
650               ihld  = iord[j];
651               for(i1 = 0; i1 < j; i1++)
652                 {
653                   if(itest == 0 && z[iord[i1]] < z[ihld])
654                     {
655                       itest = 1;
656                       for(i2 = j-1; i2 >= i1; i2--)
657                         {
658                           iord[i2+1] = iord[i2];
659                         }
660                       iord[i1] = ihld;
661                     }
662                 }
663             }
664
665           if (fClusParam == 1)
666             {
667               // Clustering algorithm returns from here for PP collisions
668               // for pp, only the output of crude clusterng is taken
669               // sigx and sigy are not calculated at this moment
670
671               Double_t supx=0., supy=0., supz=0.;
672           
673               for(j = 0;j <= ncl[i]; j++)
674                 {
675                   supx += x[iord[j]]*z[iord[j]];
676                   supy += y[iord[j]]*z[iord[j]];
677                   supz += z[iord[j]];
678                   if(j < 19)
679                     {
680                       clxy[j] = t[iord[j]];
681                     }
682                 }
683               
684               if( ncl[i] + 1 < 19)
685                 {
686                   for(Int_t ncel = ncl[i] + 1; ncel < kNmaxCell; ncel ++ )
687                     {
688                       clxy[ncel] = -1;
689                     }
690                 }
691               clusdata[0] = supx/supz;
692               clusdata[1] = supy/supz;
693               clusdata[2] = supz;
694               clusdata[3] = ncl[i]+1;
695               clusdata[4] = 0.5;
696               clusdata[5] = 0.0;
697               pmdcludata  = new AliPMDcludata(clusdata,clxy);
698               fPMDclucont->Add(pmdcludata);
699             }
700
701           /* MODIFICATION PART STARTS (Tapan July 2008)
702              iord[0] is the cell with highest ADC in the crude-cluster
703              ig is the number of local maxima in the crude-cluster
704              For the higest peak we make ig=0 which means first local maximum.
705              Next we go down in terms of the ADC sequence and find out if any
706              more of the cells form local maxima. The definition of local 
707              maxima is that all its neighbours are of less ADC compared to it. 
708           */
709
710           if (fClusParam == 2)
711             {
712               // This part is to split the supercluster
713               //
714               ig = 0;
715               xc[ig] = x[iord[0]];
716               yc[ig] = y[iord[0]];
717               zc[ig] = z[iord[0]];
718               tc[ig] = t[iord[0]];
719               Int_t ivalid = 0,  icount = 0;
720           
721               for(j=1;j<=ncl[i];j++)
722                 {
723                   x1  = x[iord[j]];
724                   y1  = y[iord[j]]; 
725                   z1  = z[iord[j]];
726                   t1  = t[iord[j]];
727                   rr=Distance(x1,y1,xc[ig],yc[ig]);
728                   
729                   // Check the cells which are outside the neighbours (rr>1.2)
730                   if(rr>1.2 ) 
731                     {
732                       ivalid=0;
733                       icount=0;
734                       for(Int_t j1=1;j1<j;j1++)
735                         {
736                           icount++;
737                           Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
738                           if(rr1>1.2) ivalid++;
739                         }
740                       if(ivalid == icount && z1>0.5*zc[ig])
741                         {
742                           ig++;
743                           xc[ig]=x1; 
744                           yc[ig]=y1; 
745                           zc[ig]=z1;
746                           tc[ig]=t1;
747                         }
748                     }     
749                 }
750           
751               icl=icl+ig+1;
752           
753               //  We use simple Gaussian weighting. (Tapan Jan 2005)
754               // compute the number of cells belonging to each cluster.
755               // cell can be shared between several clusters  
756               // in the ratio of cluster energy deposition
757               // To calculate: 
758               // (1) number of cells belonging to a cluster (ig) and 
759               // (2) total ADC of the cluster (ig) 
760               // (3) x and y positions of the cluster
761           
762               
763               Int_t *cellCount;
764               Int_t **cellXY;
765               
766               Int_t    *status;
767               Double_t *totaladc, *totaladc2, *ncell,*weight;
768               Double_t *xclust, *yclust, *sigxclust, *sigyclust;
769               Double_t *ax, *ay, *ax2, *ay2;
770           
771           
772               status    = new Int_t [ncl[i]+1];
773               cellXY    = new Int_t *[ncl[i]+1];
774           
775               cellCount = new Int_t [ig+1];
776               totaladc  = new Double_t [ig+1];
777               totaladc2 = new Double_t [ig+1];
778               ncell     = new Double_t [ig+1];
779               weight    = new Double_t [ig+1];
780               xclust    = new Double_t [ig+1];
781               yclust    = new Double_t [ig+1];
782               sigxclust = new Double_t [ig+1];
783               sigyclust = new Double_t [ig+1];
784               ax        = new Double_t [ig+1];
785               ay        = new Double_t [ig+1];
786               ax2       = new Double_t [ig+1];
787               ay2       = new Double_t [ig+1];
788               
789               for(j = 0; j < ncl[i]+1; j++) 
790                 {
791                   status[j]     = 0;
792                   cellXY[j] = new Int_t[ig+1];
793                 }
794               //initialization
795               for(Int_t kcl = 0; kcl < ig+1; kcl++)
796                 {
797                   cellCount[kcl] = 0;
798                   totaladc[kcl]  = 0.;
799                   totaladc2[kcl] = 0.;
800                   ncell[kcl]     = 0.;
801                   weight[kcl]    = 0.;  
802                   xclust[kcl]    = 0.; 
803                   yclust[kcl]    = 0.;
804                   sigxclust[kcl] = 0.; 
805                   sigyclust[kcl] = 0.;
806                   ax[kcl]        = 0.;      
807                   ay[kcl]        = 0.;      
808                   ax2[kcl]       = 0.;      
809                   ay2[kcl]       = 0.;    
810                   for(j = 0; j < ncl[i]+1; j++)
811                     {
812                       cellXY[j][kcl] = 0;
813                     }
814                 }
815               Double_t sumweight, gweight; 
816               
817               for(j = 0;j <= ncl[i]; j++)
818                 {
819                   x1 = x[iord[j]];
820                   y1 = y[iord[j]];
821                   z1 = z[iord[j]];
822                   t1 = t[iord[j]];
823                   
824                   for(Int_t kcl=0; kcl<=ig; kcl++)
825                     {
826                       x2 = xc[kcl];
827                       y2 = yc[kcl];
828                       rr = Distance(x1,y1,x2,y2);
829                       t2 = tc[kcl];               
830                       
831                       if(rr==0)
832                         {
833                           ncell[kcl]     = 1.;
834                           totaladc[kcl]  = z1;
835                           totaladc2[kcl] = z1*z1;
836                           ax[kcl]        = x1 * z1;
837                           ay[kcl]        = y1 * z1;
838                           ax2[kcl]       = 0.;
839                           ay2[kcl]       = 0.;
840                           status[j]      = 1;
841                         }
842                     }
843                 }
844           
845               for(j = 0; j <= ncl[i]; j++)
846                 {
847                   Int_t   maxweight = 0;
848                   Double_t     max  = 0.;
849                   
850                   if(status[j] == 0)
851                     { 
852                       x1 = x[iord[j]]; 
853                       y1 = y[iord[j]];
854                       z1 = z[iord[j]];
855                       t1 = t[iord[j]];
856                       sumweight = 0.;
857
858                       for(Int_t kcl = 0; kcl <= ig; kcl++)
859                         {
860                           x2 = xc[kcl]; 
861                           y2 = yc[kcl]; 
862                           rr = Distance(x1,y1,x2,y2);
863                           gweight     = exp(-(rr*rr)/(2*(1.2*1.2)));
864                           weight[kcl] = zc[kcl] * gweight;
865                           sumweight   = sumweight + weight[kcl];
866                           
867                           if(weight[kcl] > max)
868                             {
869                               max       =  weight[kcl];
870                               maxweight =  kcl;
871                             }
872                         }
873                       
874                       cellXY[cellCount[maxweight]][maxweight] = iord[j];
875                       
876                       cellCount[maxweight]++;
877                       
878                       x2 = xc[maxweight];
879                       y2 = yc[maxweight];
880                       totaladc[maxweight]  +=  z1;
881                       ax[maxweight]        +=  x1*z1;
882                       ay[maxweight]        +=  y1*z1;
883                       totaladc2[maxweight] +=  z1*z1;
884                       ax2[maxweight]       +=  z1*(x1-x2)*(x1-x2);
885                       ay2[maxweight]       +=  z1*(y1-y2)*(y1-y2);
886                       ncell[maxweight]++;
887                       
888                     }
889                 }
890           
891               for(Int_t kcl = 0; kcl <= ig; kcl++)
892                 {
893                   if(totaladc[kcl] > 0.)
894                     {
895                       xclust[kcl] = (ax[kcl])/ totaladc[kcl];
896                       yclust[kcl] = (ay[kcl])/ totaladc[kcl];
897                       
898                       //natasha
899                       Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
900                       if(totaladc2[kcl] >= sqtotadc)
901                         {
902                           sigxclust[kcl] = 0.25;
903                           sigyclust[kcl] = 0.25;
904                         }
905                       else
906                         {
907                           sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
908                           sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
909                         }       
910                     }
911                   
912                   for(j = 0; j < cellCount[kcl]; j++) clno++; 
913                   
914                   if (clno >= 4608) 
915                     {
916                       AliWarning("RefClust: Too many clusters! more than 4608");
917                       return;
918                     }
919                   clusdata[0] = xclust[kcl];
920                   clusdata[1] = yclust[kcl];
921                   clusdata[2] = totaladc[kcl];
922                   clusdata[3] = ncell[kcl];
923                   
924                   if(sigxclust[kcl] > sigyclust[kcl]) 
925                     {
926                       clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
927                       clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
928                     }
929                   else
930                     {
931                       clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
932                       clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
933                     }
934                   
935                   clxy[0] = tc[kcl];
936                   
937                   Int_t Ncell=1;
938                   for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
939                     {
940                       if(ii<18) 
941                         {       
942                           clxy[Ncell] = t[cellXY[ii][kcl]];
943                           Ncell++;
944                         }
945                     } 
946                   
947                   pmdcludata = new AliPMDcludata(clusdata,clxy);
948                   fPMDclucont->Add(pmdcludata);
949                 }
950               delete [] cellCount;
951               for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
952               
953               delete [] status;
954               delete [] totaladc;
955               delete [] totaladc2;
956               delete [] ncell;
957               delete [] xclust;
958               delete [] yclust;
959               delete [] sigxclust;
960               delete [] sigyclust;
961               delete [] ax;
962               delete [] ay;
963               delete [] ax2;
964               delete [] ay2;
965               delete [] weight;
966
967             }
968
969           delete [] iord;
970           delete [] tc;   
971           delete [] t;
972           delete [] x;
973           delete [] y;
974           delete [] z;
975           delete [] xc;
976           delete [] yc;
977           delete [] zc;
978
979
980         }
981     }
982   delete [] ncl;
983   delete [] clxy;
984 }
985 // ------------------------------------------------------------------------ //
986 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1, 
987                                       Double_t x2, Double_t y2)
988 {
989   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
990 }
991 // ------------------------------------------------------------------------ //
992 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
993 {
994   fCutoff = decut;
995 }
996 // ------------------------------------------------------------------------ //
997 void AliPMDClusteringV1::SetClusteringParam(Int_t cluspar)
998 {
999   fClusParam = cluspar;
1000 }
1001 // ------------------------------------------------------------------------ //