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 array fEdepCell[kNDIMX][kNDIMY] and cluster information is in array
  32    fClusters[5][5000]. integer fClno gives total number of clusters in the
  33    supermodule.
  34
  35    fEdepCell, fClno  and fClusters are the only global ( public ) variables.
  36    Others are local ( private ) to the code.
  37    At the moment, the data is read for whole detector ( all supermodules
  38    and pmd as well as cpv. This will have to be modify later )
  39    LAST UPDATE  :  October 23, 2002
  40 -----------------------------------------------------------------------*/
  41
  42 #include <Riostream.h>
  43 #include <TMath.h>
  44 #include <TObjArray.h>
  45 #include <stdio.h>
  46
  47 #include "AliPMDcluster.h"
  48 #include "AliPMDClustering.h"
  49 #include "AliPMDClusteringV2.h"
  50 #include "AliLog.h"
  51
  52 ClassImp(AliPMDClusteringV2)
  53
  54 const Double_t AliPMDClusteringV2::fgkSqroot3by2=0.8660254;  // sqrt(3.)/2.
  55
  56 AliPMDClusteringV2::AliPMDClusteringV2():
  57   fClno(0),
  58   fCutoff(0.0)
  59 {
  60   for(int i = 0; i < kNDIMX; i++)
  61     {
  62       for(int j = 0; j < kNDIMY; j++)
  63         {
  64           fCoord[0][i][j] = i+j/2.;
  65           fCoord[1][i][j] = fgkSqroot3by2*j;
  66           fEdepCell[i][j] = 0;
  67         }
  68     }
  69 }
  70 // ------------------------------------------------------------------------ //
  71 AliPMDClusteringV2::~AliPMDClusteringV2()
  72 {
  73
  74 }
  75 // ------------------------------------------------------------------------ //
  76 void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96], TObjArray *pmdcont)
  77 {
  78   // main function to call other necessary functions to do clustering
  79   //
  80   AliPMDcluster *pmdcl = 0;
  81
  82   Int_t    i, i1, i2, j, nmx1, incr, id, jd;
  83   Int_t    celldataX[15], celldataY[15];
  84   Float_t  clusdata[6];
  85   Double_t cutoff, ave;
  86
  87   const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
  88
  89   Int_t ndimXr =0;
  90   Int_t ndimYr =0;
  91
  92   if (ismn < 12)
  93     {
  94       ndimXr = 96;
  95       ndimYr = 48;
  96     }
  97   else if (ismn >= 12 && ismn <= 23)
  98     {
  99       ndimXr = 48;
 100       ndimYr = 96;
 101     }
 102
 103   for (Int_t i =0; i < kNDIMX; i++)
 104     {
 105       for (Int_t j =0; j < kNDIMY; j++)
 106         {
 107           fEdepCell[i][j] = 0;
 108         }
 109     }
 110
 111
 112   for (id = 0; id < ndimXr; id++)
 113     {
 114       for (jd = 0; jd < ndimYr; jd++)
 115         {
 116           j=jd;
 117           i=id+(ndimYr/2-1)-(jd/2);
 118
 119           if (ismn < 12)
 120             {
 121               fEdepCell[i][j] = celladc[jd][id];
 122             }
 123           else if (ismn >= 12 && ismn <= 23)
 124             {
 125               fEdepCell[i][j] = celladc[id][jd];
 126             }
 127
 128         }
 129     }
 130
 131   Order();          // order the data
 132   cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
 133   ave=0.;
 134   nmx1=-1;
 135
 136   for(j=0;j<kNMX; j++)
 137     {
 138       i1 = fIord[0][j];
 139       i2 = fIord[1][j];
 140       if (fEdepCell[i1][i2] > 0.) {ave = ave + fEdepCell[i1][i2];}
 141       if (fEdepCell[i1][i2] > cutoff ) nmx1 = nmx1 + 1;
 142     }
 143   // nmx1 --- number of cells having ener dep >= cutoff
 144
 145   AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
 146
 147   if (nmx1 == 0) nmx1 = 1;
 148   ave=ave/nmx1;
 149
 150   AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
 151                   kNMX,ave));
 152
 153   incr = CrClust(ave, cutoff, nmx1);
 154   RefClust(incr);
 155
 156   AliDebug(1,Form("Detector Plane = %d  Serial Module No = %d Number of clusters = %d",idet, ismn, fClno));
 157
 158   for(i1=0; i1<=fClno; i1++)
 159     {
 160       Float_t cluXC    = (Float_t) fClusters[0][i1];
 161       Float_t cluYC    = (Float_t) fClusters[1][i1];
 162       Float_t cluADC   = (Float_t) fClusters[2][i1];
 163       Float_t cluCELLS = (Float_t) fClusters[3][i1];
 164       Float_t sigmaX   = (Float_t) fClusters[4][i1];
 165       Float_t sigmaY   = (Float_t) fClusters[5][i1];
 166       Float_t cluY0    = ktwobysqrt3*cluYC;
 167       Float_t cluX0    = cluXC - cluY0/2.;
 168       //
 169       // Cluster X centroid is back transformed
 170       //
 171       if (ismn < 12)
 172         {
 173           clusdata[0] = cluX0 - (24-1) + cluY0/2.;
 174         }
 175       else if (ismn >= 12 && ismn <= 23)
 176         {
 177           clusdata[0] = cluX0 - (48-1) + cluY0/2.;
 178         }
 179
 180       clusdata[1]      = cluY0;
 181       clusdata[2]      = cluADC;
 182       clusdata[3]      = cluCELLS;
 183       clusdata[4]      = sigmaX;
 184       clusdata[5]      = sigmaY;
 185
 186       //
 187       // Cells associated with a cluster
 188       //
 189       for (Int_t ihit = 0; ihit < 15; ihit++)
 190         {
 191           celldataX[ihit] = 1;  // dummy nos. -- will be changed
 192           celldataY[ihit] = 1;  // dummy nos. -- will be changed
 193         }
 194
 195       pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
 196       pmdcont->Add(pmdcl);
 197     }
 198 }
 199 // ------------------------------------------------------------------------ //
 200 void AliPMDClusteringV2::Order()
 201 {
 202   // Sorting algorithm
 203   // sorts the ADC values from higher to lower
 204   //
 205   double dd[kNMX];
 206   // matrix fEdepCell converted into
 207   // one dimensional array dd. adum a place holder for double
 208   int i, j, i1, i2, iord1[kNMX];
 209   // information of
 210   // ordering is stored in iord1, original array not ordered
 211   //
 212   // define arrays dd and iord1
 213   for(i1=0; i1 < kNDIMX; i1++)
 214     {
 215       for(i2=0; i2 < kNDIMY; i2++)
 216         {
 217           i        = i1 + i2*kNDIMX;
 218           iord1[i] = i;
 219           dd[i]    = fEdepCell[i1][i2];
 220         }
 221     }
 222   // sort and store sorting information in iord1
 223
 224   TMath::Sort(kNMX,dd,iord1);
 225
 226   // store the sorted information in fIord for later use
 227   for(i=0; i<kNMX; i++)
 228     {
 229       j  = iord1[i];
 230       i2 = j/kNDIMX;
 231       i1 = j-i2*kNDIMX;
 232       fIord[0][i]=i1;
 233       fIord[1][i]=i2;
 234     }
 235 }
 236 // ------------------------------------------------------------------------ //
 237 Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1)
 238 {
 239   // Does crude clustering
 240   // Finds out only the big patch by just searching the
 241   // connected cells
 242   //
 243
 244   int i,j,k,id1,id2,icl, numcell;
 245   int jd1,jd2, icell, cellcount;
 246   int clust[2][5000];
 247   static int neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
 248
 249   // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
 250   // cell. There are six neighbours.
 251   // cellcount --- total number of cells having nonzero ener dep
 252   // numcell --- number of cells in a given supercluster
 253   // ofstream ofl0("cells_loc",ios::out);
 254   // initialize fInfocl[2][kNDIMX][kNDIMY]
 255
 256   AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
 257
 258   for (j=0; j < kNDIMX; j++){
 259     for(k=0; k < kNDIMY; k++){
 260       fInfocl[0][j][k] = 0;
 261       fInfocl[1][j][k] = 0;
 262     }
 263   }
 264   for(i=0; i < kNMX; i++){
 265     fInfcl[0][i] = -1;
 266     id1=fIord[0][i];
 267     id2=fIord[1][i];
 268     if(fEdepCell[id1][id2] <= cutoff){fInfocl[0][id1][id2]=-1;}
 269   }
 270   // ---------------------------------------------------------------
 271   // crude clustering begins. Start with cell having largest adc
 272   // count and loop over the cells in descending order of adc count
 273   // ---------------------------------------------------------------
 274   icl=-1;
 275   cellcount=-1;
 276   for(icell=0; icell <= nmx1; icell++){
 277     id1=fIord[0][icell];
 278     id2=fIord[1][icell];
 279     if(fInfocl[0][id1][id2] == 0 ){
 280       // ---------------------------------------------------------------
 281       // icl -- cluster #, numcell -- # of cells in it, clust -- stores
 282       // coordinates of the cells in a cluster, fInfocl[0][i1][i2] is 1 for
 283       // primary and 2 for secondary cells,
 284       // fInfocl[1][i1][i2] stores cluster #
 285       // ---------------------------------------------------------------
 286       icl=icl+1;
 287       numcell=0;
 288       cellcount = cellcount + 1;
 289       fInfocl[0][id1][id2]=1;
 290       fInfocl[1][id1][id2]=icl;
 291       fInfcl[0][cellcount]=icl;
 292       fInfcl[1][cellcount]=id1;
 293       fInfcl[2][cellcount]=id2;
 294
 295       clust[0][numcell]=id1;
 296       clust[1][numcell]=id2;
 297       for(i=1; i<5000; i++)clust[0][i] = -1;
 298       // ---------------------------------------------------------------
 299       // check for adc count in neib. cells. If ne 0 put it in this clust
 300       // ---------------------------------------------------------------
 301       for(i=0; i<6; i++){
 302         jd1=id1+neibx[i];
 303         jd2=id2+neiby[i];
 304         if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
 305             fInfocl[0][jd1][jd2] == 0){
 306           numcell=numcell+1;
 307           fInfocl[0][jd1][jd2]=2;
 308           fInfocl[1][jd1][jd2]=icl;
 309           clust[0][numcell]=jd1;
 310           clust[1][numcell]=jd2;
 311           cellcount=cellcount+1;
 312           fInfcl[0][cellcount]=icl;
 313           fInfcl[1][cellcount]=jd1;
 314           fInfcl[2][cellcount]=jd2;
 315         }
 316       }
 317       // ---------------------------------------------------------------
 318       // check adc count for neighbour's neighbours recursively and
 319       // if nonzero, add these to the cluster.
 320       // ---------------------------------------------------------------
 321       for(i=1;i < 5000;i++){
 322         if(clust[0][i] != -1){
 323           id1=clust[0][i];
 324           id2=clust[1][i];
 325           for(j=0; j<6 ; j++){
 326             jd1=id1+neibx[j];
 327             jd2=id2+neiby[j];
 328             if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
 329                 fInfocl[0][jd1][jd2] == 0 ){
 330               fInfocl[0][jd1][jd2] = 2;
 331               fInfocl[1][jd1][jd2] = icl;
 332               numcell              = numcell + 1;
 333               clust[0][numcell]    = jd1;
 334               clust[1][numcell]    = jd2;
 335               cellcount            = cellcount+1;
 336               fInfcl[0][cellcount] = icl;
 337               fInfcl[1][cellcount] = jd1;
 338               fInfcl[2][cellcount] = jd2;
 339             }
 340           }
 341         }
 342       }
 343     }
 344   }
 345   //  for(icell=0; icell<=cellcount; icell++){
 346   //    ofl0 << fInfcl[0][icell] << " " << fInfcl[1][icell] << " " <<
 347   //      fInfcl[2][icell] << endl;
 348   //  }
 349   return cellcount;
 350 }
 351 // ------------------------------------------------------------------------ //
 352 void AliPMDClusteringV2::RefClust(Int_t incr)
 353 {
 354   // Does the refining of clusters
 355   // Takes the big patch and does gaussian fitting and
 356   // finds out the more refined clusters
 357   //
 358
 359   const Int_t kndim = 4500;
 360
 361   int i, j, k, i1, i2, id, icl, itest;
 362   int ihld;
 363   int ig, nsupcl;
 364   int ncl[kndim], iord[kndim];
 365
 366   double x1, y1, z1, x2, y2, z2;
 367   double rr;
 368
 369   double x[kndim], y[kndim], z[kndim];
 370   double xc[kndim], yc[kndim], zc[kndim], cells[kndim];
 371   double rcl[kndim], rcs[kndim];
 372
 373   // fClno counts the final clusters
 374   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
 375   // x, y and z store (x,y) coordinates of and energy deposited in a cell
 376   // xc, yc store (x,y) coordinates of the cluster center
 377   // zc stores the energy deposited in a cluster
 378   // rc is cluster radius
 379   // finally the cluster information is put in 2-dimensional array clusters
 380   // ofstream ofl1("checking.5",ios::app);
 381
 382   fClno  = -1;
 383   nsupcl = -1;
 384   for(i=0; i<4500; i++){ncl[i]=-1;}
 385   for(i=0; i<incr; i++){
 386     if(fInfcl[0][i] != nsupcl){ nsupcl=nsupcl+1; }
 387     if (nsupcl > 4500) {
 388       AliWarning("RefClust: Too many superclusters!");
 389       nsupcl = 4500;
 390       break;
 391     }
 392     ncl[nsupcl]=ncl[nsupcl]+1;
 393   }
 394
 395   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
 396                   incr+1,nsupcl+1));
 397
 398   id=-1;
 399   icl=-1;
 400   for(i=0; i<nsupcl; i++){
 401     if(ncl[i] == 0){
 402       id++;
 403       icl++;
 404       // one  cell super-clusters --> single cluster
 405       // cluster center at the centyer of the cell
 406       // cluster radius = half cell dimension
 407       if (fClno >= 5000) {
 408         AliWarning("RefClust: Too many clusters! more than 5000");
 409         return;
 410       }
 411       fClno++;
 412       i1 = fInfcl[1][id];
 413       i2 = fInfcl[2][id];
 414       fClusters[0][fClno] = fCoord[0][i1][i2];
 415       fClusters[1][fClno] = fCoord[1][i1][i2];
 416       fClusters[2][fClno] = fEdepCell[i1][i2];
 417       fClusters[3][fClno] = 1.;
 418       fClusters[4][fClno] = 0.0;
 419       fClusters[5][fClno] = 0.0;
 420       //ofl1 << icl << " " << fCoord[0][i1][i2] << " " << fCoord[1][i1][i2] <<
 421       //" " << fEdepCell[i1][i2] << " " << fClusters[3][fClno] <<endl;
 422     }else if(ncl[i] == 1){
 423       // two cell super-cluster --> single cluster
 424       // cluster center is at ener. dep.-weighted mean of two cells
 425       // cluster radius == half cell dimension
 426       id++;
 427       icl++;
 428       if (fClno >= 5000) {
 429         AliWarning("RefClust: Too many clusters! more than 5000");
 430         return;
 431       }
 432       fClno++;
 433       i1   = fInfcl[1][id];
 434       i2   = fInfcl[2][id];
 435       x1   = fCoord[0][i1][i2];
 436       y1   = fCoord[1][i1][i2];
 437       z1   = fEdepCell[i1][i2];
 438
 439       id++;
 440       i1   = fInfcl[1][id];
 441       i2   = fInfcl[2][id];
 442       x2   = fCoord[0][i1][i2];
 443       y2   = fCoord[1][i1][i2];
 444       z2   = fEdepCell[i1][i2];
 445
 446       fClusters[0][fClno] = (x1*z1+x2*z2)/(z1+z2);
 447       fClusters[1][fClno] = (y1*z1+y2*z2)/(z1+z2);
 448       fClusters[2][fClno] = z1+z2;
 449       fClusters[3][fClno] = 2.;
 450       fClusters[4][fClno] = sqrt(z1*z2)/(z1+z2);
 451       fClusters[5][fClno] = 0;  // sigma large nonzero, sigma small zero
 452
 453       //ofl1 << icl << " " << fClusters[0][fClno] << " " << fClusters[1][fClno]
 454       //   << " " << fClusters[2][fClno] << " " <<fClusters[3][fClno] <<endl;
 455     }
 456     else{
 457       id      = id + 1;
 458       iord[0] = 0;
 459       // super-cluster of more than two cells - broken up into smaller
 460       // clusters gaussian centers computed. (peaks separated by > 1 cell)
 461       // Begin from cell having largest energy deposited This is first
 462       // cluster center
 463       // *****************************************************************
 464       // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
 465       // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
 466       // SINCE WE EXPECT THE SUPERCLUSTER
 467       // TO BE A SINGLE CLUSTER
 468       //*******************************************************************
 469
 470       i1      = fInfcl[1][id];
 471       i2      = fInfcl[2][id];
 472       x[0]    = fCoord[0][i1][i2];
 473       y[0]    = fCoord[1][i1][i2];
 474       z[0]    = fEdepCell[i1][i2];
 475       iord[0] = 0;
 476       for(j=1;j<=ncl[i];j++){
 477
 478         id      = id + 1;
 479         i1      = fInfcl[1][id];
 480         i2      = fInfcl[2][id];
 481         iord[j] = j;
 482         x[j]    = fCoord[0][i1][i2];
 483         y[j]    = fCoord[1][i1][i2];
 484         z[j]    = fEdepCell[i1][i2];
 485       }
 486       // arranging cells within supercluster in decreasing order
 487       for(j=1;j<=ncl[i];j++)
 488         {
 489           itest = 0;
 490           ihld  = iord[j];
 491           for(i1=0; i1<j; i1++)
 492             {
 493               if(itest == 0 && z[iord[i1]] < z[ihld])
 494                 {
 495                   itest = 1;
 496                   for(i2=j-1;i2>=i1;i2--)
 497                     {
 498                       iord[i2+1] = iord[i2];
 499                     }
 500                   iord[i1] = ihld;
 501                 }
 502             }
 503         }
 504
 505       // compute the number of clusters and their centers ( first
 506       // guess )
 507       // centers must be separated by cells having smaller ener. dep.
 508       // neighbouring centers should be either strong or well-separated
 509       ig     = 0;
 510       xc[ig] = x[iord[0]];
 511       yc[ig] = y[iord[0]];
 512       zc[ig] = z[iord[0]];
 513       for(j=1;j<=ncl[i];j++){
 514         itest = -1;
 515         x1    = x[iord[j]];
 516         y1    = y[iord[j]];
 517         for(k=0;k<=ig;k++){
 518           x2 = xc[k];
 519           y2 = yc[k];
 520           rr = Distance(x1,y1,x2,y2);
 521           //***************************************************************
 522           // finetuning cluster splitting
 523           // the numbers zc/4 and zc/10 may need to be changed.
 524           // Also one may need to add one more layer because our
 525           // cells are smaller in absolute scale
 526           //****************************************************************
 527
 528
 529           if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)
 530             itest++;
 531           if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)
 532             itest++;
 533           if( rr >= 2.1)itest++;
 534         }
 535         if(itest == ig){
 536           ig++;
 537           xc[ig] = x1;
 538           yc[ig] = y1;
 539           zc[ig] = z[iord[j]];
 540         }
 541       }
 542
 543       ClustDetails(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0],
 544                    rcl[0], rcs[0], cells[0]);
 545
 546       icl = icl + ig + 1;
 547
 548       for(j=0; j<=ig; j++)
 549         {
 550           if (fClno >= 5000)
 551             {
 552               AliWarning("RefClust: Too many clusters! more than 5000");
 553               return;
 554             }
 555           fClno++;
 556           fClusters[0][fClno] = xc[j];
 557           fClusters[1][fClno] = yc[j];
 558           fClusters[2][fClno] = zc[j];
 559           fClusters[4][fClno] = rcl[j];
 560           fClusters[5][fClno] = rcs[j];
 561           if(ig == 0)
 562             {
 563               fClusters[3][fClno] = ncl[i];
 564             }
 565           else
 566             {
 567               fClusters[3][fClno] = cells[j];
 568             }
 569         }
 570
 571
 572     }
 573   }
 574 }
 575
 576
 577 // ------------------------------------------------------------------------ //
 578
 579 void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust,
 580                                       Double_t &x, Double_t &y, Double_t &z,
 581                                       Double_t &xc, Double_t &yc, Double_t &zc,
 582                                       Double_t &rcl, Double_t &rcs,
 583                                       Double_t &cells)
 584 {
 585   // function begins
 586   //
 587
 588   const Int_t kndim1 = 4500;
 589   const Int_t kndim2 = 10;
 590   const Int_t kndim3 = 100;
 591
 592   int i, j, k, i1, i2;
 593   int cluster[kndim1][kndim2];
 594
 595   double x1, y1, x2, y2, rr;
 596   double sumx, sumy, sumxy, sumxx;
 597   double sum, sum1, sumyy;
 598   double b, c, r1, r2;
 599
 600   double xx[kndim1], yy[kndim1], zz[kndim1];
 601   double xxc[kndim1], yyc[kndim1];
 602
 603   double str[kndim1];
 604
 605   double str1[kndim1];
 606   double xcl[kndim1], ycl[kndim1], cln[kndim1];
 607   double clustcell[kndim1][kndim3];
 608
 609   for(i=0; i<=nclust; i++){
 610    xxc[i]=*(&xc+i);
 611    yyc[i]=*(&yc+i);
 612    str[i]=0.;
 613    str1[i]=0.;
 614   }
 615   for(i=0; i<=ncell; i++){
 616     xx[i]=*(&x+i);
 617     yy[i]=*(&y+i);
 618     zz[i]=*(&z+i);
 619   }
 620   // INITIALIZE
 621   for(i=0; i<4500; i++){
 622     for(j=0; j<100; j++){
 623       clustcell[i][j]=0.;
 624     }
 625   }
 626
 627   // INITIALIZE
 628   for(i=0;i<4500;i++){
 629     for(j=0;j<10;j++){
 630       cluster[i][j]=0;
 631     }
 632   }
 633
 634
 635   if(nclust > 0){
 636     // more than one cluster
 637     // checking cells shared between several  clusters.
 638     // First check if the cell is within
 639     // one cell unit ( nearest neighbour). Else,
 640     // if it is within 1.74 cell units ( next nearest )
 641     // Else if it is upto 2 cell units etc.
 642
 643     for (i=0; i<=ncell; i++){
 644       x1            = xx[i];
 645       y1            = yy[i];
 646       cluster[i][0] = 0;
 647       // distance <= 1 cell unit
 648       for(j=0; j<=nclust; j++)
 649         {
 650           x2 = xxc[j];
 651           y2 = yyc[j];
 652           rr = Distance(x1, y1, x2, y2);
 653           if(rr <= 1.)
 654             {
 655               cluster[i][0]++;
 656               i1             = cluster[i][0];
 657               cluster[i][i1] = j;
 658             }
 659         }
 660       // next nearest neighbour
 661       if(cluster[i][0] == 0)
 662         {
 663           for(j=0; j<=nclust; j++)
 664             {
 665               x2 = xxc[j];
 666               y2 = yyc[j];
 667               rr = Distance(x1, y1, x2, y2);
 668               if(rr <= sqrt(3.))
 669                 {
 670                   cluster[i][0]++;
 671                   i1             = cluster[i][0];
 672                   cluster[i][i1] = j;
 673                 }
 674             }
 675         }
 676       // next-to-next nearest neighbour
 677       if(cluster[i][0] == 0)
 678         {
 679           for(j=0; j<=nclust; j++)
 680             {
 681               x2 = xxc[j];
 682               y2 = yyc[j];
 683               rr = Distance(x1, y1, x2, y2);
 684               if(rr <= 2.)
 685                 {
 686                   cluster[i][0]++;
 687                   i1 = cluster[i][0];
 688                   cluster[i][i1] = j;
 689                 }
 690             }
 691         }
 692       // one more
 693       if(cluster[i][0] == 0)
 694         {
 695           for(j=0; j<=nclust; j++)
 696             {
 697               x2 = xxc[j];
 698               y2 = yyc[j];
 699               rr = Distance(x1, y1, x2, y2);
 700               if(rr <= 2.7)
 701                 {
 702                   cluster[i][0]++;
 703                   i1 = cluster[i][0];
 704                   cluster[i][i1] = j;
 705                 }
 706             }
 707         }
 708     }
 709
 710
 711     // computing cluster strength. Some cells are shared.
 712     for(i=0; i<=ncell; i++){
 713       if(cluster[i][0] != 0){
 714         i1 = cluster[i][0];
 715         for(j=1; j<=i1; j++){
 716           i2      = cluster[i][j];
 717           str[i2] = str[i2]+zz[i]/i1;
 718         }
 719       }
 720     }
 721
 722     for(k=0; k<5; k++)
 723       {
 724         for(i=0; i<=ncell; i++)
 725           {
 726             if(cluster[i][0] != 0)
 727               {
 728                 i1=cluster[i][0];
 729                 sum=0.;
 730                 for(j=1; j<=i1; j++)
 731                   {
 732                     sum=sum+str[cluster[i][j]];
 733                   }
 734
 735                 for(j=1; j<=i1; j++)
 736                   {
 737                     i2 = cluster[i][j];
 738                     str1[i2]         = str1[i2] + zz[i]*str[i2]/sum;
 739                     clustcell[i2][i] = zz[i]*str[i2]/sum;
 740                   }
 741               }
 742           }
 743
 744
 745           for(j=0; j<=nclust; j++)
 746             {
 747               str[j]=str1[j];
 748               str1[j]=0.;
 749             }
 750       }
 751
 752     for(i=0; i<=nclust; i++){
 753       sumx = 0.;
 754       sumy = 0.;
 755       sum  = 0.;
 756       sum1 = 0.;
 757       for(j=0; j<=ncell; j++){
 758         if(clustcell[i][j] != 0){
 759           sumx = sumx+clustcell[i][j]*xx[j];
 760           sumy = sumy+clustcell[i][j]*yy[j];
 761           sum  = sum+clustcell[i][j];
 762           sum1 = sum1+clustcell[i][j]/zz[j];
 763         }
 764       }
 765       //***** xcl and ycl are cluster centroid positions ( center of gravity )
 766
 767       xcl[i] = sumx/sum;
 768       ycl[i] = sumy/sum;
 769       cln[i] = sum1;
 770       sumxx = 0.;
 771       sumyy = 0.;
 772       sumxy = 0.;
 773       for(j=0; j<=ncell; j++){
 774         sumxx = sumxx+clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
 775         sumyy = sumyy+clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
 776         sumxy = sumxy+clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
 777       }
 778       b = sumxx+sumyy;
 779       c = sumxx*sumyy-sumxy*sumxy;
 780       // ******************r1 and r2 are major and minor axes ( r1 > r2 ).
 781       r1 = b/2.+sqrt(b*b/4.-c);
 782       r2 = b/2.-sqrt(b*b/4.-c);
 783       // final assignments to proper external variables
 784       *(&xc + i) = xcl[i];
 785       *(&yc + i) = ycl[i];
 786       *(&zc + i) = str[i];
 787       *(&cells + i) = cln[i];
 788       *(&rcl+i) = r1;
 789       *(&rcs+i) = r2;
 790     }
 791   }else{
 792     sumx = 0.;
 793     sumy = 0.;
 794     sum  = 0.;
 795     sum1 = 0.;
 796     i    = 0;
 797     for(j=0; j<=ncell; j++){
 798       sumx = sumx+zz[j]*xx[j];
 799       sumy = sumy+zz[j]*yy[j];
 800       sum  = sum+zz[j];
 801       sum1 = sum1+1.;
 802     }
 803     xcl[i] = sumx/sum;
 804     ycl[i] = sumy/sum;
 805     cln[i] = sum1;
 806     sumxx  = 0.;
 807     sumyy  = 0.;
 808     sumxy  = 0.;
 809     for(j=0; j<=ncell; j++){
 810       sumxx = sumxx+clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
 811       sumyy = sumyy+clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
 812       sumxy = sumxy+clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
 813     }
 814     b  = sumxx+sumyy;
 815     c  = sumxx*sumyy-sumxy*sumxy;
 816     r1 = b/2.+sqrt(b*b/4.-c);
 817     r2 = b/2.-sqrt(b*b/4.-c);
 818     // final assignments
 819     *(&xc + i)    = xcl[i];
 820     *(&yc + i)    = ycl[i];
 821     *(&zc + i)    = str[i];
 822     *(&cells + i) = cln[i];
 823     *(&rcl+i)     = r1;
 824     *(&rcs+i)     = r2;
 825   }
 826 }
 827
 828 // ------------------------------------------------------------------------ //
 829 Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
 830                                       Double_t x2, Double_t y2)
 831 {
 832   return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
 833 }
 834 // ------------------------------------------------------------------------ //
 835 void AliPMDClusteringV2::SetEdepCut(Float_t decut)
 836 {
 837   fCutoff = decut;
 838 }
 839 // ------------------------------------------------------------------------ //