CMS 3D CMS Logo

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, int iTimeAlgo, double iThEnergeticPulses, 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, const float gain)
 
 ~MahiFit ()
 

Public Attributes

const HcalTimeSlewhcalTimeSlewDelay_ = nullptr
 
float thEnergeticPulses_
 
float thEnergeticPulsesFC_
 

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
 
float ccTime (const float itQ) 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_
 
int timeAlgo_
 
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

◆ Index

typedef BXVector::Index MahiFit::Index

Definition at line 136 of file MahiFit.h.

◆ ShapeWithId

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

Definition at line 143 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 99 of file MahiFit.h.

99 {};

Member Function Documentation

◆ calculateArrivalTime()

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

Definition at line 419 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().

419  {
420  if (nnlsWork_.nPulseTot > 1) {
421  SamplePulseMatrix pulseDerivMatTMP = nnlsWork_.pulseDerivMat;
422  for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
423  nnlsWork_.pulseDerivMat.col(iBX) = pulseDerivMatTMP.col(nnlsWork_.bxs.coeff(iBX) + nnlsWork_.bxOffset);
424  }
425  }
426 
427  for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
428  nnlsWork_.pulseDerivMat.col(iBX) *= nnlsWork_.ampVec.coeff(iBX);
429  }
430 
431  //FIXME: avoid solve full equation for one element
433  PulseVector solution = nnlsWork_.pulseDerivMat.colPivHouseholderQr().solve(residuals);
434  float t = std::clamp((float)solution.coeff(itIndex), -timeLimit_, timeLimit_);
435 
436  return t;
437 }
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:166
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:173
Eigen::Matrix< double, SampleVectorSize, 1 > SampleVector
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49
BXVector bxs
Definition: MahiFit.h:28

◆ calculateChiSq()

float MahiFit::calculateChiSq ( ) const
private

Definition at line 543 of file MahiFit.cc.

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

Referenced by minimize().

543  {
545  .squaredNorm();
546 }
SampleVector amplitudes
Definition: MahiFit.h:31
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:166
PulseVector ampVec
Definition: MahiFit.h:56
SampleDecompLLT covDecomp
Definition: MahiFit.h:62
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49

◆ ccTime()

float MahiFit::ccTime ( const float  itQ) const
private

Definition at line 354 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, applyTimeSlew_, HcalSpecialTimes::DEFAULT_ccTIME, HcalTimeSlew::delay(), dumpMFGeometry_cfg::delta, HLT_2023v12_cff::distance, f, FitterFuncs::PulseShapeFunctor::getPulseShape(), hcalTimeSlewDelay_, meanTime_, nnlsWork_, bTagMiniDQMDeepCSV::numerator, psfPtr_, FitterFuncs::PulseShapeFunctor::singlePulseShapeFuncMahi(), slewFlavor_, mathSSE::sqrt(), FrontierCondition_GT_autoExpress_cfi::t0, thEnergeticPulsesFC_, tsDelay1GeV_, MahiNnlsWorkspace::tsOffset, MahiNnlsWorkspace::tsSize, and geometryCSVtoXML::xx.

Referenced by doFit().

