#include <MahiFit.h>

Public Types
typedef BXVector::Index	Index

Public Member Functions
void	doFit (std::array< float, 3 > &correctedOutput, const int nbx) const

	MahiFit ()

void	phase1Apply (const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &reconstructedTime, bool &useTriple, float &chi2) const

void	phase1Debug (const HBHEChannelInfo &channelData, MahiDebugInfo &mdi) const

void	resetPulseShapeTemplate (const HcalPulseShapes::Shape &ps, bool hasTimeInfo, unsigned int nSamples)

void	setParameters (bool iDynamicPed, double iTS4Thresh, double chiSqSwitch, bool iApplyTimeSlew, HcalTimeSlew::BiasSetting slewFlavor, bool iCalculateArrivalTime, double iMeanTime, double iTimeSigmaHPD, double iTimeSigmaSiPM, const std::vector< int > &iActiveBXs, int iNMaxItersMin, int iNMaxItersNNLS, double iDeltaChiSqThresh, double iNnlsThresh)

void	setPulseShapeTemplate (const HcalPulseShapes::Shape &ps, bool hasTimeInfo, const HcalTimeSlew *hcalTimeSlewDelay, unsigned int nSamples)

	~MahiFit ()

Public Attributes
const HcalPulseShapes::Shape *	currentPulseShape_ = nullptr

const HcalTimeSlew *	hcalTimeSlewDelay_ = nullptr

Private Member Functions
float	calculateArrivalTime (unsigned int iBX) const

float	calculateChiSq () const

const float	minimize () const

void	nnls () const

void	nnlsConstrainParameter (Index minratioidx) const

void	nnlsUnconstrainParameter (Index idxp) const

void	onePulseMinimize () const

void	resetWorkspace () const

void	solveSubmatrix (PulseMatrix &mat, PulseVector &invec, PulseVector &outvec, unsigned nP) const

void	updateCov (const SampleMatrix &invCovMat) const

void	updatePulseShape (const float itQ, FullSampleVector &pulseShape, FullSampleVector &pulseDeriv, FullSampleMatrix &pulseCov) const

Private Attributes
std::vector< int >	activeBXs_

bool	applyTimeSlew_

int	bxOffsetConf_

unsigned int	bxSizeConf_

bool	calculateArrivalTime_

float	chiSqSwitch_

int	cntsetPulseShape_

float	deltaChiSqThresh_

bool	dynamicPed_

float	meanTime_

int	nMaxItersMin_

int	nMaxItersNNLS_

float	nnlsThresh_

MahiNnlsWorkspace	nnlsWork_

float	norm_ = (1.f / std::sqrt(12))

std::unique_ptr< ROOT::Math::Functor >	pfunctor_

std::unique_ptr< FitterFuncs::PulseShapeFunctor >	psfPtr_

HcalTimeSlew::BiasSetting	slewFlavor_

float	timeSigmaHPD_

float	timeSigmaSiPM_

float	ts4Thresh_

float	tsDelay1GeV_ = 0.f

Static Private Attributes
static constexpr int	pedestalBX_ = 100

static constexpr float	timeLimit_ = 12.5f

Detailed Description

Definition at line 89 of file MahiFit.h.

Member Typedef Documentation

◆ Index

typedef BXVector::Index MahiFit::Index

Definition at line 125 of file MahiFit.h.

Constructor & Destructor Documentation

◆ MahiFit()

MahiFit::MahiFit ( )

Definition at line 5 of file MahiFit.cc.

5 {}

◆ ~MahiFit()

MahiFit::~MahiFit ( )

inline

Definition at line 92 of file MahiFit.h.

92 {};

Member Function Documentation

◆ calculateArrivalTime()

float MahiFit::calculateArrivalTime ( unsigned int iBX ) const

private

Definition at line 319 of file MahiFit.cc.

                                                               {
   if (nnlsWork_.nPulseTot > 1) {
     SamplePulseMatrix pulseDerivMatTMP = nnlsWork_.pulseDerivMat;
     for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
       nnlsWork_.pulseDerivMat.col(iBX) = pulseDerivMatTMP.col(nnlsWork_.bxs.coeff(iBX) + nnlsWork_.bxOffset);
     }
   }
  
   for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
     nnlsWork_.pulseDerivMat.col(iBX) *= nnlsWork_.ampVec.coeff(iBX);
   }
  
   //FIXME: avoid solve full equation for one element
   SampleVector residuals = nnlsWork_.pulseMat * nnlsWork_.ampVec - nnlsWork_.amplitudes;
   PulseVector solution = nnlsWork_.pulseDerivMat.colPivHouseholderQr().solve(residuals);
   float t = std::clamp((float)solution.coeff(itIndex), -timeLimit_, timeLimit_);
  
   return t;
 }

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, nnlsWork_, MahiNnlsWorkspace::nPulseTot, MahiNnlsWorkspace::pulseDerivMat, MahiNnlsWorkspace::pulseMat, OrderedSet::t, and timeLimit_.

Referenced by doFit().

◆ calculateChiSq()

float MahiFit::calculateChiSq ( ) const

private

Definition at line 441 of file MahiFit.cc.

                                     {
   return (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat * nnlsWork_.ampVec - nnlsWork_.amplitudes))
       .squaredNorm();
 }

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::covDecomp, nnlsWork_, and MahiNnlsWorkspace::pulseMat.

Referenced by minimize().

◆ doFit()

void MahiFit::doFit	(	std::array< float, 3 > &	correctedOutput,
		const int	nbx
	)		const

