PMD/AliPMDClusteringV1.cxx

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