Another histos for lumi
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV1.cxx
index 2a15ac3..d9a4f9a 100644 (file)
@@ -13,6 +13,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 //-----------------------------------------------------//
 //                                                     //
 //  Source File : PMDClusteringV1.cxx, Version 00      //
    Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
    ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
    builds up superclusters and breaks them into clusters. The input is
-   in array fEdepCell[kNDIMX][kNDIMY] and cluster information is in array
-   fClusters[5][5000]. integer fClno gives total number of clusters in the
+   in array edepcell[kNMX] and cluster information is in a
+   TObjarray. Integer clno gives total number of clusters in the
    supermodule.
 
-   fEdepCell, fClno  and fClusters are the only global ( public ) variables.
+   fClusters is the only global ( public ) variables.
    Others are local ( private ) to the code.
    At the moment, the data is read for whole detector ( all supermodules
    and pmd as well as cpv. This will have to be modify later )
    LAST UPDATE  :  October 23, 2002
 -----------------------------------------------------------------------*/
 
-#include "Riostream.h"
+#include <Riostream.h>
+#include <TMath.h>
 #include <TNtuple.h>
 #include <TObjArray.h>
+#include "TRandom.h"
 #include <stdio.h>
 
+#include "AliPMDcludata.h"
 #include "AliPMDcluster.h"
 #include "AliPMDClustering.h"
 #include "AliPMDClusteringV1.h"
 #include "AliLog.h"
 
+
 ClassImp(AliPMDClusteringV1)
 
 const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254;  // sqrt(3.)/2.
 
 AliPMDClusteringV1::AliPMDClusteringV1():
-  fCutoff(0.0)
+  fPMDclucont(new TObjArray()),
+  fCutoff(0.0),
+  fClusParam(0)
 {
-  for(int i = 0; i < kNDIMX; i++)
+  for(Int_t i = 0; i < kNDIMX; i++)
     {
-      for(int j = 0; j < kNDIMY; j++)
+      for(Int_t j = 0; j < kNDIMY; j++)
        {
          fCoord[0][i][j] = i+j/2.;
          fCoord[1][i][j] = fgkSqroot3by2*j;
-         fEdepCell[i][j] = 0;
        }
     }
 }
 // ------------------------------------------------------------------------ //
+AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
+  AliPMDClustering(pmdclv1),
+  fPMDclucont(0),
+  fCutoff(0),
+  fClusParam(0)
+{
+  // copy constructor
+  AliError("Copy constructor not allowed ");
+  
+}
+// ------------------------------------------------------------------------ //
+AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
+{
+  // copy constructor
+  AliError("Assignment operator not allowed ");
+  return *this;
+}
+// ------------------------------------------------------------------------ //
 AliPMDClusteringV1::~AliPMDClusteringV1()
 {
-
+  delete fPMDclucont;
 }
 // ------------------------------------------------------------------------ //
-void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96], TObjArray *pmdcont)
+void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
+                                Int_t celltrack[48][96],
+                                Int_t cellpid[48][96],
+                                Double_t celladc[48][96],
+                                TObjArray *pmdcont)
 {
   // main function to call other necessary functions to do clustering
   //
-  AliPMDcluster *pmdcl = 0;
-  /*
-    int id and jd defined to read the input data.
-    It is assumed that for data we have 0 <= id <= 48
-    and 0 <= jd <=96
-  */
-
-  int i, i1, i2, j, nmx1, incr, id, jd;
-  Int_t   celldataX[15], celldataY[15];
-  Float_t clusdata[6];
 
-  Double_t  cutoff, ave;
+  AliPMDcluster *pmdcl = 0;
 
   const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
+  const Int_t kNmaxCell   = 19;     // # of cells surrounding a cluster center
 
+  Int_t    i = 0,  j = 0, nmx1 = 0;
+  Int_t    incr = 0, id = 0, jd = 0;
+  Int_t    celldataX[kNmaxCell], celldataY[kNmaxCell];
+  Int_t    celldataTr[kNmaxCell], celldataPid[kNmaxCell];
+  Float_t  celldataAdc[kNmaxCell];
+  Float_t  clusdata[6] = {0.,0.,0.,0.,0.,0.};
+  Double_t cutoff, ave;
+  Double_t edepcell[kNMX];
+  Double_t cellenergy[kNMX];
+  
   // ndimXr and ndimYr are different because of different module size
 
-  Int_t ndimXr =0;
-  Int_t ndimYr =0;
+  Int_t ndimXr = 0;
+  Int_t ndimYr = 0;
 
   if (ismn < 12)
     {
@@ -107,12 +137,18 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96
       ndimYr = 96;
     }
 
-  for (Int_t i =0; i < kNDIMX; i++)
+  for (i = 0; i < kNMX; i++)
+  {
+      cellenergy[i] = 0.;
+  }
+
+  Int_t kk = 0;
+  for (i = 0; i < kNDIMX; i++)
     {
-      for (Int_t j =0; j < kNDIMY; j++)
+      for (j = 0; j < kNDIMY; j++)
        {
-         fEdepCell[i][j] = 0;
-         fCellTrNo[i][j] = -1;
+         edepcell[kk] = 0.;
+         kk++;
        }
     }
 
@@ -120,68 +156,82 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96
     {
       for (jd = 0; jd < ndimYr; jd++)
        {
-         j=jd;
-         i=id+(ndimYr/2-1)-(jd/2);
+         j = jd;
+         i = id+(ndimYr/2-1)-(jd/2);
 
+         Int_t ij = i + j*kNDIMX;
+         
          if (ismn < 12)
            {
-             fEdepCell[i][j] = celladc[jd][id];
-             fCellTrNo[i][j] = jd*10000+id; /* for association */
+             cellenergy[ij]    = celladc[jd][id];
            }
          else if (ismn >= 12 && ismn <= 23)
            {
-             fEdepCell[i][j] = celladc[id][jd];
-             fCellTrNo[i][j] = id*10000+jd; /* for association */
+             cellenergy[ij]    = celladc[id][jd];
            }
-
        }
     }
-  Order(); // order the data
-  cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
-  ave=0.;
-  nmx1=-1;
-
-  for(j=0;j<kNMX; j++)
+  
+  for (i = 0; i < kNMX; i++)
     {
-      i1 = fIord[0][j];
-      i2 = fIord[1][j];
-      if (fEdepCell[i1][i2] > 0.) {ave = ave + fEdepCell[i1][i2];}
-      if (fEdepCell[i1][i2] > cutoff ) nmx1 = nmx1 + 1;
+      edepcell[i] = cellenergy[i];
     }
 
+  Bool_t jsort = true;
+  // the dimension of iord1 is increased twice
+  Int_t iord1[2*kNMX];
+  TMath::Sort((Int_t)kNMX,edepcell,iord1,jsort);// order the data
+  cutoff = fCutoff;                             // cutoff to discard cells
+  ave  = 0.;
+  nmx1 = -1;
+  for(i = 0;i < kNMX; i++)
+    {
+      if(edepcell[i] > 0.) 
+       {
+         ave += edepcell[i];
+       }
+      if(edepcell[i] > cutoff )
+       {
+         nmx1++;
+       }
+    }
+  
   AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
 
   if (nmx1 == 0) nmx1 = 1;
-  ave=ave/nmx1;
-
+  ave = ave/nmx1;
   AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
                  kNMX,ave));
