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 (Int_t 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
 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           Int_t 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           //Float_t cellY    = (Float_t) (ktwobysqrt3*celldumY/10);
 484           //Float_t cellX    = (Float_t) (celldumX/10 - (celldumY/2.)/10);
 485
 486           clY = (Int_t)((ktwobysqrt3*fCoord[1][i1][i2])*10);
 487           clX = (Int_t)((fCoord[0][i1][i2] - clY/20.)*10);
 488           clxy[0] = clX*10000 + clY ;
 489
 490           //clX = (Int_t) fCoord[0][i1][i2]*10;
 491           //clY = (Int_t) fCoord[1][i1][i2]*10;
 492           //clxy[0] = clX*10000 + clY;
 493
 494
 495           for(Int_t icltr = 1; icltr < 15; icltr++)
 496             {
 497               clxy[icltr] = -1;
 498             }
 499           pmdcludata  = new AliPMDcludata(clusdata,clxy);
 500           fPMDclucont->Add(pmdcludata);
 501
 502
 503         }
 504       else if(ncl[i] == 1)
 505         {
 506           // two cell super-cluster --> single cluster
 507           // cluster center is at ener. dep.-weighted mean of two cells
 508           // cluster radius == half cell dimension
 509           id++;
 510           icl++;
 511           if (clno >= 5000)
 512             {
 513               AliWarning("RefClust: Too many clusters! more than 5000");
 514               return;
 515             }
 516           clno++;
 517           i1   = fInfcl[1][id];
 518           i2   = fInfcl[2][id];
 519           Int_t i12 = i1 + i2*kNDIMX;
 520
 521           x1   = fCoord[0][i1][i2];
 522           y1   = fCoord[1][i1][i2];
 523           z1   = edepcell[i12];
 524
 525           id++;
 526           i1   = fInfcl[1][id];
 527           i2   = fInfcl[2][id];
 528           i12  = i1 + i2*kNDIMX;
 529
 530           x2   = fCoord[0][i1][i2];
 531           y2   = fCoord[1][i1][i2];
 532           z2   = edepcell[i12];
 533
 534           clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
 535           clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
 536           clusdata[2] = z1+z2;
 537           clusdata[3] = 2.;
 538           clusdata[4] = (TMath::Sqrt(z1*z2))/(z1+z2);
 539           clusdata[5] = 0.0;
 540
 541           clY = (Int_t)((ktwobysqrt3*y1)*10);
 542           clX = (Int_t)((x1 - clY/20.)*10);
 543           clxy[0] = clX*10000 + clY ;
 544
 545           //clX = (Int_t) x1*10;
 546           //clY = (Int_t) y1*10;
 547           //clxy[0] = clX*10000 + clY;
 548
 549           clY = (Int_t)((ktwobysqrt3*y2)*10);
 550           clX = (Int_t)((x2 - clY/20.)*10);
 551           clxy[1] = clX*10000 + clY ;
 552
 553           //clX = (Int_t) x2*10;
 554           //clY = (Int_t) y2*10;
 555           //clxy[1] = clX*10000 + clY;
 556
 557           for(Int_t icltr = 2; icltr < 15; icltr++)
 558             {
 559               clxy[icltr] = -1;
 560             }
 561           pmdcludata  = new AliPMDcludata(clusdata, clxy);
 562           fPMDclucont->Add(pmdcludata);
 563         }
 564       else{
 565         id++;
 566         iord[0] = 0;
 567         // super-cluster of more than two cells - broken up into smaller
 568         // clusters gaussian centers computed. (peaks separated by > 1 cell)
 569         // Begin from cell having largest energy deposited This is first
 570         // cluster center
 571         // *****************************************************************
 572         // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
 573         // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
 574         // SINCE WE EXPECT THE SUPERCLUSTER
 575         // TO BE A SINGLE CLUSTER
 576         //*******************************************************************
 577
 578         i1      = fInfcl[1][id];
 579         i2      = fInfcl[2][id];
 580         Int_t i12 = i1 + i2*kNDIMX;
 581
 582         x[0]    = fCoord[0][i1][i2];
 583         y[0]    = fCoord[1][i1][i2];
 584         z[0]    = edepcell[i12];
 585
 586         iord[0] = 0;
 587         for(j = 1; j <= ncl[i]; j++)
 588           {
 589
 590             id++;
 591             i1      = fInfcl[1][id];
 592             i2      = fInfcl[2][id];
 593             Int_t i12 = i1 + i2*kNDIMX;
 594             iord[j] = j;
 595             x[j]    = fCoord[0][i1][i2];
 596             y[j]    = fCoord[1][i1][i2];
 597             z[j]    = edepcell[i12];
 598           }
 599
 600         // arranging cells within supercluster in decreasing order
 601         for(j = 1; j <= ncl[i];j++)
 602           {
 603             itest = 0;
 604             ihld  = iord[j];
 605             for(i1 = 0; i1 < j; i1++)
 606               {
 607                 if(itest == 0 && z[iord[i1]] < z[ihld])
 608                   {
 609                     itest = 1;
 610                     for(i2 = j-1;i2 >= i1;i2--)
 611                       {
 612                         iord[i2+1] = iord[i2];
 613                       }
 614                     iord[i1] = ihld;
 615                   }
 616               }
 617           }
 618
 619
 620         // compute the number of clusters and their centers ( first
 621         // guess )
 622         // centers must be separated by cells having smaller ener. dep.
 623         // neighbouring centers should be either strong or well-separated
 624         ig     = 0;
 625         xc[ig] = x[iord[0]];
 626         yc[ig] = y[iord[0]];
 627         zc[ig] = z[iord[0]];
 628         for(j = 1; j <= ncl[i]; j++)
 629           {
 630             itest = -1;
 631             x1    = x[iord[j]];
 632             y1    = y[iord[j]];
 633             for(k = 0; k <= ig; k++)
 634               {
 635                 x2 = xc[k];
 636                 y2 = yc[k];
 637                 rr = Distance(x1,y1,x2,y2);
 638                 //************************************************************
 639                 // finetuning cluster splitting
 640                 // the numbers zc/4 and zc/10 may need to be changed.
 641                 // Also one may need to add one more layer because our
 642                 // cells are smaller in absolute scale
 643                 //************************************************************
 644
 645
 646                 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.) itest++;
 647                 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.) itest++;
 648                 if( rr >= 2.1)itest++;
 649               }
 650
 651             if(itest == ig)
 652               {
 653                 ig++;
 654                 xc[ig] = x1;
 655                 yc[ig] = y1;
 656                 zc[ig] = z[iord[j]];
 657               }
 658           }
 659         ClustDetails(ncl[i], ig, x, y ,z, xc, yc, zc, rcl, rcs, cells,
 660                      testncl, testindex);
 661
 662         Int_t pp = 0;
 663         for(j = 0; j <= ig; j++)
 664           {
 665             clno++;
 666             if (clno >= 5000)
 667               {
 668                 AliWarning("RefClust: Too many clusters! more than 5000");
 669                 return;
 670               }
 671             clusdata[0] = xc[j];
 672             clusdata[1] = yc[j];
 673             clusdata[2] = zc[j];
 674             clusdata[4] = rcl[j];
 675             clusdata[5] = rcs[j];
 676             if(ig == 0)
 677               {
 678                 clusdata[3] = ncl[i] + 1;//ajay
 679                 //clusdata[3] = ncl[i] ;
 680               }
 681             else
 682               {
 683                 clusdata[3] = cells[j];
 684               }
 685             // cell information
 686             Int_t ncellcls =  testncl[j];
 687             if( ncellcls < 15 )
 688               {
 689                 for(Int_t kk = 1; kk <= ncellcls; kk++)
 690                   {
 691                     Int_t ll =  testindex[pp];
 692                     clY = (Int_t)((ktwobysqrt3*y[ll])*10);
 693                     clX = (Int_t)((x[ll] - clY/20.)*10);
 694                     clxy[kk-1] = clX*10000 + clY ;
 695
 696
 697                     //clX = (Int_t) x[ll]*10;
 698                     //clY = (Int_t) y[ll]*10;
 699                     //clxy[kk-1] = (Int_t) clX*10000 + clY ;
 700                     pp++;
 701                   }
 702                 for(Int_t icltr = ncellcls ; icltr < 15; icltr++)
 703                   {
 704                     clxy[icltr] = -1;
 705                   }
 706               }
 707             pmdcludata = new AliPMDcludata(clusdata, clxy);
 708             fPMDclucont->Add(pmdcludata);
 709           }
 710         testncl.Set(0);
 711         testindex.Set(0);
 712       }
 713     }
 714   delete [] ncl;
 715   delete [] iord;
 716   delete [] x;
 717   delete [] y;
 718   delete [] z;
 719   delete [] xc;
 720   delete [] yc;
 721   delete [] zc;
 722   delete [] cells;
 723   delete [] rcl;
 724   delete [] rcs;
 725 }
 726 // ------------------------------------------------------------------------ //
 727 void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust, Double_t x[],
 728                                       Double_t y[], Double_t z[],Double_t xc[],
 729                                       Double_t yc[], Double_t zc[],
 730                                       Double_t rcl[], Double_t rcs[],
 731                                       Double_t cells[], TArrayI &testncl,
 732                                       TArrayI &testindex)
 733 {
 734   // function begins
 735   //
 736
 737   Int_t kndim1 = ncell + 1;//ncell
 738   Int_t kndim2 = 20;
 739   Int_t kndim3 = nclust + 1;//nclust
 740
 741   Int_t    i, j, k, i1, i2;
 742   Double_t x1, y1, x2, y2, rr, b, c, r1, r2;
 743   Double_t sumx, sumy, sumxy, sumxx, sum, sum1, sumyy;
 744
 745   Double_t  *str, *str1, *xcl, *ycl, *cln;
 746   Int_t    **cell;
 747   Int_t    ** cluster;
 748   Double_t **clustcell;
 749   str  = new Double_t [kndim3];
 750   str1 = new Double_t [kndim3];
 751   xcl  = new Double_t [kndim3];
 752   ycl  = new Double_t [kndim3];
 753   cln  = new Double_t [kndim3];
 754
 755   clustcell = new Double_t *[kndim3];
 756   cell      = new Int_t    *[kndim3];
 757   cluster   = new Int_t    *[kndim1];
 758   for(i = 0; i < kndim1; i++)
 759     {
 760       cluster[i] = new Int_t [kndim2];
 761     }
 762
 763   for(i = 0; i < kndim3; i++)
 764     {
 765       str[i]  = 0;
 766       str1[i] = 0;
 767       xcl[i]  = 0;
 768       ycl[i]  = 0;
 769       cln[i]  = 0;
 770
 771       cell[i]    = new Int_t [kndim2];
 772       clustcell[i] = new Double_t [kndim1];
 773       for(j = 0; j < kndim1; j++)
 774         {
 775           clustcell[i][j] = 0;
 776         }
 777       for(j = 0; j < kndim2; j++)
 778         {
 779           cluster[i][j] = 0;
 780           cell[i][j] = 0;
 781         }
 782     }
 783
 784   if(nclust > 0)
 785     {
 786       // more than one cluster
 787       // checking cells shared between several  clusters.
 788       // First check if the cell is within
 789       // one cell unit ( nearest neighbour). Else,
 790       // if it is within 1.74 cell units ( next nearest )
 791       // Else if it is upto 2 cell units etc.
 792
 793       for (i = 0; i <= ncell; i++)
 794         {
 795           x1            = x[i];
 796           y1            = y[i];
 797           cluster[i][0] = 0;
 798           // distance <= 1 cell unit
 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 <= 1.)
 805                 {
 806                   cluster[i][0]++;
 807                   i1             = cluster[i][0];
 808                   cluster[i][i1] = j;
 809                 }
 810             }
 811           // next nearest neighbour
 812           if(cluster[i][0] == 0)
 813             {
 814               for(j=0; j<=nclust; j++)
 815                 {
 816                   x2 = xc[j];
 817                   y2 = yc[j];
 818                   rr = Distance(x1, y1, x2, y2);
 819                   if(rr <= TMath::Sqrt(3.))
 820                     {
 821                       cluster[i][0]++;
 822                       i1             = cluster[i][0];
 823                       cluster[i][i1] = j;
 824                     }
 825                 }
 826             }
 827           // next-to-next nearest neighbour
 828           if(cluster[i][0] == 0)
 829             {
 830               for(j=0; j<=nclust; j++)
 831                 {
 832                   x2 = xc[j];
 833                   y2 = yc[j];
 834                   rr = Distance(x1, y1, x2, y2);
 835                   if(rr <= 2.)
 836                     {
 837                       cluster[i][0]++;
 838                       i1             = cluster[i][0];
 839                       cluster[i][i1] = j;
 840                     }
 841                 }
 842             }
 843           // one more
 844           if(cluster[i][0] == 0)
 845             {
 846               for(j = 0; j <= nclust; j++)
 847                 {
 848                   x2 = xc[j];
 849                   y2 = yc[j];
 850                   rr = Distance(x1, y1, x2, y2);
 851                   if(rr <= 2.7)
 852                     {
 853                       cluster[i][0]++;
 854                       i1             = cluster[i][0];
 855                       cluster[i][i1] = j;
 856                     }
 857                 }
 858             }
 859         }
 860
 861       // computing cluster strength. Some cells are shared.
 862       for(i = 0; i <= ncell; i++)
 863         {
 864           if(cluster[i][0] != 0)
 865             {
 866               i1 = cluster[i][0];
 867               for(j = 1; j <= i1; j++)
 868                 {
 869                   i2       = cluster[i][j];
 870                   str[i2] += z[i]/i1;
 871                 }
 872             }
 873         }
 874
 875       for(k = 0; k < 5; k++)
 876         {
 877           for(i = 0; i <= ncell; i++)
 878             {
 879               if(cluster[i][0] != 0)
 880                 {
 881                   i1=cluster[i][0];
 882                   sum=0.;
 883                   for(j = 1; j <= i1; j++)
 884                     {
 885                       sum += str[cluster[i][j]];
 886                     }
 887
 888                   for(j = 1; j <= i1; j++)
 889                     {
 890                       i2               = cluster[i][j];
 891                       str1[i2]        +=  z[i]*str[i2]/sum;
 892                       clustcell[i2][i] = z[i]*str[i2]/sum;
 893                     }
 894                 }
 895             }
 896
 897
 898           for(j = 0; j <= nclust; j++)
 899             {
 900               str[j]  = str1[j];
 901               str1[j] = 0.;
 902             }
 903         }
 904
 905       for(i = 0; i <= nclust; i++)
 906         {
 907           sumx = 0.;
 908           sumy = 0.;
 909           sum  = 0.;
 910           sum1 = 0.;
 911           for(j = 0; j <= ncell; j++)
 912             {
 913               if(clustcell[i][j] != 0)
 914                 {
 915                   sumx  +=  clustcell[i][j]*x[j];
 916                   sumy  +=  clustcell[i][j]*y[j];
 917                   sum   +=  clustcell[i][j];
 918                   sum1  +=  clustcell[i][j]/z[j];
 919                 }
 920             }
 921           //** xcl and ycl are cluster centroid positions ( center of gravity )
 922
 923           xcl[i] = sumx/sum;
 924           ycl[i] = sumy/sum;
 925           cln[i] = sum1;
 926           sumxx = 0.;
 927           sumyy = 0.;
 928           sumxy = 0.;
 929           for(j = 0; j <= ncell; j++)
 930             {
 931               sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
 932               sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
 933               sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
 934             }
 935           b = sumxx+sumyy;
 936           c = sumxx*sumyy-sumxy*sumxy;
 937           // ******************r1 and r2 are major and minor axes ( r1 > r2 ).
 938           r1 = b/2.+TMath::Sqrt(b*b/4.-c);
 939           r2 = b/2.-TMath::Sqrt(b*b/4.-c);
 940           // final assignments to proper external variables
 941           xc[i]    = xcl[i];
 942           yc[i]    = ycl[i];
 943           zc[i]    = str[i];
 944           cells[i] = cln[i];
 945           rcl[i]   = r1;
 946           rcs[i]   = r2;
 947
 948         }
 949
 950       //To get the cell position in a cluster
 951
 952       for(Int_t ii=0; ii<= ncell; ii++)
 953         {
 954           Int_t jj = cluster[ii][0];
 955           for(Int_t kk=1; kk<= jj; kk++)
 956             {
 957               Int_t ll = cluster[ii][kk];
 958               cell[ll][0]++;
 959               cell[ll][cell[ll][0]] = ii;
 960             }
 961         }
 962
 963       testncl.Set(nclust+1);
 964       Int_t counter = 0;
 965
 966       for(Int_t ii=0; ii <= nclust; ii++)
 967         {
 968           testncl[ii] =  cell[ii][0];
 969           counter += testncl[ii];
 970         }
 971       testindex.Set(counter);
 972       Int_t ll = 0;
 973       for(Int_t ii=0; ii<= nclust; ii++)
 974         {
 975           for(Int_t jj = 1; jj<= testncl[ii]; jj++)
 976             {
 977               Int_t kk = cell[ii][jj];
 978               testindex[ll] = kk;
 979               ll++;
 980             }
 981         }
 982
 983     }
 984   else if(nclust == 0)
 985     {
 986       sumx = 0.;
 987       sumy = 0.;
 988       sum  = 0.;
 989       sum1 = 0.;
 990       i    = 0;
 991       for(j = 0; j <= ncell; j++)
 992         {
 993           sumx += z[j]*x[j];
 994           sumy += z[j]*y[j];
 995           sum  += z[j];
 996           sum1++;
 997         }
 998       xcl[i] = sumx/sum;
 999       ycl[i] = sumy/sum;
1000       cln[i] = sum1;
1001       sumxx  = 0.;
1002       sumyy  = 0.;
1003       sumxy  = 0.;
1004       for(j = 0; j <= ncell; j++)
1005         {
1006           sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
1007           sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
1008           sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
1009         }
1010       b  = sumxx+sumyy;
1011       c  = sumxx*sumyy-sumxy*sumxy;
1012       r1 = b/2.+ TMath::Sqrt(b*b/4.-c);
1013       r2 = b/2.- TMath::Sqrt(b*b/4.-c);
1014
1015       // To get the cell position in a cluster
1016       testncl.Set(nclust+1);
1017       //testindex.Set(ncell);
1018       testindex.Set(ncell+1);
1019       cell[0][0] = ncell + 1;//ajay
1020       //cell[0][0] = ncell;
1021       testncl[0] = cell[0][0];
1022       Int_t ll   = 0;
1023       for(Int_t ii = 1; ii <= ncell; ii++)
1024         {
1025           cell[0][ii]=ii;
1026           //clustcell[0][ii]=1.;
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]    = str[i];//ajay
1035       zc[i]    = sum;
1036       cells[i] = cln[i];
1037       rcl[i]   = r1;
1038       rcs[i]   = r2;
1039     }
1040   for(i = 0; i < kndim3; i++)
1041     {
1042       delete [] clustcell[i];
1043       delete [] cell[i];
1044     }
1045   delete [] clustcell;
1046   delete [] cell;
1047   for(i = 0; i <kndim1 ; i++)
1048     {
1049       delete [] cluster[i];
1050     }
1051   delete [] cluster;
1052   delete [] str;
1053   delete [] str1;
1054   delete [] xcl;
1055   delete [] ycl;
1056   delete [] cln;
1057 }
1058
1059 // ------------------------------------------------------------------------ //
1060 Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
1061                                       Double_t x2, Double_t y2)
1062 {
1063   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
1064 }
1065 // ------------------------------------------------------------------------ //
1066 void AliPMDClusteringV2::SetEdepCut(Float_t decut)
1067 {
1068   fCutoff = decut;
1069 }
1070 // ------------------------------------------------------------------------ //