diff --git a/RecoLocalCalo/HcalRecAlgos/interface/EigenMatrixTypes.h b/RecoLocalCalo/HcalRecAlgos/interface/EigenMatrixTypes.h index 9ba6a1aa0b7f0..16223083884ac 100644 --- a/RecoLocalCalo/HcalRecAlgos/interface/EigenMatrixTypes.h +++ b/RecoLocalCalo/HcalRecAlgos/interface/EigenMatrixTypes.h @@ -11,8 +11,8 @@ typedef Eigen::Matrix SampleVector; typedef Eigen::Matrix PulseVector; typedef Eigen::Matrix BXVector; -typedef Eigen::Matrix FullSampleVector; -typedef Eigen::Matrix FullSampleMatrix; +typedef Eigen::Matrix FullSampleVector; +typedef Eigen::Matrix FullSampleMatrix; typedef Eigen::Matrix SampleMatrix; typedef Eigen::Matrix PulseMatrix; diff --git a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h index fc89401966b85..a5fdcdf61b8b9 100644 --- a/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h +++ b/RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h @@ -21,6 +21,7 @@ struct MahiNnlsWorkspace { unsigned int tsOffset; unsigned int fullTSOffset; int bxOffset; + int maxoffset; double dt; //holds active bunch crossings @@ -66,7 +67,6 @@ struct MahiNnlsWorkspace { unsigned int nP; PulseVector ampVec; - PulseVector errVec; PulseVector ampvecpermtest; SamplePulseMatrix invcovp; @@ -75,8 +75,6 @@ struct MahiNnlsWorkspace { PulseVector updateWork; // w (vector) SampleDecompLLT covDecomp; - SampleMatrix covDecompLinv; - PulseMatrix topleft_work; PulseDecompLDLT pulseDecomp; }; diff --git a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc index 75a155375929b..76a36ae213692 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc @@ -62,7 +62,7 @@ void MahiFit::phase1Apply(const HBHEChannelInfo& channelData, channelData.tsPedestalWidth(2)*channelData.tsPedestalWidth(2)+ channelData.tsPedestalWidth(3)*channelData.tsPedestalWidth(3) ); - nnlsWork_.pedConstraint = pedVal*SampleMatrix::Ones(nnlsWork_.tsSize, nnlsWork_.tsSize); + nnlsWork_.pedConstraint.setConstant(nnlsWork_.tsSize, nnlsWork_.tsSize, pedVal); nnlsWork_.amplitudes.resize(nnlsWork_.tsSize); nnlsWork_.noiseTerms.resize(nnlsWork_.tsSize); @@ -137,7 +137,10 @@ void MahiFit::doFit(std::array &correctedOutput, int nbx) const { nnlsWork_.bxOffset = bxOffsetConf_; } - nnlsWork_.bxs.resize(bxSize); + nnlsWork_.nPulseTot = bxSize; + + if (dynamicPed_) nnlsWork_.nPulseTot++; + nnlsWork_.bxs.setZero(nnlsWork_.nPulseTot); if (nbx==1) { nnlsWork_.bxs.coeffRef(0) = 0; @@ -148,28 +151,25 @@ void MahiFit::doFit(std::array &correctedOutput, int nbx) const { } } - nnlsWork_.nPulseTot = bxSize; + nnlsWork_.maxoffset = nnlsWork_.bxs.coeffRef(bxSize-1); + if (dynamicPed_) nnlsWork_.bxs[nnlsWork_.nPulseTot-1] = pedestalBX_; - if (dynamicPed_) { - nnlsWork_.nPulseTot++; - nnlsWork_.bxs.resize(nnlsWork_.nPulseTot); - nnlsWork_.bxs[nnlsWork_.nPulseTot-1] = pedestalBX_; - } - - nnlsWork_.pulseMat.resize(nnlsWork_.tsSize,nnlsWork_.nPulseTot); - nnlsWork_.ampVec = PulseVector::Zero(nnlsWork_.nPulseTot); - nnlsWork_.errVec = PulseVector::Zero(nnlsWork_.nPulseTot); + nnlsWork_.pulseMat.setZero(nnlsWork_.tsSize,nnlsWork_.nPulseTot); + nnlsWork_.invcovp.setZero(nnlsWork_.tsSize,nnlsWork_.nPulseTot); + nnlsWork_.ampVec.setZero(nnlsWork_.nPulseTot); + nnlsWork_.ampvecpermtest.setZero(nnlsWork_.nPulseTot); int offset=0; for (unsigned int iBX=0; iBX &correctedOutput, int nbx) const { nnlsWork_.pulseDerivArray[iBX], nnlsWork_.pulseCovArray[iBX]); - nnlsWork_.ampVec.coeffRef(iBX)=0; - - nnlsWork_.pulseMat.col(iBX) = nnlsWork_.pulseShapeArray[iBX].segment(nnlsWork_.fullTSOffset - offset, nnlsWork_.tsSize); + nnlsWork_.pulseMat.col(iBX) = nnlsWork_.pulseShapeArray[iBX].segment(nnlsWork_.maxoffset - offset, nnlsWork_.tsSize); } } - nnlsWork_.pulseMat.col(nnlsWork_.nPulseTot-1) = SampleVector::Ones(nnlsWork_.tsSize); - - nnlsWork_.aTaMat.resize(nnlsWork_.nPulseTot, nnlsWork_.nPulseTot); - nnlsWork_.aTbVec.resize(nnlsWork_.nPulseTot); + nnlsWork_.aTaMat.setZero(nnlsWork_.nPulseTot, nnlsWork_.nPulseTot); + nnlsWork_.aTbVec.setZero(nnlsWork_.nPulseTot); + nnlsWork_.updateWork.setZero(nnlsWork_.nPulseTot); double chiSq = minimize(); - nnlsWork_.residuals = nnlsWork_.pulseMat*nnlsWork_.ampVec - nnlsWork_.amplitudes; - bool foundintime = false; unsigned int ipulseintime = 0; @@ -282,23 +277,23 @@ void MahiFit::updatePulseShape(double itQ, FullSampleVector &pulseShape, FullSam //with previous SOI=TS4 case assumed by psfPtr_->getPulseShape() int delta =nnlsWork_. tsOffset == 3 ? 1 : 0; - for (unsigned int iTS=nnlsWork_.fullTSOffset; iTS