File indexing completed on 2021-11-28 23:48:59
0001 #define EIGEN_NO_DEBUG
0002 #include "RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 MahiFit::MahiFit() {}
0006
0007 void MahiFit::setParameters(bool iDynamicPed,
0008 double iTS4Thresh,
0009 double chiSqSwitch,
0010 bool iApplyTimeSlew,
0011 HcalTimeSlew::BiasSetting slewFlavor,
0012 bool iCalculateArrivalTime,
0013 double iMeanTime,
0014 double iTimeSigmaHPD,
0015 double iTimeSigmaSiPM,
0016 const std::vector<int>& iActiveBXs,
0017 int iNMaxItersMin,
0018 int iNMaxItersNNLS,
0019 double iDeltaChiSqThresh,
0020 double iNnlsThresh) {
0021 dynamicPed_ = iDynamicPed;
0022
0023 ts4Thresh_ = iTS4Thresh;
0024 chiSqSwitch_ = chiSqSwitch;
0025
0026 applyTimeSlew_ = iApplyTimeSlew;
0027 slewFlavor_ = slewFlavor;
0028
0029 calculateArrivalTime_ = iCalculateArrivalTime;
0030 meanTime_ = iMeanTime;
0031 timeSigmaHPD_ = iTimeSigmaHPD;
0032 timeSigmaSiPM_ = iTimeSigmaSiPM;
0033
0034 activeBXs_ = iActiveBXs;
0035
0036 nMaxItersMin_ = iNMaxItersMin;
0037 nMaxItersNNLS_ = iNMaxItersNNLS;
0038
0039 deltaChiSqThresh_ = iDeltaChiSqThresh;
0040 nnlsThresh_ = iNnlsThresh;
0041
0042 bxOffsetConf_ = -(*std::min_element(activeBXs_.begin(), activeBXs_.end()));
0043 bxSizeConf_ = activeBXs_.size();
0044 }
0045
0046 void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
0047 float& reconstructedEnergy,
0048 float& soiPlusOneEnergy,
0049 float& reconstructedTime,
0050 bool& useTriple,
0051 float& chi2) const {
0052 const unsigned nSamples = channelData.nSamples();
0053 const unsigned soi = channelData.soi();
0054
0055
0056 assert(nSamples == 8 || nSamples == 10);
0057 assert(nnlsWork_.tsSize <= nSamples);
0058 assert(soi + 1U < nnlsWork_.tsSize);
0059
0060 resetWorkspace();
0061
0062 nnlsWork_.tsOffset = soi;
0063
0064
0065 std::array<float, 4> reconstructedVals{{0.f, -9999.f, -9999.f, 0.f}};
0066
0067 double tsTOT = 0, tstrig = 0;
0068 for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
0069 auto const amplitude = channelData.tsRawCharge(iTS) - channelData.tsPedestal(iTS);
0070
0071 nnlsWork_.amplitudes.coeffRef(iTS) = amplitude;
0072
0073
0074 auto const noiseADC = norm_ * channelData.tsDFcPerADC(iTS);
0075
0076
0077 auto const pedWidth = channelData.tsPedestalWidth(iTS);
0078
0079
0080 auto const noisePhotoSq = (amplitude > pedWidth) ? (amplitude * channelData.fcByPE()) : 0.f;
0081
0082
0083 nnlsWork_.noiseTerms.coeffRef(iTS) = noiseADC * noiseADC + noisePhotoSq + pedWidth * pedWidth;
0084
0085 nnlsWork_.pedVals.coeffRef(iTS) = pedWidth;
0086
0087 tsTOT += amplitude;
0088 if (iTS == soi)
0089 tstrig += amplitude;
0090 }
0091
0092
0093
0094 const float gain0 = channelData.tsGain(0);
0095 tsTOT *= gain0;
0096 tstrig *= gain0;
0097
0098 useTriple = false;
0099 if (tstrig > ts4Thresh_ && tsTOT > 0) {
0100 nnlsWork_.noisecorr = channelData.noisecorr();
0101
0102 if (nnlsWork_.noisecorr != 0) {
0103 nnlsWork_.pedVal = 0.f;
0104 } else {
0105
0106 nnlsWork_.pedVal = 0.25f * (channelData.tsPedestalWidth(0) * channelData.tsPedestalWidth(0) +
0107 channelData.tsPedestalWidth(1) * channelData.tsPedestalWidth(1) +
0108 channelData.tsPedestalWidth(2) * channelData.tsPedestalWidth(2) +
0109 channelData.tsPedestalWidth(3) * channelData.tsPedestalWidth(3));
0110 }
0111
0112
0113 if (chiSqSwitch_ > 0) {
0114 doFit(reconstructedVals, 1);
0115 if (reconstructedVals[2] > chiSqSwitch_) {
0116 doFit(reconstructedVals, 0);
0117 useTriple = true;
0118 }
0119 } else {
0120 doFit(reconstructedVals, 0);
0121 useTriple = true;
0122 }
0123 }
0124
0125 reconstructedEnergy = reconstructedVals[0] * gain0;
0126 soiPlusOneEnergy = reconstructedVals[3] * gain0;
0127 reconstructedTime = reconstructedVals[1];
0128 chi2 = reconstructedVals[2];
0129 }
0130
0131 void MahiFit::doFit(std::array<float, 4>& correctedOutput, int nbx) const {
0132 unsigned int bxSize = 1;
0133
0134 if (nbx == 1) {
0135 nnlsWork_.bxOffset = 0;
0136 } else {
0137 bxSize = bxSizeConf_;
0138 nnlsWork_.bxOffset = static_cast<int>(nnlsWork_.tsOffset) >= bxOffsetConf_ ? bxOffsetConf_ : nnlsWork_.tsOffset;
0139 }
0140
0141 nnlsWork_.nPulseTot = bxSize;
0142
0143 if (dynamicPed_)
0144 nnlsWork_.nPulseTot++;
0145 nnlsWork_.bxs.setZero(nnlsWork_.nPulseTot);
0146
0147 if (nbx == 1) {
0148 nnlsWork_.bxs.coeffRef(0) = 0;
0149 } else {
0150 for (unsigned int iBX = 0; iBX < bxSize; ++iBX) {
0151 nnlsWork_.bxs.coeffRef(iBX) =
0152 activeBXs_[iBX] -
0153 ((static_cast<int>(nnlsWork_.tsOffset) + activeBXs_[0]) >= 0 ? 0 : (nnlsWork_.tsOffset + activeBXs_[0]));
0154 }
0155 }
0156
0157 nnlsWork_.maxoffset = nnlsWork_.bxs.coeff(bxSize - 1);
0158 if (dynamicPed_)
0159 nnlsWork_.bxs[nnlsWork_.nPulseTot - 1] = pedestalBX_;
0160
0161 nnlsWork_.pulseMat.setZero(nnlsWork_.tsSize, nnlsWork_.nPulseTot);
0162 nnlsWork_.pulseDerivMat.setZero(nnlsWork_.tsSize, nnlsWork_.nPulseTot);
0163
0164 FullSampleVector pulseShapeArray;
0165 FullSampleVector pulseDerivArray;
0166 FullSampleMatrix pulseCov;
0167
0168 int offset = 0;
0169 for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0170 offset = nnlsWork_.bxs.coeff(iBX);
0171
0172 if (offset == pedestalBX_) {
0173 nnlsWork_.pulseMat.col(iBX) = SampleVector::Ones(nnlsWork_.tsSize);
0174 nnlsWork_.pulseDerivMat.col(iBX) = SampleVector::Zero(nnlsWork_.tsSize);
0175 } else {
0176 pulseShapeArray.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
0177 if (calculateArrivalTime_)
0178 pulseDerivArray.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
0179 pulseCov.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset,
0180 nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
0181 nnlsWork_.pulseCovArray[iBX].setZero(nnlsWork_.tsSize, nnlsWork_.tsSize);
0182
0183 updatePulseShape(
0184 nnlsWork_.amplitudes.coeff(nnlsWork_.tsOffset + offset), pulseShapeArray, pulseDerivArray, pulseCov);
0185
0186 nnlsWork_.pulseMat.col(iBX) = pulseShapeArray.segment(nnlsWork_.maxoffset - offset, nnlsWork_.tsSize);
0187 if (calculateArrivalTime_)
0188 nnlsWork_.pulseDerivMat.col(iBX) = pulseDerivArray.segment(nnlsWork_.maxoffset - offset, nnlsWork_.tsSize);
0189 nnlsWork_.pulseCovArray[iBX] = pulseCov.block(
0190 nnlsWork_.maxoffset - offset, nnlsWork_.maxoffset - offset, nnlsWork_.tsSize, nnlsWork_.tsSize);
0191 }
0192 }
0193
0194 const float chiSq = minimize();
0195
0196 bool foundintime = false;
0197 unsigned int ipulseintime = 0;
0198
0199 for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0200 if (nnlsWork_.bxs.coeff(iBX) == 0) {
0201 ipulseintime = iBX;
0202 foundintime = true;
0203 break;
0204 }
0205 }
0206
0207 if (foundintime) {
0208 correctedOutput.at(0) = nnlsWork_.ampVec.coeff(ipulseintime);
0209 correctedOutput.at(3) = nnlsWork_.ampVec.coeff(ipulseintime + 1U);
0210 if (correctedOutput.at(0) != 0) {
0211
0212 float arrivalTime = 0.f;
0213 if (calculateArrivalTime_)
0214 arrivalTime = calculateArrivalTime(ipulseintime);
0215 correctedOutput.at(1) = arrivalTime;
0216 } else
0217 correctedOutput.at(1) = -9999.f;
0218
0219 correctedOutput.at(2) = chiSq;
0220 }
0221 }
0222
0223 const float MahiFit::minimize() const {
0224 nnlsWork_.invcovp.setZero(nnlsWork_.tsSize, nnlsWork_.nPulseTot);
0225 nnlsWork_.ampVec.setZero(nnlsWork_.nPulseTot);
0226
0227 SampleMatrix invCovMat;
0228 invCovMat.setConstant(nnlsWork_.tsSize, nnlsWork_.tsSize, nnlsWork_.pedVal);
0229 invCovMat += nnlsWork_.noiseTerms.asDiagonal();
0230
0231
0232 if (nnlsWork_.noisecorr != 0) {
0233 for (unsigned int i = 1; i < nnlsWork_.tsSize; ++i) {
0234 auto const noiseCorrTerm = nnlsWork_.noisecorr * nnlsWork_.pedVals.coeff(i - 1) * nnlsWork_.pedVals.coeff(i);
0235 invCovMat(i - 1, i) += noiseCorrTerm;
0236 invCovMat(i, i - 1) += noiseCorrTerm;
0237 }
0238 }
0239
0240 float oldChiSq = 9999;
0241 float chiSq = oldChiSq;
0242
0243 for (int iter = 1; iter < nMaxItersMin_; ++iter) {
0244 updateCov(invCovMat);
0245
0246 if (nnlsWork_.nPulseTot > 1) {
0247 nnls();
0248 } else {
0249 onePulseMinimize();
0250 }
0251
0252 const float newChiSq = calculateChiSq();
0253 const float deltaChiSq = newChiSq - chiSq;
0254
0255 if (newChiSq == oldChiSq && newChiSq < chiSq) {
0256 break;
0257 }
0258 oldChiSq = chiSq;
0259 chiSq = newChiSq;
0260
0261 if (std::abs(deltaChiSq) < deltaChiSqThresh_)
0262 break;
0263 }
0264
0265 return chiSq;
0266 }
0267
0268 void MahiFit::updatePulseShape(const float itQ,
0269 FullSampleVector& pulseShape,
0270 FullSampleVector& pulseDeriv,
0271 FullSampleMatrix& pulseCov) const {
0272
0273 if (itQ <= 0.f)
0274 return;
0275
0276 float t0 = meanTime_;
0277
0278 if (applyTimeSlew_) {
0279 if (itQ <= 1.f)
0280 t0 += tsDelay1GeV_;
0281 else
0282 t0 += hcalTimeSlewDelay_->delay(float(itQ), slewFlavor_);
0283 }
0284
0285 std::array<double, hcal::constants::maxSamples> pulseN;
0286 std::array<double, hcal::constants::maxSamples> pulseM;
0287 std::array<double, hcal::constants::maxSamples> pulseP;
0288
0289 const float xx = t0;
0290 const float xxm = -nnlsWork_.dt + t0;
0291 const float xxp = nnlsWork_.dt + t0;
0292
0293 psfPtr_->singlePulseShapeFuncMahi(&xx);
0294 psfPtr_->getPulseShape(pulseN);
0295
0296 psfPtr_->singlePulseShapeFuncMahi(&xxm);
0297 psfPtr_->getPulseShape(pulseM);
0298
0299 psfPtr_->singlePulseShapeFuncMahi(&xxp);
0300 psfPtr_->getPulseShape(pulseP);
0301
0302
0303
0304 int delta = 4 - nnlsWork_.tsOffset;
0305
0306 auto invDt = 0.5 / nnlsWork_.dt;
0307
0308 for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
0309 pulseShape(iTS + nnlsWork_.maxoffset) = pulseN[iTS + delta];
0310 if (calculateArrivalTime_)
0311 pulseDeriv(iTS + nnlsWork_.maxoffset) = (pulseM[iTS + delta] - pulseP[iTS + delta]) * invDt;
0312
0313 pulseM[iTS + delta] -= pulseN[iTS + delta];
0314 pulseP[iTS + delta] -= pulseN[iTS + delta];
0315 }
0316
0317 for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
0318 for (unsigned int jTS = 0; jTS < iTS + 1; ++jTS) {
0319 auto const tmp = 0.5 * (pulseP[iTS + delta] * pulseP[jTS + delta] + pulseM[iTS + delta] * pulseM[jTS + delta]);
0320
0321 pulseCov(jTS + nnlsWork_.maxoffset, iTS + nnlsWork_.maxoffset) = tmp;
0322 if (iTS != jTS)
0323 pulseCov(iTS + nnlsWork_.maxoffset, jTS + nnlsWork_.maxoffset) = tmp;
0324 }
0325 }
0326 }
0327
0328 void MahiFit::updateCov(const SampleMatrix& samplecov) const {
0329 SampleMatrix invCovMat = samplecov;
0330
0331 for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0332 auto const amp = nnlsWork_.ampVec.coeff(iBX);
0333 if (amp == 0)
0334 continue;
0335
0336 int offset = nnlsWork_.bxs.coeff(iBX);
0337
0338 if (offset == pedestalBX_)
0339 continue;
0340 else {
0341 invCovMat.noalias() += amp * amp * nnlsWork_.pulseCovArray[offset + nnlsWork_.bxOffset];
0342 }
0343 }
0344
0345 nnlsWork_.covDecomp.compute(invCovMat);
0346 }
0347
0348 float MahiFit::calculateArrivalTime(const unsigned int itIndex) const {
0349 if (nnlsWork_.nPulseTot > 1) {
0350 SamplePulseMatrix pulseDerivMatTMP = nnlsWork_.pulseDerivMat;
0351 for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0352 nnlsWork_.pulseDerivMat.col(iBX) = pulseDerivMatTMP.col(nnlsWork_.bxs.coeff(iBX) + nnlsWork_.bxOffset);
0353 }
0354 }
0355
0356 for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0357 nnlsWork_.pulseDerivMat.col(iBX) *= nnlsWork_.ampVec.coeff(iBX);
0358 }
0359
0360
0361 SampleVector residuals = nnlsWork_.pulseMat * nnlsWork_.ampVec - nnlsWork_.amplitudes;
0362 PulseVector solution = nnlsWork_.pulseDerivMat.colPivHouseholderQr().solve(residuals);
0363 float t = std::clamp((float)solution.coeff(itIndex), -timeLimit_, timeLimit_);
0364
0365 return t;
0366 }
0367
0368 void MahiFit::nnls() const {
0369 const unsigned int npulse = nnlsWork_.nPulseTot;
0370 const unsigned int nsamples = nnlsWork_.tsSize;
0371
0372 PulseVector updateWork;
0373 PulseVector ampvecpermtest;
0374
0375 nnlsWork_.invcovp = nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat);
0376 nnlsWork_.aTaMat.noalias() = nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.invcovp);
0377 nnlsWork_.aTbVec.noalias() =
0378 nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes));
0379
0380 int iter = 0;
0381 Index idxwmax = 0;
0382 float wmax = 0.0f;
0383 float threshold = nnlsThresh_;
0384
0385 while (true) {
0386 if (iter > 0 || nnlsWork_.nP == 0) {
0387 if (nnlsWork_.nP == std::min(npulse, nsamples))
0388 break;
0389
0390 const unsigned int nActive = npulse - nnlsWork_.nP;
0391
0392 if (nActive == 0)
0393 break;
0394
0395 updateWork.noalias() = nnlsWork_.aTbVec - nnlsWork_.aTaMat.lazyProduct(nnlsWork_.ampVec);
0396
0397 Index idxwmaxprev = idxwmax;
0398 float wmaxprev = wmax;
0399 wmax = updateWork.tail(nActive).maxCoeff(&idxwmax);
0400
0401 if (wmax < threshold || (idxwmax == idxwmaxprev && wmax == wmaxprev)) {
0402 break;
0403 }
0404
0405 if (iter >= nMaxItersNNLS_) {
0406 break;
0407 }
0408
0409
0410 idxwmax += nnlsWork_.nP;
0411 nnlsUnconstrainParameter(idxwmax);
0412 }
0413
0414 while (true) {
0415 if (nnlsWork_.nP == 0)
0416 break;
0417
0418 ampvecpermtest.noalias() = nnlsWork_.ampVec.head(nnlsWork_.nP);
0419
0420 solveSubmatrix(nnlsWork_.aTaMat, nnlsWork_.aTbVec, ampvecpermtest, nnlsWork_.nP);
0421
0422
0423 if (ampvecpermtest.head(nnlsWork_.nP).minCoeff() > 0.f) {
0424 nnlsWork_.ampVec.head(nnlsWork_.nP) = ampvecpermtest.head(nnlsWork_.nP);
0425 break;
0426 }
0427
0428
0429 Index minratioidx = 0;
0430
0431
0432 float minratio = std::numeric_limits<float>::max();
0433 for (unsigned int ipulse = 0; ipulse < nnlsWork_.nP; ++ipulse) {
0434 if (ampvecpermtest.coeff(ipulse) <= 0.f) {
0435 const float c_ampvec = nnlsWork_.ampVec.coeff(ipulse);
0436 const float ratio = c_ampvec / (c_ampvec - ampvecpermtest.coeff(ipulse));
0437 if (ratio < minratio) {
0438 minratio = ratio;
0439 minratioidx = ipulse;
0440 }
0441 }
0442 }
0443 nnlsWork_.ampVec.head(nnlsWork_.nP) +=
0444 minratio * (ampvecpermtest.head(nnlsWork_.nP) - nnlsWork_.ampVec.head(nnlsWork_.nP));
0445
0446
0447 nnlsWork_.ampVec.coeffRef(minratioidx) = 0.f;
0448
0449 nnlsConstrainParameter(minratioidx);
0450 }
0451
0452 ++iter;
0453
0454
0455
0456 if (iter % 10 == 0) {
0457 threshold *= 10.;
0458 }
0459 }
0460 }
0461
0462 void MahiFit::onePulseMinimize() const {
0463 nnlsWork_.invcovp = nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat);
0464
0465 float aTaCoeff = (nnlsWork_.invcovp.transpose() * nnlsWork_.invcovp).coeff(0);
0466 float aTbCoeff =
0467 (nnlsWork_.invcovp.transpose() * (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes))).coeff(0);
0468
0469 nnlsWork_.ampVec.coeffRef(0) = std::max(0.f, aTbCoeff / aTaCoeff);
0470 }
0471
0472 float MahiFit::calculateChiSq() const {
0473 return (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat * nnlsWork_.ampVec - nnlsWork_.amplitudes))
0474 .squaredNorm();
0475 }
0476
0477 void MahiFit::setPulseShapeTemplate(const int pulseShapeId,
0478 const HcalPulseShapes& ps,
0479 const bool hasTimeInfo,
0480 const HcalTimeSlew* hcalTimeSlewDelay,
0481 const unsigned int nSamples) {
0482 if (hcalTimeSlewDelay != hcalTimeSlewDelay_) {
0483 assert(hcalTimeSlewDelay);
0484 hcalTimeSlewDelay_ = hcalTimeSlewDelay;
0485 tsDelay1GeV_ = hcalTimeSlewDelay->delay(1.0, slewFlavor_);
0486 }
0487
0488 if (pulseShapeId != currentPulseShapeId_) {
0489 resetPulseShapeTemplate(pulseShapeId, ps, nSamples);
0490 }
0491
0492
0493 nnlsWork_.dt = hasTimeInfo ? timeSigmaSiPM_ : timeSigmaHPD_;
0494
0495 if (nnlsWork_.tsSize != nSamples) {
0496 nnlsWork_.tsSize = nSamples;
0497 nnlsWork_.amplitudes.resize(nSamples);
0498 nnlsWork_.noiseTerms.resize(nSamples);
0499 nnlsWork_.pedVals.resize(nSamples);
0500 }
0501 }
0502
0503 void MahiFit::resetPulseShapeTemplate(const int pulseShapeId,
0504 const HcalPulseShapes& ps,
0505 const unsigned int ) {
0506 ++cntsetPulseShape_;
0507
0508 psfPtr_ = nullptr;
0509 for (auto& elem : knownPulseShapes_) {
0510 if (elem.first == pulseShapeId) {
0511 psfPtr_ = &*elem.second;
0512 break;
0513 }
0514 }
0515
0516 if (!psfPtr_) {
0517
0518
0519 auto p = std::make_shared<FitterFuncs::PulseShapeFunctor>(
0520 ps.getShape(pulseShapeId), false, false, false, 1, 0, 0, hcal::constants::maxSamples);
0521 knownPulseShapes_.emplace_back(pulseShapeId, p);
0522 psfPtr_ = &*p;
0523 }
0524
0525 currentPulseShapeId_ = pulseShapeId;
0526 }
0527
0528 void MahiFit::nnlsUnconstrainParameter(Index idxp) const {
0529 if (idxp != nnlsWork_.nP) {
0530 nnlsWork_.aTaMat.col(nnlsWork_.nP).swap(nnlsWork_.aTaMat.col(idxp));
0531 nnlsWork_.aTaMat.row(nnlsWork_.nP).swap(nnlsWork_.aTaMat.row(idxp));
0532 nnlsWork_.pulseMat.col(nnlsWork_.nP).swap(nnlsWork_.pulseMat.col(idxp));
0533 Eigen::numext::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP), nnlsWork_.aTbVec.coeffRef(idxp));
0534 Eigen::numext::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP), nnlsWork_.ampVec.coeffRef(idxp));
0535 Eigen::numext::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP), nnlsWork_.bxs.coeffRef(idxp));
0536 }
0537 ++nnlsWork_.nP;
0538 }
0539
0540 void MahiFit::nnlsConstrainParameter(Index minratioidx) const {
0541 if (minratioidx != (nnlsWork_.nP - 1)) {
0542 nnlsWork_.aTaMat.col(nnlsWork_.nP - 1).swap(nnlsWork_.aTaMat.col(minratioidx));
0543 nnlsWork_.aTaMat.row(nnlsWork_.nP - 1).swap(nnlsWork_.aTaMat.row(minratioidx));
0544 nnlsWork_.pulseMat.col(nnlsWork_.nP - 1).swap(nnlsWork_.pulseMat.col(minratioidx));
0545 Eigen::numext::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP - 1), nnlsWork_.aTbVec.coeffRef(minratioidx));
0546 Eigen::numext::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP - 1), nnlsWork_.ampVec.coeffRef(minratioidx));
0547 Eigen::numext::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP - 1), nnlsWork_.bxs.coeffRef(minratioidx));
0548 }
0549 --nnlsWork_.nP;
0550 }
0551
0552 void MahiFit::phase1Debug(const HBHEChannelInfo& channelData, MahiDebugInfo& mdi) const {
0553 float recoEnergy, soiPlus1Energy, recoTime, chi2;
0554 bool use3;
0555 phase1Apply(channelData, recoEnergy, soiPlus1Energy, recoTime, use3, chi2);
0556
0557 mdi.nSamples = channelData.nSamples();
0558 mdi.soi = channelData.soi();
0559
0560 mdi.use3 = use3;
0561
0562 mdi.inTimeConst = nnlsWork_.dt;
0563 mdi.inPedAvg = 0.25 * (channelData.tsPedestalWidth(0) * channelData.tsPedestalWidth(0) +
0564 channelData.tsPedestalWidth(1) * channelData.tsPedestalWidth(1) +
0565 channelData.tsPedestalWidth(2) * channelData.tsPedestalWidth(2) +
0566 channelData.tsPedestalWidth(3) * channelData.tsPedestalWidth(3));
0567 mdi.inGain = channelData.tsGain(0);
0568
0569 for (unsigned int iTS = 0; iTS < channelData.nSamples(); ++iTS) {
0570 double charge = channelData.tsRawCharge(iTS);
0571 double ped = channelData.tsPedestal(iTS);
0572
0573 mdi.inNoiseADC[iTS] = norm_ * channelData.tsDFcPerADC(iTS);
0574
0575 if ((charge - ped) > channelData.tsPedestalWidth(iTS)) {
0576 mdi.inNoisePhoto[iTS] = sqrt((charge - ped) * channelData.fcByPE());
0577 } else {
0578 mdi.inNoisePhoto[iTS] = 0;
0579 }
0580
0581 mdi.inPedestal[iTS] = channelData.tsPedestalWidth(iTS);
0582 mdi.totalUCNoise[iTS] = nnlsWork_.noiseTerms.coeffRef(iTS);
0583
0584 if (channelData.hasTimeInfo()) {
0585 mdi.inputTDC[iTS] = channelData.tsRiseTime(iTS);
0586 } else {
0587 mdi.inputTDC[iTS] = -1;
0588 }
0589 }
0590
0591 mdi.arrivalTime = recoTime;
0592 mdi.chiSq = chi2;
0593
0594 for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0595 if (nnlsWork_.bxs.coeff(iBX) == 0) {
0596 mdi.mahiEnergy = nnlsWork_.ampVec.coeff(iBX);
0597 for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
0598 mdi.count[iTS] = iTS;
0599 mdi.inputTS[iTS] = nnlsWork_.amplitudes.coeff(iTS);
0600 mdi.itPulse[iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
0601 }
0602 } else if (nnlsWork_.bxs.coeff(iBX) == pedestalBX_) {
0603 mdi.pedEnergy = nnlsWork_.ampVec.coeff(iBX);
0604 } else if (nnlsWork_.bxs.coeff(iBX) >= -3 && nnlsWork_.bxs.coeff(iBX) <= 4) {
0605 int ootIndex = nnlsWork_.bxs.coeff(iBX);
0606 if (ootIndex > 0)
0607 ootIndex += 2;
0608 else
0609 ootIndex += 3;
0610 mdi.ootEnergy[ootIndex] = nnlsWork_.ampVec.coeff(iBX);
0611 for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
0612 mdi.ootPulse[ootIndex][iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
0613 }
0614 }
0615 }
0616 }
0617
0618 void MahiFit::solveSubmatrix(PulseMatrix& mat, PulseVector& invec, PulseVector& outvec, unsigned nP) const {
0619 using namespace Eigen;
0620 switch (nP) {
0621 case 10: {
0622 Matrix<float, 10, 10> temp = mat;
0623 outvec.head<10>() = temp.ldlt().solve(invec.head<10>());
0624 } break;
0625 case 9: {
0626 Matrix<float, 9, 9> temp = mat.topLeftCorner<9, 9>();
0627 outvec.head<9>() = temp.ldlt().solve(invec.head<9>());
0628 } break;
0629 case 8: {
0630 Matrix<float, 8, 8> temp = mat.topLeftCorner<8, 8>();
0631 outvec.head<8>() = temp.ldlt().solve(invec.head<8>());
0632 } break;
0633 case 7: {
0634 Matrix<float, 7, 7> temp = mat.topLeftCorner<7, 7>();
0635 outvec.head<7>() = temp.ldlt().solve(invec.head<7>());
0636 } break;
0637 case 6: {
0638 Matrix<float, 6, 6> temp = mat.topLeftCorner<6, 6>();
0639 outvec.head<6>() = temp.ldlt().solve(invec.head<6>());
0640 } break;
0641 case 5: {
0642 Matrix<float, 5, 5> temp = mat.topLeftCorner<5, 5>();
0643 outvec.head<5>() = temp.ldlt().solve(invec.head<5>());
0644 } break;
0645 case 4: {
0646 Matrix<float, 4, 4> temp = mat.topLeftCorner<4, 4>();
0647 outvec.head<4>() = temp.ldlt().solve(invec.head<4>());
0648 } break;
0649 case 3: {
0650 Matrix<float, 3, 3> temp = mat.topLeftCorner<3, 3>();
0651 outvec.head<3>() = temp.ldlt().solve(invec.head<3>());
0652 } break;
0653 case 2: {
0654 Matrix<float, 2, 2> temp = mat.topLeftCorner<2, 2>();
0655 outvec.head<2>() = temp.ldlt().solve(invec.head<2>());
0656 } break;
0657 case 1: {
0658 Matrix<float, 1, 1> temp = mat.topLeftCorner<1, 1>();
0659 outvec.head<1>() = temp.ldlt().solve(invec.head<1>());
0660 } break;
0661 default:
0662 throw cms::Exception("HcalMahiWeirdState")
0663 << "Weird number of pulses encountered in Mahi, module is configured incorrectly!";
0664 }
0665 }
0666
0667 void MahiFit::resetWorkspace() const {
0668 nnlsWork_.nPulseTot = 0;
0669 nnlsWork_.tsOffset = 0;
0670 nnlsWork_.bxOffset = 0;
0671 nnlsWork_.maxoffset = 0;
0672 nnlsWork_.nP = 0;
0673
0674
0675 nnlsWork_.amplitudes.setZero();
0676 nnlsWork_.noiseTerms.setZero();
0677 nnlsWork_.pedVals.setZero();
0678 }