CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
List of all members | Public Types | Public Member Functions | Public Attributes | Private Types | Private Member Functions | Private Attributes | Static Private Attributes
MahiFit Class Reference

#include <MahiFit.h>

Public Types

typedef BXVector::Index Index
 

Public Member Functions

void doFit (std::array< float, 4 > &correctedOutput, const int nbx) const
 
 MahiFit ()
 
void phase1Apply (const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &soiPlusOneEnergy, float &reconstructedTime, bool &useTriple, float &chi2) const
 
void phase1Debug (const HBHEChannelInfo &channelData, MahiDebugInfo &mdi) const
 
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 (int pulseShapeId, const HcalPulseShapes &ps, bool hasTimeInfo, const HcalTimeSlew *hcalTimeSlewDelay, unsigned int nSamples)
 
 ~MahiFit ()
 

Public Attributes

const HcalTimeSlewhcalTimeSlewDelay_ = nullptr
 

Private Types

typedef std::pair< int,
std::shared_ptr
< FitterFuncs::PulseShapeFunctor > > 
ShapeWithId
 

Private Member Functions

float calculateArrivalTime (const 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 resetPulseShapeTemplate (int pulseShapeId, const HcalPulseShapes &ps, unsigned int nSamples)
 
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_ = 0
 
int currentPulseShapeId_ = INT_MIN
 
float deltaChiSqThresh_
 
bool dynamicPed_
 
std::vector< ShapeWithIdknownPulseShapes_
 
float meanTime_
 
int nMaxItersMin_
 
int nMaxItersNNLS_
 
float nnlsThresh_
 
MahiNnlsWorkspace nnlsWork_
 
float norm_ = (1.f / std::sqrt(12))
 
FitterFuncs::PulseShapeFunctorpsfPtr_ = nullptr
 
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 96 of file MahiFit.h.

Member Typedef Documentation

typedef BXVector::Index MahiFit::Index

Definition at line 133 of file MahiFit.h.

typedef std::pair<int, std::shared_ptr<FitterFuncs::PulseShapeFunctor> > MahiFit::ShapeWithId
private

Definition at line 137 of file MahiFit.h.

Constructor & Destructor Documentation

MahiFit::MahiFit ( )

Definition at line 5 of file MahiFit.cc.

5 {}
MahiFit::~MahiFit ( )
inline

Definition at line 99 of file MahiFit.h.

99 {};

Member Function Documentation

float MahiFit::calculateArrivalTime ( const unsigned int  iBX) const
private

Definition at line 348 of file MahiFit.cc.

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

Referenced by doFit().

348  {
349  if (nnlsWork_.nPulseTot > 1) {
350  SamplePulseMatrix pulseDerivMatTMP = nnlsWork_.pulseDerivMat;
351  for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
352  nnlsWork_.pulseDerivMat.col(iBX) = pulseDerivMatTMP.col(nnlsWork_.bxs.coeff(iBX) + nnlsWork_.bxOffset);
353  }
354  }
355 
356  for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
357  nnlsWork_.pulseDerivMat.col(iBX) *= nnlsWork_.ampVec.coeff(iBX);
358  }
359 
360  //FIXME: avoid solve full equation for one element
362  PulseVector solution = nnlsWork_.pulseDerivMat.colPivHouseholderQr().solve(residuals);
363  float t = std::clamp((float)solution.coeff(itIndex), -timeLimit_, timeLimit_);
364 
365  return t;
366 }
unsigned int nPulseTot
Definition: MahiFit.h:20
Eigen::Matrix< double, SampleVectorSize, Eigen::Dynamic, 0, SampleVectorSize, PulseVectorSize > SamplePulseMatrix
SampleVector amplitudes
Definition: MahiFit.h:31
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
PulseVector ampVec
Definition: MahiFit.h:56
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, PulseVectorSize, 1 > PulseVector
SamplePulseMatrix pulseDerivMat
Definition: MahiFit.h:52
static constexpr float timeLimit_
Definition: MahiFit.h:165
Eigen::Matrix< double, SampleVectorSize, 1 > SampleVector
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49
BXVector bxs
Definition: MahiFit.h:28
float MahiFit::calculateChiSq ( ) const
private

Definition at line 472 of file MahiFit.cc.

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

Referenced by minimize().

472  {
474  .squaredNorm();
475 }
SampleVector amplitudes
Definition: MahiFit.h:31
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
PulseVector ampVec
Definition: MahiFit.h:56
SampleDecompLLT covDecomp
Definition: MahiFit.h:62
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49
void MahiFit::doFit ( std::array< float, 4 > &  correctedOutput,
const int  nbx 
) const

Definition at line 131 of file MahiFit.cc.

References activeBXs_, MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::bxOffset, bxOffsetConf_, MahiNnlsWorkspace::bxs, bxSizeConf_, calculateArrivalTime(), calculateArrivalTime_, dynamicPed_, MahiNnlsWorkspace::maxoffset, minimize(), nnlsWork_, MahiNnlsWorkspace::nPulseTot, hltrates_dqm_sourceclient-live_cfg::offset, pedestalBX_, MahiNnlsWorkspace::pulseCovArray, MahiNnlsWorkspace::pulseDerivMat, MahiNnlsWorkspace::pulseMat, MahiNnlsWorkspace::tsOffset, MahiNnlsWorkspace::tsSize, and updatePulseShape().

Referenced by phase1Apply().

