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