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