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