354  {
355  // those conditions are now on data time slices, can be done on the fitted pulse i.e. using nlsWork_.ampVec.coeff(itIndex);
356 
357  // Selecting energetic hits - Fitted Energy > 20 GeV
358  if (itQ < thEnergeticPulsesFC_)
360 
361  // Rejecting late hits Energy in TS[3] > (Energy in TS[4] and TS[5])
362  // With small OOTPU (Energy in TS[0] ,TS[1] and TS[2]) < 5 GeV
363  // To speed up check around the fitted time (?)
364 
365  // distanze as in formula of page 6
366  // https://indico.cern.ch/event/1142347/contributions/4793749/attachments/2412936/4129323/HCAL%20timing%20update.pdf
367 
368  float t0 = meanTime_;
369 
370  if (applyTimeSlew_) {
371  if (itQ <= 1.f)
372  t0 += tsDelay1GeV_;
373  else
374  t0 += hcalTimeSlewDelay_->delay(float(itQ), slewFlavor_);
375  }
376 
377  float ccTime = 0.f;
378  float distance_delta_max = 0.f;
379 
380  std::array<double, hcal::constants::maxSamples> pulseN;
381 
382  // making 8 TS out of the template i.e. 200 points
383  int TS_SOIandAfter = 25 * (nnlsWork_.tsSize - nnlsWork_.tsOffset);
384  int TS_beforeSOI = -25 * nnlsWork_.tsOffset;
385 
386  for (int deltaNS = TS_beforeSOI; deltaNS < TS_SOIandAfter; ++deltaNS) { // from -75ns and + 125ns
387  const float xx = t0 + deltaNS;
388 
390  psfPtr_->getPulseShape(pulseN);
391 
392  float pulse2 = 0;
393  float norm2 = 0;
394  float numerator = 0;
395  //
396  int delta = 4 - nnlsWork_.tsOffset; // like in updatePulseShape
397 
398  // rm TS0 and TS7, to speed up and reduce noise
399  for (unsigned int iTS = 1; iTS < (nnlsWork_.tsSize - 1); ++iTS) {
400  //pulseN[iTS] is the area of the template
401  float norm = nnlsWork_.amplitudes.coeffRef(iTS);
402 
403  // Finding the distance after each iteration.
404  numerator += norm * pulseN[iTS + delta];
405  pulse2 += pulseN[iTS + delta] * pulseN[iTS + delta];
406  norm2 += norm * norm;
407  }
408 
409  float distance = numerator / sqrt(pulse2 * norm2);
410  if (distance > distance_delta_max) {
411  distance_delta_max = distance;
412  ccTime = deltaNS;
413  }
414  }
415 
416  return ccTime;
417 }
float thEnergeticPulsesFC_
Definition: MahiFit.h:140
HcalTimeSlew::BiasSetting slewFlavor_
Definition: MahiFit.h:183
SampleVector amplitudes
Definition: MahiFit.h:31
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:166
void singlePulseShapeFuncMahi(const float *x)
bool applyTimeSlew_
Definition: MahiFit.h:182
float meanTime_
Definition: MahiFit.h:188
void getPulseShape(std::array< double, hcal::constants::maxSamples > &fillPulseShape)
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
T sqrt(T t)
Definition: SSEVec.h:19
double f[11][100]
unsigned int tsOffset
Definition: MahiFit.h:22
float ccTime(const float itQ) const
Definition: MahiFit.cc:354
float tsDelay1GeV_
Definition: MahiFit.h:184
const HcalTimeSlew * hcalTimeSlewDelay_
Definition: MahiFit.h:137
unsigned int tsSize
Definition: MahiFit.h:21
constexpr float DEFAULT_ccTIME
FitterFuncs::PulseShapeFunctor * psfPtr_
Definition: MahiFit.h:206

◆ doFit()

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

Definition at line 135 of file MahiFit.cc.

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

Referenced by phase1Apply().

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

◆ minimize()

const float MahiFit::minimize ( ) const
private

Definition at line 229 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().

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

◆ nnls()

void MahiFit::nnls ( ) const
private

Definition at line 439 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, particleFlowDisplacedVertex_cfi::ratio, solveSubmatrix(), remoteMonitoring_LASER_era2018_cfg::threshold, MahiNnlsWorkspace::tsSize, and cms::cuda::wmax.

Referenced by minimize().

439  {
440  const unsigned int npulse = nnlsWork_.nPulseTot;
441  const unsigned int nsamples = nnlsWork_.tsSize;
442 
443  PulseVector updateWork;
444  PulseVector ampvecpermtest;
445 
447  nnlsWork_.aTaMat.noalias() = nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.invcovp);
448  nnlsWork_.aTbVec.noalias() =
449  nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes));
450 
451  int iter = 0;
452  Index idxwmax = 0;
453  float wmax = 0.0f;
454  float threshold = nnlsThresh_;
455 
456  while (true) {
457  if (iter > 0 || nnlsWork_.nP == 0) {
458  if (nnlsWork_.nP == std::min(npulse, nsamples))
459  break;
460 
461  const unsigned int nActive = npulse - nnlsWork_.nP;
462  // exit if there are no more pulses to constrain
463  if (nActive == 0)
464  break;
465 
466  updateWork.noalias() = nnlsWork_.aTbVec - nnlsWork_.aTaMat.lazyProduct(nnlsWork_.ampVec);
467 
468  Index idxwmaxprev = idxwmax;
469  float wmaxprev = wmax;
470  wmax = updateWork.tail(nActive).maxCoeff(&idxwmax);
471 
472  if (wmax < threshold || (idxwmax == idxwmaxprev && wmax == wmaxprev)) {
473  break;
474  }
475 
476  if (iter >= nMaxItersNNLS_) {
477  break;
478  }
479 
480  //unconstrain parameter
481  idxwmax += nnlsWork_.nP;
482  nnlsUnconstrainParameter(idxwmax);
483  }
484 
485  while (true) {
486  if (nnlsWork_.nP == 0)
487  break;
488 
489  ampvecpermtest.noalias() = nnlsWork_.ampVec.head(nnlsWork_.nP);
490 
492 
493  //check solution
494  if (ampvecpermtest.head(nnlsWork_.nP).minCoeff() > 0.f) {
495  nnlsWork_.ampVec.head(nnlsWork_.nP) = ampvecpermtest.head(nnlsWork_.nP);
496  break;
497  }
498 
499  //update parameter vector
500  Index minratioidx = 0;
501 
502  // no realizable optimization here (because it autovectorizes!)
503  float minratio = std::numeric_limits<float>::max();
504  for (unsigned int ipulse = 0; ipulse < nnlsWork_.nP; ++ipulse) {
505  if (ampvecpermtest.coeff(ipulse) <= 0.f) {
506  const float c_ampvec = nnlsWork_.ampVec.coeff(ipulse);
507  const float ratio = c_ampvec / (c_ampvec - ampvecpermtest.coeff(ipulse));
508  if (ratio < minratio) {
509  minratio = ratio;
510  minratioidx = ipulse;
511  }
512  }
513  }
514  nnlsWork_.ampVec.head(nnlsWork_.nP) +=
515  minratio * (ampvecpermtest.head(nnlsWork_.nP) - nnlsWork_.ampVec.head(nnlsWork_.nP));
516 
517  //avoid numerical problems with later ==0. check
518  nnlsWork_.ampVec.coeffRef(minratioidx) = 0.f;
519 
520  nnlsConstrainParameter(minratioidx);
521  }
522 
523  ++iter;
524 
525  //adaptive convergence threshold to avoid infinite loops but still
526  //ensure best value is used
527  if (iter % 10 == 0) {
528  threshold *= 10.;
529  }
530  }
531 }
unsigned int nPulseTot
Definition: MahiFit.h:20
void nnlsUnconstrainParameter(Index idxp) const
Definition: MahiFit.cc:603
SamplePulseMatrix invcovp
Definition: MahiFit.h:58
SampleVector amplitudes
Definition: MahiFit.h:31
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:166
void nnlsConstrainParameter(Index minratioidx) const
Definition: MahiFit.cc:615
BXVector::Index Index
Definition: MahiFit.h:136
PulseVector ampVec
Definition: MahiFit.h:56
SampleDecompLLT covDecomp
Definition: MahiFit.h:62
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, PulseVectorSize, 1 > PulseVector
PulseVector aTbVec
Definition: MahiFit.h:60
float nnlsThresh_
Definition: MahiFit.h:198
void solveSubmatrix(PulseMatrix &mat, PulseVector &invec, PulseVector &outvec, unsigned nP) const
Definition: MahiFit.cc:693
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:195
__host__ __device__ V V wmax

