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
 655           if (fClusParam == 1)
 656             {
 657               // Clustering algorithm returns from here for PP collisions
 658               // for pp, only the output of crude clusterng is taken
 659               // sigx and sigy are not calculated at this moment
 660
 661               Double_t supx=0., supy=0., supz=0.;
 662
 663               for(j = 0;j <= ncl[i]; j++)
 664                 {
 665                   supx += x[iord[j]]*z[iord[j]];
 666                   supy += y[iord[j]]*z[iord[j]];
 667                   supz += z[iord[j]];
 668                   if(j < 19)
 669                     {
 670                       clxy[j] = t[iord[j]];
 671                     }
 672                 }
 673
 674               if( ncl[i] + 1 < 19)
 675                 {
 676                   for(Int_t ncel = ncl[i] + 1; ncel < kNmaxCell; ncel ++ )
 677                     {
 678                       clxy[ncel] = -1;
 679                     }
 680                 }
 681               clusdata[0] = supx/supz;
 682               clusdata[1] = supy/supz;
 683               clusdata[2] = supz;
 684               clusdata[3] = ncl[i]+1;
 685               clusdata[4] = 0.5;
 686               clusdata[5] = 0.0;
 687               pmdcludata  = new AliPMDcludata(clusdata,clxy);
 688               fPMDclucont->Add(pmdcludata);
 689             }
 690
 691           /* MODIFICATION PART STARTS (Tapan July 2008)
 692              iord[0] is the cell with highest ADC in the crude-cluster
 693              ig is the number of local maxima in the crude-cluster
 694              For the higest peak we make ig=0 which means first local maximum.
 695              Next we go down in terms of the ADC sequence and find out if any
 696              more of the cells form local maxima. The definition of local
 697              maxima is that all its neighbours are of less ADC compared to it.
 698           */
 699
 700           if (fClusParam == 2)
 701             {
 702               // This part is to split the supercluster
 703               //
 704               ig = 0;
 705               xc[ig] = x[iord[0]];
 706               yc[ig] = y[iord[0]];
 707               zc[ig] = z[iord[0]];
 708               tc[ig] = t[iord[0]];
 709               Int_t ivalid = 0,  icount = 0;
 710
 711               for(j=1;j<=ncl[i];j++)
 712                 {
 713                   x1  = x[iord[j]];
 714                   y1  = y[iord[j]];
 715                   z1  = z[iord[j]];
 716                   t1  = t[iord[j]];
 717                   rr=Distance(x1,y1,xc[ig],yc[ig]);
 718
 719                   // Check the cells which are outside the neighbours (rr>1.2)
 720                   if(rr>1.2 )
 721                     {
 722                       ivalid=0;
 723                       icount=0;
 724                       for(Int_t j1=1;j1<j;j1++)
 725                         {
 726                           icount++;
 727                           Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
 728                           if(rr1>1.2) ivalid++;
 729                         }
 730                       if(ivalid == icount && z1>0.5*zc[ig])
 731                         {
 732                           ig++;
 733                           xc[ig]=x1;
 734                           yc[ig]=y1;
 735                           zc[ig]=z1;
 736                           tc[ig]=t1;
 737                         }
 738                     }
 739                 }
 740
 741               icl=icl+ig+1;
 742
 743               //  We use simple Gaussian weighting. (Tapan Jan 2005)
 744               // compute the number of cells belonging to each cluster.
 745               // cell can be shared between several clusters
 746               // in the ratio of cluster energy deposition
 747               // To calculate:
 748               // (1) number of cells belonging to a cluster (ig) and
 749               // (2) total ADC of the cluster (ig)
 750               // (3) x and y positions of the cluster
 751
 752
 753               Int_t *cellCount;
 754               Int_t **cellXY;
 755
 756               Int_t    *status;
 757               Double_t *totaladc, *totaladc2, *ncell,*weight;
 758               Double_t *xclust, *yclust, *sigxclust, *sigyclust;
 759               Double_t *ax, *ay, *ax2, *ay2;
 760
 761
 762               status    = new Int_t [ncl[i]+1];
 763               cellXY    = new Int_t *[ncl[i]+1];
 764
 765               cellCount = new Int_t [ig+1];
 766               totaladc  = new Double_t [ig+1];
 767               totaladc2 = new Double_t [ig+1];
 768               ncell     = new Double_t [ig+1];
 769               weight    = new Double_t [ig+1];
 770               xclust    = new Double_t [ig+1];
 771               yclust    = new Double_t [ig+1];
 772               sigxclust = new Double_t [ig+1];
 773               sigyclust = new Double_t [ig+1];
 774               ax        = new Double_t [ig+1];
 775               ay        = new Double_t [ig+1];
 776               ax2       = new Double_t [ig+1];
 777               ay2       = new Double_t [ig+1];
 778
 779               for(j = 0; j < ncl[i]+1; j++)
 780                 {
 781                   status[j]     = 0;
 782                   cellXY[j] = new Int_t[ig+1];
 783                 }
 784               //initialization
 785               for(Int_t kcl = 0; kcl < ig+1; kcl++)
 786                 {
 787                   cellCount[kcl] = 0;
 788                   totaladc[kcl]  = 0.;
 789                   totaladc2[kcl] = 0.;
 790                   ncell[kcl]     = 0.;
 791                   weight[kcl]    = 0.;
 792                   xclust[kcl]    = 0.;
 793                   yclust[kcl]    = 0.;
 794                   sigxclust[kcl] = 0.;
 795                   sigyclust[kcl] = 0.;
 796                   ax[kcl]        = 0.;
 797                   ay[kcl]        = 0.;
 798                   ax2[kcl]       = 0.;
 799                   ay2[kcl]       = 0.;
 800                   for(j = 0; j < ncl[i]+1; j++)
 801                     {
 802                       cellXY[j][kcl] = 0;
 803                     }
 804                 }
 805               Double_t sumweight, gweight;
 806
 807               for(j = 0;j <= ncl[i]; j++)
 808                 {
 809                   x1 = x[iord[j]];
 810                   y1 = y[iord[j]];
 811                   z1 = z[iord[j]];
 812                   t1 = t[iord[j]];
 813
 814                   for(Int_t kcl=0; kcl<=ig; kcl++)
 815                     {
 816                       x2 = xc[kcl];
 817                       y2 = yc[kcl];
 818                       rr = Distance(x1,y1,x2,y2);
 819                       t2 = tc[kcl];
 820
 821                       if(rr==0)
 822                         {
 823                           ncell[kcl]     = 1.;
 824                           totaladc[kcl]  = z1;
 825                           totaladc2[kcl] = z1*z1;
 826                           ax[kcl]        = x1 * z1;
 827                           ay[kcl]        = y1 * z1;
 828                           ax2[kcl]       = 0.;
 829                           ay2[kcl]       = 0.;
 830                           status[j]      = 1;
 831                         }
 832                     }
 833                 }
 834
 835               for(j = 0; j <= ncl[i]; j++)
 836                 {
 837                   Int_t   maxweight = 0;
 838                   Double_t     max  = 0.;
 839
 840                   if(status[j] == 0)
 841                     {
 842                       x1 = x[iord[j]];
 843                       y1 = y[iord[j]];
 844                       z1 = z[iord[j]];
 845                       t1 = t[iord[j]];
 846                       sumweight = 0.;
 847
 848                       for(Int_t kcl = 0; kcl <= ig; kcl++)
 849                         {
 850                           x2 = xc[kcl];
 851                           y2 = yc[kcl];
 852                           rr = Distance(x1,y1,x2,y2);
 853                           gweight     = exp(-(rr*rr)/(2*(1.2*1.2)));
 854                           weight[kcl] = zc[kcl] * gweight;
 855                           sumweight   = sumweight + weight[kcl];
 856
 857                           if(weight[kcl] > max)
 858                             {
 859                               max       =  weight[kcl];
 860                               maxweight =  kcl;
 861                             }
 862                         }
 863
 864                       cellXY[cellCount[maxweight]][maxweight] = iord[j];
 865
 866                       cellCount[maxweight]++;
 867
 868                       x2 = xc[maxweight];
 869                       y2 = yc[maxweight];
 870                       totaladc[maxweight]  +=  z1;
 871                       ax[maxweight]        +=  x1*z1;
 872                       ay[maxweight]        +=  y1*z1;
 873                       totaladc2[maxweight] +=  z1*z1;
 874                       ax2[maxweight]       +=  z1*(x1-x2)*(x1-x2);
 875                       ay2[maxweight]       +=  z1*(y1-y2)*(y1-y2);
 876                       ncell[maxweight]++;
 877
 878                     }
 879                 }
 880
 881               for(Int_t kcl = 0; kcl <= ig; kcl++)
 882                 {
 883                   if(totaladc[kcl] > 0.)
 884                     {
 885                       xclust[kcl] = (ax[kcl])/ totaladc[kcl];
 886                       yclust[kcl] = (ay[kcl])/ totaladc[kcl];
 887
 888                       //natasha
 889                       Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
 890                       if(totaladc2[kcl] >= sqtotadc)
 891                         {
 892                           sigxclust[kcl] = 0.25;
 893                           sigyclust[kcl] = 0.25;
 894                         }
 895                       else
 896                         {
 897                           sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
 898                           sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
 899                         }
 900                     }
 901
 902                   for(j = 0; j < cellCount[kcl]; j++) clno++;
 903
 904                   if (clno >= 4608)
 905                     {
 906                       AliWarning("RefClust: Too many clusters! more than 4608");
 907                       return;
 908                     }
 909                   clusdata[0] = xclust[kcl];
 910                   clusdata[1] = yclust[kcl];
 911                   clusdata[2] = totaladc[kcl];
 912                   clusdata[3] = ncell[kcl];
 913
 914                   if(sigxclust[kcl] > sigyclust[kcl])
 915                     {
 916                       clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
 917                       clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
 918                     }
 919                   else
 920                     {
 921                       clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
 922                       clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
 923                     }
 924
 925                   clxy[0] = tc[kcl];
 926
 927                   Int_t Ncell=1;
 928                   for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
 929                     {
 930                       if(ii<18)
 931                         {
 932                           clxy[Ncell] = t[cellXY[ii][kcl]];
 933                           Ncell++;
 934                         }
 935                     }
 936
 937                   pmdcludata = new AliPMDcludata(clusdata,clxy);
 938                   fPMDclucont->Add(pmdcludata);
 939                 }
 940               delete [] cellCount;
 941               for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
 942
 943               delete [] status;
 944               delete [] totaladc;
 945               delete [] totaladc2;
 946               delete [] ncell;
 947               delete [] xclust;
 948               delete [] yclust;
 949               delete [] sigxclust;
 950               delete [] sigyclust;
 951               delete [] ax;
 952               delete [] ay;
 953               delete [] ax2;
 954               delete [] ay2;
 955               delete [] weight;
 956
 957             }
 958
 959           delete [] iord;
 960           delete [] tc;
 961           delete [] t;
 962           delete [] x;
 963           delete [] y;
 964           delete [] z;
 965           delete [] xc;
 966           delete [] yc;
 967           delete [] zc;
 968
 969
 970         }
 971     }
 972   delete [] ncl;
 973   delete [] clxy;
 974 }
 975 // ------------------------------------------------------------------------ //
 976 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
 977                                       Double_t x2, Double_t y2)
 978 {
 979   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
 980 }
 981 // ------------------------------------------------------------------------ //
 982 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
 983 {
 984   fCutoff = decut;
 985 }
 986 // ------------------------------------------------------------------------ //
 987 void AliPMDClusteringV1::SetClusteringParam(Int_t cluspar)
 988 {
 989   fClusParam = cluspar;
 990 }
 991 // ------------------------------------------------------------------------ //