Breaking Superclusters Into Smaller Clusters

The superclusters constructed in function crclust generally have large number of cells. This is particularly so in case of large multiplicities. If the particle multiplicity is small, each particle would produce a cluster in the detector and these clusters will, in general, be non-overlapping. In that case, each supercluster would consist of few cells and each supercluster may corespond to a particle. In such a case, there is no need of breaking the superclusters further. But when the multiplicity is large, the clusters produced by each particle may overlap, thus giving rise to large superclusters. The aim here is to devise a scheme to identify individual clusters when one has overlapping clusters forming a supercluster.

Apart from identifing clusters we also want to know certain other properties of the cluster. These are

1. The position of the cluster center. This will consist of x and y coordinates in a coordinate system in which the origin is at one of the corners of the supermodule.
2. The strength of the cluster which is the sum of edep's of all the cells belonging to the cluster.
3. The number of cells belonging to the cluster.
4. The width or extent or size of the cluster. Since the cluster is in a plane, its shape may be elongated. So this would mean there would be two sizes, a long one and a short one.

When the supercluster consists of one cell, its coordinates are taken to be the center of the cell. Its strength is just the edep of the cell. The number of cells of the cluster are ( obviously ) 1 and the width is taken to be zero ( both long and short ).

For a supercluster having two cells, the cluster center is the center of gravity of the two cells. That is, if the positions of the two cells are $(x_1, y_1)$ and $(x_2, y_2)$ and their adc's are $z_1$ and $z_2$ respectively, the cluster centers $(x_c, y_c)$ is given by $x_c = \frac{z1 x1 + z2 x2}{z1+z2}$ and $y_c = \frac{z1 y1 + z2 y2}{z1+z2}$. The cluster strength is $z_1+z_2$. Number of cells is 2, the width along the line joining the two cells is $\frac{z_1 z_2}{(z_1+z_2)^2}$ and the width perpendicular to the line joining the two cells is zero.

If the supercluster has more than two cells it may be possible to break it into a number of clusters. For breaking up the supercluster let us assume that the supercluster consists of overlapping clusters. Then, to begin with, the maxima ( peaks ) in edep's may be identified with the centers of the clusters. This ansatz is a reasonable starting point but needs to be refined. Further, we don't want to identify each maximum in edep with a cluster. For example, a weak peak very close to a strong peak may be due to fluctuation in edep. In any case, we should not be having two cluster centers very close ( distance $\sim$ 1 cell unit ). Thus we begin the determination of cluster centers from the cell with largest edep ( obvious cluster center ) and work downwards. Several conditions are imposed for a cell to be accepted as a cluster center. These are

1. The proposed center should be at least one cell unit away from previously determined cluster centers. That is, neighbouring cells cannot be cluster centers.
2. If the distance between the proposed center and the previously accpted center is between 1 and 2, the strength of the ( new ) cell should be larger than 25% of the previously accepted center cell. This is to ensure that flustuations do not give rise to clusters. The number ( 25% ) is ad hoc and needs to be fine tuned. Note that this distance condition implies the two cells are next nearest neighbours.
3. If the distance between the proposed center and the previously accpted center is 2 cell units ( next-to-next nearest neighbour ), the cell strength should be larger than 10% of the strength of the previously accpted cell.
4. If the distance is larger than 2 cell units, it is accepted as new cell center.

It must be mentioned that the procedure outlined above is adhoc and needs to be fine tuned keeping the results of test beam in mind. Further more, this procedure determines only the cluster centers and even these need to be refined. We still need to compute the strength of each cluster, the number of cells belonging to it and its size. This is done in another function. Below is the listing of the function crclust which does the computations outlined above.

void refclust(int incr, int supmod, int idet){
  int i, j, k, i1, i2, id, icl, ncl[4500], iord[4500], itest; 
  int ihld, ncell, ipause;
  int ig, clustno, nsupcl, idum, neib[200], idd, lev1[20], lev2[20];
  double x[4500], y[4500], z[4500], x1, y1, z1, x2, y2, z2, dst, dist;
  double xc[4500], yc[4500], zc[4500], cells[4500], sum, rcl[4500], rcs[4500], rr;
  // 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
  clno=-1;
  nsupcl=-1;
  for(i=0; i<4500; i++){ncl[i]=-1;}
  for(i=0; i<=incr; i++){
    if(infcl[0][i] != nsupcl){ nsupcl=nsupcl+1; }
    ncl[nsupcl]=ncl[nsupcl]+1;
  }
  id=-1;
  icl=-1;
  for(i=0; i<=nsupcl; i++){
    if(ncl[i] == 0){ 
      // one  cell super-clusters --> single cluster
      // cluster center at the center of the cell
      // cluster radius = zero
      id=id+1;  icl=icl+1; clno=clno+1;  i1=infcl[1][id];  i2=infcl[2][id];
      clusters[0][clno]=coord[0][i1][i2];  clusters[1][clno]=coord[1][i1][i2];
      clusters[2][clno]=d[i1][i2];  clusters[3][clno]=1.;  
	clusters[4][clno]=0.; clusters[5][clno]=0.; // for one-cell cluster, sigma = 0.
    }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; clno=clno+1; 
      i1=infcl[1][id];  i2=infcl[2][id]; 
      x1=coord[0][i1][i2]; y1=coord[1][i1][i2];  z1=d[i1][i2];
      id=id+1; i1=infcl[1][id];  i2=infcl[2][id];
      x2=coord[0][i1][i2];  y2=coord[1][i1][i2];  z2=d[i1][i2];
      clusters[0][clno]=(x1*z1+x2*z2)/(z1+z2);  
      clusters[1][clno]=(y1*z1+y2*z2)/(z1+z2);
      clusters[2][clno]=z1+z2;  clusters[3][clno]=2.;  clusters[4][clno]=sqrt(z1*z2)/(z1+z2);
	clusters[5][clno]=0.; // sigma large nonzero, sigma small zero.
    }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
      i1=infcl[1][id];  i2=infcl[2][id];
      x[0]=coord[0][i1][i2];  y[0]=coord[1][i1][i2];  z[0]=d[i1][i2];
      iord[0]=0;
      for(j=1;j<=ncl[i];j++){
	id=id+1; i1=infcl[1][id];  i2=infcl[2][id]; iord[j]=j;
	x[j]=coord[0][i1][i2];  y[j]=coord[1][i1][i2]; z[j]=d[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=Dist(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]];
	}
      }
      centroidCal(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++){
	clno=clno+1; 
	clusters[0][clno]=xc[j]; 
	clusters[1][clno]=yc[j]; 
	clusters[2][clno]=zc[j];
	clusters[3][clno]=cells[j];
	clusters[4][clno]=rcl[j];
	clusters[5][clno]=rcs[j];
	if(clusters[3][clno] == 1){clusters[4][clno]=0.; clusters[5][clno]=0.;}
      }
    }
  }
}