121void TileCorrelation::sum(std::vector<double> &digits,
int ros,
int drawer,
int channel,
int gain,
int &dignum) {
126 dignum = digits.size();
128 m_N[ros][drawer][channel][gain]++;
129 N_d = double(
m_N[ros][drawer][channel][gain]);
130 for (
int i = 0; i < dignum; ++i) {
131 m_S[ros][drawer][channel][gain][i] += digits[i];
132 for (
int j = 0; j < dignum; ++j)
133 m_SS[ros][drawer][channel][gain][i][j] += digits[i] * digits[j];
138 <<
" drawer=" << drawer
139 <<
" channel=" << channel
141 <<
" N=" <<
m_N[ros][drawer][channel][gain]
142 <<
" Sum[1]=" <<
m_S[ros][drawer][channel][gain][1]
143 <<
" Sum[2]=" <<
m_S[ros][drawer][channel][gain][2]
144 <<
" Sum[1][2]=" <<
m_SS[ros][drawer][channel][gain][1][2]
145 <<
" Sum[1][1]=" <<
m_SS[ros][drawer][channel][gain][1][1]
146 <<
" Sum[2][2]=" <<
m_SS[ros][drawer][channel][gain][2][2]
147 <<
" B[1][2]=" <<
m_SS[ros][drawer][channel][gain][1][2] / N_d
148 -
m_S[ros][drawer][channel][gain][1] / N_d *
m_S[ros][drawer][channel][gain][2] / N_d
149 <<
" Correlation[1][2]=" << (N_d *
m_SS[ros][drawer][channel][gain][1][2]
150 -
m_S[ros][drawer][channel][gain][1] *
m_S[ros][drawer][channel][gain][2])
151 / sqrt((N_d *
m_SS[ros][drawer][channel][gain][1][1]
152 -
m_S[ros][drawer][channel][gain][1]
153 *
m_S[ros][drawer][channel][gain][1])
154 * (N_d *
m_SS[ros][drawer][channel][gain][2][2]
155 -
m_S[ros][drawer][channel][gain][2]
156 *
m_S[ros][drawer][channel][gain][2]))
213 int &dignum,
int chthres) {
217 dignum = digits.size();
221 m_N[ros][drawer][channel][gain]++;
225 if (ros == 1 && drawer == 1 && channel == 0 && gain == 1)
229 for (
int i = 0; i < dignum -
m_lag; ++i) {
230 m_S1[ros][drawer][channel][gain][
m_lag - 1] += digits[i];
232 m_S12[ros][drawer][channel][gain][
m_lag - 1] += digits[i] * digits[i +
m_lag];
233 m_S11[ros][drawer][channel][gain][
m_lag - 1] += digits[i] * digits[i];
237 if (
m_lag == 1 && ros == 1 && drawer == 1 && channel == 0 && gain == 1)
240 <<
" S1=" <<
m_S1[ros][drawer][channel][gain][
m_lag - 1]);
246 -
m_S1[ros][drawer][channel][gain][
m_lag - 1] *
m_S2[ros][drawer][channel][gain][
m_lag - 1])
249 -
m_S1[ros][drawer][channel][gain][
m_lag - 1] *
m_S1[ros][drawer][channel][gain][
m_lag - 1])
251 -
m_S2[ros][drawer][channel][gain][
m_lag - 1] *
m_S2[ros][drawer][channel][gain][
m_lag - 1]));
253 if (
m_lag == 1 && ros == 1 && drawer == 1 && channel == 0 && gain == 1)
255 <<
" corr_sum=" <<
m_corrSum[ros][drawer][channel][gain][
m_lag - 1]
256 <<
" corr_sum_sq=" <<
m_corrSum2[ros][drawer][channel][gain][
m_lag - 1]);
264 if (
m_lag == 1 && ros == 1 && drawer == 1 && channel == 0 && gain == 1)
266 <<
" jentry-chthres=" <<
m_jEntry - chthres
267 <<
" lag=1, ros=1, drawer=1, channel=0, gain=1"
270 <<
" sum corr_mean=" <<
m_corrSum[ros][drawer][channel][gain][
m_lag - 1]
287 for (
int ros = 0; ros < 4; ++ros)
295 if (is7to9 && dignum == 7) {
296 for (
int i = 0; i < 9; ++i)
297 m_R[ros][drawer][channel][gain][i][i] = 1.;
300 for (
int i = 0; i < 9 -
m_lag; ++i) {
307 m_R[ros][drawer][channel][gain][i][i +
m_lag] = 0.;
308 m_R[ros][drawer][channel][gain][i +
m_lag][i] = 0.;
311 if (-1. >
m_R[ros][drawer][channel][gain][i][i +
m_lag]
312 ||
m_R[ros][drawer][channel][gain][i][i +
m_lag] > 1.)
313 m_R[ros][drawer][channel][gain][i][i +
m_lag] = 0.;
314 if (-1. >
m_R[ros][drawer][channel][gain][i +
m_lag][i]
315 ||
m_R[ros][drawer][channel][gain][i +
m_lag][i] > 1.)
316 m_R[ros][drawer][channel][gain][i +
m_lag][i] = 0.;
320 for (
int i = 0; i < dignum; i++)
321 m_R[ros][drawer][channel][gain][i][i] = 1.;
324 for (
int i = 0; i < dignum -
m_lag; ++i) {
329 if (-1. >
m_R[ros][drawer][channel][gain][i][i +
m_lag]
330 ||
m_R[ros][drawer][channel][gain][i][i +
m_lag] > 1.)
331 m_R[ros][drawer][channel][gain][i][i +
m_lag] = 0.;
332 if (-1. >
m_R[ros][drawer][channel][gain][i +
m_lag][i]
333 ||
m_R[ros][drawer][channel][gain][i +
m_lag][i] > 1.)
334 m_R[ros][drawer][channel][gain][i +
m_lag][i] = 0.;
369 const TileHWID *tileHWID,
int dignum) {
373 CLHEP::HepMatrix correlation(dignum, 1, 0);
375 std::fstream *correlationFile =
new std::fstream(correlationSummOptFilterFile.c_str(), std::fstream::out);
376 if (correlationFile->is_open())
ATH_MSG_INFO(correlationSummOptFilterFile <<
" file open");
378 if (deltaCorrelation) {
380 for (
int j = 0; j < dignum; ++j) {
382 if (
m_R[0][0][0][0][i][j] > -100000. &&
m_R[0][0][0][0][i][j] < 100000.)
383 correlation[i][j] =
m_R[0][0][0][0][i][j];
385 correlation[i][j] = 0.0;
388 *correlationFile << correlation.T() << std::endl;
390 for (
int ros = 0; ros < 4; ++ros)
392 int frag = tileHWID->
frag(ros + 1, drawer);
396 <<
" drawer " << drawer << MSG::hex
397 <<
" frag0x " << frag << MSG::dec
398 <<
" channel " << channel
400 <<
" N " <<
m_N[ros][drawer][channel][gain]);
402 if (
m_N[ros][drawer][channel][gain] > 0) {
404 for (
int j = 0; j < dignum; ++j) {
406 if (
m_R[ros][drawer][channel][gain][i][j] > -100000. &&
m_R[ros][drawer][channel][gain][i][j] < 100000.)
407 correlation[i][j] =
m_R[ros][drawer][channel][gain][i][j];
409 correlation[i][j] = 0.0;
412 *correlationFile <<
"ros " << ros
413 <<
" drawer " << drawer << std::hex
414 <<
" frag0x " << frag << std::dec
415 <<
" channel " << channel
417 <<
" N " <<
m_N[ros][drawer][channel][gain]
428 correlationFile->close();
433 const TileHWID *tileHWID,
int dignum) {
437 CLHEP::HepMatrix correlation(dignum, dignum, 0);
439 std::fstream *correlationFile =
new std::fstream(correlationMatrixOptFilterFile.c_str(), std::fstream::out);
440 if (correlationFile->is_open())
ATH_MSG_INFO(correlationMatrixOptFilterFile <<
" file open");
442 if (deltaCorrelation) {
443 for (
int i = 0; i < dignum; ++i)
444 for (
int j = 0; j < dignum; ++j) {
445 if (
m_R[0][0][0][0][i][j] > -100000. &&
m_R[0][0][0][0][i][j] < 100000.)
446 correlation[i][j] =
m_R[0][0][0][0][i][j];
448 correlation[i][j] = 0.0;
451 *correlationFile << correlation << std::endl;
453 for (
int ros = 0; ros < 4; ++ros)
455 int frag = tileHWID->
frag(ros + 1, drawer);
459 <<
" drawer " << drawer << MSG::hex
460 <<
" frag0x " << frag << MSG::dec
461 <<
" channel " << channel
463 <<
" N " <<
m_N[ros][drawer][channel][gain]);
466 if (
m_N[ros][drawer][channel][gain] > 0) {
467 for (
int i = 0; i < dignum; ++i)
468 for (
int j = 0; j < dignum; ++j) {
469 if (
m_R[ros][drawer][channel][gain][i][j] > -100000.
470 &&
m_R[ros][drawer][channel][gain][i][j] < 100000.)
471 correlation[i][j] =
m_R[ros][drawer][channel][gain][i][j];
473 correlation[i][j] = 0.0;
476 *correlationFile <<
"ros " << ros
477 <<
" drawer " << drawer << std::hex
478 <<
" frag0x " << frag << std::dec
479 <<
" channel " << channel
481 <<
" N " <<
m_N[ros][drawer][channel][gain]
482 << correlation << std::endl;
491 correlationFile->close();
496 (
bool deltaCorrelation,
497 const std::vector<double>& shapeFormLG,
498 const std::vector<double>& shapeFormHG,
499 const std::vector<double>& shapeFormDerivativeLG,
500 const std::vector<double>& shapeFormDerivativeHG,
501 const std::string& aiLoOptFilterFile,
502 const std::string& biLoOptFilterFile,
503 const std::string& aiHiOptFilterFile,
504 const std::string& biHiOptFilterFile,
505 const TileHWID *tileHWID,
int dignum) {
509 CLHEP::HepMatrix correlation(dignum, dignum, 0);
510 CLHEP::HepMatrix inverse(dignum, dignum, 0);
511 CLHEP::HepMatrix
zero(dignum, dignum, 0);
512 CLHEP::HepMatrix pulseShape(dignum, 1, 0);
513 CLHEP::HepMatrix pulseShapeDerivative(dignum, 1, 0);
514 CLHEP::HepMatrix
a(dignum, 1, 0);
515 CLHEP::HepMatrix b(dignum, 1, 0);
517 std::fstream *aiLoFile =
new std::fstream(aiLoOptFilterFile.c_str(), std::fstream::out);
518 std::fstream *biLoFile =
new std::fstream(biLoOptFilterFile.c_str(), std::fstream::out);
519 std::fstream *aiHiFile =
new std::fstream(aiHiOptFilterFile.c_str(), std::fstream::out);
520 std::fstream *biHiFile =
new std::fstream(biHiOptFilterFile.c_str(), std::fstream::out);
523 if (aiLoFile->is_open() && aiLoFile->is_open() && aiLoFile->is_open() && aiLoFile->is_open()) {
552 if (
msgLvl(MSG::VERBOSE)) {
553 msg(MSG::VERBOSE) <<
"shapeFormLG, shapeFormDerivativeLG, shapeFormHG, shapeFormDerivativeHG" <<
endmsg;
554 for (
int i = 0; i < int(shapeFormLG.size()); ++i)
555 msg(MSG::VERBOSE) << i <<
" " << std::setw(18) << std::setprecision(10)
556 << shapeFormLG[i] <<
" " << std::setw(18) << std::setprecision(10)
557 << shapeFormDerivativeLG[i] <<
" " << std::setw(18) << std::setprecision(10)
558 << shapeFormHG[i] <<
" " << std::setw(18) << std::setprecision(10)
559 << shapeFormDerivativeHG[i] <<
" " <<
endmsg;
570 double Q1, Q2, Q3, Delta;
573 if (deltaCorrelation) {
575 for (
int pha = -12; pha < 13; ++pha) {
579 for (
int i = 0; i < dignum; ++i)
580 for (
int j = 0; j < dignum; ++j)
581 correlation[i][j] =
m_R[0][0][0][0][i][j];
583 inverse = correlation.inverse(ierr);
585 for (
int i = 0; i < dignum; ++i) {
587 pulseShape[i][0] = shapeFormLG[i * 25 + 12 + pha];
588 pulseShapeDerivative[i][0] = shapeFormDerivativeLG[i * 25 + 12 + pha];
590 pulseShape[i][0] = shapeFormHG[i * 25 + 12 + pha];
591 pulseShapeDerivative[i][0] = shapeFormDerivativeHG[i * 25 + 12 + pha];
595 Q1 = ((pulseShape.T()) * inverse * pulseShape).determinant();
596 Q2 = ((pulseShapeDerivative.T()) * inverse * pulseShapeDerivative).determinant();
597 Q3 = ((pulseShapeDerivative.T()) * inverse * pulseShape).determinant();
598 Delta = Q1 * Q2 - Q3 * Q3;
600 a = Q2 / Delta * inverse * pulseShape - Q3 / Delta * inverse * pulseShapeDerivative;
601 b = Q3 / Delta * inverse * pulseShape - Q1 / Delta * inverse * pulseShapeDerivative;
604 *aiLoFile << std::setw(6) << pha;
605 for (
int i = 0; i < dignum; ++i)
606 *aiLoFile << std::setw(18) << std::setprecision(10) <<
a(i + 1, 1);
607 *aiLoFile << std::endl;
609 *biLoFile << std::setw(6) << pha;
610 for (
int i = 0; i < dignum; ++i)
611 *biLoFile << std::setw(18) << std::setprecision(10) << b(i + 1, 1);
612 *biLoFile << std::endl;
614 *aiHiFile << std::setw(6) << pha;
615 for (
int i = 0; i < dignum; ++i)
616 *aiHiFile << std::setw(18) << std::setprecision(10) <<
a(i + 1, 1);
617 *aiHiFile << std::endl;
619 *biHiFile << std::setw(6) << pha;
620 for (
int i = 0; i < dignum; ++i)
621 *biHiFile << std::setw(18) << std::setprecision(10) << b(i + 1, 1);
622 *biHiFile << std::endl;
627 for (
int ros = 0; ros < 4; ++ros)
629 int frag = tileHWID->
frag(ros + 1, drawer);
632 if (
m_N[ros][drawer][channel][gain] > 0) {
634 *aiLoFile <<
"ros " << ros
635 <<
" drawer " << drawer << std::hex
636 <<
" frag0x " << frag << std::dec
637 <<
" channel " << channel
638 <<
" N " <<
m_N[ros][drawer][channel][0] << std::endl;
640 *biLoFile <<
"ros " << ros
641 <<
" drawer " << drawer << std::hex
642 <<
" frag0x " << frag << std::dec
643 <<
" channel " << channel
644 <<
" N " <<
m_N[ros][drawer][channel][0] << std::endl;
647 *aiHiFile <<
"ros " << ros
648 <<
" drawer " << drawer << std::hex
649 <<
" frag0x " << frag << std::dec
650 <<
" channel " << channel
651 <<
" N " <<
m_N[ros][drawer][channel][1] << std::endl;
653 *biHiFile <<
"ros " << ros
654 <<
" drawer " << drawer << std::hex
655 <<
" frag0x " << frag << std::dec
656 <<
" channel " << channel
657 <<
" N " <<
m_N[ros][drawer][channel][1] << std::endl;
660 for (
int pha = -12; pha < 13; ++pha) {
664 for (
int i = 0; i < dignum; ++i)
665 for (
int j = 0; j < dignum; ++j)
666 correlation[i][j] =
m_R[ros][drawer][channel][gain][i][j];
668 inverse = correlation.inverse(ierr);
670 for (
int i = 0; i < dignum; ++i) {
672 pulseShape[i][0] = shapeFormLG[i * 25 + 12 + pha];
673 pulseShapeDerivative[i][0] = shapeFormDerivativeLG[i * 25 + 12 + pha];
675 pulseShape[i][0] = shapeFormHG[i * 25 + 12 + pha];
676 pulseShapeDerivative[i][0] = shapeFormDerivativeHG[i * 25 + 12 + pha];
695 Q1 = ((pulseShape.T()) * inverse * pulseShape).determinant();
696 Q2 = ((pulseShapeDerivative.T()) * inverse * pulseShapeDerivative).determinant();
697 Q3 = ((pulseShapeDerivative.T()) * inverse * pulseShape).determinant();
698 Delta = Q1 * Q2 - Q3 * Q3;
701 a = Q2 / Delta * inverse * pulseShape - Q3 / Delta * inverse * pulseShapeDerivative;
702 b = Q3 / Delta * inverse * pulseShape - Q1 / Delta * inverse * pulseShapeDerivative;
708 *aiLoFile << std::setw(6) << pha;
709 for (
int i = 0; i < dignum; ++i)
710 *aiLoFile << std::setw(18) << std::setprecision(10) <<
a(i + 1, 1);
711 *aiLoFile << std::endl;
713 *biLoFile << std::setw(6) << pha;
714 for (
int i = 0; i < dignum; ++i)
715 *biLoFile << std::setw(18) << std::setprecision(10) << b(i + 1, 1);
716 *biLoFile << std::endl;
718 *aiHiFile << std::setw(6) << pha;
719 for (
int i = 0; i < dignum; ++i)
720 *aiHiFile << std::setw(18) << std::setprecision(10) <<
a(i + 1, 1);
721 *aiHiFile << std::endl;
723 *biHiFile << std::setw(6) << pha;
724 for (
int i = 0; i < dignum; ++i)
725 *biHiFile << std::setw(18) << std::setprecision(10) << b(i + 1, 1);
726 *biHiFile << std::endl;
742 std::vector<double> &pulseShapeT,
int dignum) {
747 pulseShape.resize(dignum * 25);
748 ATH_MSG_DEBUG(
"Set dimension of m_pulseShape to dignum*25=" << dignum * 25);
752 double tmin = 10000., tmax = -10000.;
753 size = pulseShapeT.size();
754 for (
int i = 0; i < size; ++i) {
755 if (pulseShapeT[i] < tmin) tmin = pulseShapeT[i];
756 if (pulseShapeT[i] > tmax) tmax = pulseShapeT[i];
757 if (pulseShapeT[i] == 0) nt0 = i;
763 <<
" central point=" << nt0
764 <<
" pulseShapeT[nt0]=" << pulseShapeT[nt0]
765 <<
" pulseShapeY[nt0]=" << pulseShapeY[nt0]);
771 double minn, minp, tdist;
772 pulseShape[dignum * 25 / 2] = pulseShapeY[nt0];
773 for (
int i = 1; i < dignum * 25 / 2 + 1; ++i) {
776 pulseShape[dignum * 25 / 2 - i] = 0.;
783 for (
int j = 0; j < nt0 + 1 && !exact; ++j) {
784 if (pulseShapeT[j] ==
double(-i)) {
785 pulseShape[dignum * 25 / 2 - i] = pulseShapeY[j];
788 tdist = pulseShapeT[j] - double(-i);
789 if (tdist < 0. && tdist > minn) {
793 if (tdist > 0. && tdist < minp) {
802 <<
" pulseShape=" << pulseShape[dignum * 25 / 2 - i]);
807 <<
" nminn=" << nminn
808 <<
" pulseShapeT=" << pulseShapeT[nminn]
809 <<
" pulseShapeY=" << pulseShapeY[nminn] << std::endl
810 <<
" nminp=" << nminp
811 <<
" pulseShapeT=" << pulseShapeT[nminp]
812 <<
" pulseShapeY=" << pulseShapeY[nminp]);
815 pulseShape[dignum * 25 / 2 - i] = pulseShapeY[nminn]
816 + (pulseShapeY[nminp] - pulseShapeY[nminn]) / (pulseShapeT[nminp] - pulseShapeT[nminn])
817 * (-i - pulseShapeT[nminn]);
824 pulseShape[dignum * 25 / 2 + i] = 0.;
831 for (
int j = nt0; j < size && !exact; ++j) {
832 if (pulseShapeT[j] ==
double(i)) {
833 pulseShape[dignum * 25 / 2 + i] = pulseShapeY[j];
836 tdist = pulseShapeT[j] - double(i);
837 if (tdist < 0)
if (tdist > minn) {
841 if (tdist > 0)
if (tdist < minp) {
849 <<
" pulseShape=" << pulseShape[dignum * 25 / 2 + i]);
853 <<
" nminn=" << nminn
854 <<
" pulseShapeT=" << pulseShapeT[nminn]
855 <<
" pulseShapeY=" << pulseShapeY[nminn] << std::endl
856 <<
" nminp=" << nminp
857 <<
" pulseShapeT=" << pulseShapeT[nminp]
858 <<
" pulseShapeY=" << pulseShapeY[nminp]);
861 pulseShape[dignum * 25 / 2 + i] = pulseShapeY[nminn]
862 + (pulseShapeY[nminp] - pulseShapeY[nminn]) / (pulseShapeT[nminp] - pulseShapeT[nminn])
863 * (i - pulseShapeT[nminn]);