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 // ------------------------------------------------------------------------ //