131  {
132  unsigned int bxSize = 1;
133 
134  if (nbx == 1) {
135  nnlsWork_.bxOffset = 0;
136  } else {
137  bxSize = bxSizeConf_;
139  }
140 
141  nnlsWork_.nPulseTot = bxSize;
142 
143  if (dynamicPed_)
145  nnlsWork_.bxs.setZero(nnlsWork_.nPulseTot);
146 
147  if (nbx == 1) {
148  nnlsWork_.bxs.coeffRef(0) = 0;
149  } else {
150  for (unsigned int iBX = 0; iBX < bxSize; ++iBX) {
151  nnlsWork_.bxs.coeffRef(iBX) =
152  activeBXs_[iBX] -
153  ((static_cast<int>(nnlsWork_.tsOffset) + activeBXs_[0]) >= 0 ? 0 : (nnlsWork_.tsOffset + activeBXs_[0]));
154  }
155  }
156 
157  nnlsWork_.maxoffset = nnlsWork_.bxs.coeff(bxSize - 1);
158  if (dynamicPed_)
160 
163 
164  FullSampleVector pulseShapeArray;
165  FullSampleVector pulseDerivArray;
166  FullSampleMatrix pulseCov;
167 
168  int offset = 0;
169  for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
170  offset = nnlsWork_.bxs.coeff(iBX);
171 
172  if (offset == pedestalBX_) {
173  nnlsWork_.pulseMat.col(iBX) = SampleVector::Ones(nnlsWork_.tsSize);
174  nnlsWork_.pulseDerivMat.col(iBX) = SampleVector::Zero(nnlsWork_.tsSize);
175  } else {
176  pulseShapeArray.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
178  pulseDerivArray.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset);
179  pulseCov.setZero(nnlsWork_.tsSize + nnlsWork_.maxoffset + nnlsWork_.bxOffset,
182 
184  nnlsWork_.amplitudes.coeff(nnlsWork_.tsOffset + offset), pulseShapeArray, pulseDerivArray, pulseCov);
185 
186  nnlsWork_.pulseMat.col(iBX) = pulseShapeArray.segment(nnlsWork_.maxoffset - offset, nnlsWork_.tsSize);
188  nnlsWork_.pulseDerivMat.col(iBX) = pulseDerivArray.segment(nnlsWork_.maxoffset - offset, nnlsWork_.tsSize);
189  nnlsWork_.pulseCovArray[iBX] = pulseCov.block(
191  }
192  }
193 
194  const float chiSq = minimize();
195 
196  bool foundintime = false;
197  unsigned int ipulseintime = 0;
198 
199  for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
200  if (nnlsWork_.bxs.coeff(iBX) == 0) {
201  ipulseintime = iBX;
202  foundintime = true;
203  break;
204  }
205  }
206 
207  if (foundintime) {
208  correctedOutput.at(0) = nnlsWork_.ampVec.coeff(ipulseintime); //charge
209  correctedOutput.at(3) = nnlsWork_.ampVec.coeff(ipulseintime + 1U); //charge for SOI+1
210  if (correctedOutput.at(0) != 0) {
211  // fixME store the timeslew
212  float arrivalTime = 0.f;
214  arrivalTime = calculateArrivalTime(ipulseintime);
215  correctedOutput.at(1) = arrivalTime; //time
216  } else
217  correctedOutput.at(1) = -9999.f; //time
218 
219  correctedOutput.at(2) = chiSq; //chi2
220  }
221 }
unsigned int nPulseTot
Definition: MahiFit.h:20
SampleVector amplitudes
Definition: MahiFit.h:31
void updatePulseShape(const float itQ, FullSampleVector &pulseShape, FullSampleVector &pulseDeriv, FullSampleMatrix &pulseCov) const
Definition: MahiFit.cc:268
Eigen::Matrix< double, FullSampleVectorSize, FullSampleVectorSize > FullSampleMatrix
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
Eigen::Matrix< double, FullSampleVectorSize, 1 > FullSampleVector
static constexpr int pedestalBX_
Definition: MahiFit.h:161
PulseVector ampVec
Definition: MahiFit.h:56
std::array< SampleMatrix, MaxPVSize > pulseCovArray
Definition: MahiFit.h:46
SamplePulseMatrix pulseDerivMat
Definition: MahiFit.h:52
bool calculateArrivalTime_
Definition: MahiFit.h:177
int bxOffsetConf_
Definition: MahiFit.h:191
unsigned int bxSizeConf_
Definition: MahiFit.h:190
const float minimize() const
Definition: MahiFit.cc:223
float calculateArrivalTime(const unsigned int iBX) const
Definition: MahiFit.cc:348
std::vector< int > activeBXs_
Definition: MahiFit.h:182
unsigned int tsOffset
Definition: MahiFit.h:22
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49
bool dynamicPed_
Definition: MahiFit.h:168
unsigned int tsSize
Definition: MahiFit.h:21
BXVector bxs
Definition: MahiFit.h:28
const float MahiFit::minimize ( ) const
private

Definition at line 223 of file MahiFit.cc.

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

Referenced by doFit().

223  {
226 
227  SampleMatrix invCovMat;
228  invCovMat.setConstant(nnlsWork_.tsSize, nnlsWork_.tsSize, nnlsWork_.pedVal);
229  invCovMat += nnlsWork_.noiseTerms.asDiagonal();
230 
231  //Add off-Diagonal components up to first order
232  if (nnlsWork_.noisecorr != 0) {
233  for (unsigned int i = 1; i < nnlsWork_.tsSize; ++i) {
234  auto const noiseCorrTerm = nnlsWork_.noisecorr * nnlsWork_.pedVals.coeff(i - 1) * nnlsWork_.pedVals.coeff(i);
235  invCovMat(i - 1, i) += noiseCorrTerm;
236  invCovMat(i, i - 1) += noiseCorrTerm;
237  }
238  }
239 
240  float oldChiSq = 9999;
241  float chiSq = oldChiSq;
242 
243  for (int iter = 1; iter < nMaxItersMin_; ++iter) {
244  updateCov(invCovMat);
245 
246  if (nnlsWork_.nPulseTot > 1) {
247  nnls();
248  } else {
250  }
251 
252  const float newChiSq = calculateChiSq();
253  const float deltaChiSq = newChiSq - chiSq;
254 
255  if (newChiSq == oldChiSq && newChiSq < chiSq) {
256  break;
257  }
258  oldChiSq = chiSq;
259  chiSq = newChiSq;
260 
261  if (std::abs(deltaChiSq) < deltaChiSqThresh_)
262  break;
263  }
264 
265  return chiSq;
266 }
unsigned int nPulseTot
Definition: MahiFit.h:20
float noisecorr
Definition: MahiFit.h:42
void nnls() const
Definition: MahiFit.cc:368
SamplePulseMatrix invcovp
Definition: MahiFit.h:58
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
PulseVector ampVec
Definition: MahiFit.h:56
void onePulseMinimize() const
Definition: MahiFit.cc:462
int nMaxItersMin_
Definition: MahiFit.h:184
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void updateCov(const SampleMatrix &invCovMat) const
Definition: MahiFit.cc:328
Eigen::Matrix< double, SampleVectorSize, SampleVectorSize > SampleMatrix
SampleVector noiseTerms
Definition: MahiFit.h:34
float calculateChiSq() const
Definition: MahiFit.cc:472
unsigned int tsSize
Definition: MahiFit.h:21
float deltaChiSqThresh_
Definition: MahiFit.h:187
SampleVector pedVals
Definition: MahiFit.h:37
void MahiFit::nnls ( ) const
private

Definition at line 368 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::covDecomp, MahiNnlsWorkspace::invcovp, SiStripPI::max, SiStripPI::min, nMaxItersNNLS_, nnlsConstrainParameter(), nnlsThresh_, nnlsUnconstrainParameter(), nnlsWork_, MahiNnlsWorkspace::nP, MahiNnlsWorkspace::nPulseTot, MahiNnlsWorkspace::pulseMat, solveSubmatrix(), dtDQMClient_cfg::threshold, MahiNnlsWorkspace::tsSize, and cms::cuda::wmax.

