PMD/AliPMDClusteringV2.cxx

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