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