Referenced by minimize().

368  {
369  const unsigned int npulse = nnlsWork_.nPulseTot;
370  const unsigned int nsamples = nnlsWork_.tsSize;
371 
372  PulseVector updateWork;
373  PulseVector ampvecpermtest;
374 
376  nnlsWork_.aTaMat.noalias() = nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.invcovp);
377  nnlsWork_.aTbVec.noalias() =
378  nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes));
379 
380  int iter = 0;
381  Index idxwmax = 0;
382  float wmax = 0.0f;
383  float threshold = nnlsThresh_;
384 
385  while (true) {
386  if (iter > 0 || nnlsWork_.nP == 0) {
387  if (nnlsWork_.nP == std::min(npulse, nsamples))
388  break;
389 
390  const unsigned int nActive = npulse - nnlsWork_.nP;
391  // exit if there are no more pulses to constrain
392  if (nActive == 0)
393  break;
394 
395  updateWork.noalias() = nnlsWork_.aTbVec - nnlsWork_.aTaMat.lazyProduct(nnlsWork_.ampVec);
396 
397  Index idxwmaxprev = idxwmax;
398  float wmaxprev = wmax;
399  wmax = updateWork.tail(nActive).maxCoeff(&idxwmax);
400 
401  if (wmax < threshold || (idxwmax == idxwmaxprev && wmax == wmaxprev)) {
402  break;
403  }
404 
405  if (iter >= nMaxItersNNLS_) {
406  break;
407  }
408 
409  //unconstrain parameter
410  idxwmax += nnlsWork_.nP;
411  nnlsUnconstrainParameter(idxwmax);
412  }
413 
414  while (true) {
415  if (nnlsWork_.nP == 0)
416  break;
417 
418  ampvecpermtest.noalias() = nnlsWork_.ampVec.head(nnlsWork_.nP);
419 
421 
422  //check solution
423  if (ampvecpermtest.head(nnlsWork_.nP).minCoeff() > 0.f) {
424  nnlsWork_.ampVec.head(nnlsWork_.nP) = ampvecpermtest.head(nnlsWork_.nP);
425  break;
426  }
427 
428  //update parameter vector
429  Index minratioidx = 0;
430 
431  // no realizable optimization here (because it autovectorizes!)
432  float minratio = std::numeric_limits<float>::max();
433  for (unsigned int ipulse = 0; ipulse < nnlsWork_.nP; ++ipulse) {
434  if (ampvecpermtest.coeff(ipulse) <= 0.f) {
435  const float c_ampvec = nnlsWork_.ampVec.coeff(ipulse);
436  const float ratio = c_ampvec / (c_ampvec - ampvecpermtest.coeff(ipulse));
437  if (ratio < minratio) {
438  minratio = ratio;
439  minratioidx = ipulse;
440  }
441  }
442  }
443  nnlsWork_.ampVec.head(nnlsWork_.nP) +=
444  minratio * (ampvecpermtest.head(nnlsWork_.nP) - nnlsWork_.ampVec.head(nnlsWork_.nP));
445 
446  //avoid numerical problems with later ==0. check
447  nnlsWork_.ampVec.coeffRef(minratioidx) = 0.f;
448 
449  nnlsConstrainParameter(minratioidx);
450  }
451 
452  ++iter;
453 
454  //adaptive convergence threshold to avoid infinite loops but still
455  //ensure best value is used
456  if (iter % 10 == 0) {
457  threshold *= 10.;
458  }
459  }
460 }
unsigned int nPulseTot
Definition: MahiFit.h:20
SamplePulseMatrix invcovp
Definition: MahiFit.h:58
void nnlsConstrainParameter(Index minratioidx) const
Definition: MahiFit.cc:540
SampleVector amplitudes
Definition: MahiFit.h:31
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
BXVector::Index Index
Definition: MahiFit.h:133
void solveSubmatrix(PulseMatrix &mat, PulseVector &invec, PulseVector &outvec, unsigned nP) const
Definition: MahiFit.cc:618
PulseVector ampVec
Definition: MahiFit.h:56
SampleDecompLLT covDecomp
Definition: MahiFit.h:62
void nnlsUnconstrainParameter(Index idxp) const
Definition: MahiFit.cc:528
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, PulseVectorSize, 1 > PulseVector
PulseVector aTbVec
Definition: MahiFit.h:60
float nnlsThresh_
Definition: MahiFit.h:188
PulseMatrix aTaMat
Definition: MahiFit.h:59
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49
unsigned int tsSize
Definition: MahiFit.h:21
unsigned int nP
Definition: MahiFit.h:55
int nMaxItersNNLS_
Definition: MahiFit.h:185
__host__ __device__ V V wmax
void MahiFit::nnlsConstrainParameter ( Index  minratioidx) const
private

Definition at line 540 of file MahiFit.cc.

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

Referenced by nnls().

540  {
541  if (minratioidx != (nnlsWork_.nP - 1)) {
542  nnlsWork_.aTaMat.col(nnlsWork_.nP - 1).swap(nnlsWork_.aTaMat.col(minratioidx));
543  nnlsWork_.aTaMat.row(nnlsWork_.nP - 1).swap(nnlsWork_.aTaMat.row(minratioidx));
544  nnlsWork_.pulseMat.col(nnlsWork_.nP - 1).swap(nnlsWork_.pulseMat.col(minratioidx));
545  Eigen::numext::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP - 1), nnlsWork_.aTbVec.coeffRef(minratioidx));
546  Eigen::numext::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP - 1), nnlsWork_.ampVec.coeffRef(minratioidx));
547  Eigen::numext::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP - 1), nnlsWork_.bxs.coeffRef(minratioidx));
548  }
549  --nnlsWork_.nP;
550 }
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
PulseVector ampVec
Definition: MahiFit.h:56
void swap(Association< C > &lhs, Association< C > &rhs)
Definition: Association.h:117
PulseVector aTbVec
Definition: MahiFit.h:60
PulseMatrix aTaMat
Definition: MahiFit.h:59
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49
unsigned int nP
Definition: MahiFit.h:55
BXVector bxs
Definition: MahiFit.h:28
void MahiFit::nnlsUnconstrainParameter ( Index  idxp) const
private

Definition at line 528 of file MahiFit.cc.

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

Referenced by nnls().

