- 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]
- for (j=0; j < 72; j++){
- for(k=0; k < 72; 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) &&
- 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) &&
- 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;
-}
-
-double AliPMDClustering::ranmar()
-{
- /* C==========================C*/
- /*===================================C==========================*/
- /* Universal random number generator proposed by Marsaglia and Zaman
- in report FSU-SCRI-87-50 */