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