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