◆ nnlsConstrainParameter()

void MahiFit::nnlsConstrainParameter ( Index  minratioidx) const
private

Definition at line 615 of file MahiFit.cc.

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

Referenced by nnls().

615  {
616  if (minratioidx != (nnlsWork_.nP - 1)) {
617  nnlsWork_.aTaMat.col(nnlsWork_.nP - 1).swap(nnlsWork_.aTaMat.col(minratioidx));
618  nnlsWork_.aTaMat.row(nnlsWork_.nP - 1).swap(nnlsWork_.aTaMat.row(minratioidx));
619  nnlsWork_.pulseMat.col(nnlsWork_.nP - 1).swap(nnlsWork_.pulseMat.col(minratioidx));
620  Eigen::numext::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP - 1), nnlsWork_.aTbVec.coeffRef(minratioidx));
621  Eigen::numext::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP - 1), nnlsWork_.ampVec.coeffRef(minratioidx));
622  Eigen::numext::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP - 1), nnlsWork_.bxs.coeffRef(minratioidx));
623  }
624  --nnlsWork_.nP;
625 }
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:166
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

◆ nnlsUnconstrainParameter()

void MahiFit::nnlsUnconstrainParameter ( Index  idxp) const
private

Definition at line 603 of file MahiFit.cc.

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

Referenced by nnls().

603  {
604  if (idxp != nnlsWork_.nP) {
605  nnlsWork_.aTaMat.col(nnlsWork_.nP).swap(nnlsWork_.aTaMat.col(idxp));
606  nnlsWork_.aTaMat.row(nnlsWork_.nP).swap(nnlsWork_.aTaMat.row(idxp));
607  nnlsWork_.pulseMat.col(nnlsWork_.nP).swap(nnlsWork_.pulseMat.col(idxp));
608  Eigen::numext::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP), nnlsWork_.aTbVec.coeffRef(idxp));
609  Eigen::numext::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP), nnlsWork_.ampVec.coeffRef(idxp));
610  Eigen::numext::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP), nnlsWork_.bxs.coeffRef(idxp));
611  }
612  ++nnlsWork_.nP;
613 }
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:166
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

◆ onePulseMinimize()

void MahiFit::onePulseMinimize ( ) const
private

Definition at line 533 of file MahiFit.cc.

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

Referenced by minimize().

533  {
535 
536  float aTaCoeff = (nnlsWork_.invcovp.transpose() * nnlsWork_.invcovp).coeff(0);
537  float aTbCoeff =
538  (nnlsWork_.invcovp.transpose() * (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes))).coeff(0);
539 
540  nnlsWork_.ampVec.coeffRef(0) = std::max(0.f, aTbCoeff / aTaCoeff);
541 }
SamplePulseMatrix invcovp
Definition: MahiFit.h:58
SampleVector amplitudes
Definition: MahiFit.h:31
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:166
PulseVector ampVec
Definition: MahiFit.h:56
SampleDecompLLT covDecomp
Definition: MahiFit.h:62
double f[11][100]
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49

