Calculate energy, time and chi2 for one channel using fitted pulse shape.
Method uses HFITV.
627 if (
rms > 0.0)
break;
631 if (
rms > 0.0)
break;
635 if (
rms > 0.0)
break;
640 }
else if (
igain == 1) {
644 if (
rms > 0.0)
break;
648 if (
rms > 0.0)
break;
652 if (
rms > 0.0)
break;
678 <<
" ifit1=" << ifit1
679 <<
" ifit2=" << ifit1 + nfit - 1
686 std::vector<float> samples =
digit->samples();
687 double maxsamp = 0.0;
693 bool no_signal =
true;
694 pedestal = samples[0];
695 const double delta = 1.e-6;
697 for (
int isamp = 0; isamp < nfit; ++isamp) {
698 int j = isamp + ifit1;
700 yvec0[isamp] = samples[j];
701 if (no_signal) no_signal = (fabs(samples[j] - pedestal) < delta);
709 ATH_MSG_VERBOSE (
"Saturated ADC value yvec0[" << isamp <<
"]=" << yvec0[isamp]
713 ATH_MSG_VERBOSE (
"Zero ADC value yvec0[" << isamp <<
"]=" << yvec0[isamp]
714 <<
" RMS=" << eyvec[isamp] );
718 if (yvec0[isamp] > maxsamp) {
721 maxsamp = yvec0[isamp];
724 if (yvec0[isamp] < minsamp) {
725 minsamp = yvec0[isamp];
729 <<
", xvec=" << xvec[isamp]
730 <<
", yvec0=" << yvec0[isamp]
731 <<
", eyvec=" << eyvec[isamp] );
740 double pedg = (yvec0[0] + yvec0[nfit - 1]) / 2.0;
763 <<
" ipeakMax=" << ipeakMax
764 <<
" ipeakMin=" << ipeakMin
765 <<
" fixedTime=" << ((fixedTime) ?
"true" :
"false") );
767 const std::vector<double>* tpulse = &
m_dummy;
768 const std::vector<double>* ypulse = &
m_dummy;
769 const std::vector<double>* tdpulse = &
m_dummy;
770 const std::vector<double>* dpulse = &
m_dummy;
771 const std::vector<double>* tleak = &
m_dummy;
772 const std::vector<double>* yleak = &
m_dummy;
773 const std::vector<double>* tdleak = &
m_dummy;
774 const std::vector<double>* dleak = &
m_dummy;
847 double gval, gpval,
sy, syg, sygp, sg, sgp, sgg, sgpgp, sggp, serr, err2;
848 double dgg0, dgg, dggp, dgpgp, dyg, dygp,
dg, dc,
xd;
849 double sllp, sylp, slplp, dleakage, leakage;
850 double fixtau = 0.0, fixped = 0.0, fixampl = 0.0, fixchi2 =
MAX_CHI2;
851 double leaktau = 0.0, leakped = 0.0, leakampl = 0.0, leakchi2 =
MAX_CHI2;
852 double cistau = 0.0, cisped = 0.0, cisampl = 0.0, cisatau = 0.0, cischi2 =
MAX_CHI2;
861 if (tempChi2 < oldchi2) oldchi2 = tempChi2;
878 for (
int isamp = 0; isamp < nfit; ++isamp) {
879 ATH_MSG_VERBOSE (
"Lo gain leakage xvec[" << isamp <<
"]=" << xvec[isamp] );
881 leakage =
pulse(xvec[isamp], tleak, yleak);
882 dleakage =
derivative(xvec[isamp], tdleak, dleak);
885 yvec[isamp] = yvec0[isamp] - pedg;
888 <<
" yvec0[" << isamp <<
"]=" << yvec0[isamp]
889 <<
" pedg=" << pedg );
891 sllp += leakage * dleakage;
892 sylp +=
yvec[isamp] * dleakage;
893 slplp += dleakage * dleakage;
896 if (fabs(slplp) >
EPS_DG && !fixedTime) {
897 leaktau = (sllp - sylp) / slplp;
908 <<
" slplp=" << slplp
909 <<
" leaktau=" << leaktau );
919 for (
int isamp = 0; isamp < nfit; ++isamp) {
920 leakage =
pulse(xvec[isamp], tleak, yleak);
922 yvec[isamp] = yvec0[isamp] - leakage;
925 <<
" yvec0[" << isamp <<
"]=" << yvec0[isamp]
926 <<
" leakage=" << leakage );
928 err2 = eyvec[isamp] * eyvec[isamp];
931 syg +=
yvec[isamp] * gval / err2;
932 sgg += gval * gval / err2;
935 dgg0 = sg * sg - serr * sgg;
936 if (fabs(dgg0) >
EPS_DG) {
937 leakampl = (
sy * sg - serr * syg) / dgg0;
938 leakped = (syg * sg -
sy * sgg) / dgg0;
951 for (
int isamp = 0; isamp < nfit; ++isamp) {
953 leakage =
pulse(xvec[isamp], tleak, yleak);
954 xd = yvec0[isamp] - (gval + leakage);
955 leakchi2 = leakchi2 + (
xd *
xd) / (eyvec[isamp] * eyvec[isamp]);
958 <<
" yvec0[" << isamp <<
"]=" << yvec0[isamp]
960 <<
", leakage=" << leakage
963 leakchi2 = leakchi2 / (nfit - 3.0);
966 <<
" leakped=" << leakped
967 <<
" leakampl=" << leakampl
968 <<
" leakchi2=" << leakchi2 );
975 for (
int isamp = 0; isamp < nfit; ++isamp) {
976 leakage =
pulse(xvec[isamp], tleak, yleak);
979 yvec[isamp] = yvec0[isamp] - leakage;
982 <<
" yvec0[" << isamp <<
"]=" << yvec0[isamp]
983 <<
" leakage=" << leakage
984 <<
" yvec[" << isamp <<
"]=" <<
yvec[isamp] );
998 for (
int isamp = 0; isamp < nfit; ++isamp) {
1000 gpval =
derivative(xvec[isamp], tdpulse, dpulse);
1001 err2 = eyvec[isamp] * eyvec[isamp];
1002 sy +=
yvec[isamp] / err2;
1004 sgp += gpval / err2;
1005 syg +=
yvec[isamp] * gval / err2;
1006 sygp +=
yvec[isamp] * gpval / err2;
1007 sgg += gval * gval / err2;
1008 sggp += gval * gpval / err2;
1009 sgpgp += gpval * gpval / err2;
1012 dgg = sgg - sg * sg / serr;
1013 dggp = sggp - sg * sgp / serr;
1014 dgpgp = sgpgp - sgp * sgp / serr;
1015 dyg = syg -
sy * sg / serr;
1016 dygp = sygp -
sy * sgp / serr;
1017 dg = dgg * dgpgp - dggp * dggp;
1020 cisampl = (dyg * dgpgp - dygp * dggp) /
dg;
1021 cisatau = (dyg * dggp - dygp * dgg) /
dg;
1023 - (dyg * dgpgp * sg - dygp * dggp * sg + dyg * dggp * sgp - dygp * dgg * sgp) /
dg)
1026 if (fabs(cisampl) >
EPS_DG) {
1027 cistau = cisatau / cisampl;
1049 <<
" sgpgp=" << sgpgp <<
endmsg;
1053 <<
" sgpgp=" << sgpgp
1059 <<
" cisped=" << cisped
1060 <<
" cisampl=" << cisampl <<
endmsg;
1068 for (
int isamp = 0; isamp < nfit; ++isamp) {
1070 leakage =
pulse(xvec[isamp], tleak, yleak);
1072 yvec[isamp] = yvec0[isamp] - leakage;
1074 cischi2 = cischi2 + (
xd *
xd) / (eyvec[isamp] * eyvec[isamp]);
1077 <<
" yvec[" << isamp <<
"]=" <<
yvec[isamp]
1078 <<
" leakage=" << leakage
1082 cischi2 = cischi2 / (nfit - 3.0);
1088 if ((cischi2 < leakchi2) && !fixedTime) {
1097 tempChi2 = leakchi2;
1115 for (
int isamp = 0; isamp < nfit; ++isamp) {
1118 int jsamp = (
int) xvec[isamp] - delta_peak;
1121 else if (jsamp >= nfit) jsamp = nfit - 1;
1131 err2 = eyvec[isamp] * eyvec[isamp];
1132 sy += yvec0[isamp] / err2;
1134 syg += yvec0[isamp] * gval / err2;
1135 sgg += gval * gval / err2;
1139 dgg0 = sg * sg - serr * sgg;
1140 if (fabs(dgg0) >
EPS_DG) {
1141 fixampl = (
sy * sg - serr * syg) / dgg0;
1142 fixped = (syg * sg -
sy * sgg) / dgg0;
1144 <<
" fixampl = " << fixampl
1145 <<
" fixped = " << fixped );
1150 <<
" small dgg0 = " << dgg0
1151 <<
", fixampl = " << fixampl
1152 <<
" fixped = " << fixped );
1158 for (
int isamp = 0; isamp < nfit; ++isamp) {
1159 dc = yvec0[isamp] -
scaledPulse(xvec[isamp], tpulse, ypulse);
1160 fixchi2 = fixchi2 + (dc * dc) / (eyvec[isamp] * eyvec[isamp]);
1162 <<
" yvec0[" << isamp <<
"]= " << yvec0[isamp]
1163 <<
" eyvec[" << isamp <<
"]= " << eyvec[isamp]
1164 <<
" fixchi2= " << fixchi2 );
1166 fixchi2 = fixchi2 / (nfit - 2.0);
1167 ATH_MSG_VERBOSE (
" fixchi2/(nfit-2.0)=" << fixchi2 <<
" nfit=" << nfit );
1182 tempChi2 = oldchi2 = -fabs(fixchi2);
1195 for (
int isamp = 0; isamp < nfit; ++isamp) {
1199 int jsamp = (
int) xvec[isamp] - delta_peak;
1202 else if (jsamp >= nfit) jsamp = nfit - 1;
1214 gpval =
derivative(xvec[isamp], tdpulse, dpulse);
1217 err2 = eyvec[isamp] * eyvec[isamp];
1218 sy += yvec0[isamp] / err2;
1220 sgp += gpval / err2;
1221 syg += yvec0[isamp] * gval / err2;
1222 sygp += yvec0[isamp] * gpval / err2;
1223 sgg += gval * gval / err2;
1224 sggp += gval * gpval / err2;
1225 sgpgp += gpval * gpval / err2;
1231 <<
" gpval=" << gpval
1232 <<
" sgp=" << sgp );
1235 dgg = sgg - sg * sg / serr;
1236 dggp = sggp - sg * sgp / serr;
1237 dgpgp = sgpgp - sgp * sgp / serr;
1238 dyg = syg -
sy * sg / serr;
1239 dygp = sygp -
sy * sgp / serr;
1240 dg = dgg * dgpgp - dggp * dggp;
1245 ampl = (dyg * dgpgp - dygp * dggp) /
dg;
1247 atau = (dyg * dggp - dygp * dgg) /
dg;
1249 ped = (
sy - ((dyg * dgpgp - dygp * dggp) * sg + (dyg * dggp - dygp * dgg) * sgp) /
dg)
1252 if (fabs(ampl) >
EPS_DG) {
1272 <<
" sgp=" << sgp <<
endmsg;
1276 <<
" sgg=" << sgg <<
endmsg;
1279 <<
" sgpgp=" << sgpgp <<
endmsg;
1284 <<
" dgpgp=" << dgpgp
1285 <<
" dyg*dgpgp=" << (dyg * dgpgp) <<
endmsg;
1289 <<
" dygp*dggp=" << (dygp * dggp) <<
endmsg;
1293 <<
" dyg*dggp=" << (dyg * dggp) <<
endmsg;
1297 <<
" dygp*dgg=" << (dygp * dgg) <<
endmsg;
1301 <<
" tau=" << tau <<
endmsg;
1310 for (
int isamp = 0; isamp < nfit; ++isamp) {
1311 dc = yvec0[isamp] -
scaledPulse(xvec[isamp], tpulse, ypulse);
1312 tempChi2 = tempChi2 + (dc * dc) / (eyvec[isamp] * eyvec[isamp]);
1314 <<
" yvec0[" << isamp <<
"]=" << yvec0[isamp]
1315 <<
" eyvec[" << isamp <<
"]=" << eyvec[isamp]
1316 <<
" dc=" << dc <<
" chi2=" << tempChi2 );
1318 tempChi2 = tempChi2 / (nfit - 3.0);
1319 ATH_MSG_VERBOSE (
" chi2/(nfit-3.0)=" << tempChi2 <<
" nfit=" << nfit );
1352 if (tempChi2 <
chi2) {
1373 <<
" chi2=" << tempChi2 );
1396 <<
" Ped=" << pedestal
1397 <<
" Amplitude=" << amplitude
1398 <<
" Chi2=" <<
chi2);
1421 ATH_MSG_VERBOSE (
"Lo gain leakage xvec[" << isamp <<
"]=" << xvec[isamp] );
1423 leakage =
pulse(xvec[isamp], tleak, yleak);
1424 dleakage =
derivative(xvec[isamp], tdleak, dleak);
1427 yvec[isamp] = yvec0[isamp] - pedg;
1430 <<
" yvec0[" << isamp <<
"]=" << yvec0[isamp]
1431 <<
" pedg=" << pedg );
1433 sllp += leakage * dleakage;
1434 sylp +=
yvec[isamp] * dleakage;
1435 slplp += dleakage * dleakage;
1438 if (fabs(slplp) >
EPS_DG && !fixedTime) {
1439 leaktau = (sllp - sylp) / slplp;
1450 <<
" slplp=" << slplp
1451 <<
" leaktau=" << leaktau );
1463 leakage =
pulse(xvec[isamp], tleak, yleak);
1465 yvec[isamp] = yvec0[isamp] - leakage;
1468 <<
" yvec0[" << isamp <<
"]=" << yvec0[isamp]
1469 <<
" leakage=" << leakage );
1471 err2 = eyvec[isamp] * eyvec[isamp];
1472 sy +=
yvec[isamp] / err2;
1474 syg +=
yvec[isamp] * gval / err2;
1475 sgg += gval * gval / err2;
1478 dgg0 = sg * sg - serr * sgg;
1479 if (fabs(dgg0) >
EPS_DG) {
1480 leakampl = (
sy * sg - serr * syg) / dgg0;
1481 leakped = (syg * sg -
sy * sgg) / dgg0;
1484 leakped =
sy / serr;
1499 leakage =
pulse(xvec[isamp], tleak, yleak);
1500 xd = yvec0[isamp] - (gval + leakage);
1501 leakchi2 = leakchi2 + (
xd *
xd) / (eyvec[isamp] * eyvec[isamp]);
1504 <<
" yvec0[" << isamp <<
"]=" << yvec0[isamp]
1506 <<
", leakage=" << leakage
1509 leakchi2 = leakchi2 / (nfit_real - 3.0);
1512 <<
" leakped=" << leakped
1513 <<
" leakampl=" << leakampl
1514 <<
" leakchi2=" << leakchi2 );
1523 leakage =
pulse(xvec[isamp], tleak, yleak);
1526 yvec[isamp] = yvec0[isamp] - leakage;
1529 <<
" yvec0[" << isamp <<
"]=" << yvec0[isamp]
1530 <<
" leakage=" << leakage
1531 <<
" yvec[" << isamp <<
"]=" <<
yvec[isamp] );
1548 gpval =
derivative(xvec[isamp], tdpulse, dpulse);
1549 err2 = eyvec[isamp] * eyvec[isamp];
1550 sy +=
yvec[isamp] / err2;
1552 sgp += gpval / err2;
1553 syg +=
yvec[isamp] * gval / err2;
1554 sygp +=
yvec[isamp] * gpval / err2;
1555 sgg += gval * gval / err2;
1556 sggp += gval * gpval / err2;
1557 sgpgp += gpval * gpval / err2;
1560 dgg = sgg - sg * sg / serr;
1561 dggp = sggp - sg * sgp / serr;
1562 dgpgp = sgpgp - sgp * sgp / serr;
1563 dyg = syg -
sy * sg / serr;
1564 dygp = sygp -
sy * sgp / serr;
1565 dg = dgg * dgpgp - dggp * dggp;
1568 cisampl = (dyg * dgpgp - dygp * dggp) /
dg;
1569 cisatau = (dyg * dggp - dygp * dgg) /
dg;
1571 - (dyg * dgpgp * sg - dygp * dggp * sg + dyg * dggp * sgp - dygp * dgg * sgp) /
dg)
1574 if (fabs(cisampl) >
EPS_DG) {
1575 cistau = cisatau / cisampl;
1597 <<
" sgpgp=" << sgpgp <<
endmsg;
1601 <<
" sgpgp=" << sgpgp
1607 <<
" cisped=" << cisped
1608 <<
" cisampl=" << cisampl <<
endmsg;
1621 leakage =
pulse(xvec[isamp], tleak, yleak);
1623 yvec[isamp] = yvec0[isamp] - leakage;
1625 cischi2 = cischi2 + (
xd *
xd) / (eyvec[isamp] * eyvec[isamp]);
1628 <<
" yvec[" << isamp <<
"]=" <<
yvec[isamp]
1629 <<
" leakage=" << leakage
1633 cischi2 = cischi2 / (nfit_real - 3.0);
1638 if ((cischi2 < leakchi2) && !fixedTime) {
1647 tempChi2 = leakchi2;
1673 gpval =
derivative(xvec[isamp], tdpulse, dpulse);
1675 err2 = eyvec[isamp] * eyvec[isamp];
1676 sy += yvec0[isamp] / err2;
1678 sgp += gpval / err2;
1679 syg += yvec0[isamp] * gval / err2;
1680 sygp += yvec0[isamp] * gpval / err2;
1681 sgg += gval * gval / err2;
1682 sggp += gval * gpval / err2;
1683 sgpgp += gpval * gpval / err2;
1689 <<
" gpval=" << gpval
1690 <<
" sgp=" << sgp );
1693 dgg = sgg - sg * sg / serr;
1694 dggp = sggp - sg * sgp / serr;
1695 dgpgp = sgpgp - sgp * sgp / serr;
1696 dyg = syg -
sy * sg / serr;
1697 dygp = sygp -
sy * sgp / serr;
1698 dg = dgg * dgpgp - dggp * dggp;
1703 ampl = (dyg * dgpgp - dygp * dggp) /
dg;
1705 atau = (dyg * dggp - dygp * dgg) /
dg;
1707 ped = (
sy - ((dyg * dgpgp - dygp * dggp) * sg + (dyg * dggp - dygp * dgg) * sgp) /
dg)
1710 if (fabs(ampl) >
EPS_DG) {
1730 <<
" sgp=" << sgp <<
endmsg;
1734 <<
" sgg=" << sgg <<
endmsg;
1737 <<
" sgpgp=" << sgpgp <<
endmsg;
1741 <<
" dgpgp=" << dgpgp
1742 <<
" dyg*dgpgp=" << (dyg * dgpgp) <<
endmsg;
1746 <<
" dygp*dggp=" << (dygp * dggp) <<
endmsg;
1750 <<
" dyg*dggp=" << (dyg * dggp) <<
endmsg;
1754 <<
" dygp*dgg=" << (dygp * dgg) <<
endmsg;
1758 <<
" tau=" << tau <<
endmsg;
1771 dc = yvec0[isamp] -
scaledPulse(xvec[isamp], tpulse, ypulse);
1772 tempChi2 = tempChi2 + (dc * dc) / (eyvec[isamp] * eyvec[isamp]);
1774 <<
" yvec0[" << isamp <<
"]=" << yvec0[isamp]
1775 <<
" eyvec[" << isamp <<
"]=" << eyvec[isamp]
1777 <<
" chi2=" << tempChi2 );
1779 tempChi2 = tempChi2 / (nfit_real - 3.0);
1781 <<
" nfit_real=" << nfit_real );
1811 <<
" Ped=" << pedestal
1812 <<
" Amplitude=" << amplitude
1813 <<
" Chi2=" <<
chi2 );