PMD/AliPMDClusteringV1.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 : PMDClusteringV1.cxx, Version 00      //
  21 //                                                     //
  22 //  Date   : September 26 2002                         //
  23 //                                                     //
  24 //  clustering code for alice pmd                      //
  25 //                                                     //
  26 //-----------------------------------------------------//
  27
  28 /* --------------------------------------------------------------------
  29    Code developed by S. C. Phatak, Institute of Physics,
  30    Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
  31    ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
  32    builds up superclusters and breaks them into clusters. The input is
  33    in array edepcell[kNMX] and cluster information is in a
  34    TObjarray. Integer clno gives total number of clusters in the
  35    supermodule.
  36
  37    fClusters is the only global ( public ) variables.
  38    Others are local ( private ) to the code.
  39    At the moment, the data is read for whole detector ( all supermodules
  40    and pmd as well as cpv. This will have to be modify later )
  41    LAST UPDATE  :  October 23, 2002
  42 -----------------------------------------------------------------------*/
  43
  44 #include <Riostream.h>
  45 #include <TMath.h>
  46 #include <TNtuple.h>
  47 #include <TObjArray.h>
  48 #include "TRandom.h"
  49 #include <stdio.h>
  50
  51 #include "AliPMDcludata.h"
  52 #include "AliPMDcluster.h"
  53 #include "AliPMDClustering.h"
  54 #include "AliPMDClusteringV1.h"
  55 #include "AliLog.h"
  56
  57 ClassImp(AliPMDClusteringV1)
  58
  59 const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254;  // sqrt(3.)/2.
  60
  61 AliPMDClusteringV1::AliPMDClusteringV1():
  62   fPMDclucont(new TObjArray()),
  63   fCutoff(0.0)
  64 {
  65   for(Int_t i = 0; i < kNDIMX; i++)
  66     {
  67       for(Int_t j = 0; j < kNDIMY; j++)
  68         {
  69           fCoord[0][i][j] = i+j/2.;
  70           fCoord[1][i][j] = fgkSqroot3by2*j;
  71         }
  72     }
  73 }
  74 // ------------------------------------------------------------------------ //
  75 AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
  76   AliPMDClustering(pmdclv1),
  77   fPMDclucont(0),
  78   fCutoff(0)
  79 {
  80   // copy constructor
  81   AliError("Copy constructor not allowed ");
  82
  83 }
  84 // ------------------------------------------------------------------------ //
  85 AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
  86 {
  87   // copy constructor
  88   AliError("Assignment operator not allowed ");
  89   return *this;
  90 }
  91 // ------------------------------------------------------------------------ //
  92 AliPMDClusteringV1::~AliPMDClusteringV1()
  93 {
  94   delete fPMDclucont;
  95 }
  96 // ------------------------------------------------------------------------ //
  97 void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
  98                                  Double_t celladc[48][96], TObjArray *pmdcont)
  99 {
 100   // main function to call other necessary functions to do clustering
 101   //
 102
 103   AliPMDcluster *pmdcl = 0;
 104
 105   Int_t    i,  j, nmx1, incr, id, jd;
 106   Int_t    celldataX[15], celldataY[15];
 107   Float_t  clusdata[6];
 108   Double_t cutoff, ave;
 109   Double_t edepcell[kNMX];
 110
 111   const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
 112
 113   // ndimXr and ndimYr are different because of different module size
 114
 115   Int_t ndimXr = 0;
 116   Int_t ndimYr = 0;
 117
 118   if (ismn < 12)
 119     {
 120       ndimXr = 96;
 121       ndimYr = 48;
 122     }
 123   else if (ismn >= 12 && ismn <= 23)
 124     {
 125       ndimXr = 48;
 126       ndimYr = 96;
 127     }
 128
 129   Int_t kk = 0;
 130   for (Int_t i = 0; i < kNDIMX; i++)
 131     {
 132       for (Int_t j = 0; j < kNDIMY; j++)
 133         {
 134           fCellTrNo[i][j] = -1;
 135           edepcell[kk] = 0.;
 136           kk++;
 137         }
 138     }
 139
 140   for (id = 0; id < ndimXr; id++)
 141     {
 142       for (jd = 0; jd < ndimYr; jd++)
 143         {
 144           j = jd;
 145           i = id+(ndimYr/2-1)-(jd/2);
 146
 147           //PH    Int_t ij = i + j*kNDIMX;
 148           Int_t ij = i + j*ndimXr;
 149
 150           if (ismn < 12)
 151             {
 152               edepcell[ij]    = celladc[jd][id];
 153               fCellTrNo[i][j] = jd*10000+id;  // for association
 154             }
 155           else if (ismn >= 12 && ismn <= 23)
 156             {
 157               edepcell[ij]    = celladc[id][jd];
 158               fCellTrNo[i][j] = id*10000+jd;  // for association
 159             }
 160
 161         }
 162     }
 163
 164   Int_t iord1[kNMX];
 165   TMath::Sort(kNMX,edepcell,iord1);// order the data
 166   cutoff = fCutoff;                // cutoff to discard cells
 167   ave  = 0.;
 168   nmx1 = -1;
 169   for(i = 0;i < kNMX; i++)
 170     {
 171       if(edepcell[i] > 0.)
 172         {
 173           ave += edepcell[i];
 174         }
 175       if(edepcell[i] > cutoff )
 176         {
 177           nmx1++;
 178         }
 179     }
 180
 181   AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
 182
 183   if (nmx1 == 0) nmx1 = 1;
 184   ave = ave/nmx1;
 185
 186   AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
 187                   kNMX,ave));
 188
 189   incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
 190   RefClust(incr,edepcell);
 191   Int_t nentries1 = fPMDclucont->GetEntries();
 192   AliDebug(1,Form("Detector Plane = %d  Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
 193   AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
 194
 195   for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 196     {
 197       AliPMDcludata *pmdcludata =
 198         (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
 199       Float_t cluXC    = pmdcludata->GetClusX();
 200       Float_t cluYC    = pmdcludata->GetClusY();
 201       Float_t cluADC   = pmdcludata->GetClusADC();
 202       Float_t cluCELLS = pmdcludata->GetClusCells();
 203       Float_t cluSIGX  = pmdcludata->GetClusSigmaX();
 204       Float_t cluSIGY  = pmdcludata->GetClusSigmaY();
 205
 206       Float_t cluY0    = ktwobysqrt3*cluYC;
 207       Float_t cluX0    = cluXC - cluY0/2.;
 208
 209       //
 210       // Cluster X centroid is back transformed
 211       //
 212       if (ismn < 12)
 213         {
 214           clusdata[0] = cluX0 - (24-1) + cluY0/2.;
 215         }
 216       else if ( ismn >= 12 && ismn <= 23)
 217         {
 218           clusdata[0] = cluX0 - (48-1) + cluY0/2.;
 219         }
 220
 221       clusdata[1]     = cluY0;
 222       clusdata[2]     = cluADC;
 223       clusdata[3]     = cluCELLS;
 224       clusdata[4]     = cluSIGX;
 225       clusdata[5]     = cluSIGY;
 226
 227       //
 228       // Cells associated with a cluster
 229       //
 230
 231       for (Int_t ihit = 0; ihit < 15; ihit++)
 232         {
 233           if (ismn < 12)
 234             {
 235               celldataX[ihit] = pmdcludata->GetCellXY(ihit)%10000;
 236               celldataY[ihit] = pmdcludata->GetCellXY(ihit)/10000;
 237             }
 238           else if (ismn >= 12 && ismn <= 23)
 239             {
 240               celldataX[ihit] = pmdcludata->GetCellXY(ihit)/10000;
 241               celldataY[ihit] = pmdcludata->GetCellXY(ihit)%10000;
 242             }
 243         }
 244       pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
 245       pmdcont->Add(pmdcl);
 246     }
 247
 248   fPMDclucont->Clear();
 249
 250 }
 251 // ------------------------------------------------------------------------ //
 252 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
 253                                   Int_t iord1[], Double_t edepcell[])
 254 {
 255   // Does crude clustering
 256   // Finds out only the big patch by just searching the
 257   // connected cells
 258   //
 259   Int_t i,j,k,id1,id2,icl, numcell, clust[2][5000];
 260   Int_t jd1,jd2, icell, cellcount;
 261   static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
 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   for(i=0; i < kNMX; i++)
 274     {
 275       fInfcl[0][i] = -1;
 276
 277       j  = iord1[i];
 278       id2 = j/kNDIMX;
 279       id1 = j-id2*kNDIMX;
 280
 281       if(edepcell[j] <= cutoff)
 282         {
 283           fInfocl[0][id1][id2] = -1;
 284         }
 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
 292   icl       = -1;
 293   cellcount = -1;
 294
 295   for(icell = 0; icell <= nmx1; icell++)
 296     {
 297       j  = iord1[icell];
 298       id2 = j/kNDIMX;
 299       id1 = j-id2*kNDIMX;
 300
 301       if(fInfocl[0][id1][id2] == 0 )
 302         {
 303           icl++;
 304           numcell = 0;
 305           cellcount++;
 306           fInfocl[0][id1][id2] = 1;
 307           fInfocl[1][id1][id2] = icl;
 308           fInfcl[0][cellcount] = icl;
 309           fInfcl[1][cellcount] = id1;
 310           fInfcl[2][cellcount] = id2;
 311
 312           clust[0][numcell] = id1;
 313           clust[1][numcell] = id2;
 314
 315           for(i = 1; i < 5000; i++)
 316             {
 317               clust[0][i]=0;
 318             }
 319           // ---------------------------------------------------------------
 320           // check for adc count in neib. cells. If ne 0 put it in this clust
 321           // ---------------------------------------------------------------
 322           for(i = 0; i < 6; i++)
 323             {
 324               jd1 = id1 + neibx[i];
 325               jd2 = id2 + neiby[i];
 326               if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
 327                   fInfocl[0][jd1][jd2] == 0)
 328                 {
 329                   numcell++;
 330                   fInfocl[0][jd1][jd2] = 2;
 331                   fInfocl[1][jd1][jd2] = icl;
 332                   clust[0][numcell]    = jd1;
 333                   clust[1][numcell]    = jd2;
 334                   cellcount++;
 335                   fInfcl[0][cellcount] = icl;
 336                   fInfcl[1][cellcount] = jd1;
 337                   fInfcl[2][cellcount] = jd2;
 338                 }
 339             }
 340           // ---------------------------------------------------------------
 341           // check adc count for neighbour's neighbours recursively and
 342           // if nonzero, add these to the cluster.
 343           // ---------------------------------------------------------------
 344           for(i = 1; i < 5000;i++)
 345             {
 346               if(clust[0][i] != 0)
 347                 {
 348                   id1 = clust[0][i];
 349                   id2 = clust[1][i];
 350                   for(j = 0; j < 6 ; j++)
 351                     {
 352                       jd1 = id1 + neibx[j];
 353                       jd2 = id2 + neiby[j];
 354                       if( (jd1 >= 0 && jd1 < kNDIMX) &&
 355                           (jd2 >= 0 && jd2 < kNDIMY) &&
 356                           fInfocl[0][jd1][jd2] == 0 )
 357                         {
 358                           fInfocl[0][jd1][jd2] = 2;
 359                           fInfocl[1][jd1][jd2] = icl;
 360                           numcell++;
 361                           clust[0][numcell]    = jd1;
 362                           clust[1][numcell]    = jd2;
 363                           cellcount++;
 364                           fInfcl[0][cellcount] = icl;
 365                           fInfcl[1][cellcount] = jd1;
 366                           fInfcl[2][cellcount] = jd2;
 367                         }
 368                     }
 369                 }
 370             }
 371         }
 372     }
 373   return cellcount;
 374 }
 375 // ------------------------------------------------------------------------ //
 376 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
 377 {
 378   // Does the refining of clusters
 379   // Takes the big patch and does gaussian fitting and
 380   // finds out the more refined clusters
 381   //
 382
 383
 384
 385   AliPMDcludata *pmdcludata = 0;
 386
 387   Int_t *cellCount = 0x0;
 388   Int_t **cellXY = 0x0;
 389   const Int_t kdim = 4500;
 390
 391   Int_t    i, j, k, i1, i2, id, icl,  itest,ihld, ig, nsupcl,clno;
 392   Int_t    t[kdim];
 393   Int_t    ncl[kdim], iord[kdim], lev1[kdim], lev2[kdim];
 394   Int_t    clxy[15];
 395   Float_t  clusdata[6];
 396   Double_t x1, y1, z1, x2, y2, z2, dist,rr,sum;
 397   Double_t x[kdim], y[kdim], z[kdim];
 398   Double_t xc[kdim], yc[kdim], zc[kdim], cells[kdim], rc[kdim];
 399
 400   // Initialisation
 401   for(i = 0; i<kdim; i++)
 402     {
 403       t[i]         = -1;
 404       ncl[i]       = -1;
 405       if (i < 6) clusdata[i] = 0.;
 406       if (i < 15) clxy[i] = 0;
 407     }
 408
 409   // clno counts the final clusters
 410   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
 411   // x, y and z store (x,y) coordinates of and energy deposited in a cell
 412   // xc, yc store (x,y) coordinates of the cluster center
 413   // zc stores the energy deposited in a cluster
 414   // rc is cluster radius
 415
 416   clno  = -1;
 417   nsupcl = -1;
 418
 419   for(i = 0; i <= incr; i++)
 420     {
 421       if(fInfcl[0][i] != nsupcl)
 422         {
 423           nsupcl++;
 424         }
 425       if (nsupcl > kdim)
 426         {
 427           AliWarning("RefClust: Too many superclusters!");
 428           nsupcl = kdim;
 429           break;
 430         }
 431       ncl[nsupcl]++;
 432     }
 433
 434   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
 435                   incr+1,nsupcl+1));
 436   id  = -1;
 437   icl = -1;
 438
 439   for(i = 0; i <= nsupcl; i++)
 440     {
 441       if(ncl[i] == 0)
 442         {
 443           id++;
 444           icl++;
 445           if (clno >= 5000)
 446             {
 447               AliWarning("RefClust: Too many clusters! more than 5000");
 448               return;
 449             }
 450           clno++;
 451           i1 = fInfcl[1][id];
 452           i2 = fInfcl[2][id];
 453
 454           Int_t i12 = i1 + i2*kNDIMX;
 455
 456           clusdata[0] = fCoord[0][i1][i2];
 457           clusdata[1] = fCoord[1][i1][i2];
 458           clusdata[2] = edepcell[i12];
 459           clusdata[3] = 1.;
 460           clusdata[4] = 0.5;
 461           clusdata[5] = 0.0;
 462
 463           clxy[0] = fCellTrNo[i1][i2];            //association
 464           for(Int_t icltr = 1; icltr < 15; icltr++)
 465             {
 466               clxy[icltr] = -1;
 467             }
 468           pmdcludata  = new AliPMDcludata(clusdata,clxy);
 469           fPMDclucont->Add(pmdcludata);
 470         }
 471       else if(ncl[i] == 1)
 472         {
 473           id++;
 474           icl++;
 475           if (clno >= 5000)
 476             {
 477               AliWarning("RefClust: Too many clusters! more than 5000");
 478               return;
 479             }
 480           clno++;
 481           i1   = fInfcl[1][id];
 482           i2   = fInfcl[2][id];
 483           Int_t i12 = i1 + i2*kNDIMX;
 484
 485           x1   = fCoord[0][i1][i2];
 486           y1   = fCoord[1][i1][i2];
 487           z1   = edepcell[i12];
 488           clxy[0] = fCellTrNo[i1][i2];    //asso
 489           id++;
 490           i1   = fInfcl[1][id];
 491           i2   = fInfcl[2][id];
 492
 493           Int_t i22 = i1 + i2*kNDIMX;
 494           x2   = fCoord[0][i1][i2];
 495           y2   = fCoord[1][i1][i2];
 496           z2   = edepcell[i22];
 497           clxy[1] = fCellTrNo[i1][i2];            //asso
 498           for(Int_t icltr = 2; icltr < 15; icltr++)
 499             {
 500               clxy[icltr] = -1;
 501             }
 502
 503           clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
 504           clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
 505           clusdata[2] = z1+z2;
 506           clusdata[3] = 2.;
 507           clusdata[4] = 0.5;
 508           clusdata[5] = 0.0;
 509           pmdcludata  = new AliPMDcludata(clusdata,clxy);
 510           fPMDclucont->Add(pmdcludata);
 511         }
 512       else
 513         {
 514           id++;
 515           iord[0] = 0;
 516           // super-cluster of more than two cells - broken up into smaller
 517           // clusters gaussian centers computed. (peaks separated by > 1 cell)
 518           // Begin from cell having largest energy deposited This is first
 519           // cluster center
 520           i1      = fInfcl[1][id];
 521           i2      = fInfcl[2][id];
 522           Int_t i12 = i1 + i2*kNDIMX;
 523
 524           x[0]    = fCoord[0][i1][i2];
 525           y[0]    = fCoord[1][i1][i2];
 526           z[0]    = edepcell[i12];
 527
 528           t[0] = fCellTrNo[i1][i2];       //asso
 529
 530           iord[0] = 0;
 531           for(j = 1; j <= ncl[i]; j++)
 532             {
 533               id++;
 534               i1      = fInfcl[1][id];
 535               i2      = fInfcl[2][id];
 536               Int_t i12 = i1 + i2*kNDIMX;
 537
 538               iord[j] = j;
 539               x[j]    = fCoord[0][i1][i2];
 540               y[j]    = fCoord[1][i1][i2];
 541               z[j]    = edepcell[i12];
 542
 543               t[j]    = fCellTrNo[i1][i2];            //asso
 544             }
 545
 546           // arranging cells within supercluster in decreasing order
 547
 548           for(j = 1;j <= ncl[i]; j++)
 549             {
 550               itest = 0;
 551               ihld  = iord[j];
 552               for(i1 = 0; i1 < j; i1++)
 553                 {
 554                   if(itest == 0 && z[iord[i1]] < z[ihld])
 555                     {
 556                       itest = 1;
 557                       for(i2 = j-1; i2 >= i1; i2--)
 558                         {
 559                           iord[i2+1] = iord[i2];
 560                         }
 561                       iord[i1] = ihld;
 562                     }
 563                 }
 564             }
 565           // compute the number of Gaussians and their centers ( first
 566           // guess )
 567           // centers must be separated by cells having smaller ener. dep.
 568           // neighbouring centers should be either strong or well-separated
 569           ig = 0;
 570           xc[ig] = x[iord[0]];
 571           yc[ig] = y[iord[0]];
 572           zc[ig] = z[iord[0]];
 573           for(j = 1; j <= ncl[i]; j++)
 574             {
 575               itest = -1;
 576               x1    = x[iord[j]];
 577               y1    = y[iord[j]];
 578               for(k = 0; k <= ig; k++)
 579                 {
 580                   x2 = xc[k];
 581                   y2 = yc[k];
 582                   rr = Distance(x1,y1,x2,y2);
 583                   if(rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)itest++;
 584                   if(rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)itest++;
 585                   if( rr >= 2.1)itest++;
 586                 }
 587               if(itest == ig)
 588                 {
 589                   ig++;
 590                   xc[ig] = x1;
 591                   yc[ig] = y1;
 592                   zc[ig] = z[iord[j]];
 593                 }
 594             }
 595           GaussFit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]);
 596           icl += ig+1;
 597           // compute the number of cells belonging to each cluster.
 598           // cell is shared between several clusters ( if they are equidistant
 599           // from it ) in the ratio of cluster energy deposition
 600
 601           Int_t jj = 15;
 602           cellCount = new Int_t [ig+1];
 603           cellXY = new Int_t *[jj];
 604           for(Int_t ij = 0; ij < 15; ij++) cellXY[ij] = new Int_t [ig+1];
 605
 606           for(j = 0; j <= ig; j++)
 607             {
 608               cellCount[j] = 0;
 609               cells[j]     = 0.;
 610             }
 611
 612           if(ig > 0)
 613             {
 614               for(j = 0; j <= ncl[i]; j++)
 615                 {
 616                   lev1[j] = 0;
 617                   lev2[j] = 0;
 618                   for(k = 0; k <= ig; k++)
 619                     {
 620                       dist = Distance(x[j], y[j], xc[k], yc[k]);
 621                       if(dist < TMath::Sqrt(3.) )
 622                         {
 623                           //asso
 624                           if (cellCount[k] < 15)
 625                             {
 626                               cellXY[cellCount[k]][k] = t[j];
 627                             }
 628                           cellCount[k]++;
 629                           //
 630                           lev1[0]++;
 631                           i1       = lev1[0];
 632                           lev1[i1] = k;
 633                         }
 634                       else
 635                         {
 636                           if(dist < 2.1)
 637                             {
 638                               lev2[0]++;
 639                               i1       = lev2[0];
 640                               lev2[i1] = k;
 641                             }
 642                         }
 643                     }
 644                   if(lev1[0] != 0)
 645                     {
 646                       if(lev1[0] == 1)
 647                         {
 648                           cells[lev1[1]]++;
 649                         }
 650                       else
 651                         {
 652                           sum=0.;
 653                           for(k = 1; k <= lev1[0]; k++)
 654                             {
 655                               sum  += zc[lev1[k]];
 656                             }
 657                           for(k = 1; k <= lev1[0]; k++)
 658                             {
 659                               cells[lev1[k]] += zc[lev1[k]]/sum;
 660                             }
 661                         }
 662                     }
 663                   else
 664                     {
 665                       if(lev2[0] == 0)
 666                         {
 667                           cells[lev2[1]]++;
 668                         }
 669                       else
 670                         {
 671                           sum=0.;
 672                           for( k = 1; k <= lev2[0]; k++)
 673                             {
 674                               sum += zc[lev2[k]];
 675                             }
 676                           for(k = 1; k <= lev2[0]; k++)
 677                             {
 678                               cells[lev2[k]] +=  zc[lev2[k]]/sum;
 679                             }
 680                         }
 681                     }
 682                 }
 683             }
 684
 685           // zero rest of the cell array
 686           //asso
 687           for( k = 0; k <= ig; k++)
 688             {
 689               for(Int_t icltr = cellCount[k]; icltr < 15; icltr++)
 690                 {
 691                   cellXY[icltr][k] = -1;
 692                 }
 693             }
 694           //
 695
 696           for(j = 0; j <= ig; j++)
 697             {
 698               clno++;
 699               if (clno >= 5000)
 700                 {
 701                   AliWarning("RefClust: Too many clusters! more than 5000");
 702                   return;
 703                 }
 704               clusdata[0] = xc[j];
 705               clusdata[1] = yc[j];
 706               clusdata[2] = zc[j];
 707               clusdata[4] = rc[j];
 708               clusdata[5] = 0.0;
 709               if(ig == 0)
 710                 {
 711                   clusdata[3] = ncl[i];
 712                 }
 713               else
 714                 {
 715                   clusdata[3] = cells[j];
 716                 }
 717
 718
 719               for (Int_t ii=0; ii < 15; ii++)
 720                 {
 721                   clxy[ii] = cellXY[ii][j];
 722                 }
 723               pmdcludata = new AliPMDcludata(clusdata,clxy);
 724               fPMDclucont->Add(pmdcludata);
 725             }
 726           delete [] cellCount;
 727           for(Int_t jj = 0; jj < 15; jj++)delete [] cellXY[jj];
 728           delete [] cellXY;
 729         }
 730     }
 731 }
 732 // ------------------------------------------------------------------------ //
 733 void AliPMDClusteringV1::GaussFit(Int_t ncell, Int_t nclust, Double_t &x,
 734                                   Double_t &y ,Double_t &z, Double_t &xc,
 735                                   Double_t &yc, Double_t &zc, Double_t &rc)
 736 {
 737   // Does gaussian fitting
 738   //
 739
 740   const Int_t kdim = 4500;
 741   Int_t i, j, i1, i2, novar, idd, jj;
 742   Int_t neib[kdim][50];
 743
 744   Double_t sum, dx, dy, str, str1, aint, sum1, rr, dum;
 745   Double_t x1, x2, y1, y2;
 746   Double_t xx[kdim], yy[kdim], zz[kdim], xxc[kdim], yyc[kdim];
 747   Double_t a[kdim], b[kdim], c[kdim], d[kdim], ha[kdim], hb[kdim];
 748   Double_t hc[kdim], hd[kdim], zzc[kdim], rrc[kdim];
 749
 750   TRandom rnd;
 751
 752   str   = 0.;
 753   str1  = 0.;
 754   rr    = 0.3;
 755   novar = 0;
 756   j     = 0;
 757
 758   for(i = 0; i <= ncell; i++)
 759     {
 760       xx[i] = *(&x+i);
 761       yy[i] = *(&y+i);
 762       zz[i] = *(&z+i);
 763       str  += zz[i];
 764     }
 765   for(i=0; i<=nclust; i++)
 766     {
 767       xxc[i] = *(&xc+i);
 768       yyc[i] = *(&yc+i);
 769       zzc[i] = *(&zc+i);
 770       str1  += zzc[i];
 771       rrc[i] = 0.5;
 772     }
 773   for(i = 0; i <= nclust; i++)
 774     {
 775       zzc[i] = str/str1*zzc[i];
 776       ha[i]  = xxc[i];
 777       hb[i]  = yyc[i];
 778       hc[i]  = zzc[i];
 779       hd[i]  = rrc[i];
 780       x1     = xxc[i];
 781       y1     = yyc[i];
 782     }
 783
 784   for(i = 0; i <= ncell; i++)
 785     {
 786       idd = 0;
 787       x1  = xx[i];
 788       y1  = yy[i];
 789       for(j = 0; j <= nclust; j++)
 790         {
 791           x2 = xxc[j];
 792           y2 = yyc[j];
 793           if(Distance(x1,y1,x2,y2) <= 3.)
 794             {
 795               idd++;
 796               neib[i][idd] = j;
 797             }
 798         }
 799       neib[i][0] = idd;
 800     }
 801   sum = 0.;
 802   for(i1 = 0; i1 <= ncell; i1++)
 803     {
 804       aint = 0.;
 805       idd = neib[i1][0];
 806       for(i2 = 1; i2 <= idd; i2++)
 807         {
 808           jj    = neib[i1][i2];
 809           dx    = xx[i1] - xxc[jj];
 810           dy    = yy[i1] - yyc[jj];
 811           dum   = rrc[j]*rrc[jj] + rr*rr;
 812           aint += exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum;
 813         }
 814       sum += (aint - zz[i1])*(aint - zz[i1])/str;
 815     }
 816   str1 = 0.;
 817
 818   for(i = 0; i <= nclust; i++)
 819     {
 820       a[i]  = xxc[i] + 0.6*(rnd.Uniform() - 0.5);
 821       b[i]  = yyc[i] + 0.6*(rnd.Uniform() - 0.5);
 822       c[i]  = zzc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.2);
 823       str1 += zzc[i];
 824       d[i]  = rrc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.1);
 825
 826       if(d[i] < 0.25)
 827         {
 828           d[i]=0.25;
 829         }
 830     }
 831   for(i = 0; i <= nclust; i++)
 832     {
 833       c[i] = c[i]*str/str1;
 834     }
 835   sum1=0.;
 836
 837   for(i1 = 0; i1 <= ncell; i1++)
 838     {
 839       aint = 0.;
 840       idd = neib[i1][0];
 841       for(i2 = 1; i2 <= idd; i2++)
 842         {
 843           jj    = neib[i1][i2];
 844           dx    = xx[i1] - a[jj];
 845           dy    = yy[i1] - b[jj];
 846           dum   = d[jj]*d[jj]+rr*rr;
 847           aint += exp(-(dx*dx+dy*dy)/dum)*c[i2]*rr*rr/dum;
 848         }
 849       sum1 += (aint - zz[i1])*(aint - zz[i1])/str;
 850     }
 851
 852     if(sum1 < sum)
 853       {
 854         for(i2 = 0; i2 <= nclust; i2++)
 855         {
 856           xxc[i2] = a[i2];
 857           yyc[i2] = b[i2];
 858           zzc[i2] = c[i2];
 859           rrc[i2] = d[i2];
 860           sum     = sum1;
 861         }
 862       }
 863     for(j = 0; j <= nclust; j++)
 864       {
 865         *(&xc+j) = xxc[j];
 866         *(&yc+j) = yyc[j];
 867         *(&zc+j) = zzc[j];
 868         *(&rc+j) = rrc[j];
 869       }
 870 }
 871 // ------------------------------------------------------------------------ //
 872 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
 873                                       Double_t x2, Double_t y2)
 874 {
 875   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
 876 }
 877 // ------------------------------------------------------------------------ //
 878 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
 879 {
 880   fCutoff = decut;
 881 }
 882 // ------------------------------------------------------------------------ //