Definition at line 115 of file MahiFit.cc.

                                                                       {
   unsigned int bxSize = 1;
  
   if (nbx == 1) {
     nnlsWork_.bxOffset = 0;
   } else {
     bxSize = bxSizeConf_;
     nnlsWork_.bxOffset = static_cast<int>(nnlsWork_.tsOffset) >= bxOffsetConf_ ? bxOffsetConf_ : nnlsWork_.tsOffset;
   }
  
   nnlsWork_.nPulseTot = bxSize;
  
   if (dynamicPed_)
     nnlsWork_.nPulseTot++;
   nnlsWork_.bxs.setZero(nnlsWork_.nPulseTot);
  
   if (nbx == 1) {
     nnlsWork_.bxs.coeffRef(0) = 0;
   } else {
     for (unsigned int iBX = 0; iBX < bxSize; ++iBX) {
       nnlsWork_.bxs.coeffRef(iBX) =
           activeBXs_[iBX] -
           ((static_cast<int>(nnlsWork_.tsOffset) + activeBXs_[0]) >= 0 ? 0 : (nnlsWork_.tsOffset + activeBXs_[0]));
     }
   }
  
   nnlsWork_.maxoffset = nnlsWork_.bxs.coeff(bxSize - 1);
   if (dynamicPed_)
     nnlsWork_.bxs[nnlsWork_.nPulseTot - 1] = pedestalBX_;
  
   nnlsWork_.pulseMat.setZero(nnlsWork_.tsSize, nnlsWork_.nPulseTot);
   nnlsWork_.pulseDerivMat.setZero(nnlsWork_.tsSize, nnlsWork_.nPulseTot);
  
   FullSampleVector pulseShapeArray;
   FullSampleVector pulseDerivArray;
   FullSampleMatrix pulseCov;
  
   int offset = 0;
   for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
     offset = nnlsWork_.bxs.coeff(iBX);
  
     if (offset == pedestalBX_) {
       nnlsWork_.pulseMat.col(iBX) = SampleVector::Ones(nnlsWork_.tsSize);
       nnlsWork_.pulseDerivMat.col(iBX) = SampleVector::Zero(nnlsWork_.tsSize);
     } else {
       pulseShapeArray.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
       if (calculateArrivalTime_)
         pulseDerivArray.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
       pulseCov.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset,
                        nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
       nnlsWork_.pulseCovArray[iBX].setZero(nnlsWork_.tsSize, nnlsWork_.tsSize);
  
       updatePulseShape(
           nnlsWork_.amplitudes.coeff(nnlsWork_.tsOffset + offset), pulseShapeArray, pulseDerivArray, pulseCov);
  
       nnlsWork_.pulseMat.col(iBX) = pulseShapeArray.segment(nnlsWork_.maxoffset - offset, nnlsWork_.tsSize);
       if (calculateArrivalTime_)
         nnlsWork_.pulseDerivMat.col(iBX) = pulseDerivArray.segment(nnlsWork_.maxoffset - offset, nnlsWork_.tsSize);
       nnlsWork_.pulseCovArray[iBX] = pulseCov.block(
           nnlsWork_.maxoffset - offset, nnlsWork_.maxoffset - offset, nnlsWork_.tsSize, nnlsWork_.tsSize);
     }
   }
  
   const float chiSq = minimize();
  
   bool foundintime = false;
   unsigned int ipulseintime = 0;
  
   for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
     if (nnlsWork_.bxs.coeff(iBX) == 0) {
       ipulseintime = iBX;
       foundintime = true;
       break;
     }
   }
  
   if (foundintime) {
     correctedOutput.at(0) = nnlsWork_.ampVec.coeff(ipulseintime);  //charge
     if (correctedOutput.at(0) != 0) {
       // fixME store the timeslew
       float arrivalTime = 0.f;
       if (calculateArrivalTime_)
         arrivalTime = calculateArrivalTime(ipulseintime);
       correctedOutput.at(1) = arrivalTime;  //time
     } else
       correctedOutput.at(1) = -9999.f;  //time
  
     correctedOutput.at(2) = chiSq;  //chi2
   }
 }

Referenced by phase1Apply().

◆ minimize()

const float MahiFit::minimize ( ) const

private

Definition at line 206 of file MahiFit.cc.

                                     {
   nnlsWork_.invcovp.setZero(nnlsWork_.tsSize, nnlsWork_.nPulseTot);
   nnlsWork_.ampVec.setZero(nnlsWork_.nPulseTot);
  
   SampleMatrix invCovMat;
   invCovMat.setConstant(nnlsWork_.tsSize, nnlsWork_.tsSize, nnlsWork_.pedVal);
   invCovMat += nnlsWork_.noiseTerms.asDiagonal();
  
   float oldChiSq = 9999;
   float chiSq = oldChiSq;
  
   for (int iter = 1; iter < nMaxItersMin_; ++iter) {
     updateCov(invCovMat);
  
     if (nnlsWork_.nPulseTot > 1) {
       nnls();
     } else {
       onePulseMinimize();
     }
  
     const float newChiSq = calculateChiSq();
     const float deltaChiSq = newChiSq - chiSq;
  
     if (newChiSq == oldChiSq && newChiSq < chiSq) {
       break;
     }
     oldChiSq = chiSq;
     chiSq = newChiSq;
  
     if (std::abs(deltaChiSq) < deltaChiSqThresh_)
       break;
   }
  
   return chiSq;
 }

References funct::abs(), MahiNnlsWorkspace::ampVec, calculateChiSq(), deltaChiSqThresh_, MahiNnlsWorkspace::invcovp, nMaxItersMin_, nnls(), nnlsWork_, MahiNnlsWorkspace::noiseTerms, MahiNnlsWorkspace::nPulseTot, onePulseMinimize(), MahiNnlsWorkspace::pedVal, MahiNnlsWorkspace::tsSize, and updateCov().

Referenced by doFit().

◆ nnls()

void MahiFit::nnls ( ) const