-          
-  incr = CrClust(ave, cutoff, nmx1);
-  RefClust(incr);
-
-  AliDebug(1,Form("Detector Plane = %d  Serial Module No = %d Number of clusters = %d",idet, ismn, fClno));
-
   
-  for(i1=0; i1<=fClno; i1++)
+  incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
+  RefClust(incr,edepcell);
+  Int_t nentries1 = fPMDclucont->GetEntries();
+  AliDebug(1,Form("Detector Plane = %d  Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
+  AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
+  
+  for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
     {
-      Float_t cluXC    = (Float_t) fClusters[0][i1];
-      Float_t cluYC    = (Float_t) fClusters[1][i1];
-      Float_t cluADC   = (Float_t) fClusters[2][i1];
-      Float_t cluCELLS = (Float_t) fClusters[3][i1];
-      Float_t cluRAD   = (Float_t) fClusters[4][i1];
+      AliPMDcludata *pmdcludata = 
+       (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
+      Float_t cluXC    = pmdcludata->GetClusX();
+      Float_t cluYC    = pmdcludata->GetClusY();
+      Float_t cluADC   = pmdcludata->GetClusADC();
+      Float_t cluCELLS = pmdcludata->GetClusCells();
+      Float_t cluSIGX  = pmdcludata->GetClusSigmaX();
+      Float_t cluSIGY  = pmdcludata->GetClusSigmaY();
+      
       Float_t cluY0    = ktwobysqrt3*cluYC;
       Float_t cluX0    = cluXC - cluY0/2.;
-
-
+      
       // 
       // Cluster X centroid is back transformed
       //
+
       if (ismn < 12)
        {
          clusdata[0] = cluX0 - (24-1) + cluY0/2.;
        }
-      else if (ismn >= 12 && ismn <= 23)
+      else if ( ismn >= 12 && ismn <= 23)
        {
          clusdata[0] = cluX0 - (48-1) + cluY0/2.;
        }         
@@ -189,675 +239,799 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96
       clusdata[1]     = cluY0;
       clusdata[2]     = cluADC;
       clusdata[3]     = cluCELLS;
-      clusdata[4]     = cluRAD;
-      clusdata[5]     = 0.;
+      clusdata[4]     = cluSIGX;
+      clusdata[5]     = cluSIGY;
 
       //
       // Cells associated with a cluster
       //
-      for (Int_t ihit = 0; ihit < 15; ihit++)
+
+      for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
        {
+         Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
+         Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
 
          if (ismn < 12)
            {
-             celldataX[ihit] = fClTr[ihit][i1]%10000;
-             celldataY[ihit] = fClTr[ihit][i1]/10000;
+             celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
            }
          else if (ismn >= 12 && ismn <= 23)
            {
-             celldataX[ihit] = fClTr[ihit][i1]/10000;
-             celldataY[ihit] = fClTr[ihit][i1]%10000;
+             celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
            }
+         
+         celldataY[ihit] = cellcol;
+         
+         Int_t irow = celldataX[ihit];
+         Int_t icol = celldataY[ihit];
+
+         if (ismn < 12)
+           {
+             if ((irow >= 0 && irow < 96) && (icol >= 0 && icol < 48))
+               {
+                 celldataTr[ihit]  = celltrack[icol][irow];
+                 celldataPid[ihit] = cellpid[icol][irow];
+                 celldataAdc[ihit] = (Float_t) celladc[icol][irow];
+               }
+             else
+               {
+                 celldataTr[ihit]  = -1;
+                 celldataPid[ihit] = -1;
+                 celldataAdc[ihit] = -1;
+               }
+           }
+         else if (ismn >= 12 && ismn < 24)
+           {
+             if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
+               {
+                 celldataTr[ihit]  = celltrack[irow][icol];
+                 celldataPid[ihit] = cellpid[irow][icol];
+                 celldataAdc[ihit] = (Float_t) celladc[irow][icol];
+
+               }
+             else
+               {
+                 celldataTr[ihit]  = -1;
+                 celldataPid[ihit] = -1;
+                 celldataAdc[ihit] = -1;
+               }
+           }
+         
        }
-      //printf("%d %f %f\n",idet,cluXC,cluYC );
-      pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
+      
+      pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
+                               celldataTr, celldataPid, celldataAdc);
       pmdcont->Add(pmdcl);
     }
+  
+  fPMDclucont->Delete();
 }
 // ------------------------------------------------------------------------ //
-void AliPMDClusteringV1::Order()
+Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
+                                 Int_t iord1[], Double_t edepcell[])
 {
-  // Sorting algorithm
-  // sorts the ADC values from higher to lower
-  //
-  double dd[kNMX];
-  // matrix fEdepCell converted into
-  // one dimensional array dd. adum a place holder for double
-  int i, j, i1, i2, iord1[kNMX];
-  // information of
-  // ordering is stored in iord1, original array not ordered
+  // Does crude clustering 
+  // Finds out only the big patch by just searching the
+  // connected cells
   //
-  // define arrays dd and iord1
-  for(i1=0; i1 < kNDIMX; i1++)
+  const Int_t kndim = 4609;
+  Int_t i=0,j=0,k=0,id1=0,id2=0,icl=0, numcell=0;
+  Int_t jd1=0,jd2=0, icell=0, cellcount=0;
+  Int_t clust[2][kndim];
+  static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
+
+  AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
+  
+  for (j = 0; j < kNDIMX; j++)
     {
-      for(i2=0; i2 < kNDIMY; i2++)
+      for(k = 0; k < kNDIMY; k++)
        {
-         i        = i1 + i2*kNDIMX;
-         iord1[i] = i;
-         dd[i]    = fEdepCell[i1][i2];
+         fInfocl[0][j][k] = 0;
+         fInfocl[1][j][k] = 0;
        }
     }
-  
-  TMath::Sort(kNMX,dd,iord1); //PH Using much better algorithm...
-  // store the sorted information in fIord for later use
-  for(i=0; i<kNMX; i++)
+  for(i=0; i < kNMX; i++)
     {
+      fInfcl[0][i] = -1;
+      
       j  = iord1[i];
-      i2 = j/kNDIMX;
-      i1 = j-i2*kNDIMX;
-      fIord[0][i]=i1;
-      fIord[1][i]=i2;
-    }
-}
-// ------------------------------------------------------------------------ //
-int AliPMDClusteringV1::CrClust(double ave, double cutoff, int nmx1)
-{
-  // Does crude clustering
-  // Finds out only the big patch by just searching the
-  // connected cells
-  //
-  int i,j,k,id1,id2,icl, numcell, clust[2][5000];
-  int jd1,jd2, icell, cellcount;
-  static int neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
-  // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
-  // cell. There are six neighbours.
-  // cellcount --- total number of cells having nonzero ener dep
-  // numcell --- number of cells in a given supercluster
-  // ofstream ofl0("cells_loc",ios::out);
-  // initialize fInfocl[2][kNDIMX][kNDIMY]
+      id2 = j/kNDIMX;
+      id1 = j-id2*kNDIMX;
 
-  AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
-  
-  for (j=0; j < kNDIMX; j++){
-    for(k=0; k < kNDIMY; k++){
-      fInfocl[0][j][k] = 0;
-      fInfocl[1][j][k] = 0;
+      if(edepcell[j] <= cutoff)
+       {
+         fInfocl[0][id1][id2] = -1;
+       }
     }
-  }
-  for(i=0; i < kNMX; i++){
-    fInfcl[0][i] = -1;
-    id1=fIord[0][i];
-    id2=fIord[1][i];
-    if(fEdepCell[id1][id2] <= cutoff){fInfocl[0][id1][id2]=-1;}
-  }
+
   // ---------------------------------------------------------------
   // crude clustering begins. Start with cell having largest adc
   // count and loop over the cells in descending order of adc count
   // ---------------------------------------------------------------
-  icl=-1;
-  cellcount=-1;
-  for(icell=0; icell <= nmx1; icell++){
-    id1=fIord[0][icell];
-    id2=fIord[1][icell];
-    if(fInfocl[0][id1][id2] == 0 ){
-      // ---------------------------------------------------------------
-      // icl -- cluster #, numcell -- # of cells in it, clust -- stores
-      // coordinates of the cells in a cluster, fInfocl[0][i1][i2] is 1 for
-      // primary and 2 for secondary cells,
-      // fInfocl[1][i1][i2] stores cluster #
-      // ---------------------------------------------------------------
-      icl=icl+1;
-      numcell=0;
-      cellcount = cellcount + 1;
-      fInfocl[0][id1][id2]=1;
-      fInfocl[1][id1][id2]=icl;
-      fInfcl[0][cellcount]=icl;
-      fInfcl[1][cellcount]=id1;
-      fInfcl[2][cellcount]=id2;
-
-      clust[0][numcell]=id1;
-      clust[1][numcell]=id2;
-      for(i=1; i<5000; i++)clust[0][i]=0;
-      // ---------------------------------------------------------------
-      // check for adc count in neib. cells. If ne 0 put it in this clust
-      // ---------------------------------------------------------------
-      for(i=0; i<6; i++){
-       jd1=id1+neibx[i];
-       jd2=id2+neiby[i];
-       if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
-           fInfocl[0][jd1][jd2] == 0){
-         numcell=numcell+1;
-         fInfocl[0][jd1][jd2]=2;
-         fInfocl[1][jd1][jd2]=icl;
-         clust[0][numcell]=jd1;
-         clust[1][numcell]=jd2;
-         cellcount=cellcount+1;
-         fInfcl[0][cellcount]=icl;
-         fInfcl[1][cellcount]=jd1;
-         fInfcl[2][cellcount]=jd2;
-       }
-      }
-      // ---------------------------------------------------------------
-      // check adc count for neighbour's neighbours recursively and
-      // if nonzero, add these to the cluster.
-      // ---------------------------------------------------------------
-      for(i=1;i < 5000;i++){
-       if(clust[0][i] != 0){
-         id1=clust[0][i];
-         id2=clust[1][i];
-         for(j=0; j<6 ; j++){
-           jd1=id1+neibx[j];
-           jd2=id2+neiby[j];
-           if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
-               fInfocl[0][jd1][jd2] == 0 ){
-             fInfocl[0][jd1][jd2] = 2;
-             fInfocl[1][jd1][jd2] = icl;
-             numcell              = numcell + 1;
-             clust[0][numcell]    = jd1;
-             clust[1][numcell]    = jd2;
-             cellcount            = cellcount+1;
-             fInfcl[0][cellcount] = icl;
-             fInfcl[1][cellcount] = jd1;
-             fInfcl[2][cellcount] = jd2;
+
+  icl       = -1;
+  cellcount = -1;
+
+  for(icell = 0; icell <= nmx1; icell++)
+    {
+      j  = iord1[icell];
+      id2 = j/kNDIMX;
+      id1 = j-id2*kNDIMX;
+
+      if(fInfocl[0][id1][id2] == 0 )
+       {
+         icl++;
+         numcell = 0;
+         cellcount++; 
+         fInfocl[0][id1][id2] = 1;
+         fInfocl[1][id1][id2] = icl;
+         fInfcl[0][cellcount] = icl;
+         fInfcl[1][cellcount] = id1;
+         fInfcl[2][cellcount] = id2;
+
+         clust[0][numcell] = id1;
+         clust[1][numcell] = id2;
+         
+         for(i = 1; i < kndim; i++)
+           {
+             clust[0][i]=0;
+           }
+         // ---------------------------------------------------------------
+         // check for adc count in neib. cells. If ne 0 put it in this clust
+         // ---------------------------------------------------------------
+         for(i = 0; i < 6; i++)
+           {
+             jd1 = id1 + neibx[i];
+             jd2 = id2 + neiby[i];
+             if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
+                 fInfocl[0][jd1][jd2] == 0)
+               {
+                 numcell++;
+                 fInfocl[0][jd1][jd2] = 2;
+                 fInfocl[1][jd1][jd2] = icl;
+                 clust[0][numcell]    = jd1;
+                 clust[1][numcell]    = jd2;
+                 cellcount++;
+                 fInfcl[0][cellcount] = icl;
+                 fInfcl[1][cellcount] = jd1;
+                 fInfcl[2][cellcount] = jd2;
+               }
+           }
+         // ---------------------------------------------------------------
+         // check adc count for neighbour's neighbours recursively and
+         // if nonzero, add these to the cluster.
+         // ---------------------------------------------------------------
+         for(i = 1; i < kndim;i++)
+           {
+             if(clust[0][i] != 0)
+               {
+                 id1 = clust[0][i];
+                 id2 = clust[1][i];
+                 for(j = 0; j < 6 ; j++)
+                   {
+                     jd1 = id1 + neibx[j];
+                     jd2 = id2 + neiby[j];
+                     if( (jd1 >= 0 && jd1 < kNDIMX) && 
+                         (jd2 >= 0 && jd2 < kNDIMY) &&
+                         fInfocl[0][jd1][jd2] == 0 )
+                       {
+                         fInfocl[0][jd1][jd2] = 2;
+                         fInfocl[1][jd1][jd2] = icl;
+                         numcell++;
+                         clust[0][numcell]    = jd1;
+                         clust[1][numcell]    = jd2;
+                         cellcount++;
+                         fInfcl[0][cellcount] = icl;
+                         fInfcl[1][cellcount] = jd1;
+                         fInfcl[2][cellcount] = jd2;
+                       }
+                   }
+               }
            }
-         }
        }
-      }
     }
-  }
-
-  //  for(icell=0; icell<=cellcount; icell++){
-  //    ofl0 << fInfcl[0][icell] << " " << fInfcl[1][icell] << " " <<
-  //      fInfcl[2][icell] << endl;
-  //}
-
   return cellcount;
 }
 // ------------------------------------------------------------------------ //
