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