private

Definition at line 339 of file MahiFit.cc.

                          {
   const unsigned int npulse = nnlsWork_.nPulseTot;
   const unsigned int nsamples = nnlsWork_.tsSize;
  
   PulseVector updateWork;
  
   nnlsWork_.invcovp = nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat);
   nnlsWork_.aTaMat = nnlsWork_.invcovp.transpose() * nnlsWork_.invcovp;
   nnlsWork_.aTbVec = nnlsWork_.invcovp.transpose() * (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes));
  
   int iter = 0;
   Index idxwmax = 0;
   float wmax = 0.0f;
   float threshold = nnlsThresh_;
  
   while (true) {
     if (iter > 0 || nnlsWork_.nP == 0) {
       if (nnlsWork_.nP == std::min(npulse, nsamples))
         break;
  
       const unsigned int nActive = npulse - nnlsWork_.nP;
       // exit if there are no more pulses to constrain
       if (nActive == 0)
         break;
  
       updateWork = nnlsWork_.aTbVec - nnlsWork_.aTaMat * nnlsWork_.ampVec;
  
       Index idxwmaxprev = idxwmax;
       float wmaxprev = wmax;
       wmax = updateWork.tail(nActive).maxCoeff(&idxwmax);
  
       if (wmax < threshold || (idxwmax == idxwmaxprev && wmax == wmaxprev)) {
         break;
       }
  
       if (iter >= nMaxItersNNLS_) {
         break;
       }
  
       //unconstrain parameter
       idxwmax += nnlsWork_.nP;
       nnlsUnconstrainParameter(idxwmax);
     }
  
     while (true) {
       if (nnlsWork_.nP == 0)
         break;
  
       PulseVector ampvecpermtest = nnlsWork_.ampVec.head(nnlsWork_.nP);
  
       solveSubmatrix(nnlsWork_.aTaMat, nnlsWork_.aTbVec, ampvecpermtest, nnlsWork_.nP);
  
       //check solution
       if (ampvecpermtest.head(nnlsWork_.nP).minCoeff() > 0.f) {
         nnlsWork_.ampVec.head(nnlsWork_.nP) = ampvecpermtest.head(nnlsWork_.nP);
         break;
       }
  
       //update parameter vector
       Index minratioidx = 0;
  
       // no realizable optimization here (because it autovectorizes!)
       float minratio = std::numeric_limits<float>::max();
       for (unsigned int ipulse = 0; ipulse < nnlsWork_.nP; ++ipulse) {
         if (ampvecpermtest.coeff(ipulse) <= 0.f) {
           const float c_ampvec = nnlsWork_.ampVec.coeff(ipulse);
           const float ratio = c_ampvec / (c_ampvec - ampvecpermtest.coeff(ipulse));
           if (ratio < minratio) {
             minratio = ratio;
             minratioidx = ipulse;
           }
         }
       }
       nnlsWork_.ampVec.head(nnlsWork_.nP) +=
           minratio * (ampvecpermtest.head(nnlsWork_.nP) - nnlsWork_.ampVec.head(nnlsWork_.nP));
  
       //avoid numerical problems with later ==0. check
       nnlsWork_.ampVec.coeffRef(minratioidx) = 0.f;
  
       nnlsConstrainParameter(minratioidx);
     }
  
     ++iter;
  
     //adaptive convergence threshold to avoid infinite loops but still
     //ensure best value is used
     if (iter % 10 == 0) {
       threshold *= 10.;
     }
   }
 }

Referenced by minimize().

◆ nnlsConstrainParameter()

void MahiFit::nnlsConstrainParameter ( Index minratioidx ) const

private

Definition at line 486 of file MahiFit.cc.

                                                             {
   if (minratioidx != (nnlsWork_.nP - 1)) {
     nnlsWork_.aTaMat.col(nnlsWork_.nP - 1).swap(nnlsWork_.aTaMat.col(minratioidx));
     nnlsWork_.aTaMat.row(nnlsWork_.nP - 1).swap(nnlsWork_.aTaMat.row(minratioidx));
     nnlsWork_.pulseMat.col(nnlsWork_.nP - 1).swap(nnlsWork_.pulseMat.col(minratioidx));
     Eigen::numext::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP - 1), nnlsWork_.aTbVec.coeffRef(minratioidx));
     Eigen::numext::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP - 1), nnlsWork_.ampVec.coeffRef(minratioidx));
     Eigen::numext::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP - 1), nnlsWork_.bxs.coeffRef(minratioidx));
   }
   --nnlsWork_.nP;
 }

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::bxs, nnlsWork_, MahiNnlsWorkspace::nP, MahiNnlsWorkspace::pulseMat, and edm::swap().

Referenced by nnls().

◆ nnlsUnconstrainParameter()

void MahiFit::nnlsUnconstrainParameter ( Index idxp ) const

private

Definition at line 474 of file MahiFit.cc.

                                                        {
   if (idxp != nnlsWork_.nP) {
     nnlsWork_.aTaMat.col(nnlsWork_.nP).swap(nnlsWork_.aTaMat.col(idxp));
     nnlsWork_.aTaMat.row(nnlsWork_.nP).swap(nnlsWork_.aTaMat.row(idxp));
     nnlsWork_.pulseMat.col(nnlsWork_.nP).swap(nnlsWork_.pulseMat.col(idxp));
     Eigen::numext::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP), nnlsWork_.aTbVec.coeffRef(idxp));
     Eigen::numext::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP), nnlsWork_.ampVec.coeffRef(idxp));
     Eigen::numext::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP), nnlsWork_.bxs.coeffRef(idxp));
   }
   ++nnlsWork_.nP;
 }

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::bxs, nnlsWork_, MahiNnlsWorkspace::nP, MahiNnlsWorkspace::pulseMat, and edm::swap().

