17 #include "CLHEP/Matrix/Matrix.h"
18 #include <CLHEP/Matrix/Vector.h>
24 #include <TRobustEstimator.h>
28 using CLHEP::HepMatrix;
29 using CLHEP::HepVector;
88 for (
int i=0;
i<dignum;
i++){
90 for (
int j=0;j<dignum;j++){
127 for (
int i=0;
i<dignum;
i++)
128 for (
int j=0;j<dignum;j++)
141 dignum=digits.size();
149 for (
int i=0;
i<dignum;
i++){
156 <<
" TileOFCorrelation::Sum, ros="<<
ros
182 TMatrixDSym TempMatrix(dignum);
184 double *Row =
new double[dignum];
195 for (
int j=0; j<dignum;j++){
198 REstimator.AddRow(Row);
200 REstimator.Evaluate();
201 REstimator.GetCovariance(TempMatrix);
202 for (
int i=0;
i<dignum;
i++)
203 for (
int j=0;j<dignum;j++)
207 if(dignum==7 && flag_7to9)
208 for (
int j=0;j<9;j++)
211 for (
int i=0;
i<dignum;
i++){
212 for (
int j=0;j<dignum;j++){
228 if(dignum==7 && flag_7to9)
229 for (
int j=0;j<9;j++)
232 for (
int i=0;
i<dignum;
i++){
233 for (
int j=0;j<dignum;j++){
254 vector < vector <float> >
x;
265 dignum=digits.size();
319 <<
" jentry="<<
m_jentry<<
" jentry-chthres="<<
m_jentry-chthres<<
" m_lag=1, ros=1, drawer=1, channel=0, gain=1"
344 if (flag_7to9 && dignum==7){
345 for (
int i=0;
i<9;
i++)
365 for (
int i=0;
i<dignum;
i++)
388 for(
int i=0;
i<dignum;
i++)
390 for(
int j=0;j<dignum;j++)
401 cout<<
" TileOFCorrelation::PrintCorrelation()..."<<endl;
403 cout<<
" ros="<<
ros<<endl;
405 cout<<
" drawer="<<
drawer<<endl;
407 cout<<
" channel="<<
channel<<endl;
409 cout<<
" gain="<<
gain<<endl;
410 for (
int i=0;
i<dignum;
i++){
411 for (
int j=0;j<dignum;j++){
429 const string& OptFilterFile_CorrelationSumm,
436 HepMatrix M_correlation(dignum,1,0);
438 fstream *f_correlation =
new fstream(OptFilterFile_CorrelationSumm.c_str(),
fstream::out);
439 if (f_correlation->is_open())
log<<MSG::INFO<<OptFilterFile_CorrelationSumm<<
" file open"<<
endmsg;
441 if (deltaCorrelation){
443 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 M_correlation[
i][j]=
m_R[0][0][0][0][
i][j];
448 M_correlation[
i][j]=0.0;
452 *f_correlation<<M_correlation.T()<<endl;
461 <<
" drawer "<<
drawer<<MSG::hex
462 <<
" frag0x "<<frag<<MSG::dec
470 for (
int j=0;j<dignum;j++){
476 M_correlation[
i][j]=0.0;
481 *f_correlation<<
"ros "<<
ros
482 <<
" drawer "<<
drawer<<std::hex
483 <<
" frag0x "<<frag<<std::dec
493 f_correlation->close();
505 if(!deltaCorrelation){
508 for (
int i=1;
i<dignum;
i++){
525 const string& OptFilterFile_CorrelationMatrix,
532 HepMatrix M_correlation(dignum,dignum,0);
534 fstream *f_correlation =
new fstream(OptFilterFile_CorrelationMatrix.c_str(),
fstream::out);
535 if (f_correlation->is_open())
log<<MSG::INFO<<OptFilterFile_CorrelationMatrix<<
" file open"<<
endmsg;
537 if (deltaCorrelation){
538 for (
int i=0;
i<dignum;
i++)
539 for (
int j=0;j<dignum;j++){
540 if (
m_R[0][0][0][0][
i][j]>-100000. &&
m_R[0][0][0][0][
i][j]<100000.)
541 M_correlation[
i][j]=
m_R[0][0][0][0][
i][j];
543 M_correlation[
i][j]=0.0;
547 *f_correlation<<M_correlation<<endl;
556 <<
" drawer "<<
drawer<<MSG::hex
557 <<
" frag0x "<<frag<<MSG::dec
564 for (
int i=0;
i<dignum;
i++)
565 for (
int j=0;j<dignum;j++){
570 M_correlation[
i][j]=0.0;
575 *f_correlation<<
"ros "<<
ros
576 <<
" drawer "<<
drawer<<std::hex
577 <<
" frag0x "<<frag<<std::dec
581 <<M_correlation<<endl;
586 f_correlation->close();
593 bool deltaCorrelation,
594 const vector<double>& LshapeForm,
595 const vector<double>& HshapeForm,
596 const vector<double>& LdshapeForm,
597 const vector<double>& HdshapeForm,
607 HepMatrix Correlation(dignum,dignum,0), Inverse(dignum,dignum,0),
Zero(dignum,dignum,0);
608 HepMatrix PulseShape(dignum,1,0), DPulseShape(dignum,1,0);
609 HepMatrix
a(dignum,1,0),
b(dignum,1,0), WeightZero(dignum,1,0);
614 log<<MSG::INFO<<
" Calculating OF2 weights ";
616 log<<MSG::INFO<<
" Calculating OF1 weights ";
617 if (deltaCorrelation)
618 log<<MSG::INFO<<
"with Delta correlation matrix "<<
endmsg;
620 log<<MSG::INFO<<
"with correlation matrix obtained from data "<<
endmsg;
626 if (deltaCorrelation){
628 for (
int i=0;
i<dignum;
i++)
629 for (
int j=0;j<dignum;j++)
630 Correlation[
i][j]=
m_R[0][0][0][0][
i][j];
632 Inverse=Correlation.inverse(ierr);
634 for (
int pha=-100;pha<101;pha++){
637 for (
int i=0;
i<dignum;
i++){
639 PulseShape[
i][0]=LshapeForm[
i*25+100+pha];
640 DPulseShape[
i][0]=LdshapeForm[
i*25+100+pha];
642 PulseShape[
i][0]=HshapeForm[
i*25+100+pha];
643 DPulseShape[
i][0]=HdshapeForm[
i*25+100+pha];
645 g[pha+100][
i]=PulseShape[
i][0];
651 HepMatrix SystemMatrix(dignum+2,dignum+2,0);
652 HepVector
Result(dignum + 2,0);
653 HepVector IndependTermsAmp(dignum+2,0);
654 HepVector IndependTermsTime(dignum+2,0);
655 for (
int i=0;
i<dignum;
i++){
657 for (
int j=0; j<dignum; j++){
658 SystemMatrix[
i][j]=Correlation[
i][j];
660 SystemMatrix[dignum][
i] = PulseShape[
i][0];
661 SystemMatrix[
i][dignum] = PulseShape[
i][0];
663 SystemMatrix[dignum + 1][
i] = DPulseShape[
i][0];
664 SystemMatrix[
i][dignum + 1] = DPulseShape[
i][0];
668 IndependTermsAmp[dignum] = 1.;
669 IndependTermsTime[dignum+1] = -1.;
671 Result = solve(SystemMatrix,IndependTermsAmp);
673 for (
int ismp=0; ismp<dignum; ismp++) {
677 Result = solve(SystemMatrix,IndependTermsTime);
679 for (
int ismp=0; ismp<dignum; ismp++){
684 HepMatrix SystemMatrix(dignum+3,dignum+3,0);
685 HepVector
Result(dignum+3,0);
686 HepVector IndependTermsAmp(dignum+3,0);
687 HepVector IndependTermsTime(dignum+3,0);
688 HepVector IndependTermsPed(dignum+3,0);
690 for (
int i=0;
i<dignum;
i++){
691 for (
int j=0; j<dignum; j++){
692 SystemMatrix[
i][j]=Correlation[
i][j];
694 SystemMatrix[dignum][
i] = PulseShape[
i][0];
695 SystemMatrix[
i][dignum] = PulseShape[
i][0];
697 SystemMatrix[dignum + 1][
i] = DPulseShape[
i][0];
698 SystemMatrix[
i][dignum + 1] = DPulseShape[
i][0];
701 SystemMatrix[dignum + 2][
i] = 1.;
702 SystemMatrix[
i][dignum + 2] = 1.;
705 IndependTermsAmp[dignum] = 1.;
706 IndependTermsTime[dignum+1] = -1.;
707 IndependTermsPed[dignum+2] = 1.;
709 Result = solve(SystemMatrix,IndependTermsAmp);
711 for (
int ismp=0; ismp<dignum; ismp++) {
715 Result = solve(SystemMatrix,IndependTermsTime);
717 for (
int ismp=0; ismp<dignum; ismp++){
721 Result = solve(SystemMatrix,IndependTermsPed);
723 for (
int ismp=0; ismp<dignum; ismp++){
729 for (
int ismp=0; ismp<dignum; ismp++){
730 w_a[pha+100][ismp] = 0.;
731 w_b[pha+100][ismp] = 0.;
732 w_c[pha+100][ismp] = 0.;
737 printf(
"ros=%d drawer=%d channel=%d gain=%d N=%d \n",
741 for (
int i=0;
i<dignum;
i++)
742 for (
int j=0;j<dignum;j++)
745 Inverse=Correlation.inverse(ierr);
748 log<<MSG::WARNING<<
" Correlation matrix cannot be inverted for ros="<<
ros<<
" drawer="<<
drawer
750 log<<MSG::WARNING<<
" Weights set to zero for this channel..."<<
endmsg;
751 for (
int i=0;
i<dignum;
i++){
752 for (
int j=0; j<dignum; j++){
759 for (
int pha=-100;pha<101;pha++){
762 for (
int i=0;
i<dignum;
i++){
764 PulseShape[
i][0]=LshapeForm[
i*25+100+pha];
765 DPulseShape[
i][0]=LdshapeForm[
i*25+100+pha];
767 PulseShape[
i][0]=HshapeForm[
i*25+100+pha];
768 DPulseShape[
i][0]=HdshapeForm[
i*25+100+pha];
770 g[pha+100][
i]=PulseShape[
i][0];
775 HepMatrix SystemMatrix(dignum+2,dignum+2,0);
776 HepVector
Result(dignum + 2,0);
777 HepVector IndependTermsAmp(dignum+2,0);
778 HepVector IndependTermsTime(dignum+2,0);
780 for (
int i=0;
i<dignum;
i++){
781 for (
int j=0; j<dignum; j++){
782 SystemMatrix[
i][j]=Correlation[
i][j];
784 SystemMatrix[dignum][
i] = PulseShape[
i][0];
785 SystemMatrix[
i][dignum] = PulseShape[
i][0];
787 SystemMatrix[dignum + 1][
i] = DPulseShape[
i][0];
788 SystemMatrix[
i][dignum + 1] = DPulseShape[
i][0];
791 IndependTermsAmp[dignum] = 1.;
792 IndependTermsTime[dignum+1] = -1.;
794 Result = solve(SystemMatrix,IndependTermsAmp);
796 for (
int ismp=0; ismp<dignum; ismp++) {
800 Result = solve(SystemMatrix,IndependTermsTime);
802 for (
int ismp=0; ismp<dignum; ismp++){
807 HepMatrix SystemMatrix(dignum+3,dignum+3,0);
808 HepVector
Result(dignum + 3,0);
809 HepVector IndependTermsAmp(dignum+3,0);
810 HepVector IndependTermsTime(dignum+3,0);
811 HepVector IndependTermsPed(dignum+3,0);
813 for (
int i=0;
i<dignum;
i++){
814 for (
int j=0; j<dignum; j++){
815 SystemMatrix[
i][j]=Correlation[
i][j];
817 SystemMatrix[dignum][
i] = PulseShape[
i][0];
818 SystemMatrix[
i][dignum] = PulseShape[
i][0];
820 SystemMatrix[dignum + 1][
i] = DPulseShape[
i][0];
821 SystemMatrix[
i][dignum + 1] = DPulseShape[
i][0];
823 SystemMatrix[dignum + 2][
i] = 1.;
824 SystemMatrix[
i][dignum + 2] = 1.;
827 IndependTermsAmp[dignum] = 1.;
828 IndependTermsTime[dignum+1] = -1.;
829 IndependTermsPed[dignum+2] = 1.;
831 Result = solve(SystemMatrix,IndependTermsAmp);
833 for (
int ismp=0; ismp<dignum; ismp++) {
837 Result = solve(SystemMatrix,IndependTermsTime);
839 for (
int ismp=0; ismp<dignum; ismp++){
843 Result = solve(SystemMatrix,IndependTermsPed);
845 for (
int ismp=0; ismp<dignum; ismp++){
850 printf(
"PHASE: %d deter=%f\n",pha,SystemMatrix.determinant());
851 for (
int i=0;
i<dignum+2;
i++){
852 for (
int j=0; j<dignum+2; j++){
853 printf(
"%4.4f ",SystemMatrix[
i][j]);
857 for (
int i=0;
i<dignum;
i++){
858 for (
int j=0; j<dignum; j++){
863 for (
int ismp=0; ismp<dignum; ismp++)
864 printf(
"samp: %d ->w_a: %f w_b:%f\n",
865 ismp,
w_a[pha+100][ismp],
w_b[pha+100][ismp]);
870 for (
int ismp=0; ismp<dignum; ismp++){
871 w_a[pha+100][ismp] = 0.;
872 w_b[pha+100][ismp] = 0.;
873 w_c[pha+100][ismp] = 0.;
890 vector<double> &pulseShapeX,
891 vector<double> &pulseShapeT,
898 int shape_size=(dignum-1)*25+200;
899 pulseShape.resize(shape_size);
904 double tmin=10000., tmax=-10000.;
905 int ntmin=0, ntmax=0;
906 size=pulseShapeT.size();
908 if (pulseShapeT[
i]<tmin) {
912 if (pulseShapeT[
i]>tmax){
916 if (pulseShapeT[
i]==0) nt0=
i;
918 log<<
MSG::DEBUG<<
"pulseShapeX & pulseShapeT size ="<<
size<<
", tmin="<<tmin<<
", tmax="<<tmax<<
" nt0="<<nt0<<
" pulseShapeT[nt0]="<<pulseShapeT[nt0]<<
" pulseShapeX[nt0]="<<pulseShapeX[nt0]<<
endmsg;
923 double minn, minp, tdist;
924 pulseShape[(shape_size)/2]=pulseShapeX[nt0];
926 for (
int i=1;
i<(shape_size)/2;
i++){
929 if (-
i<tmin) pulseShape[(shape_size)/2-
i]=pulseShapeX[ntmin];
936 for (
int j=0;j<nt0+1&&!exact;j++){
937 if (pulseShapeT[j]==
double(-
i)){
938 pulseShape[(shape_size)/2-
i]=pulseShapeX[j];
941 tdist=pulseShapeT[j]-
double(-
i);
942 if (tdist < 0. && tdist>minn){
946 if (tdist>0. && tdist<minp){
958 <<
" nminn="<<nminn<<
" pulseShapeT="<<pulseShapeT[nminn]<<
" pulseShapeX="<<pulseShapeX[nminn]<<std::endl
959 <<
" nminp="<<nminp<<
" pulseShapeT="<<pulseShapeT[nminp]<<
" pulseShapeX="<<pulseShapeX[nminp]<<
endmsg;
960 pulseShape[(shape_size)/2-
i]=pulseShapeX[nminn]+(pulseShapeX[nminp]-pulseShapeX[nminn])/(pulseShapeT[nminp]-pulseShapeT[nminn])*(-
i-pulseShapeT[nminn]);
967 if (
i>tmax) pulseShape[(shape_size)/2+
i]=pulseShapeX[ntmax];
974 for (
int j=nt0;j<
size&&!exact;j++){
975 if (pulseShapeT[j]==
double(
i)){
976 pulseShape[(shape_size)/2+
i]=pulseShapeX[j];
979 tdist=pulseShapeT[j]-
double(
i);
998 <<
" nminn="<<nminn<<
" pulseShapeT="<<pulseShapeT[nminn]<<
" pulseShapeX="<<pulseShapeX[nminn]<<std::endl
999 <<
" nminp="<<nminp<<
" pulseShapeT="<<pulseShapeT[nminp]<<
" pulseShapeX="<<pulseShapeX[nminp]<<
endmsg;
1001 pulseShape[(shape_size)/2+
i]=pulseShapeX[nminn]+(pulseShapeX[nminp]-pulseShapeX[nminn])/(pulseShapeT[nminp]-pulseShapeT[nminn])*(
i-pulseShapeT[nminn]);
1011 const string& OptFilterFile_ai_lo,
1012 const string& OptFilterFile_bi_lo,
1013 const string& OptFilterFile_ai_hi,
1014 const string& OptFilterFile_bi_hi,
1023 fstream *f_ai_lo =
new fstream(OptFilterFile_ai_lo.c_str(), fstream::app|
fstream::out);
1024 fstream *f_bi_lo =
new fstream(OptFilterFile_bi_lo.c_str(), fstream::app|
fstream::out);
1025 fstream *f_ai_hi =
new fstream(OptFilterFile_ai_hi.c_str(), fstream::app|
fstream::out);
1026 fstream *f_bi_hi =
new fstream(OptFilterFile_bi_hi.c_str(), fstream::app|
fstream::out);
1029 if (f_ai_lo->is_open()&&f_ai_lo->is_open()&&f_ai_lo->is_open()&&f_ai_lo->is_open())
log<<MSG::INFO<<
" Weights files open"<<
endmsg;
1030 else log<<MSG::INFO<<
" Weights files didn't open succesfully"<<
endmsg;
1035 std::cout<<
"ros "<<
ros
1036 <<
" drawer "<<
drawer<<std::hex
1037 <<
" frag0x "<<frag<<std::dec
1041 *f_ai_lo<<
"ros "<<
ros
1042 <<
" drawer "<<
drawer<<std::hex
1043 <<
" frag0x "<<frag<<std::dec
1048 *f_bi_lo<<
"ros "<<
ros
1049 <<
" drawer "<<
drawer<<std::hex
1050 <<
" frag0x "<<frag<<std::dec
1056 *f_ai_hi<<
"ros "<<
ros
1057 <<
" drawer "<<
drawer<<std::hex
1058 <<
" frag0x "<<frag<<std::dec
1063 *f_bi_hi<<
"ros "<<
ros
1064 <<
" drawer "<<
drawer<<std::hex
1065 <<
" frag0x "<<frag<<std::dec
1074 for (
int pha=-100;pha<101;pha++){
1076 *f_ai_lo<<std::setw(6)<<pha;
1077 for (
int i=0;
i<dignum;
i++) *f_ai_lo<<std::setw(18)<<std::setprecision(10)<<
w_a[pha+100][
i];
1080 *f_bi_lo<<std::setw(6)<<pha;
1081 for (
int i=0;
i<dignum;
i++) *f_bi_lo<<std::setw(18)<<std::setprecision(10)<<
w_b[pha+100][
i];
1084 *f_ai_hi<<std::setw(6)<<pha;
1085 for (
int i=0;
i<dignum;
i++) *f_ai_hi<<std::setw(18)<<std::setprecision(10)<<
w_a[pha+100][
i];
1088 *f_bi_hi<<std::setw(6)<<pha;
1089 for (
int i=0;
i<dignum;
i++) *f_bi_hi<<std::setw(18)<<std::setprecision(10)<<
w_b[pha+100][
i];
1106 const string& OptFilterFile_bi_lo,
1107 const string& OptFilterFile_ai_hi,
1108 const string& OptFilterFile_bi_hi,
1112 fstream *f_ai_lo =
new fstream(OptFilterFile_ai_lo.c_str(), fstream::trunc|
fstream::out);
1113 fstream *f_bi_lo =
new fstream(OptFilterFile_bi_lo.c_str(), fstream::trunc|
fstream::out);
1114 fstream *f_ai_hi =
new fstream(OptFilterFile_ai_hi.c_str(), fstream::trunc|
fstream::out);
1115 fstream *f_bi_hi =
new fstream(OptFilterFile_bi_hi.c_str(), fstream::trunc|
fstream::out);
1117 f_ai_lo->open(OptFilterFile_ai_hi.c_str(), fstream::trunc|
fstream::out);
1118 f_bi_lo->open(OptFilterFile_ai_hi.c_str(), fstream::trunc|
fstream::out);
1119 f_ai_hi->open(OptFilterFile_ai_hi.c_str(), fstream::trunc|
fstream::out);
1120 f_bi_hi->open(OptFilterFile_ai_hi.c_str(), fstream::trunc|
fstream::out);