+
+
+int AliPMDClustering::crclust(double ave, double cutoff, int nmx1)
+{
+ 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 infocl[2][ndimx][ndimy]
+
+ if (fDebug == 1)
+ {
+ printf(" *** Inside crclust ** nmx = %d nmx1 = %d ndimx = %d ndimy = %d ave = %f cutoff = %f\n",
+ nmx,nmx1,ndimx,ndimy,ave,cutoff);
+ }
+ for (j=0; j < ndimx; j++){
+ for(k=0; k < ndimy; k++){
+ infocl[0][j][k] = 0;
+ infocl[1][j][k] = 0;
+ }
+ }
+ for(i=0; i < nmx; i++){
+ infcl[0][i] = -1;
+ id1=iord[0][i];
+ id2=iord[1][i];
+ if(d[id1][id2] <= cutoff){infocl[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=iord[0][icell];
+ id2=iord[1][icell];
+ if(infocl[0][id1][id2] == 0 ){
+ // ---------------------------------------------------------------
+ // icl -- cluster #, numcell -- # of cells in it, clust -- stores
+ // coordinates of the cells in a cluster, infocl[0][i1][i2] is 1 for
+ // primary and 2 for secondary cells,
+ // infocl[1][i1][i2] stores cluster #
+ // ---------------------------------------------------------------
+ icl=icl+1;
+ numcell=0;
+ cellcount=cellcount+1;
+ infocl[0][id1][id2]=1;
+ infocl[1][id1][id2]=icl;
+ infcl[0][cellcount]=icl;
+ infcl[1][cellcount]=id1;
+ infcl[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 < 72) && (jd2 >= 0 && jd2 < 72) &&
+ if( (jd1 >= 0 && jd1 < ndimx) && (jd2 >= 0 && jd2 < ndimy) &&
+ infocl[0][jd1][jd2] == 0){
+ numcell=numcell+1;
+ infocl[0][jd1][jd2]=2;
+ infocl[1][jd1][jd2]=icl;
+ clust[0][numcell]=jd1;
+ clust[1][numcell]=jd2;
+ cellcount=cellcount+1;
+ infcl[0][cellcount]=icl;
+ infcl[1][cellcount]=jd1;
+ infcl[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 < 72) && (jd2 >= 0 && jd2 < 72) &&
+ if( (jd1 >= 0 && jd1 < ndimx) && (jd2 >= 0 && jd2 < ndimy) &&
+ infocl[0][jd1][jd2] == 0 ){
+ infocl[0][jd1][jd2]=2;
+ infocl[1][jd1][jd2]=icl;
+ numcell=numcell+1;
+ clust[0][numcell]=jd1;
+ clust[1][numcell]=jd2;
+ cellcount=cellcount+1;
+ infcl[0][cellcount]=icl;
+ infcl[1][cellcount]=jd1;
+ infcl[2][cellcount]=jd2;
+ }
+ }
+ }
+ }
+ }
+ }
+ // for(icell=0; icell<=cellcount; icell++){
+ // ofl0 << infcl[0][icell] << " " << infcl[1][icell] << " " <<
+ // infcl[2][icell] << endl;
+ // }
+ return cellcount;
+}
+
+void AliPMDClustering::refclust(int incr)