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
  58 ClassImp(AliPMDClusteringV1)
  59
  60 const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254;  // sqrt(3.)/2.
  61
  62 AliPMDClusteringV1::AliPMDClusteringV1():
  63   fPMDclucont(new TObjArray()),
  64   fCutoff(0.0)
  65 {
  66   for(Int_t i = 0; i < kNDIMX; i++)
  67     {
  68       for(Int_t j = 0; j < kNDIMY; j++)
  69         {
  70           fCoord[0][i][j] = i+j/2.;
  71           fCoord[1][i][j] = fgkSqroot3by2*j;
  72         }
  73     }
  74 }
  75 // ------------------------------------------------------------------------ //
  76 AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
  77   AliPMDClustering(pmdclv1),
  78   fPMDclucont(0),
  79   fCutoff(0)
  80 {
  81   // copy constructor
  82   AliError("Copy constructor not allowed ");
  83
  84 }
  85 // ------------------------------------------------------------------------ //
  86 AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
  87 {
  88   // copy constructor
  89   AliError("Assignment operator not allowed ");
  90   return *this;
  91 }
  92 // ------------------------------------------------------------------------ //
  93 AliPMDClusteringV1::~AliPMDClusteringV1()
  94 {
  95   delete fPMDclucont;
  96 }
  97 // ------------------------------------------------------------------------ //
  98 void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
  99                                  Double_t celladc[48][96], TObjArray *pmdcont)
 100 {
 101   // main function to call other necessary functions to do clustering
 102   //
 103
 104   AliPMDcluster *pmdcl = 0;
 105
 106   const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
 107   const Int_t kNmaxCell   = 19;     // # of cells surrounding a cluster center
 108
 109   Int_t    i,  j, nmx1, incr, id, jd;
 110   Int_t    celldataX[kNmaxCell], celldataY[kNmaxCell];
 111   Float_t  clusdata[6];
 112   Double_t cutoff, ave;
 113   Double_t edepcell[kNMX];
 114
 115
 116   Double_t *cellenergy = new Double_t [11424];
 117
 118
 119   // ndimXr and ndimYr are different because of different module size
 120
 121   Int_t ndimXr = 0;
 122   Int_t ndimYr = 0;
 123
 124   if (ismn < 12)
 125     {
 126       ndimXr = 96;
 127       ndimYr = 48;
 128     }
 129   else if (ismn >= 12 && ismn <= 23)
 130     {
 131       ndimXr = 48;
 132       ndimYr = 96;
 133     }
 134
 135   for (i =0; i < 11424; i++)
 136   {
 137       cellenergy[i] = 0.;
 138   }
 139
 140
 141   Int_t kk = 0;
 142   for (i = 0; i < kNDIMX; i++)
 143     {
 144       for (j = 0; j < kNDIMY; j++)
 145         {
 146           edepcell[kk] = 0.;
 147           kk++;
 148         }
 149     }
 150
 151   for (id = 0; id < ndimXr; id++)
 152     {
 153       for (jd = 0; jd < ndimYr; jd++)
 154         {
 155           j = jd;
 156           i = id+(ndimYr/2-1)-(jd/2);
 157
 158           Int_t ij = i + j*kNDIMX;
 159
 160           if (ismn < 12)
 161             {
 162               cellenergy[ij]    = celladc[jd][id];
 163             }
 164           else if (ismn >= 12 && ismn <= 23)
 165             {
 166               cellenergy[ij]    = celladc[id][jd];
 167             }
 168         }
 169     }
 170
 171   for (i = 0; i < kNMX; i++)
 172   {
 173     edepcell[i] = cellenergy[i];
 174   }
 175
 176   delete [] cellenergy;
 177
 178   Int_t iord1[kNMX];
 179   TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
 180   cutoff = fCutoff;                       // cutoff to discard cells
 181   ave  = 0.;
 182   nmx1 = -1;
 183   for(i = 0;i < kNMX; i++)
 184     {
 185       if(edepcell[i] > 0.)
 186         {
 187           ave += edepcell[i];
 188         }
 189       if(edepcell[i] > cutoff )
 190         {
 191           nmx1++;
 192         }
 193     }
 194
 195   AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
 196
 197   if (nmx1 == 0) nmx1 = 1;
 198   ave = ave/nmx1;
 199   AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
 200                   kNMX,ave));
 201
 202   incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
 203   RefClust(incr,edepcell);
 204   Int_t nentries1 = fPMDclucont->GetEntries();
 205   AliDebug(1,Form("Detector Plane = %d  Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
 206   AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
 207
 208   for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 209     {
 210       AliPMDcludata *pmdcludata =
 211         (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
 212       Float_t cluXC    = pmdcludata->GetClusX();
 213       Float_t cluYC    = pmdcludata->GetClusY();
 214       Float_t cluADC   = pmdcludata->GetClusADC();
 215       Float_t cluCELLS = pmdcludata->GetClusCells();
 216       Float_t cluSIGX  = pmdcludata->GetClusSigmaX();
 217       Float_t cluSIGY  = pmdcludata->GetClusSigmaY();
 218
 219       Float_t cluY0    = ktwobysqrt3*cluYC;
 220       Float_t cluX0    = cluXC - cluY0/2.;
 221
 222       //
 223       // Cluster X centroid is back transformed
 224       //
 225
 226       if (ismn < 12)
 227         {
 228           clusdata[0] = cluX0 - (24-1) + cluY0/2.;
 229         }
 230       else if ( ismn >= 12 && ismn <= 23)
 231         {
 232           clusdata[0] = cluX0 - (48-1) + cluY0/2.;
 233         }
 234
 235
 236       clusdata[1]     = cluY0;
 237       clusdata[2]     = cluADC;
 238       clusdata[3]     = cluCELLS;
 239       clusdata[4]     = cluSIGX;
 240       clusdata[5]     = cluSIGY;
 241
 242       //
 243       // Cells associated with a cluster
 244       //
 245
 246       for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
 247         {
 248           Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
 249           Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
 250
 251           if (ismn < 12)
 252             {
 253               celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
 254             }
 255           else if (ismn >= 12 && ismn <= 23)
 256             {
 257               celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
 258             }
 259
 260           celldataY[ihit] = cellcol;
 261         }
 262
 263
 264       pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
 265       pmdcont->Add(pmdcl);
 266     }
 267
 268   fPMDclucont->Delete();
 269 }
 270 // ------------------------------------------------------------------------ //
 271 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
 272                                   Int_t iord1[], Double_t edepcell[])
 273 {
 274   // Does crude clustering
 275   // Finds out only the big patch by just searching the
 276   // connected cells
 277   //
 278   const Int_t kndim = 4609;
 279   Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
 280   Int_t jd1,jd2, icell, cellcount;
 281   static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
 282
 283   AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
 284
 285   for (j = 0; j < kNDIMX; j++)
 286     {
 287       for(k = 0; k < kNDIMY; k++)
 288         {
 289           fInfocl[0][j][k] = 0;
 290           fInfocl[1][j][k] = 0;
 291         }
 292     }
 293   for(i=0; i < kNMX; i++)
 294     {
 295       fInfcl[0][i] = -1;
 296
 297       j  = iord1[i];
 298       id2 = j/kNDIMX;
 299       id1 = j-id2*kNDIMX;
 300
 301       if(edepcell[j] <= cutoff)
 302         {
 303           fInfocl[0][id1][id2] = -1;
 304         }
 305     }
 306
 307   // ---------------------------------------------------------------
 308   // crude clustering begins. Start with cell having largest adc
 309   // count and loop over the cells in descending order of adc count
 310   // ---------------------------------------------------------------
 311
 312   icl       = -1;
 313   cellcount = -1;
 314
 315   for(icell = 0; icell <= nmx1; icell++)
 316     {
 317       j  = iord1[icell];
 318       id2 = j/kNDIMX;
 319       id1 = j-id2*kNDIMX;
 320
 321       if(fInfocl[0][id1][id2] == 0 )
 322         {
 323           icl++;
 324           numcell = 0;
 325           cellcount++;
 326           fInfocl[0][id1][id2] = 1;
 327           fInfocl[1][id1][id2] = icl;
 328           fInfcl[0][cellcount] = icl;
 329           fInfcl[1][cellcount] = id1;
 330           fInfcl[2][cellcount] = id2;
 331
 332           clust[0][numcell] = id1;
 333           clust[1][numcell] = id2;
 334
 335           for(i = 1; i < kndim; i++)
 336             {
 337               clust[0][i]=0;
 338             }
 339           // ---------------------------------------------------------------
 340           // check for adc count in neib. cells. If ne 0 put it in this clust
 341           // ---------------------------------------------------------------
 342           for(i = 0; i < 6; i++)
 343             {
 344               jd1 = id1 + neibx[i];
 345               jd2 = id2 + neiby[i];
 346               if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
 347                   fInfocl[0][jd1][jd2] == 0)
 348                 {
 349                   numcell++;
 350                   fInfocl[0][jd1][jd2] = 2;
 351                   fInfocl[1][jd1][jd2] = icl;
 352                   clust[0][numcell]    = jd1;
 353                   clust[1][numcell]    = jd2;
 354                   cellcount++;
 355                   fInfcl[0][cellcount] = icl;
 356                   fInfcl[1][cellcount] = jd1;
 357                   fInfcl[2][cellcount] = jd2;
 358                 }
 359             }
 360           // ---------------------------------------------------------------
 361           // check adc count for neighbour's neighbours recursively and
 362           // if nonzero, add these to the cluster.
 363           // ---------------------------------------------------------------
 364           for(i = 1; i < kndim;i++)
 365             {
 366               if(clust[0][i] != 0)
 367                 {
 368                   id1 = clust[0][i];
 369                   id2 = clust[1][i];
 370                   for(j = 0; j < 6 ; j++)
 371                     {
 372                       jd1 = id1 + neibx[j];
 373                       jd2 = id2 + neiby[j];
 374                       if( (jd1 >= 0 && jd1 < kNDIMX) &&
 375                           (jd2 >= 0 && jd2 < kNDIMY) &&
 376                           fInfocl[0][jd1][jd2] == 0 )
 377                         {
 378                           fInfocl[0][jd1][jd2] = 2;
 379                           fInfocl[1][jd1][jd2] = icl;
 380                           numcell++;
 381                           clust[0][numcell]    = jd1;
 382                           clust[1][numcell]    = jd2;
 383                           cellcount++;
 384                           fInfcl[0][cellcount] = icl;
 385                           fInfcl[1][cellcount] = jd1;
 386                           fInfcl[2][cellcount] = jd2;
 387                         }
 388                     }
 389                 }
 390             }
 391         }
 392     }
 393   return cellcount;
 394 }
 395 // ------------------------------------------------------------------------ //
 396 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
 397 {
 398   // Does the refining of clusters
 399   // Takes the big patch and does gaussian fitting and
 400   // finds out the more refined clusters
 401   //
 402
 403   AliPMDcludata *pmdcludata = 0;
 404
 405   const Int_t kNmaxCell   = 19;    // # of cells surrounding a cluster center
 406
 407   Int_t ndim = incr + 1;
 408
 409   Int_t    *ncl  = 0x0;
 410   Int_t    *clxy = 0x0;
 411   Int_t    i12, i22;
 412   Int_t    i, j, k, i1, i2, id, icl,  itest,ihld, ig, nsupcl,clno, t1, t2;
 413   Float_t  clusdata[6];
 414   Double_t x1, y1, z1, x2, y2, z2, rr;
 415
 416   ncl   = new Int_t [ndim];
 417   clxy  = new Int_t [kNmaxCell];
 418
 419   // Initialisation
 420   for(i = 0; i<ndim; i++)
 421     {
 422       ncl[i]  = -1;
 423       if (i < 6) clusdata[i] = 0.;
 424       if (i < kNmaxCell) clxy[i]    = 0;
 425     }
 426
 427   // clno counts the final clusters
 428   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
 429   // x, y and z store (x,y) coordinates of and energy deposited in a cell
 430   // xc, yc store (x,y) coordinates of the cluster center
 431   // zc stores the energy deposited in a cluster
 432   // rc is cluster radius
 433
 434   clno  = -1;
 435   nsupcl = -1;
 436
 437   for(i = 0; i <= incr; i++)
 438     {
 439       if(fInfcl[0][i] != nsupcl)
 440         {
 441           nsupcl++;
 442         }
 443       if (nsupcl > ndim)
 444         {
 445           AliWarning("RefClust: Too many superclusters!");
 446           nsupcl = ndim;
 447           break;
 448         }
 449       ncl[nsupcl]++;
 450     }
 451
 452   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
 453                   incr+1,nsupcl+1));
 454   id  = -1;
 455   icl = -1;
 456
 457   for(i = 0; i <= nsupcl; i++)
 458     {
 459       if(ncl[i] == 0)
 460         {
 461           id++;
 462           icl++;
 463           if (clno >= 4608)
 464             {
 465               AliWarning("RefClust: Too many clusters! more than 4608");
 466               return;
 467             }
 468           clno++;
 469           i1 = fInfcl[1][id];
 470           i2 = fInfcl[2][id];
 471
 472           i12 = i1 + i2*kNDIMX;
 473
 474           clusdata[0] = fCoord[0][i1][i2];
 475           clusdata[1] = fCoord[1][i1][i2];
 476           clusdata[2] = edepcell[i12];
 477           clusdata[3] = 1.;
 478           clusdata[4] = 0.5;
 479           clusdata[5] = 0.0;
 480
 481           clxy[0] = i1*10000 + i2;
 482
 483           for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
 484             {
 485               clxy[icltr] = -1;
 486             }
 487           pmdcludata  = new AliPMDcludata(clusdata,clxy);
 488           fPMDclucont->Add(pmdcludata);
 489         }
 490       else if(ncl[i] == 1)
 491         {
 492           id++;
 493           icl++;
 494           if (clno >= 4608)
 495             {
 496               AliWarning("RefClust: Too many clusters! more than 4608");
 497               return;
 498             }
 499           clno++;
 500           i1   = fInfcl[1][id];
 501           i2   = fInfcl[2][id];
 502           i12  = i1 + i2*kNDIMX;
 503
 504           x1   = fCoord[0][i1][i2];
 505           y1   = fCoord[1][i1][i2];
 506           z1   = edepcell[i12];
 507
 508           clxy[0] = i1*10000 + i2;
 509
 510           id++;
 511           i1   = fInfcl[1][id];
 512           i2   = fInfcl[2][id];
 513
 514           i22  = i1 + i2*kNDIMX;
 515           x2   = fCoord[0][i1][i2];
 516           y2   = fCoord[1][i1][i2];
 517           z2   = edepcell[i22];
 518
 519           clxy[1] = i1*10000 + i2;
 520
 521
 522           for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
 523             {
 524               clxy[icltr] = -1;
 525             }
 526
 527           clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
 528           clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
 529           clusdata[2] = z1+z2;
 530           clusdata[3] = 2.;
 531           clusdata[4] = 0.5;
 532           clusdata[5] = 0.0;
 533           pmdcludata  = new AliPMDcludata(clusdata,clxy);
 534           fPMDclucont->Add(pmdcludata);
 535         }
 536       else
 537         {
 538
 539           Int_t    *iord, *tc, *t;
 540           Double_t *x, *y, *z, *xc, *yc, *zc;
 541
 542           iord = new Int_t [ncl[i]+1];
 543           tc   = new Int_t [ncl[i]+1];
 544           t    = new Int_t [ncl[i]+1];
 545
 546           x    = new Double_t [ncl[i]+1];
 547           y    = new Double_t [ncl[i]+1];
 548           z    = new Double_t [ncl[i]+1];
 549           xc   = new Double_t [ncl[i]+1];
 550           yc   = new Double_t [ncl[i]+1];
 551           zc   = new Double_t [ncl[i]+1];
 552
 553           for( k = 0; k < ncl[i]+1; k++)
 554             {
 555               iord[k] = -1;
 556               t[k]    = -1;
 557               tc[k]   = -1;
 558               x[k]    = -1;
 559               y[k]    = -1;
 560               z[k]    = -1;
 561               xc[k]   = -1;
 562               yc[k]   = -1;
 563               zc[k]   = -1;
 564             }
 565           id++;
 566           // super-cluster of more than two cells - broken up into smaller
 567           // clusters gaussian centers computed. (peaks separated by > 1 cell)
 568           // Begin from cell having largest energy deposited This is first
 569           // cluster center
 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           t[0]    = i1*10000 + i2;
 578
 579
 580           iord[0] = 0;
 581           for(j = 1; j <= ncl[i]; j++)
 582             {
 583               id++;
 584               i1      = fInfcl[1][id];
 585               i2      = fInfcl[2][id];
 586               i12     = i1 + i2*kNDIMX;
 587
 588               iord[j] = j;
 589               x[j]    = fCoord[0][i1][i2];
 590               y[j]    = fCoord[1][i1][i2];
 591               z[j]    = edepcell[i12];
 592               t[j]    = i1*10000 + i2;
 593
 594             }
 595
 596           // arranging cells within supercluster in decreasing order
 597
 598           for(j = 1;j <= ncl[i]; j++)
 599             {
 600               itest = 0;
 601               ihld  = iord[j];
 602               for(i1 = 0; i1 < j; i1++)
 603                 {
 604                   if(itest == 0 && z[iord[i1]] < z[ihld])
 605                     {
 606                       itest = 1;
 607                       for(i2 = j-1; i2 >= i1; i2--)
 608                         {
 609                           iord[i2+1] = iord[i2];
 610                         }
 611                       iord[i1] = ihld;
 612                     }
 613                 }
 614             }
 615           /* MODIFICATION PART STARTS (Tapan July 2008)
 616              iord[0] is the cell with highest ADC in the crude-cluster
 617              ig is the number of local maxima in the crude-cluster
 618              For the higest peak we make ig=0 which means first local maximum.
 619              Next we go down in terms of the ADC sequence and find out if any
 620              more of the cells form local maxima. The definition of local
 621              maxima is that all its neighbours are of less ADC compared to it.
 622           */
 623           ig = 0;
 624           xc[ig] = x[iord[0]];
 625           yc[ig] = y[iord[0]];
 626           zc[ig] = z[iord[0]];
 627           tc[ig] = t[iord[0]];
 628           Int_t ivalid = 0,  icount = 0;
 629
 630           for(j=1;j<=ncl[i];j++)
 631             {
 632               x1  = x[iord[j]];
 633               y1  = y[iord[j]];
 634               z1  = z[iord[j]];
 635               t1  = t[iord[j]];
 636               rr=Distance(x1,y1,xc[ig],yc[ig]);
 637
 638               // Check the cells which are outside the neighbours (rr>1.2)
 639               if(rr>1.2 )
 640                 {
 641                   ivalid=0;
 642                   icount=0;
 643                   for(Int_t j1=1;j1<j;j1++)
 644                     {
 645                       icount++;
 646                       Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
 647                       if(rr1>1.2) ivalid++;
 648                     }
 649                   if(ivalid == icount && z1>0.5*zc[ig])
 650                     {
 651                       ig++;
 652                       xc[ig]=x1;
 653                       yc[ig]=y1;
 654                       zc[ig]=z1;
 655                       tc[ig]=t1;
 656                     }
 657                 }
 658             }
 659
 660           icl=icl+ig+1;
 661
 662           //  We use simple Gaussian weighting. (Tapan Jan 2005)
 663           // compute the number of cells belonging to each cluster.
 664           // cell can be shared between several clusters
 665           // in the ratio of cluster energy deposition
 666           // To calculate:
 667           // (1) number of cells belonging to a cluster (ig) and
 668           // (2) total ADC of the cluster (ig)
 669           // (3) x and y positions of the cluster
 670
 671
 672           Int_t *cellCount;
 673           Int_t **cellXY;
 674
 675           Int_t    *status;
 676           Double_t *totaladc, *totaladc2, *ncell,*weight;
 677           Double_t *xclust, *yclust, *sigxclust, *sigyclust;
 678           Double_t *ax, *ay, *ax2, *ay2;
 679
 680
 681           status    = new Int_t [ncl[i]+1];
 682           cellXY    = new Int_t *[ncl[i]+1];
 683
 684           cellCount = new Int_t [ig+1];
 685           totaladc  = new Double_t [ig+1];
 686           totaladc2 = new Double_t [ig+1];
 687           ncell     = new Double_t [ig+1];
 688           weight    = new Double_t [ig+1];
 689           xclust    = new Double_t [ig+1];
 690           yclust    = new Double_t [ig+1];
 691           sigxclust = new Double_t [ig+1];
 692           sigyclust = new Double_t [ig+1];
 693           ax        = new Double_t [ig+1];
 694           ay        = new Double_t [ig+1];
 695           ax2       = new Double_t [ig+1];
 696           ay2       = new Double_t [ig+1];
 697
 698           for(j = 0; j < ncl[i]+1; j++)
 699             {
 700               status[j]     = 0;
 701               cellXY[j] = new Int_t[ig+1];
 702             }
 703           //initialization
 704           for(Int_t kcl = 0; kcl < ig+1; kcl++)
 705             {
 706               cellCount[kcl] = 0;
 707               totaladc[kcl]  = 0.;
 708               totaladc2[kcl] = 0.;
 709               ncell[kcl]     = 0.;
 710               weight[kcl]    = 0.;
 711               xclust[kcl]    = 0.;
 712               yclust[kcl]    = 0.;
 713               sigxclust[kcl] = 0.;
 714               sigyclust[kcl] = 0.;
 715               ax[kcl]        = 0.;
 716               ay[kcl]        = 0.;
 717               ax2[kcl]       = 0.;
 718               ay2[kcl]       = 0.;
 719               for(j = 0; j < ncl[i]+1; j++)
 720                 {
 721                   cellXY[j][kcl] = 0;
 722                 }
 723             }
 724           Double_t sumweight, gweight;
 725
 726           for(j = 0;j <= ncl[i]; j++)
 727             {
 728               x1 = x[iord[j]];
 729               y1 = y[iord[j]];
 730               z1 = z[iord[j]];
 731               t1 = t[iord[j]];
 732
 733               for(Int_t kcl=0; kcl<=ig; kcl++)
 734                 {
 735                   x2 = xc[kcl];
 736                   y2 = yc[kcl];
 737                   rr = Distance(x1,y1,x2,y2);
 738                   t2 = tc[kcl];
 739
 740                   if(rr==0)
 741                     {
 742                       ncell[kcl] = 1.;
 743                       totaladc[kcl]  = z1;
 744                       totaladc2[kcl]  = pow(z1,2);
 745                       ax[kcl] =  x1 * z1;
 746                       ay[kcl] =  y1 * z1;
 747                       ax2[kcl]=  0.;
 748                       ay2[kcl]=  0.;
 749                       status[j] = 1;
 750                     }
 751                 }
 752             }
 753
 754           for(j = 0; j <= ncl[i]; j++)
 755             {
 756               Int_t   maxweight = 0;
 757               Double_t     max  = 0.;
 758
 759               if(status[j] == 0)
 760                 {
 761                   x1 = x[iord[j]];
 762                   y1 = y[iord[j]];
 763                   z1 = z[iord[j]];
 764                   t1 = t[iord[j]];
 765                   sumweight = 0.;
 766
 767                   for(Int_t kcl = 0; kcl <= ig; kcl++)
 768                     {
 769                       x2 = xc[kcl];
 770                       y2 = yc[kcl];
 771                       rr = Distance(x1,y1,x2,y2);
 772                       gweight     = exp(-(rr*rr)/(2*(1.2*1.2)));
 773                       weight[kcl] = zc[kcl] * gweight;
 774                       sumweight   = sumweight + weight[kcl];
 775
 776                       if(weight[kcl] > max)
 777                         {
 778                           max       =  weight[kcl];
 779                           maxweight =  kcl;
 780                         }
 781                     }
 782
 783                   cellXY[cellCount[maxweight]][maxweight] = iord[j];
 784
 785                   cellCount[maxweight]++;
 786
 787                   for(Int_t kcl = 0; kcl <= ig; kcl++)
 788                     {
 789                       x2 = xc[kcl];
 790                       y2 = yc[kcl];
 791                       rr = Distance(x1,y1,x2,y2);
 792                       t2 = tc[kcl];
 793                       // For calculating weighted mean:
 794                       // Weighted_mean = (\sum w_i x_i) / (\sum w_i)
 795
 796                       if(sumweight>0.)
 797                         {
 798                           totaladc[kcl] = totaladc[kcl] + z1*(weight[kcl]/sumweight);
 799                           ncell[kcl]    = ncell[kcl]    + (weight[kcl]/sumweight);
 800                           ax[kcl]       = ax[kcl]       + (x1 * z1 *weight[kcl]/sumweight);
 801                           ay[kcl]       = ay[kcl]       + (y1 * z1 *weight[kcl]/sumweight);
 802
 803                           // For calculating weighted sigma:
 804                           // Wieghted_sigma= ((\sum w_i)/((\sum w_i)^2 - \sum (w_i^2))) \sum w_i (x_i - \Hat\mu)^2
 805                           // I assume here x1,y1 are cluster centers, otherwise life becomes too difficult
 806                           totaladc2[kcl] = totaladc2[kcl] +  pow((z1*(weight[kcl]/sumweight)),2);
 807                           ax2[kcl] = ax2[kcl] + (z1 *weight[kcl]/sumweight) * pow((x1-x2),2);
 808                           ay2[kcl] = ay2[kcl] + (z1 *weight[kcl]/sumweight) * pow((y1-y2),2);
 809                         }
 810                     }
 811                 }
 812             }
 813
 814           for(Int_t kcl = 0; kcl <= ig; kcl++)
 815             {
 816
 817               if(totaladc[kcl]>0){
 818                 if(totaladc[kcl]>0.)xclust[kcl] = (ax[kcl])/ totaladc[kcl];
 819                 if(totaladc[kcl]>0.)yclust[kcl] = (ay[kcl])/ totaladc[kcl];
 820
 821                 if(totaladc[kcl]>0.)sigxclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ax2[kcl];
 822                 if(totaladc[kcl]>0.)sigyclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ay2[kcl];
 823               }
 824
 825               for(j = 0; j < cellCount[kcl]; j++) clno++;
 826
 827               if (clno >= 4608)
 828                 {
 829                   AliWarning("RefClust: Too many clusters! more than 4608");
 830                   return;
 831                 }
 832               clusdata[0] = xclust[kcl];
 833               clusdata[1] = yclust[kcl];
 834               clusdata[2] = totaladc[kcl];
 835               clusdata[3] = ncell[kcl];
 836               if(sigxclust[kcl] > sigyclust[kcl])
 837                 {
 838                   clusdata[4] = pow(sigxclust[kcl],0.5);
 839                   clusdata[5] = pow(sigyclust[kcl],0.5);
 840                 }
 841               else
 842                 {
 843                   clusdata[4] = pow(sigyclust[kcl],0.5);
 844                   clusdata[5] = pow(sigxclust[kcl],0.5);
 845                 }
 846
 847               clxy[0] = tc[kcl];
 848
 849               Int_t Ncell=1;
 850               for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
 851                 {
 852                   Float_t xx = x[cellXY[ii][kcl]];
 853                   Float_t yy = y[cellXY[ii][kcl]];
 854
 855                   rr=Distance(xclust[kcl],yclust[kcl],xx,yy);
 856                   // Natasha
 857                   // We store only the nearest and nest-nearest neighbours
 858                   if(rr<2.2)
 859                     {
 860                       clxy[Ncell] = t[cellXY[ii][kcl]];
 861                       Ncell++;
 862                     }
 863                 }
 864
 865               pmdcludata = new AliPMDcludata(clusdata,clxy);
 866               fPMDclucont->Add(pmdcludata);
 867             }
 868
 869           delete [] iord;
 870           delete [] tc;
 871           delete [] t;
 872           delete [] x;
 873           delete [] y;
 874           delete [] z;
 875           delete [] xc;
 876           delete [] yc;
 877           delete [] zc;
 878
 879           delete [] cellCount;
 880           for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
 881
 882           delete [] status;
 883           delete [] totaladc;
 884           delete [] totaladc2;
 885           delete [] ncell;
 886           delete [] xclust;
 887           delete [] yclust;
 888           delete [] sigxclust;
 889           delete [] sigyclust;
 890           delete [] ax;
 891           delete [] ay;
 892           delete [] ax2;
 893           delete [] ay2;
 894           delete [] weight;
 895         }
 896     }
 897   delete [] ncl;
 898   delete [] clxy;
 899 }
 900 // ------------------------------------------------------------------------ //
 901 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
 902                                       Double_t x2, Double_t y2)
 903 {
 904   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
 905 }
 906 // ------------------------------------------------------------------------ //
 907 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
 908 {
 909   fCutoff = decut;
 910 }
 911 // ------------------------------------------------------------------------ //