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