Referenced by nnls().

◆ onePulseMinimize()

void MahiFit::onePulseMinimize ( ) const

private

Definition at line 431 of file MahiFit.cc.

                                      {
   nnlsWork_.invcovp = nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat);
  
   float aTaCoeff = (nnlsWork_.invcovp.transpose() * nnlsWork_.invcovp).coeff(0);
   float aTbCoeff =
       (nnlsWork_.invcovp.transpose() * (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes))).coeff(0);
  
   nnlsWork_.ampVec.coeffRef(0) = std::max(0.f, aTbCoeff / aTaCoeff);
 }

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::covDecomp, f, MahiNnlsWorkspace::invcovp, SiStripPI::max, nnlsWork_, and MahiNnlsWorkspace::pulseMat.

Referenced by minimize().

◆ phase1Apply()

void MahiFit::phase1Apply	(	const HBHEChannelInfo &	channelData,
		float &	reconstructedEnergy,
		float &	reconstructedTime,
		bool &	useTriple,
		float &	chi2
	)		const

Definition at line 46 of file MahiFit.cc.

                                              {
   assert(channelData.nSamples() == 8 || channelData.nSamples() == 10);
  
   resetWorkspace();
  
   nnlsWork_.tsOffset = channelData.soi();
  
   std::array<float, 3> reconstructedVals{{0.0f, -9999.f, -9999.f}};
  
   double tsTOT = 0, tstrig = 0;  // in GeV
   for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
     auto const amplitude = channelData.tsRawCharge(iTS) - channelData.tsPedestal(iTS);
  
     nnlsWork_.amplitudes.coeffRef(iTS) = amplitude;
  
     //ADC granularity
     auto const noiseADC = norm_ * channelData.tsDFcPerADC(iTS);
  
     //Electronic pedestal
     auto const pedWidth = channelData.tsPedestalWidth(iTS);
  
     //Photostatistics
     auto const noisePhoto = (amplitude > pedWidth) ? std::sqrt(amplitude * channelData.fcByPE()) : 0.f;
  
     //Total uncertainty from all sources
     nnlsWork_.noiseTerms.coeffRef(iTS) = noiseADC * noiseADC + noisePhoto * noisePhoto + pedWidth * pedWidth;
  
     tsTOT += amplitude;
     if (iTS == nnlsWork_.tsOffset)
       tstrig += amplitude;
   }
  
   tsTOT *= channelData.tsGain(0);
   tstrig *= channelData.tsGain(0);
  
   useTriple = false;
   if (tstrig >= ts4Thresh_ && tsTOT > 0) {
     //Average pedestal width (for covariance matrix constraint)
     nnlsWork_.pedVal = 0.25f * (channelData.tsPedestalWidth(0) * channelData.tsPedestalWidth(0) +
                                 channelData.tsPedestalWidth(1) * channelData.tsPedestalWidth(1) +
                                 channelData.tsPedestalWidth(2) * channelData.tsPedestalWidth(2) +
                                 channelData.tsPedestalWidth(3) * channelData.tsPedestalWidth(3));
  
     // only do pre-fit with 1 pulse if chiSq threshold is positive
     if (chiSqSwitch_ > 0) {
       doFit(reconstructedVals, 1);
       if (reconstructedVals[2] > chiSqSwitch_) {
         doFit(reconstructedVals, 0);  //nbx=0 means use configured BXs
         useTriple = true;
       }
     } else {
       doFit(reconstructedVals, 0);
       useTriple = true;
     }
   } else {
     reconstructedVals.at(0) = 0.f;      //energy
     reconstructedVals.at(1) = -9999.f;  //time
     reconstructedVals.at(2) = -9999.f;  //chi2
   }
  
   reconstructedEnergy = reconstructedVals[0] * channelData.tsGain(0);
   reconstructedTime = reconstructedVals[1];
   chi2 = reconstructedVals[2];
 }

References CustomPhysics_cfi::amplitude, MahiNnlsWorkspace::amplitudes, cms::cuda::assert(), hltPixelTracks_cff::chi2, chiSqSwitch_, doFit(), f, HBHEChannelInfo::fcByPE(), nnlsWork_, MahiNnlsWorkspace::noiseTerms, norm_, HBHEChannelInfo::nSamples(), MahiNnlsWorkspace::pedVal, resetWorkspace(), HBHEChannelInfo::soi(), mathSSE::sqrt(), ts4Thresh_, HBHEChannelInfo::tsDFcPerADC(), HBHEChannelInfo::tsGain(), MahiNnlsWorkspace::tsOffset, HBHEChannelInfo::tsPedestal(), HBHEChannelInfo::tsPedestalWidth(), HBHEChannelInfo::tsRawCharge(), and MahiNnlsWorkspace::tsSize.

Referenced by phase1Debug(), and SimpleHBHEPhase1Algo::reconstruct().

◆ phase1Debug()

void MahiFit::phase1Debug	(	const HBHEChannelInfo &	channelData,
		MahiDebugInfo &	mdi
	)		const

