26 #include "CLHEP/Matrix/Matrix.h"
75 for (
int i = 0;
i < dignum; ++
i) {
77 for (
int j = 0; j < dignum; ++j) {
111 for (
int i = 0;
i < dignum; ++
i)
112 for (
int j = 0; j < dignum; ++j)
126 dignum = digits.size();
130 for (
int i = 0;
i < dignum; ++
i) {
132 for (
int j = 0; j < dignum; ++j)
178 for (
int i = 0;
i < dignum; ++
i)
180 for (
int j = 0; j < dignum; ++j) {
213 int &dignum,
int chthres) {
217 dignum = digits.size();
229 for (
int i = 0;
i < dignum -
m_lag; ++
i) {
266 <<
" jentry-chthres=" <<
m_jEntry - chthres
267 <<
" lag=1, ros=1, drawer=1, channel=0, gain=1"
295 if (is7to9 && dignum == 7) {
296 for (
int i = 0;
i < 9; ++
i)
300 for (
int i = 0;
i < 9 -
m_lag; ++
i) {
320 for (
int i = 0;
i < dignum;
i++)
324 for (
int i = 0;
i < dignum -
m_lag; ++
i) {
344 std::cout <<
" TileCorrelation::PrintCorrelation()..." << std::endl;
346 std::cout <<
" ros=" <<
ros << std::endl;
348 std::cout <<
" drawer=" <<
drawer << std::endl;
350 std::cout <<
" channel=" <<
channel << std::endl;
352 std::cout <<
" gain=" <<
gain << std::endl;
353 for (
int i = 0;
i < dignum; ++
i) {
354 for (
int j = 0; j < dignum; ++j) {
357 std::cout << std::endl;
359 std::cout << std::endl;
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;
396 <<
" drawer " <<
drawer << MSG::hex
397 <<
" frag0x " << frag << MSG::dec
404 for (
int j = 0; j < dignum; ++j) {
409 correlation[
i][j] = 0.0;
412 *correlationFile <<
"ros " <<
ros
413 <<
" drawer " <<
drawer << std::hex
414 <<
" frag0x " << frag << std::dec
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;
459 <<
" drawer " <<
drawer << MSG::hex
460 <<
" frag0x " << frag << MSG::dec
467 for (
int i = 0;
i < dignum; ++
i)
468 for (
int j = 0; j < dignum; ++j) {
473 correlation[
i][j] = 0.0;
476 *correlationFile <<
"ros " <<
ros
477 <<
" drawer " <<
drawer << std::hex
478 <<
" frag0x " << frag << std::dec
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()) {
553 msg(
MSG::VERBOSE) <<
"shapeFormLG, shapeFormDerivativeLG, shapeFormHG, shapeFormDerivativeHG" <<
endmsg;
554 for (
int i = 0;
i <
int(shapeFormLG.size()); ++
i)
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;
634 *aiLoFile <<
"ros " <<
ros
635 <<
" drawer " <<
drawer << std::hex
636 <<
" frag0x " << frag << std::dec
640 *biLoFile <<
"ros " <<
ros
641 <<
" drawer " <<
drawer << std::hex
642 <<
" frag0x " << frag << std::dec
647 *aiHiFile <<
"ros " <<
ros
648 <<
" drawer " <<
drawer << std::hex
649 <<
" frag0x " << frag << std::dec
653 *biHiFile <<
"ros " <<
ros
654 <<
" drawer " <<
drawer << std::hex
655 <<
" frag0x " << frag << std::dec
660 for (
int pha = -12; pha < 13; ++pha) {
664 for (
int i = 0;
i < dignum; ++
i)
665 for (
int j = 0; j < dignum; ++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]);