◆ phase1Apply()

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

Definition at line 50 of file MahiFit.cc.

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

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

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

◆ phase1Debug()

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

Definition at line 627 of file MahiFit.cc.

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.

627  {
628  float recoEnergy, soiPlus1Energy, recoTime, chi2;
629  bool use3;
630  phase1Apply(channelData, recoEnergy, soiPlus1Energy, recoTime, use3, chi2);
631 
632  mdi.nSamples = channelData.nSamples();
633  mdi.soi = channelData.soi();
634 
635  mdi.use3 = use3;
636 
637  mdi.inTimeConst = nnlsWork_.dt;
638  mdi.inPedAvg = 0.25 * (channelData.tsPedestalWidth(0) * channelData.tsPedestalWidth(0) +
639  channelData.tsPedestalWidth(1) * channelData.tsPedestalWidth(1) +
640  channelData.tsPedestalWidth(2) * channelData.tsPedestalWidth(2) +
641  channelData.tsPedestalWidth(3) * channelData.tsPedestalWidth(3));
642  mdi.inGain = channelData.tsGain(0);
643 
644  for (unsigned int iTS = 0; iTS < channelData.nSamples(); ++iTS) {
645  double charge = channelData.tsRawCharge(iTS);
646  double ped = channelData.tsPedestal(iTS);
647 
648  mdi.inNoiseADC[iTS] = norm_ * channelData.tsDFcPerADC(iTS);
649 
650  if ((charge - ped) > channelData.tsPedestalWidth(iTS)) {
651  mdi.inNoisePhoto[iTS] = sqrt((charge - ped) * channelData.fcByPE());
652  } else {
653  mdi.inNoisePhoto[iTS] = 0;
654  }
655 
656  mdi.inPedestal[iTS] = channelData.tsPedestalWidth(iTS);
657  mdi.totalUCNoise[iTS] = nnlsWork_.noiseTerms.coeffRef(iTS);
658 
659  if (channelData.hasTimeInfo()) {
660  mdi.inputTDC[iTS] = channelData.tsRiseTime(iTS);
661  } else {
662  mdi.inputTDC[iTS] = -1;
663  }
664  }
665 
666  mdi.arrivalTime = recoTime;
667  mdi.chiSq = chi2;
668 
669  for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
670  if (nnlsWork_.bxs.coeff(iBX) == 0) {
671  mdi.mahiEnergy = nnlsWork_.ampVec.coeff(iBX);
672  for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
673  mdi.count[iTS] = iTS;
674  mdi.inputTS[iTS] = nnlsWork_.amplitudes.coeff(iTS);
675  mdi.itPulse[iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
676  }
677  } else if (nnlsWork_.bxs.coeff(iBX) == pedestalBX_) {
678  mdi.pedEnergy = nnlsWork_.ampVec.coeff(iBX);
679  } else if (nnlsWork_.bxs.coeff(iBX) >= -3 && nnlsWork_.bxs.coeff(iBX) <= 4) {
680  int ootIndex = nnlsWork_.bxs.coeff(iBX);
681  if (ootIndex > 0)
682  ootIndex += 2;
683  else
684  ootIndex += 3;
685  mdi.ootEnergy[ootIndex] = nnlsWork_.ampVec.coeff(iBX);
686  for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
687  mdi.ootPulse[ootIndex][iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
688  }
689  }
690  }
691 }
unsigned int nPulseTot
Definition: MahiFit.h:20
float ootEnergy[7]
Definition: MahiFit.h:87
constexpr double fcByPE() const
float itPulse[MaxSVSize]
Definition: MahiFit.h:92
float ootPulse[7][MaxSVSize]
Definition: MahiFit.h:93
SampleVector amplitudes
Definition: MahiFit.h:31
constexpr float tsDFcPerADC(const unsigned ts) const
constexpr double tsRawCharge(const unsigned ts) const
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:166
constexpr float tsRiseTime(const unsigned ts) const
static constexpr int pedestalBX_
Definition: MahiFit.h:169
PulseVector ampVec
Definition: MahiFit.h:56
float norm_
Definition: MahiFit.h:185
constexpr double tsPedestalWidth(const unsigned ts) const
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
void phase1Apply(const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &soiPlusOneEnergy, float &reconstructedTime, bool &useTriple, float &chi2) const
Definition: MahiFit.cc:50
float inPedestal[MaxSVSize]
Definition: MahiFit.h:79
T sqrt(T t)
Definition: SSEVec.h:19
constexpr double tsGain(const unsigned ts) const
int nSamples
Definition: MahiFit.h:66
float chiSq
Definition: MahiFit.h:84
constexpr unsigned soi() const
bool use3
Definition: MahiFit.h:69
float inGain
Definition: MahiFit.h:74
float arrivalTime
Definition: MahiFit.h:85
float inPedAvg
Definition: MahiFit.h:73
float mahiEnergy
Definition: MahiFit.h:83
SamplePulseMatrix pulseMat
Definition: MahiFit.h:49
constexpr bool hasTimeInfo() const
constexpr unsigned nSamples() const
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 double tsPedestal(const unsigned ts) const
unsigned int tsSize
Definition: MahiFit.h:21
BXVector bxs
Definition: MahiFit.h:28
float inTimeConst
Definition: MahiFit.h:71