528  {
529  if (idxp != nnlsWork_.nP) {
530  nnlsWork_.aTaMat.col(nnlsWork_.nP).swap(nnlsWork_.aTaMat.col(idxp));
531  nnlsWork_.aTaMat.row(nnlsWork_.nP).swap(nnlsWork_.aTaMat.row(idxp));
532  nnlsWork_.pulseMat.col(nnlsWork_.nP).swap(nnlsWork_.pulseMat.col(idxp));
533  Eigen::numext::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP), nnlsWork_.aTbVec.coeffRef(idxp));
534  Eigen::numext::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP), nnlsWork_.ampVec.coeffRef(idxp));
535  Eigen::numext::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP), nnlsWork_.bxs.coeffRef(idxp));
536  }
537  ++nnlsWork_.nP;
538 }
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
PulseVector ampVec
Definition: MahiFit.h:56
void swap(Association< C > &lhs, Association< C > &rhs)
Definition: Association.h:117
PulseVector aTbVec
Definition: MahiFit.h:60
PulseMatrix aTaMat
Definition: MahiFit.h:59
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49
unsigned int nP
Definition: MahiFit.h:55
BXVector bxs
Definition: MahiFit.h:28
void MahiFit::onePulseMinimize ( ) const
private

Definition at line 462 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::covDecomp, validate-o2o-wbm::f, MahiNnlsWorkspace::invcovp, SiStripPI::max, nnlsWork_, and MahiNnlsWorkspace::pulseMat.

Referenced by minimize().

462  {
464 
465  float aTaCoeff = (nnlsWork_.invcovp.transpose() * nnlsWork_.invcovp).coeff(0);
466  float aTbCoeff =
467  (nnlsWork_.invcovp.transpose() * (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes))).coeff(0);
468 
469  nnlsWork_.ampVec.coeffRef(0) = std::max(0.f, aTbCoeff / aTaCoeff);
470 }
SamplePulseMatrix invcovp
Definition: MahiFit.h:58
SampleVector amplitudes
Definition: MahiFit.h:31
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
PulseVector ampVec
Definition: MahiFit.h:56
SampleDecompLLT covDecomp
Definition: MahiFit.h:62
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49
void MahiFit::phase1Apply ( const HBHEChannelInfo channelData,
float &  reconstructedEnergy,
float &  soiPlusOneEnergy,
float &  reconstructedTime,
bool &  useTriple,
float &  chi2 
) const

Definition at line 46 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, cms::cuda::assert(), chiSqSwitch_, doFit(), validate-o2o-wbm::f, HBHEChannelInfo::fcByPE(), nnlsWork_, MahiNnlsWorkspace::noisecorr, HBHEChannelInfo::noisecorr(), MahiNnlsWorkspace::noiseTerms, norm_, hgc_digi::nSamples, HBHEChannelInfo::nSamples(), MahiNnlsWorkspace::pedVal, MahiNnlsWorkspace::pedVals, resetWorkspace(), HBHEChannelInfo::soi(), ts4Thresh_, HBHEChannelInfo::tsDFcPerADC(), HBHEChannelInfo::tsGain(), MahiNnlsWorkspace::tsOffset, HBHEChannelInfo::tsPedestal(), HBHEChannelInfo::tsPedestalWidth(), HBHEChannelInfo::tsRawCharge(), and MahiNnlsWorkspace::tsSize.

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

51  {
52  const unsigned nSamples = channelData.nSamples();
53  const unsigned soi = channelData.soi();
54 
55  // Check some of the assumptions used by the subsequent code
56  assert(nSamples == 8 || nSamples == 10);
57  assert(nnlsWork_.tsSize <= nSamples);
58  assert(soi + 1U < nnlsWork_.tsSize);
59 
61 
62  nnlsWork_.tsOffset = soi;
63 
64  // Packing energy, time, chiSq, soiPlus1Energy, in that order
65  std::array<float, 4> reconstructedVals{{0.f, -9999.f, -9999.f, 0.f}};
66 
67  double tsTOT = 0, tstrig = 0; // in GeV
68  for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
69  auto const amplitude = channelData.tsRawCharge(iTS) - channelData.tsPedestal(iTS);
70 
71  nnlsWork_.amplitudes.coeffRef(iTS) = amplitude;
72 
73  //ADC granularity
74  auto const noiseADC = norm_ * channelData.tsDFcPerADC(iTS);
75 
76  //Electronic pedestal
77  auto const pedWidth = channelData.tsPedestalWidth(iTS);
78 
79  //Photostatistics
80  auto const noisePhotoSq = (amplitude > pedWidth) ? (amplitude * channelData.fcByPE()) : 0.f;
81 
82  //Total uncertainty from all sources
83  nnlsWork_.noiseTerms.coeffRef(iTS) = noiseADC * noiseADC + noisePhotoSq + pedWidth * pedWidth;
84 
85  nnlsWork_.pedVals.coeffRef(iTS) = pedWidth;
86 
87  tsTOT += amplitude;
88  if (iTS == soi)
89  tstrig += amplitude;
90  }
91 
92  // This preserves the original Mahi approach but will have to
93  // change in case we eventually calibrate gains per capId
94  const float gain0 = channelData.tsGain(0);
95  tsTOT *= gain0;
96  tstrig *= gain0;
97 
98  useTriple = false;
99  if (tstrig > ts4Thresh_ && tsTOT > 0) {
100  nnlsWork_.noisecorr = channelData.noisecorr();
101 
102  if (nnlsWork_.noisecorr != 0) {
103  nnlsWork_.pedVal = 0.f;
104  } else {
105  //Average pedestal width (for covariance matrix constraint)
106  nnlsWork_.pedVal = 0.25f * (channelData.tsPedestalWidth(0) * channelData.tsPedestalWidth(0) +
107  channelData.tsPedestalWidth(1) * channelData.tsPedestalWidth(1) +
108  channelData.tsPedestalWidth(2) * channelData.tsPedestalWidth(2) +
109  channelData.tsPedestalWidth(3) * channelData.tsPedestalWidth(3));
110  }
111 
112  // only do pre-fit with 1 pulse if chiSq threshold is positive
113  if (chiSqSwitch_ > 0) {
114  doFit(reconstructedVals, 1);
115  if (reconstructedVals[2] > chiSqSwitch_) {
116  doFit(reconstructedVals, 0); //nbx=0 means use configured BXs
117  useTriple = true;
118  }
119  } else {
120  doFit(reconstructedVals, 0);
121  useTriple = true;
122  }
123  }
124 
125  reconstructedEnergy = reconstructedVals[0] * gain0;
126  soiPlusOneEnergy = reconstructedVals[3] * gain0;
127  reconstructedTime = reconstructedVals[1];
128  chi2 = reconstructedVals[2];
129 }
float noisecorr
Definition: MahiFit.h:42
constexpr unsigned nSamples() const
float chiSqSwitch_
Definition: MahiFit.h:170
SampleVector amplitudes
Definition: MahiFit.h:31
constexpr double tsGain(const unsigned ts) const
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
float norm_
Definition: MahiFit.h:175
assert(be >=bs)
constexpr double tsPedestalWidth(const unsigned ts) const
void doFit(std::array< float, 4 > &correctedOutput, const int nbx) const
Definition: MahiFit.cc:131
constexpr double noisecorr() const
constexpr size_t nSamples
constexpr unsigned soi() const
constexpr double tsPedestal(const unsigned ts) const
constexpr float tsDFcPerADC(const unsigned ts) const
unsigned int tsOffset
Definition: MahiFit.h:22
void resetWorkspace() const
Definition: MahiFit.cc:667
constexpr double fcByPE() const
SampleVector noiseTerms
Definition: MahiFit.h:34
float ts4Thresh_
Definition: MahiFit.h:169
unsigned int tsSize
Definition: MahiFit.h:21
SampleVector pedVals
Definition: MahiFit.h:37
constexpr double tsRawCharge(const unsigned ts) const
void MahiFit::phase1Debug ( const HBHEChannelInfo channelData,
MahiDebugInfo mdi 
) const

