Calculate energy, time and chi2 for one channel using fitted pulse shape.
Method uses HFITV.
286 {
287
288 amplitude = 0.0;
290 pedestal = 0.0;
292
293 const HWIdentifier adcId =
digit->adc_HWID();
297
300
302
303
305 int noise_channel = (
ros < 3) ? channel :
channel + 48;
306
307 if (igain == 0) {
309 case 3:
311 if (rms > 0.0) break;
312
313 case 2:
315 if (rms > 0.0) break;
316
317 case 1:
319 if (rms > 0.0) break;
320
321 default:
323 }
324 } else if (igain == 1) {
326 case 3:
328 if (rms > 0.0) break;
329
330 case 2:
332 if (rms > 0.0) break;
333
334 case 1:
336 if (rms > 0.0) break;
337
338 default:
340 }
341 } else {
342
344 return;
345 }
346
348 << " ROS=" << ros
349 << " Channel=" << channel
350 << " gain=" << igain
351 << " RMS=" << rms
353
354
355 int ifit1 = 0;
357
360
362 << " ifit1=" << ifit1
363 << " ifit2=" << ifit1 + nfit - 1
364 << " nfit=" << nfit
369 << " demoCh=" << ((demo)?"true":"false") );
370
371 std::vector<float> samples =
digit->samples();
372 samples.erase(samples.begin(),samples.begin()+
m_firstSample);
374 double maxsamp = 0.0;
376
379
380 bool no_signal = true;
381 pedestal = samples[0];
382 const double delta = 1.e-6;
383
384 for (int isamp = 0; isamp < nfit; ++isamp) {
385 int j = isamp + ifit1;
386 xvec[isamp] = j;
387 yvec0[isamp] = samples[j];
388 if (no_signal) no_signal = (fabs(samples[j] - pedestal) < delta);
389
390
393 } else {
396 ATH_MSG_VERBOSE(
"Saturated ADC value yvec0[" << isamp <<
"]=" << yvec0[isamp]
398
399 } else {
401 ATH_MSG_VERBOSE(
"Zero ADC value yvec0[" << isamp <<
"]=" << yvec0[isamp]
402 << " RMS=" << eyvec[isamp] );
403 }
404 }
405
406 if (yvec0[isamp] > maxsamp) {
407
408
409 maxsamp = yvec0[isamp];
410 ipeakMax = j;
411 }
412 if (yvec0[isamp] < minsamp) {
413 minsamp = yvec0[isamp];
414 ipeakMin = j;
415 }
417 << ", xvec=" << xvec[isamp]
418 << ", yvec0=" << yvec0[isamp]
419 << ", eyvec=" << eyvec[isamp] );
420 }
421
422 if (no_signal) {
424 return;
425 }
426
427
428 double pedg = (yvec0[0] + yvec0[nfit - 1]) / 2.0;
429
430
431 int delta_peak = 0;
432
434
436
437 if (!fixedTime) {
439
440 pedg = yvec0[0];
444 } else {
445 fixedTime = true;
447 }
448 } else {
455 } else {
456 fixedTime = true;
458 }
459 }
460 }
461
462 double expectedTime = 0.;
465 delta_peak = std::round(expectedTime /
DTIME);
467 }
468
471 << " ipeakMax=" << ipeakMax
472 << " ipeakMin=" << ipeakMin
473 << " fixedTime=" << ((fixedTime) ? "true" : "false") );
474
475 const std::vector<double>* tpulse = &
m_dummy;
476 const std::vector<double>* ypulse = &
m_dummy;
477 const std::vector<double>* tdpulse = &
m_dummy;
478 const std::vector<double>* dpulse = &
m_dummy;
479 const std::vector<double>* tleak = &
m_dummy;
480 const std::vector<double>* yleak = &
m_dummy;
481 const std::vector<double>* tdleak = &
m_dummy;
482 const std::vector<double>* dleak = &
m_dummy;
483
484
487 if (demo) {
490 }
491
493 if (igain == 0) {
503 } else {
512 }
513 } else {
523 } else {
532 }
533 }
534 } else {
535 if (dolas) {
536 if (igain == 0) {
541 } else {
546 }
547 } else {
548 if (igain == 0) {
553 } else {
558 }
559 }
560 }
561
562 if (demo) {
563 docis = false;
565 }
566
567
568 double gval, gpval,
sy, syg, sygp, sg, sgp, sgg, sgpgp, sggp, serr, err2;
569 double dgg0, dgg, dggp, dgpgp, dyg, dygp,
dg, dc,
xd;
570 double sllp, sylp, slplp, dleakage, leakage;
571 double fixtau = 0.0, fixped = 0.0, fixampl = 0.0, fixchi2 =
MAX_CHI2;
572 double leaktau = 0.0, leakped = 0.0, leakampl = 0.0, leakchi2 =
MAX_CHI2;
573 double cistau = 0.0, cisped = 0.0, cisampl = 0.0, cisatau = 0.0, cischi2 =
MAX_CHI2;
575
576
577 int niter = 0;
578 do {
579 ++niter;
581
582 if (tempChi2 < oldchi2) oldchi2 = tempChi2;
583
584
585
587
589
591
592
595
596 sllp = 0.0;
597 sylp = 0.0;
598 slplp = 0.0;
599 for (int isamp = 0; isamp < nfit; ++isamp) {
600 ATH_MSG_VERBOSE (
"Lo gain leakage xvec[" << isamp <<
"]=" << xvec[isamp] );
601
602 leakage =
pulse(xvec[isamp], tleak, yleak);
603 dleakage =
derivative(xvec[isamp], tdleak, dleak);
604
605
606 yvec[isamp] = yvec0[isamp] - pedg;
607
609 << " yvec0[" << isamp << "]=" << yvec0[isamp]
610 << " pedg=" << pedg);
611
612 sllp += leakage * dleakage;
613 sylp +=
yvec[isamp] * dleakage;
614 slplp += dleakage * dleakage;
615 }
616
617 if (fabs(slplp) >
EPS_DG && !fixedTime) {
618 leaktau = (sllp - sylp) / slplp;
619
623 } else {
624 leaktau = 0.0;
625 }
626
628 << " sylp=" << sylp
629 << " slplp=" << slplp
630 << " leaktau=" << leaktau);
631
632
633
636 sg = 0.0;
637 syg = 0.0;
638 sgg = 0.0;
639 serr = 0.0;
640 for (int isamp = 0; isamp < nfit; ++isamp) {
641 leakage =
pulse(xvec[isamp], tleak, yleak);
643 yvec[isamp] = yvec0[isamp] - leakage;
644
646 << " yvec0[" << isamp << "]=" << yvec0[isamp]
647 << " leakage=" << leakage );
648
649 err2 = eyvec[isamp] * eyvec[isamp];
651 sg += gval / err2;
652 syg +=
yvec[isamp] * gval / err2;
653 sgg += gval * gval / err2;
654 serr += 1.0 / err2;
655 }
656
657 dgg0 = sg * sg - serr * sgg;
658 if (fabs(dgg0) >
EPS_DG) {
659 leakampl = (
sy * sg - serr * syg) / dgg0;
660 leakped = (syg * sg -
sy * sgg) / dgg0;
661 } else {
662 leakampl = 0.0;
664 }
665
666
668
669 leakchi2 = 0.0;
673
674 for (int isamp = 0; isamp < nfit; ++isamp) {
676 leakage =
pulse(xvec[isamp], tleak, yleak);
677 xd = yvec0[isamp] - (gval + leakage);
678 leakchi2 = leakchi2 + (
xd *
xd) / (eyvec[isamp] * eyvec[isamp]);
679
681 << " yvec0[" << isamp << "]=" << yvec0[isamp]
682 << " gval=" << gval
683 << ", leakage=" << leakage
684 << ", xd=" << xd );
685 }
686
687 leakchi2 = leakchi2/(nfit-3.0);
688
690 << " leakped=" << leakped
691 << " leakampl=" << leakampl
692 << " leakchi2=" << leakchi2 );
693
694
695 if (!fixedTime) {
699
700 for (int isamp = 0; isamp < nfit; ++isamp) {
701 leakage =
pulse(xvec[isamp], tleak, yleak);
702
703
704 yvec[isamp] = yvec0[isamp] - leakage;
705
707 << " yvec0[" << isamp << "]=" << yvec0[isamp]
708 << " leakage=" << leakage
709 << " yvec[" << isamp << "]=" << yvec[isamp] );
710 }
711
716 sg = 0.0;
717 sgp = 0.0;
718 syg = 0.0;
719 sygp = 0.0;
720 sgg = 0.0;
721 sggp = 0.0;
722 sgpgp = 0.0;
723 serr = 0.0;
724
725 for (int isamp = 0; isamp < nfit; ++isamp) {
727 gpval =
derivative(xvec[isamp], tdpulse, dpulse);
728 err2 = eyvec[isamp] * eyvec[isamp];
730 sg += gval / err2;
731 sgp += gpval / err2;
732 syg +=
yvec[isamp] * gval / err2;
733 sygp +=
yvec[isamp] * gpval / err2;
734 sgg += gval * gval / err2;
735 sggp += gval * gpval / err2;
736 sgpgp += gpval * gpval / err2;
737 serr += 1.0 / err2;
738 }
739
740 dgg = sgg - sg * sg / serr;
741 dggp = sggp - sg * sgp / serr;
742 dgpgp = sgpgp - sgp * sgp / serr;
743 dyg = syg -
sy * sg / serr;
744 dygp = sygp -
sy * sgp / serr;
745 dg = dgg * dgpgp - dggp * dggp;
746
748 cisampl = (dyg * dgpgp - dygp * dggp) / dg;
749 cisatau = (dyg * dggp - dygp * dgg) / dg;
751 - (dyg * dgpgp * sg - dygp * dggp * sg + dyg * dggp * sgp - dygp * dgg * sgp) /
dg)
752 / serr;
753
754 if (fabs(cisampl) >
EPS_DG) {
755 cistau = cisatau / cisampl;
759 } else {
760 cistau = 0.0;
761 }
762 } else {
763 cisampl = 0.0;
764 cisatau = 0.0;
765 cistau = 0.0;
767 }
768
769 if (
msgLvl(MSG::VERBOSE)) {
770 msg(MSG::VERBOSE) <<
" sy=" <<
sy
771 << " sg=" << sg
772 << " sgp=" << sgp
773 << " syg=" << syg
774 << " sygp=" << sygp
775 << " sgg=" << sgg
776 << " sggp=" << sggp
777 <<
" sgpgp=" << sgpgp <<
endmsg;
778
779 msg(MSG::VERBOSE) <<
" dgg=" << dgg
780 << " dggp=" << dggp
781 << " sgpgp=" << sgpgp
782 << " dyg=" << dyg
783 << " dygp=" << dygp
785
786 msg(MSG::VERBOSE) <<
" cistau=" << cistau
787 << " cisped=" << cisped
788 <<
" cisampl=" << cisampl <<
endmsg;
789 }
790
791
792 cischi2 = 0.0;
796
797 for (int isamp = 0; isamp < nfit; ++isamp) {
799 leakage =
pulse(xvec[isamp], tleak, yleak);
800
801 yvec[isamp] = yvec0[isamp] - leakage;
803 cischi2 = cischi2 + (
xd *
xd) / (eyvec[isamp] * eyvec[isamp]);
804
806 << " yvec[" << isamp << "]=" << yvec[isamp]
807 << " leakage=" << leakage
808 << " gval=" << gval
809 << " xd=" << xd );
810 }
811
812 cischi2 = cischi2 / (nfit - 3.0);
813
815 }
816
817
818 if ((cischi2 < leakchi2) && !fixedTime) {
819 tau = cistau;
821 ampl = cisampl;
822 tempChi2 = cischi2;
823 } else {
824 tau = leaktau;
826 ampl = leakampl;
827 tempChi2 = leakchi2;
828 }
829
830 } else {
831
832
833 if (niter == 1) {
834
835
838
840 sg = 0.0;
841 sgg = 0.0;
842 syg = 0.0;
843 serr = 0.0;
844
845 for (int isamp = 0; isamp < nfit; ++isamp) {
846 if (!dolas) {
847
848 int jsamp = (
int) xvec[isamp] - delta_peak;
849 if (jsamp < 0)
850 jsamp = 0;
851 else if (jsamp >= nfit) jsamp = nfit - 1;
852
853 if (igain == 0) {
855 } else {
857 }
858 } else {
860 }
861 err2 = eyvec[isamp] * eyvec[isamp];
862 sy += yvec0[isamp] / err2;
863 sg += gval / err2;
864 syg += yvec0[isamp] * gval / err2;
865 sgg += gval * gval / err2;
866 serr += 1.0 / err2;
867 }
868
869 fixtau = 0.0;
870 dgg0 = sg * sg - serr * sgg;
871 if (fabs(dgg0) >
EPS_DG) {
872 fixampl = (
sy * sg - serr * syg) / dgg0;
873 fixped = (syg * sg -
sy * sgg) / dgg0;
875 << " fixped = " << fixped );
876 } else {
877 fixampl = 0.0;
880 << ", fixampl = " << fixampl
881 << " fixped = " << fixped );
882 }
883
887 fixchi2 = 0.0;
888 for (int isamp = 0; isamp < nfit; ++isamp) {
889 dc = yvec0[isamp] -
scaledPulse(xvec[isamp], tpulse, ypulse);
890 fixchi2 = fixchi2 + (dc * dc) / (eyvec[isamp] * eyvec[isamp]);
892 << " yvec0[" << isamp << "]= " << yvec0[isamp]
893 << " eyvec[" << isamp << "]= " << eyvec[isamp]
894 << " fixchi2= " << fixchi2 );
895 }
896 fixchi2 = fixchi2 / (nfit - 2.0);
898 << " nfit=" << nfit );
899
901
902
906 }
907
908 if (fixedTime) {
910 tau = fixtau;
912 ampl = fixampl;
913 tempChi2 = oldchi2 = -fabs(fixchi2);
914 } else {
915
917 sg = 0.0;
918 sgp = 0.0;
919 syg = 0.0;
920 sygp = 0.0;
921 sgg = 0.0;
922 sggp = 0.0;
923 sgpgp = 0.0;
924 serr = 0.0;
925
926 for (int isamp = 0; isamp < nfit; ++isamp) {
927 if ((niter == 1) && (!dolas)) {
928
929
930 int jsamp = (
int) xvec[isamp] - delta_peak;
931 if (jsamp < 0)
932 jsamp = 0;
933 else if (jsamp >= nfit) jsamp = nfit - 1;
934
935 if (igain == 0) {
938 } else {
941 }
942 } else {
943
945 gpval =
derivative(xvec[isamp], tdpulse, dpulse);
946 }
947
948 err2 = eyvec[isamp] * eyvec[isamp];
949 sy += yvec0[isamp] / err2;
950 sg += gval / err2;
951 sgp += gpval / err2;
952 syg += yvec0[isamp] * gval / err2;
953 sygp += yvec0[isamp] * gpval / err2;
954 sgg += gval * gval / err2;
955 sggp += gval * gpval / err2;
956 sgpgp += gpval * gpval / err2;
957 serr += 1.0 / err2;
958
960 << " gval=" << gval
961 << " sg=" << sg
962 << " gpval=" << gpval
963 << " sgp=" << sgp );
964 }
965
966 dgg = sgg - sg * sg / serr;
967 dggp = sggp - sg * sgp / serr;
968 dgpgp = sgpgp - sgp * sgp / serr;
969 dyg = syg -
sy * sg / serr;
970 dygp = sygp -
sy * sgp / serr;
971 dg = dgg * dgpgp - dggp * dggp;
972
973
975
976 ampl = (dyg * dgpgp - dygp * dggp) / dg;
977
978 atau = (dyg * dggp - dygp * dgg) / dg;
979
980 ped = (
sy - ((dyg * dgpgp - dygp * dggp) * sg + (dyg * dggp - dygp * dgg) * sgp) /
dg)
981 / serr;
982
983 if (fabs(ampl) >
EPS_DG) {
984
985 tau = atau / ampl;
989 } else {
990 tau = 0.0;
991 }
992 } else {
993 ampl = 0.0;
994 atau = 0.0;
995 tau = 0.0;
997 }
998
999 if (
msgLvl(MSG::VERBOSE)) {
1001 msg(MSG::VERBOSE) <<
" sy=" <<
sy
1002 << " sg=" << sg
1003 <<
" sgp=" << sgp <<
endmsg;
1004
1005 msg(MSG::VERBOSE) <<
" syg=" << syg
1006 << " sygp=" << sygp
1007 <<
" sgg=" << sgg <<
endmsg;
1008
1009 msg(MSG::VERBOSE) <<
" sggp=" << sggp
1010 <<
" sgpgp=" << sgpgp <<
endmsg;
1011
1012 msg(MSG::VERBOSE) <<
" ampl = (dyg*dgpgp - dygp*dggp)= " << ampl <<
endmsg;
1013
1014 msg(MSG::VERBOSE) <<
" dyg=" << dyg
1015 << " dgpgp=" << dgpgp
1016 <<
" dyg*dgpgp=" << (dyg * dgpgp) <<
endmsg;
1017
1018 msg(MSG::VERBOSE) <<
" dygp=" << dygp
1019 << " dggp=" << dggp
1020 <<
" dygp*dggp=" << (dygp * dggp) <<
endmsg;
1021
1022 msg(MSG::VERBOSE) <<
" dyg=" << dyg
1023 << " dggp=" << dggp
1024 <<
" dyg*dggp=" << (dyg * dggp) <<
endmsg;
1025
1026 msg(MSG::VERBOSE) <<
" dygp=" << dygp
1027 << " dgg=" << dgg
1028 <<
" dygp*dgg=" << (dygp * dgg) <<
endmsg;
1029
1030 msg(MSG::VERBOSE) <<
" dg=" <<
dg
1031 << " atau=" << atau
1032 <<
" tau=" << tau <<
endmsg;
1033 }
1034
1038
1039 tempChi2 = 0;
1040
1041 for (int isamp = 0; isamp < nfit; ++isamp) {
1042 dc = yvec0[isamp] -
scaledPulse(xvec[isamp], tpulse, ypulse);
1043 tempChi2 = tempChi2 + (dc * dc) / (eyvec[isamp] * eyvec[isamp]);
1045 << " yvec0[" << isamp << "]=" << yvec0[isamp]
1046 << " eyvec[" << isamp << "]=" << eyvec[isamp]
1047 << " dc=" << dc
1048 << " chi2=" << tempChi2 );
1049 }
1050
1051 tempChi2 = tempChi2 / (nfit - 3.0);
1053 << " nfit=" << nfit );
1054 }
1055 }
1056
1057 if (
msgLvl(MSG::VERBOSE))
1058 msg(MSG::VERBOSE) <<
" t0fit: " <<
m_t0Fit << ((tau < 0.0) ?
" - " :
" + ") << fabs(tau);
1059
1060
1063 if (
msgLvl(MSG::VERBOSE))
1070 } else {
1071
1072
1075
1076
1083 }
1084 }
1085
1086
1087 if (tempChi2 <
chi2) {
1090 amplitude = ampl;
1092 }
1093
1097 else
1099
1101 }
1102
1105 << " phase=" << tau
1106 << " ped=" << ped
1107 << " ampl=" << ampl
1108 << " chi2=" << tempChi2 );
1109
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1128
1130 << " Time=" << time
1131 << " Ped=" << pedestal
1132 << " Amplitude=" << amplitude
1133 <<
" Chi2=" <<
chi2 );
1134
1136 int nfit_real;
1137
1138
1139
1141
1143
1145
1146
1148 if (!fixedTime) {
1150
1151 sllp = 0.0;
1152 sylp = 0.0;
1153 slplp = 0.0;
1155
1156 ATH_MSG_VERBOSE (
"Lo gain leakage xvec[" << isamp <<
"]=" << xvec[isamp] );
1157
1158 leakage =
pulse(xvec[isamp], tleak, yleak);
1159 dleakage =
derivative(xvec[isamp], tdleak, dleak);
1160
1161
1162 yvec[isamp] = yvec0[isamp] - pedg;
1163
1165 << " yvec0[" << isamp << "]=" << yvec0[isamp]
1166 << " pedg=" << pedg );
1167
1168 sllp += leakage * dleakage;
1169 sylp +=
yvec[isamp] * dleakage;
1170 slplp += dleakage * dleakage;
1171 }
1172
1173 if (fabs(slplp) >
EPS_DG && !fixedTime) {
1174 leaktau = (sllp - sylp) / slplp;
1175
1179 } else {
1180 leaktau = 0.0;
1181 }
1182
1184 << " sylp=" << sylp
1185 << " slplp=" << slplp
1186 << " leaktau=" << leaktau );
1187
1188
1189
1192 sg = 0.0;
1193 syg = 0.0;
1194 sgg = 0.0;
1195 serr = 0.0;
1196
1197 for (int isamp = 0;
1199
1200 leakage =
pulse(xvec[isamp], tleak, yleak);
1202 yvec[isamp] = yvec0[isamp] - leakage;
1203
1205 << " yvec0[" << isamp << "]=" << yvec0[isamp]
1206 << " leakage=" << leakage );
1207
1208 err2 = eyvec[isamp] * eyvec[isamp];
1209 sy +=
yvec[isamp] / err2;
1210 sg += gval / err2;
1211 syg +=
yvec[isamp] * gval / err2;
1212 sgg += gval * gval / err2;
1213 serr += 1.0 / err2;
1214 }
1215
1216 dgg0 = sg * sg - serr * sgg;
1217 if (fabs(dgg0) >
EPS_DG) {
1218 leakampl = (
sy * sg - serr * syg) / dgg0;
1219 leakped = (syg * sg -
sy * sgg) / dgg0;
1220 } else {
1221 leakampl = 0.0;
1222 leakped =
sy / serr;
1223 }
1224
1225
1227
1228 leakchi2 = 0.0;
1229 nfit_real = 0;
1233
1234 for (int isamp = 0;
1236
1237 ++nfit_real;
1239 leakage =
pulse(xvec[isamp], tleak, yleak);
1240 xd = yvec0[isamp] - (gval + leakage);
1241 leakchi2 = leakchi2 + (
xd *
xd) / (eyvec[isamp] * eyvec[isamp]);
1242
1244 << " yvec0[" << isamp << "]=" << yvec0[isamp]
1245 << " gval=" << gval
1246 << ", leakage=" << leakage
1247 << ", xd=" << xd );
1248 }
1249
1250 leakchi2 = leakchi2 / (nfit_real - 3.0);
1251
1253 << " leakped=" << leakped
1254 << " leakampl=" << leakampl
1255 << " leakchi2=" << leakchi2 );
1256
1257
1261
1262 for (int isamp = 0;
1264
1265 leakage =
pulse(xvec[isamp], tleak, yleak);
1266
1267
1268 yvec[isamp] = yvec0[isamp] - leakage;
1269
1271 << " yvec0[" << isamp << "]=" << yvec0[isamp]
1272 << " leakage=" << leakage
1273 << " yvec[" << isamp << "]=" << yvec[isamp] );
1274 }
1275
1280 sg = 0.0;
1281 sgp = 0.0;
1282 syg = 0.0;
1283 sygp = 0.0;
1284 sgg = 0.0;
1285 sggp = 0.0;
1286 sgpgp = 0.0;
1287 serr = 0.0;
1288
1289 for (int isamp = 0;
1291
1293 gpval =
derivative(xvec[isamp], tdpulse, dpulse);
1294 err2 = eyvec[isamp] * eyvec[isamp];
1295 sy +=
yvec[isamp] / err2;
1296 sg += gval / err2;
1297 sgp += gpval / err2;
1298 syg +=
yvec[isamp] * gval / err2;
1299 sygp +=
yvec[isamp] * gpval / err2;
1300 sgg += gval * gval / err2;
1301 sggp += gval * gpval / err2;
1302 sgpgp += gpval * gpval / err2;
1303 serr += 1.0 / err2;
1304 }
1305
1306 if (serr == 0) serr = 1;
1307 dgg = sgg - sg * sg / serr;
1308 dggp = sggp - sg * sgp / serr;
1309 dgpgp = sgpgp - sgp * sgp / serr;
1310 dyg = syg -
sy * sg / serr;
1311 dygp = sygp -
sy * sgp / serr;
1312 dg = dgg * dgpgp - dggp * dggp;
1313
1315 cisampl = (dyg * dgpgp - dygp * dggp) / dg;
1316 cisatau = (dyg * dggp - dygp * dgg) / dg;
1318 - (dyg * dgpgp * sg - dygp * dggp * sg + dyg * dggp * sgp - dygp * dgg * sgp) /
dg)
1319 / serr;
1320
1321 if (fabs(cisampl) >
EPS_DG) {
1322 cistau = cisatau / cisampl;
1326 } else {
1327 cistau = 0.0;
1328 }
1329 } else {
1330 cisampl = 0.0;
1331 cisatau = 0.0;
1332 cistau = 0.0;
1334 }
1335
1336 if (
msgLvl(MSG::VERBOSE)) {
1337 msg(MSG::VERBOSE) <<
" sy=" <<
sy
1338 << " sg=" << sg
1339 << " sgp=" << sgp
1340 << " syg=" << syg
1341 << " sygp=" << sygp
1342 << " sgg=" << sgg
1343 << " sggp=" << sggp
1344 <<
" sgpgp=" << sgpgp <<
endmsg;
1345
1346 msg(MSG::VERBOSE) <<
" dgg=" << dgg
1347 << " dggp=" << dggp
1348 << " sgpgp=" << sgpgp
1349 << " dyg=" << dyg
1350 << " dygp=" << dygp
1352
1353 msg(MSG::VERBOSE) <<
" cistau=" << cistau
1354 << " cisped=" << cisped
1355 <<
" cisampl=" << cisampl <<
endmsg;
1356 }
1357
1358
1359 cischi2 = 0.0;
1360 nfit_real = 0;
1364
1365 for (int isamp = 0;
1367
1368 ++nfit_real;
1370 leakage =
pulse(xvec[isamp], tleak, yleak);
1371
1372 yvec[isamp] = yvec0[isamp] - leakage;
1374 cischi2 = cischi2 + (
xd *
xd) / (eyvec[isamp] * eyvec[isamp]);
1375
1377 << " yvec[" << isamp << "]=" << yvec[isamp]
1378 << " leakage=" << leakage
1379 << " gval=" << gval
1380 << " xd=" << xd );
1381 }
1382
1383 cischi2=cischi2/(nfit_real-3.0);
1384
1386
1387
1388 if ((cischi2 < leakchi2) && !fixedTime) {
1389 tau = cistau;
1391 ampl = cisampl;
1392 tempChi2 = cischi2;
1393 } else {
1394 tau = leaktau;
1396 ampl = leakampl;
1397 tempChi2 = leakchi2;
1398 }
1399
1400 }
1401 } else {
1402 if (!fixedTime) {
1403
1404
1408
1410 sg = 0.0;
1411 sgp = 0.0;
1412 syg = 0.0;
1413 sygp = 0.0;
1414 sgg = 0.0;
1415 sggp = 0.0;
1416 sgpgp = 0.0;
1417 serr = 0.0;
1418
1419 for (int isamp = 0;
1421
1423 gpval =
derivative(xvec[isamp], tdpulse, dpulse);
1424
1425 err2 = eyvec[isamp] * eyvec[isamp];
1426 sy += yvec0[isamp] / err2;
1427 sg += gval / err2;
1428 sgp += gpval / err2;
1429 syg += yvec0[isamp] * gval / err2;
1430 sygp += yvec0[isamp] * gpval / err2;
1431 sgg += gval * gval / err2;
1432 sggp += gval * gpval / err2;
1433 sgpgp += gpval * gpval / err2;
1434 serr += 1.0 / err2;
1435
1437 << " gval=" << gval
1438 << " sg=" << sg
1439 << " gpval=" << gpval
1440 << " sgp=" << sgp );
1441 }
1442
1443 dgg = sgg - sg * sg / serr;
1444 dggp = sggp - sg * sgp / serr;
1445 dgpgp = sgpgp - sgp * sgp / serr;
1446 dyg = syg -
sy * sg / serr;
1447 dygp = sygp -
sy * sgp / serr;
1448 dg = dgg * dgpgp - dggp * dggp;
1449
1450
1452
1453 ampl = (dyg * dgpgp - dygp * dggp) / dg;
1454
1455 atau = (dyg * dggp - dygp * dgg) / dg;
1456
1457 ped = (
sy - ((dyg * dgpgp - dygp * dggp) * sg + (dyg * dggp - dygp * dgg) * sgp) /
dg)
1458 / serr;
1459
1460 if (fabs(ampl) >
EPS_DG) {
1461
1462 tau = atau / ampl;
1466 } else {
1467 tau = 0.0;
1468 }
1469 } else {
1470 ampl = 0.0;
1471 atau = 0.0;
1472 tau = 0.0;
1474 }
1475
1476 if (
msgLvl(MSG::VERBOSE)) {
1478 msg(MSG::VERBOSE) <<
" sy=" <<
sy
1479 << " sg=" << sg
1480 <<
" sgp=" << sgp <<
endmsg;
1481
1482 msg(MSG::VERBOSE) <<
" syg=" << syg
1483 << " sygp=" << sygp
1484 <<
" sgg=" << sgg <<
endmsg;
1485
1486 msg(MSG::VERBOSE) <<
" sggp=" << sggp
1487 <<
" sgpgp=" << sgpgp <<
endmsg;
1488
1489 msg(MSG::VERBOSE) <<
" ampl = (dyg*dgpgp - dygp*dggp)= " << ampl <<
endmsg;
1490 msg(MSG::VERBOSE) <<
" dyg=" << dyg
1491 << " dgpgp=" << dgpgp
1492 <<
" dyg*dgpgp=" << (dyg * dgpgp) <<
endmsg;
1493
1494 msg(MSG::VERBOSE) <<
" dygp=" << dygp
1495 << " dggp=" << dggp
1496 <<
" dygp*dggp=" << (dygp * dggp) <<
endmsg;
1497
1498 msg(MSG::VERBOSE) <<
" dyg=" << dyg
1499 << " dggp=" << dggp
1500 <<
" dyg*dggp=" << (dyg * dggp) <<
endmsg;
1501
1502 msg(MSG::VERBOSE) <<
" dygp=" << dygp
1503 << " dgg=" << dgg
1504 <<
" dygp*dgg=" << (dygp * dgg) <<
endmsg;
1505
1506 msg(MSG::VERBOSE) <<
" dg=" <<
dg
1507 << " atau=" << atau
1508 <<
" tau=" << tau <<
endmsg;
1509 }
1510
1514
1515 tempChi2 = 0;
1516 nfit_real = 0;
1517
1518 for (int isamp = 0;
1520
1521 ++nfit_real;
1522 dc = yvec0[isamp] -
scaledPulse(xvec[isamp], tpulse, ypulse);
1523 tempChi2 = tempChi2 + (dc * dc) / (eyvec[isamp] * eyvec[isamp]);
1525 << " yvec0[" << isamp << "]=" << yvec0[isamp]
1526 << " eyvec[" << isamp << "]=" << eyvec[isamp]
1527 << " dc=" << dc
1528 << " chi2=" << tempChi2 );
1529 }
1530
1531 tempChi2 = tempChi2 / (nfit_real - 3.0);
1533 << " nfit_real=" << nfit_real );
1534 }
1535 }
1536
1537 if (
msgLvl(MSG::VERBOSE))
1538 msg(MSG::VERBOSE) <<
" t0fit: " <<
m_t0Fit << ((tau < 0.0) ?
" - " :
" + ") << fabs(tau);
1539
1540
1542 if (
msgLvl(MSG::VERBOSE) )
1544
1545
1552 }
1553
1557 amplitude = ampl;
1559 }
1560
1561 }
1562
1564 << " Ped=" << pedestal
1565 << " Amplitude=" << amplitude
1566 <<
" Chi2=" <<
chi2 );
1567}
static unsigned int getDrawerIdx(unsigned int ros, unsigned int drawer)
Returns a drawer hash.
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)