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   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   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 (ismn < 12)
 274             {
 275               if ((irow >= 0 && irow < 96) && (icol >= 0 && icol < 48))
 276                 {
 277                   celldataTr[ihit]  = celltrack[icol][irow];
 278                   celldataPid[ihit] = cellpid[icol][irow];
 279                   celldataAdc[ihit] = (Float_t) celladc[icol][irow];
 280                 }
 281               else
 282                 {
 283                   celldataTr[ihit]  = -1;
 284                   celldataPid[ihit] = -1;
 285                   celldataAdc[ihit] = -1;
 286                 }
 287             }
 288           else if (ismn >= 12 && ismn < 24)
 289             {
 290               if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
 291                 {
 292                   celldataTr[ihit]  = celltrack[irow][icol];
 293                   celldataPid[ihit] = cellpid[irow][icol];
 294                   celldataAdc[ihit] = (Float_t) celladc[irow][icol];
 295
 296                 }
 297               else
 298                 {
 299                   celldataTr[ihit]  = -1;
 300                   celldataPid[ihit] = -1;
 301                   celldataAdc[ihit] = -1;
 302                 }
 303             }
 304
 305         }
 306
 307       pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
 308                                 celldataTr, celldataPid, celldataAdc);
 309       pmdcont->Add(pmdcl);
 310     }
 311
 312   fPMDclucont->Delete();
 313 }
 314 // ------------------------------------------------------------------------ //
 315 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
 316                                   Int_t iord1[], Double_t edepcell[])
 317 {
 318   // Does crude clustering
 319   // Finds out only the big patch by just searching the
 320   // connected cells
 321   //
 322   const Int_t kndim = 4609;
 323   Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
 324   Int_t jd1,jd2, icell, cellcount;
 325   static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
 326
 327   AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
 328
 329   for (j = 0; j < kNDIMX; j++)
 330     {
 331       for(k = 0; k < kNDIMY; k++)
 332         {
 333           fInfocl[0][j][k] = 0;
 334           fInfocl[1][j][k] = 0;
 335         }
 336     }
 337   for(i=0; i < kNMX; i++)
 338     {
 339       fInfcl[0][i] = -1;
 340
 341       j  = iord1[i];
 342       id2 = j/kNDIMX;
 343       id1 = j-id2*kNDIMX;
 344
 345       if(edepcell[j] <= cutoff)
 346         {
 347           fInfocl[0][id1][id2] = -1;
 348         }
 349     }
 350
 351   // ---------------------------------------------------------------
 352   // crude clustering begins. Start with cell having largest adc
 353   // count and loop over the cells in descending order of adc count
 354   // ---------------------------------------------------------------
 355
 356   icl       = -1;
 357   cellcount = -1;
 358
 359   for(icell = 0; icell <= nmx1; icell++)
 360     {
 361       j  = iord1[icell];
 362       id2 = j/kNDIMX;
 363       id1 = j-id2*kNDIMX;
 364
 365       if(fInfocl[0][id1][id2] == 0 )
 366         {
 367           icl++;
 368           numcell = 0;
 369           cellcount++;
 370           fInfocl[0][id1][id2] = 1;
 371           fInfocl[1][id1][id2] = icl;
 372           fInfcl[0][cellcount] = icl;
 373           fInfcl[1][cellcount] = id1;
 374           fInfcl[2][cellcount] = id2;
 375
 376           clust[0][numcell] = id1;
 377           clust[1][numcell] = id2;
 378
 379           for(i = 1; i < kndim; i++)
 380             {
 381               clust[0][i]=0;
 382             }
 383           // ---------------------------------------------------------------
 384           // check for adc count in neib. cells. If ne 0 put it in this clust
 385           // ---------------------------------------------------------------
 386           for(i = 0; i < 6; i++)
 387             {
 388               jd1 = id1 + neibx[i];
 389               jd2 = id2 + neiby[i];
 390               if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
 391                   fInfocl[0][jd1][jd2] == 0)
 392                 {
 393                   numcell++;
 394                   fInfocl[0][jd1][jd2] = 2;
 395                   fInfocl[1][jd1][jd2] = icl;
 396                   clust[0][numcell]    = jd1;
 397                   clust[1][numcell]    = jd2;
 398                   cellcount++;
 399                   fInfcl[0][cellcount] = icl;
 400                   fInfcl[1][cellcount] = jd1;
 401                   fInfcl[2][cellcount] = jd2;
 402                 }
 403             }
 404           // ---------------------------------------------------------------
 405           // check adc count for neighbour's neighbours recursively and
 406           // if nonzero, add these to the cluster.
 407           // ---------------------------------------------------------------
 408           for(i = 1; i < kndim;i++)
 409             {
 410               if(clust[0][i] != 0)
 411                 {
 412                   id1 = clust[0][i];
 413                   id2 = clust[1][i];
 414                   for(j = 0; j < 6 ; j++)
 415                     {
 416                       jd1 = id1 + neibx[j];
 417                       jd2 = id2 + neiby[j];
 418                       if( (jd1 >= 0 && jd1 < kNDIMX) &&
 419                           (jd2 >= 0 && jd2 < kNDIMY) &&
 420                           fInfocl[0][jd1][jd2] == 0 )
 421                         {
 422                           fInfocl[0][jd1][jd2] = 2;
 423                           fInfocl[1][jd1][jd2] = icl;
 424                           numcell++;
 425                           clust[0][numcell]    = jd1;
 426                           clust[1][numcell]    = jd2;
 427                           cellcount++;
 428                           fInfcl[0][cellcount] = icl;
 429                           fInfcl[1][cellcount] = jd1;
 430                           fInfcl[2][cellcount] = jd2;
 431                         }
 432                     }
 433                 }
 434             }
 435         }
 436     }
 437   return cellcount;
 438 }
 439 // ------------------------------------------------------------------------ //
 440 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
 441 {
 442   // Does the refining of clusters
 443   // Takes the big patch and does gaussian fitting and
 444   // finds out the more refined clusters
 445   //
 446
 447   AliPMDcludata *pmdcludata = 0;
 448
 449   const Int_t kNmaxCell   = 19;    // # of cells surrounding a cluster center
 450
 451   Int_t ndim = incr + 1;
 452
 453   Int_t    *ncl  = 0x0;
 454   Int_t    *clxy = 0x0;
 455   Int_t    i12, i22;
 456   Int_t    i, j, k, i1, i2, id, icl,  itest,ihld, ig, nsupcl,clno, t1, t2;
 457   Float_t  clusdata[6];
 458   Double_t x1, y1, z1, x2, y2, z2, rr;
 459
 460   ncl   = new Int_t [ndim];
 461   clxy  = new Int_t [kNmaxCell];
 462
 463   // Initialisation
 464   for(i = 0; i<ndim; i++)
 465     {
 466       ncl[i]  = -1;
 467       if (i < 6) clusdata[i] = 0.;
 468       if (i < kNmaxCell) clxy[i]    = 0;
 469     }
 470
 471   // clno counts the final clusters
 472   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
 473   // x, y and z store (x,y) coordinates of and energy deposited in a cell
 474   // xc, yc store (x,y) coordinates of the cluster center
 475   // zc stores the energy deposited in a cluster
 476   // rc is cluster radius
 477
 478   clno  = -1;
 479   nsupcl = -1;
 480
 481   for(i = 0; i <= incr; i++)
 482     {
 483       if(fInfcl[0][i] != nsupcl)
 484         {
 485           nsupcl++;
 486         }
 487       if (nsupcl > ndim)
 488         {
 489           AliWarning("RefClust: Too many superclusters!");
 490           nsupcl = ndim;
 491           break;
 492         }
 493       ncl[nsupcl]++;
 494     }
 495
 496   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
 497                   incr+1,nsupcl+1));
 498   id  = -1;
 499   icl = -1;
 500
 501   for(i = 0; i <= nsupcl; i++)
 502     {
 503       if(ncl[i] == 0)
 504         {
 505           id++;
 506           icl++;
 507           if (clno >= 4608)
 508             {
 509               AliWarning("RefClust: Too many clusters! more than 4608");
 510               return;
 511             }
 512           clno++;
 513           i1 = fInfcl[1][id];
 514           i2 = fInfcl[2][id];
 515
 516           i12 = i1 + i2*kNDIMX;
 517
 518           clusdata[0] = fCoord[0][i1][i2];
 519           clusdata[1] = fCoord[1][i1][i2];
 520           clusdata[2] = edepcell[i12];
 521           clusdata[3] = 1.;
 522           clusdata[4] = 0.5;
 523           clusdata[5] = 0.0;
 524
 525           clxy[0] = i1*10000 + i2;
 526
 527           for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
 528             {
 529               clxy[icltr] = -1;
 530             }
 531
 532           pmdcludata  = new AliPMDcludata(clusdata,clxy);
 533           fPMDclucont->Add(pmdcludata);
 534         }
 535       else if(ncl[i] == 1)
 536         {
 537           id++;
 538           icl++;
 539           if (clno >= 4608)
 540             {
 541               AliWarning("RefClust: Too many clusters! more than 4608");
 542               return;
 543             }
 544           clno++;
 545           i1   = fInfcl[1][id];
 546           i2   = fInfcl[2][id];
 547           i12  = i1 + i2*kNDIMX;
 548
 549           x1   = fCoord[0][i1][i2];
 550           y1   = fCoord[1][i1][i2];
 551           z1   = edepcell[i12];
 552
 553           clxy[0] = i1*10000 + i2;
 554
 555           id++;
 556           i1   = fInfcl[1][id];
 557           i2   = fInfcl[2][id];
 558
 559           i22  = i1 + i2*kNDIMX;
 560           x2   = fCoord[0][i1][i2];
 561           y2   = fCoord[1][i1][i2];
 562           z2   = edepcell[i22];
 563
 564           clxy[1] = i1*10000 + i2;
 565
 566
 567           for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
 568             {
 569               clxy[icltr] = -1;
 570             }
 571
 572           clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
 573           clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
 574           clusdata[2] = z1+z2;
 575           clusdata[3] = 2.;
 576           clusdata[4] = 0.5;
 577           clusdata[5] = 0.0;
 578           pmdcludata  = new AliPMDcludata(clusdata,clxy);
 579           fPMDclucont->Add(pmdcludata);
 580         }
 581       else
 582         {
 583
 584           Int_t    *iord, *tc, *t;
 585           Double_t *x, *y, *z, *xc, *yc, *zc;
 586
 587           iord = new Int_t [ncl[i]+1];
 588           tc   = new Int_t [ncl[i]+1];
 589           t    = new Int_t [ncl[i]+1];
 590
 591           x    = new Double_t [ncl[i]+1];
 592           y    = new Double_t [ncl[i]+1];
 593           z    = new Double_t [ncl[i]+1];
 594           xc   = new Double_t [ncl[i]+1];
 595           yc   = new Double_t [ncl[i]+1];
 596           zc   = new Double_t [ncl[i]+1];
 597
 598           for( k = 0; k < ncl[i]+1; k++)
 599             {
 600               iord[k] = -1;
 601               t[k]    = -1;
 602               tc[k]   = -1;
 603               x[k]    = -1;
 604               y[k]    = -1;
 605               z[k]    = -1;
 606               xc[k]   = -1;
 607               yc[k]   = -1;
 608               zc[k]   = -1;
 609             }
 610           id++;
 611           // super-cluster of more than two cells - broken up into smaller
 612           // clusters gaussian centers computed. (peaks separated by > 1 cell)
 613           // Begin from cell having largest energy deposited This is first
 614           // cluster center
 615           i1      = fInfcl[1][id];
 616           i2      = fInfcl[2][id];
 617           i12     = i1 + i2*kNDIMX;
 618
 619           x[0]    = fCoord[0][i1][i2];
 620           y[0]    = fCoord[1][i1][i2];
 621           z[0]    = edepcell[i12];
 622           t[0]    = i1*10000 + i2;
 623
 624
 625           iord[0] = 0;
 626           for(j = 1; j <= ncl[i]; j++)
 627             {
 628               id++;
 629               i1      = fInfcl[1][id];
 630               i2      = fInfcl[2][id];
 631               i12     = i1 + i2*kNDIMX;
 632
 633               iord[j] = j;
 634               x[j]    = fCoord[0][i1][i2];
 635               y[j]    = fCoord[1][i1][i2];
 636               z[j]    = edepcell[i12];
 637               t[j]    = i1*10000 + i2;
 638
 639             }
 640
 641           // arranging cells within supercluster in decreasing order
 642
 643           for(j = 1;j <= ncl[i]; j++)
 644             {
 645               itest = 0;
 646               ihld  = iord[j];
 647               for(i1 = 0; i1 < j; i1++)
 648                 {
 649                   if(itest == 0 && z[iord[i1]] < z[ihld])
 650                     {
 651                       itest = 1;
 652                       for(i2 = j-1; i2 >= i1; i2--)
 653                         {
 654                           iord[i2+1] = iord[i2];
 655                         }
 656                       iord[i1] = ihld;
 657                     }
 658                 }
 659             }
 660           /* MODIFICATION PART STARTS (Tapan July 2008)
 661              iord[0] is the cell with highest ADC in the crude-cluster
 662              ig is the number of local maxima in the crude-cluster
 663              For the higest peak we make ig=0 which means first local maximum.
 664              Next we go down in terms of the ADC sequence and find out if any
 665              more of the cells form local maxima. The definition of local
 666              maxima is that all its neighbours are of less ADC compared to it.
 667           */
 668           ig = 0;
 669           xc[ig] = x[iord[0]];
 670           yc[ig] = y[iord[0]];
 671           zc[ig] = z[iord[0]];
 672           tc[ig] = t[iord[0]];
 673           Int_t ivalid = 0,  icount = 0;
 674
 675           for(j=1;j<=ncl[i];j++)
 676             {
 677               x1  = x[iord[j]];
 678               y1  = y[iord[j]];
 679               z1  = z[iord[j]];
 680               t1  = t[iord[j]];
 681               rr=Distance(x1,y1,xc[ig],yc[ig]);
 682
 683               // Check the cells which are outside the neighbours (rr>1.2)
 684               if(rr>1.2 )
 685                 {
 686                   ivalid=0;
 687                   icount=0;
 688                   for(Int_t j1=1;j1<j;j1++)
 689                     {
 690                       icount++;
 691                       Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
 692                       if(rr1>1.2) ivalid++;
 693                     }
 694                   if(ivalid == icount && z1>0.5*zc[ig])
 695                     {
 696                       ig++;
 697                       xc[ig]=x1;
 698                       yc[ig]=y1;
 699                       zc[ig]=z1;
 700                       tc[ig]=t1;
 701                     }
 702                 }
 703             }
 704
 705           icl=icl+ig+1;
 706
 707           //  We use simple Gaussian weighting. (Tapan Jan 2005)
 708           // compute the number of cells belonging to each cluster.
 709           // cell can be shared between several clusters
 710           // in the ratio of cluster energy deposition
 711           // To calculate:
 712           // (1) number of cells belonging to a cluster (ig) and
 713           // (2) total ADC of the cluster (ig)
 714           // (3) x and y positions of the cluster
 715
 716
 717           Int_t *cellCount;
 718           Int_t **cellXY;
 719
 720           Int_t    *status;
 721           Double_t *totaladc, *totaladc2, *ncell,*weight;
 722           Double_t *xclust, *yclust, *sigxclust, *sigyclust;
 723           Double_t *ax, *ay, *ax2, *ay2;
 724
 725
 726           status    = new Int_t [ncl[i]+1];
 727           cellXY    = new Int_t *[ncl[i]+1];
 728
 729           cellCount = new Int_t [ig+1];
 730           totaladc  = new Double_t [ig+1];
 731           totaladc2 = new Double_t [ig+1];
 732           ncell     = new Double_t [ig+1];
 733           weight    = new Double_t [ig+1];
 734           xclust    = new Double_t [ig+1];
 735           yclust    = new Double_t [ig+1];
 736           sigxclust = new Double_t [ig+1];
 737           sigyclust = new Double_t [ig+1];
 738           ax        = new Double_t [ig+1];
 739           ay        = new Double_t [ig+1];
 740           ax2       = new Double_t [ig+1];
 741           ay2       = new Double_t [ig+1];
 742
 743           for(j = 0; j < ncl[i]+1; j++)
 744             {
 745               status[j]     = 0;
 746               cellXY[j] = new Int_t[ig+1];
 747             }
 748           //initialization
 749           for(Int_t kcl = 0; kcl < ig+1; kcl++)
 750             {
 751               cellCount[kcl] = 0;
 752               totaladc[kcl]  = 0.;
 753               totaladc2[kcl] = 0.;
 754               ncell[kcl]     = 0.;
 755               weight[kcl]    = 0.;
 756               xclust[kcl]    = 0.;
 757               yclust[kcl]    = 0.;
 758               sigxclust[kcl] = 0.;
 759               sigyclust[kcl] = 0.;
 760               ax[kcl]        = 0.;
 761               ay[kcl]        = 0.;
 762               ax2[kcl]       = 0.;
 763               ay2[kcl]       = 0.;
 764               for(j = 0; j < ncl[i]+1; j++)
 765                 {
 766                   cellXY[j][kcl] = 0;
 767                 }
 768             }
 769           Double_t sumweight, gweight;
 770
 771           for(j = 0;j <= ncl[i]; j++)
 772             {
 773               x1 = x[iord[j]];
 774               y1 = y[iord[j]];
 775               z1 = z[iord[j]];
 776               t1 = t[iord[j]];
 777
 778               for(Int_t kcl=0; kcl<=ig; kcl++)
 779                 {
 780                   x2 = xc[kcl];
 781                   y2 = yc[kcl];
 782                   rr = Distance(x1,y1,x2,y2);
 783                   t2 = tc[kcl];
 784
 785                   if(rr==0)
 786                     {
 787                       ncell[kcl] = 1.;
 788                       totaladc[kcl]  = z1;
 789                       totaladc2[kcl]  = pow(z1,2);
 790                       ax[kcl] =  x1 * z1;
 791                       ay[kcl] =  y1 * z1;
 792                       ax2[kcl]=  0.;
 793                       ay2[kcl]=  0.;
 794                       status[j] = 1;
 795                     }
 796                 }
 797             }
 798
 799           for(j = 0; j <= ncl[i]; j++)
 800             {
 801               Int_t   maxweight = 0;
 802               Double_t     max  = 0.;
 803
 804               if(status[j] == 0)
 805                 {
 806                   x1 = x[iord[j]];
 807                   y1 = y[iord[j]];
 808                   z1 = z[iord[j]];
 809                   t1 = t[iord[j]];
 810                   sumweight = 0.;
 811
 812                   for(Int_t kcl = 0; kcl <= ig; kcl++)
 813                     {
 814                       x2 = xc[kcl];
 815                       y2 = yc[kcl];
 816                       rr = Distance(x1,y1,x2,y2);
 817                       gweight     = exp(-(rr*rr)/(2*(1.2*1.2)));
 818                       weight[kcl] = zc[kcl] * gweight;
 819                       sumweight   = sumweight + weight[kcl];
 820
 821                       if(weight[kcl] > max)
 822                         {
 823                           max       =  weight[kcl];
 824                           maxweight =  kcl;
 825                         }
 826                     }
 827
 828                   cellXY[cellCount[maxweight]][maxweight] = iord[j];
 829
 830                   cellCount[maxweight]++;
 831
 832                   x2 = xc[maxweight];
 833                   y2 = yc[maxweight];
 834                   totaladc[maxweight]  +=  z1;
 835                   ax[maxweight]        +=  x1 * z1;
 836                   ay[maxweight]        +=  y1 * z1;
 837                   totaladc2[maxweight] +=  pow(z1,2);
 838                   ax2[maxweight]       +=  z1 * pow((x1-x2),2);
 839                   ay2[maxweight]       +=  z1 * pow((y1-y2),2);
 840                   ncell[maxweight]++;
 841
 842                 }
 843             }
 844
 845           for(Int_t kcl = 0; kcl <= ig; kcl++)
 846             {
 847
 848               if(totaladc[kcl]>0){
 849                 if(totaladc[kcl]>0.)xclust[kcl] = (ax[kcl])/ totaladc[kcl];
 850                 if(totaladc[kcl]>0.)yclust[kcl] = (ay[kcl])/ totaladc[kcl];
 851
 852                 //natasha
 853                 if(totaladc2[kcl] >= pow(totaladc[kcl],2))
 854                   {
 855                     sigxclust[kcl] = 0.25;
 856                     sigyclust[kcl] = 0.25;
 857                   }
 858                 else
 859                   {
 860                     sigxclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ax2[kcl];
 861                     sigyclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-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] = pow(sigxclust[kcl],0.5);
 881                   clusdata[5] = pow(sigyclust[kcl],0.5);
 882                 }
 883               else
 884                 {
 885                   clusdata[4] = pow(sigyclust[kcl],0.5);
 886                   clusdata[5] = pow(sigxclust[kcl],0.5);
 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                   Float_t adc = (Float_t) celladc[id1][jd1];
 971                   if(adc == 0.)
 972                     {
 973                       isocount++;
 974                       if(isocount == kCellNeighbour)
 975                         {
 976                           Float_t cadc = (Float_t) celladc[irow][icol];
 977
 978                           isocell = new AliPMDisocell(idet,ismn,irow,icol,cadc);
 979                           pmdisocell->Add(isocell);
 980
 981                         }
 982                     }
 983                 }  // neigh cell cond.
 984             }
 985         }
 986     }
 987
 988
 989 }
 990 // ------------------------------------------------------------------------ //
 991 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
 992 {
 993   fCutoff = decut;
 994 }
 995 // ------------------------------------------------------------------------ //