Definition at line 552 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiDebugInfo::arrivalTime, MahiNnlsWorkspace::bxs, RecoTauCleanerPlugins::charge, HLT_FULL_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.

552  {
553  float recoEnergy, soiPlus1Energy, recoTime, chi2;
554  bool use3;
555  phase1Apply(channelData, recoEnergy, soiPlus1Energy, recoTime, use3, chi2);
556 
557  mdi.nSamples = channelData.nSamples();
558  mdi.soi = channelData.soi();
559 
560  mdi.use3 = use3;
561 
562  mdi.inTimeConst = nnlsWork_.dt;
563  mdi.inPedAvg = 0.25 * (channelData.tsPedestalWidth(0) * channelData.tsPedestalWidth(0) +
564  channelData.tsPedestalWidth(1) * channelData.tsPedestalWidth(1) +
565  channelData.tsPedestalWidth(2) * channelData.tsPedestalWidth(2) +
566  channelData.tsPedestalWidth(3) * channelData.tsPedestalWidth(3));
567  mdi.inGain = channelData.tsGain(0);
568 
569  for (unsigned int iTS = 0; iTS < channelData.nSamples(); ++iTS) {
570  double charge = channelData.tsRawCharge(iTS);
571  double ped = channelData.tsPedestal(iTS);
572 
573  mdi.inNoiseADC[iTS] = norm_ * channelData.tsDFcPerADC(iTS);
574 
575  if ((charge - ped) > channelData.tsPedestalWidth(iTS)) {
576  mdi.inNoisePhoto[iTS] = sqrt((charge - ped) * channelData.fcByPE());
577  } else {
578  mdi.inNoisePhoto[iTS] = 0;
579  }
580 
581  mdi.inPedestal[iTS] = channelData.tsPedestalWidth(iTS);
582  mdi.totalUCNoise[iTS] = nnlsWork_.noiseTerms.coeffRef(iTS);
583 
584  if (channelData.hasTimeInfo()) {
585  mdi.inputTDC[iTS] = channelData.tsRiseTime(iTS);
586  } else {
587  mdi.inputTDC[iTS] = -1;
588  }
589  }
590 
591  mdi.arrivalTime = recoTime;
592  mdi.chiSq = chi2;
593 
594  for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
595  if (nnlsWork_.bxs.coeff(iBX) == 0) {
596  mdi.mahiEnergy = nnlsWork_.ampVec.coeff(iBX);
597  for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
598  mdi.count[iTS] = iTS;
599  mdi.inputTS[iTS] = nnlsWork_.amplitudes.coeff(iTS);
600  mdi.itPulse[iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
601  }
602  } else if (nnlsWork_.bxs.coeff(iBX) == pedestalBX_) {
603  mdi.pedEnergy = nnlsWork_.ampVec.coeff(iBX);
604  } else if (nnlsWork_.bxs.coeff(iBX) >= -3 && nnlsWork_.bxs.coeff(iBX) <= 4) {
605  int ootIndex = nnlsWork_.bxs.coeff(iBX);
606  if (ootIndex > 0)
607  ootIndex += 2;
608  else
609  ootIndex += 3;
610  mdi.ootEnergy[ootIndex] = nnlsWork_.ampVec.coeff(iBX);
611  for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
612  mdi.ootPulse[ootIndex][iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
613  }
614  }
615  }
616 }
unsigned int nPulseTot
Definition: MahiFit.h:20
constexpr unsigned nSamples() const
float ootEnergy[7]
Definition: MahiFit.h:87
float itPulse[MaxSVSize]
Definition: MahiFit.h:92
float ootPulse[7][MaxSVSize]
Definition: MahiFit.h:93
constexpr bool hasTimeInfo() const
SampleVector amplitudes
Definition: MahiFit.h:31
void phase1Apply(const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &soiPlusOneEnergy, float &reconstructedTime, bool &useTriple, float &chi2) const
Definition: MahiFit.cc:46
constexpr double tsGain(const unsigned ts) const
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
static constexpr int pedestalBX_
Definition: MahiFit.h:161
PulseVector ampVec
Definition: MahiFit.h:56
float norm_
Definition: MahiFit.h:175
float count[MaxSVSize]
Definition: MahiFit.h:89
float totalUCNoise[MaxSVSize]
Definition: MahiFit.h:81
float inNoiseADC[MaxSVSize]
Definition: MahiFit.h:76
float pedEnergy
Definition: MahiFit.h:86
float inPedestal[MaxSVSize]
Definition: MahiFit.h:79
constexpr double tsPedestalWidth(const unsigned ts) const
T sqrt(T t)
Definition: SSEVec.h:19
int nSamples
Definition: MahiFit.h:66
float chiSq
Definition: MahiFit.h:84
constexpr unsigned soi() const
constexpr double tsPedestal(const unsigned ts) const
bool use3
Definition: MahiFit.h:69
constexpr float tsDFcPerADC(const unsigned ts) const
float inGain
Definition: MahiFit.h:74
float arrivalTime
Definition: MahiFit.h:85
float inPedAvg
Definition: MahiFit.h:73
float mahiEnergy
Definition: MahiFit.h:83
constexpr double fcByPE() const
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49
float inputTS[MaxSVSize]
Definition: MahiFit.h:90
int inputTDC[MaxSVSize]
Definition: MahiFit.h:91
float inNoisePhoto[MaxSVSize]
Definition: MahiFit.h:78
SampleVector noiseTerms
Definition: MahiFit.h:34
constexpr float tsRiseTime(const unsigned ts) const
unsigned int tsSize
Definition: MahiFit.h:21
BXVector bxs
Definition: MahiFit.h:28
float inTimeConst
Definition: MahiFit.h:71
constexpr double tsRawCharge(const unsigned ts) const
void MahiFit::resetPulseShapeTemplate ( int  pulseShapeId,
const HcalPulseShapes ps,
unsigned int  nSamples 
)
private

