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