◆ resetPulseShapeTemplate()

void MahiFit::resetPulseShapeTemplate ( int  pulseShapeId,
const HcalPulseShapes ps,
unsigned int  nSamples 
)
private

Definition at line 578 of file MahiFit.cc.

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

Referenced by setPulseShapeTemplate().

580  {
582 
583  psfPtr_ = nullptr;
584  for (auto& elem : knownPulseShapes_) {
585  if (elem.first == pulseShapeId) {
586  psfPtr_ = &*elem.second;
587  break;
588  }
589  }
590 
591  if (!psfPtr_) {
592  // only the pulse shape itself from PulseShapeFunctor is used for Mahi
593  // the uncertainty terms calculated inside PulseShapeFunctor are used for Method 2 only
594  auto p = std::make_shared<FitterFuncs::PulseShapeFunctor>(
595  ps.getShape(pulseShapeId), false, false, false, 1, 0, 0, hcal::constants::maxSamples);
596  knownPulseShapes_.emplace_back(pulseShapeId, p);
597  psfPtr_ = &*p;
598  }
599 
600  currentPulseShapeId_ = pulseShapeId;
601 }
int currentPulseShapeId_
Definition: MahiFit.h:204
int cntsetPulseShape_
Definition: MahiFit.h:205
const Shape & getShape(int shapeType) const
constexpr int maxSamples
Definition: HcalConstants.h:6
FitterFuncs::PulseShapeFunctor * psfPtr_
Definition: MahiFit.h:206
std::vector< ShapeWithId > knownPulseShapes_
Definition: MahiFit.h:207

◆ resetWorkspace()

void MahiFit::resetWorkspace ( ) const
private

Definition at line 742 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().

742  {
743  nnlsWork_.nPulseTot = 0;
744  nnlsWork_.tsOffset = 0;
745  nnlsWork_.bxOffset = 0;
746  nnlsWork_.maxoffset = 0;
747  nnlsWork_.nP = 0;
748 
749  // NOT SURE THIS IS NEEDED
750  nnlsWork_.amplitudes.setZero();
751  nnlsWork_.noiseTerms.setZero();
752  nnlsWork_.pedVals.setZero();
753 }
unsigned int nPulseTot
Definition: MahiFit.h:20
SampleVector amplitudes
Definition: MahiFit.h:31
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:166
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

◆ setParameters()