Definition at line 498 of file MahiFit.cc.

                                                                                       {
   float recoEnergy, recoTime, chi2;
   bool use3;
   phase1Apply(channelData, recoEnergy, recoTime, use3, chi2);
  
   mdi.nSamples = channelData.nSamples();
   mdi.soi = channelData.soi();
  
   mdi.use3 = use3;
  
   mdi.inTimeConst = nnlsWork_.dt;
   mdi.inPedAvg = 0.25 * (channelData.tsPedestalWidth(0) * channelData.tsPedestalWidth(0) +
                          channelData.tsPedestalWidth(1) * channelData.tsPedestalWidth(1) +
                          channelData.tsPedestalWidth(2) * channelData.tsPedestalWidth(2) +
                          channelData.tsPedestalWidth(3) * channelData.tsPedestalWidth(3));
   mdi.inGain = channelData.tsGain(0);
  
   for (unsigned int iTS = 0; iTS < channelData.nSamples(); ++iTS) {
     double charge = channelData.tsRawCharge(iTS);
     double ped = channelData.tsPedestal(iTS);
  
     mdi.inNoiseADC[iTS] = norm_ * channelData.tsDFcPerADC(iTS);
  
     if ((charge - ped) > channelData.tsPedestalWidth(iTS)) {
       mdi.inNoisePhoto[iTS] = sqrt((charge - ped) * channelData.fcByPE());
     } else {
       mdi.inNoisePhoto[iTS] = 0;
     }
  
     mdi.inPedestal[iTS] = channelData.tsPedestalWidth(iTS);
     mdi.totalUCNoise[iTS] = nnlsWork_.noiseTerms.coeffRef(iTS);
  
     if (channelData.hasTimeInfo()) {
       mdi.inputTDC[iTS] = channelData.tsRiseTime(iTS);
     } else {
       mdi.inputTDC[iTS] = -1;
     }
   }
  
   mdi.arrivalTime = recoTime;
   mdi.chiSq = chi2;
  
   for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
     if (nnlsWork_.bxs.coeff(iBX) == 0) {
       mdi.mahiEnergy = nnlsWork_.ampVec.coeff(iBX);
       for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
         mdi.count[iTS] = iTS;
         mdi.inputTS[iTS] = nnlsWork_.amplitudes.coeff(iTS);
         mdi.itPulse[iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
       }
     } else if (nnlsWork_.bxs.coeff(iBX) == pedestalBX_) {
       mdi.pedEnergy = nnlsWork_.ampVec.coeff(iBX);
     } else if (nnlsWork_.bxs.coeff(iBX) >= -3 && nnlsWork_.bxs.coeff(iBX) <= 4) {
       int ootIndex = nnlsWork_.bxs.coeff(iBX);
       if (ootIndex > 0)
         ootIndex += 2;
       else
         ootIndex += 3;
       mdi.ootEnergy[ootIndex] = nnlsWork_.ampVec.coeff(iBX);
       for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
         mdi.ootPulse[ootIndex][iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
       }
     }
   }
 }

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiDebugInfo::arrivalTime, MahiNnlsWorkspace::bxs, ALCARECOTkAlJpsiMuMu_cff::charge, hltPixelTracks_cff::chi2, MahiDebugInfo::chiSq, MahiDebugInfo::count, MahiNnlsWorkspace::dt, HBHEChannelInfo::fcByPE(), HBHEChannelInfo::hasTimeInfo(), MahiDebugInfo::inGain, MahiDebugInfo::inNoiseADC, MahiDebugInfo::inNoisePhoto, MahiDebugInfo::inPedAvg, MahiDebugInfo::inPedestal, MahiDebugInfo::inputTDC, MahiDebugInfo::inputTS, MahiDebugInfo::inTimeConst, MahiDebugInfo::itPulse, MahiDebugInfo::mahiEnergy, nnlsWork_, MahiNnlsWorkspace::noiseTerms, norm_, MahiNnlsWorkspace::nPulseTot, MahiDebugInfo::nSamples, HBHEChannelInfo::nSamples(), MahiDebugInfo::ootEnergy, MahiDebugInfo::ootPulse, MahiDebugInfo::pedEnergy, pedestalBX_, phase1Apply(), MahiNnlsWorkspace::pulseMat, MahiDebugInfo::soi, HBHEChannelInfo::soi(), mathSSE::sqrt(), MahiDebugInfo::totalUCNoise, HBHEChannelInfo::tsDFcPerADC(), HBHEChannelInfo::tsGain(), HBHEChannelInfo::tsPedestal(), HBHEChannelInfo::tsPedestalWidth(), HBHEChannelInfo::tsRawCharge(), HBHEChannelInfo::tsRiseTime(), MahiNnlsWorkspace::tsSize, and MahiDebugInfo::use3.

◆ resetPulseShapeTemplate()

void MahiFit::resetPulseShapeTemplate	(	const HcalPulseShapes::Shape &	ps,
		bool	hasTimeInfo,
		unsigned int	nSamples
	)

Definition at line 459 of file MahiFit.cc.

                                                                                                              {
   ++cntsetPulseShape_;
  
   // only the pulse shape itself from PulseShapeFunctor is used for Mahi
   // the uncertainty terms calculated inside PulseShapeFunctor are used for Method 2 only
   psfPtr_.reset(new FitterFuncs::PulseShapeFunctor(ps, false, false, false, 1, 0, 0, HcalConst::maxSamples));
  
   // 1 sigma time constraint
   nnlsWork_.dt = hasTimeInfo ? timeSigmaSiPM_ : timeSigmaHPD_;
  
   nnlsWork_.tsSize = nSamples;
   nnlsWork_.amplitudes.resize(nnlsWork_.tsSize);
   nnlsWork_.noiseTerms.resize(nnlsWork_.tsSize);
 }

References MahiNnlsWorkspace::amplitudes, cntsetPulseShape_, MahiNnlsWorkspace::dt, HcalConst::maxSamples, nnlsWork_, MahiNnlsWorkspace::noiseTerms, PresampleTask_cfi::nSamples, psfPtr_, timeSigmaHPD_, timeSigmaSiPM_, and MahiNnlsWorkspace::tsSize.

Referenced by setPulseShapeTemplate().

◆ resetWorkspace()

void MahiFit::resetWorkspace ( ) const

private

