Calculate energy, time and chi2 for one channel using fitted pulse shape.
Method uses HFITV.
604 {
605
606 amplitude = 0.0;
608 pedestal = 0.0;
610
611 const HWIdentifier adcId =
digit->adc_HWID();
615
618
619
621 int noise_channel = (
ros < 3) ? channel :
channel + 48;
622
623 if (igain == 0) {
625 case 3:
627 if (rms > 0.0) break;
628
629 case 2:
631 if (rms > 0.0) break;
632
633 case 1:
635 if (rms > 0.0) break;
636
637 default:
639 }
640 } else if (igain == 1) {
642 case 3:
644 if (rms > 0.0) break;
645
646 case 2:
648 if (rms > 0.0) break;
649
650 case 1:
652 if (rms > 0.0) break;
653
654 default:
656 }
657 } else {
658
660 return;
661 }
662
664 << " ROS=" << ros
665 << " Channel=" << channel
666 << " gain=" << igain
667 << " RMS=" << rms
669
670
671 int ifit1 = 0;
673
676
678 << " ifit1=" << ifit1
679 << " ifit2=" << ifit1 + nfit - 1
680 << " nfit=" << nfit
685
686 std::vector<float> samples =
digit->samples();
687 double maxsamp = 0.0;
689
692
693 bool no_signal = true;
694 pedestal = samples[0];
695 const double delta = 1.e-6;
696
697 for (int isamp = 0; isamp < nfit; ++isamp) {
698 int j = isamp + ifit1;
699 xvec[isamp] = j;
700 yvec0[isamp] = samples[j];
701 if (no_signal) no_signal = (fabs(samples[j] - pedestal) < delta);
702
703
706 } else {
709 ATH_MSG_VERBOSE (
"Saturated ADC value yvec0[" << isamp <<
"]=" << yvec0[isamp]
711 } else {
713 ATH_MSG_VERBOSE (
"Zero ADC value yvec0[" << isamp <<
"]=" << yvec0[isamp]
714 << " RMS=" << eyvec[isamp] );
715 }
716 }
717
718 if (yvec0[isamp] > maxsamp) {
719
720
721 maxsamp = yvec0[isamp];
722 ipeakMax = j;
723 }
724 if (yvec0[isamp] < minsamp) {
725 minsamp = yvec0[isamp];
726 ipeakMin = j;
727 }
729 << ", xvec=" << xvec[isamp]
730 << ", yvec0=" << yvec0[isamp]
731 << ", eyvec=" << eyvec[isamp] );
732 }
733
734 if (no_signal) {
736 return;
737 }
738
739
740 double pedg = (yvec0[0] + yvec0[nfit - 1]) / 2.0;
741
742
743 int delta_peak = 0;
744
746
748
749 if (!fixedTime) {
756 } else {
757 fixedTime = true;
759 }
760 }
761
763 << " ipeakMax=" << ipeakMax
764 << " ipeakMin=" << ipeakMin
765 << " fixedTime=" << ((fixedTime) ? "true" : "false") );
766
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;
775
777 if (igain == 0) {
787 } else {
796 }
797 } else {
807 } else {
816 }
817 }
818 } else {
820 if (igain == 0) {
825 } else {
830 }
831 } else {
832 if (igain == 0) {
837 } else {
842 }
843 }
844 }
845
846
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;
854
855
856 int niter = 0;
857 do {
858 ++niter;
860
861 if (tempChi2 < oldchi2) oldchi2 = tempChi2;
862
863
864
866
868
870
871
874
875 sllp = 0.0;
876 sylp = 0.0;
877 slplp = 0.0;
878 for (int isamp = 0; isamp < nfit; ++isamp) {
879 ATH_MSG_VERBOSE (
"Lo gain leakage xvec[" << isamp <<
"]=" << xvec[isamp] );
880
881 leakage =
pulse(xvec[isamp], tleak, yleak);
882 dleakage =
derivative(xvec[isamp], tdleak, dleak);
883
884
885 yvec[isamp] = yvec0[isamp] - pedg;
886
888 << " yvec0[" << isamp << "]=" << yvec0[isamp]
889 << " pedg=" << pedg );
890
891 sllp += leakage * dleakage;
892 sylp +=
yvec[isamp] * dleakage;
893 slplp += dleakage * dleakage;
894 }
895
896 if (fabs(slplp) >
EPS_DG && !fixedTime) {
897 leaktau = (sllp - sylp) / slplp;
898
902 } else {
903 leaktau = 0.0;
904 }
905
907 << " sylp=" << sylp
908 << " slplp=" << slplp
909 << " leaktau=" << leaktau );
910
911
912
915 sg = 0.0;
916 syg = 0.0;
917 sgg = 0.0;
918 serr = 0.0;
919 for (int isamp = 0; isamp < nfit; ++isamp) {
920 leakage =
pulse(xvec[isamp], tleak, yleak);
922 yvec[isamp] = yvec0[isamp] - leakage;
923
925 << " yvec0[" << isamp << "]=" << yvec0[isamp]
926 << " leakage=" << leakage );
927
928 err2 = eyvec[isamp] * eyvec[isamp];
930 sg += gval / err2;
931 syg +=
yvec[isamp] * gval / err2;
932 sgg += gval * gval / err2;
933 serr += 1.0 / err2;
934 }
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;
939 } else {
940 leakampl = 0.0;
942 }
943
944
946
947 leakchi2 = 0.0;
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]);
956
958 << " yvec0[" << isamp << "]=" << yvec0[isamp]
959 << " gval=" << gval
960 << ", leakage=" << leakage
961 << ", xd=" << xd );
962 }
963 leakchi2 = leakchi2 / (nfit - 3.0);
964
966 << " leakped=" << leakped
967 << " leakampl=" << leakampl
968 << " leakchi2=" << leakchi2 );
969
970
971 if (!fixedTime) {
975 for (int isamp = 0; isamp < nfit; ++isamp) {
976 leakage =
pulse(xvec[isamp], tleak, yleak);
977
978
979 yvec[isamp] = yvec0[isamp] - leakage;
980
982 << " yvec0[" << isamp << "]=" << yvec0[isamp]
983 << " leakage=" << leakage
984 << " yvec[" << isamp << "]=" << yvec[isamp] );
985 }
990 sg = 0.0;
991 sgp = 0.0;
992 syg = 0.0;
993 sygp = 0.0;
994 sgg = 0.0;
995 sggp = 0.0;
996 sgpgp = 0.0;
997 serr = 0.0;
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;
1003 sg += gval / 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;
1010 serr += 1.0 / err2;
1011 }
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;
1018
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)
1024 / serr;
1025
1026 if (fabs(cisampl) >
EPS_DG) {
1027 cistau = cisatau / cisampl;
1031 } else {
1032 cistau = 0.0;
1033 }
1034 } else {
1035 cisampl = 0.0;
1036 cisatau = 0.0;
1037 cistau = 0.0;
1039 }
1040
1041 if (
msgLvl(MSG::VERBOSE)) {
1042 msg(MSG::VERBOSE) <<
" sy=" <<
sy
1043 << " sg=" << sg
1044 << " sgp=" << sgp
1045 << " syg=" << syg
1046 << " sygp=" << sygp
1047 << " sgg=" << sgg
1048 << " sggp=" << sggp
1049 <<
" sgpgp=" << sgpgp <<
endmsg;
1050
1051 msg(MSG::VERBOSE) <<
" dgg=" << dgg
1052 << " dggp=" << dggp
1053 << " sgpgp=" << sgpgp
1054 << " dyg=" << dyg
1055 << " dygp=" << dygp
1057
1058 msg(MSG::VERBOSE) <<
" cistau=" << cistau
1059 << " cisped=" << cisped
1060 <<
" cisampl=" << cisampl <<
endmsg;
1061 }
1062
1063
1064 cischi2 = 0.0;
1068 for (int isamp = 0; isamp < nfit; ++isamp) {
1070 leakage =
pulse(xvec[isamp], tleak, yleak);
1071
1072 yvec[isamp] = yvec0[isamp] - leakage;
1074 cischi2 = cischi2 + (
xd *
xd) / (eyvec[isamp] * eyvec[isamp]);
1075
1077 << " yvec[" << isamp << "]=" << yvec[isamp]
1078 << " leakage=" << leakage
1079 << " gval=" << gval
1080 << " xd=" << xd );
1081 }
1082 cischi2 = cischi2 / (nfit - 3.0);
1083
1085 }
1086
1087
1088 if ((cischi2 < leakchi2) && !fixedTime) {
1089 tau = cistau;
1091 ampl = cisampl;
1092 tempChi2 = cischi2;
1093 } else {
1094 tau = leaktau;
1096 ampl = leakampl;
1097 tempChi2 = leakchi2;
1098 }
1099 }
1100
1101
1102 else {
1103
1104 if (niter == 1) {
1105
1108
1110 sg = 0.0;
1111 sgg = 0.0;
1112 syg = 0.0;
1113 serr = 0.0;
1114
1115 for (int isamp = 0; isamp < nfit; ++isamp) {
1117
1118 int jsamp = (
int) xvec[isamp] - delta_peak;
1119 if (jsamp < 0)
1120 jsamp = 0;
1121 else if (jsamp >= nfit) jsamp = nfit - 1;
1122
1123 if (igain == 0) {
1125 } else {
1127 }
1128 } else {
1130 }
1131 err2 = eyvec[isamp] * eyvec[isamp];
1132 sy += yvec0[isamp] / err2;
1133 sg += gval / err2;
1134 syg += yvec0[isamp] * gval / err2;
1135 sgg += gval * gval / err2;
1136 serr += 1.0 / err2;
1137 }
1138 fixtau = 0.0;
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 );
1146 } else {
1147 fixampl = 0.0;
1150 << " small dgg0 = " << dgg0
1151 << ", fixampl = " << fixampl
1152 << " fixped = " << fixped );
1153 }
1157 fixchi2 = 0.0;
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 );
1165 }
1166 fixchi2 = fixchi2 / (nfit - 2.0);
1167 ATH_MSG_VERBOSE (
" fixchi2/(nfit-2.0)=" << fixchi2 <<
" nfit=" << nfit );
1168
1170
1171
1175 }
1176
1177 if (fixedTime) {
1179 tau = fixtau;
1181 ampl = fixampl;
1182 tempChi2 = oldchi2 = -fabs(fixchi2);
1183 } else {
1184
1186 sg = 0.0;
1187 sgp = 0.0;
1188 syg = 0.0;
1189 sygp = 0.0;
1190 sgg = 0.0;
1191 sggp = 0.0;
1192 sgpgp = 0.0;
1193 serr = 0.0;
1194
1195 for (int isamp = 0; isamp < nfit; ++isamp) {
1197
1198
1199 int jsamp = (
int) xvec[isamp] - delta_peak;
1200 if (jsamp < 0)
1201 jsamp = 0;
1202 else if (jsamp >= nfit) jsamp = nfit - 1;
1203
1204 if (igain == 0) {
1207 } else {
1210 }
1211 } else {
1212
1214 gpval =
derivative(xvec[isamp], tdpulse, dpulse);
1215 }
1216
1217 err2 = eyvec[isamp] * eyvec[isamp];
1218 sy += yvec0[isamp] / err2;
1219 sg += gval / 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;
1226 serr += 1.0 / err2;
1227
1229 << " gval=" << gval
1230 << " sg=" << sg
1231 << " gpval=" << gpval
1232 << " sgp=" << sgp );
1233 }
1234
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;
1241
1242
1244
1245 ampl = (dyg * dgpgp - dygp * dggp) / dg;
1246
1247 atau = (dyg * dggp - dygp * dgg) / dg;
1248
1249 ped = (
sy - ((dyg * dgpgp - dygp * dggp) * sg + (dyg * dggp - dygp * dgg) * sgp) /
dg)
1250 / serr;
1251
1252 if (fabs(ampl) >
EPS_DG) {
1253
1254 tau = atau / ampl;
1258 } else {
1259 tau = 0.0;
1260 }
1261 } else {
1262 ampl = 0.0;
1263 atau = 0.0;
1264 tau = 0.0;
1266 }
1267
1268 if (
msgLvl(MSG::VERBOSE)) {
1270 msg(MSG::VERBOSE) <<
" sy=" <<
sy
1271 << " sg=" << sg
1272 <<
" sgp=" << sgp <<
endmsg;
1273
1274 msg(MSG::VERBOSE) <<
" syg=" << syg
1275 << " sygp=" << sygp
1276 <<
" sgg=" << sgg <<
endmsg;
1277
1278 msg(MSG::VERBOSE) <<
" sggp=" << sggp
1279 <<
" sgpgp=" << sgpgp <<
endmsg;
1280
1281 msg(MSG::VERBOSE) <<
" ampl = (dyg*dgpgp - dygp*dggp)= " << ampl <<
endmsg;
1282
1283 msg(MSG::VERBOSE) <<
" dyg=" << dyg
1284 << " dgpgp=" << dgpgp
1285 <<
" dyg*dgpgp=" << (dyg * dgpgp) <<
endmsg;
1286
1287 msg(MSG::VERBOSE) <<
" dygp=" << dygp
1288 << " dggp=" << dggp
1289 <<
" dygp*dggp=" << (dygp * dggp) <<
endmsg;
1290
1291 msg(MSG::VERBOSE) <<
" dyg=" << dyg
1292 << " dggp=" << dggp
1293 <<
" dyg*dggp=" << (dyg * dggp) <<
endmsg;
1294
1295 msg(MSG::VERBOSE) <<
" dygp=" << dygp
1296 << " dgg=" << dgg
1297 <<
" dygp*dgg=" << (dygp * dgg) <<
endmsg;
1298
1299 msg(MSG::VERBOSE) <<
" dg=" <<
dg
1300 << " atau=" << atau
1301 <<
" tau=" << tau <<
endmsg;
1302 }
1303
1307
1308 tempChi2 = 0;
1309
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 );
1317 }
1318 tempChi2 = tempChi2 / (nfit - 3.0);
1319 ATH_MSG_VERBOSE (
" chi2/(nfit-3.0)=" << tempChi2 <<
" nfit=" << nfit );
1320 }
1321 }
1322
1323 if (
msgLvl(MSG::VERBOSE))
1324 msg(MSG::VERBOSE) <<
" t0fit: " <<
m_t0Fit << ((tau < 0.0) ?
" - " :
" + ") << fabs(tau);
1325
1326
1335 } else {
1336
1337
1340
1341
1348 }
1349 }
1350
1351
1352 if (tempChi2 <
chi2) {
1355 amplitude = ampl;
1357 }
1358
1362 else
1364
1366 }
1367
1370 << " phase=" << tau
1371 << " ped=" << ped
1372 << " ampl=" << ampl
1373 << " chi2=" << tempChi2 );
1374
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1393
1395 << " Time=" << time
1396 << " Ped=" << pedestal
1397 << " Amplitude=" << amplitude
1398 <<
" Chi2=" <<
chi2);
1399
1401 int nfit_real;
1402
1403
1404
1406
1408
1410
1411
1413 if (!fixedTime) {
1415
1416 sllp = 0.0;
1417 sylp = 0.0;
1418 slplp = 0.0;
1419 for (int isamp = 0;
1421 ATH_MSG_VERBOSE (
"Lo gain leakage xvec[" << isamp <<
"]=" << xvec[isamp] );
1422
1423 leakage =
pulse(xvec[isamp], tleak, yleak);
1424 dleakage =
derivative(xvec[isamp], tdleak, dleak);
1425
1426
1427 yvec[isamp] = yvec0[isamp] - pedg;
1428
1430 << " yvec0[" << isamp << "]=" << yvec0[isamp]
1431 << " pedg=" << pedg );
1432
1433 sllp += leakage * dleakage;
1434 sylp +=
yvec[isamp] * dleakage;
1435 slplp += dleakage * dleakage;
1436 }
1437
1438 if (fabs(slplp) >
EPS_DG && !fixedTime) {
1439 leaktau = (sllp - sylp) / slplp;
1440
1444 } else {
1445 leaktau = 0.0;
1446 }
1447
1449 << " sylp=" << sylp
1450 << " slplp=" << slplp
1451 << " leaktau=" << leaktau );
1452
1453
1454
1457 sg = 0.0;
1458 syg = 0.0;
1459 sgg = 0.0;
1460 serr = 0.0;
1461 for (int isamp = 0;
1463 leakage =
pulse(xvec[isamp], tleak, yleak);
1465 yvec[isamp] = yvec0[isamp] - leakage;
1466
1468 << " yvec0[" << isamp << "]=" << yvec0[isamp]
1469 << " leakage=" << leakage );
1470
1471 err2 = eyvec[isamp] * eyvec[isamp];
1472 sy +=
yvec[isamp] / err2;
1473 sg += gval / err2;
1474 syg +=
yvec[isamp] * gval / err2;
1475 sgg += gval * gval / err2;
1476 serr += 1.0 / err2;
1477 }
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;
1482 } else {
1483 leakampl = 0.0;
1484 leakped =
sy / serr;
1485 }
1486
1487
1489
1490 leakchi2 = 0.0;
1491 nfit_real = 0;
1495 for (int isamp = 0;
1497 ++nfit_real;
1499 leakage =
pulse(xvec[isamp], tleak, yleak);
1500 xd = yvec0[isamp] - (gval + leakage);
1501 leakchi2 = leakchi2 + (
xd *
xd) / (eyvec[isamp] * eyvec[isamp]);
1502
1504 << " yvec0[" << isamp << "]=" << yvec0[isamp]
1505 << " gval=" << gval
1506 << ", leakage=" << leakage
1507 << ", xd=" << xd );
1508 }
1509 leakchi2 = leakchi2 / (nfit_real - 3.0);
1510
1512 << " leakped=" << leakped
1513 << " leakampl=" << leakampl
1514 << " leakchi2=" << leakchi2 );
1515
1516
1520 for (int isamp = 0;
1522
1523 leakage =
pulse(xvec[isamp], tleak, yleak);
1524
1525
1526 yvec[isamp] = yvec0[isamp] - leakage;
1527
1529 << " yvec0[" << isamp << "]=" << yvec0[isamp]
1530 << " leakage=" << leakage
1531 << " yvec[" << isamp << "]=" << yvec[isamp] );
1532 }
1537 sg = 0.0;
1538 sgp = 0.0;
1539 syg = 0.0;
1540 sygp = 0.0;
1541 sgg = 0.0;
1542 sggp = 0.0;
1543 sgpgp = 0.0;
1544 serr = 0.0;
1545 for (int isamp = 0;
1548 gpval =
derivative(xvec[isamp], tdpulse, dpulse);
1549 err2 = eyvec[isamp] * eyvec[isamp];
1550 sy +=
yvec[isamp] / err2;
1551 sg += gval / 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;
1558 serr += 1.0 / err2;
1559 }
1560 if (serr == 0) serr = 1;
1561 dgg = sgg - sg * sg / serr;
1562 dggp = sggp - sg * sgp / serr;
1563 dgpgp = sgpgp - sgp * sgp / serr;
1564 dyg = syg -
sy * sg / serr;
1565 dygp = sygp -
sy * sgp / serr;
1566 dg = dgg * dgpgp - dggp * dggp;
1567
1569 cisampl = (dyg * dgpgp - dygp * dggp) / dg;
1570 cisatau = (dyg * dggp - dygp * dgg) / dg;
1572 - (dyg * dgpgp * sg - dygp * dggp * sg + dyg * dggp * sgp - dygp * dgg * sgp) /
dg)
1573 / serr;
1574
1575 if (fabs(cisampl) >
EPS_DG) {
1576 cistau = cisatau / cisampl;
1580 } else {
1581 cistau = 0.0;
1582 }
1583 } else {
1584 cisampl = 0.0;
1585 cisatau = 0.0;
1586 cistau = 0.0;
1588 }
1589
1590 if (
msgLvl(MSG::VERBOSE)) {
1591 msg(MSG::VERBOSE) <<
" sy=" <<
sy
1592 << " sg=" << sg
1593 << " sgp=" << sgp
1594 << " syg=" << syg
1595 << " sygp=" << sygp
1596 << " sgg=" << sgg
1597 << " sggp=" << sggp
1598 <<
" sgpgp=" << sgpgp <<
endmsg;
1599
1600 msg(MSG::VERBOSE) <<
" dgg=" << dgg
1601 << " dggp=" << dggp
1602 << " sgpgp=" << sgpgp
1603 << " dyg=" << dyg
1604 << " dygp=" << dygp
1606
1607 msg(MSG::VERBOSE) <<
" cistau=" << cistau
1608 << " cisped=" << cisped
1609 <<
" cisampl=" << cisampl <<
endmsg;
1610 }
1611
1612
1613 cischi2 = 0.0;
1614 nfit_real = 0;
1618 for (int isamp = 0;
1620 ++nfit_real;
1622 leakage =
pulse(xvec[isamp], tleak, yleak);
1623
1624 yvec[isamp] = yvec0[isamp] - leakage;
1626 cischi2 = cischi2 + (
xd *
xd) / (eyvec[isamp] * eyvec[isamp]);
1627
1629 << " yvec[" << isamp << "]=" << yvec[isamp]
1630 << " leakage=" << leakage
1631 << " gval=" << gval
1632 << " xd=" << xd );
1633 }
1634 cischi2 = cischi2 / (nfit_real - 3.0);
1635
1637
1638
1639 if ((cischi2 < leakchi2) && !fixedTime) {
1640 tau = cistau;
1642 ampl = cisampl;
1643 tempChi2 = cischi2;
1644 } else {
1645 tau = leaktau;
1647 ampl = leakampl;
1648 tempChi2 = leakchi2;
1649 }
1650
1651 }
1652 } else {
1653 if (!fixedTime) {
1654
1655
1659
1661 sg = 0.0;
1662 sgp = 0.0;
1663 syg = 0.0;
1664 sygp = 0.0;
1665 sgg = 0.0;
1666 sggp = 0.0;
1667 sgpgp = 0.0;
1668 serr = 0.0;
1669
1670 for (int isamp = 0;
1672
1674 gpval =
derivative(xvec[isamp], tdpulse, dpulse);
1675
1676 err2 = eyvec[isamp] * eyvec[isamp];
1677 sy += yvec0[isamp] / err2;
1678 sg += gval / err2;
1679 sgp += gpval / err2;
1680 syg += yvec0[isamp] * gval / err2;
1681 sygp += yvec0[isamp] * gpval / err2;
1682 sgg += gval * gval / err2;
1683 sggp += gval * gpval / err2;
1684 sgpgp += gpval * gpval / err2;
1685 serr += 1.0 / err2;
1686
1688 << " gval=" << gval
1689 << " sg=" << sg
1690 << " gpval=" << gpval
1691 << " sgp=" << sgp );
1692 }
1693
1694 dgg = sgg - sg * sg / serr;
1695 dggp = sggp - sg * sgp / serr;
1696 dgpgp = sgpgp - sgp * sgp / serr;
1697 dyg = syg -
sy * sg / serr;
1698 dygp = sygp -
sy * sgp / serr;
1699 dg = dgg * dgpgp - dggp * dggp;
1700
1701
1703
1704 ampl = (dyg * dgpgp - dygp * dggp) / dg;
1705
1706 atau = (dyg * dggp - dygp * dgg) / dg;
1707
1708 ped = (
sy - ((dyg * dgpgp - dygp * dggp) * sg + (dyg * dggp - dygp * dgg) * sgp) /
dg)
1709 / serr;
1710
1711 if (fabs(ampl) >
EPS_DG) {
1712
1713 tau = atau / ampl;
1717 } else {
1718 tau = 0.0;
1719 }
1720 } else {
1721 ampl = 0.0;
1722 atau = 0.0;
1723 tau = 0.0;
1725 }
1726
1727 if (
msgLvl(MSG::VERBOSE)) {
1729 msg(MSG::VERBOSE) <<
" sy=" <<
sy
1730 << " sg=" << sg
1731 <<
" sgp=" << sgp <<
endmsg;
1732
1733 msg(MSG::VERBOSE) <<
" syg=" << syg
1734 << " sygp=" << sygp
1735 <<
" sgg=" << sgg <<
endmsg;
1736
1737 msg(MSG::VERBOSE) <<
" sggp=" << sggp
1738 <<
" sgpgp=" << sgpgp <<
endmsg;
1739
1740 msg(MSG::VERBOSE) <<
" ampl = (dyg*dgpgp - dygp*dggp)= " << ampl <<
endmsg;
1741 msg(MSG::VERBOSE) <<
" dyg=" << dyg
1742 << " dgpgp=" << dgpgp
1743 <<
" dyg*dgpgp=" << (dyg * dgpgp) <<
endmsg;
1744
1745 msg(MSG::VERBOSE) <<
" dygp=" << dygp
1746 << " dggp=" << dggp
1747 <<
" dygp*dggp=" << (dygp * dggp) <<
endmsg;
1748
1749 msg(MSG::VERBOSE) <<
" dyg=" << dyg
1750 << " dggp=" << dggp
1751 <<
" dyg*dggp=" << (dyg * dggp) <<
endmsg;
1752
1753 msg(MSG::VERBOSE) <<
" dygp=" << dygp
1754 << " dgg=" << dgg
1755 <<
" dygp*dgg=" << (dygp * dgg) <<
endmsg;
1756
1757 msg(MSG::VERBOSE) <<
" dg=" <<
dg
1758 << " atau=" << atau
1759 <<
" tau=" << tau <<
endmsg;
1760 }
1761
1765
1766 tempChi2 = 0;
1767 nfit_real = 0;
1768
1769 for (int isamp = 0;
1771 ++nfit_real;
1772 dc = yvec0[isamp] -
scaledPulse(xvec[isamp], tpulse, ypulse);
1773 tempChi2 = tempChi2 + (dc * dc) / (eyvec[isamp] * eyvec[isamp]);
1775 << " yvec0[" << isamp << "]=" << yvec0[isamp]
1776 << " eyvec[" << isamp << "]=" << eyvec[isamp]
1777 << " dc=" << dc
1778 << " chi2=" << tempChi2 );
1779 }
1780 tempChi2 = tempChi2 / (nfit_real - 3.0);
1782 << " nfit_real=" << nfit_real );
1783 }
1784 }
1785
1786 if (
msgLvl(MSG::VERBOSE))
1787 msg(MSG::VERBOSE) <<
" t0fit: " <<
m_t0Fit << ((tau < 0.0) ?
" - " :
" + ") << fabs(tau);
1788
1791
1792
1799 }
1800
1804 amplitude = ampl;
1806 }
1807
1808 }
1809
1811 << " Time=" << time
1812 << " Ped=" << pedestal
1813 << " Amplitude=" << amplitude
1814 <<
" Chi2=" <<
chi2 );
1815}
#define SATURATED_ADC_VALUE
std::vector< double > m_dummy
double derivative(double x, const std::vector< double > *xvec, const std::vector< double > *yvec) const
double chi2(TH1 *h0, TH1 *h1)
time(flags, cells_name, *args, **kw)