void MahiFit::setParameters ( bool  iDynamicPed,
double  iTS4Thresh,
double  chiSqSwitch,
bool  iApplyTimeSlew,
HcalTimeSlew::BiasSetting  slewFlavor,
bool  iCalculateArrivalTime,
int  iTimeAlgo,
double  iThEnergeticPulses,
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_, HLT_2023v12_cff::chiSqSwitch, chiSqSwitch_, deltaChiSqThresh_, dynamicPed_, meanTime_, nMaxItersMin_, nMaxItersNNLS_, nnlsThresh_, slewFlavor_, thEnergeticPulses_, HBHEMahiParameters_cfi::timeAlgo, timeAlgo_, timeSigmaHPD_, timeSigmaSiPM_, and ts4Thresh_.

22  {
23  dynamicPed_ = iDynamicPed;
24 
25  ts4Thresh_ = iTS4Thresh;
27 
28  applyTimeSlew_ = iApplyTimeSlew;
29  slewFlavor_ = slewFlavor;
30 
31  calculateArrivalTime_ = iCalculateArrivalTime;
33  thEnergeticPulses_ = iThEnergeticPulses;
34  meanTime_ = iMeanTime;
35  timeSigmaHPD_ = iTimeSigmaHPD;
36  timeSigmaSiPM_ = iTimeSigmaSiPM;
37 
38  activeBXs_ = iActiveBXs;
39 
40  nMaxItersMin_ = iNMaxItersMin;
41  nMaxItersNNLS_ = iNMaxItersNNLS;
42 
43  deltaChiSqThresh_ = iDeltaChiSqThresh;
44  nnlsThresh_ = iNnlsThresh;
45 
46  bxOffsetConf_ = -(*std::min_element(activeBXs_.begin(), activeBXs_.end()));
47  bxSizeConf_ = activeBXs_.size();
48 }
float chiSqSwitch_
Definition: MahiFit.h:180
HcalTimeSlew::BiasSetting slewFlavor_
Definition: MahiFit.h:183
bool applyTimeSlew_
Definition: MahiFit.h:182
float meanTime_
Definition: MahiFit.h:188
float thEnergeticPulses_
Definition: MahiFit.h:139
bool calculateArrivalTime_
Definition: MahiFit.h:187
int bxOffsetConf_
Definition: MahiFit.h:201
float nnlsThresh_
Definition: MahiFit.h:198
unsigned int bxSizeConf_
Definition: MahiFit.h:200
int nMaxItersMin_
Definition: MahiFit.h:194
std::vector< int > activeBXs_
Definition: MahiFit.h:192
bool dynamicPed_
Definition: MahiFit.h:178
float timeSigmaSiPM_
Definition: MahiFit.h:190
float ts4Thresh_
Definition: MahiFit.h:179
float deltaChiSqThresh_
Definition: MahiFit.h:197
float timeSigmaHPD_
Definition: MahiFit.h:189
int nMaxItersNNLS_
Definition: MahiFit.h:195
int timeAlgo_
Definition: MahiFit.h:176

◆ setPulseShapeTemplate()

void MahiFit::setPulseShapeTemplate ( int  pulseShapeId,
const HcalPulseShapes ps,
bool  hasTimeInfo,
const HcalTimeSlew hcalTimeSlewDelay,
unsigned int  nSamples,
const float  gain 
)

Definition at line 548 of file MahiFit.cc.

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

553  {
554  if (hcalTimeSlewDelay != hcalTimeSlewDelay_) {
555  assert(hcalTimeSlewDelay);
556  hcalTimeSlewDelay_ = hcalTimeSlewDelay;
557  tsDelay1GeV_ = hcalTimeSlewDelay->delay(1.0, slewFlavor_);
558  }
559 
560  if (pulseShapeId != currentPulseShapeId_) {
561  resetPulseShapeTemplate(pulseShapeId, ps, nSamples);
562  }
563 
564  // 1 sigma time constraint
565  nnlsWork_.dt = hasTimeInfo ? timeSigmaSiPM_ : timeSigmaHPD_;
566 
567  if (nnlsWork_.tsSize != nSamples) {
569  nnlsWork_.amplitudes.resize(nSamples);
570  nnlsWork_.noiseTerms.resize(nSamples);
571  nnlsWork_.pedVals.resize(nSamples);
572  }
573 
574  // threshold in GeV for ccTime
576 }
float thEnergeticPulsesFC_
Definition: MahiFit.h:140
HcalTimeSlew::BiasSetting slewFlavor_
Definition: MahiFit.h:183
SampleVector amplitudes
Definition: MahiFit.h:31
int currentPulseShapeId_
Definition: MahiFit.h:204
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:166
assert(be >=bs)
float thEnergeticPulses_
Definition: MahiFit.h:139
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
float tsDelay1GeV_
Definition: MahiFit.h:184
void resetPulseShapeTemplate(int pulseShapeId, const HcalPulseShapes &ps, unsigned int nSamples)
Definition: MahiFit.cc:578
float timeSigmaSiPM_
Definition: MahiFit.h:190
SampleVector noiseTerms
Definition: MahiFit.h:34
const HcalTimeSlew * hcalTimeSlewDelay_
Definition: MahiFit.h:137
unsigned int tsSize
Definition: MahiFit.h:21
float timeSigmaHPD_
Definition: MahiFit.h:189
SampleVector pedVals
Definition: MahiFit.h:37

◆ solveSubmatrix()

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

Definition at line 693 of file MahiFit.cc.

References Exception, and groupFilesInBlocks::temp.

Referenced by nnls().

693  {
694  using namespace Eigen;
695  switch (nP) { // pulse matrix is always square.
696  case 10: {
697  Matrix<float, 10, 10> temp = mat;
698  outvec.head<10>() = temp.ldlt().solve(invec.head<10>());
699  } break;
700  case 9: {
701  Matrix<float, 9, 9> temp = mat.topLeftCorner<9, 9>();
702  outvec.head<9>() = temp.ldlt().solve(invec.head<9>());
703  } break;
704  case 8: {
705  Matrix<float, 8, 8> temp = mat.topLeftCorner<8, 8>();
706  outvec.head<8>() = temp.ldlt().solve(invec.head<8>());
707  } break;
708  case 7: {
709  Matrix<float, 7, 7> temp = mat.topLeftCorner<7, 7>();
710  outvec.head<7>() = temp.ldlt().solve(invec.head<7>());
711  } break;
712  case 6: {
713  Matrix<float, 6, 6> temp = mat.topLeftCorner<6, 6>();
714  outvec.head<6>() = temp.ldlt().solve(invec.head<6>());
715  } break;
716  case 5: {
717  Matrix<float, 5, 5> temp = mat.topLeftCorner<5, 5>();
718  outvec.head<5>() = temp.ldlt().solve(invec.head<5>());
719  } break;
720  case 4: {
721  Matrix<float, 4, 4> temp = mat.topLeftCorner<4, 4>();
722  outvec.head<4>() = temp.ldlt().solve(invec.head<4>());
723  } break;
724  case 3: {
725  Matrix<float, 3, 3> temp = mat.topLeftCorner<3, 3>();
726  outvec.head<3>() = temp.ldlt().solve(invec.head<3>());
727  } break;
728  case 2: {
729  Matrix<float, 2, 2> temp = mat.topLeftCorner<2, 2>();
730  outvec.head<2>() = temp.ldlt().solve(invec.head<2>());
731  } break;
732  case 1: {
733  Matrix<float, 1, 1> temp = mat.topLeftCorner<1, 1>();
734  outvec.head<1>() = temp.ldlt().solve(invec.head<1>());
735  } break;
736  default:
737  throw cms::Exception("HcalMahiWeirdState")
738  << "Weird number of pulses encountered in Mahi, module is configured incorrectly!";
739  }
740 }

◆ updateCov()

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

Definition at line 334 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().

334  {
335  SampleMatrix invCovMat = samplecov;
336 
337  for (unsigned int iBX = 0; iBX < nnlsWork_.nPulseTot; ++iBX) {
338  auto const amp = nnlsWork_.ampVec.coeff(iBX);
339  if (amp == 0)
340  continue;
341 
342  int offset = nnlsWork_.bxs.coeff(iBX);
343 
344  if (offset == pedestalBX_)
345  continue;
346  else {
347  invCovMat.noalias() += amp * amp * nnlsWork_.pulseCovArray[offset + nnlsWork_.bxOffset];
348  }
349  }
350 
351  nnlsWork_.covDecomp.compute(invCovMat);
352 }
unsigned int nPulseTot
Definition: MahiFit.h:20
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:166
static constexpr int pedestalBX_
Definition: MahiFit.h:169
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

◆ updatePulseShape()

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

Definition at line 274 of file MahiFit.cc.

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

Referenced by doFit().

277  {
278  // set a null pulse shape for negative / or null TS
279  if (itQ <= 0.f)
280  return;
281 
282  float t0 = meanTime_;
283 
284  if (applyTimeSlew_) {
285  if (itQ <= 1.f)
286  t0 += tsDelay1GeV_;
287  else
288  t0 += hcalTimeSlewDelay_->delay(float(itQ), slewFlavor_);
289  }
290 
291  std::array<double, hcal::constants::maxSamples> pulseN;
292  std::array<double, hcal::constants::maxSamples> pulseM;
293  std::array<double, hcal::constants::maxSamples> pulseP;
294 
295  const float xx = t0;
296  const float xxm = -nnlsWork_.dt + t0;
297  const float xxp = nnlsWork_.dt + t0;
298 
300  psfPtr_->getPulseShape(pulseN);
301 
303  psfPtr_->getPulseShape(pulseM);
304 
306  psfPtr_->getPulseShape(pulseP);
307 
308  //in the 2018+ case where the sample of interest (SOI) is in TS3, add an extra offset to align
309  //with previous SOI=TS4 case assumed by psfPtr_->getPulseShape()
310  int delta = 4 - nnlsWork_.tsOffset;
311 
312  auto invDt = 0.5 / nnlsWork_.dt;
313 
314  for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
315  pulseShape(iTS + nnlsWork_.maxoffset) = pulseN[iTS + delta];
317  pulseDeriv(iTS + nnlsWork_.maxoffset) = (pulseM[iTS + delta] - pulseP[iTS + delta]) * invDt;
318 
319  pulseM[iTS + delta] -= pulseN[iTS + delta];
320  pulseP[iTS + delta] -= pulseN[iTS + delta];
321  }
322 
323  for (unsigned int iTS = 0; iTS < nnlsWork_.tsSize; ++iTS) {
324  for (unsigned int jTS = 0; jTS < iTS + 1; ++jTS) {
325  auto const tmp = 0.5 * (pulseP[iTS + delta] * pulseP[jTS + delta] + pulseM[iTS + delta] * pulseM[jTS + delta]);
326 
327  pulseCov(jTS + nnlsWork_.maxoffset, iTS + nnlsWork_.maxoffset) = tmp;
328  if (iTS != jTS)
329  pulseCov(iTS + nnlsWork_.maxoffset, jTS + nnlsWork_.maxoffset) = tmp;
330  }
331  }
332 }
HcalTimeSlew::BiasSetting slewFlavor_
Definition: MahiFit.h:183
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:166
void singlePulseShapeFuncMahi(const float *x)
bool applyTimeSlew_
Definition: MahiFit.h:182
float meanTime_
Definition: MahiFit.h:188
void getPulseShape(std::array< double, hcal::constants::maxSamples > &fillPulseShape)
bool calculateArrivalTime_
Definition: MahiFit.h:187
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
double f[11][100]
unsigned int tsOffset
Definition: MahiFit.h:22
float tsDelay1GeV_
Definition: MahiFit.h:184
const HcalTimeSlew * hcalTimeSlewDelay_
Definition: MahiFit.h:137
unsigned int tsSize
Definition: MahiFit.h:21
tmp
align.sh
Definition: createJobs.py:716
FitterFuncs::PulseShapeFunctor * psfPtr_
Definition: MahiFit.h:206

