26 fModules =
new std::vector<EmcModule>;
53 std::cout <<
"Warning from BEmcRec::LoadProfile(): No acton defined for shower profile evaluation; should be defined in a detector specific module " <<
Name() << std::endl;
58 std::ofstream outfile(fname);
59 if (!outfile.is_open())
61 std::cout <<
"Error in BEmcRec::PrintTowerGeometry(): Failed to open file "
62 << fname << std::endl;
65 outfile <<
"Number of bins:" << std::endl;
66 outfile <<
fNx <<
" " <<
fNy << std::endl;
67 outfile <<
"ix iy x y z dx0 dy0 dz0 dx1 dy1 dz1" << std::endl;
70 std::map<int, TowerGeom>::iterator
it;
71 for (
int iy = 0; iy <
fNy; iy++)
73 for (
int ix = 0; ix <
fNx; ix++)
80 outfile << ix <<
" " << iy <<
" " << geom.
Xcenter <<
" "
82 << geom.
dY[0] <<
" " << geom.
dZ[0] <<
" " << geom.
dX[1] <<
" "
83 << geom.
dY[1] <<
" " << geom.
dZ[1] << std::endl;
92 if (ix < 0 || ix >=
fNx || iy < 0 || iy >=
fNy)
return false;
94 int ich = iy *
fNx + ix;
95 std::map<int, TowerGeom>::iterator
it =
fTowerGeom.find(ich);
104 if (ix < 0 || ix >=
fNx || iy < 0 || iy >=
fNy)
return false;
110 geom.
dX[0] = geom.
dX[1] = 0;
111 geom.
dY[0] = geom.
dY[1] = 0;
112 geom.
dZ[0] = geom.
dZ[1] = 0;
114 int ich = iy *
fNx + ix;
124 std::cout <<
"Error in BEmcRec::CalculateTowerSize(): Tower geometry not well setup (NX = "
125 <<
fNx <<
")" << std::endl;
130 int idx[nb] = {0, 1, 0, -1, -1, 1, 1, -1};
131 int idy[nb] = {-1, 0, 1, 0, -1, -1, 1, 1};
133 std::map<int, TowerGeom>::iterator
it;
145 while (inx < nb && (idx[inx] == 0 || !
GetTowerGeometry(ix + idx[inx], iy + idy[inx], geomx))) inx++;
148 std::cout <<
"Error in BEmcRec::CompleteTowerGeometry(): Error when locating neighbour for (ix,iy)=("
149 << ix <<
"," << iy <<
")" << std::endl;
156 while (iny < nb && (idy[iny] == 0 || !
GetTowerGeometry(ix + idx[iny], iy + idy[iny], geomy))) iny++;
159 std::cout <<
"Error in BEmcRec::CompleteTowerGeometry(): Error when locating neighbour for (ix,iy)=("
160 << ix <<
"," << iy <<
")" << std::endl;
179 float& xA,
float& yA,
float& zA)
188 if (ix < 0 || ix >=
fNx)
190 std::cout <<
m_ThisName <<
" Error in BEmcRec::Tower2Global: wrong input x: " << ix << std::endl;
195 if (iy < 0 || iy >=
fNy)
197 std::cout <<
"Error in BEmcRec::Tower2Global: wrong input y: " << iy << std::endl;
207 int idx[4] = {1, 0, -1, 0};
208 int idy[4] = {0, 1, 0, -1};
210 while (ii < 4 && !
GetTowerGeometry(ix + idx[ii], iy + idy[ii], geom0)) ii++;
213 std::cout <<
"Error in BEmcRec::Tower2Global: can not identify neighbour for tower ("
214 << ix <<
"," << iy <<
")" << std::endl;
217 float Xc = geom0.
Xcenter - idx[ii] * geom0.
dX[0] - idy[ii] * geom0.
dX[1];
218 float Yc = geom0.
Ycenter - idx[ii] * geom0.
dY[0] - idy[ii] * geom0.
dY[1];
219 float Zc = geom0.
Zcenter - idx[ii] * geom0.
dZ[0] - idy[ii] * geom0.
dZ[1];
225 float xt = geom0.
Xcenter + (xC - ix) * geom0.
dX[0] + (yC - iy) * geom0.
dX[1];
226 float yt = geom0.
Ycenter + (xC - ix) * geom0.
dY[0] + (yC - iy) * geom0.
dY[1];
227 float zt = geom0.
Zcenter + (xC - ix) * geom0.
dZ[0] + (yC - iy) * geom0.
dZ[1];
240 int idist = ix2 - ix1;
243 int idistr =
fNx -
abs(idist);
244 if (idistr <
abs(idist))
258 float dist = x2 -
x1;
261 float distr =
fNx - fabs(dist);
262 if (distr <
abs(dist))
282 int next, ib, ie, iab, iae, last, LastCl, leng, ich;
288 std::vector<EmcModule>::iterator ph;
289 std::vector<EmcModule> hl;
291 (*fClusters).clear();
292 nhit = (*fModules).size();
294 if (nhit <= 0)
return 0;
303 LenCl =
new int[MaxLen];
309 ph = (*fModules).begin();
311 while (ph != (*fModules).end()) *vv++ = *ph++;
317 for (ich = 1; ich < nhit + 1; ich++)
319 if (ich < nhit) ia = vhit[ich].
ich;
323 if ((ia - vhit[ich - 1].ich > 1)
325 || (ia - ia /
fNx *
fNx == 0))
336 int* LenCltmp =
new int[MaxLen];
339 LenCl =
new int[MaxLen * 2];
347 LenCl[nCl - 1] = next - ib;
356 for (
int iCl = LastCl; iCl >= 0; iCl--)
360 if (iab - vhit[last].ich >
fNx)
goto new_ich;
361 for (
int ichc = last; ichc >= last - leng + 1; ichc--)
366 if ((vhit[ichc].ich +
fNx <= iae && vhit[ichc].ich +
fNx >= iab) || (
bCYL && (iae %
fNx ==
fNx - 1) && (iae - vhit[ichc].ich ==
fNx - 1))
371 CopyVector(&vhit[last + 1], &vhit[last + 1 - leng], ib - 1 - last);
375 for (
int i = iCl; i < nCl - 2; i++) LenCl[i] = LenCl[i + 1];
377 LenCl[nCl - 2] = LenCl[nCl - 1] + leng;
398 for (
int iCl = 0; iCl < nCl; iCl++)
402 for (ich = 0; ich < leng; ich++) hl.push_back(vhit[ib + ich]);
418 float& py,
float& pxx,
float& pyy,
float& pyx,
423 float a,
x,
y,
e,
xx, yy, yx;
424 std::vector<EmcModule>::iterator ph;
432 if (phit->empty())
return;
439 while (ph != phit->end())
449 if (emax <= 0)
return;
451 int iymax = ichmax /
fNx;
452 int ixmax = ichmax - iymax *
fNx;
463 while (ph != phit->end())
468 int iy = ph->ich /
fNx;
469 int ix = ph->ich - iy *
fNx;
471 int idy = iy - iymax;
494 while (x < -0.5) x += float(fNx);
495 while (x >= fNx - 0.5) x -= float(fNx);
523 float fPpar1, fPpar2, fPpar3, fPpar4;
558 dx = fabs(xc - fPshiftx);
559 dy = fabs(yc - fPshifty);
560 r2 = dx * dx + dy *
dy;
563 double e = fPpar1 * exp(-r3 / fPpar2) + fPpar3 * exp(-r1 / fPpar4);
574 while (xcg < -0.5) xcg += float(
fNx);
575 while (xcg >=
fNx - 0.5) xcg -= float(
fNx);
577 int ixcg =
int(xcg + 0.5);
578 int iycg =
int(ycg + 0.5);
579 float ddx = fabs(xcg - ixcg);
580 float ddy = fabs(ycg - iycg);
589 if (xcg - ixcg < 0) isx = -1;
591 if (ycg - iycg < 0) isy = -1;
594 int idy = (iy - iycg) * isy;
597 if (idx == 0 && idy == 0)
599 else if (idx == 1 && idy == 0)
601 else if (idx == 1 && idy == 1)
603 else if (idx == 0 && idy == 1)
609 float dy = fabs(iy - ycg);
610 float rr = sqrt(dx * dx + dy * dy);
616 for (
int ip = 0; ip < 4; ip++)
624 eout = (ep[1] + ep[2]) / 2. + ep[3];
626 eout = (ep[0] - ep[2]) / 2. - ep[3];
628 eout = (ep[0] - ep[1]) / 2. - ep[3];
633 if (eout < 0) eout = 1
e-6;
642 int nn = plist->size();
643 if (nn <= 0)
return 0;
645 for (
int i = 0; i <
nn; i++)
647 int ich = (*plist)[i].ich;
650 if (iy == iyt && iz == izt)
652 return (*plist)[i].amp;
659 float BEmcRec::GetProb(std::vector<EmcModule> HitList,
float en,
float xg,
float yg,
float zg,
float& chi2,
int& ndf)
676 int nn = HitList.size();
677 if (nn <= 0)
return -1;
687 Momenta(&HitList, etot, zcg, ycg, zz, yy, yz, thresh);
689 int iz0cg =
int(zcg + 0.5);
690 int iy0cg =
int(ycg + 0.5);
691 float ddz = fabs(zcg - iz0cg);
692 float ddy = fabs(ycg - iy0cg);
695 if (zcg - iz0cg < 0) isz = -1;
697 if (ycg - iy0cg < 0) isy = -1;
707 if (e1 < thresh) e1 = 0;
708 if (e2 < thresh) e2 = 0;
709 if (e3 < thresh) e3 = 0;
710 if (e4 < thresh) e4 = 0;
712 float e1t = (e1 + e2 + e3 +
e4) / etot;
713 float e2t = (e1 + e2 - e3 -
e4) / etot;
714 float e3t = (e1 - e2 - e3 +
e4) / etot;
715 float e4t = (
e3) / etot;
722 for (
int ip = 0; ip <
NP; ip++)
731 err[ip] = sqrt(err[ip] * err[ip] + 4 * enoise * enoise / etot / etot);
735 err[ip] = sqrt(err[ip] * err[ip] + 1 * enoise * enoise / etot / etot);
740 chi2 += (ep[0] - e1t) * (ep[0] - e1t) / err[0] / err[0];
741 chi2 += (ep[1] - e2t) * (ep[1] - e2t) / err[1] / err[1];
742 chi2 += (ep[2] - e3t) * (ep[2] - e3t) / err[2] / err[2];
743 chi2 += (ep[3] - e4t) * (ep[3] - e4t) / err[3] / err[3];
748 float prob = TMath::Prob(chi2, ndf);
758 return (static_cast<const EmcModule*>(h1)->ich - static_cast<const EmcModule*>(h2)->ich);
765 float amp1 =
static_cast<const EmcModule*
>(
h1)->amp;
766 float amp2 =
static_cast<const EmcModule*
>(
h2)->amp;
767 return (amp1 < amp2) ? 1 : (amp1 > amp2) ? -1
776 for (
int i = 0; i <
N; i++) *p++ = 0;
784 for (
int i = 0; i <
N; i++) *p++ = 0;
792 for (
int i = 0; i <
N; i++)
805 for (
int i = 0; i <
N; i++) to[i] = from[i];
816 for (
int i = 0; i <
N; i++)