-void AliPMDClusteringV1::RefClust(int incr)
+void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
 {
   // Does the refining of clusters
   // Takes the big patch and does gaussian fitting and
   // finds out the more refined clusters
   //
-  int i, j, k, i1, i2, id, icl, ncl[4500], iord[4500], itest;
-  int ihld;
-  int ig, nsupcl, lev1[20], lev2[20];
-  double x[4500], y[4500], z[4500], x1, y1, z1, x2, y2, z2, dist;
-  double xc[4500], yc[4500], zc[4500], cells[4500], sum, rc[4500], rr;
   
+  AliPMDcludata *pmdcludata = 0;
   
-  //asso
-  Int_t t[4500],cellCount[4500];
-  for(i=0; i<4500; i++)
-    {
-      t[i]=-1;
-      cellCount[i]=0;
-    }
+  const Int_t kNmaxCell   = 19;    // # of cells surrounding a cluster center
+  
+  Int_t ndim = incr + 1;
   
+  Int_t    *ncl  = 0x0;
+  Int_t    *clxy = 0x0;  
+  Int_t    i12 = 0, i22 = 0;
+  Int_t    i = 0, j = 0, k = 0;
+  Int_t    i1 = 0, i2 = 0, id = 0, icl = 0;
+  Int_t    itest = 0, ihld = 0, ig = 0;
+  Int_t    nsupcl = 0, clno = 0, t1 = 0, t2 = 0;
+  Float_t  clusdata[6];
+  Double_t x1 = 0, y1 = 0, z1 = 0, x2 = 0, y2 = 0, z2 = 0, rr = 0;
   
-  // fClno counts the final clusters
+  ncl   = new Int_t [ndim];
+  clxy  = new Int_t [kNmaxCell];
+  
+  // Initialisation  
+  for(i = 0; i<ndim; i++)
+    {
+      ncl[i]  = -1; 
+    }
+  for(i = 0; i<6; i++)
+    {
+      clusdata[i] = 0.;
+    }
+  for(i = 0; i<19; i++)
+    {
+      clxy[i] = 0;
+    }
+
+  // clno counts the final clusters
   // nsupcl =  # of superclusters; ncl[i]= # of cells in supercluster i
   // x, y and z store (x,y) coordinates of and energy deposited in a cell
   // xc, yc store (x,y) coordinates of the cluster center
   // zc stores the energy deposited in a cluster
   // rc is cluster radius
-  // finally the cluster information is put in 2-dimensional array clusters
-  //ofstream ofl1("checking.5",ios::app);
-  fClno  = -1;
+  
+  clno  = -1;
   nsupcl = -1;
-  for(i=0; i<4500; i++){ncl[i]=-1;}
-  for(i=0; i<= incr; i++){
-    if(fInfcl[0][i] != nsupcl){ nsupcl=nsupcl+1; }
-    if (nsupcl > 4500) {
-      AliWarning("RefClust: Too many superclusters!");
-      nsupcl = 4500;
-      break;
-    }
-     
-    ncl[nsupcl]=ncl[nsupcl]+1;
-  }
 
+  for(i = 0; i <= incr; i++)
+    {
+      if(fInfcl[0][i] != nsupcl)
+       {
+         nsupcl++;
+       }
+      if (nsupcl >= ndim) 
+       {
+         AliWarning("RefClust: Too many superclusters!");
+         nsupcl = ndim - 1;
+         break;
+       }
+      ncl[nsupcl]++;
+    }
+  
   AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
                  incr+1,nsupcl+1));