Definition at line 503 of file MahiFit.cc.

References cntsetPulseShape_, currentPulseShapeId_, HcalPulseShapes::getShape(), knownPulseShapes_, hcal::constants::maxSamples, AlCaHLTBitMon_ParallelJobs::p, and psfPtr_.

Referenced by setPulseShapeTemplate().

505  {
507 
508  psfPtr_ = nullptr;
509  for (auto& elem : knownPulseShapes_) {
510  if (elem.first == pulseShapeId) {
511  psfPtr_ = &*elem.second;
512  break;
513  }
514  }
515 
516  if (!psfPtr_) {
517  // only the pulse shape itself from PulseShapeFunctor is used for Mahi
518  // the uncertainty terms calculated inside PulseShapeFunctor are used for Method 2 only
519  auto p = std::make_shared<FitterFuncs::PulseShapeFunctor>(
520  ps.getShape(pulseShapeId), false, false, false, 1, 0, 0, hcal::constants::maxSamples);
521  knownPulseShapes_.emplace_back(pulseShapeId, p);
522  psfPtr_ = &*p;
523  }
524 
525  currentPulseShapeId_ = pulseShapeId;
526 }
const Shape & getShape(int shapeType) const
int currentPulseShapeId_
Definition: MahiFit.h:194
int cntsetPulseShape_
Definition: MahiFit.h:195
constexpr int maxSamples
Definition: HcalConstants.h:6
FitterFuncs::PulseShapeFunctor * psfPtr_
Definition: MahiFit.h:196
std::vector< ShapeWithId > knownPulseShapes_
Definition: MahiFit.h:197
void MahiFit::resetWorkspace ( ) const
private

Definition at line 667 of file MahiFit.cc.

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

Referenced by phase1Apply().

667  {
668  nnlsWork_.nPulseTot = 0;
669  nnlsWork_.tsOffset = 0;
670  nnlsWork_.bxOffset = 0;
671  nnlsWork_.maxoffset = 0;
672  nnlsWork_.nP = 0;
673 
674  // NOT SURE THIS IS NEEDED
675  nnlsWork_.amplitudes.setZero();
676  nnlsWork_.noiseTerms.setZero();
677  nnlsWork_.pedVals.setZero();
678 }
unsigned int nPulseTot
Definition: MahiFit.h:20
SampleVector amplitudes
Definition: MahiFit.h:31
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
unsigned int tsOffset
Definition: MahiFit.h:22
SampleVector noiseTerms
Definition: MahiFit.h:34
unsigned int nP
Definition: MahiFit.h:55
SampleVector pedVals
Definition: MahiFit.h:37
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.

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

20  {
21  dynamicPed_ = iDynamicPed;
22 
23  ts4Thresh_ = iTS4Thresh;
24  chiSqSwitch_ = chiSqSwitch;
25 
26  applyTimeSlew_ = iApplyTimeSlew;
27  slewFlavor_ = slewFlavor;
28 
29  calculateArrivalTime_ = iCalculateArrivalTime;
30  meanTime_ = iMeanTime;
31  timeSigmaHPD_ = iTimeSigmaHPD;
32  timeSigmaSiPM_ = iTimeSigmaSiPM;
33 
34  activeBXs_ = iActiveBXs;
35 
36  nMaxItersMin_ = iNMaxItersMin;
37  nMaxItersNNLS_ = iNMaxItersNNLS;
38 
39  deltaChiSqThresh_ = iDeltaChiSqThresh;
40  nnlsThresh_ = iNnlsThresh;
41 
42  bxOffsetConf_ = -(*std::min_element(activeBXs_.begin(), activeBXs_.end()));
43  bxSizeConf_ = activeBXs_.size();
44 }
float chiSqSwitch_
Definition: MahiFit.h:170
HcalTimeSlew::BiasSetting slewFlavor_
Definition: MahiFit.h:173
bool applyTimeSlew_
Definition: MahiFit.h:172
float meanTime_
Definition: MahiFit.h:178
bool calculateArrivalTime_
Definition: MahiFit.h:177
int bxOffsetConf_
Definition: MahiFit.h:191
float nnlsThresh_
Definition: MahiFit.h:188
unsigned int bxSizeConf_
Definition: MahiFit.h:190
int nMaxItersMin_
Definition: MahiFit.h:184
std::vector< int > activeBXs_
Definition: MahiFit.h:182
bool dynamicPed_
Definition: MahiFit.h:168
float timeSigmaSiPM_
Definition: MahiFit.h:180
float ts4Thresh_
Definition: MahiFit.h:169
float deltaChiSqThresh_
Definition: MahiFit.h:187
float timeSigmaHPD_
Definition: MahiFit.h:179
int nMaxItersNNLS_
Definition: MahiFit.h:185
void MahiFit::setPulseShapeTemplate ( int  pulseShapeId,
const HcalPulseShapes ps,
bool  hasTimeInfo,
const HcalTimeSlew hcalTimeSlewDelay,
unsigned int  nSamples 
)

Definition at line 477 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, cms::cuda::assert(), currentPulseShapeId_, HcalTimeSlew::delay(), MahiNnlsWorkspace::dt, hcalTimeSlewDelay_, nnlsWork_, MahiNnlsWorkspace::noiseTerms, hgc_digi::nSamples, MahiNnlsWorkspace::pedVals, resetPulseShapeTemplate(), slewFlavor_, timeSigmaHPD_, timeSigmaSiPM_, tsDelay1GeV_, and MahiNnlsWorkspace::tsSize.

