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                             int timeAlgo,
0014                             double iThEnergeticPulses,
0015                             double iMeanTime,
0016                             double iTimeSigmaHPD,
0017                             double iTimeSigmaSiPM,
0018                             const std::vector<int>& iActiveBXs,
0019                             int iNMaxItersMin,
0020                             int iNMaxItersNNLS,
0021                             double iDeltaChiSqThresh,
0022                             double iNnlsThresh) {
0023   dynamicPed_ = iDynamicPed;
0024
0025   ts4Thresh_ = iTS4Thresh;
0026   chiSqSwitch_ = chiSqSwitch;
0027
0028   applyTimeSlew_ = iApplyTimeSlew;
0029   slewFlavor_ = slewFlavor;
0030
0031   calculateArrivalTime_ = iCalculateArrivalTime;
0032   timeAlgo_ = timeAlgo;
0033   thEnergeticPulses_ = iThEnergeticPulses;
0034   meanTime_ = iMeanTime;
0035   timeSigmaHPD_ = iTimeSigmaHPD;
0036   timeSigmaSiPM_ = iTimeSigmaSiPM;
0037
0038   activeBXs_ = iActiveBXs;
0039
0040   nMaxItersMin_ = iNMaxItersMin;
0041   nMaxItersNNLS_ = iNMaxItersNNLS;
0042
0043   deltaChiSqThresh_ = iDeltaChiSqThresh;
0044   nnlsThresh_ = iNnlsThresh;
0045
0046   bxOffsetConf_ = -(*std::min_element(activeBXs_.begin(), activeBXs_.end()));
0047   bxSizeConf_ = activeBXs_.size();
0048 }
0049
0050 void MahiFit::phase1Apply(const HBHEChannelInfo& channelData,
0051                           float& reconstructedEnergy,
0052                           float& soiPlusOneEnergy,
0053                           float& reconstructedTime,
0054                           bool& useTriple,
0055                           float& chi2) const {
0056   const unsigned nSamples = channelData.nSamples();
0057   const unsigned soi = channelData.soi();
0058
0059   // Check some of the assumptions used by the subsequent code
0060   assert(nSamples == 8 || nSamples == 10);
0061   assert(nnlsWork_.tsSize <= nSamples);
0062   assert(soi + 1U < nnlsWork_.tsSize);
0063
0064   resetWorkspace();
0065
0066   nnlsWork_.tsOffset = soi;
0067
0068   // Packing energy, time, chiSq, soiPlus1Energy, in that order
0069   std::array<float, 4> reconstructedVals{{0.f, -9999.f, -9999.f, 0.f}};
0070
0071   double tsTOT = 0, tstrig = 0;  // in GeV
0072   for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
0073     auto const amplitude = channelData.tsRawCharge(iTS) - channelData.tsPedestal(iTS);
0074
0075     nnlsWork_.amplitudes.coeffRef(iTS) = amplitude;
0076
0077     //ADC granularity
0078     auto const noiseADC = norm_ * channelData.tsDFcPerADC(iTS);
0079
0080     //Electronic pedestal
0081     auto const pedWidth = channelData.tsPedestalWidth(iTS);
0082
0083     //Photostatistics
0084     auto const noisePhotoSq = (amplitude > pedWidth) ? (amplitude * channelData.fcByPE()) : 0.f;
0085
0086     //Total uncertainty from all sources
0087     nnlsWork_.noiseTerms.coeffRef(iTS) = noiseADC * noiseADC + noisePhotoSq + pedWidth * pedWidth;
0088
0089     nnlsWork_.pedVals.coeffRef(iTS) = pedWidth;
0090
0091     tsTOT += amplitude;
0092     if (iTS == soi)
0093       tstrig += amplitude;
0094   }
0095
0096   // This preserves the original Mahi approach but will have to
0097   // change in case we eventually calibrate gains per capId
0098   const float gain0 = channelData.tsGain(0);
0099   tsTOT *= gain0;
0100   tstrig *= gain0;
0101
0102   useTriple = false;
0103   if (tstrig > ts4Thresh_ && tsTOT > 0) {
0104     nnlsWork_.noisecorr = channelData.noisecorr();
0105
0106     if (nnlsWork_.noisecorr != 0) {
0107       nnlsWork_.pedVal = 0.f;
0108     } else {
0109       //Average pedestal width (for covariance matrix constraint)
0110       nnlsWork_.pedVal = 0.25f * (channelData.tsPedestalWidth(0) * channelData.tsPedestalWidth(0) +
0111                                   channelData.tsPedestalWidth(1) * channelData.tsPedestalWidth(1) +
0112                                   channelData.tsPedestalWidth(2) * channelData.tsPedestalWidth(2) +
0113                                   channelData.tsPedestalWidth(3) * channelData.tsPedestalWidth(3));
0114     }
0115
0116     // only do pre-fit with 1 pulse if chiSq threshold is positive
0117     if (chiSqSwitch_ > 0) {
0118       doFit(reconstructedVals, 1);
0119       if (reconstructedVals[2] > chiSqSwitch_) {
0120         doFit(reconstructedVals, 0);  //nbx=0 means use configured BXs
0121         useTriple = true;
0122       }
0123     } else {
0124       doFit(reconstructedVals, 0);
0125       useTriple = true;
0126     }
0127   }
0128
0129   reconstructedEnergy = reconstructedVals[0] * gain0;
0130   soiPlusOneEnergy = reconstructedVals[3] * gain0;
0131   reconstructedTime = reconstructedVals[1];
0132   chi2 = reconstructedVals[2];
0133 }
0134
0135 void MahiFit::doFit(std::array<float, 4>& correctedOutput, int nbx) const {
0136   unsigned int bxSize = 1;
0137
0138   if (nbx == 1) {
0139     nnlsWork_.bxOffset = 0;
0140   } else {
0141     bxSize = bxSizeConf_;
0142     nnlsWork_.bxOffset = static_cast<int>(nnlsWork_.tsOffset) >= bxOffsetConf_ ? bxOffsetConf_ : nnlsWork_.tsOffset;
0143   }
0144
0145   nnlsWork_.nPulseTot = bxSize;
0146
0147   if (dynamicPed_)
0148     nnlsWork_.nPulseTot++;
0149   nnlsWork_.bxs.setZero(nnlsWork_.nPulseTot);
0150
0151   if (nbx == 1) {
0152     nnlsWork_.bxs.coeffRef(0) = 0;
0153   } else {
0154     for (unsigned int iBX = 0; iBX < bxSize; ++iBX) {
0155       nnlsWork_.bxs.coeffRef(iBX) =
0156           activeBXs_[iBX] -
0157           ((static_cast<int>(nnlsWork_.tsOffset) + activeBXs_[0]) >= 0 ? 0 : (nnlsWork_.tsOffset + activeBXs_[0]));
0158     }
0159   }
0160
0161   nnlsWork_.maxoffset = nnlsWork_.bxs.coeff(bxSize - 1);
0162   if (dynamicPed_)
0163     nnlsWork_.bxs[nnlsWork_.nPulseTot - 1] = pedestalBX_;
0164
0165   nnlsWork_.pulseMat.setZero(nnlsWork_.tsSize, nnlsWork_.nPulseTot);
0166   nnlsWork_.pulseDerivMat.setZero(nnlsWork_.tsSize, nnlsWork_.nPulseTot);
0167
0168   FullSampleVector pulseShapeArray;
0169   FullSampleVector pulseDerivArray;
0170   FullSampleMatrix pulseCov;
0171
0172   int offset = 0;
0173   for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0174     offset = nnlsWork_.bxs.coeff(iBX);
0175
0176     if (offset == pedestalBX_) {
0177       nnlsWork_.pulseMat.col(iBX) = SampleVector::Ones(nnlsWork_.tsSize);
0178       nnlsWork_.pulseDerivMat.col(iBX) = SampleVector::Zero(nnlsWork_.tsSize);
0179     } else {
0180       pulseShapeArray.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
0181       if (calculateArrivalTime_)
0182         pulseDerivArray.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
0183       pulseCov.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset,
0184                        nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
0185       nnlsWork_.pulseCovArray[iBX].setZero(nnlsWork_.tsSize, nnlsWork_.tsSize);
0186
0187       updatePulseShape(
0188           nnlsWork_.amplitudes.coeff(nnlsWork_.tsOffset + offset), pulseShapeArray, pulseDerivArray, pulseCov);
0189
0190       nnlsWork_.pulseMat.col(iBX) = pulseShapeArray.segment(nnlsWork_.maxoffset - offset, nnlsWork_.tsSize);
0191       if (calculateArrivalTime_)
0192         nnlsWork_.pulseDerivMat.col(iBX) = pulseDerivArray.segment(nnlsWork_.maxoffset - offset, nnlsWork_.tsSize);
0193       nnlsWork_.pulseCovArray[iBX] = pulseCov.block(
0194           nnlsWork_.maxoffset - offset, nnlsWork_.maxoffset - offset, nnlsWork_.tsSize, nnlsWork_.tsSize);
0195     }
0196   }
0197
0198   const float chiSq = minimize();
0199
0200   bool foundintime = false;
0201   unsigned int ipulseintime = 0;
0202
0203   for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0204     if (nnlsWork_.bxs.coeff(iBX) == 0) {
0205       ipulseintime = iBX;
0206       foundintime = true;
0207       break;
0208     }
0209   }
0210
0211   if (foundintime) {
0212     correctedOutput.at(0) = nnlsWork_.ampVec.coeff(ipulseintime);       //charge
0213     correctedOutput.at(3) = nnlsWork_.ampVec.coeff(ipulseintime + 1U);  //charge for SOI+1
0214     if (correctedOutput.at(0) != 0) {
0215       // fixME store the timeslew
0216       float arrivalTime = 0.f;
0217       if (calculateArrivalTime_ && timeAlgo_ == 1)
0218         arrivalTime = calculateArrivalTime(ipulseintime);
0219       else if (calculateArrivalTime_ && timeAlgo_ == 2)
0220         arrivalTime = ccTime(nnlsWork_.ampVec.coeff(ipulseintime));
0221       correctedOutput.at(1) = arrivalTime;  //time
0222     } else
0223       correctedOutput.at(1) = -9999.f;  //time
0224
0225     correctedOutput.at(2) = chiSq;  //chi2
0226   }
0227 }
0228
0229 const float MahiFit::minimize() const {
0230   nnlsWork_.invcovp.setZero(nnlsWork_.tsSize, nnlsWork_.nPulseTot);
0231   nnlsWork_.ampVec.setZero(nnlsWork_.nPulseTot);
0232
0233   SampleMatrix invCovMat;
0234   invCovMat.setConstant(nnlsWork_.tsSize, nnlsWork_.tsSize, nnlsWork_.pedVal);
0235   invCovMat += nnlsWork_.noiseTerms.asDiagonal();
0236
0237   //Add off-Diagonal components up to first order
0238   if (nnlsWork_.noisecorr != 0) {
0239     for (unsigned int i = 1; i < nnlsWork_.tsSize; ++i) {
0240       auto const noiseCorrTerm = nnlsWork_.noisecorr * nnlsWork_.pedVals.coeff(i - 1) * nnlsWork_.pedVals.coeff(i);
0241       invCovMat(i - 1, i) += noiseCorrTerm;
0242       invCovMat(i, i - 1) += noiseCorrTerm;
0243     }
0244   }
0245
0246   float oldChiSq = 9999;
0247   float chiSq = oldChiSq;
0248
0249   for (int iter = 1; iter < nMaxItersMin_; ++iter) {
0250     updateCov(invCovMat);
0251
0252     if (nnlsWork_.nPulseTot > 1) {
0253       nnls();
0254     } else {
0255       onePulseMinimize();
0256     }
0257
0258     const float newChiSq = calculateChiSq();
0259     const float deltaChiSq = newChiSq - chiSq;
0260
0261     if (newChiSq == oldChiSq && newChiSq < chiSq) {
0262       break;
0263     }
0264     oldChiSq = chiSq;
0265     chiSq = newChiSq;
0266
0267     if (std::abs(deltaChiSq) < deltaChiSqThresh_)
0268       break;
0269   }
0270
0271   return chiSq;
0272 }
0273
0274 void MahiFit::updatePulseShape(const float itQ,
0275                                FullSampleVector& pulseShape,
0276                                FullSampleVector& pulseDeriv,
0277                                FullSampleMatrix& pulseCov) const {
0278   // set a null pulse shape for negative / or null TS
0279   if (itQ <= 0.f)
0280     return;
0281
0282   float t0 = meanTime_;
0283
0284   if (applyTimeSlew_) {
0285     if (itQ <= 1.f)
0286       t0 += tsDelay1GeV_;
0287     else
0288       t0 += hcalTimeSlewDelay_->delay(float(itQ), slewFlavor_);
0289   }
0290
0291   std::array<double, hcal::constants::maxSamples> pulseN;
0292   std::array<double, hcal::constants::maxSamples> pulseM;
0293   std::array<double, hcal::constants::maxSamples> pulseP;
0294
0295   const float xx = t0;
0296   const float xxm = -nnlsWork_.dt + t0;
0297   const float xxp = nnlsWork_.dt + t0;
0298
0299   psfPtr_->singlePulseShapeFuncMahi(&xx);
0300   psfPtr_->getPulseShape(pulseN);
0301
0302   psfPtr_->singlePulseShapeFuncMahi(&xxm);
0303   psfPtr_->getPulseShape(pulseM);
0304
0305   psfPtr_->singlePulseShapeFuncMahi(&xxp);
0306   psfPtr_->getPulseShape(pulseP);
0307
0308   //in the 2018+ case where the sample of interest (SOI) is in TS3, add an extra offset to align
0309   //with previous SOI=TS4 case assumed by psfPtr_->getPulseShape()
0310   int delta = 4 - nnlsWork_.tsOffset;
0311
0312   auto invDt = 0.5 / nnlsWork_.dt;
0313
0314   for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
0315     pulseShape(iTS + nnlsWork_.maxoffset) = pulseN[iTS + delta];
0316     if (calculateArrivalTime_)
0317       pulseDeriv(iTS + nnlsWork_.maxoffset) = (pulseM[iTS + delta] - pulseP[iTS + delta]) * invDt;
0318
0319     pulseM[iTS + delta] -= pulseN[iTS + delta];
0320     pulseP[iTS + delta] -= pulseN[iTS + delta];
0321   }
0322
0323   for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
0324     for (unsigned int jTS = 0; jTS < iTS + 1; ++jTS) {
0325       auto const tmp = 0.5 * (pulseP[iTS + delta] * pulseP[jTS + delta] + pulseM[iTS + delta] * pulseM[jTS + delta]);
0326
0327       pulseCov(jTS + nnlsWork_.maxoffset, iTS + nnlsWork_.maxoffset) = tmp;
0328       if (iTS != jTS)
0329         pulseCov(iTS + nnlsWork_.maxoffset, jTS + nnlsWork_.maxoffset) = tmp;
0330     }
0331   }
0332 }
0333
0334 void MahiFit::updateCov(const SampleMatrix& samplecov) const {
0335   SampleMatrix invCovMat = samplecov;
0336
0337   for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0338     auto const amp = nnlsWork_.ampVec.coeff(iBX);
0339     if (amp == 0)
0340       continue;
0341
0342     int offset = nnlsWork_.bxs.coeff(iBX);
0343
0344     if (offset == pedestalBX_)
0345       continue;
0346     else {
0347       invCovMat.noalias() += amp * amp * nnlsWork_.pulseCovArray[offset + nnlsWork_.bxOffset];
0348     }
0349   }
0350
0351   nnlsWork_.covDecomp.compute(invCovMat);
0352 }
0353
0354 float MahiFit::ccTime(const float itQ) const {
0355   // those conditions are now on data time slices, can be done on the fitted pulse i.e. using nlsWork_.ampVec.coeff(itIndex);
0356
0357   // Selecting energetic hits - Fitted Energy > 20 GeV
0358   if (itQ < thEnergeticPulsesFC_)
0359     return HcalSpecialTimes::DEFAULT_ccTIME;
0360
0361   // Rejecting late hits  Energy in TS[3] > (Energy in TS[4] and TS[5])
0362   // With small OOTPU (Energy in TS[0] ,TS[1] and TS[2]) < 5 GeV
0363   // To speed up check around the fitted time (?)
0364
0365   // distanze as in formula of page 6
0366   // https://indico.cern.ch/event/1142347/contributions/4793749/attachments/2412936/4129323/HCAL%20timing%20update.pdf
0367
0368   float t0 = meanTime_;
0369
0370   if (applyTimeSlew_) {
0371     if (itQ <= 1.f)
0372       t0 += tsDelay1GeV_;
0373     else
0374       t0 += hcalTimeSlewDelay_->delay(float(itQ), slewFlavor_);
0375   }
0376
0377   float ccTime = 0.f;
0378   float distance_delta_max = 0.f;
0379
0380   std::array<double, hcal::constants::maxSamples> pulseN;
0381
0382   // making 8 TS out of the template i.e. 200 points
0383   int TS_SOIandAfter = 25 * (nnlsWork_.tsSize - nnlsWork_.tsOffset);
0384   int TS_beforeSOI = -25 * nnlsWork_.tsOffset;
0385
0386   for (int deltaNS = TS_beforeSOI; deltaNS < TS_SOIandAfter; ++deltaNS) {  // from -75ns and + 125ns
0387     const float xx = t0 + deltaNS;
0388
0389     psfPtr_->singlePulseShapeFuncMahi(&xx);
0390     psfPtr_->getPulseShape(pulseN);
0391
0392     float pulse2 = 0;
0393     float norm2 = 0;
0394     float numerator = 0;
0395     //
0396     int delta = 4 - nnlsWork_.tsOffset;  // like in updatePulseShape
0397
0398     // rm TS0 and TS7, to speed up and reduce noise
0399     for (unsigned int iTS = 1; iTS < (nnlsWork_.tsSize - 1); ++iTS) {
0400       //pulseN[iTS] is the area of the template
0401       float norm = nnlsWork_.amplitudes.coeffRef(iTS);
0402
0403       //  Finding the distance after each iteration.
0404       numerator += norm * pulseN[iTS + delta];
0405       pulse2 += pulseN[iTS + delta] * pulseN[iTS + delta];
0406       norm2 += norm * norm;
0407     }
0408
0409     float distance = numerator / sqrt(pulse2 * norm2);
0410     if (distance > distance_delta_max) {
0411       distance_delta_max = distance;
0412       ccTime = deltaNS;
0413     }
0414   }
0415
0416   return ccTime;
0417 }
0418
0419 float MahiFit::calculateArrivalTime(const unsigned int itIndex) const {
0420   if (nnlsWork_.nPulseTot > 1) {
0421     SamplePulseMatrix pulseDerivMatTMP = nnlsWork_.pulseDerivMat;
0422     for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0423       nnlsWork_.pulseDerivMat.col(iBX) = pulseDerivMatTMP.col(nnlsWork_.bxs.coeff(iBX) + nnlsWork_.bxOffset);
0424     }
0425   }
0426
0427   for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0428     nnlsWork_.pulseDerivMat.col(iBX) *= nnlsWork_.ampVec.coeff(iBX);
0429   }
0430
0431   //FIXME: avoid solve full equation for one element
0432   SampleVector residuals = nnlsWork_.pulseMat * nnlsWork_.ampVec - nnlsWork_.amplitudes;
0433   PulseVector solution = nnlsWork_.pulseDerivMat.colPivHouseholderQr().solve(residuals);
0434   float t = std::clamp((float)solution.coeff(itIndex), -timeLimit_, timeLimit_);
0435
0436   return t;
0437 }
0438
0439 void MahiFit::nnls() const {
0440   const unsigned int npulse = nnlsWork_.nPulseTot;
0441   const unsigned int nsamples = nnlsWork_.tsSize;
0442
0443   PulseVector updateWork;
0444   PulseVector ampvecpermtest;
0445
0446   nnlsWork_.invcovp = nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat);
0447   nnlsWork_.aTaMat.noalias() = nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.invcovp);
0448   nnlsWork_.aTbVec.noalias() =
0449       nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes));
0450
0451   int iter = 0;
0452   Index idxwmax = 0;
0453   float wmax = 0.0f;
0454   float threshold = nnlsThresh_;
0455
0456   while (true) {
0457     if (iter > 0 || nnlsWork_.nP == 0) {
0458       if (nnlsWork_.nP == std::min(npulse, nsamples))
0459         break;
0460
0461       const unsigned int nActive = npulse - nnlsWork_.nP;
0462       // exit if there are no more pulses to constrain
0463       if (nActive == 0)
0464         break;
0465
0466       updateWork.noalias() = nnlsWork_.aTbVec - nnlsWork_.aTaMat.lazyProduct(nnlsWork_.ampVec);
0467
0468       Index idxwmaxprev = idxwmax;
0469       float wmaxprev = wmax;
0470       wmax = updateWork.tail(nActive).maxCoeff(&idxwmax);
0471
0472       if (wmax < threshold || (idxwmax == idxwmaxprev && wmax == wmaxprev)) {
0473         break;
0474       }
0475
0476       if (iter >= nMaxItersNNLS_) {
0477         break;
0478       }
0479
0480       //unconstrain parameter
0481       idxwmax += nnlsWork_.nP;
0482       nnlsUnconstrainParameter(idxwmax);
0483     }
0484
0485     while (true) {
0486       if (nnlsWork_.nP == 0)
0487         break;
0488
0489       ampvecpermtest.noalias() = nnlsWork_.ampVec.head(nnlsWork_.nP);
0490
0491       solveSubmatrix(nnlsWork_.aTaMat, nnlsWork_.aTbVec, ampvecpermtest, nnlsWork_.nP);
0492
0493       //check solution
0494       if (ampvecpermtest.head(nnlsWork_.nP).minCoeff() > 0.f) {
0495         nnlsWork_.ampVec.head(nnlsWork_.nP) = ampvecpermtest.head(nnlsWork_.nP);
0496         break;
0497       }
0498
0499       //update parameter vector
0500       Index minratioidx = 0;
0501
0502       // no realizable optimization here (because it autovectorizes!)
0503       float minratio = std::numeric_limits<float>::max();
0504       for (unsigned int ipulse = 0; ipulse < nnlsWork_.nP; ++ipulse) {
0505         if (ampvecpermtest.coeff(ipulse) <= 0.f) {
0506           const float c_ampvec = nnlsWork_.ampVec.coeff(ipulse);
0507           const float ratio = c_ampvec / (c_ampvec - ampvecpermtest.coeff(ipulse));
0508           if (ratio < minratio) {
0509             minratio = ratio;
0510             minratioidx = ipulse;
0511           }
0512         }
0513       }
0514       nnlsWork_.ampVec.head(nnlsWork_.nP) +=
0515           minratio * (ampvecpermtest.head(nnlsWork_.nP) - nnlsWork_.ampVec.head(nnlsWork_.nP));
0516
0517       //avoid numerical problems with later ==0. check
0518       nnlsWork_.ampVec.coeffRef(minratioidx) = 0.f;
0519
0520       nnlsConstrainParameter(minratioidx);
0521     }
0522
0523     ++iter;
0524
0525     //adaptive convergence threshold to avoid infinite loops but still
0526     //ensure best value is used
0527     if (iter % 10 == 0) {
0528       threshold *= 10.;
0529     }
0530   }
0531 }
0532
0533 void MahiFit::onePulseMinimize() const {
0534   nnlsWork_.invcovp = nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat);
0535
0536   float aTaCoeff = (nnlsWork_.invcovp.transpose() * nnlsWork_.invcovp).coeff(0);
0537   float aTbCoeff =
0538       (nnlsWork_.invcovp.transpose() * (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes))).coeff(0);
0539
0540   nnlsWork_.ampVec.coeffRef(0) = std::max(0.f, aTbCoeff / aTaCoeff);
0541 }
0542
0543 float MahiFit::calculateChiSq() const {
0544   return (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat * nnlsWork_.ampVec - nnlsWork_.amplitudes))
0545       .squaredNorm();
0546 }
0547
0548 void MahiFit::setPulseShapeTemplate(const int pulseShapeId,
0549                                     const HcalPulseShapes& ps,
0550                                     const bool hasTimeInfo,
0551                                     const HcalTimeSlew* hcalTimeSlewDelay,
0552                                     const unsigned int nSamples,
0553                                     const float gain0) {
0554   if (hcalTimeSlewDelay != hcalTimeSlewDelay_) {
0555     assert(hcalTimeSlewDelay);
0556     hcalTimeSlewDelay_ = hcalTimeSlewDelay;
0557     tsDelay1GeV_ = hcalTimeSlewDelay->delay(1.0, slewFlavor_);
0558   }
0559
0560   if (pulseShapeId != currentPulseShapeId_) {
0561     resetPulseShapeTemplate(pulseShapeId, ps, nSamples);
0562   }
0563
0564   // 1 sigma time constraint
0565   nnlsWork_.dt = hasTimeInfo ? timeSigmaSiPM_ : timeSigmaHPD_;
0566
0567   if (nnlsWork_.tsSize != nSamples) {
0568     nnlsWork_.tsSize = nSamples;
0569     nnlsWork_.amplitudes.resize(nSamples);
0570     nnlsWork_.noiseTerms.resize(nSamples);
0571     nnlsWork_.pedVals.resize(nSamples);
0572   }
0573
0574   // threshold in GeV for ccTime
0575   thEnergeticPulsesFC_ = thEnergeticPulses_ / gain0;
0576 }
0577
0578 void MahiFit::resetPulseShapeTemplate(const int pulseShapeId,
0579                                       const HcalPulseShapes& ps,
0580                                       const unsigned int /* nSamples */) {
0581   ++cntsetPulseShape_;
0582
0583   psfPtr_ = nullptr;
0584   for (auto& elem : knownPulseShapes_) {
0585     if (elem.first == pulseShapeId) {
0586       psfPtr_ = &*elem.second;
0587       break;
0588     }
0589   }
0590
0591   if (!psfPtr_) {
0592     // only the pulse shape itself from PulseShapeFunctor is used for Mahi
0593     // the uncertainty terms calculated inside PulseShapeFunctor are used for Method 2 only
0594     auto p = std::make_shared<FitterFuncs::PulseShapeFunctor>(
0595         ps.getShape(pulseShapeId), false, false, false, 1, 0, 0, hcal::constants::maxSamples);
0596     knownPulseShapes_.emplace_back(pulseShapeId, p);
0597     psfPtr_ = &*p;
0598   }
0599
0600   currentPulseShapeId_ = pulseShapeId;
0601 }
0602
0603 void MahiFit::nnlsUnconstrainParameter(Index idxp) const {
0604   if (idxp != nnlsWork_.nP) {
0605     nnlsWork_.aTaMat.col(nnlsWork_.nP).swap(nnlsWork_.aTaMat.col(idxp));
0606     nnlsWork_.aTaMat.row(nnlsWork_.nP).swap(nnlsWork_.aTaMat.row(idxp));
0607     nnlsWork_.pulseMat.col(nnlsWork_.nP).swap(nnlsWork_.pulseMat.col(idxp));
0608     Eigen::numext::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP), nnlsWork_.aTbVec.coeffRef(idxp));
0609     Eigen::numext::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP), nnlsWork_.ampVec.coeffRef(idxp));
0610     Eigen::numext::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP), nnlsWork_.bxs.coeffRef(idxp));
0611   }
0612   ++nnlsWork_.nP;
0613 }
0614
0615 void MahiFit::nnlsConstrainParameter(Index minratioidx) const {
0616   if (minratioidx != (nnlsWork_.nP - 1)) {
0617     nnlsWork_.aTaMat.col(nnlsWork_.nP - 1).swap(nnlsWork_.aTaMat.col(minratioidx));
0618     nnlsWork_.aTaMat.row(nnlsWork_.nP - 1).swap(nnlsWork_.aTaMat.row(minratioidx));
0619     nnlsWork_.pulseMat.col(nnlsWork_.nP - 1).swap(nnlsWork_.pulseMat.col(minratioidx));
0620     Eigen::numext::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP - 1), nnlsWork_.aTbVec.coeffRef(minratioidx));
0621     Eigen::numext::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP - 1), nnlsWork_.ampVec.coeffRef(minratioidx));
0622     Eigen::numext::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP - 1), nnlsWork_.bxs.coeffRef(minratioidx));
0623   }
0624   --nnlsWork_.nP;
0625 }
0626
0627 void MahiFit::phase1Debug(const HBHEChannelInfo& channelData, MahiDebugInfo& mdi) const {
0628   float recoEnergy, soiPlus1Energy, recoTime, chi2;
0629   bool use3;
0630   phase1Apply(channelData, recoEnergy, soiPlus1Energy, recoTime, use3, chi2);
0631
0632   mdi.nSamples = channelData.nSamples();
0633   mdi.soi = channelData.soi();
0634
0635   mdi.use3 = use3;
0636
0637   mdi.inTimeConst = nnlsWork_.dt;
0638   mdi.inPedAvg = 0.25 * (channelData.tsPedestalWidth(0) * channelData.tsPedestalWidth(0) +
0639                          channelData.tsPedestalWidth(1) * channelData.tsPedestalWidth(1) +
0640                          channelData.tsPedestalWidth(2) * channelData.tsPedestalWidth(2) +
0641                          channelData.tsPedestalWidth(3) * channelData.tsPedestalWidth(3));
0642   mdi.inGain = channelData.tsGain(0);
0643
0644   for (unsigned int iTS = 0; iTS < channelData.nSamples(); ++iTS) {
0645     double charge = channelData.tsRawCharge(iTS);
0646     double ped = channelData.tsPedestal(iTS);
0647
0648     mdi.inNoiseADC[iTS] = norm_ * channelData.tsDFcPerADC(iTS);
0649
0650     if ((charge - ped) > channelData.tsPedestalWidth(iTS)) {
0651       mdi.inNoisePhoto[iTS] = sqrt((charge - ped) * channelData.fcByPE());
0652     } else {
0653       mdi.inNoisePhoto[iTS] = 0;
0654     }
0655
0656     mdi.inPedestal[iTS] = channelData.tsPedestalWidth(iTS);
0657     mdi.totalUCNoise[iTS] = nnlsWork_.noiseTerms.coeffRef(iTS);
0658
0659     if (channelData.hasTimeInfo()) {
0660       mdi.inputTDC[iTS] = channelData.tsRiseTime(iTS);
0661     } else {
0662       mdi.inputTDC[iTS] = -1;
0663     }
0664   }
0665
0666   mdi.arrivalTime = recoTime;
0667   mdi.chiSq = chi2;
0668
0669   for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
0670     if (nnlsWork_.bxs.coeff(iBX) == 0) {
0671       mdi.mahiEnergy = nnlsWork_.ampVec.coeff(iBX);
0672       for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
0673         mdi.count[iTS] = iTS;
0674         mdi.inputTS[iTS] = nnlsWork_.amplitudes.coeff(iTS);
0675         mdi.itPulse[iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
0676       }
0677     } else if (nnlsWork_.bxs.coeff(iBX) == pedestalBX_) {
0678       mdi.pedEnergy = nnlsWork_.ampVec.coeff(iBX);
0679     } else if (nnlsWork_.bxs.coeff(iBX) >= -3 && nnlsWork_.bxs.coeff(iBX) <= 4) {
0680       int ootIndex = nnlsWork_.bxs.coeff(iBX);
0681       if (ootIndex > 0)
0682         ootIndex += 2;
0683       else
0684         ootIndex += 3;
0685       mdi.ootEnergy[ootIndex] = nnlsWork_.ampVec.coeff(iBX);
0686       for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
0687         mdi.ootPulse[ootIndex][iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
0688       }
0689     }
0690   }
0691 }
0692
0693 void MahiFit::solveSubmatrix(PulseMatrix& mat, PulseVector& invec, PulseVector& outvec, unsigned nP) const {
0694   using namespace Eigen;
0695   switch (nP) {  // pulse matrix is always square.
0696     case 10: {
0697       Matrix<float, 10, 10> temp = mat;
0698       outvec.head<10>() = temp.ldlt().solve(invec.head<10>());
0699     } break;
0700     case 9: {
0701       Matrix<float, 9, 9> temp = mat.topLeftCorner<9, 9>();
0702       outvec.head<9>() = temp.ldlt().solve(invec.head<9>());
0703     } break;
0704     case 8: {
0705       Matrix<float, 8, 8> temp = mat.topLeftCorner<8, 8>();
0706       outvec.head<8>() = temp.ldlt().solve(invec.head<8>());
0707     } break;
0708     case 7: {
0709       Matrix<float, 7, 7> temp = mat.topLeftCorner<7, 7>();
0710       outvec.head<7>() = temp.ldlt().solve(invec.head<7>());
0711     } break;
0712     case 6: {
0713       Matrix<float, 6, 6> temp = mat.topLeftCorner<6, 6>();
0714       outvec.head<6>() = temp.ldlt().solve(invec.head<6>());
0715     } break;
0716     case 5: {
0717       Matrix<float, 5, 5> temp = mat.topLeftCorner<5, 5>();
0718       outvec.head<5>() = temp.ldlt().solve(invec.head<5>());
0719     } break;
0720     case 4: {
0721       Matrix<float, 4, 4> temp = mat.topLeftCorner<4, 4>();
0722       outvec.head<4>() = temp.ldlt().solve(invec.head<4>());
0723     } break;
0724     case 3: {
0725       Matrix<float, 3, 3> temp = mat.topLeftCorner<3, 3>();
0726       outvec.head<3>() = temp.ldlt().solve(invec.head<3>());
0727     } break;
0728     case 2: {
0729       Matrix<float, 2, 2> temp = mat.topLeftCorner<2, 2>();
0730       outvec.head<2>() = temp.ldlt().solve(invec.head<2>());
0731     } break;
0732     case 1: {
0733       Matrix<float, 1, 1> temp = mat.topLeftCorner<1, 1>();
0734       outvec.head<1>() = temp.ldlt().solve(invec.head<1>());
0735     } break;
0736     default:
0737       throw cms::Exception("HcalMahiWeirdState")
0738           << "Weird number of pulses encountered in Mahi, module is configured incorrectly!";
0739   }
0740 }
0741
0742 void MahiFit::resetWorkspace() const {
0743   nnlsWork_.nPulseTot = 0;
0744   nnlsWork_.tsOffset = 0;
0745   nnlsWork_.bxOffset = 0;
0746   nnlsWork_.maxoffset = 0;
0747   nnlsWork_.nP = 0;
0748
0749   // NOT SURE THIS IS NEEDED
0750   nnlsWork_.amplitudes.setZero();
0751   nnlsWork_.noiseTerms.setZero();
0752   nnlsWork_.pedVals.setZero();
0753 }