-  id=-1;
-  icl=-1;
-  for(i=0; i<=nsupcl; i++) {
-    if(ncl[i] == 0){
-      id=id+1;
-      icl=icl+1;
-      // one  cell super-clusters --> single cluster
-      // cluster center at the centyer of the cell
-      // cluster radius = half cell dimension
-      if (fClno >= 5000) {
-       AliWarning("RefClust: Too many clusters! more than 5000");
-       return;
-      }
-      fClno = fClno + 1;
-      i1 = fInfcl[1][id];
-      i2 = fInfcl[2][id];
-      fClusters[0][fClno] = fCoord[0][i1][i2];
-      fClusters[1][fClno] = fCoord[1][i1][i2];
-      fClusters[2][fClno] = fEdepCell[i1][i2];
-      fClusters[3][fClno] = 1.;
-      fClusters[4][fClno] = 0.5;
-
-      //association
-
-      fClTr[0][fClno]=fCellTrNo[i1][i2];
-      for(Int_t icltr=1;icltr<14;icltr++)
+  id  = -1;
+  icl = -1;
+  
+  for(i = 0; i <= nsupcl; i++) 
+    {
+      if(ncl[i] == 0)
        {
-         fClTr[icltr][fClno]=-1;
-       }
-      
-      //ofl1 << icl << " " << fCoord[0][i1][i2] << " " << fCoord[1][i1][i2] <<
-      //" " << fEdepCell[i1][i2] << " " << fClusters[3][fClno] <<endl;
-      
-    }else if(ncl[i] == 1){
-      // two cell super-cluster --> single cluster
-      // cluster center is at ener. dep.-weighted mean of two cells
-      // cluster radius == half cell dimension
-      id   = id + 1;
-      icl  = icl+1;
-      if (fClno >= 5000) {
-       AliWarning("RefClust: Too many clusters! more than 5000");
-       return;
-      }
-      fClno = fClno+1;
-      i1   = fInfcl[1][id];
-      i2   = fInfcl[2][id];
-      x1   = fCoord[0][i1][i2];
-      y1   = fCoord[1][i1][i2];
-      z1   = fEdepCell[i1][i2];
-
-      //asso
-      fClTr[0][fClno]=fCellTrNo[i1][i2];
-      //
-
-      id   = id+1;
-      i1   = fInfcl[1][id];
-      i2   = fInfcl[2][id];
-      x2   = fCoord[0][i1][i2];
-      y2   = fCoord[1][i1][i2];
-      z2   = fEdepCell[i1][i2];
+         id++;
+         icl++;
+         if (clno >= 4608) 
+           {
+             AliWarning("RefClust: Too many clusters! more than 4608");
+             delete [] ncl;
+             delete [] clxy;
+             return;
+           }
+         clno++;
+         i1 = fInfcl[1][id];
+         i2 = fInfcl[2][id];
+         
+         i12 = i1 + i2*kNDIMX;
+         
+         clusdata[0] = fCoord[0][i1][i2];
+         clusdata[1] = fCoord[1][i1][i2];
+         clusdata[2] = edepcell[i12];
+         clusdata[3] = 1.;
+         clusdata[4] = 999.5;
+         clusdata[5] = 999.5;
 
-      //asso
+         clxy[0] = i1*10000 + i2;
+         
+         for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
+           {
+             clxy[icltr] = -1;
+           }
 
-      fClTr[1][fClno]=fCellTrNo[i1][i2];
-      for(Int_t icltr=2;icltr<14;icltr++)
-       {
-         fClTr[icltr][fClno] = -1;
+         pmdcludata  = new AliPMDcludata(clusdata,clxy);
+         fPMDclucont->Add(pmdcludata);
        }
-      //
+      else if(ncl[i] == 1) 
+       {
+         id++;
+         icl++;
+         if (clno >= 4608) 
+           {
+             AliWarning("RefClust: Too many clusters! more than 4608");
+             delete [] ncl;
+             delete [] clxy;
 
-      fClusters[0][fClno] = (x1*z1+x2*z2)/(z1+z2);
-      fClusters[1][fClno] = (y1*z1+y2*z2)/(z1+z2);
-      fClusters[2][fClno] = z1+z2;
-      fClusters[3][fClno] = 2.;
-      fClusters[4][fClno] = 0.5;
+             return;
+           }
+         clno++;
+         i1   = fInfcl[1][id];
+         i2   = fInfcl[2][id];
+         i12  = i1 + i2*kNDIMX;
 
+         x1   = fCoord[0][i1][i2];
+         y1   = fCoord[1][i1][i2];
+         z1   = edepcell[i12];
 
-      //ofl1 << icl << " " << fClusters[0][fClno] << " " << fClusters[1][fClno]
-      //  << " " << fClusters[2][fClno] << " " <<fClusters[3][fClno] <<endl;
-    }
-    else{
-      
-      //asso
-      for(Int_t icg=0;icg<4500;icg++)
-       {
-         cellCount[icg]=0;
-       }
-      //
+         clxy[0] = i1*10000 + i2;
+         
+         id++;
+         i1   = fInfcl[1][id];
+         i2   = fInfcl[2][id];
 
-      id      = id + 1;
-      iord[0] = 0;
-      // super-cluster of more than two cells - broken up into smaller
-      // clusters gaussian centers computed. (peaks separated by > 1 cell)
-      // Begin from cell having largest energy deposited This is first
-      // cluster center
-      i1      = fInfcl[1][id];
-      i2      = fInfcl[2][id];
-      x[0]    = fCoord[0][i1][i2];
-      y[0]    = fCoord[1][i1][i2];
-      z[0]    = fEdepCell[i1][i2];
-      
-      //asso
-      t[0]=fCellTrNo[i1][i2];
-      //
+         i22  = i1 + i2*kNDIMX;
+         x2   = fCoord[0][i1][i2];
+         y2   = fCoord[1][i1][i2];
+         z2   = edepcell[i22];
+
+         clxy[1] = i1*10000 + i2;
+         
 
-      iord[0] = 0;
-      for(j=1;j<=ncl[i];j++){
-
-       id      = id + 1;
-       i1      = fInfcl[1][id];
-       i2      = fInfcl[2][id];
-       iord[j] = j;
-       x[j]    = fCoord[0][i1][i2];
-       y[j]    = fCoord[1][i1][i2];
-       z[j]    = fEdepCell[i1][i2];
-
-       //asso
-       t[j]=fCellTrNo[i1][i2];
-       //
-
-
-      }
-      // arranging cells within supercluster in decreasing order
-      for(j=1;j<=ncl[i];j++){
-       itest=0;
-       ihld=iord[j];
-       for(i1=0;i1<j;i1++){
-         if(itest == 0 && z[iord[i1]] < z[ihld]){
-           itest=1;
-           for(i2=j-1;i2>=i1;i2--){
-             iord[i2+1]=iord[i2];
+         for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
+           {
+             clxy[icltr] = -1;
            }
-           iord[i1]=ihld;
-         }
+         
+         clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
+         clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
+         clusdata[2] = z1+z2;
+         clusdata[3] = 2.;
+         clusdata[4] = 0.5;
+         clusdata[5] = 0.0;
+         pmdcludata  = new AliPMDcludata(clusdata,clxy);
+         fPMDclucont->Add(pmdcludata);
        }
-      }
-
-      // compute the number of Gaussians and their centers ( first
-      // guess )
-      // centers must be separated by cells having smaller ener. dep.
-      // neighbouring centers should be either strong or well-separated
-      ig=0;
-      xc[ig]=x[iord[0]];
-      yc[ig]=y[iord[0]];
-      zc[ig]=z[iord[0]];
-      for(j=1;j<=ncl[i];j++){
-       itest=-1;
-       x1=x[iord[j]];
-       y1=y[iord[j]];
-       for(k=0;k<=ig;k++){
-         x2=xc[k]; y2=yc[k];
-         rr=Distance(x1,y1,x2,y2);
-         if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)
-           itest=itest+1;
-         if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)
-           itest=itest+1;
-         if( rr >= 2.1)itest=itest+1;
-       }
-       if(itest == ig){
-         ig=ig+1;
-         xc[ig]=x1;
-         yc[ig]=y1;
-         zc[ig]=z[iord[j]];
-       }
-      }
-      // for(j=0; j<=ig; j++){
-      //ofl1 << icl+j+1 << " " << xc[j] << " " <<yc[j] <<" "<<zc[j]<<endl;
-      //}
-      // GaussFit to adjust cluster parameters to minimize
-      GaussFit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]);
-      icl=icl+ig+1;
-      // compute the number of cells belonging to each cluster.
-      // cell is shared between several clusters ( if they are equidistant
-      // from it ) in the ratio of cluster energy deposition
-      for(j=0; j<=ig; j++){
-       cells[j]=0.;
-      }
-      if(ig > 0){
-       for(j=0; j<=ncl[i]; j++){
-         lev1[j]=0;
-         lev2[j]=0;
-         for(k=0; k<=ig; k++){
-           dist=Distance(x[j], y[j], xc[k], yc[k]);
-           if(dist < sqrt(3.) ){
-
-             //asso
-             fClTr[cellCount[k]][fClno+k+1]=t[j];
-             cellCount[k]++;
-             //
+      else
+       {
+         Int_t    *iord, *tc, *t;
+         Double_t *x, *y, *z, *xc, *yc, *zc;
 
-             lev1[0]++;
-             i1=lev1[0];
-             lev1[i1]=k;
-           }else{
-             if(dist < 2.1){
-               lev2[0]++;
-               i1=lev2[0];
-               lev2[i1]=k;
-             }
+         iord = new Int_t [ncl[i]+1];
+         tc   = new Int_t [ncl[i]+1];
+         t    = new Int_t [ncl[i]+1];
+         x    = new Double_t [ncl[i]+1];
+         y    = new Double_t [ncl[i]+1];
+         z    = new Double_t [ncl[i]+1];
+         xc   = new Double_t [ncl[i]+1];
+         yc   = new Double_t [ncl[i]+1];
+         zc   = new Double_t [ncl[i]+1];
+         
+         for( k = 0; k < ncl[i]+1; k++)
+           {
+             iord[k] = -1;
+             t[k]    = -1;
+             tc[k]   = -1;
+             x[k]    = -1;
+             y[k]    = -1;
+             z[k]    = -1;
+             xc[k]   = -1;
+             yc[k]   = -1;
+             zc[k]   = -1;
            }
-         }
-         if(lev1[0] != 0){
-           if(lev1[0] == 1){cells[lev1[1]]=cells[lev1[1]]+1.;}
-           else{
-             sum=0.;
-             for(k=1; k<=lev1[0]; k++){
-               sum=sum+zc[lev1[k]];
-             }
-             for(k=1; k<=lev1[0]; k++){
-               cells[lev1[k]]=cells[lev1[k]]+zc[lev1[k]]/sum;
-             }
+         id++;
+         // super-cluster of more than two cells - broken up into smaller
+         // clusters gaussian centers computed. (peaks separated by > 1 cell)
+         // Begin from cell having largest energy deposited This is first
+         // cluster center
+         i1      = fInfcl[1][id];
+         i2      = fInfcl[2][id];
+         i12     = i1 + i2*kNDIMX;
+         
+         x[0]    = fCoord[0][i1][i2];
+         y[0]    = fCoord[1][i1][i2];
+         z[0]    = edepcell[i12];
+         t[0]    = i1*10000 + i2;
+         
+
+         iord[0] = 0;
+         for(j = 1; j <= ncl[i]; j++)
+           {
+             id++;
+             i1      = fInfcl[1][id];
+             i2      = fInfcl[2][id];
+             i12     = i1 + i2*kNDIMX;
+             iord[j] = j;
+             x[j]    = fCoord[0][i1][i2];
+             y[j]    = fCoord[1][i1][i2];
+             z[j]    = edepcell[i12];
+             t[j]    = i1*10000 + i2;
            }
-         }else{
-           if(lev2[0] == 0){cells[lev2[1]]=cells[lev2[1]]+1.;}
-           else{
-             sum=0.;
-             for(k=1; k<=lev2[0]; k++){
-               sum=sum+zc[lev2[k]];
-             }
-             for(k=1; k<=lev2[0]; k++){
-               cells[lev2[k]]=cells[lev2[k]]+zc[lev2[k]]/sum;
-             }
+         
+         // arranging cells within supercluster in decreasing order
+         
+         for(j = 1;j <= ncl[i]; j++)
+           {
+             itest = 0;
+             ihld  = iord[j];
+             for(i1 = 0; i1 < j; i1++)
+               {
+                 if(itest == 0 && z[iord[i1]] < z[ihld])
+                   {
+                     itest = 1;
+                     for(i2 = j-1; i2 >= i1; i2--)
+                       {
+                         iord[i2+1] = iord[i2];
+                       }
+                     iord[i1] = ihld;
+                   }
+               }
            }
-         }
-       }
-      }
 
-      // zero rest of the cell array
-      //asso
-      for(k=0; k<=ig; k++)
-       {
-         for(Int_t icltr=cellCount[k];icltr<14;icltr++)
+         if (fClusParam == 1)
+           {
+             // Clustering algorithm returns from here for PP collisions
+             // for pp, only the output of crude clusterng is taken
+             // sigx and sigy are not calculated at this moment
+
+             Double_t supx=0., supy=0., supz=0.;
+         
+             for(j = 0;j <= ncl[i]; j++)
+               {
+                 supx += x[iord[j]]*z[iord[j]];
+                 supy += y[iord[j]]*z[iord[j]];
+                 supz += z[iord[j]];
+                 if(j < 19)
+                   {
+                     clxy[j] = t[iord[j]];
+                   }
+               }
+             
+             if( ncl[i] + 1 < 19)
+               {
+                 for(Int_t ncel = ncl[i] + 1; ncel < kNmaxCell; ncel ++ )
+                   {
+                     clxy[ncel] = -1;
+                   }
+               }
+             clusdata[0] = supx/supz;
+             clusdata[1] = supy/supz;
+             clusdata[2] = supz;
+             clusdata[3] = ncl[i]+1;
+             clusdata[4] = 0.5;
+             clusdata[5] = 0.0;
+             pmdcludata  = new AliPMDcludata(clusdata,clxy);
+             fPMDclucont->Add(pmdcludata);
+           }
+
+         /* MODIFICATION PART STARTS (Tapan July 2008)
+            iord[0] is the cell with highest ADC in the crude-cluster
+            ig is the number of local maxima in the crude-cluster
+            For the higest peak we make ig=0 which means first local maximum.
+            Next we go down in terms of the ADC sequence and find out if any
+            more of the cells form local maxima. The definition of local 
+            maxima is that all its neighbours are of less ADC compared to it. 
+         */
+
+         if (fClusParam == 2)
            {
-             fClTr[icltr][fClno]=-1;
+             // This part is to split the supercluster
+             //
+             ig = 0;
+             xc[ig] = x[iord[0]];
+             yc[ig] = y[iord[0]];
+             zc[ig] = z[iord[0]];
+             tc[ig] = t[iord[0]];
+             Int_t ivalid = 0,  icount = 0;
+         
+             for(j=1;j<=ncl[i];j++)
+               {
+                 x1  = x[iord[j]];
+                 y1  = y[iord[j]]; 
+                 z1  = z[iord[j]];
+                 t1  = t[iord[j]];
+                 rr=Distance(x1,y1,xc[ig],yc[ig]);
+                 
+                 // Check the cells which are outside the neighbours (rr>1.2)
+                 if(rr>1.2 ) 
+                   {
+                     ivalid=0;
+                     icount=0;
+                     for(Int_t j1=1;j1<j;j1++)
+                       {
+                         icount++;
+                         Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
+                         if(rr1>1.2) ivalid++;
+                       }
+                     if(ivalid == icount && z1>0.5*zc[ig])
+                       {
+                         ig++;
+                         xc[ig]=x1; 
+                         yc[ig]=y1; 
+                         zc[ig]=z1;
+                         tc[ig]=t1;
+                       }
+                   }     
+               }
+         
+             icl=icl+ig+1;
+         
+             //  We use simple Gaussian weighting. (Tapan Jan 2005)
+             // compute the number of cells belonging to each cluster.
+             // cell can be shared between several clusters  
+             // in the ratio of cluster energy deposition
+             // To calculate: 
+             // (1) number of cells belonging to a cluster (ig) and 
+             // (2) total ADC of the cluster (ig) 
+             // (3) x and y positions of the cluster
+         
+             
+             Int_t *cellCount;
+             Int_t **cellXY;
+             
+             Int_t    *status;
+             Double_t *totaladc, *totaladc2, *ncell,*weight;
+             Double_t *xclust, *yclust, *sigxclust, *sigyclust;
+             Double_t *ax, *ay, *ax2, *ay2;
+         
+         
+             status    = new Int_t [ncl[i]+1];
+             cellXY    = new Int_t *[ncl[i]+1];
+         
+             cellCount = new Int_t [ig+1];
+             totaladc  = new Double_t [ig+1];
+             totaladc2 = new Double_t [ig+1];
+             ncell     = new Double_t [ig+1];
+             weight    = new Double_t [ig+1];
+             xclust    = new Double_t [ig+1];
+             yclust    = new Double_t [ig+1];
+             sigxclust = new Double_t [ig+1];
+             sigyclust = new Double_t [ig+1];
+             ax        = new Double_t [ig+1];
+             ay        = new Double_t [ig+1];
+             ax2       = new Double_t [ig+1];
+             ay2       = new Double_t [ig+1];
+             
+             for(j = 0; j < ncl[i]+1; j++) 
+               {
+                 status[j]     = 0;
+                 cellXY[j] = new Int_t[ig+1];
+               }
+             //initialization
+             for(Int_t kcl = 0; kcl < ig+1; kcl++)
+               {
+                 cellCount[kcl] = 0;
+                 totaladc[kcl]  = 0.;
+                 totaladc2[kcl] = 0.;
+                 ncell[kcl]     = 0.;
+                 weight[kcl]    = 0.;  
+                 xclust[kcl]    = 0.; 
+                 yclust[kcl]    = 0.;
+                 sigxclust[kcl] = 0.; 
+                 sigyclust[kcl] = 0.;
+                 ax[kcl]        = 0.;      
+                 ay[kcl]        = 0.;      
+                 ax2[kcl]       = 0.;      
+                 ay2[kcl]       = 0.;    
+                 for(j = 0; j < ncl[i]+1; j++)
+                   {
+                     cellXY[j][kcl] = 0;
+                   }
+               }
+             Double_t sumweight, gweight; 
+             
+             for(j = 0;j <= ncl[i]; j++)
+               {
+                 x1 = x[iord[j]];
+                 y1 = y[iord[j]];
+                 z1 = z[iord[j]];
+                 t1 = t[iord[j]];
+                 
+                 for(Int_t kcl=0; kcl<=ig; kcl++)
+                   {
+                     x2 = xc[kcl];
+                     y2 = yc[kcl];
+                     rr = Distance(x1,y1,x2,y2);
+                     t2 = tc[kcl];               
+                     
+                     if(rr==0)
+                       {
+                         ncell[kcl]     = 1.;
+                         totaladc[kcl]  = z1;
+                         totaladc2[kcl] = z1*z1;
+                         ax[kcl]        = x1 * z1;
+                         ay[kcl]        = y1 * z1;
+                         ax2[kcl]       = 0.;
+                         ay2[kcl]       = 0.;
+                         status[j]      = 1;
+                       }
+                   }
+               }
+         
+             for(j = 0; j <= ncl[i]; j++)
+               {
+                 Int_t   maxweight = 0;
+                 Double_t     max  = 0.;
+                 
+                 if(status[j] == 0)
+                   { 
+                     x1 = x[iord[j]]; 
+                     y1 = y[iord[j]];
+                     z1 = z[iord[j]];
+                     t1 = t[iord[j]];
+                     sumweight = 0.;
+
+                     for(Int_t kcl = 0; kcl <= ig; kcl++)
+                       {
+                         x2 = xc[kcl]; 
+                         y2 = yc[kcl]; 
+                         rr = Distance(x1,y1,x2,y2);
+                         gweight     = exp(-(rr*rr)/(2*(1.2*1.2)));
+                         weight[kcl] = zc[kcl] * gweight;
+                         sumweight   = sumweight + weight[kcl];
+                         
+                         if(weight[kcl] > max)
+                           {
+                             max       =  weight[kcl];
+                             maxweight =  kcl;
+                           }
+                       }
+                     
+                     cellXY[cellCount[maxweight]][maxweight] = iord[j];
+                     
+                     cellCount[maxweight]++;
+                     
+                     x2 = xc[maxweight];
+                     y2 = yc[maxweight];
+                     totaladc[maxweight]  +=  z1;
+                     ax[maxweight]        +=  x1*z1;
+                     ay[maxweight]        +=  y1*z1;
+                     totaladc2[maxweight] +=  z1*z1;
+                     ax2[maxweight]       +=  z1*(x1-x2)*(x1-x2);
+                     ay2[maxweight]       +=  z1*(y1-y2)*(y1-y2);
+                     ncell[maxweight]++;
+                     
+                   }
+               }
+         
+             for(Int_t kcl = 0; kcl <= ig; kcl++)
+               {
+                 if(totaladc[kcl] > 0.)
+                   {
+                     xclust[kcl] = (ax[kcl])/ totaladc[kcl];
+                     yclust[kcl] = (ay[kcl])/ totaladc[kcl];
+                     
+                     //natasha
+                     Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
+                     if(totaladc2[kcl] >= sqtotadc)
+                       {
+                         sigxclust[kcl] = 0.25;
+                         sigyclust[kcl] = 0.25;
+                       }
+                     else
+                       {
+                         sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
+                         sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
+                       }       
+                   }
+                 
+                 for(j = 0; j < cellCount[kcl]; j++) clno++; 
+                 
+                 if (clno >= 4608) 
+                   {
+                     AliWarning("RefClust: Too many clusters! more than 4608");
+
+                     delete [] cellCount;
+                     for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
+                     delete [] cellXY;
+
+                     delete [] status;
+                     delete [] totaladc;
+                     delete [] totaladc2;
+                     delete [] ncell;
+                     delete [] xclust;
+                     delete [] yclust;
+                     delete [] sigxclust;
+                     delete [] sigyclust;
+                     delete [] ax;
+                     delete [] ay;
+                     delete [] ax2;
+                     delete [] ay2;
+                     delete [] weight;
+
+                     delete [] iord;
+                     delete [] tc;       
+                     delete [] t;
+                     delete [] x;
+                     delete [] y;
+                     delete [] z;
+                     delete [] xc;
+                     delete [] yc;
+                     delete [] zc;
+
+
+                     delete [] ncl;
+                     delete [] clxy;
+
+                     return;
+                   }
+                 clusdata[0] = xclust[kcl];
+                 clusdata[1] = yclust[kcl];
+                 clusdata[2] = totaladc[kcl];
+                 clusdata[3] = ncell[kcl];
+                 
+                 if(sigxclust[kcl] > sigyclust[kcl]) 
+                   {
+                     clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
+                     clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
+                   }
+                 else
+                   {
+                     clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
+                     clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
+                   }
+                 
+                 clxy[0] = tc[kcl];
+                 
+                 Int_t Ncell=1;
+                 for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
+                   {
+                     if(ii<18) 
+                       {       
+                         clxy[Ncell] = t[cellXY[ii][kcl]];
+                         Ncell++;
+                       }
+                   } 
+                 
+                 pmdcludata = new AliPMDcludata(clusdata,clxy);
+                 fPMDclucont->Add(pmdcludata);
+               }
+             delete [] cellCount;
+             for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
+             delete [] cellXY;
+
+             delete [] status;
+             delete [] totaladc;
+             delete [] totaladc2;
+             delete [] ncell;
+             delete [] xclust;
+             delete [] yclust;
+             delete [] sigxclust;
+             delete [] sigyclust;
+             delete [] ax;
+             delete [] ay;
+             delete [] ax2;
+             delete [] ay2;
+             delete [] weight;
+
            }
-       }
-      //
 
+         delete [] iord;
+         delete [] tc;   
+         delete [] t;
+         delete [] x;
+         delete [] y;
+         delete [] z;
+         delete [] xc;
+         delete [] yc;
+         delete [] zc;
 
 
-      for(j=0; j<=ig; j++){
-       if (fClno >= 5000) {
-         AliWarning("RefClust: Too many clusters! more than 5000");
-         return;
-       }
-       fClno               = fClno + 1;
-       fClusters[0][fClno] = xc[j];
-       fClusters[1][fClno] = yc[j];
-       fClusters[2][fClno] = zc[j];
-       fClusters[4][fClno] = rc[j];
-       if(ig == 0){
-         fClusters[3][fClno] = ncl[i];
-       }else{
-         fClusters[3][fClno] = cells[j];
        }
-      }
-    }
-  }
-}
-// ------------------------------------------------------------------------ //
-void AliPMDClusteringV1::GaussFit(Int_t ncell, Int_t nclust, Double_t &x, Double_t &y ,Double_t &z, Double_t &xc, Double_t &yc, Double_t &zc, Double_t &rc)
-{
-  // Does gaussian fitting
-  //
-  int i, j, i1, i2, novar, idd, jj;
-  double xx[4500], yy[4500], zz[4500], xxc[4500], yyc[4500];
-  double a[4500], b[4500], c[4500], d[4500], ha[4500], hb[4500];
-  double hc[4500], hd[4500], zzc[4500], rrc[4500];
-  int neib[4500][50];
-  double sum, dx, dy, str, str1, aint, sum1, rr, dum;
-  double x1, x2, y1, y2;
-  str   = 0.;
-  str1  = 0.;
-  rr    = 0.3;
-  novar = 0;
-  j = 0;  // Just put not to see the compiler warning, BKN
-
-  for(i=0; i<=ncell; i++)
-    {
-      xx[i] = *(&x+i);
-      yy[i] = *(&y+i);
-      zz[i] = *(&z+i);
-      str   = str + zz[i];
-    }
-  for(i=0; i<=nclust; i++)
-    {
-      xxc[i] = *(&xc+i);
-      yyc[i] = *(&yc+i);
-      zzc[i] = *(&zc+i);
-      str1   = str1 + zzc[i];
-      rrc[i] = 0.5;
-    }
-  for(i=0; i<=nclust; i++)
-    {
-      zzc[i] = str/str1*zzc[i];
-      ha[i]  = xxc[i];
-      hb[i]  = yyc[i];
-      hc[i]  = zzc[i];
-      hd[i]  = rrc[i];
-      x1     = xxc[i];
-      y1     = yyc[i];
-    }
-  for(i=0; i<=ncell; i++){
-    idd=0;
-    x1=xx[i];
-    y1=yy[i];
-    for(j=0; j<=nclust; j++){
-      x2=xxc[j];
-      y2=yyc[j];
-      if(Distance(x1,y1,x2,y2) <= 3.){ idd=idd+1; neib[i][idd]=j; }
     }
-    neib[i][0]=idd;
-  }
-  sum=0.;
-  for(i1=0; i1<=ncell; i1++){
-    aint=0.;
-    idd=neib[i1][0];
-    for(i2=1; i2<=idd; i2++){
-      jj=neib[i1][i2];
-      dx=xx[i1]-xxc[jj];
-      dy=yy[i1]-yyc[jj];
-      dum=rrc[j]*rrc[jj]+rr*rr;
-      aint=aint+exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum;
-    }
-    sum=sum+(aint-zz[i1])*(aint-zz[i1])/str;
-  }
-//   jmax=nclust*1000;
-//   if(nclust > 20)jmax=20000;
-//   for(j=0; j<jmax; j++){
-    str1=0.;
-    for(i=0; i<=nclust; i++){
-      a[i]=xxc[i]+0.6*(Ranmar()-0.5);
-      b[i]=yyc[i]+0.6*(Ranmar()-0.5);
-      c[i]=zzc[i]*(1.+(Ranmar()-0.5)*0.2);
-      str1=str1+zzc[i];
-      d[i]=rrc[i]*(1.+(Ranmar()-0.5)*0.1);
-      if(d[i] < 0.25)d[i]=0.25;
-    }
-    for(i=0; i<=nclust; i++){ c[i]=c[i]*str/str1; }
-    sum1=0.;
-    for(i1=0; i1<=ncell; i1++){
-      aint=0.;
-      idd=neib[i1][0];
-      for(i2=1; i2<=idd; i2++){
-       jj=neib[i1][i2];
-       dx=xx[i1]-a[jj];
-       dy=yy[i1]-b[jj];
-       dum=d[jj]*d[jj]+rr*rr;
-       aint=aint+exp(-(dx*dx+dy*dy)/dum)*c[i2]*rr*rr/dum;
-      }
-      sum1=sum1+(aint-zz[i1])*(aint-zz[i1])/str;
-    }
-
-    if(sum1 < sum){
-      for(i2=0; i2<=nclust; i2++){
-       xxc[i2]=a[i2];
-       yyc[i2]=b[i2];
-       zzc[i2]=c[i2];
-       rrc[i2]=d[i2];
-       sum=sum1;
-      }
-    }
-//   }
-  for(j=0; j<=nclust; j++){
-    *(&xc+j)=xxc[j];
-    *(&yc+j)=yyc[j];
-    *(&zc+j)=zzc[j];
-    *(&rc+j)=rrc[j];
-  }
+  delete [] ncl;
+  delete [] clxy;
 }
 // ------------------------------------------------------------------------ //
-double AliPMDClusteringV1::Distance(double x1, double y1, double x2, double y2)
+Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1, 
+                                     Double_t x2, Double_t y2)
 {
-  return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
+  return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
 }
 // ------------------------------------------------------------------------ //
-double AliPMDClusteringV1::Ranmar() const
+void AliPMDClusteringV1::SetEdepCut(Float_t decut)
 {
-  //  Universal random number generator proposed by Marsaglia and Zaman
-  //  in report FSU-SCRI-87-50
-
-  //  clock_t start;
-  int ii, jj;
-  static int i=96, j=32, itest=0, i1, i2, i3, i4, i5;
-  static double u[97], c, cd, cm, s, t;
-  static double uni;
-  int count1,count2,idum;
-  /*    $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$  */
-  if (itest == 0) {
-    //*******************************************************
-    // following three lines if the seed to be provided by computer
-    // start = time(NULL);
-    // ii=start;
-    // jj=start;
-    //*******************************************************
-    //following two lines for fixed seed ( during testing only. Else
-    //use preceeing three lines
-    ii=8263;
-    jj=5726;
-    if(ii > 31328 ) ii = ii - ( ii / 31328 ) * 31328;
-    if(jj > 30081 ) jj = jj - ( jj / 30081 ) * 30081;
-    itest=itest+1;
-    if((( ii > 0 ) &&  ( ii <= 31328 )) && (( jj > 0 ) &&
-                                           ( jj <= 30081 ))){
-      i1=ii/177+2; i2=ii-(i1-2)*177+2; i3=jj/169+1; i4=jj-(i3-1)*169;
-      i4 = jj - (i3-1)*169;
-      count1=0;
-      while ( count1 < 97 ){
-       s=0.;
-       t=0.5;
-       count2=0;
-       while( count2 < 24 ){
-         idum=i1*i2/179;
-         idum=( i1*i2 - (i1*i2/179)*179 ) * i3;
-         i5=idum-(idum/179)*179;
-         i1=i2; i2=i3; i3=i5; idum=53*i4+1; i4=idum-(idum/169)*169;
-         if( i4*i5-((i4*i5)/64)*64 >= 32 ) s=s+t;
-         t=0.5*t;
-         count2=count2+1;
-       }
-       u[count1] = s;
-       count1 = count1 +1;
-      }
-      c = 362436./16777216.;  cd = 7654321./16777216.;
-      cm = 16777213./16777216.;
-    }
-    else{
-      AliWarning("Wrong initialization");
-    }
-  }
-  else{
-    uni = u[i] - u[j];
-    if( uni < 0.) uni = uni + 1;
-    u[i] = uni;
-    i = i -1;
-    if( i < 0 ) i = 96;
-    j = j - 1;
-    if ( j < 0 ) j = 96;
-    c = c - cd;
-    if( c < 0. ) c = c+cm;
-    uni = uni-c ;
-    if( uni < 0. )uni = uni+1.;
-  }
-  return uni;
+  fCutoff = decut;
 }
 // ------------------------------------------------------------------------ //
-void AliPMDClusteringV1::SetEdepCut(Float_t decut)
+void AliPMDClusteringV1::SetClusteringParam(Int_t cluspar)
 {
-  fCutoff = decut;
+  fClusParam = cluspar;
 }
 // ------------------------------------------------------------------------ //