481  {
482  if (hcalTimeSlewDelay != hcalTimeSlewDelay_) {
483  assert(hcalTimeSlewDelay);
484  hcalTimeSlewDelay_ = hcalTimeSlewDelay;
485  tsDelay1GeV_ = hcalTimeSlewDelay->delay(1.0, slewFlavor_);
486  }
487 
488  if (pulseShapeId != currentPulseShapeId_) {
489  resetPulseShapeTemplate(pulseShapeId, ps, nSamples);
490  }
491 
492  // 1 sigma time constraint
493  nnlsWork_.dt = hasTimeInfo ? timeSigmaSiPM_ : timeSigmaHPD_;
494 
495  if (nnlsWork_.tsSize != nSamples) {
497  nnlsWork_.amplitudes.resize(nSamples);
498  nnlsWork_.noiseTerms.resize(nSamples);
499  nnlsWork_.pedVals.resize(nSamples);
500  }
501 }
HcalTimeSlew::BiasSetting slewFlavor_
Definition: MahiFit.h:173
SampleVector amplitudes
Definition: MahiFit.h:31
int currentPulseShapeId_
Definition: MahiFit.h:194
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
assert(be >=bs)
constexpr size_t nSamples
float tsDelay1GeV_
Definition: MahiFit.h:174
void resetPulseShapeTemplate(int pulseShapeId, const HcalPulseShapes &ps, unsigned int nSamples)
Definition: MahiFit.cc:503
float timeSigmaSiPM_
Definition: MahiFit.h:180
SampleVector noiseTerms
Definition: MahiFit.h:34
const HcalTimeSlew * hcalTimeSlewDelay_
Definition: MahiFit.h:134
unsigned int tsSize
Definition: MahiFit.h:21
float timeSigmaHPD_
Definition: MahiFit.h:179
float delay(float fC, BiasSetting bias=Medium) const
Returns the amount (ns) by which a pulse of the given number of fC will be delayed by the timeslew ef...
Definition: HcalTimeSlew.cc:20
SampleVector pedVals
Definition: MahiFit.h:37
void MahiFit::solveSubmatrix ( PulseMatrix mat,
PulseVector invec,
PulseVector outvec,
unsigned  nP 
) const
private

Definition at line 618 of file MahiFit.cc.

References Exception, and groupFilesInBlocks::temp.

Referenced by nnls().

618  {
619  using namespace Eigen;
620  switch (nP) { // pulse matrix is always square.
621  case 10: {
622  Matrix<float, 10, 10> temp = mat;
623  outvec.head<10>() = temp.ldlt().solve(invec.head<10>());
624  } break;
625  case 9: {
626  Matrix<float, 9, 9> temp = mat.topLeftCorner<9, 9>();
627  outvec.head<9>() = temp.ldlt().solve(invec.head<9>());
628  } break;
629  case 8: {
630  Matrix<float, 8, 8> temp = mat.topLeftCorner<8, 8>();
631  outvec.head<8>() = temp.ldlt().solve(invec.head<8>());
632  } break;
633  case 7: {
634  Matrix<float, 7, 7> temp = mat.topLeftCorner<7, 7>();
635  outvec.head<7>() = temp.ldlt().solve(invec.head<7>());
636  } break;
637  case 6: {
638  Matrix<float, 6, 6> temp = mat.topLeftCorner<6, 6>();
639  outvec.head<6>() = temp.ldlt().solve(invec.head<6>());
640  } break;
641  case 5: {
642  Matrix<float, 5, 5> temp = mat.topLeftCorner<5, 5>();
643  outvec.head<5>() = temp.ldlt().solve(invec.head<5>());
644  } break;
645  case 4: {
646  Matrix<float, 4, 4> temp = mat.topLeftCorner<4, 4>();
647  outvec.head<4>() = temp.ldlt().solve(invec.head<4>());
648  } break;
649  case 3: {
650  Matrix<float, 3, 3> temp = mat.topLeftCorner<3, 3>();
651  outvec.head<3>() = temp.ldlt().solve(invec.head<3>());
652  } break;
653  case 2: {
654  Matrix<float, 2, 2> temp = mat.topLeftCorner<2, 2>();
655  outvec.head<2>() = temp.ldlt().solve(invec.head<2>());
656  } break;
657  case 1: {
658  Matrix<float, 1, 1> temp = mat.topLeftCorner<1, 1>();
659  outvec.head<1>() = temp.ldlt().solve(invec.head<1>());
660  } break;
661  default:
662  throw cms::Exception("HcalMahiWeirdState")
663  << "Weird number of pulses encountered in Mahi, module is configured incorrectly!";
664  }
665 }
void MahiFit::updateCov ( const SampleMatrix invCovMat) const
private

Definition at line 328 of file MahiFit.cc.

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().

328  {
329  SampleMatrix invCovMat = samplecov;
330 
331  for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
332  auto const amp = nnlsWork_.ampVec.coeff(iBX);
333  if (amp == 0)
334  continue;
335 
336  int offset = nnlsWork_.bxs.coeff(iBX);
337 
338  if (offset == pedestalBX_)
339  continue;
340  else {
341  invCovMat.noalias() += amp * amp * nnlsWork_.pulseCovArray[offset + nnlsWork_.bxOffset];
342  }
343  }
344 
345  nnlsWork_.covDecomp.compute(invCovMat);
346 }
unsigned int nPulseTot
Definition: MahiFit.h:20
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
static constexpr int pedestalBX_
Definition: MahiFit.h:161
PulseVector ampVec
Definition: MahiFit.h:56
std::array< SampleMatrix, MaxPVSize > pulseCovArray
Definition: MahiFit.h:46
SampleDecompLLT covDecomp
Definition: MahiFit.h:62
Eigen::Matrix< double, SampleVectorSize, SampleVectorSize > SampleMatrix
BXVector bxs
Definition: MahiFit.h:28
void MahiFit::updatePulseShape ( const float  itQ,
FullSampleVector pulseShape,
FullSampleVector pulseDeriv,
FullSampleMatrix pulseCov 
) const
private

Definition at line 268 of file MahiFit.cc.

References applyTimeSlew_, calculateArrivalTime_, HcalTimeSlew::delay(), CommonMethods::delta(), MahiNnlsWorkspace::dt, validate-o2o-wbm::f, FitterFuncs::PulseShapeFunctor::getPulseShape(), hcalTimeSlewDelay_, MahiNnlsWorkspace::maxoffset, meanTime_, nnlsWork_, psfPtr_, FitterFuncs::PulseShapeFunctor::singlePulseShapeFuncMahi(), slewFlavor_, FrontierCondition_GT_autoExpress_cfi::t0, createJobs::tmp, tsDelay1GeV_, MahiNnlsWorkspace::tsOffset, and MahiNnlsWorkspace::tsSize.

Referenced by doFit().

