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