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