47 float e,
x,
y,
xx, yy, xy;
69 vector<EmcModule>::iterator ph;
74 if ((*ph).ich == ich)
return (*ph).amp;
85 vector<EmcModule>::iterator ph;
98 vector<EmcModule>::iterator ph;
103 if ((*ph).ich == ich)
return (*ph).tof;
114 vector<EmcModule>::iterator ph;
132 float e,
x,
y,
xx, yy, xy;
133 float ecore, ecorecorr;
143 const float thresh = 0.01;
145 vector<EmcModule>::iterator ph;
146 float xcg, ycg,
xx, xy, yy;
165 if (et > thresh) es += (*ph).amp;
176 float et, xcg, ycg,
xx, xy, yy;
178 int ix0, iy0, isx, isy;
181 ix0 =
int(xcg + 0.5);
182 iy0 =
int(ycg + 0.5);
185 if (xcg - ix0 < 0) isx = -1;
187 if (ycg - iy0 < 0) isy = -1;
194 return (e1 + e2 + e3 + e4);
201 float et, xcg, ycg,
xx, xy, yy;
202 int ich, ix0, iy0, nhit;
206 if (nhit <= 0)
return 0;
209 ix0 =
int(xcg + 0.5);
210 iy0 =
int(ycg + 0.5);
221 vector<EmcModule>::iterator ph;
236 if (
abs(idx) <= 1 &&
abs(idy) <= 1) es += (*ph).amp;
247 vector<EmcModule>::iterator ph;
259 if ((*ph).amp > emax)
305 int ixypk, ixpk, iypk,
in, nh, ic;
315 vector<EmcModule>::iterator ph;
316 vector<EmcModule> hl;
323 if (nhit <= 0)
return 0;
329 while (ph !=
fHitList.end()) *vv++ = *ph++;
338 for (ic = 0; ic < nhit; ic++)
340 float amp = hlist[ic].
amp;
343 int ich = hlist[ic].
ich;
352 while ((in < nhit) && (hlist[in].ich <= ichmax))
361 while ((in >= 0) && (hlist[
in].
ich >= ichmin))
371 cout <<
"!!! Error in EmcCluster::GetSubClusters(): too many peaks in a cluster (>"
373 <<
"). May need tower energy threshold increase for clustering." << endl;
395 for (
int ich = 0; ich < nhit; ich++) hl.push_back(hlist[ich]);
397 PkList->push_back(peak);
400 ppeaks->push_back(hlist[PeakCh[0]]);
410 for (ipk = 0; ipk < npk; ipk++)
412 Energy[ipk] =
new float[nhit];
414 totEnergy =
new float[nhit];
415 tmpEnergy =
new float[nhit];
416 for (
int i = 0; i < nhit; ++i)
428 for (ipk = 0; ipk < npk; ipk++)
431 if (iter > 0) ratio = Energy[ipk][ic] / totEnergy[ic];
433 ixypk = hlist[ic].
ich;
446 for (in = ic + 1; in < nhit; in++)
453 if (ixy > ichmax)
break;
459 if (iter > 0) ratio = Energy[ipk][
in] / totEnergy[
in];
462 xpk[ipk] += eg *
idx;
463 ypk[ipk] += eg * idy;
471 for (in = ic - 1; in >= 0; in--)
478 if (ixy < ichmin)
break;
484 if (iter > 0) ratio = Energy[ipk][
in] / totEnergy[
in];
487 xpk[ipk] += eg *
idx;
488 ypk[ipk] += eg * idy;
493 xpk[ipk] = xpk[ipk] / epk[ipk] + ixpk;
494 ypk[ipk] = ypk[ipk] / epk[ipk] + iypk;
497 for (in = 0; in < nhit; in++)
507 if (
ABS(dx) < 2.5 &&
ABS(dy) < 2.5)
516 float flim = 0.000001;
517 for (in = 0; in < nhit; in++)
519 if (tmpEnergy[in] > flim)
521 totEnergy[
in] = tmpEnergy[
in];
525 totEnergy[
in] = flim;
533 for (ipk = 0; ipk < npk; ipk++)
538 for (in = 0; in < nhit; in++)
540 if (tmpEnergy[in] > 0)
543 a = hlist[
in].
amp * Energy[ipk][
in] / tmpEnergy[
in];
583 for (ipk = 0; ipk < npk; ipk++)
591 for (in = 0; in < nhit; in++)
594 for (ig = igmpk1[ipk]; ig <= igmpk2[ipk]; ig++)
603 Energy[ipk][
in] +=
a;
611 for (ipk = 0; ipk < npk; ipk++)
614 for (in = 0; in < nhit; in++)
616 if (tmpEnergy[in] > 0)
619 a = hlist[
in].
amp * Energy[ipk][
in] / tmpEnergy[
in];
633 ppeaks->push_back(hlist[PeakCh[ipk]]);
635 for (in = 0; in < nh; in++) hl.push_back(phit[in]);
637 PkList->push_back(peak);
644 for (ipk = 0; ipk < npk; ipk++)
delete[] Energy[ipk];