Definition at line 613 of file MahiFit.cc.

                                    {
   nnlsWork_.nPulseTot = 0;
   nnlsWork_.tsOffset = 0;
   nnlsWork_.bxOffset = 0;
   nnlsWork_.maxoffset = 0;
   nnlsWork_.nP = 0;
  
   // NOT SURE THIS IS NEEDED
   nnlsWork_.amplitudes.setZero();
   nnlsWork_.noiseTerms.setZero();
 }

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::maxoffset, nnlsWork_, MahiNnlsWorkspace::noiseTerms, MahiNnlsWorkspace::nP, MahiNnlsWorkspace::nPulseTot, and MahiNnlsWorkspace::tsOffset.

Referenced by phase1Apply().

◆ setParameters()

void MahiFit::setParameters	(	bool	iDynamicPed,
		double	iTS4Thresh,
		double	chiSqSwitch,
		bool	iApplyTimeSlew,
		HcalTimeSlew::BiasSetting	slewFlavor,
		bool	iCalculateArrivalTime,
		double	iMeanTime,
		double	iTimeSigmaHPD,
		double	iTimeSigmaSiPM,
		const std::vector< int > &	iActiveBXs,
		int	iNMaxItersMin,
		int	iNMaxItersNNLS,
		double	iDeltaChiSqThresh,
		double	iNnlsThresh
	)

Definition at line 7 of file MahiFit.cc.

                                                 {
   dynamicPed_ = iDynamicPed;
  
   ts4Thresh_ = iTS4Thresh;
   chiSqSwitch_ = chiSqSwitch;
  
   applyTimeSlew_ = iApplyTimeSlew;
   slewFlavor_ = slewFlavor;
  
   calculateArrivalTime_ = iCalculateArrivalTime;
   meanTime_ = iMeanTime;
   timeSigmaHPD_ = iTimeSigmaHPD;
   timeSigmaSiPM_ = iTimeSigmaSiPM;
  
   activeBXs_ = iActiveBXs;
  
   nMaxItersMin_ = iNMaxItersMin;
   nMaxItersNNLS_ = iNMaxItersNNLS;
  
   deltaChiSqThresh_ = iDeltaChiSqThresh;
   nnlsThresh_ = iNnlsThresh;
  
   bxOffsetConf_ = -(*std::min_element(activeBXs_.begin(), activeBXs_.end()));
   bxSizeConf_ = activeBXs_.size();
 }

References activeBXs_, applyTimeSlew_, bxOffsetConf_, bxSizeConf_, calculateArrivalTime_, HLT_2018_cff::chiSqSwitch, chiSqSwitch_, deltaChiSqThresh_, dynamicPed_, meanTime_, nMaxItersMin_, nMaxItersNNLS_, nnlsThresh_, slewFlavor_, timeSigmaHPD_, timeSigmaSiPM_, and ts4Thresh_.

◆ setPulseShapeTemplate()

void MahiFit::setPulseShapeTemplate	(	const HcalPulseShapes::Shape &	ps,
		bool	hasTimeInfo,
		const HcalTimeSlew *	hcalTimeSlewDelay,
		unsigned int	nSamples
	)

Definition at line 446 of file MahiFit.cc.

                                                            {
   if (!(&ps == currentPulseShape_)) {
     hcalTimeSlewDelay_ = hcalTimeSlewDelay;
     tsDelay1GeV_ = hcalTimeSlewDelay->delay(1.0, slewFlavor_);
  
     resetPulseShapeTemplate(ps, hasTimeInfo, nSamples);
     currentPulseShape_ = &ps;
   }
 }

References currentPulseShape_, HcalTimeSlew::delay(), hcalTimeSlewDelay_, PresampleTask_cfi::nSamples, resetPulseShapeTemplate(), slewFlavor_, and tsDelay1GeV_.

◆ solveSubmatrix()

void MahiFit::solveSubmatrix	(	PulseMatrix &	mat,
		PulseVector &	invec,
		PulseVector &	outvec,
		unsigned	nP
	)		const

private

Definition at line 564 of file MahiFit.cc.

                                                                                                          {
   using namespace Eigen;
   switch (nP) {  // pulse matrix is always square.
     case 10: {
       Matrix<float, 10, 10> temp = mat;
       outvec.head<10>() = temp.ldlt().solve(invec.head<10>());
     } break;
     case 9: {
       Matrix<float, 9, 9> temp = mat.topLeftCorner<9, 9>();
       outvec.head<9>() = temp.ldlt().solve(invec.head<9>());
     } break;
     case 8: {
       Matrix<float, 8, 8> temp = mat.topLeftCorner<8, 8>();
       outvec.head<8>() = temp.ldlt().solve(invec.head<8>());
     } break;
     case 7: {
       Matrix<float, 7, 7> temp = mat.topLeftCorner<7, 7>();
       outvec.head<7>() = temp.ldlt().solve(invec.head<7>());
     } break;
     case 6: {
       Matrix<float, 6, 6> temp = mat.topLeftCorner<6, 6>();
       outvec.head<6>() = temp.ldlt().solve(invec.head<6>());
     } break;
     case 5: {
       Matrix<float, 5, 5> temp = mat.topLeftCorner<5, 5>();
       outvec.head<5>() = temp.ldlt().solve(invec.head<5>());
     } break;
     case 4: {
       Matrix<float, 4, 4> temp = mat.topLeftCorner<4, 4>();
       outvec.head<4>() = temp.ldlt().solve(invec.head<4>());
     } break;
     case 3: {
       Matrix<float, 3, 3> temp = mat.topLeftCorner<3, 3>();
       outvec.head<3>() = temp.ldlt().solve(invec.head<3>());
     } break;
     case 2: {
       Matrix<float, 2, 2> temp = mat.topLeftCorner<2, 2>();
       outvec.head<2>() = temp.ldlt().solve(invec.head<2>());
     } break;
     case 1: {
       Matrix<float, 1, 1> temp = mat.topLeftCorner<1, 1>();
       outvec.head<1>() = temp.ldlt().solve(invec.head<1>());
     } break;
     default:
       throw cms::Exception("HcalMahiWeirdState")
           << "Weird number of pulses encountered in Mahi, module is configured incorrectly!";
   }
 }

