3.
Search for clusters: convert to ÔcmÕ, search
for 2-dim clusters, produces sorted list of clusters with energy above Õcluster_thresholdÕ
for(edge=0; edge<3; edge++) {
for(id=0;
id<npeak[edge]; id++) {
peak[edge][id].coord = peak[edge][id].coord * geom->edge[edge] /
36.;
peak[edge][id].width = peak[edge][id].width * geom->edge[edge] /
36.;
peak[edge][id].tmp = 0.0;
} }
npsble = -1;
for(u=0; u<npeak[0]; u++) {
for(v=0; v<npeak[1];
v++) {
for(w=0;
w<npeak[2]; w++) {
dltz = peak[0][u].coord/geom->edge[0] +
peak[1][v].coord/geom->edge[1] + peak[2][w].coord/geom->edge[2];
spread = ( peak[0][u].width + peak[1][v].width + peak[2][w].width )
*2.*sqrt(12.);
sum_edge = geom->edge[0] + geom->edge[1] + geom->edge[2];
if(ABS(dltz-2.)*sum_edge < spread) {
npsble++;
if(npsble ==
NHIT) return(-1);
hit[npsble].sector = sector + 1; /* in hit[] sector = 1..6 */
hit[npsble].layer = io + 1; /* in hit[] io = 1..2 */
hit[npsble].i = geom->h *
( peak[0][u].coord/geom->edge[0] - peak[1][v].coord/geom->edge[1]
- peak[2][w].coord/geom->edge[2] ) / 2. + geom->h2;
hit[npsble].j = geom->edge[1] *
( peak[2][w].coord/geom->edge[2] - peak[1][v].coord/geom->edge[1]
) / 2.;
hit[npsble].k = geom->d;
hit[npsble].di = sqrt( peak[0][u].width * peak[0][u].width /
(geom->edge[0] * geom->edge[0]) +
peak[1][v].width * peak[1][v].width / (geom->edge[1] *
geom->edge[1]) +
peak[2][w].width * peak[2][w].width / (geom->edge[2] *
geom->edge[2]) ) * geom->h / 2.;
hit[npsble].dj = sqrt( peak[2][w].width * peak[2][w].width /
(geom->edge[2] * geom->edge[2]) +
peak[1][v].width * peak[1][v].width / (geom->edge[1] *
geom->edge[1]) ) * geom->edge[1] / 2.;
/*hit[npsble].width = sqrt( hit[npsble].di*hit[npsble].di +
hit[npsble].dj*hit[npsble].dj );*/
hit[npsble].energy = peak[0][u].energy + peak[1][v].energy +
peak[2][w].energy;
hit[npsble].peak1[0] = u;
hit[npsble].peak1[1] = v;
hit[npsble].peak1[2] = w;
hit[npsble].peakn[0] = 1; hit[npsble].peakn[1] = 1; hit[npsble].peakn[2]
= 1; /* temporary one peak per layer only */
peak[0][u].tmp = peak[0][u].tmp + hit[npsble].energy;
peak[1][v].tmp = peak[1][v].tmp + hit[npsble].energy;
peak[2][w].tmp = peak[2][w].tmp + hit[npsble].energy;
}
} } }
qsort((void *)hit, npsble,
sizeof(ECHit), (int (*) (const void *, const void *))hit_compare);
nhit = 0; while(hit[nhit].energy >
th && nhit < NHIT && nhit < npsble) nhit++;