new version clustering algo
authorbnandi <bnandi@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 13 Jul 2006 12:07:51 +0000 (12:07 +0000)
committerbnandi <bnandi@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 13 Jul 2006 12:07:51 +0000 (12:07 +0000)
PMD/AliPMDClusteringV2.cxx [new file with mode: 0644]
PMD/AliPMDClusteringV2.h [new file with mode: 0644]

diff --git a/PMD/AliPMDClusteringV2.cxx b/PMD/AliPMDClusteringV2.cxx
new file mode 100644 (file)
index 0000000..d8f446d
--- /dev/null
@@ -0,0 +1,835 @@
+/***************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------//
+//                                                     //
+//  Source File : PMDClusteringV2.cxx                  //
+//                                                     //
+//  clustering code for alice pmd                      //
+//                                                     //
+//-----------------------------------------------------//
+
+/* --------------------------------------------------------------------
+   Code developed by S. C. Phatak, Institute of Physics,
+   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
+   supermodule.
+
+   fEdepCell, fClno  and fClusters are 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 <TObjArray.h>
+#include <stdio.h>
+
+#include "AliPMDcluster.h"
+#include "AliPMDClustering.h"
+#include "AliPMDClusteringV2.h"
+#include "AliLog.h"
+
+ClassImp(AliPMDClusteringV2)
+
+const Double_t AliPMDClusteringV2::fgkSqroot3by2=0.8660254;  // sqrt(3.)/2.
+
+AliPMDClusteringV2::AliPMDClusteringV2():
+  fCutoff(0.0)
+{
+  for(int i = 0; i < kNDIMX; i++)
+    {
+      for(int j = 0; j < kNDIMY; j++)
+       {
+         fCoord[0][i][j] = i+j/2.;
+         fCoord[1][i][j] = fgkSqroot3by2*j;
+         fEdepCell[i][j] = 0;
+       }
+    }
+}
+// ------------------------------------------------------------------------ //
+AliPMDClusteringV2::~AliPMDClusteringV2()
+{
+
+}
+// ------------------------------------------------------------------------ //
+void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96], TObjArray *pmdcont)
+{
+  // main function to call other necessary functions to do clustering
+  //
+  AliPMDcluster *pmdcl = 0;
+
+  Int_t    i, i1, i2, j, nmx1, incr, id, jd;
+  Int_t    celldataX[15], celldataY[15];
+  Float_t  clusdata[6];
+  Double_t cutoff, ave;
+
+  const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
+
+  Int_t ndimXr =0;
+  Int_t ndimYr =0;
+
+  if (ismn < 12)
+    {
+      ndimXr = 96;
+      ndimYr = 48;
+    }
+  else if (ismn >= 12 && ismn <= 23)
+    {
+      ndimXr = 48;
+      ndimYr = 96;
+    }
+
+  for (Int_t i =0; i < kNDIMX; i++)
+    {
+      for (Int_t j =0; j < kNDIMY; j++)
+       {
+         fEdepCell[i][j] = 0;
+       }
+    }
+
+
+  for (id = 0; id < ndimXr; id++)
+    {
+      for (jd = 0; jd < ndimYr; jd++)
+       {
+         j=jd;
+         i=id+(ndimYr/2-1)-(jd/2);
+         
+         if (ismn < 12)
+           {
+             fEdepCell[i][j] = celladc[jd][id];
+           }
+         else if (ismn >= 12 && ismn <= 23)
+           {
+             fEdepCell[i][j] = 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++)
+    {
+      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;
+    }
+  // nmx1 --- number of cells having ener dep >= cutoff
+
+  AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
+
+  if (nmx1 == 0) nmx1 = 1;
+  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++)
+    {
+      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 sigmaX   = (Float_t) fClusters[4][i1];
+      Float_t sigmaY   = (Float_t) fClusters[5][i1];
+      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)
+       {
+         clusdata[0] = cluX0 - (48-1) + cluY0/2.;
+       }         
+
+      clusdata[1]      = cluY0;
+      clusdata[2]      = cluADC;
+      clusdata[3]      = cluCELLS;
+      clusdata[4]      = sigmaX;
+      clusdata[5]      = sigmaY;
+
+      //
+      // Cells associated with a cluster
+      //
+      for (Int_t ihit = 0; ihit < 15; ihit++)
+       {
+         celldataX[ihit] = 1;  // dummy nos. -- will be changed
+         celldataY[ihit] = 1;  // dummy nos. -- will be changed
+       }
+
+      pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
+      pmdcont->Add(pmdcl);
+    }
+}
+// ------------------------------------------------------------------------ //
+void AliPMDClusteringV2::Order()
+{
+  // 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
+  //
+  // define arrays dd and iord1
+  for(i1=0; i1 < kNDIMX; i1++)
+    {
+      for(i2=0; i2 < kNDIMY; i2++)
+       {
+         i        = i1 + i2*kNDIMX;
+         iord1[i] = i;
+         dd[i]    = fEdepCell[i1][i2];
+       }
+    }
+  // sort and store sorting information in iord1
+
+  TMath::Sort(kNMX,dd,iord1);
+
+  // store the sorted information in fIord for later use
+  for(i=0; i<kNMX; i++)
+    {
+      j  = iord1[i];
+      i2 = j/kNDIMX;
+      i1 = j-i2*kNDIMX;
+      fIord[0][i]=i1;
+      fIord[1][i]=i2;
+    }
+}
+// ------------------------------------------------------------------------ //
+Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1)
+{
+  // Does crude clustering
+  // Finds out only the big patch by just searching the
+  // connected cells
+  //
+
+  int i,j,k,id1,id2,icl, numcell;
+  int jd1,jd2, icell, cellcount;
+  int clust[2][5000];
+  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]
+
+  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;
+    }
+  }
+  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] = -1;
+      // ---------------------------------------------------------------
+      // 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] != -1){
+         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;
+           }
+         }
+       }
+      }
+    }
+  }
+  //  for(icell=0; icell<=cellcount; icell++){
+  //    ofl0 << fInfcl[0][icell] << " " << fInfcl[1][icell] << " " <<
+  //      fInfcl[2][icell] << endl;
+  //  }
+  return cellcount;
+}
+// ------------------------------------------------------------------------ //
+void AliPMDClusteringV2::RefClust(Int_t incr)
+{
+  // Does the refining of clusters
+  // Takes the big patch and does gaussian fitting and
+  // finds out the more refined clusters
+  //
+
+  const Int_t kndim = 4500;
+
+  int i, j, k, i1, i2, id, icl, itest;
+  int ihld;
+  int ig, nsupcl;
+  int ncl[kndim], iord[kndim];
+
+  double x1, y1, z1, x2, y2, z2;
+  double rr;
+
+  double x[kndim], y[kndim], z[kndim];
+  double xc[kndim], yc[kndim], zc[kndim], cells[kndim];
+  double rcl[kndim], rcs[kndim];
+
+  // fClno 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;
+  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;
+  }
+
+  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++;
+      icl++;
+      // 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++;
+      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.0;
+      fClusters[5][fClno] = 0.0;
+      //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++;
+      icl++;
+      if (fClno >= 5000) {
+       AliWarning("RefClust: Too many clusters! more than 5000");
+       return;
+      }
+      fClno++;
+      i1   = fInfcl[1][id];
+      i2   = fInfcl[2][id];
+      x1   = fCoord[0][i1][i2];
+      y1   = fCoord[1][i1][i2];
+      z1   = fEdepCell[i1][i2];
+
+      id++;
+      i1   = fInfcl[1][id];
+      i2   = fInfcl[2][id];
+      x2   = fCoord[0][i1][i2];
+      y2   = fCoord[1][i1][i2];
+      z2   = fEdepCell[i1][i2];
+
+      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] = sqrt(z1*z2)/(z1+z2);
+      fClusters[5][fClno] = 0;  // sigma large nonzero, sigma small zero
+
+      //ofl1 << icl << " " << fClusters[0][fClno] << " " << fClusters[1][fClno]
+      //   << " " << fClusters[2][fClno] << " " <<fClusters[3][fClno] <<endl;
+    }
+    else{
+      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
+      // *****************************************************************
+      // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
+      // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
+      // SINCE WE EXPECT THE SUPERCLUSTER 
+      // TO BE A SINGLE CLUSTER
+      //*******************************************************************
+
+      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];
+      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];
+      }
+      // 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;
+               }
+           }
+       }
+
+      // compute the number of clusters 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);
+         //***************************************************************
+         // finetuning cluster splitting
+         // the numbers zc/4 and zc/10 may need to be changed. 
+         // Also one may need to add one more layer because our 
+         // cells are smaller in absolute scale
+         //****************************************************************
+
+
+         if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)
+           itest++;
+         if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)
+           itest++;
+         if( rr >= 2.1)itest++;
+       }
+       if(itest == ig){
+         ig++;
+         xc[ig] = x1;
+         yc[ig] = y1;
+         zc[ig] = z[iord[j]];
+       }
+      }
+
+      ClustDetails(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0],
+                  rcl[0], rcs[0], cells[0]);
+
+      icl = icl + ig + 1;
+
+      for(j=0; j<=ig; j++)
+       {
+         if (fClno >= 5000)
+           {
+             AliWarning("RefClust: Too many clusters! more than 5000");
+             return;
+           }
+         fClno++;
+         fClusters[0][fClno] = xc[j];
+         fClusters[1][fClno] = yc[j];
+         fClusters[2][fClno] = zc[j];
+         fClusters[4][fClno] = rcl[j];
+         fClusters[5][fClno] = rcs[j];
+         if(ig == 0)
+           {
+             fClusters[3][fClno] = ncl[i];
+           }
+         else
+           {
+             fClusters[3][fClno] = cells[j];
+           }
+       }
+
+
+    }
+  }
+}
+
+
+// ------------------------------------------------------------------------ //
+
+void AliPMDClusteringV2::ClustDetails(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 &rcl, Double_t &rcs,
+                                     Double_t &cells)
+{
+  // function begins
+  //
+  
+  const Int_t kndim1 = 4500;
+  const Int_t kndim2 = 10;
+  const Int_t kndim3 = 100;
+
+  int i, j, k, i1, i2;
+  int cluster[kndim1][kndim2];
+  
+  double x1, y1, x2, y2, rr;
+  double sumx, sumy, sumxy, sumxx;
+  double sum, sum1, sumyy;
+  double b, c, r1, r2;
+
+  double xx[kndim1], yy[kndim1], zz[kndim1];
+  double xxc[kndim1], yyc[kndim1];
+
+  double str[kndim1];
+
+  double str1[kndim1];
+  double xcl[kndim1], ycl[kndim1], cln[kndim1];
+  double clustcell[kndim1][kndim3];
+
+  for(i=0; i<=nclust; i++){
+   xxc[i]=*(&xc+i); 
+   yyc[i]=*(&yc+i); 
+   str[i]=0.; 
+   str1[i]=0.;
+  }
+  for(i=0; i<=ncell; i++){
+    xx[i]=*(&x+i); 
+    yy[i]=*(&y+i); 
+    zz[i]=*(&z+i);
+  }
+  // INITIALIZE 
+  for(i=0; i<4500; i++){
+    for(j=0; j<100; j++){
+      clustcell[i][j]=0.;
+    }
+  }
+
+  // INITIALIZE
+  for(i=0;i<4500;i++){
+    for(j=0;j<10;j++){
+      cluster[i][j]=0;
+    }
+  }
+
+
+  if(nclust > 0){
+    // more than one cluster
+    // checking cells shared between several  clusters.
+    // First check if the cell is within
+    // one cell unit ( nearest neighbour). Else, 
+    // if it is within 1.74 cell units ( next nearest )
+    // Else if it is upto 2 cell units etc.
+
+    for (i=0; i<=ncell; i++){
+      x1            = xx[i];
+      y1            = yy[i];
+      cluster[i][0] = 0;
+      // distance <= 1 cell unit
+      for(j=0; j<=nclust; j++)
+       {
+         x2 = xxc[j];
+         y2 = yyc[j];
+         rr = Distance(x1, y1, x2, y2);
+         if(rr <= 1.)
+           {
+             cluster[i][0]++;
+             i1             = cluster[i][0];
+             cluster[i][i1] = j;
+           }
+       }
+      // next nearest neighbour
+      if(cluster[i][0] == 0)
+       {
+         for(j=0; j<=nclust; j++)
+           {
+             x2 = xxc[j];
+             y2 = yyc[j];
+             rr = Distance(x1, y1, x2, y2);
+             if(rr <= sqrt(3.))
+               {
+                 cluster[i][0]++;
+                 i1             = cluster[i][0];
+                 cluster[i][i1] = j;
+               }
+           }
+       }
+      // next-to-next nearest neighbour
+      if(cluster[i][0] == 0)
+       {
+         for(j=0; j<=nclust; j++)
+           {
+             x2 = xxc[j];
+             y2 = yyc[j];
+             rr = Distance(x1, y1, x2, y2);
+             if(rr <= 2.)
+               {
+                 cluster[i][0]++;
+                 i1 = cluster[i][0];
+                 cluster[i][i1] = j;
+               }
+           }
+       }
+      // one more
+      if(cluster[i][0] == 0)
+       {
+         for(j=0; j<=nclust; j++)
+           {
+             x2 = xxc[j];
+             y2 = yyc[j];
+             rr = Distance(x1, y1, x2, y2);
+             if(rr <= 2.7)
+               {
+                 cluster[i][0]++;
+                 i1 = cluster[i][0];
+                 cluster[i][i1] = j;
+               }
+           }
+       }
+    }
+
+
+    // computing cluster strength. Some cells are shared.
+    for(i=0; i<=ncell; i++){
+      if(cluster[i][0] != 0){
+       i1 = cluster[i][0];
+       for(j=1; j<=i1; j++){
+         i2      = cluster[i][j];
+         str[i2] = str[i2]+zz[i]/i1;
+       }
+      }
+    }
+
+    for(k=0; k<5; k++)
+      {
+       for(i=0; i<=ncell; i++)
+         {
+           if(cluster[i][0] != 0)
+             {
+               i1=cluster[i][0];
+               sum=0.;
+               for(j=1; j<=i1; j++)
+                 {
+                   sum=sum+str[cluster[i][j]];
+                 }
+
+               for(j=1; j<=i1; j++)
+                 {
+                   i2 = cluster[i][j]; 
+                   str1[i2]         = str1[i2] + zz[i]*str[i2]/sum;
+                   clustcell[i2][i] = zz[i]*str[i2]/sum;
+                 }
+             }
+         }
+
+
+         for(j=0; j<=nclust; j++)
+           {
+             str[j]=str1[j];
+             str1[j]=0.;
+           }
+      }
+
+    for(i=0; i<=nclust; i++){
+      sumx = 0.;
+      sumy = 0.;
+      sum  = 0.;
+      sum1 = 0.;
+      for(j=0; j<=ncell; j++){
+       if(clustcell[i][j] != 0){
+         sumx = sumx+clustcell[i][j]*xx[j];
+         sumy = sumy+clustcell[i][j]*yy[j];
+         sum  = sum+clustcell[i][j];
+         sum1 = sum1+clustcell[i][j]/zz[j];
+       }
+      }
+      //***** xcl and ycl are cluster centroid positions ( center of gravity )
+
+      xcl[i] = sumx/sum;
+      ycl[i] = sumy/sum;
+      cln[i] = sum1;
+      sumxx = 0.;
+      sumyy = 0.;
+      sumxy = 0.;
+      for(j=0; j<=ncell; j++){
+       sumxx = sumxx+clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
+       sumyy = sumyy+clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
+       sumxy = sumxy+clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
+      }
+      b = sumxx+sumyy;
+      c = sumxx*sumyy-sumxy*sumxy;
+      // ******************r1 and r2 are major and minor axes ( r1 > r2 ). 
+      r1 = b/2.+sqrt(b*b/4.-c);
+      r2 = b/2.-sqrt(b*b/4.-c);
+      // final assignments to proper external variables
+      *(&xc + i) = xcl[i];
+      *(&yc + i) = ycl[i];
+      *(&zc + i) = str[i];
+      *(&cells + i) = cln[i];
+      *(&rcl+i) = r1;
+      *(&rcs+i) = r2;
+    }
+  }else{
+    sumx = 0.;
+    sumy = 0.;
+    sum  = 0.;
+    sum1 = 0.;
+    i    = 0;
+    for(j=0; j<=ncell; j++){
+      sumx = sumx+zz[j]*xx[j];
+      sumy = sumy+zz[j]*yy[j];
+      sum  = sum+zz[j];
+      sum1 = sum1+1.;
+    }
+    xcl[i] = sumx/sum;
+    ycl[i] = sumy/sum;
+    cln[i] = sum1;
+    sumxx  = 0.;
+    sumyy  = 0.;
+    sumxy  = 0.;
+    for(j=0; j<=ncell; j++){
+      sumxx = sumxx+clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
+      sumyy = sumyy+clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
+      sumxy = sumxy+clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
+    }
+    b  = sumxx+sumyy;
+    c  = sumxx*sumyy-sumxy*sumxy;
+    r1 = b/2.+sqrt(b*b/4.-c);
+    r2 = b/2.-sqrt(b*b/4.-c);
+    // final assignments
+    *(&xc + i)    = xcl[i];
+    *(&yc + i)    = ycl[i];
+    *(&zc + i)    = str[i];
+    *(&cells + i) = cln[i];
+    *(&rcl+i)     = r1;
+    *(&rcs+i)     = r2;
+  }
+}
+
+// ------------------------------------------------------------------------ //
+Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
+                                     Double_t x2, Double_t y2)
+{
+  return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
+}
+// ------------------------------------------------------------------------ //
+void AliPMDClusteringV2::SetEdepCut(Float_t decut)
+{
+  fCutoff = decut;
+}
+// ------------------------------------------------------------------------ //
diff --git a/PMD/AliPMDClusteringV2.h b/PMD/AliPMDClusteringV2.h
new file mode 100644 (file)
index 0000000..9b7bf9e
--- /dev/null
@@ -0,0 +1,105 @@
+#ifndef ALIPMDCLUSTERINGV2_H
+#define ALIPMDCLUSTERINGV2_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//-----------------------------------------------------//
+//                                                     //
+//  Header File : PMDClusteringV2.h,                   //
+//                                                     //
+//  clustering code for alice pmd                      //
+//                                                     //
+//-----------------------------------------------------//
+/* --------------------------------------------------------------------
+   Code developed by S. C. Phatak, Institute of Physics,
+   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 d[ndimx][ndimy] and cluster information is in array
+   clusters[5][5000]. integer clno gives total number of clusters in the
+   supermodule.
+   d, clno  and clusters are 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 "Rtypes.h"
+
+class TObjArray;
+class AliPMDcluster;
+class AliPMDClusteringV2 : public AliPMDClustering
+{
+
+ public:
+  AliPMDClusteringV2();
+  virtual ~AliPMDClusteringV2();
+
+  void     DoClust(Int_t idet, Int_t ismn, Double_t celladc[][96],
+                  TObjArray *pmdcont);
+  void     Order();
+  
+  Int_t    CrClust(Double_t ave, Double_t cutoff, Int_t nmx1);
+  void     RefClust(Int_t incr);
+
+  void     ClustDetails(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 &rcl, Double_t &rcs, Double_t &cells);
+  Double_t Distance(Double_t x1, Double_t y1,
+                   Double_t x2, Double_t y2);
+  void     SetEdepCut(Float_t decut);
+  
+ protected:
+
+  static const Double_t fgkSqroot3by2;  // fgkSqroot3by2 = sqrt(3.)/2.
+
+
+  enum {
+    kNMX    = 11424,
+    kNDIMX  = 119,
+    kNDIMY  = 96
+  };
+  /*
+    kNMX   : # of cells in a supermodule
+    kNDIMX : maximum number of cells along x direction (origin at one corner)
+    kNDIMY : maximum number of cells along axis at 60 degrees with x axis
+  */
+
+  Double_t fEdepCell[kNDIMX][kNDIMY]; //energy(ADC) in each cell of the supermodule
+  Double_t fClusters[6][5000]; // Cluster informations
+  Int_t    fClno;   // number of clusters in a supermodule
+
+  /*
+    clusters[0][i] --- x position of the cluster center
+    clusters[1][i] --- y position of the cluster center
+    clusters[2][i] --- total energy in the cluster
+    clusters[3][i] --- number of cells forming the cluster
+                       ( possibly fractional )
+    clusters[4][i] --- cluster sigma x
+    clusters[5][i] --- cluster sigma y
+  */
+
+  Int_t    fIord[2][kNMX];             // ordered list of i and j according to
+                                       // decreasing energy dep.
+  Int_t    fInfocl[2][kNDIMX][kNDIMY]; // cellwise information on the 
+                                       // cluster to which the cell
+  Int_t    fInfcl[3][kNMX];            // cluster information [0][i]
+                                       // -- cluster number
+  Double_t fCoord[2][kNDIMX][kNDIMY];
+
+  /*
+    fIord --- ordered list of i and j according to decreasing energy dep.
+    fInfocl --- cellwise information on the cluster to which the cell
+    belongs and whether it has largest energy dep. or not
+    ( now redundant - probably )
+    fInfcl ---  cluster information [0][i] -- cluster number
+    [1][i] -- i of the cell
+    [2][i] -- j of the cell
+    coord --- x and y coordinates of center of each cell
+  */
+
+  Float_t fCutoff; // Energy(ADC) cutoff per cell before clustering
+
+  ClassDef(AliPMDClusteringV2,0) // Does clustering for PMD
+};
+#endif