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           Int_t ij = i + j*kNDIMX;
 148
 149           if (ismn < 12)
 150             {
 151               edepcell[ij]    = celladc[jd][id];
 152               fCellTrNo[i][j] = jd*10000+id;  // for association
 153             }
 154           else if (ismn >= 12 && ismn <= 23)
 155             {
 156               edepcell[ij]    = celladc[id][jd];
 157               fCellTrNo[i][j] = id*10000+jd;  // for association
 158             }
 159
 160         }
 161     }
 162
 163   Int_t iord1[kNMX];
 164   TMath::Sort(kNMX,edepcell,iord1);// order the data
 165   cutoff = fCutoff;                // cutoff to discard cells
 166   ave  = 0.;
 167   nmx1 = -1;
 168   for(i = 0;i < kNMX; i++)
 169     {
 170       if(edepcell[i] > 0.)
 171         {
 172           ave += edepcell[i];
 173         }
 174       if(edepcell[i] > cutoff )
 175         {
 176           nmx1++;
 177         }
 178     }
 179
 180   AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
 181
 182   if (nmx1 == 0) nmx1 = 1;
 183   ave = ave/nmx1;
 184
 185   AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
 186                   kNMX,ave));
 187
 188   incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
 189
 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   AliPMDcludata *pmdcludata = 0;
 384
 385   const Int_t kdim = 4500;
 386
 387   Int_t i, j, k, i1, i2, id, icl,  itest,ihld, ig, nsupcl,clno;
 388   Int_t lev1[20], lev2[20];
 389   Int_t ncl[kdim], iord[kdim], t[kdim];
 390   Int_t clxy[15];
 391
 392   Int_t *cellCount;
 393   Int_t **cellXY;
 394
 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
 401   // Initialisation
 402   for(i = 0; i<kdim; i++)
 403     {
 404       t[i]         = -1;
 405       ncl[i]       = -1;
 406       if (i < 6) clusdata[i] = 0.;
 407       if (i < 15) clxy[i] = 0;
 408     }
 409
 410   // clno counts the final clusters
 411   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
 412   // x, y and z store (x,y) coordinates of and energy deposited in a cell
 413   // xc, yc store (x,y) coordinates of the cluster center
 414   // zc stores the energy deposited in a cluster
 415   // rc is cluster radius
 416
 417   clno  = -1;
 418   nsupcl = -1;
 419
 420   for(i = 0; i <= incr; i++)
 421     {
 422       if(fInfcl[0][i] != nsupcl)
 423         {
 424           nsupcl++;
 425         }
 426       if (nsupcl > kdim)
 427         {
 428           AliWarning("RefClust: Too many superclusters!");
 429           nsupcl = kdim;
 430           break;
 431         }
 432       ncl[nsupcl]++;
 433     }
 434
 435   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
 436                   incr+1,nsupcl+1));
 437   id  = -1;
 438   icl = -1;
 439
 440   for(i = 0; i <= nsupcl; i++)
 441     {
 442       if(ncl[i] == 0)
 443         {
 444           id++;
 445           icl++;
 446           if (clno >= 5000)
 447             {
 448               AliWarning("RefClust: Too many clusters! more than 5000");
 449               return;
 450             }
 451           clno++;
 452           i1 = fInfcl[1][id];
 453           i2 = fInfcl[2][id];
 454
 455           Int_t i12 = i1 + i2*kNDIMX;
 456
 457           clusdata[0] = fCoord[0][i1][i2];
 458           clusdata[1] = fCoord[1][i1][i2];
 459           clusdata[2] = edepcell[i12];
 460           clusdata[3] = 1.;
 461           clusdata[4] = 0.5;
 462           clusdata[5] = 0.0;
 463
 464           clxy[0] = fCellTrNo[i1][i2];            //association
 465           for(Int_t icltr = 1; icltr < 15; icltr++)
 466             {
 467               clxy[icltr] = -1;
 468             }
 469           pmdcludata  = new AliPMDcludata(clusdata,clxy);
 470           fPMDclucont->Add(pmdcludata);
 471         }
 472       else if(ncl[i] == 1)
 473         {
 474           id++;
 475           icl++;
 476           if (clno >= 5000)
 477             {
 478               AliWarning("RefClust: Too many clusters! more than 5000");
 479               return;
 480             }
 481           clno++;
 482           i1   = fInfcl[1][id];
 483           i2   = fInfcl[2][id];
 484           Int_t i12 = i1 + i2*kNDIMX;
 485
 486           x1   = fCoord[0][i1][i2];
 487           y1   = fCoord[1][i1][i2];
 488           z1   = edepcell[i12];
 489           clxy[0] = fCellTrNo[i1][i2];    //asso
 490           id++;
 491           i1   = fInfcl[1][id];
 492           i2   = fInfcl[2][id];
 493
 494           Int_t i22 = i1 + i2*kNDIMX;
 495           x2   = fCoord[0][i1][i2];
 496           y2   = fCoord[1][i1][i2];
 497           z2   = edepcell[i22];
 498           clxy[1] = fCellTrNo[i1][i2];            //asso
 499           for(Int_t icltr = 2; icltr < 15; icltr++)
 500             {
 501               clxy[icltr] = -1;
 502             }
 503
 504           clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
 505           clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
 506           clusdata[2] = z1+z2;
 507           clusdata[3] = 2.;
 508           clusdata[4] = 0.5;
 509           clusdata[5] = 0.0;
 510           pmdcludata  = new AliPMDcludata(clusdata,clxy);
 511           fPMDclucont->Add(pmdcludata);
 512         }
 513       else
 514         {
 515           id++;
 516           iord[0] = 0;
 517           // super-cluster of more than two cells - broken up into smaller
 518           // clusters gaussian centers computed. (peaks separated by > 1 cell)
 519           // Begin from cell having largest energy deposited This is first
 520           // cluster center
 521           i1      = fInfcl[1][id];
 522           i2      = fInfcl[2][id];
 523           Int_t i12 = i1 + i2*kNDIMX;
 524
 525           x[0]    = fCoord[0][i1][i2];
 526           y[0]    = fCoord[1][i1][i2];
 527           z[0]    = edepcell[i12];
 528
 529           t[0] = fCellTrNo[i1][i2];       //asso
 530
 531           iord[0] = 0;
 532           for(j = 1; j <= ncl[i]; j++)
 533             {
 534               id++;
 535               i1      = fInfcl[1][id];
 536               i2      = fInfcl[2][id];
 537               Int_t i12 = i1 + i2*kNDIMX;
 538
 539               iord[j] = j;
 540               x[j]    = fCoord[0][i1][i2];
 541               y[j]    = fCoord[1][i1][i2];
 542               z[j]    = edepcell[i12];
 543
 544               t[j]    = fCellTrNo[i1][i2];            //asso
 545             }
 546
 547           // arranging cells within supercluster in decreasing order
 548
 549           for(j = 1;j <= ncl[i]; j++)
 550             {
 551               itest = 0;
 552               ihld  = iord[j];
 553               for(i1 = 0; i1 < j; i1++)
 554                 {
 555                   if(itest == 0 && z[iord[i1]] < z[ihld])
 556                     {
 557                       itest = 1;
 558                       for(i2 = j-1; i2 >= i1; i2--)
 559                         {
 560                           iord[i2+1] = iord[i2];
 561                         }
 562                       iord[i1] = ihld;
 563                     }
 564                 }
 565             }
 566           // compute the number of Gaussians and their centers ( first
 567           // guess )
 568           // centers must be separated by cells having smaller ener. dep.
 569           // neighbouring centers should be either strong or well-separated
 570           ig = 0;
 571           xc[ig] = x[iord[0]];
 572           yc[ig] = y[iord[0]];
 573           zc[ig] = z[iord[0]];
 574           for(j = 1; j <= ncl[i]; j++)
 575             {
 576               itest = -1;
 577               x1    = x[iord[j]];
 578               y1    = y[iord[j]];
 579               for(k = 0; k <= ig; k++)
 580                 {
 581                   x2 = xc[k];
 582                   y2 = yc[k];
 583                   rr = Distance(x1,y1,x2,y2);
 584                   if(rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)itest++;
 585                   if(rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)itest++;
 586                   if( rr >= 2.1)itest++;
 587                 }
 588               if(itest == ig)
 589                 {
 590                   ig++;
 591                   xc[ig] = x1;
 592                   yc[ig] = y1;
 593                   zc[ig] = z[iord[j]];
 594                 }
 595             }
 596           GaussFit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]);
 597           icl += ig+1;
 598           // compute the number of cells belonging to each cluster.
 599           // cell is shared between several clusters ( if they are equidistant
 600           // from it ) in the ratio of cluster energy deposition
 601
 602           Int_t jj = 15;
 603           cellCount = new Int_t [ig+1];
 604           cellXY = new Int_t *[jj];
 605           for(Int_t ij = 0; ij < 15; ij++) cellXY[ij] = new Int_t [ig+1];
 606
 607           for(j = 0; j <= ig; j++)
 608             {
 609               cellCount[j] = 0;
 610               cells[j]     = 0.;
 611             }
 612
 613           if(ig > 0)
 614             {
 615               for(j = 0; j <= ncl[i]; j++)
 616                 {
 617                   lev1[j] = 0;
 618                   lev2[j] = 0;
 619                   for(k = 0; k <= ig; k++)
 620                     {
 621                       dist = Distance(x[j], y[j], xc[k], yc[k]);
 622                       if(dist < TMath::Sqrt(3.) )
 623                         {
 624                           //asso
 625                           if (cellCount[k] < 15)
 626                             {
 627                               cellXY[cellCount[k]][k] = t[j];
 628                             }
 629                           cellCount[k]++;
 630                           //
 631                           lev1[0]++;
 632                           i1       = lev1[0];
 633                           lev1[i1] = k;
 634                         }
 635                       else
 636                         {
 637                           if(dist < 2.1)
 638                             {
 639                               lev2[0]++;
 640                               i1       = lev2[0];
 641                               lev2[i1] = k;
 642                             }
 643                         }
 644                     }
 645                   if(lev1[0] != 0)
 646                     {
 647                       if(lev1[0] == 1)
 648                         {
 649                           cells[lev1[1]]++;
 650                         }
 651                       else
 652                         {
 653                           sum=0.;
 654                           for(k = 1; k <= lev1[0]; k++)
 655                             {
 656                               sum  += zc[lev1[k]];
 657                             }
 658                           for(k = 1; k <= lev1[0]; k++)
 659                             {
 660                               cells[lev1[k]] += zc[lev1[k]]/sum;
 661                             }
 662                         }
 663                     }
 664                   else
 665                     {
 666                       if(lev2[0] == 0)
 667                         {
 668                           cells[lev2[1]]++;
 669                         }
 670                       else
 671                         {
 672                           sum=0.;
 673                           for( k = 1; k <= lev2[0]; k++)
 674                             {
 675                               sum += zc[lev2[k]];
 676                             }
 677                           for(k = 1; k <= lev2[0]; k++)
 678                             {
 679                               cells[lev2[k]] +=  zc[lev2[k]]/sum;
 680                             }
 681                         }
 682                     }
 683                 }
 684             }
 685
 686           // zero rest of the cell array
 687           //asso
 688           for( k = 0; k <= ig; k++)
 689             {
 690               for(Int_t icltr = cellCount[k]; icltr < 15; icltr++)
 691                 {
 692                   cellXY[icltr][k] = -1;
 693                 }
 694             }
 695           //
 696
 697           for(j = 0; j <= ig; j++)
 698             {
 699               clno++;
 700               if (clno >= 5000)
 701                 {
 702                   AliWarning("RefClust: Too many clusters! more than 5000");
 703                   return;
 704                 }
 705               clusdata[0] = xc[j];
 706               clusdata[1] = yc[j];
 707               clusdata[2] = zc[j];
 708               clusdata[4] = rc[j];
 709               clusdata[5] = 0.0;
 710               if(ig == 0)
 711                 {
 712                   clusdata[3] = ncl[i];
 713                 }
 714               else
 715                 {
 716                   clusdata[3] = cells[j];
 717                 }
 718
 719
 720               for (Int_t ii=0; ii < 15; ii++)
 721                 {
 722                   clxy[ii] = cellXY[ii][j];
 723                 }
 724               pmdcludata = new AliPMDcludata(clusdata,clxy);
 725               fPMDclucont->Add(pmdcludata);
 726             }
 727           delete [] cellCount;
 728           for(Int_t jj = 0; jj < 15; jj++)delete [] cellXY[jj];
 729           delete [] cellXY;
 730         }
 731     }
 732 }
 733 // ------------------------------------------------------------------------ //
 734 void AliPMDClusteringV1::GaussFit(Int_t ncell, Int_t nclust, Double_t &x,
 735                                   Double_t &y ,Double_t &z, Double_t &xc,
 736                                   Double_t &yc, Double_t &zc, Double_t &rc)
 737 {
 738   // Does gaussian fitting
 739   //
 740
 741   const Int_t kdim = 4500;
 742   Int_t i, j, i1, i2, novar, idd, jj;
 743   Int_t neib[kdim][50];
 744
 745   Double_t sum, dx, dy, str, str1, aint, sum1, rr, dum;
 746   Double_t x1, x2, y1, y2;
 747   Double_t xx[kdim], yy[kdim], zz[kdim], xxc[kdim], yyc[kdim];
 748   Double_t a[kdim], b[kdim], c[kdim], d[kdim], ha[kdim], hb[kdim];
 749   Double_t hc[kdim], hd[kdim], zzc[kdim], rrc[kdim];
 750
 751   TRandom rnd;
 752
 753   str   = 0.;
 754   str1  = 0.;
 755   rr    = 0.3;
 756   novar = 0;
 757   j     = 0;
 758
 759   for(i = 0; i <= ncell; i++)
 760     {
 761       xx[i] = *(&x+i);
 762       yy[i] = *(&y+i);
 763       zz[i] = *(&z+i);
 764       str  += zz[i];
 765     }
 766   for(i=0; i<=nclust; i++)
 767     {
 768       xxc[i] = *(&xc+i);
 769       yyc[i] = *(&yc+i);
 770       zzc[i] = *(&zc+i);
 771       str1  += zzc[i];
 772       rrc[i] = 0.5;
 773     }
 774   for(i = 0; i <= nclust; i++)
 775     {
 776       zzc[i] = str/str1*zzc[i];
 777       ha[i]  = xxc[i];
 778       hb[i]  = yyc[i];
 779       hc[i]  = zzc[i];
 780       hd[i]  = rrc[i];
 781       x1     = xxc[i];
 782       y1     = yyc[i];
 783     }
 784
 785   for(i = 0; i <= ncell; i++)
 786     {
 787       idd = 0;
 788       x1  = xx[i];
 789       y1  = yy[i];
 790       for(j = 0; j <= nclust; j++)
 791         {
 792           x2 = xxc[j];
 793           y2 = yyc[j];
 794           if(Distance(x1,y1,x2,y2) <= 3.)
 795             {
 796               idd++;
 797               neib[i][idd] = j;
 798             }
 799         }
 800       neib[i][0] = idd;
 801     }
 802   sum = 0.;
 803   for(i1 = 0; i1 <= ncell; i1++)
 804     {
 805       aint = 0.;
 806       idd = neib[i1][0];
 807       for(i2 = 1; i2 <= idd; i2++)
 808         {
 809           jj    = neib[i1][i2];
 810           dx    = xx[i1] - xxc[jj];
 811           dy    = yy[i1] - yyc[jj];
 812           dum   = rrc[j]*rrc[jj] + rr*rr;
 813           aint += exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum;
 814         }
 815       sum += (aint - zz[i1])*(aint - zz[i1])/str;
 816     }
 817   str1 = 0.;
 818
 819   for(i = 0; i <= nclust; i++)
 820     {
 821       a[i]  = xxc[i] + 0.6*(rnd.Uniform() - 0.5);
 822       b[i]  = yyc[i] + 0.6*(rnd.Uniform() - 0.5);
 823       c[i]  = zzc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.2);
 824       str1 += zzc[i];
 825       d[i]  = rrc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.1);
 826
 827       if(d[i] < 0.25)
 828         {
 829           d[i]=0.25;
 830         }
 831     }
 832   for(i = 0; i <= nclust; i++)
 833     {
 834       c[i] = c[i]*str/str1;
 835     }
 836   sum1=0.;
 837
 838   for(i1 = 0; i1 <= ncell; i1++)
 839     {
 840       aint = 0.;
 841       idd = neib[i1][0];
 842       for(i2 = 1; i2 <= idd; i2++)
 843         {
 844           jj    = neib[i1][i2];
 845           dx    = xx[i1] - a[jj];
 846           dy    = yy[i1] - b[jj];
 847           dum   = d[jj]*d[jj]+rr*rr;
 848           aint += exp(-(dx*dx+dy*dy)/dum)*c[i2]*rr*rr/dum;
 849         }
 850       sum1 += (aint - zz[i1])*(aint - zz[i1])/str;
 851     }
 852
 853     if(sum1 < sum)
 854       {
 855         for(i2 = 0; i2 <= nclust; i2++)
 856         {
 857           xxc[i2] = a[i2];
 858           yyc[i2] = b[i2];
 859           zzc[i2] = c[i2];
 860           rrc[i2] = d[i2];
 861           sum     = sum1;
 862         }
 863       }
 864     for(j = 0; j <= nclust; j++)
 865       {
 866         *(&xc+j) = xxc[j];
 867         *(&yc+j) = yyc[j];
 868         *(&zc+j) = zzc[j];
 869         *(&rc+j) = rrc[j];
 870       }
 871 }
 872 // ------------------------------------------------------------------------ //
 873 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
 874                                       Double_t x2, Double_t y2)
 875 {
 876   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
 877 }
 878 // ------------------------------------------------------------------------ //
 879 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
 880 {
 881   fCutoff = decut;
 882 }
 883 // ------------------------------------------------------------------------ //