References Exception, and groupFilesInBlocks::temp.

Referenced by nnls().

◆ updateCov()

void MahiFit::updateCov ( const SampleMatrix & invCovMat ) const

private

Definition at line 298 of file MahiFit.cc.

                                                            {
   SampleMatrix invCovMat = samplecov;
  
   for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
     auto const amp = nnlsWork_.ampVec.coeff(iBX);
     if (amp == 0)
       continue;
  
     int offset = nnlsWork_.bxs.coeff(iBX);
  
     if (offset == pedestalBX_)
       continue;
     else {
       auto const ampsq = amp * amp;
       invCovMat += ampsq * nnlsWork_.pulseCovArray[offset + nnlsWork_.bxOffset];
     }
   }
  
   nnlsWork_.covDecomp.compute(invCovMat);
 }

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, MahiNnlsWorkspace::covDecomp, nnlsWork_, MahiNnlsWorkspace::nPulseTot, hltrates_dqm_sourceclient-live_cfg::offset, pedestalBX_, and MahiNnlsWorkspace::pulseCovArray.

Referenced by minimize().

◆ updatePulseShape()

void MahiFit::updatePulseShape	(	const float	itQ,
		FullSampleVector &	pulseShape,
		FullSampleVector &	pulseDeriv,
		FullSampleMatrix &	pulseCov
	)		const

private

Definition at line 242 of file MahiFit.cc.

                                                                  {
   float t0 = meanTime_;
  
   if (applyTimeSlew_) {
     if (itQ <= 1.f)
       t0 += tsDelay1GeV_;
     else
       t0 += hcalTimeSlewDelay_->delay(float(itQ), slewFlavor_);
   }
  
   std::array<double, HcalConst::maxSamples> pulseN;
   std::array<double, HcalConst::maxSamples> pulseM;
   std::array<double, HcalConst::maxSamples> pulseP;
  
   const float xx = t0;
   const float xxm = -nnlsWork_.dt + t0;
   const float xxp = nnlsWork_.dt + t0;
  
   psfPtr_->singlePulseShapeFuncMahi(&xx);
   psfPtr_->getPulseShape(pulseN);
  
   psfPtr_->singlePulseShapeFuncMahi(&xxm);
   psfPtr_->getPulseShape(pulseM);
  
   psfPtr_->singlePulseShapeFuncMahi(&xxp);
   psfPtr_->getPulseShape(pulseP);
  
   //in the 2018+ case where the sample of interest (SOI) is in TS3, add an extra offset to align
   //with previous SOI=TS4 case assumed by psfPtr_->getPulseShape()
   int delta = 4 - nnlsWork_.tsOffset;
  
   auto invDt = 0.5 / nnlsWork_.dt;
  
   for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
     pulseShape(iTS + nnlsWork_.maxoffset) = pulseN[iTS + delta];
     if (calculateArrivalTime_)
       pulseDeriv(iTS + nnlsWork_.maxoffset) = (pulseM[iTS + delta] - pulseP[iTS + delta]) * invDt;
  
     pulseM[iTS + delta] -= pulseN[iTS + delta];
     pulseP[iTS + delta] -= pulseN[iTS + delta];
   }
  
   for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
     for (unsigned int jTS = 0; jTS < iTS + 1; ++jTS) {
       auto const tmp = 0.5 * (pulseP[iTS + delta] * pulseP[jTS + delta] + pulseM[iTS + delta] * pulseM[jTS + delta]);
  
       pulseCov(jTS + nnlsWork_.maxoffset, iTS + nnlsWork_.maxoffset) = tmp;
       if (iTS != jTS)
         pulseCov(iTS + nnlsWork_.maxoffset, jTS + nnlsWork_.maxoffset) = tmp;
     }
   }
 }

References applyTimeSlew_, calculateArrivalTime_, HcalTimeSlew::delay(), dumpMFGeometry_cfg::delta, MahiNnlsWorkspace::dt, f, hcalTimeSlewDelay_, MahiNnlsWorkspace::maxoffset, meanTime_, nnlsWork_, psfPtr_, slewFlavor_, FrontierCondition_GT_autoExpress_cfi::t0, createJobs::tmp, tsDelay1GeV_, MahiNnlsWorkspace::tsOffset, MahiNnlsWorkspace::tsSize, and geometryCSVtoXML::xx.

Referenced by doFit().

Member Data Documentation

◆ activeBXs_

std::vector<int> MahiFit::activeBXs_

private

Definition at line 172 of file MahiFit.h.

Referenced by doFit(), and setParameters().

◆ applyTimeSlew_

bool MahiFit::applyTimeSlew_

private

Definition at line 162 of file MahiFit.h.

Referenced by setParameters(), and updatePulseShape().

◆ bxOffsetConf_

int MahiFit::bxOffsetConf_

private

Definition at line 181 of file MahiFit.h.

Referenced by doFit(), and setParameters().

◆ bxSizeConf_

unsigned int MahiFit::bxSizeConf_

private

Definition at line 180 of file MahiFit.h.

Referenced by doFit(), and setParameters().

◆ calculateArrivalTime_

bool MahiFit::calculateArrivalTime_

private

Definition at line 167 of file MahiFit.h.

Referenced by doFit(), setParameters(), and updatePulseShape().

◆ chiSqSwitch_

float MahiFit::chiSqSwitch_

