PMD/AliPMDClusteringV2.cxx

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