271  {
272  // set a null pulse shape for negative / or null TS
273  if (itQ <= 0.f)
274  return;
275 
276  float t0 = meanTime_;
277 
278  if (applyTimeSlew_) {
279  if (itQ <= 1.f)
280  t0 += tsDelay1GeV_;
281  else
282  t0 += hcalTimeSlewDelay_->delay(float(itQ), slewFlavor_);
283  }
284 
285  std::array<double, hcal::constants::maxSamples> pulseN;
286  std::array<double, hcal::constants::maxSamples> pulseM;
287  std::array<double, hcal::constants::maxSamples> pulseP;
288 
289  const float xx = t0;
290  const float xxm = -nnlsWork_.dt + t0;
291  const float xxp = nnlsWork_.dt + t0;
292 
294  psfPtr_->getPulseShape(pulseN);
295 
297  psfPtr_->getPulseShape(pulseM);
298 
300  psfPtr_->getPulseShape(pulseP);
301 
302  //in the 2018+ case where the sample of interest (SOI) is in TS3, add an extra offset to align
303  //with previous SOI=TS4 case assumed by psfPtr_->getPulseShape()
304  int delta = 4 - nnlsWork_.tsOffset;
305 
306  auto invDt = 0.5 / nnlsWork_.dt;
307 
308  for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
309  pulseShape(iTS + nnlsWork_.maxoffset) = pulseN[iTS + delta];
311  pulseDeriv(iTS + nnlsWork_.maxoffset) = (pulseM[iTS + delta] - pulseP[iTS + delta]) * invDt;
312 
313  pulseM[iTS + delta] -= pulseN[iTS + delta];
314  pulseP[iTS + delta] -= pulseN[iTS + delta];
315  }
316 
317  for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
318  for (unsigned int jTS = 0; jTS < iTS + 1; ++jTS) {
319  auto const tmp = 0.5 * (pulseP[iTS + delta] * pulseP[jTS + delta] + pulseM[iTS + delta] * pulseM[jTS + delta]);
320 
321  pulseCov(jTS + nnlsWork_.maxoffset, iTS + nnlsWork_.maxoffset) = tmp;
322  if (iTS != jTS)
323  pulseCov(iTS + nnlsWork_.maxoffset, jTS + nnlsWork_.maxoffset) = tmp;
324  }
325  }
326 }
HcalTimeSlew::BiasSetting slewFlavor_
Definition: MahiFit.h:173
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:158
void singlePulseShapeFuncMahi(const float *x)
bool applyTimeSlew_
Definition: MahiFit.h:172
float meanTime_
Definition: MahiFit.h:178
void getPulseShape(std::array< double, hcal::constants::maxSamples > &fillPulseShape)
bool calculateArrivalTime_
Definition: MahiFit.h:177
unsigned int tsOffset
Definition: MahiFit.h:22
float tsDelay1GeV_
Definition: MahiFit.h:174
const HcalTimeSlew * hcalTimeSlewDelay_
Definition: MahiFit.h:134
unsigned int tsSize
Definition: MahiFit.h:21
float delay(float fC, BiasSetting bias=Medium) const
Returns the amount (ns) by which a pulse of the given number of fC will be delayed by the timeslew ef...
Definition: HcalTimeSlew.cc:20
tmp
align.sh
Definition: createJobs.py:716
FitterFuncs::PulseShapeFunctor * psfPtr_
Definition: MahiFit.h:196

Member Data Documentation

std::vector<int> MahiFit::activeBXs_
private

Definition at line 182 of file MahiFit.h.

Referenced by doFit(), and setParameters().

bool MahiFit::applyTimeSlew_
private

Definition at line 172 of file MahiFit.h.

Referenced by setParameters(), and updatePulseShape().

int MahiFit::bxOffsetConf_
private

Definition at line 191 of file MahiFit.h.

Referenced by doFit(), and setParameters().

unsigned int MahiFit::bxSizeConf_
private

Definition at line 190 of file MahiFit.h.

Referenced by doFit(), and setParameters().

bool MahiFit::calculateArrivalTime_
private

Definition at line 177 of file MahiFit.h.

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

float MahiFit::chiSqSwitch_
private

Definition at line 170 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

int MahiFit::cntsetPulseShape_ = 0
private

Definition at line 195 of file MahiFit.h.

Referenced by resetPulseShapeTemplate().

int MahiFit::currentPulseShapeId_ = INT_MIN
private

Definition at line 194 of file MahiFit.h.

Referenced by resetPulseShapeTemplate(), and setPulseShapeTemplate().

float MahiFit::deltaChiSqThresh_
private

Definition at line 187 of file MahiFit.h.

Referenced by minimize(), and setParameters().

bool MahiFit::dynamicPed_
private

Definition at line 168 of file MahiFit.h.

Referenced by doFit(), and setParameters().

const HcalTimeSlew* MahiFit::hcalTimeSlewDelay_ = nullptr

Definition at line 134 of file MahiFit.h.

Referenced by setPulseShapeTemplate(), and updatePulseShape().

std::vector<ShapeWithId> MahiFit::knownPulseShapes_
private

Definition at line 197 of file MahiFit.h.

Referenced by resetPulseShapeTemplate().

float MahiFit::meanTime_
private

Definition at line 178 of file MahiFit.h.

Referenced by setParameters(), and updatePulseShape().

int MahiFit::nMaxItersMin_
private

Definition at line 184 of file MahiFit.h.

Referenced by minimize(), and setParameters().

int MahiFit::nMaxItersNNLS_
private

Definition at line 185 of file MahiFit.h.

Referenced by nnls(), and setParameters().

float MahiFit::nnlsThresh_
private

Definition at line 188 of file MahiFit.h.

Referenced by nnls(), and setParameters().

MahiNnlsWorkspace MahiFit::nnlsWork_
mutableprivate
float MahiFit::norm_ = (1.f / std::sqrt(12))
private

Definition at line 175 of file MahiFit.h.

Referenced by phase1Apply(), and phase1Debug().

constexpr int MahiFit::pedestalBX_ = 100
staticprivate

Definition at line 161 of file MahiFit.h.

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

FitterFuncs::PulseShapeFunctor* MahiFit::psfPtr_ = nullptr
private

Definition at line 196 of file MahiFit.h.

Referenced by resetPulseShapeTemplate(), and updatePulseShape().

HcalTimeSlew::BiasSetting MahiFit::slewFlavor_
private

Definition at line 173 of file MahiFit.h.

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

constexpr float MahiFit::timeLimit_ = 12.5f
staticprivate

Definition at line 165 of file MahiFit.h.

Referenced by calculateArrivalTime().

float MahiFit::timeSigmaHPD_
private

Definition at line 179 of file MahiFit.h.

Referenced by setParameters(), and setPulseShapeTemplate().

float MahiFit::timeSigmaSiPM_
private

Definition at line 180 of file MahiFit.h.

Referenced by setParameters(), and setPulseShapeTemplate().

float MahiFit::ts4Thresh_
private

Definition at line 169 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

float MahiFit::tsDelay1GeV_ = 0.f
private

Definition at line 174 of file MahiFit.h.

Referenced by setPulseShapeTemplate(), and updatePulseShape().