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