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               else
 297                 {
 298                   celldataTr[ihit]  = -1;
 299                   celldataPid[ihit] = -1;
 300                   celldataAdc[ihit] = -1;
 301                 }
 302             }
 303
 304         }
 305
 306       pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
 307                                 celldataTr, celldataPid, celldataAdc);
 308       pmdcont->Add(pmdcl);
 309     }
 310
 311   fPMDclucont->Delete();
 312 }
 313 // ------------------------------------------------------------------------ //
 314 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
 315                                   Int_t iord1[], Double_t edepcell[])
 316 {
 317   // Does crude clustering
 318   // Finds out only the big patch by just searching the
 319   // connected cells
 320   //
 321   const Int_t kndim = 4609;
 322   Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
 323   Int_t jd1,jd2, icell, cellcount;
 324   static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
 325
 326   AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
 327
 328   for (j = 0; j < kNDIMX; j++)
 329     {
 330       for(k = 0; k < kNDIMY; k++)
 331         {
 332           fInfocl[0][j][k] = 0;
 333           fInfocl[1][j][k] = 0;
 334         }
 335     }
 336   for(i=0; i < kNMX; i++)
 337     {
 338       fInfcl[0][i] = -1;
 339
 340       j  = iord1[i];
 341       id2 = j/kNDIMX;
 342       id1 = j-id2*kNDIMX;
 343
 344       if(edepcell[j] <= cutoff)
 345         {
 346           fInfocl[0][id1][id2] = -1;
 347         }
 348     }
 349
 350   // ---------------------------------------------------------------
 351   // crude clustering begins. Start with cell having largest adc
 352   // count and loop over the cells in descending order of adc count
 353   // ---------------------------------------------------------------
 354
 355   icl       = -1;
 356   cellcount = -1;
 357
 358   for(icell = 0; icell <= nmx1; icell++)
 359     {
 360       j  = iord1[icell];
 361       id2 = j/kNDIMX;
 362       id1 = j-id2*kNDIMX;
 363
 364       if(fInfocl[0][id1][id2] == 0 )
 365         {
 366           icl++;
 367           numcell = 0;
 368           cellcount++;
 369           fInfocl[0][id1][id2] = 1;
 370           fInfocl[1][id1][id2] = icl;
 371           fInfcl[0][cellcount] = icl;
 372           fInfcl[1][cellcount] = id1;
 373           fInfcl[2][cellcount] = id2;
 374
 375           clust[0][numcell] = id1;
 376           clust[1][numcell] = id2;
 377
 378           for(i = 1; i < kndim; i++)
 379             {
 380               clust[0][i]=0;
 381             }
 382           // ---------------------------------------------------------------
 383           // check for adc count in neib. cells. If ne 0 put it in this clust
 384           // ---------------------------------------------------------------
 385           for(i = 0; i < 6; i++)
 386             {
 387               jd1 = id1 + neibx[i];
 388               jd2 = id2 + neiby[i];
 389               if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
 390                   fInfocl[0][jd1][jd2] == 0)
 391                 {
 392                   numcell++;
 393                   fInfocl[0][jd1][jd2] = 2;
 394                   fInfocl[1][jd1][jd2] = icl;
 395                   clust[0][numcell]    = jd1;
 396                   clust[1][numcell]    = jd2;
 397                   cellcount++;
 398                   fInfcl[0][cellcount] = icl;
 399                   fInfcl[1][cellcount] = jd1;
 400                   fInfcl[2][cellcount] = jd2;
 401                 }
 402             }
 403           // ---------------------------------------------------------------
 404           // check adc count for neighbour's neighbours recursively and
 405           // if nonzero, add these to the cluster.
 406           // ---------------------------------------------------------------
 407           for(i = 1; i < kndim;i++)
 408             {
 409               if(clust[0][i] != 0)
 410                 {
 411                   id1 = clust[0][i];
 412                   id2 = clust[1][i];
 413                   for(j = 0; j < 6 ; j++)
 414                     {
 415                       jd1 = id1 + neibx[j];
 416                       jd2 = id2 + neiby[j];
 417                       if( (jd1 >= 0 && jd1 < kNDIMX) &&
 418                           (jd2 >= 0 && jd2 < kNDIMY) &&
 419                           fInfocl[0][jd1][jd2] == 0 )
 420                         {
 421                           fInfocl[0][jd1][jd2] = 2;
 422                           fInfocl[1][jd1][jd2] = icl;
 423                           numcell++;
 424                           clust[0][numcell]    = jd1;
 425                           clust[1][numcell]    = jd2;
 426                           cellcount++;
 427                           fInfcl[0][cellcount] = icl;
 428                           fInfcl[1][cellcount] = jd1;
 429                           fInfcl[2][cellcount] = jd2;
 430                         }
 431                     }
 432                 }
 433             }
 434         }
 435     }
 436   return cellcount;
 437 }
 438 // ------------------------------------------------------------------------ //
 439 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
 440 {
 441   // Does the refining of clusters
 442   // Takes the big patch and does gaussian fitting and
 443   // finds out the more refined clusters
 444   //
 445
 446   AliPMDcludata *pmdcludata = 0;
 447
 448   const Int_t kNmaxCell   = 19;    // # of cells surrounding a cluster center
 449
 450   Int_t ndim = incr + 1;
 451
 452   Int_t    *ncl  = 0x0;
 453   Int_t    *clxy = 0x0;
 454   Int_t    i12, i22;
 455   Int_t    i, j, k, i1, i2, id, icl,  itest,ihld, ig, nsupcl,clno, t1, t2;
 456   Float_t  clusdata[6];
 457   Double_t x1, y1, z1, x2, y2, z2, rr;
 458
 459   ncl   = new Int_t [ndim];
 460   clxy  = new Int_t [kNmaxCell];
 461
 462   // Initialisation
 463   for(i = 0; i<ndim; i++)
 464     {
 465       ncl[i]  = -1;
 466       if (i < 6) clusdata[i] = 0.;
 467       if (i < kNmaxCell) clxy[i]    = 0;
 468     }
 469
 470   // clno counts the final clusters
 471   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
 472   // x, y and z store (x,y) coordinates of and energy deposited in a cell
 473   // xc, yc store (x,y) coordinates of the cluster center
 474   // zc stores the energy deposited in a cluster
 475   // rc is cluster radius
 476
 477   clno  = -1;
 478   nsupcl = -1;
 479
 480   for(i = 0; i <= incr; i++)
 481     {
 482       if(fInfcl[0][i] != nsupcl)
 483         {
 484           nsupcl++;
 485         }
 486       if (nsupcl > ndim)
 487         {
 488           AliWarning("RefClust: Too many superclusters!");
 489           nsupcl = ndim;
 490           break;
 491         }
 492       ncl[nsupcl]++;
 493     }
 494
 495   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
 496                   incr+1,nsupcl+1));
 497   id  = -1;
 498   icl = -1;
 499
 500   for(i = 0; i <= nsupcl; i++)
 501     {
 502       if(ncl[i] == 0)
 503         {
 504           id++;
 505           icl++;
 506           if (clno >= 4608)
 507             {
 508               AliWarning("RefClust: Too many clusters! more than 4608");
 509               return;
 510             }
 511           clno++;
 512           i1 = fInfcl[1][id];
 513           i2 = fInfcl[2][id];
 514
 515           i12 = i1 + i2*kNDIMX;
 516
 517           clusdata[0] = fCoord[0][i1][i2];
 518           clusdata[1] = fCoord[1][i1][i2];
 519           clusdata[2] = edepcell[i12];
 520           clusdata[3] = 1.;
 521           clusdata[4] = 0.5;
 522           clusdata[5] = 0.0;
 523
 524           clxy[0] = i1*10000 + i2;
 525
 526           for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
 527             {
 528               clxy[icltr] = -1;
 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]  = pow(z1,2);
 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] +=  pow(z1,2);
 836                   ax2[maxweight]       +=  z1 * pow((x1-x2),2);
 837                   ay2[maxweight]       +=  z1 * pow((y1-y2),2);
 838                   ncell[maxweight]++;
 839
 840                 }
 841             }
 842
 843           for(Int_t kcl = 0; kcl <= ig; kcl++)
 844             {
 845
 846               if(totaladc[kcl]>0){
 847                 if(totaladc[kcl]>0.)xclust[kcl] = (ax[kcl])/ totaladc[kcl];
 848                 if(totaladc[kcl]>0.)yclust[kcl] = (ay[kcl])/ totaladc[kcl];
 849
 850                 if(totaladc[kcl]>0.)sigxclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ax2[kcl];
 851                 if(totaladc[kcl]>0.)sigyclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ay2[kcl];
 852               }
 853
 854               for(j = 0; j < cellCount[kcl]; j++) clno++;
 855
 856               if (clno >= 4608)
 857                 {
 858                   AliWarning("RefClust: Too many clusters! more than 4608");
 859                   return;
 860                 }
 861               clusdata[0] = xclust[kcl];
 862               clusdata[1] = yclust[kcl];
 863               clusdata[2] = totaladc[kcl];
 864               clusdata[3] = ncell[kcl];
 865               if(sigxclust[kcl] > sigyclust[kcl])
 866                 {
 867                   clusdata[4] = pow(sigxclust[kcl],0.5);
 868                   clusdata[5] = pow(sigyclust[kcl],0.5);
 869                 }
 870               else
 871                 {
 872                   clusdata[4] = pow(sigyclust[kcl],0.5);
 873                   clusdata[5] = pow(sigxclust[kcl],0.5);
 874                 }
 875
 876               clxy[0] = tc[kcl];
 877
 878               Int_t Ncell=1;
 879               for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
 880                 {
 881                   if(ii<18)
 882                     {
 883                       clxy[Ncell] = t[cellXY[ii][kcl]];
 884                       Ncell++;
 885                     }
 886                 }
 887
 888               pmdcludata = new AliPMDcludata(clusdata,clxy);
 889               fPMDclucont->Add(pmdcludata);
 890             }
 891
 892           delete [] iord;
 893           delete [] tc;
 894           delete [] t;
 895           delete [] x;
 896           delete [] y;
 897           delete [] z;
 898           delete [] xc;
 899           delete [] yc;
 900           delete [] zc;
 901
 902           delete [] cellCount;
 903           for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
 904
 905           delete [] status;
 906           delete [] totaladc;
 907           delete [] totaladc2;
 908           delete [] ncell;
 909           delete [] xclust;
 910           delete [] yclust;
 911           delete [] sigxclust;
 912           delete [] sigyclust;
 913           delete [] ax;
 914           delete [] ay;
 915           delete [] ax2;
 916           delete [] ay2;
 917           delete [] weight;
 918         }
 919     }
 920   delete [] ncl;
 921   delete [] clxy;
 922 }
 923 // ------------------------------------------------------------------------ //
 924 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
 925                                       Double_t x2, Double_t y2)
 926 {
 927   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
 928 }
 929 // ------------------------------------------------------------------------ //
 930 void AliPMDClusteringV1::FindIsoCell(Int_t idet, Int_t ismn, Double_t celladc[][96], TObjArray *pmdisocell)
 931 {
 932   // Does isolated cell search for offline calibration
 933
 934   AliPMDisocell *isocell = 0;
 935
 936   const Int_t kMaxRow = 48;
 937   const Int_t kMaxCol = 96;
 938   const Int_t kCellNeighbour = 6;
 939
 940   Int_t id1, jd1;
 941
 942   Int_t neibx[6] = {1,0,-1,-1,0,1};
 943   Int_t neiby[6] = {0,1,1,0,-1,-1};
 944
 945
 946   for(Int_t irow = 0; irow < kMaxRow; irow++)
 947     {
 948       for(Int_t icol = 0; icol < kMaxCol; icol++)
 949         {
 950           if(celladc[irow][icol] > 0)
 951             {
 952               Int_t isocount = 0;
 953               for(Int_t ii = 0; ii < kCellNeighbour; ii++)
 954                 {
 955                   id1 = irow + neibx[ii];
 956                   jd1 = icol + neiby[ii];
 957                   Float_t adc = (Float_t) celladc[id1][jd1];
 958                   if(adc == 0.)
 959                     {
 960                       isocount++;
 961                       if(isocount == kCellNeighbour)
 962                         {
 963                           Float_t cadc = (Float_t) celladc[irow][icol];
 964
 965                           isocell = new AliPMDisocell(idet,ismn,irow,icol,cadc);
 966                           pmdisocell->Add(isocell);
 967
 968                         }
 969                     }
 970                 }  // neigh cell cond.
 971             }
 972         }
 973     }
 974
 975
 976 }
 977 // ------------------------------------------------------------------------ //
 978 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
 979 {
 980   fCutoff = decut;
 981 }
 982 // ------------------------------------------------------------------------ //