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