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   Int_t iord1[kNMX];
 182   TMath::Sort((Int_t)kNMX,edepcell,iord1,jsort);// 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=0,j=0,k=0,id1=0,id2=0,icl=0, numcell=0;
 318   Int_t jd1=0,jd2=0, icell=0, cellcount=0;
 319   Int_t clust[2][kndim];
 320   static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
 321
 322   AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
 323
 324   for (j = 0; j < kNDIMX; j++)
 325     {
 326       for(k = 0; k < kNDIMY; k++)
 327         {
 328           fInfocl[0][j][k] = 0;
 329           fInfocl[1][j][k] = 0;
 330         }
 331     }
 332   for(i=0; i < kNMX; i++)
 333     {
 334       fInfcl[0][i] = -1;
 335
 336       j  = iord1[i];
 337       id2 = j/kNDIMX;
 338       id1 = j-id2*kNDIMX;
 339
 340       if(edepcell[j] <= cutoff)
 341         {
 342           fInfocl[0][id1][id2] = -1;
 343         }
 344     }
 345
 346   // ---------------------------------------------------------------
 347   // crude clustering begins. Start with cell having largest adc
 348   // count and loop over the cells in descending order of adc count
 349   // ---------------------------------------------------------------
 350
 351   icl       = -1;
 352   cellcount = -1;
 353
 354   for(icell = 0; icell <= nmx1; icell++)
 355     {
 356       j  = iord1[icell];
 357       id2 = j/kNDIMX;
 358       id1 = j-id2*kNDIMX;
 359
 360       if(fInfocl[0][id1][id2] == 0 )
 361         {
 362           icl++;
 363           numcell = 0;
 364           cellcount++;
 365           fInfocl[0][id1][id2] = 1;
 366           fInfocl[1][id1][id2] = icl;
 367           fInfcl[0][cellcount] = icl;
 368           fInfcl[1][cellcount] = id1;
 369           fInfcl[2][cellcount] = id2;
 370
 371           clust[0][numcell] = id1;
 372           clust[1][numcell] = id2;
 373
 374           for(i = 1; i < kndim; i++)
 375             {
 376               clust[0][i]=0;
 377             }
 378           // ---------------------------------------------------------------
 379           // check for adc count in neib. cells. If ne 0 put it in this clust
 380           // ---------------------------------------------------------------
 381           for(i = 0; i < 6; i++)
 382             {
 383               jd1 = id1 + neibx[i];
 384               jd2 = id2 + neiby[i];
 385               if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
 386                   fInfocl[0][jd1][jd2] == 0)
 387                 {
 388                   numcell++;
 389                   fInfocl[0][jd1][jd2] = 2;
 390                   fInfocl[1][jd1][jd2] = icl;
 391                   clust[0][numcell]    = jd1;
 392                   clust[1][numcell]    = jd2;
 393                   cellcount++;
 394                   fInfcl[0][cellcount] = icl;
 395                   fInfcl[1][cellcount] = jd1;
 396                   fInfcl[2][cellcount] = jd2;
 397                 }
 398             }
 399           // ---------------------------------------------------------------
 400           // check adc count for neighbour's neighbours recursively and
 401           // if nonzero, add these to the cluster.
 402           // ---------------------------------------------------------------
 403           for(i = 1; i < kndim;i++)
 404             {
 405               if(clust[0][i] != 0)
 406                 {
 407                   id1 = clust[0][i];
 408                   id2 = clust[1][i];
 409                   for(j = 0; j < 6 ; j++)
 410                     {
 411                       jd1 = id1 + neibx[j];
 412                       jd2 = id2 + neiby[j];
 413                       if( (jd1 >= 0 && jd1 < kNDIMX) &&
 414                           (jd2 >= 0 && jd2 < kNDIMY) &&
 415                           fInfocl[0][jd1][jd2] == 0 )
 416                         {
 417                           fInfocl[0][jd1][jd2] = 2;
 418                           fInfocl[1][jd1][jd2] = icl;
 419                           numcell++;
 420                           clust[0][numcell]    = jd1;
 421                           clust[1][numcell]    = jd2;
 422                           cellcount++;
 423                           fInfcl[0][cellcount] = icl;
 424                           fInfcl[1][cellcount] = jd1;
 425                           fInfcl[2][cellcount] = jd2;
 426                         }
 427                     }
 428                 }
 429             }
 430         }
 431     }
 432   return cellcount;
 433 }
 434 // ------------------------------------------------------------------------ //
 435 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
 436 {
 437   // Does the refining of clusters
 438   // Takes the big patch and does gaussian fitting and
 439   // finds out the more refined clusters
 440   //
 441
 442   AliPMDcludata *pmdcludata = 0;
 443
 444   const Int_t kNmaxCell   = 19;    // # of cells surrounding a cluster center
 445
 446   Int_t ndim = incr + 1;
 447
 448   Int_t    *ncl  = 0x0;
 449   Int_t    *clxy = 0x0;
 450   Int_t    i12 = 0, i22 = 0;
 451   Int_t    i = 0, j = 0, k = 0;
 452   Int_t    i1 = 0, i2 = 0, id = 0, icl = 0;
 453   Int_t    itest = 0, ihld = 0, ig = 0;
 454   Int_t    nsupcl = 0, clno = 0, t1 = 0, t2 = 0;
 455   Float_t  clusdata[6];
 456   Double_t x1 = 0, y1 = 0, z1 = 0, x2 = 0, y2 = 0, z2 = 0, rr = 0;
 457
 458   ncl   = new Int_t [ndim];
 459   clxy  = new Int_t [kNmaxCell];
 460
 461   // Initialisation
 462   for(i = 0; i<ndim; i++)
 463     {
 464       ncl[i]  = -1;
 465     }
 466   for(i = 0; i<6; i++)
 467     {
 468       clusdata[i] = 0.;
 469     }
 470   for(i = 0; i<19; i++)
 471     {
 472       clxy[i] = 0;
 473     }
 474
 475   // clno counts the final clusters
 476   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
 477   // x, y and z store (x,y) coordinates of and energy deposited in a cell
 478   // xc, yc store (x,y) coordinates of the cluster center
 479   // zc stores the energy deposited in a cluster
 480   // rc is cluster radius
 481
 482   clno  = -1;
 483   nsupcl = -1;
 484
 485   for(i = 0; i <= incr; i++)
 486     {
 487       if(fInfcl[0][i] != nsupcl)
 488         {
 489           nsupcl++;
 490         }
 491       if (nsupcl > ndim)
 492         {
 493           AliWarning("RefClust: Too many superclusters!");
 494           nsupcl = ndim;
 495           break;
 496         }
 497       ncl[nsupcl]++;
 498     }
 499
 500   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
 501                   incr+1,nsupcl+1));
 502   id  = -1;
 503   icl = -1;
 504
 505   for(i = 0; i <= nsupcl; i++)
 506     {
 507       if(ncl[i] == 0)
 508         {
 509           id++;
 510           icl++;
 511           if (clno >= 4608)
 512             {
 513               AliWarning("RefClust: Too many clusters! more than 4608");
 514               return;
 515             }
 516           clno++;
 517           i1 = fInfcl[1][id];
 518           i2 = fInfcl[2][id];
 519
 520           i12 = i1 + i2*kNDIMX;
 521
 522           clusdata[0] = fCoord[0][i1][i2];
 523           clusdata[1] = fCoord[1][i1][i2];
 524           clusdata[2] = edepcell[i12];
 525           clusdata[3] = 1.;
 526           clusdata[4] = 999.5;
 527           clusdata[5] = 999.5;
 528
 529           clxy[0] = i1*10000 + i2;
 530
 531           for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
 532             {
 533               clxy[icltr] = -1;
 534             }
 535
 536           pmdcludata  = new AliPMDcludata(clusdata,clxy);
 537           fPMDclucont->Add(pmdcludata);
 538         }
 539       else if(ncl[i] == 1)
 540         {
 541           id++;
 542           icl++;
 543           if (clno >= 4608)
 544             {
 545               AliWarning("RefClust: Too many clusters! more than 4608");
 546               return;
 547             }
 548           clno++;
 549           i1   = fInfcl[1][id];
 550           i2   = fInfcl[2][id];
 551           i12  = i1 + i2*kNDIMX;
 552
 553           x1   = fCoord[0][i1][i2];
 554           y1   = fCoord[1][i1][i2];
 555           z1   = edepcell[i12];
 556
 557           clxy[0] = i1*10000 + i2;
 558
 559           id++;
 560           i1   = fInfcl[1][id];
 561           i2   = fInfcl[2][id];
 562
 563           i22  = i1 + i2*kNDIMX;
 564           x2   = fCoord[0][i1][i2];
 565           y2   = fCoord[1][i1][i2];
 566           z2   = edepcell[i22];
 567
 568           clxy[1] = i1*10000 + i2;
 569
 570
 571           for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
 572             {
 573               clxy[icltr] = -1;
 574             }
 575
 576           clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
 577           clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
 578           clusdata[2] = z1+z2;
 579           clusdata[3] = 2.;
 580           clusdata[4] = 0.5;
 581           clusdata[5] = 0.0;
 582           pmdcludata  = new AliPMDcludata(clusdata,clxy);
 583           fPMDclucont->Add(pmdcludata);
 584         }
 585       else
 586         {
 587
 588           Int_t    *iord, *tc, *t;
 589           Double_t *x, *y, *z, *xc, *yc, *zc;
 590
 591           iord = new Int_t [ncl[i]+1];
 592           tc   = new Int_t [ncl[i]+1];
 593           t    = new Int_t [ncl[i]+1];
 594
 595           x    = new Double_t [ncl[i]+1];
 596           y    = new Double_t [ncl[i]+1];
 597           z    = new Double_t [ncl[i]+1];
 598           xc   = new Double_t [ncl[i]+1];
 599           yc   = new Double_t [ncl[i]+1];
 600           zc   = new Double_t [ncl[i]+1];
 601
 602           for( k = 0; k < ncl[i]+1; k++)
 603             {
 604               iord[k] = -1;
 605               t[k]    = -1;
 606               tc[k]   = -1;
 607               x[k]    = -1;
 608               y[k]    = -1;
 609               z[k]    = -1;
 610               xc[k]   = -1;
 611               yc[k]   = -1;
 612               zc[k]   = -1;
 613             }
 614           id++;
 615           // super-cluster of more than two cells - broken up into smaller
 616           // clusters gaussian centers computed. (peaks separated by > 1 cell)
 617           // Begin from cell having largest energy deposited This is first
 618           // cluster center
 619           i1      = fInfcl[1][id];
 620           i2      = fInfcl[2][id];
 621           i12     = i1 + i2*kNDIMX;
 622
 623           x[0]    = fCoord[0][i1][i2];
 624           y[0]    = fCoord[1][i1][i2];
 625           z[0]    = edepcell[i12];
 626           t[0]    = i1*10000 + i2;
 627
 628
 629           iord[0] = 0;
 630           for(j = 1; j <= ncl[i]; j++)
 631             {
 632               id++;
 633               i1      = fInfcl[1][id];
 634               i2      = fInfcl[2][id];
 635               i12     = i1 + i2*kNDIMX;
 636
 637               iord[j] = j;
 638               x[j]    = fCoord[0][i1][i2];
 639               y[j]    = fCoord[1][i1][i2];
 640               z[j]    = edepcell[i12];
 641               t[j]    = i1*10000 + i2;
 642
 643             }
 644
 645           // arranging cells within supercluster in decreasing order
 646
 647           for(j = 1;j <= ncl[i]; j++)
 648             {
 649               itest = 0;
 650               ihld  = iord[j];
 651               for(i1 = 0; i1 < j; i1++)
 652                 {
 653                   if(itest == 0 && z[iord[i1]] < z[ihld])
 654                     {
 655                       itest = 1;
 656                       for(i2 = j-1; i2 >= i1; i2--)
 657                         {
 658                           iord[i2+1] = iord[i2];
 659                         }
 660                       iord[i1] = ihld;
 661                     }
 662                 }
 663             }
 664
 665           if (fClusParam == 1)
 666             {
 667               // Clustering algorithm returns from here for PP collisions
 668               // for pp, only the output of crude clusterng is taken
 669               // sigx and sigy are not calculated at this moment
 670
 671               Double_t supx=0., supy=0., supz=0.;
 672
 673               for(j = 0;j <= ncl[i]; j++)
 674                 {
 675                   supx += x[iord[j]]*z[iord[j]];
 676                   supy += y[iord[j]]*z[iord[j]];
 677                   supz += z[iord[j]];
 678                   if(j < 19)
 679                     {
 680                       clxy[j] = t[iord[j]];
 681                     }
 682                 }
 683
 684               if( ncl[i] + 1 < 19)
 685                 {
 686                   for(Int_t ncel = ncl[i] + 1; ncel < kNmaxCell; ncel ++ )
 687                     {
 688                       clxy[ncel] = -1;
 689                     }
 690                 }
 691               clusdata[0] = supx/supz;
 692               clusdata[1] = supy/supz;
 693               clusdata[2] = supz;
 694               clusdata[3] = ncl[i]+1;
 695               clusdata[4] = 0.5;
 696               clusdata[5] = 0.0;
 697               pmdcludata  = new AliPMDcludata(clusdata,clxy);
 698               fPMDclucont->Add(pmdcludata);
 699             }
 700
 701           /* MODIFICATION PART STARTS (Tapan July 2008)
 702              iord[0] is the cell with highest ADC in the crude-cluster
 703              ig is the number of local maxima in the crude-cluster
 704              For the higest peak we make ig=0 which means first local maximum.
 705              Next we go down in terms of the ADC sequence and find out if any
 706              more of the cells form local maxima. The definition of local
 707              maxima is that all its neighbours are of less ADC compared to it.
 708           */
 709
 710           if (fClusParam == 2)
 711             {
 712               // This part is to split the supercluster
 713               //
 714               ig = 0;
 715               xc[ig] = x[iord[0]];
 716               yc[ig] = y[iord[0]];
 717               zc[ig] = z[iord[0]];
 718               tc[ig] = t[iord[0]];
 719               Int_t ivalid = 0,  icount = 0;
 720
 721               for(j=1;j<=ncl[i];j++)
 722                 {
 723                   x1  = x[iord[j]];
 724                   y1  = y[iord[j]];
 725                   z1  = z[iord[j]];
 726                   t1  = t[iord[j]];
 727                   rr=Distance(x1,y1,xc[ig],yc[ig]);
 728
 729                   // Check the cells which are outside the neighbours (rr>1.2)
 730                   if(rr>1.2 )
 731                     {
 732                       ivalid=0;
 733                       icount=0;
 734                       for(Int_t j1=1;j1<j;j1++)
 735                         {
 736                           icount++;
 737                           Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
 738                           if(rr1>1.2) ivalid++;
 739                         }
 740                       if(ivalid == icount && z1>0.5*zc[ig])
 741                         {
 742                           ig++;
 743                           xc[ig]=x1;
 744                           yc[ig]=y1;
 745                           zc[ig]=z1;
 746                           tc[ig]=t1;
 747                         }
 748                     }
 749                 }
 750
 751               icl=icl+ig+1;
 752
 753               //  We use simple Gaussian weighting. (Tapan Jan 2005)
 754               // compute the number of cells belonging to each cluster.
 755               // cell can be shared between several clusters
 756               // in the ratio of cluster energy deposition
 757               // To calculate:
 758               // (1) number of cells belonging to a cluster (ig) and
 759               // (2) total ADC of the cluster (ig)
 760               // (3) x and y positions of the cluster
 761
 762
 763               Int_t *cellCount;
 764               Int_t **cellXY;
 765
 766               Int_t    *status;
 767               Double_t *totaladc, *totaladc2, *ncell,*weight;
 768               Double_t *xclust, *yclust, *sigxclust, *sigyclust;
 769               Double_t *ax, *ay, *ax2, *ay2;
 770
 771
 772               status    = new Int_t [ncl[i]+1];
 773               cellXY    = new Int_t *[ncl[i]+1];
 774
 775               cellCount = new Int_t [ig+1];
 776               totaladc  = new Double_t [ig+1];
 777               totaladc2 = new Double_t [ig+1];
 778               ncell     = new Double_t [ig+1];
 779               weight    = new Double_t [ig+1];
 780               xclust    = new Double_t [ig+1];
 781               yclust    = new Double_t [ig+1];
 782               sigxclust = new Double_t [ig+1];
 783               sigyclust = new Double_t [ig+1];
 784               ax        = new Double_t [ig+1];
 785               ay        = new Double_t [ig+1];
 786               ax2       = new Double_t [ig+1];
 787               ay2       = new Double_t [ig+1];
 788
 789               for(j = 0; j < ncl[i]+1; j++)
 790                 {
 791                   status[j]     = 0;
 792                   cellXY[j] = new Int_t[ig+1];
 793                 }
 794               //initialization
 795               for(Int_t kcl = 0; kcl < ig+1; kcl++)
 796                 {
 797                   cellCount[kcl] = 0;
 798                   totaladc[kcl]  = 0.;
 799                   totaladc2[kcl] = 0.;
 800                   ncell[kcl]     = 0.;
 801                   weight[kcl]    = 0.;
 802                   xclust[kcl]    = 0.;
 803                   yclust[kcl]    = 0.;
 804                   sigxclust[kcl] = 0.;
 805                   sigyclust[kcl] = 0.;
 806                   ax[kcl]        = 0.;
 807                   ay[kcl]        = 0.;
 808                   ax2[kcl]       = 0.;
 809                   ay2[kcl]       = 0.;
 810                   for(j = 0; j < ncl[i]+1; j++)
 811                     {
 812                       cellXY[j][kcl] = 0;
 813                     }
 814                 }
 815               Double_t sumweight, gweight;
 816
 817               for(j = 0;j <= ncl[i]; j++)
 818                 {
 819                   x1 = x[iord[j]];
 820                   y1 = y[iord[j]];
 821                   z1 = z[iord[j]];
 822                   t1 = t[iord[j]];
 823
 824                   for(Int_t kcl=0; kcl<=ig; kcl++)
 825                     {
 826                       x2 = xc[kcl];
 827                       y2 = yc[kcl];
 828                       rr = Distance(x1,y1,x2,y2);
 829                       t2 = tc[kcl];
 830
 831                       if(rr==0)
 832                         {
 833                           ncell[kcl]     = 1.;
 834                           totaladc[kcl]  = z1;
 835                           totaladc2[kcl] = z1*z1;
 836                           ax[kcl]        = x1 * z1;
 837                           ay[kcl]        = y1 * z1;
 838                           ax2[kcl]       = 0.;
 839                           ay2[kcl]       = 0.;
 840                           status[j]      = 1;
 841                         }
 842                     }
 843                 }
 844
 845               for(j = 0; j <= ncl[i]; j++)
 846                 {
 847                   Int_t   maxweight = 0;
 848                   Double_t     max  = 0.;
 849
 850                   if(status[j] == 0)
 851                     {
 852                       x1 = x[iord[j]];
 853                       y1 = y[iord[j]];
 854                       z1 = z[iord[j]];
 855                       t1 = t[iord[j]];
 856                       sumweight = 0.;
 857
 858                       for(Int_t kcl = 0; kcl <= ig; kcl++)
 859                         {
 860                           x2 = xc[kcl];
 861                           y2 = yc[kcl];
 862                           rr = Distance(x1,y1,x2,y2);
 863                           gweight     = exp(-(rr*rr)/(2*(1.2*1.2)));
 864                           weight[kcl] = zc[kcl] * gweight;
 865                           sumweight   = sumweight + weight[kcl];
 866
 867                           if(weight[kcl] > max)
 868                             {
 869                               max       =  weight[kcl];
 870                               maxweight =  kcl;
 871                             }
 872                         }
 873
 874                       cellXY[cellCount[maxweight]][maxweight] = iord[j];
 875
 876                       cellCount[maxweight]++;
 877
 878                       x2 = xc[maxweight];
 879                       y2 = yc[maxweight];
 880                       totaladc[maxweight]  +=  z1;
 881                       ax[maxweight]        +=  x1*z1;
 882                       ay[maxweight]        +=  y1*z1;
 883                       totaladc2[maxweight] +=  z1*z1;
 884                       ax2[maxweight]       +=  z1*(x1-x2)*(x1-x2);
 885                       ay2[maxweight]       +=  z1*(y1-y2)*(y1-y2);
 886                       ncell[maxweight]++;
 887
 888                     }
 889                 }
 890
 891               for(Int_t kcl = 0; kcl <= ig; kcl++)
 892                 {
 893                   if(totaladc[kcl] > 0.)
 894                     {
 895                       xclust[kcl] = (ax[kcl])/ totaladc[kcl];
 896                       yclust[kcl] = (ay[kcl])/ totaladc[kcl];
 897
 898                       //natasha
 899                       Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
 900                       if(totaladc2[kcl] >= sqtotadc)
 901                         {
 902                           sigxclust[kcl] = 0.25;
 903                           sigyclust[kcl] = 0.25;
 904                         }
 905                       else
 906                         {
 907                           sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
 908                           sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
 909                         }
 910                     }
 911
 912                   for(j = 0; j < cellCount[kcl]; j++) clno++;
 913
 914                   if (clno >= 4608)
 915                     {
 916                       AliWarning("RefClust: Too many clusters! more than 4608");
 917                       return;
 918                     }
 919                   clusdata[0] = xclust[kcl];
 920                   clusdata[1] = yclust[kcl];
 921                   clusdata[2] = totaladc[kcl];
 922                   clusdata[3] = ncell[kcl];
 923
 924                   if(sigxclust[kcl] > sigyclust[kcl])
 925                     {
 926                       clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
 927                       clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
 928                     }
 929                   else
 930                     {
 931                       clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
 932                       clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
 933                     }
 934
 935                   clxy[0] = tc[kcl];
 936
 937                   Int_t Ncell=1;
 938                   for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
 939                     {
 940                       if(ii<18)
 941                         {
 942                           clxy[Ncell] = t[cellXY[ii][kcl]];
 943                           Ncell++;
 944                         }
 945                     }
 946
 947                   pmdcludata = new AliPMDcludata(clusdata,clxy);
 948                   fPMDclucont->Add(pmdcludata);
 949                 }
 950               delete [] cellCount;
 951               for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
 952
 953               delete [] status;
 954               delete [] totaladc;
 955               delete [] totaladc2;
 956               delete [] ncell;
 957               delete [] xclust;
 958               delete [] yclust;
 959               delete [] sigxclust;
 960               delete [] sigyclust;
 961               delete [] ax;
 962               delete [] ay;
 963               delete [] ax2;
 964               delete [] ay2;
 965               delete [] weight;
 966
 967             }
 968
 969           delete [] iord;
 970           delete [] tc;
 971           delete [] t;
 972           delete [] x;
 973           delete [] y;
 974           delete [] z;
 975           delete [] xc;
 976           delete [] yc;
 977           delete [] zc;
 978
 979
 980         }
 981     }
 982   delete [] ncl;
 983   delete [] clxy;
 984 }
 985 // ------------------------------------------------------------------------ //
 986 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
 987                                       Double_t x2, Double_t y2)
 988 {
 989   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
 990 }
 991 // ------------------------------------------------------------------------ //
 992 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
 993 {
 994   fCutoff = decut;
 995 }
 996 // ------------------------------------------------------------------------ //
 997 void AliPMDClusteringV1::SetClusteringParam(Int_t cluspar)
 998 {
 999   fClusParam = cluspar;
1000 }
1001 // ------------------------------------------------------------------------ //