PMD/AliPMDClusteringV2.cxx

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