Member Data Documentation

◆ activeBXs_

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

Definition at line 192 of file MahiFit.h.

Referenced by doFit(), and setParameters().

◆ applyTimeSlew_

bool MahiFit::applyTimeSlew_
private

Definition at line 182 of file MahiFit.h.

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

◆ bxOffsetConf_

int MahiFit::bxOffsetConf_
private

Definition at line 201 of file MahiFit.h.

Referenced by doFit(), and setParameters().

◆ bxSizeConf_

unsigned int MahiFit::bxSizeConf_
private

Definition at line 200 of file MahiFit.h.

Referenced by doFit(), and setParameters().

◆ calculateArrivalTime_

bool MahiFit::calculateArrivalTime_
private

Definition at line 187 of file MahiFit.h.

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

◆ chiSqSwitch_

float MahiFit::chiSqSwitch_
private

Definition at line 180 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

◆ cntsetPulseShape_

int MahiFit::cntsetPulseShape_ = 0
private

Definition at line 205 of file MahiFit.h.

Referenced by resetPulseShapeTemplate().

◆ currentPulseShapeId_

int MahiFit::currentPulseShapeId_ = INT_MIN
private

Definition at line 204 of file MahiFit.h.

Referenced by resetPulseShapeTemplate(), and setPulseShapeTemplate().