private

Definition at line 160 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

◆ cntsetPulseShape_

int MahiFit::cntsetPulseShape_

private

Definition at line 184 of file MahiFit.h.

Referenced by resetPulseShapeTemplate().

◆ currentPulseShape_

const HcalPulseShapes::Shape* MahiFit::currentPulseShape_ = nullptr

Definition at line 126 of file MahiFit.h.

Referenced by setPulseShapeTemplate().

◆ deltaChiSqThresh_

float MahiFit::deltaChiSqThresh_

private

Definition at line 177 of file MahiFit.h.

Referenced by minimize(), and setParameters().

◆ dynamicPed_

bool MahiFit::dynamicPed_

private

Definition at line 158 of file MahiFit.h.

Referenced by doFit(), and setParameters().

◆ hcalTimeSlewDelay_

const HcalTimeSlew* MahiFit::hcalTimeSlewDelay_ = nullptr

Definition at line 127 of file MahiFit.h.

Referenced by setPulseShapeTemplate(), and updatePulseShape().

◆ meanTime_

float MahiFit::meanTime_

private

Definition at line 168 of file MahiFit.h.

Referenced by setParameters(), and updatePulseShape().

◆ nMaxItersMin_

int MahiFit::nMaxItersMin_

private

Definition at line 174 of file MahiFit.h.

Referenced by minimize(), and setParameters().

◆ nMaxItersNNLS_

int MahiFit::nMaxItersNNLS_

private

Definition at line 175 of file MahiFit.h.

Referenced by nnls(), and setParameters().

◆ nnlsThresh_

float MahiFit::nnlsThresh_

private

Definition at line 178 of file MahiFit.h.

Referenced by nnls(), and setParameters().

◆ nnlsWork_

MahiNnlsWorkspace MahiFit::nnlsWork_

mutableprivate

Definition at line 148 of file MahiFit.h.

Referenced by calculateArrivalTime(), calculateChiSq(), doFit(), minimize(), nnls(), nnlsConstrainParameter(), nnlsUnconstrainParameter(), onePulseMinimize(), phase1Apply(), phase1Debug(), resetPulseShapeTemplate(), resetWorkspace(), updateCov(), and updatePulseShape().

◆ norm_

float MahiFit::norm_ = (1.f / std::sqrt(12))

private

Definition at line 165 of file MahiFit.h.

Referenced by phase1Apply(), and phase1Debug().

◆ pedestalBX_

constexpr int MahiFit::pedestalBX_ = 100

staticconstexprprivate

Definition at line 151 of file MahiFit.h.

Referenced by doFit(), phase1Debug(), and updateCov().

◆ pfunctor_

std::unique_ptr<ROOT::Math::Functor> MahiFit::pfunctor_

private

Definition at line 186 of file MahiFit.h.

◆ psfPtr_

std::unique_ptr<FitterFuncs::PulseShapeFunctor> MahiFit::psfPtr_

private

Definition at line 185 of file MahiFit.h.

Referenced by resetPulseShapeTemplate(), and updatePulseShape().

◆ slewFlavor_

HcalTimeSlew::BiasSetting MahiFit::slewFlavor_

private

Definition at line 163 of file MahiFit.h.

Referenced by setParameters(), setPulseShapeTemplate(), and updatePulseShape().

◆ timeLimit_

constexpr float MahiFit::timeLimit_ = 12.5f

staticconstexprprivate

Definition at line 155 of file MahiFit.h.

Referenced by calculateArrivalTime().

◆ timeSigmaHPD_

float MahiFit::timeSigmaHPD_

private

Definition at line 169 of file MahiFit.h.

Referenced by resetPulseShapeTemplate(), and setParameters().

◆ timeSigmaSiPM_

float MahiFit::timeSigmaSiPM_

private

Definition at line 170 of file MahiFit.h.

Referenced by resetPulseShapeTemplate(), and setParameters().

◆ ts4Thresh_

float MahiFit::ts4Thresh_

private

Definition at line 159 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

◆ tsDelay1GeV_

float MahiFit::tsDelay1GeV_ = 0.f

private

Definition at line 164 of file MahiFit.h.

Referenced by setPulseShapeTemplate(), and updatePulseShape().

Public Types

Public Member Functions

Public Attributes

Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Member Typedef Documentation

◆ Index

Constructor & Destructor Documentation

◆ MahiFit()

◆ ~MahiFit()

Member Function Documentation

◆ calculateArrivalTime()

◆ calculateChiSq()

◆ doFit()

◆ minimize()

◆ nnls()

◆ nnlsConstrainParameter()

◆ nnlsUnconstrainParameter()

◆ onePulseMinimize()

◆ phase1Apply()

◆ phase1Debug()

◆ resetPulseShapeTemplate()

◆ resetWorkspace()

◆ setParameters()

◆ setPulseShapeTemplate()

◆ solveSubmatrix()

◆ updateCov()

◆ updatePulseShape()

Member Data Documentation

◆ activeBXs_

◆ applyTimeSlew_

◆ bxOffsetConf_

◆ bxSizeConf_

◆ calculateArrivalTime_

◆ chiSqSwitch_

◆ cntsetPulseShape_

◆ currentPulseShape_

◆ deltaChiSqThresh_

◆ dynamicPed_

◆ hcalTimeSlewDelay_

◆ meanTime_

◆ nMaxItersMin_

◆ nMaxItersNNLS_

◆ nnlsThresh_

◆ nnlsWork_

◆ norm_

◆ pedestalBX_

◆ pfunctor_

◆ psfPtr_

◆ slewFlavor_

◆ timeLimit_

◆ timeSigmaHPD_

◆ timeSigmaSiPM_

◆ ts4Thresh_

◆ tsDelay1GeV_