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