fixed allocation of memory converted to dynamic
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV2.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 : PMDClusteringV2.cxx                  //
21 //                                                     //
22 //  clustering code for alice pmd                      //
23 //                                                     //
24 //-----------------------------------------------------//
25
26 /* --------------------------------------------------------------------
27    Code developed by S. C. Phatak, Institute of Physics,
28    Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
29    ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
30    builds up superclusters and breaks them into clusters. The input is
31    in TObjarray  and cluster information is in TObjArray.
32    integer clno gives total number of clusters in the  supermodule.
33    fClusters is the  global ( public ) variables.
34    Others are local ( private ) to the code.
35    At the moment, the data is read for whole detector ( all supermodules
36    and pmd as well as cpv. This will have to be modify later )
37    LAST UPDATE  :  October 23, 2002
38 -----------------------------------------------------------------------*/
39
40 #include <Riostream.h>
41 #include <TMath.h>
42 #include <TObjArray.h>
43 #include <TArrayI.h>
44
45 #include "AliPMDcludata.h"
46 #include "AliPMDcluster.h"
47 #include "AliPMDClustering.h"
48 #include "AliPMDClusteringV2.h"
49 #include "AliLog.h"
50
51 ClassImp(AliPMDClusteringV2)
52
53 const Double_t AliPMDClusteringV2::fgkSqroot3by2=0.8660254;  // sqrt(3.)/2.
54
55 AliPMDClusteringV2::AliPMDClusteringV2():
56   fPMDclucont(new TObjArray()),
57   fCutoff(0.0)
58 {
59   for(int i = 0; i < kNDIMX; i++)
60     {
61       for(int j = 0; j < kNDIMY; j++)
62         {
63           fCoord[0][i][j] = i+j/2.;
64           fCoord[1][i][j] = fgkSqroot3by2*j;
65         }
66     }
67 }
68 // ------------------------------------------------------------------------ //
69
70
71 AliPMDClusteringV2::AliPMDClusteringV2(const AliPMDClusteringV2& pmdclv2):
72   AliPMDClustering(pmdclv2),
73   fPMDclucont(0),
74   fCutoff(0)
75 {
76   // copy constructor
77   AliError("Copy constructor not allowed ");
78   
79 }
80 // ------------------------------------------------------------------------ //
81 AliPMDClusteringV2 &AliPMDClusteringV2::operator=(const AliPMDClusteringV2& /*pmdclv2*/)
82 {
83   // copy constructor
84   AliError("Assignment operator not allowed ");
85   return *this;
86 }
87 // ------------------------------------------------------------------------ //
88 AliPMDClusteringV2::~AliPMDClusteringV2()
89 {
90   delete fPMDclucont;
91 }
92 // ------------------------------------------------------------------------ //
93
94 void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn, 
95                                  Double_t celladc[48][96], TObjArray *pmdcont)
96 {
97   // main function to call other necessary functions to do clustering
98   //
99   AliPMDcluster *pmdcl = 0;
100
101   const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
102
103   Int_t    i, j, nmx1, incr, id, jd;
104   Int_t    ndimXr = 0;
105   Int_t    ndimYr = 0;
106   Int_t    celldataX[15], celldataY[15];
107   Float_t  clusdata[6];  
108   Double_t cutoff, ave;
109   Double_t edepcell[kNMX];
110
111
112   if (ismn < 12)
113     {
114       ndimXr = 96;
115       ndimYr = 48;
116     }
117   else if (ismn >= 12 && ismn <= 23)
118     {
119       ndimXr = 48;
120       ndimYr = 96;
121     }
122   
123   for (Int_t i =0; i < kNMX; i++)
124     {
125      edepcell[i] = 0.;
126     }
127     
128   for (id = 0; id < ndimXr; id++)
129     {
130       for (jd = 0; jd < ndimYr; jd++)
131         {
132           j = jd;
133           i = id + (ndimYr/2-1) - (jd/2);
134           Int_t ij = i + j*kNDIMX;
135           if (ismn < 12)
136             {
137               edepcell[ij]    = celladc[jd][id];
138             }
139           else if (ismn >= 12 && ismn <= 23)
140             {
141              edepcell[ij]    = celladc[id][jd];
142             }
143
144         }
145     }
146
147   Int_t iord1[kNMX];
148   TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
149   cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
150   ave  = 0.;
151   nmx1 = -1;
152
153   for(i = 0;i < kNMX; i++)
154     {
155       if(edepcell[i] > 0.) 
156         {
157           ave += edepcell[i];
158         }
159       if(edepcell[i] > cutoff )
160         {
161           nmx1++;
162         }
163     }
164   
165   AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
166   
167   if (nmx1 == 0) 
168     {
169       nmx1 = 1;
170     }
171   ave = ave/nmx1;
172   
173   AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
174                   kNMX,ave));
175   
176   incr = CrClust(ave, cutoff, nmx1,iord1, edepcell);
177   RefClust(incr,edepcell );
178   
179   Int_t nentries1 = fPMDclucont->GetEntries();
180   AliDebug(1,Form("Detector Plane = %d  Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
181   AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
182   for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
183     {
184       AliPMDcludata *pmdcludata = 
185         (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
186       Float_t cluXC    = pmdcludata->GetClusX();
187       Float_t cluYC    = pmdcludata->GetClusY();
188       Float_t cluADC   = pmdcludata->GetClusADC();
189       Float_t cluCELLS = pmdcludata->GetClusCells();
190       Float_t cluSIGX  = pmdcludata->GetClusSigmaX();
191       Float_t cluSIGY  = pmdcludata->GetClusSigmaY();
192       
193       Float_t cluY0    = ktwobysqrt3*cluYC;
194       Float_t cluX0    = cluXC - cluY0/2.;
195       
196       // 
197       // Cluster X centroid is back transformed
198       //
199       if (ismn < 12)
200         {
201           clusdata[0] = cluX0 - (24-1) + cluY0/2.;
202         }
203       else if (ismn  >= 12 && ismn <= 23)
204         {
205           clusdata[0] = cluX0 - (48-1) + cluY0/2.;
206         }         
207
208       clusdata[1]     = cluY0;
209       clusdata[2]     = cluADC;
210       clusdata[3]     = cluCELLS;
211       clusdata[4]     = cluSIGX;
212       clusdata[5]     = cluSIGY;
213       //
214       // Cells associated with a cluster
215       //
216       for (Int_t ihit = 0; ihit < 15; ihit++)
217         {
218           Int_t dummyXY = pmdcludata->GetCellXY(ihit);
219          
220           Int_t celldumY   = dummyXY%10000;
221           Int_t celldumX   = dummyXY/10000;
222           //Float_t cellY    = (Float_t) (ktwobysqrt3*celldumY/10);
223           //Float_t cellX    = (Float_t) (celldumX/10 - (celldumY/2.)/10);
224            
225           Float_t cellY    = (Float_t) celldumY/10;
226           Float_t cellX    = (Float_t) celldumX/10;
227
228           
229           // 
230           // Cell X centroid is back transformed
231           //
232           if (ismn < 12)
233             {
234               celldataX[ihit] = (Int_t) ((cellX - (24-1) + cellY/2.) + 0.5);
235             }
236           else if (ismn  >= 12 && ismn <= 23)
237             {
238               celldataX[ihit] = (Int_t) ((cellX - (48-1) + cellY/2.) + 0.5 );
239             }     
240           celldataY[ihit] = (Int_t) (cellY + 0.5);
241         }
242
243       pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
244       pmdcont->Add(pmdcl);
245     }
246   fPMDclucont->Clear();
247 }
248 // ------------------------------------------------------------------------ //
249 Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
250                                   Int_t iord1[], Double_t edepcell[])
251 {
252   // Does crude clustering
253   // Finds out only the big patch by just searching the
254   // connected cells
255   //
256
257   Int_t i,j,k,id1,id2,icl, numcell;
258   Int_t jd1,jd2, icell, cellcount;
259   Int_t clust[2][5000];
260   static Int_t neibx[6] = {1,0,-1,-1,0,1}, neiby[6] = {0,1,1,0,-1,-1};
261
262   // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
263   // cell. There are six neighbours.
264   // cellcount --- total number of cells having nonzero ener dep
265   // numcell --- number of cells in a given supercluster
266   
267   AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
268   
269   for (j=0; j < kNDIMX; j++)
270     {
271       for(k=0; k < kNDIMY; k++)
272         {
273           fInfocl[0][j][k] = 0;
274           fInfocl[1][j][k] = 0;
275         }
276     }
277  
278   for(i=0; i < kNMX; i++)
279     {
280       fInfcl[0][i] = -1;
281       
282       j  = iord1[i];
283       id2 = j/kNDIMX;
284       id1 = j-id2*kNDIMX;
285       
286       if(edepcell[j] <= cutoff)
287         {
288           fInfocl[0][id1][id2] = -1;
289         }
290     }
291   // ---------------------------------------------------------------
292   // crude clustering begins. Start with cell having largest adc
293   // count and loop over the cells in descending order of adc count
294   // ---------------------------------------------------------------
295   icl       = -1;
296   cellcount = -1;
297   for(icell=0; icell <= nmx1; icell++)
298     {
299       j  = iord1[icell];
300       id2 = j/kNDIMX;
301       id1 = j-id2*kNDIMX;
302       if(fInfocl[0][id1][id2] == 0 )
303         {
304           // ---------------------------------------------------------------
305           // icl -- cluster #, numcell -- # of cells in it, clust -- stores
306           // coordinates of the cells in a cluster, fInfocl[0][i1][i2] is 1 for
307           // primary and 2 for secondary cells,
308           // fInfocl[1][i1][i2] stores cluster #
309           // ---------------------------------------------------------------
310           icl++;
311           numcell = 0;
312           cellcount++;
313           fInfocl[0][id1][id2]  = 1;
314           fInfocl[1][id1][id2]  = icl;
315           fInfcl[0][cellcount]  = icl;
316           fInfcl[1][cellcount]  = id1;
317           fInfcl[2][cellcount]  = id2;
318
319           clust[0][numcell]     = id1;
320           clust[1][numcell]     = id2;
321           for(i = 1; i < 5000; i++)
322             {
323               clust[0][i] = -1;
324             }
325           // ---------------------------------------------------------------
326           // check for adc count in neib. cells. If ne 0 put it in this clust
327           // ---------------------------------------------------------------
328           for(i = 0; i < 6; i++)
329             {
330             jd1 = id1 + neibx[i];
331             jd2 = id2 + neiby[i];
332             if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
333                 fInfocl[0][jd1][jd2] == 0)
334               {
335                 numcell++;
336                 fInfocl[0][jd1][jd2] = 2;
337                 fInfocl[1][jd1][jd2] = icl;
338                 clust[0][numcell]    = jd1;
339                 clust[1][numcell]    = jd2;
340                 cellcount++;
341                 fInfcl[0][cellcount] = icl;
342                 fInfcl[1][cellcount] = jd1;
343                 fInfcl[2][cellcount] = jd2;
344               }
345             }
346           // ---------------------------------------------------------------
347           // check adc count for neighbour's neighbours recursively and
348           // if nonzero, add these to the cluster.
349           // ---------------------------------------------------------------
350           for(i = 1;i < 5000; i++)
351             { 
352               if(clust[0][i] != -1)
353                 {
354                   id1 = clust[0][i];
355                   id2 = clust[1][i];
356                   for(j = 0; j < 6 ; j++)
357                     {
358                       jd1 = id1 + neibx[j];
359                       jd2 = id2 + neiby[j];
360                       if( (jd1 >= 0 && jd1 < kNDIMX) && 
361                           (jd2 >= 0 && jd2 < kNDIMY) 
362                           && fInfocl[0][jd1][jd2] == 0 )
363                         {
364                           fInfocl[0][jd1][jd2] = 2;
365                           fInfocl[1][jd1][jd2] = icl;
366                           numcell++;
367                           clust[0][numcell]    = jd1;
368                           clust[1][numcell]    = jd2;
369                           cellcount++;
370                           fInfcl[0][cellcount] = icl;
371                           fInfcl[1][cellcount] = jd1;
372                           fInfcl[2][cellcount] = jd2;
373                         }
374                     }
375                 }
376             }
377         }
378     }
379   return cellcount;
380 }
381 // ------------------------------------------------------------------------ //
382   void AliPMDClusteringV2::RefClust(Int_t incr, Double_t edepcell[])
383 {
384   // Does the refining of clusters
385   // Takes the big patch and does gaussian fitting and
386   // finds out the more refined clusters
387
388   const Float_t ktwobysqrt3 = 1.1547;
389   AliPMDcludata *pmdcludata = 0;
390
391
392   Int_t    i, j, k, i1, i2, id, icl, itest, ihld;
393   Int_t    ig, nsupcl, clno, clX,clY;
394   Int_t    clxy[15];
395
396   Float_t  clusdata[6];
397   Double_t x1, y1, z1, x2, y2, z2, rr;
398
399   Int_t kndim = incr + 1;
400
401   TArrayI testncl;
402   TArrayI testindex;
403
404   Int_t    *ncl, *iord;
405
406   Double_t *x, *y, *z, *xc, *yc, *zc, *cells, *rcl, *rcs;
407
408   ncl   = new Int_t [kndim];
409   iord  = new Int_t [kndim];
410   x     = new Double_t [kndim];
411   y     = new Double_t [kndim];
412   z     = new Double_t [kndim];
413   xc    = new Double_t [kndim];
414   yc    = new Double_t [kndim];
415   zc    = new Double_t [kndim];
416   cells = new Double_t [kndim];
417   rcl   = new Double_t [kndim];
418   rcs   = new Double_t [kndim];
419   
420   for(Int_t kk = 0; kk < 15; kk++)
421     {
422       if( kk < 6 )clusdata[kk] = 0.;
423     }
424    
425   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
426   // x, y and z store (x,y) coordinates of and energy deposited in a cell
427   // xc, yc store (x,y) coordinates of the cluster center
428   // zc stores the energy deposited in a cluster, rc is cluster radius
429
430   clno   = -1;
431   nsupcl = -1;
432
433   for(i = 0; i < kndim; i++)
434     {
435       ncl[i] = -1;
436     }
437   for(i = 0; i <= incr; i++)
438     {
439       if(fInfcl[0][i] != nsupcl)
440         {
441           nsupcl++;
442         }
443       if (nsupcl > 4500) 
444         {
445           AliWarning("RefClust: Too many superclusters!");
446           nsupcl = 4500;
447           break;
448         }
449       ncl[nsupcl]++;
450     }
451   
452   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
453                   incr+1,nsupcl+1));
454   
455   id  = -1;
456   icl = -1;
457   for(i = 0; i <= nsupcl; i++)
458     {
459       if(ncl[i] == 0)
460         {
461           id++;
462           icl++;
463           // one  cell super-clusters --> single cluster
464           // cluster center at the centyer of the cell
465           // cluster radius = half cell dimension
466           if (clno >= 5000) 
467             {
468               AliWarning("RefClust: Too many clusters! more than 5000");
469               return;
470             }
471           clno++;
472           i1          = fInfcl[1][id];
473           i2          = fInfcl[2][id];
474           Int_t i12 = i1 + i2*kNDIMX;
475           clusdata[0] = fCoord[0][i1][i2];
476           clusdata[1] = fCoord[1][i1][i2];
477           clusdata[2] = edepcell[i12];
478           clusdata[3] = 1.;
479           clusdata[4] = 0.0;
480           clusdata[5] = 0.0;
481           
482           //cell information
483           //Float_t cellY    = (Float_t) (ktwobysqrt3*celldumY/10);
484           //Float_t cellX    = (Float_t) (celldumX/10 - (celldumY/2.)/10);
485           
486           clY = (Int_t)((ktwobysqrt3*fCoord[1][i1][i2])*10);
487           clX = (Int_t)((fCoord[0][i1][i2] - clY/20.)*10);
488           clxy[0] = clX*10000 + clY ;
489
490           //clX = (Int_t) fCoord[0][i1][i2]*10;
491           //clY = (Int_t) fCoord[1][i1][i2]*10;
492           //clxy[0] = clX*10000 + clY;
493
494
495           for(Int_t icltr = 1; icltr < 15; icltr++)
496             {
497               clxy[icltr] = -1;
498             }
499           pmdcludata  = new AliPMDcludata(clusdata,clxy);
500           fPMDclucont->Add(pmdcludata);
501           
502           
503         }
504       else if(ncl[i] == 1)
505         {
506           // two cell super-cluster --> single cluster
507           // cluster center is at ener. dep.-weighted mean of two cells
508           // cluster radius == half cell dimension
509           id++;
510           icl++;
511           if (clno >= 5000) 
512             {
513               AliWarning("RefClust: Too many clusters! more than 5000");
514               return;
515             }
516           clno++;
517           i1   = fInfcl[1][id];
518           i2   = fInfcl[2][id];
519           Int_t i12 = i1 + i2*kNDIMX;
520           
521           x1   = fCoord[0][i1][i2];
522           y1   = fCoord[1][i1][i2];
523           z1   = edepcell[i12];
524           
525           id++;
526           i1   = fInfcl[1][id];
527           i2   = fInfcl[2][id];
528           i12  = i1 + i2*kNDIMX;
529           
530           x2   = fCoord[0][i1][i2];
531           y2   = fCoord[1][i1][i2];
532           z2   = edepcell[i12];
533           
534           clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
535           clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
536           clusdata[2] = z1+z2;
537           clusdata[3] = 2.;
538           clusdata[4] = (TMath::Sqrt(z1*z2))/(z1+z2);
539           clusdata[5] = 0.0;
540
541           clY = (Int_t)((ktwobysqrt3*y1)*10);
542           clX = (Int_t)((x1 - clY/20.)*10);
543           clxy[0] = clX*10000 + clY ;
544
545           //clX = (Int_t) x1*10;
546           //clY = (Int_t) y1*10;
547           //clxy[0] = clX*10000 + clY;
548           
549           clY = (Int_t)((ktwobysqrt3*y2)*10);
550           clX = (Int_t)((x2 - clY/20.)*10);
551           clxy[1] = clX*10000 + clY ;
552
553           //clX = (Int_t) x2*10;
554           //clY = (Int_t) y2*10;
555           //clxy[1] = clX*10000 + clY;
556           
557           for(Int_t icltr = 2; icltr < 15; icltr++)
558             {
559               clxy[icltr] = -1;
560             }
561           pmdcludata  = new AliPMDcludata(clusdata, clxy);
562           fPMDclucont->Add(pmdcludata);
563         }
564       else{
565         id++;
566         iord[0] = 0;
567         // super-cluster of more than two cells - broken up into smaller
568         // clusters gaussian centers computed. (peaks separated by > 1 cell)
569         // Begin from cell having largest energy deposited This is first
570         // cluster center
571         // *****************************************************************
572         // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
573         // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
574         // SINCE WE EXPECT THE SUPERCLUSTER 
575         // TO BE A SINGLE CLUSTER
576         //*******************************************************************
577         
578         i1      = fInfcl[1][id];
579         i2      = fInfcl[2][id];
580         Int_t i12 = i1 + i2*kNDIMX;
581         
582         x[0]    = fCoord[0][i1][i2];
583         y[0]    = fCoord[1][i1][i2];
584         z[0]    = edepcell[i12];
585         
586         iord[0] = 0;
587         for(j = 1; j <= ncl[i]; j++)
588           {
589             
590             id++;
591             i1      = fInfcl[1][id];
592             i2      = fInfcl[2][id];
593             Int_t i12 = i1 + i2*kNDIMX;
594             iord[j] = j;
595             x[j]    = fCoord[0][i1][i2];
596             y[j]    = fCoord[1][i1][i2];
597             z[j]    = edepcell[i12];
598           }
599         
600         // arranging cells within supercluster in decreasing order
601         for(j = 1; j <= ncl[i];j++)
602           {
603             itest = 0;
604             ihld  = iord[j];
605             for(i1 = 0; i1 < j; i1++)
606               {
607                 if(itest == 0 && z[iord[i1]] < z[ihld])
608                   {
609                     itest = 1;
610                     for(i2 = j-1;i2 >= i1;i2--)
611                       {
612                         iord[i2+1] = iord[i2];
613                       }
614                     iord[i1] = ihld;
615                   }
616               }
617           }
618         
619         
620         // compute the number of clusters and their centers ( first
621         // guess )
622         // centers must be separated by cells having smaller ener. dep.
623         // neighbouring centers should be either strong or well-separated
624         ig     = 0;
625         xc[ig] = x[iord[0]];
626         yc[ig] = y[iord[0]];
627         zc[ig] = z[iord[0]];
628         for(j = 1; j <= ncl[i]; j++)
629           {
630             itest = -1;
631             x1    = x[iord[j]];
632             y1    = y[iord[j]];
633             for(k = 0; k <= ig; k++)
634               {
635                 x2 = xc[k];
636                 y2 = yc[k];
637                 rr = Distance(x1,y1,x2,y2);
638                 //************************************************************
639                 // finetuning cluster splitting
640                 // the numbers zc/4 and zc/10 may need to be changed. 
641                 // Also one may need to add one more layer because our 
642                 // cells are smaller in absolute scale
643                 //************************************************************
644                 
645                 
646                 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.) itest++;
647                 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.) itest++;
648                 if( rr >= 2.1)itest++;
649               }
650             
651             if(itest == ig)
652               {
653                 ig++;
654                 xc[ig] = x1;
655                 yc[ig] = y1;
656                 zc[ig] = z[iord[j]];
657               }
658           }
659         ClustDetails(ncl[i], ig, x, y ,z, xc, yc, zc, rcl, rcs, cells, 
660                      testncl, testindex);
661         
662         Int_t pp = 0;
663         for(j = 0; j <= ig; j++)
664           { 
665             clno++;
666             if (clno >= 5000)
667               {
668                 AliWarning("RefClust: Too many clusters! more than 5000");
669                 return;
670               }
671             clusdata[0] = xc[j];
672             clusdata[1] = yc[j];
673             clusdata[2] = zc[j];
674             clusdata[4] = rcl[j];
675             clusdata[5] = rcs[j];
676             if(ig == 0)
677               {
678                 clusdata[3] = ncl[i] + 1;//ajay
679                 //clusdata[3] = ncl[i] ;
680               }
681             else
682               {
683                 clusdata[3] = cells[j];
684               }
685             // cell information
686             Int_t ncellcls =  testncl[j];
687             if( ncellcls < 15 )
688               {
689                 for(Int_t kk = 1; kk <= ncellcls; kk++)
690                   {
691                     Int_t ll =  testindex[pp];
692                     clY = (Int_t)((ktwobysqrt3*y[ll])*10);
693                     clX = (Int_t)((x[ll] - clY/20.)*10);
694                     clxy[kk-1] = clX*10000 + clY ;
695                     
696
697                     //clX = (Int_t) x[ll]*10;
698                     //clY = (Int_t) y[ll]*10;
699                     //clxy[kk-1] = (Int_t) clX*10000 + clY ;
700                     pp++;
701                   }
702                 for(Int_t icltr = ncellcls ; icltr < 15; icltr++)
703                   {
704                     clxy[icltr] = -1;
705                   }
706               }
707             pmdcludata = new AliPMDcludata(clusdata, clxy);
708             fPMDclucont->Add(pmdcludata);
709           }
710         testncl.Set(0);
711         testindex.Set(0);
712       }
713     }
714   delete [] ncl;
715   delete [] iord;
716   delete [] x;
717   delete [] y;
718   delete [] z;
719   delete [] xc;
720   delete [] yc;
721   delete [] zc;
722   delete [] cells;
723   delete [] rcl;
724   delete [] rcs;
725 }
726 // ------------------------------------------------------------------------ //
727 void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust, Double_t x[], 
728                                       Double_t y[], Double_t z[],Double_t xc[],
729                                       Double_t yc[], Double_t zc[],
730                                       Double_t rcl[], Double_t rcs[], 
731                                       Double_t cells[], TArrayI &testncl,
732                                       TArrayI &testindex)
733 {
734   // function begins
735   //
736
737   Int_t kndim1 = ncell + 1;//ncell
738   Int_t kndim2 = 20;
739   Int_t kndim3 = nclust + 1;//nclust
740   
741   Int_t    i, j, k, i1, i2;
742   Double_t x1, y1, x2, y2, rr, b, c, r1, r2;
743   Double_t sumx, sumy, sumxy, sumxx, sum, sum1, sumyy;
744
745   Double_t  *str, *str1, *xcl, *ycl, *cln; 
746   Int_t    **cell;
747   Int_t    ** cluster;
748   Double_t **clustcell;
749   str  = new Double_t [kndim3];
750   str1 = new Double_t [kndim3];
751   xcl  = new Double_t [kndim3];
752   ycl  = new Double_t [kndim3];
753   cln  = new Double_t [kndim3];
754
755   clustcell = new Double_t *[kndim3];
756   cell      = new Int_t    *[kndim3];
757   cluster   = new Int_t    *[kndim1];
758   for(i = 0; i < kndim1; i++)
759     {
760       cluster[i] = new Int_t [kndim2];
761     }
762   
763   for(i = 0; i < kndim3; i++)
764     {
765       str[i]  = 0;
766       str1[i] = 0;
767       xcl[i]  = 0;
768       ycl[i]  = 0;
769       cln[i]  = 0;
770       
771       cell[i]    = new Int_t [kndim2];
772       clustcell[i] = new Double_t [kndim1];       
773       for(j = 0; j < kndim1; j++)
774         {
775           clustcell[i][j] = 0;
776         }
777       for(j = 0; j < kndim2; j++)
778         {
779           cluster[i][j] = 0;
780           cell[i][j] = 0;
781         }
782     }
783   
784   if(nclust > 0)
785     {
786       // more than one cluster
787       // checking cells shared between several  clusters.
788       // First check if the cell is within
789       // one cell unit ( nearest neighbour). Else, 
790       // if it is within 1.74 cell units ( next nearest )
791       // Else if it is upto 2 cell units etc.
792       
793       for (i = 0; i <= ncell; i++)
794         {
795           x1            = x[i];
796           y1            = y[i];
797           cluster[i][0] = 0;
798           // distance <= 1 cell unit
799           for(j = 0; j <= nclust; j++)
800             {
801               x2 = xc[j];
802               y2 = yc[j];
803               rr = Distance(x1, y1, x2, y2);
804               if(rr <= 1.)
805                 {
806                   cluster[i][0]++;
807                   i1             = cluster[i][0];
808                   cluster[i][i1] = j;
809                 }
810             }
811           // next nearest neighbour
812           if(cluster[i][0] == 0)
813             {
814               for(j=0; j<=nclust; j++)
815                 {
816                   x2 = xc[j];
817                   y2 = yc[j];
818                   rr = Distance(x1, y1, x2, y2);
819                   if(rr <= TMath::Sqrt(3.))
820                     {
821                       cluster[i][0]++;
822                       i1             = cluster[i][0];
823                       cluster[i][i1] = j;
824                     }
825                 }
826             }
827           // next-to-next nearest neighbour
828           if(cluster[i][0] == 0)
829             {
830               for(j=0; j<=nclust; j++)
831                 {
832                   x2 = xc[j];
833                   y2 = yc[j];
834                   rr = Distance(x1, y1, x2, y2);
835                   if(rr <= 2.)
836                     {
837                       cluster[i][0]++;
838                       i1             = cluster[i][0];
839                       cluster[i][i1] = j;
840                     }
841                 }
842             }
843           // one more
844           if(cluster[i][0] == 0)
845             {
846               for(j = 0; j <= nclust; j++)
847                 {
848                   x2 = xc[j];
849                   y2 = yc[j];
850                   rr = Distance(x1, y1, x2, y2);
851                   if(rr <= 2.7)
852                     {
853                       cluster[i][0]++;
854                       i1             = cluster[i][0];
855                       cluster[i][i1] = j;
856                     }
857                 }
858             }
859         }
860       
861       // computing cluster strength. Some cells are shared.
862       for(i = 0; i <= ncell; i++)
863         {
864           if(cluster[i][0] != 0)
865             {
866               i1 = cluster[i][0];
867               for(j = 1; j <= i1; j++)
868                 {
869                   i2       = cluster[i][j];
870                   str[i2] += z[i]/i1;
871                 }
872             }
873         }
874       
875       for(k = 0; k < 5; k++)
876         {
877           for(i = 0; i <= ncell; i++)
878             {
879               if(cluster[i][0] != 0)
880                 {
881                   i1=cluster[i][0];
882                   sum=0.;
883                   for(j = 1; j <= i1; j++)
884                     {
885                       sum += str[cluster[i][j]];
886                     }
887                   
888                   for(j = 1; j <= i1; j++)
889                     {
890                       i2               = cluster[i][j]; 
891                       str1[i2]        +=  z[i]*str[i2]/sum;
892                       clustcell[i2][i] = z[i]*str[i2]/sum;
893                     }
894                 }
895             }
896           
897           
898           for(j = 0; j <= nclust; j++)
899             {
900               str[j]  = str1[j];
901               str1[j] = 0.;
902             }
903         }
904       
905       for(i = 0; i <= nclust; i++)
906         {
907           sumx = 0.;
908           sumy = 0.;
909           sum  = 0.;
910           sum1 = 0.;
911           for(j = 0; j <= ncell; j++)
912             {
913               if(clustcell[i][j] != 0)
914                 {
915                   sumx  +=  clustcell[i][j]*x[j];
916                   sumy  +=  clustcell[i][j]*y[j];
917                   sum   +=  clustcell[i][j];
918                   sum1  +=  clustcell[i][j]/z[j];
919                 }
920             }
921           //** xcl and ycl are cluster centroid positions ( center of gravity )
922           
923           xcl[i] = sumx/sum;
924           ycl[i] = sumy/sum;
925           cln[i] = sum1;
926           sumxx = 0.;
927           sumyy = 0.;
928           sumxy = 0.;
929           for(j = 0; j <= ncell; j++)
930             {
931               sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
932               sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
933               sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
934             }
935           b = sumxx+sumyy;
936           c = sumxx*sumyy-sumxy*sumxy;
937           // ******************r1 and r2 are major and minor axes ( r1 > r2 ). 
938           r1 = b/2.+TMath::Sqrt(b*b/4.-c);
939           r2 = b/2.-TMath::Sqrt(b*b/4.-c);
940           // final assignments to proper external variables
941           xc[i]    = xcl[i];
942           yc[i]    = ycl[i];
943           zc[i]    = str[i];
944           cells[i] = cln[i];
945           rcl[i]   = r1;
946           rcs[i]   = r2;
947
948         }
949       
950       //To get the cell position in a cluster
951       
952       for(Int_t ii=0; ii<= ncell; ii++)
953         {
954           Int_t jj = cluster[ii][0]; 
955           for(Int_t kk=1; kk<= jj; kk++)
956             {
957               Int_t ll = cluster[ii][kk];
958               cell[ll][0]++;
959               cell[ll][cell[ll][0]] = ii;
960             }
961         }
962       
963       testncl.Set(nclust+1);
964       Int_t counter = 0;
965       
966       for(Int_t ii=0; ii <= nclust; ii++)
967         {
968           testncl[ii] =  cell[ii][0];
969           counter += testncl[ii];
970         }
971       testindex.Set(counter);
972       Int_t ll = 0;
973       for(Int_t ii=0; ii<= nclust; ii++)
974         {
975           for(Int_t jj = 1; jj<= testncl[ii]; jj++)
976             {
977               Int_t kk = cell[ii][jj];
978               testindex[ll] = kk;
979               ll++;
980             }
981         }
982       
983     }
984   else if(nclust == 0)
985     {
986       sumx = 0.;
987       sumy = 0.;
988       sum  = 0.;
989       sum1 = 0.;
990       i    = 0;
991       for(j = 0; j <= ncell; j++)
992         {
993           sumx += z[j]*x[j];
994           sumy += z[j]*y[j];
995           sum  += z[j];
996           sum1++;
997         }
998       xcl[i] = sumx/sum;
999       ycl[i] = sumy/sum;
1000       cln[i] = sum1;
1001       sumxx  = 0.;
1002       sumyy  = 0.;
1003       sumxy  = 0.;
1004       for(j = 0; j <= ncell; j++)
1005         {
1006           sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
1007           sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
1008           sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
1009         }
1010       b  = sumxx+sumyy;
1011       c  = sumxx*sumyy-sumxy*sumxy;
1012       r1 = b/2.+ TMath::Sqrt(b*b/4.-c);
1013       r2 = b/2.- TMath::Sqrt(b*b/4.-c);
1014       
1015       // To get the cell position in a cluster
1016       testncl.Set(nclust+1);
1017       //testindex.Set(ncell);
1018       testindex.Set(ncell+1);
1019       cell[0][0] = ncell + 1;//ajay
1020       //cell[0][0] = ncell;
1021       testncl[0] = cell[0][0];
1022       Int_t ll   = 0;
1023       for(Int_t ii = 1; ii <= ncell; ii++)
1024         {
1025           cell[0][ii]=ii;
1026           //clustcell[0][ii]=1.;
1027           Int_t kk = cell[0][ii];
1028           testindex[ll] = kk;
1029           ll++;
1030         }
1031       // final assignments
1032       xc[i]    = xcl[i];
1033       yc[i]    = ycl[i];
1034       //zc[i]    = str[i];//ajay
1035       zc[i]    = sum;
1036       cells[i] = cln[i];
1037       rcl[i]   = r1;
1038       rcs[i]   = r2;
1039     }
1040   for(i = 0; i < kndim3; i++)
1041     {
1042       delete [] clustcell[i];
1043       delete [] cell[i];
1044     }
1045   delete [] clustcell;
1046   delete [] cell;
1047   for(i = 0; i <kndim1 ; i++)
1048     {
1049       delete [] cluster[i];
1050     }
1051   delete [] cluster;
1052   delete [] str;
1053   delete [] str1;
1054   delete [] xcl;
1055   delete [] ycl;
1056   delete [] cln;
1057 }
1058
1059 // ------------------------------------------------------------------------ //
1060 Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
1061                                       Double_t x2, Double_t y2)
1062 {
1063   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
1064 }
1065 // ------------------------------------------------------------------------ //
1066 void AliPMDClusteringV2::SetEdepCut(Float_t decut)
1067 {
1068   fCutoff = decut;
1069 }
1070 // ------------------------------------------------------------------------ //