◆ deltaChiSqThresh_

float MahiFit::deltaChiSqThresh_
private

Definition at line 197 of file MahiFit.h.

Referenced by minimize(), and setParameters().

◆ dynamicPed_

bool MahiFit::dynamicPed_
private

Definition at line 178 of file MahiFit.h.

Referenced by doFit(), and setParameters().

◆ hcalTimeSlewDelay_

const HcalTimeSlew* MahiFit::hcalTimeSlewDelay_ = nullptr

Definition at line 137 of file MahiFit.h.

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

◆ knownPulseShapes_

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

Definition at line 207 of file MahiFit.h.

Referenced by resetPulseShapeTemplate().

◆ meanTime_

float MahiFit::meanTime_
private

Definition at line 188 of file MahiFit.h.

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

◆ nMaxItersMin_

int MahiFit::nMaxItersMin_
private

Definition at line 194 of file MahiFit.h.

Referenced by minimize(), and setParameters().

◆ nMaxItersNNLS_

int MahiFit::nMaxItersNNLS_
private

Definition at line 195 of file MahiFit.h.

Referenced by nnls(), and setParameters().

◆ nnlsThresh_

float MahiFit::nnlsThresh_
private

Definition at line 198 of file MahiFit.h.

Referenced by nnls(), and setParameters().

◆ nnlsWork_

MahiNnlsWorkspace MahiFit::nnlsWork_
mutableprivate

◆ norm_

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

Definition at line 185 of file MahiFit.h.

Referenced by phase1Apply(), and phase1Debug().

◆ pedestalBX_

constexpr int MahiFit::pedestalBX_ = 100
staticprivate

Definition at line 169 of file MahiFit.h.

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

◆ psfPtr_

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

Definition at line 206 of file MahiFit.h.

Referenced by ccTime(), resetPulseShapeTemplate(), and updatePulseShape().

◆ slewFlavor_

HcalTimeSlew::BiasSetting MahiFit::slewFlavor_
private

Definition at line 183 of file MahiFit.h.

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

◆ thEnergeticPulses_

float MahiFit::thEnergeticPulses_

Definition at line 139 of file MahiFit.h.

Referenced by setParameters(), and setPulseShapeTemplate().

◆ thEnergeticPulsesFC_

float MahiFit::thEnergeticPulsesFC_

Definition at line 140 of file MahiFit.h.

Referenced by ccTime(), and setPulseShapeTemplate().

◆ timeAlgo_

int MahiFit::timeAlgo_
private

Definition at line 176 of file MahiFit.h.

Referenced by doFit(), and setParameters().

◆ timeLimit_

constexpr float MahiFit::timeLimit_ = 12.5f
staticprivate

Definition at line 173 of file MahiFit.h.

Referenced by calculateArrivalTime().

◆ timeSigmaHPD_

float MahiFit::timeSigmaHPD_
private

Definition at line 189 of file MahiFit.h.

Referenced by setParameters(), and setPulseShapeTemplate().

◆ timeSigmaSiPM_

float MahiFit::timeSigmaSiPM_
private

Definition at line 190 of file MahiFit.h.

Referenced by setParameters(), and setPulseShapeTemplate().

◆ ts4Thresh_

float MahiFit::ts4Thresh_
private

Definition at line 179 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

◆ tsDelay1GeV_

float MahiFit::tsDelay1GeV_ = 0.f
private

Definition at line 184 of file MahiFit.h.

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