HcalRecAlgos/src/MahiFit.cc

0001 #define EIGEN_NO_DEBUG  // kill throws in eigen code
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   // Check some of the assumptions used by the subsequent code
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   // Packing energy, time, chiSq, soiPlus1Energy, in that order
0065   std::array<float, 4> reconstructedVals{{0.f, -9999.f, -9999.f, 0.f}};
0066
0067   double tsTOT = 0, tstrig = 0;  // in GeV
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     //ADC granularity
0074     auto const noiseADC = norm_ * channelData.tsDFcPerADC(iTS);
0075
0076     //Electronic pedestal
0077     auto const pedWidth = channelData.tsPedestalWidth(iTS);
0078
0079     //Photostatistics
0080     auto const noisePhotoSq = (amplitude > pedWidth) ? (amplitude * channelData.fcByPE()) : 0.f;
0081
0082     //Total uncertainty from all sources
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   // This preserves the original Mahi approach but will have to
0093   // change in case we eventually calibrate gains per capId
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       //Average pedestal width (for covariance matrix constraint)
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     // only do pre-fit with 1 pulse if chiSq threshold is positive
0113     if (chiSqSwitch_ > 0) {
0114       doFit(reconstructedVals, 1);
0115       if (reconstructedVals[2] > chiSqSwitch_) {
0116         doFit(reconstructedVals, 0);  //nbx=0 means use configured BXs
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);       //charge
0209     correctedOutput.at(3) = nnlsWork_.ampVec.coeff(ipulseintime + 1U);  //charge for SOI+1
0210     if (correctedOutput.at(0) != 0) {
0211       // fixME store the timeslew
0212       float arrivalTime = 0.f;
0213       if (calculateArrivalTime_)
0214         arrivalTime = calculateArrivalTime(ipulseintime);
0215       correctedOutput.at(1) = arrivalTime;  //time
0216     } else
0217       correctedOutput.at(1) = -9999.f;  //time
0218
0219     correctedOutput.at(2) = chiSq;  //chi2
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   //Add off-Diagonal components up to first order
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   // set a null pulse shape for negative / or null TS
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   //in the 2018+ case where the sample of interest (SOI) is in TS3, add an extra offset to align
0303   //with previous SOI=TS4 case assumed by psfPtr_->getPulseShape()
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   //FIXME: avoid solve full equation for one element
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       // exit if there are no more pulses to constrain
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       //unconstrain parameter
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       //check solution
0423       if (ampvecpermtest.head(nnlsWork_.nP).minCoeff() > 0.f) {
0424         nnlsWork_.ampVec.head(nnlsWork_.nP) = ampvecpermtest.head(nnlsWork_.nP);
0425         break;
0426       }
0427
0428       //update parameter vector
0429       Index minratioidx = 0;
0430
0431       // no realizable optimization here (because it autovectorizes!)
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       //avoid numerical problems with later ==0. check
0447       nnlsWork_.ampVec.coeffRef(minratioidx) = 0.f;
0448
0449       nnlsConstrainParameter(minratioidx);
0450     }
0451
0452     ++iter;
0453
0454     //adaptive convergence threshold to avoid infinite loops but still
0455     //ensure best value is used
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   // 1 sigma time constraint
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 /* nSamples */) {
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     // only the pulse shape itself from PulseShapeFunctor is used for Mahi
0518     // the uncertainty terms calculated inside PulseShapeFunctor are used for Method 2 only
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) {  // pulse matrix is always square.
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   // NOT SURE THIS IS NEEDED
0675   nnlsWork_.amplitudes.setZero();
0676   nnlsWork_.noiseTerms.setZero();
0677   nnlsWork_.pedVals.setZero();
0678 }