CMS 3D CMS Logo

List of all members | Public Types | Public Member Functions | Public Attributes | 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, 3 > &correctedOutput, const int nbx) const
 
 MahiFit ()
 
void phase1Apply (const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &reconstructedTime, bool &useTriple, float &chi2) const
 
void phase1Debug (const HBHEChannelInfo &channelData, MahiDebugInfo &mdi) const
 
void resetPulseShapeTemplate (const HcalPulseShapes::Shape &ps)
 
void setParameters (bool iDynamicPed, double iTS4Thresh, double chiSqSwitch, bool iApplyTimeSlew, HcalTimeSlew::BiasSetting slewFlavor, double iMeanTime, double iTimeSigmaHPD, double iTimeSigmaSiPM, const std::vector< int > &iActiveBXs, int iNMaxItersMin, int iNMaxItersNNLS, double iDeltaChiSqThresh, double iNnlsThresh)
 
void setPulseShapeTemplate (const HcalPulseShapes::Shape &ps, const HcalTimeSlew *hcalTimeSlewDelay)
 
 ~MahiFit ()
 

Public Attributes

const HcalPulseShapes::ShapecurrentPulseShape_ = 0
 
const HcalTimeSlewhcalTimeSlewDelay_ = 0
 

Private Member Functions

double calculateArrivalTime () const
 
double calculateChiSq () const
 
double minimize () const
 
void nnls () const
 
void nnlsConstrainParameter (Index minratioidx) const
 
void nnlsUnconstrainParameter (Index idxp) const
 
void onePulseMinimize () const
 
void resetWorkspace () const
 
void solveSubmatrix (PulseMatrix &mat, PulseVector &invec, PulseVector &outvec, unsigned nP) const
 
void updateCov () const
 
void updatePulseShape (double itQ, FullSampleVector &pulseShape, FullSampleVector &pulseDeriv, FullSampleMatrix &pulseCov) const
 

Private Attributes

std::vector< int > activeBXs_
 
bool applyTimeSlew_
 
int bxOffsetConf_
 
unsigned int bxSizeConf_
 
float chiSqSwitch_
 
int cntsetPulseShape_
 
float deltaChiSqThresh_
 
bool dynamicPed_
 
const unsigned int fullTSofInterest_
 
const unsigned int fullTSSize_
 
float meanTime_
 
int nMaxItersMin_
 
int nMaxItersNNLS_
 
float nnlsThresh_
 
MahiNnlsWorkspace nnlsWork_
 
std::unique_ptr< ROOT::Math::Functor > pfunctor_
 
std::unique_ptr< FitterFuncs::PulseShapeFunctorpsfPtr_
 
HcalTimeSlew::BiasSetting slewFlavor_
 
float timeSigmaHPD_
 
float timeSigmaSiPM_
 
float ts4Thresh_
 
double tsDelay1GeV_ =0
 

Static Private Attributes

static int pedestalBX_ = 100
 
static float timeLimit_ = 12.5
 

Detailed Description

Definition at line 120 of file MahiFit.h.

Member Typedef Documentation

typedef BXVector::Index MahiFit::Index

Definition at line 146 of file MahiFit.h.

Constructor & Destructor Documentation

MahiFit::MahiFit ( )

Definition at line 4 of file MahiFit.cc.

4  :
5  fullTSSize_(19),
7 {}
const unsigned int fullTSofInterest_
Definition: MahiFit.h:173
const unsigned int fullTSSize_
Definition: MahiFit.h:172
MahiFit::~MahiFit ( )
inline

Definition at line 124 of file MahiFit.h.

References vertices_cff::chi2, and HBHEMahiParameters_cfi::chiSqSwitch.

124 { };

Member Function Documentation

double MahiFit::calculateArrivalTime ( ) const
private

Definition at line 330 of file MahiFit.cc.

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, MahiNnlsWorkspace::fullTSOffset, nnlsWork_, MahiNnlsWorkspace::nPulseTot, PFRecoTauDiscriminationByIsolation_cfi::offset, pedestalBX_, MahiNnlsWorkspace::pulseDerivArray, MahiNnlsWorkspace::pulseDerivMat, MahiNnlsWorkspace::residuals, lumiQTWidget::t, timeLimit_, and MahiNnlsWorkspace::tsSize.

Referenced by doFit().

330  {
331 
333 
334  int itIndex=0;
335 
336  for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
337  int offset=nnlsWork_.bxs.coeff(iBX);
338  if (offset==0) itIndex=iBX;
339 
340  if (offset==pedestalBX_) {
341  nnlsWork_.pulseDerivMat.col(iBX) = SampleVector::Zero(nnlsWork_.tsSize);
342  }
343  else {
345  }
346  }
347 
348  PulseVector solution = nnlsWork_.pulseDerivMat.colPivHouseholderQr().solve(nnlsWork_.residuals);
349  float t = solution.coeff(itIndex)/nnlsWork_.ampVec.coeff(itIndex);
350  t = (t>timeLimit_) ? timeLimit_ :
351  ((t<-timeLimit_) ? -timeLimit_ : t);
352 
353  return t;
354 
355 }
unsigned int nPulseTot
Definition: MahiFit.h:19
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
PulseVector ampVec
Definition: MahiFit.h:67
SamplePulseMatrix pulseDerivMat
Definition: MahiFit.h:60
PulseVector residuals
Definition: MahiFit.h:63
static float timeLimit_
Definition: MahiFit.h:179
std::array< FullSampleVector, MaxPVSize > pulseDerivArray
Definition: MahiFit.h:49
unsigned int fullTSOffset
Definition: MahiFit.h:22
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, PulseVectorSize, 1 > PulseVector
unsigned int tsSize
Definition: MahiFit.h:20
BXVector bxs
Definition: MahiFit.h:27
static int pedestalBX_
Definition: MahiFit.h:175
double MahiFit::calculateChiSq ( ) const
private

Definition at line 470 of file MahiFit.cc.

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

Referenced by minimize().

470  {
471 
472  return (nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.pulseMat*nnlsWork_.ampVec - nnlsWork_.amplitudes)).squaredNorm();
473 }
SampleVector amplitudes
Definition: MahiFit.h:30
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
PulseVector ampVec
Definition: MahiFit.h:67
SampleDecompLLT covDecomp
Definition: MahiFit.h:77
SamplePulseMatrix pulseMat
Definition: MahiFit.h:57
void MahiFit::doFit ( std::array< float, 3 > &  correctedOutput,
const int  nbx 
) const

Definition at line 128 of file MahiFit.cc.

References activeBXs_, MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::bxOffset, bxOffsetConf_, MahiNnlsWorkspace::bxs, bxSizeConf_, calculateArrivalTime(), dynamicPed_, MahiNnlsWorkspace::errVec, MahiNnlsWorkspace::fullTSOffset, MaxFSVSize, minimize(), nnlsWork_, MahiNnlsWorkspace::nPulseTot, PFRecoTauDiscriminationByIsolation_cfi::offset, pedestalBX_, MahiNnlsWorkspace::pulseCovArray, MahiNnlsWorkspace::pulseDerivArray, MahiNnlsWorkspace::pulseMat, MahiNnlsWorkspace::pulseShapeArray, MahiNnlsWorkspace::residuals, BXVector< T >::resize(), MahiNnlsWorkspace::tsOffset, MahiNnlsWorkspace::tsSize, and updatePulseShape().

Referenced by phase1Apply().

128  {
129 
130  unsigned int bxSize=1;
131 
132  if (nbx==1) {
133  nnlsWork_.bxOffset = 0;
134  }
135  else {
136  bxSize = bxSizeConf_;
138  }
139 
140  nnlsWork_.bxs.resize(bxSize);
141 
142  if (nbx==1) {
143  nnlsWork_.bxs.coeffRef(0) = 0;
144  }
145  else {
146  for (unsigned int iBX=0; iBX<bxSize; ++iBX) {
147  nnlsWork_.bxs.coeffRef(iBX) = activeBXs_[iBX];
148  }
149  }
150 
151  nnlsWork_.nPulseTot = bxSize;
152 
153  if (dynamicPed_) {
157  }
158 
160  nnlsWork_.ampVec = PulseVector::Zero(nnlsWork_.nPulseTot);
161  nnlsWork_.errVec = PulseVector::Zero(nnlsWork_.nPulseTot);
162 
163  int offset=0;
164  for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
165  offset=nnlsWork_.bxs.coeff(iBX);
166 
167  nnlsWork_.pulseShapeArray[iBX] = FullSampleVector::Zero(MaxFSVSize);
168  nnlsWork_.pulseDerivArray[iBX] = FullSampleVector::Zero(MaxFSVSize);
170 
171  if (offset==pedestalBX_) {
172  nnlsWork_.ampVec.coeffRef(iBX) = 0;
173  }
174  else {
175 
179  nnlsWork_.pulseCovArray[iBX]);
180 
181  nnlsWork_.ampVec.coeffRef(iBX)=0;
182 
184 
185  }
186  }
187 
188  nnlsWork_.pulseMat.col(nnlsWork_.nPulseTot-1) = SampleVector::Ones(nnlsWork_.tsSize);
189 
192 
193  double chiSq = minimize();
194 
196 
197  bool foundintime = false;
198  unsigned int ipulseintime = 0;
199 
200  for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
201  if (nnlsWork_.bxs.coeff(iBX)==0) {
202  ipulseintime = iBX;
203  foundintime = true;
204  }
205  }
206 
207  if (foundintime) {
208  correctedOutput.at(0) = nnlsWork_.ampVec.coeff(ipulseintime); //charge
209  if (correctedOutput.at(0)!=0) {
210  double arrivalTime = calculateArrivalTime();
211  correctedOutput.at(1) = arrivalTime; //time
212  }
213  else correctedOutput.at(1) = -9999;//time
214 
215  correctedOutput.at(2) = chiSq; //chi2
216 
217  }
218 
219 }
unsigned int nPulseTot
Definition: MahiFit.h:19
SampleVector amplitudes
Definition: MahiFit.h:30
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
PulseVector ampVec
Definition: MahiFit.h:67
constexpr int MaxFSVSize
std::array< FullSampleVector, MaxPVSize > pulseShapeArray
Definition: MahiFit.h:46
int bxOffsetConf_
Definition: MahiFit.h:203
PulseVector aTbVec
Definition: MahiFit.h:74
unsigned int bxSizeConf_
Definition: MahiFit.h:202
PulseVector residuals
Definition: MahiFit.h:63
double calculateArrivalTime() const
Definition: MahiFit.cc:330
std::vector< int > activeBXs_
Definition: MahiFit.h:194
std::array< FullSampleVector, MaxPVSize > pulseDerivArray
Definition: MahiFit.h:49
Polynomial< 0 > Constant
Definition: Constant.h:6
void updatePulseShape(double itQ, FullSampleVector &pulseShape, FullSampleVector &pulseDeriv, FullSampleMatrix &pulseCov) const
Definition: MahiFit.cc:254
unsigned int tsOffset
Definition: MahiFit.h:21
void resize(int bx, unsigned size)
PulseMatrix aTaMat
Definition: MahiFit.h:73
unsigned int fullTSOffset
Definition: MahiFit.h:22
SamplePulseMatrix pulseMat
Definition: MahiFit.h:57
bool dynamicPed_
Definition: MahiFit.h:182
std::array< FullSampleMatrix, MaxPVSize > pulseCovArray
Definition: MahiFit.h:43
unsigned int tsSize
Definition: MahiFit.h:20
BXVector bxs
Definition: MahiFit.h:27
PulseVector errVec
Definition: MahiFit.h:69
double minimize() const
Definition: MahiFit.cc:221
static int pedestalBX_
Definition: MahiFit.h:175
double MahiFit::minimize ( ) const
private

Definition at line 221 of file MahiFit.cc.

References funct::abs(), calculateChiSq(), deltaChiSqThresh_, nMaxItersMin_, nnls(), nnlsWork_, MahiNnlsWorkspace::nPulseTot, onePulseMinimize(), and updateCov().

Referenced by doFit().

221  {
222 
223  double oldChiSq=9999;
224  double chiSq=oldChiSq;
225 
226  for( int iter=1; iter<nMaxItersMin_ ; ++iter) {
227 
228  updateCov();
229 
230  if (nnlsWork_.nPulseTot>1) {
231  nnls();
232  }
233  else {
235  }
236 
237  double newChiSq=calculateChiSq();
238  double deltaChiSq = newChiSq - chiSq;
239 
240  if (newChiSq==oldChiSq && newChiSq<chiSq) {
241  break;
242  }
243  oldChiSq=chiSq;
244  chiSq = newChiSq;
245 
246  if (std::abs(deltaChiSq)<deltaChiSqThresh_) break;
247 
248  }
249 
250  return chiSq;
251 
252 }
unsigned int nPulseTot
Definition: MahiFit.h:19
void nnls() const
Definition: MahiFit.cc:358
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
void onePulseMinimize() const
Definition: MahiFit.cc:458
int nMaxItersMin_
Definition: MahiFit.h:196
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double calculateChiSq() const
Definition: MahiFit.cc:470
float deltaChiSqThresh_
Definition: MahiFit.h:199
void updateCov() const
Definition: MahiFit.cc:308
void MahiFit::nnls ( ) const
private

Definition at line 358 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::ampvecpermtest, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, MahiNnlsWorkspace::covDecomp, MahiNnlsWorkspace::fullTSOffset, mps_fire::i, MahiNnlsWorkspace::invcovp, SiStripPI::max, min(), nMaxItersNNLS_, nnlsConstrainParameter(), nnlsThresh_, nnlsUnconstrainParameter(), nnlsWork_, MahiNnlsWorkspace::nP, MahiNnlsWorkspace::nPulseTot, PFRecoTauDiscriminationByIsolation_cfi::offset, pedestalBX_, MahiNnlsWorkspace::pulseMat, MahiNnlsWorkspace::pulseShapeArray, particleFlowDisplacedVertex_cfi::ratio, solveSubmatrix(), electronIdCutBased_cfi::threshold, MahiNnlsWorkspace::tsSize, and MahiNnlsWorkspace::updateWork.

Referenced by minimize().

358  {
359  const unsigned int npulse = nnlsWork_.nPulseTot;
360 
361  for (unsigned int iBX=0; iBX<npulse; ++iBX) {
362  int offset=nnlsWork_.bxs.coeff(iBX);
363  if (offset==pedestalBX_) {
364  nnlsWork_.pulseMat.col(iBX) = SampleVector::Ones(nnlsWork_.tsSize);
365  }
366  else {
368  }
369  }
370 
372  nnlsWork_.aTaMat = nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.invcovp);
373  nnlsWork_.aTbVec = nnlsWork_.invcovp.transpose().lazyProduct(nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes));
374 
375  int iter = 0;
376  Index idxwmax = 0;
377  double wmax = 0.0;
378  double threshold = nnlsThresh_;
379 
380  nnlsWork_.nP=0;
381 
382  while (true) {
383  if (iter>0 || nnlsWork_.nP==0) {
384  if ( nnlsWork_.nP==std::min(npulse, nnlsWork_.tsSize)) break;
385 
386  const unsigned int nActive = npulse - nnlsWork_.nP;
388 
389  Index idxwmaxprev = idxwmax;
390  double wmaxprev = wmax;
391  wmax = nnlsWork_.updateWork.tail(nActive).maxCoeff(&idxwmax);
392 
393  if (wmax<threshold || (idxwmax==idxwmaxprev && wmax==wmaxprev)) {
394  break;
395  }
396 
397  if (iter>=nMaxItersNNLS_) {
398  break;
399  }
400 
401  //unconstrain parameter
402  Index idxp = nnlsWork_.nP + idxwmax;
404 
405  }
406 
407  while (true) {
408  if (nnlsWork_.nP==0) break;
409 
411 
413 
414  //check solution
415  bool positive = true;
416  for (unsigned int i = 0; i < nnlsWork_.nP; ++i)
417  positive &= (nnlsWork_.ampvecpermtest(i) > 0);
418  if (positive) {
420  break;
421  }
422 
423  //update parameter vector
424  Index minratioidx=0;
425 
426  // no realizable optimization here (because it autovectorizes!)
427  double minratio = std::numeric_limits<double>::max();
428  for (unsigned int ipulse=0; ipulse<nnlsWork_.nP; ++ipulse) {
429  if (nnlsWork_.ampvecpermtest.coeff(ipulse)<=0.) {
430  const double c_ampvec = nnlsWork_.ampVec.coeff(ipulse);
431  const double ratio = c_ampvec/(c_ampvec-nnlsWork_.ampvecpermtest.coeff(ipulse));
432  if (ratio<minratio) {
433  minratio = ratio;
434  minratioidx = ipulse;
435  }
436  }
437  }
439 
440  //avoid numerical problems with later ==0. check
441  nnlsWork_.ampVec.coeffRef(minratioidx) = 0.;
442 
443  nnlsConstrainParameter(minratioidx);
444  }
445 
446  ++iter;
447 
448  //adaptive convergence threshold to avoid infinite loops but still
449  //ensure best value is used
450  if (iter%10==0) {
451  threshold *= 10.;
452  }
453 
454  }
455 
456 }
unsigned int nPulseTot
Definition: MahiFit.h:19
SamplePulseMatrix invcovp
Definition: MahiFit.h:72
void nnlsConstrainParameter(Index minratioidx) const
Definition: MahiFit.cc:510
SampleVector amplitudes
Definition: MahiFit.h:30
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
BXVector::Index Index
Definition: MahiFit.h:146
void solveSubmatrix(PulseMatrix &mat, PulseVector &invec, PulseVector &outvec, unsigned nP) const
Definition: MahiFit.cc:594
PulseVector ampVec
Definition: MahiFit.h:67
SampleDecompLLT covDecomp
Definition: MahiFit.h:77
PulseVector ampvecpermtest
Definition: MahiFit.h:70
std::array< FullSampleVector, MaxPVSize > pulseShapeArray
Definition: MahiFit.h:46
void nnlsUnconstrainParameter(Index idxp) const
Definition: MahiFit.cc:500
PulseVector aTbVec
Definition: MahiFit.h:74
float nnlsThresh_
Definition: MahiFit.h:200
T min(T a, T b)
Definition: MathUtil.h:58
PulseMatrix aTaMat
Definition: MahiFit.h:73
unsigned int fullTSOffset
Definition: MahiFit.h:22
SamplePulseMatrix pulseMat
Definition: MahiFit.h:57
unsigned int tsSize
Definition: MahiFit.h:20
unsigned int nP
Definition: MahiFit.h:66
BXVector bxs
Definition: MahiFit.h:27
int nMaxItersNNLS_
Definition: MahiFit.h:197
PulseVector updateWork
Definition: MahiFit.h:75
static int pedestalBX_
Definition: MahiFit.h:175
void MahiFit::nnlsConstrainParameter ( Index  minratioidx) const
private

Definition at line 510 of file MahiFit.cc.

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

Referenced by nnls().

510  {
511  nnlsWork_.aTaMat.col(nnlsWork_.nP-1).swap(nnlsWork_.aTaMat.col(minratioidx));
512  nnlsWork_.aTaMat.row(nnlsWork_.nP-1).swap(nnlsWork_.aTaMat.row(minratioidx));
513  nnlsWork_.pulseMat.col(nnlsWork_.nP-1).swap(nnlsWork_.pulseMat.col(minratioidx));
514  std::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP-1),nnlsWork_.aTbVec.coeffRef(minratioidx));
515  std::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP-1),nnlsWork_.ampVec.coeffRef(minratioidx));
516  std::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP-1),nnlsWork_.bxs.coeffRef(minratioidx));
517  --nnlsWork_.nP;
518 
519 }
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
PulseVector ampVec
Definition: MahiFit.h:67
PulseVector aTbVec
Definition: MahiFit.h:74
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
PulseMatrix aTaMat
Definition: MahiFit.h:73
SamplePulseMatrix pulseMat
Definition: MahiFit.h:57
unsigned int nP
Definition: MahiFit.h:66
BXVector bxs
Definition: MahiFit.h:27
void MahiFit::nnlsUnconstrainParameter ( Index  idxp) const
private

Definition at line 500 of file MahiFit.cc.

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

Referenced by nnls().

500  {
501  nnlsWork_.aTaMat.col(nnlsWork_.nP).swap(nnlsWork_.aTaMat.col(idxp));
502  nnlsWork_.aTaMat.row(nnlsWork_.nP).swap(nnlsWork_.aTaMat.row(idxp));
503  nnlsWork_.pulseMat.col(nnlsWork_.nP).swap(nnlsWork_.pulseMat.col(idxp));
504  std::swap(nnlsWork_.aTbVec.coeffRef(nnlsWork_.nP),nnlsWork_.aTbVec.coeffRef(idxp));
505  std::swap(nnlsWork_.ampVec.coeffRef(nnlsWork_.nP),nnlsWork_.ampVec.coeffRef(idxp));
506  std::swap(nnlsWork_.bxs.coeffRef(nnlsWork_.nP),nnlsWork_.bxs.coeffRef(idxp));
507  ++nnlsWork_.nP;
508 }
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
PulseVector ampVec
Definition: MahiFit.h:67
PulseVector aTbVec
Definition: MahiFit.h:74
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
PulseMatrix aTaMat
Definition: MahiFit.h:73
SamplePulseMatrix pulseMat
Definition: MahiFit.h:57
unsigned int nP
Definition: MahiFit.h:66
BXVector bxs
Definition: MahiFit.h:27
void MahiFit::onePulseMinimize ( ) const
private

Definition at line 458 of file MahiFit.cc.

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

Referenced by minimize().

458  {
459 
461 
462  SingleMatrix aTamatval = nnlsWork_.invcovp.transpose()*nnlsWork_.invcovp;
463  SingleVector aTbvecval = nnlsWork_.invcovp.transpose()*nnlsWork_.covDecomp.matrixL().solve(nnlsWork_.amplitudes);
464 
465  nnlsWork_.ampVec.coeffRef(0) = std::max(0., aTbvecval.coeff(0)/aTamatval.coeff(0));
466 
467 
468 }
SamplePulseMatrix invcovp
Definition: MahiFit.h:72
SampleVector amplitudes
Definition: MahiFit.h:30
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
PulseVector ampVec
Definition: MahiFit.h:67
SampleDecompLLT covDecomp
Definition: MahiFit.h:77
Eigen::Matrix< double, 1, 1 > SingleVector
SamplePulseMatrix pulseMat
Definition: MahiFit.h:57
Eigen::Matrix< double, 1, 1 > SingleMatrix
void MahiFit::phase1Apply ( const HBHEChannelInfo channelData,
float &  reconstructedEnergy,
float &  reconstructedTime,
bool &  useTriple,
float &  chi2 
) const

Definition at line 40 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, ALCARECOTkAlJpsiMuMu_cff::charge, chiSqSwitch_, doFit(), MahiNnlsWorkspace::dt, HBHEChannelInfo::fcByPE(), MahiNnlsWorkspace::fullTSOffset, fullTSofInterest_, HBHEChannelInfo::hasTimeInfo(), nnlsWork_, MahiNnlsWorkspace::noiseTerms, HBHEChannelInfo::nSamples(), MahiNnlsWorkspace::pedConstraint, resetWorkspace(), HBHEChannelInfo::soi(), mathSSE::sqrt(), timeSigmaHPD_, timeSigmaSiPM_, ts4Thresh_, HBHEChannelInfo::tsDFcPerADC(), HBHEChannelInfo::tsGain(), MahiNnlsWorkspace::tsOffset, HBHEChannelInfo::tsPedestal(), HBHEChannelInfo::tsPedestalWidth(), HBHEChannelInfo::tsRawCharge(), and MahiNnlsWorkspace::tsSize.

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

44  {
45 
46  assert(channelData.nSamples()==8||channelData.nSamples()==10);
47 
49 
50  nnlsWork_.tsSize = channelData.nSamples();
51  nnlsWork_.tsOffset = channelData.soi();
53 
54  // 1 sigma time constraint
55  if (channelData.hasTimeInfo()) nnlsWork_.dt=timeSigmaSiPM_;
57 
58 
59  //Average pedestal width (for covariance matrix constraint)
60  float pedVal = 0.25*( channelData.tsPedestalWidth(0)*channelData.tsPedestalWidth(0)+
61  channelData.tsPedestalWidth(1)*channelData.tsPedestalWidth(1)+
62  channelData.tsPedestalWidth(2)*channelData.tsPedestalWidth(2)+
63  channelData.tsPedestalWidth(3)*channelData.tsPedestalWidth(3) );
64 
65  nnlsWork_.pedConstraint = pedVal*SampleMatrix::Ones(nnlsWork_.tsSize, nnlsWork_.tsSize);
68 
69  std::array<float,3> reconstructedVals {{ 0.0, -9999, -9999 }};
70 
71  double tsTOT = 0, tstrig = 0; // in GeV
72  for(unsigned int iTS=0; iTS<nnlsWork_.tsSize; ++iTS){
73  double charge = channelData.tsRawCharge(iTS);
74  double ped = channelData.tsPedestal(iTS);
75 
76  nnlsWork_.amplitudes.coeffRef(iTS) = charge - ped;
77 
78  //ADC granularity
79  double noiseADC = (1./sqrt(12))*channelData.tsDFcPerADC(iTS);
80 
81  //Photostatistics
82  double noisePhoto = 0;
83  if ( (charge-ped)>channelData.tsPedestalWidth(iTS)) {
84  noisePhoto = sqrt((charge-ped)*channelData.fcByPE());
85  }
86 
87  //Electronic pedestal
88  double pedWidth = channelData.tsPedestalWidth(iTS);
89 
90  //Total uncertainty from all sources
91  nnlsWork_.noiseTerms.coeffRef(iTS) = noiseADC*noiseADC + noisePhoto*noisePhoto + pedWidth*pedWidth;
92 
93  tsTOT += (charge - ped)*channelData.tsGain(0);
94  if( iTS==nnlsWork_.tsOffset ){
95  tstrig += (charge - ped)*channelData.tsGain(0);
96  }
97  }
98 
99  if(tstrig >= ts4Thresh_ && tsTOT > 0) {
100 
101  useTriple=false;
102 
103  // only do pre-fit with 1 pulse if chiSq threshold is positive
104  if (chiSqSwitch_>0) {
105  doFit(reconstructedVals,1);
106  if (reconstructedVals[2]>chiSqSwitch_) {
107  doFit(reconstructedVals,0); //nbx=0 means use configured BXs
108  useTriple=true;
109  }
110  }
111  else {
112  doFit(reconstructedVals,0);
113  useTriple=true;
114  }
115  }
116  else{
117  reconstructedVals.at(0) = 0.; //energy
118  reconstructedVals.at(1) = -9999.; //time
119  reconstructedVals.at(2) = -9999.; //chi2
120  }
121 
122  reconstructedEnergy = reconstructedVals[0]*channelData.tsGain(0);
123  reconstructedTime = reconstructedVals[1];
124  chi2 = reconstructedVals[2];
125 
126 }
float chiSqSwitch_
Definition: MahiFit.h:184
bool hasTimeInfo() const
double tsGain(const unsigned ts) const
SampleVector amplitudes
Definition: MahiFit.h:30
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
unsigned soi() const
double tsPedestal(const unsigned ts) const
double tsRawCharge(const unsigned ts) const
T sqrt(T t)
Definition: SSEVec.h:18
void doFit(std::array< float, 3 > &correctedOutput, const int nbx) const
Definition: MahiFit.cc:128
double fcByPE() const
SampleMatrix pedConstraint
Definition: MahiFit.h:39
const unsigned int fullTSofInterest_
Definition: MahiFit.h:173
unsigned int tsOffset
Definition: MahiFit.h:21
void resetWorkspace() const
Definition: MahiFit.cc:663
double tsPedestalWidth(const unsigned ts) const
unsigned int fullTSOffset
Definition: MahiFit.h:22
float timeSigmaSiPM_
Definition: MahiFit.h:192
SampleVector noiseTerms
Definition: MahiFit.h:36
float ts4Thresh_
Definition: MahiFit.h:183
unsigned int tsSize
Definition: MahiFit.h:20
float timeSigmaHPD_
Definition: MahiFit.h:191
unsigned nSamples() const
float tsDFcPerADC(const unsigned ts) const
void MahiFit::phase1Debug ( const HBHEChannelInfo channelData,
MahiDebugInfo mdi 
) const

Definition at line 521 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiDebugInfo::arrivalTime, MahiNnlsWorkspace::bxs, ALCARECOTkAlJpsiMuMu_cff::charge, vertices_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, MahiDebugInfo::nEnergy, nnlsWork_, MahiNnlsWorkspace::noiseTerms, MahiDebugInfo::nPulse, MahiNnlsWorkspace::nPulseTot, MahiDebugInfo::nSamples, HBHEChannelInfo::nSamples(), MahiDebugInfo::pedEnergy, pedestalBX_, MahiDebugInfo::pEnergy, phase1Apply(), MahiDebugInfo::pPulse, 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.

522  {
523 
524  float recoEnergy, recoTime, chi2;
525  bool use3;
526  phase1Apply(channelData, recoEnergy, recoTime, use3, chi2);
527 
528 
529  mdi.nSamples = channelData.nSamples();
530  mdi.soi = channelData.soi();
531 
532  mdi.use3 = use3;
533 
534  mdi.inTimeConst = nnlsWork_.dt;
535  mdi.inPedAvg = 0.25*( channelData.tsPedestalWidth(0)*channelData.tsPedestalWidth(0)+
536  channelData.tsPedestalWidth(1)*channelData.tsPedestalWidth(1)+
537  channelData.tsPedestalWidth(2)*channelData.tsPedestalWidth(2)+
538  channelData.tsPedestalWidth(3)*channelData.tsPedestalWidth(3) );
539  mdi.inGain = channelData.tsGain(0);
540 
541  for (unsigned int iTS=0; iTS<channelData.nSamples(); ++iTS) {
542 
543  double charge = channelData.tsRawCharge(iTS);
544  double ped = channelData.tsPedestal(iTS);
545 
546  mdi.inNoiseADC[iTS] = (1./sqrt(12))*channelData.tsDFcPerADC(iTS);
547 
548  if ( (charge-ped)>channelData.tsPedestalWidth(iTS)) {
549  mdi.inNoisePhoto[iTS] = sqrt((charge-ped)*channelData.fcByPE());
550  }
551  else { mdi.inNoisePhoto[iTS] = 0; }
552 
553  mdi.inPedestal[iTS] = channelData.tsPedestalWidth(iTS);
554  mdi.totalUCNoise[iTS] = nnlsWork_.noiseTerms.coeffRef(iTS);
555 
556  if (channelData.hasTimeInfo()) {
557  mdi.inputTDC[iTS] = channelData.tsRiseTime(iTS);
558  }
559  else { mdi.inputTDC[iTS]=-1; }
560 
561  }
562 
563  mdi.arrivalTime = recoTime;
564  mdi.chiSq = chi2;
565 
566  for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
567  if (nnlsWork_.bxs.coeff(iBX)==0) {
568  mdi.mahiEnergy=nnlsWork_.ampVec.coeff(iBX);
569  for(unsigned int iTS=0; iTS<nnlsWork_.tsSize; ++iTS){
570  mdi.count[iTS] = iTS;
571  mdi.inputTS[iTS] = nnlsWork_.amplitudes.coeff(iTS);
572  mdi.itPulse[iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
573  }
574  }
575  else if (nnlsWork_.bxs.coeff(iBX)==pedestalBX_) {
576  mdi.pedEnergy=nnlsWork_.ampVec.coeff(iBX);
577  }
578  else if (nnlsWork_.bxs.coeff(iBX)==-1) {
579  mdi.pEnergy=nnlsWork_.ampVec.coeff(iBX);
580  for(unsigned int iTS=0; iTS<nnlsWork_.tsSize; ++iTS){
581  mdi.pPulse[iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
582  }
583  }
584  else if (nnlsWork_.bxs.coeff(iBX)==1) {
585  mdi.nEnergy=nnlsWork_.ampVec.coeff(iBX);
586  for(unsigned int iTS=0; iTS<nnlsWork_.tsSize; ++iTS){
587  mdi.nPulse[iTS] = nnlsWork_.pulseMat.col(iBX).coeff(iTS);
588  }
589  }
590  }
591 }
unsigned int nPulseTot
Definition: MahiFit.h:19
void phase1Apply(const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &reconstructedTime, bool &useTriple, float &chi2) const
Definition: MahiFit.cc:40
float itPulse[MaxSVSize]
Definition: MahiFit.h:114
bool hasTimeInfo() const
double tsGain(const unsigned ts) const
SampleVector amplitudes
Definition: MahiFit.h:30
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
unsigned soi() const
PulseVector ampVec
Definition: MahiFit.h:67
float count[MaxSVSize]
Definition: MahiFit.h:111
double tsPedestal(const unsigned ts) const
float totalUCNoise[MaxSVSize]
Definition: MahiFit.h:101
float nPulse[MaxSVSize]
Definition: MahiFit.h:116
float inNoiseADC[MaxSVSize]
Definition: MahiFit.h:96
float pedEnergy
Definition: MahiFit.h:109
double tsRawCharge(const unsigned ts) const
float inPedestal[MaxSVSize]
Definition: MahiFit.h:99
T sqrt(T t)
Definition: SSEVec.h:18
double fcByPE() const
int nSamples
Definition: MahiFit.h:86
float chiSq
Definition: MahiFit.h:104
float nEnergy
Definition: MahiFit.h:108
bool use3
Definition: MahiFit.h:89
float inGain
Definition: MahiFit.h:94
float arrivalTime
Definition: MahiFit.h:105
float inPedAvg
Definition: MahiFit.h:93
float tsRiseTime(const unsigned ts) const
double tsPedestalWidth(const unsigned ts) const
float mahiEnergy
Definition: MahiFit.h:103
SamplePulseMatrix pulseMat
Definition: MahiFit.h:57
float inputTS[MaxSVSize]
Definition: MahiFit.h:112
int inputTDC[MaxSVSize]
Definition: MahiFit.h:113
float inNoisePhoto[MaxSVSize]
Definition: MahiFit.h:98
SampleVector noiseTerms
Definition: MahiFit.h:36
float pEnergy
Definition: MahiFit.h:107
float pPulse[MaxSVSize]
Definition: MahiFit.h:115
unsigned int tsSize
Definition: MahiFit.h:20
BXVector bxs
Definition: MahiFit.h:27
float inTimeConst
Definition: MahiFit.h:91
unsigned nSamples() const
static int pedestalBX_
Definition: MahiFit.h:175
float tsDFcPerADC(const unsigned ts) const
void MahiFit::resetPulseShapeTemplate ( const HcalPulseShapes::Shape ps)

Definition at line 488 of file MahiFit.cc.

References cntsetPulseShape_, pfunctor_, psfPtr_, and FitterFuncs::PulseShapeFunctor::singlePulseShapeFunc().

Referenced by setPulseShapeTemplate().

488  {
490 
491  // only the pulse shape itself from PulseShapeFunctor is used for Mahi
492  // the uncertainty terms calculated inside PulseShapeFunctor are used for Method 2 only
493  psfPtr_.reset(new FitterFuncs::PulseShapeFunctor(ps,false,false,false,
494  1,0,0,10));
495  pfunctor_ = std::unique_ptr<ROOT::Math::Functor>( new ROOT::Math::Functor(psfPtr_.get(),&FitterFuncs::PulseShapeFunctor::singlePulseShapeFunc, 3) );
496 
497 
498 }
double singlePulseShapeFunc(const double *x)
int cntsetPulseShape_
Definition: MahiFit.h:206
std::unique_ptr< FitterFuncs::PulseShapeFunctor > psfPtr_
Definition: MahiFit.h:207
std::unique_ptr< ROOT::Math::Functor > pfunctor_
Definition: MahiFit.h:208
void MahiFit::resetWorkspace ( ) const
private

Definition at line 663 of file MahiFit.cc.

References MahiNnlsWorkspace::amplitudes, MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::ampvecpermtest, MahiNnlsWorkspace::aTaMat, MahiNnlsWorkspace::aTbVec, begin, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, MahiNnlsWorkspace::covDecompLinv, MahiNnlsWorkspace::dt, end, MahiNnlsWorkspace::errVec, lumiContext::fill, MahiNnlsWorkspace::fullTSOffset, MahiNnlsWorkspace::invCovMat, MahiNnlsWorkspace::invcovp, nnlsWork_, MahiNnlsWorkspace::noiseTerms, MahiNnlsWorkspace::nPulseTot, MahiNnlsWorkspace::pedConstraint, MahiNnlsWorkspace::pulseCovArray, MahiNnlsWorkspace::pulseDerivArray, MahiNnlsWorkspace::pulseDerivMat, MahiNnlsWorkspace::pulseM, MahiNnlsWorkspace::pulseMat, MahiNnlsWorkspace::pulseN, MahiNnlsWorkspace::pulseP, MahiNnlsWorkspace::pulseShapeArray, MahiNnlsWorkspace::residuals, MahiNnlsWorkspace::topleft_work, MahiNnlsWorkspace::tsOffset, MahiNnlsWorkspace::tsSize, and MahiNnlsWorkspace::updateWork.

Referenced by phase1Apply().

663  {
664 
666  nnlsWork_.tsSize=0;
670  nnlsWork_.dt=0;
671 
675 
679 
680  nnlsWork_.amplitudes.setZero();
681  nnlsWork_.bxs.setZero();
682  nnlsWork_.invCovMat.setZero();
683  nnlsWork_.noiseTerms.setZero();
684  nnlsWork_.pedConstraint.setZero();
685  nnlsWork_.pulseMat.setZero();
686  nnlsWork_.pulseDerivMat.setZero();
687  nnlsWork_.residuals.setZero();
688  nnlsWork_.ampVec.setZero();
689  nnlsWork_.errVec.setZero();
690  nnlsWork_.ampvecpermtest.setZero();
691  nnlsWork_.invcovp.setZero();
692  nnlsWork_.aTaMat.setZero();
693  nnlsWork_.aTbVec.setZero();
694  nnlsWork_.updateWork.setZero();
695  nnlsWork_.covDecompLinv.setZero();
696  nnlsWork_.topleft_work.setZero();
697 
698 
699 
700 }
unsigned int nPulseTot
Definition: MahiFit.h:19
SamplePulseMatrix invcovp
Definition: MahiFit.h:72
SampleVector amplitudes
Definition: MahiFit.h:30
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
PulseVector ampVec
Definition: MahiFit.h:67
PulseVector ampvecpermtest
Definition: MahiFit.h:70
std::array< FullSampleVector, MaxPVSize > pulseShapeArray
Definition: MahiFit.h:46
std::array< double, MaxSVSize > pulseM
Definition: MahiFit.h:53
PulseMatrix topleft_work
Definition: MahiFit.h:79
SamplePulseMatrix pulseDerivMat
Definition: MahiFit.h:60
std::array< double, MaxSVSize > pulseN
Definition: MahiFit.h:52
PulseVector aTbVec
Definition: MahiFit.h:74
PulseVector residuals
Definition: MahiFit.h:63
#define end
Definition: vmac.h:39
SampleMatrix pedConstraint
Definition: MahiFit.h:39
std::array< FullSampleVector, MaxPVSize > pulseDerivArray
Definition: MahiFit.h:49
SampleMatrix covDecompLinv
Definition: MahiFit.h:78
SampleMatrix invCovMat
Definition: MahiFit.h:33
unsigned int tsOffset
Definition: MahiFit.h:21
PulseMatrix aTaMat
Definition: MahiFit.h:73
unsigned int fullTSOffset
Definition: MahiFit.h:22
SamplePulseMatrix pulseMat
Definition: MahiFit.h:57
std::array< FullSampleMatrix, MaxPVSize > pulseCovArray
Definition: MahiFit.h:43
SampleVector noiseTerms
Definition: MahiFit.h:36
#define begin
Definition: vmac.h:32
unsigned int tsSize
Definition: MahiFit.h:20
BXVector bxs
Definition: MahiFit.h:27
PulseVector errVec
Definition: MahiFit.h:69
std::array< double, MaxSVSize > pulseP
Definition: MahiFit.h:54
PulseVector updateWork
Definition: MahiFit.h:75
void MahiFit::setParameters ( bool  iDynamicPed,
double  iTS4Thresh,
double  chiSqSwitch,
bool  iApplyTimeSlew,
HcalTimeSlew::BiasSetting  slewFlavor,
double  iMeanTime,
double  iTimeSigmaHPD,
double  iTimeSigmaSiPM,
const std::vector< int > &  iActiveBXs,
int  iNMaxItersMin,
int  iNMaxItersNNLS,
double  iDeltaChiSqThresh,
double  iNnlsThresh 
)

Definition at line 9 of file MahiFit.cc.

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

13  {
14 
15  dynamicPed_ = iDynamicPed;
16 
17  ts4Thresh_ = iTS4Thresh;
19 
20  applyTimeSlew_ = iApplyTimeSlew;
21  slewFlavor_ = slewFlavor;
22 
23  meanTime_ = iMeanTime;
24  timeSigmaHPD_ = iTimeSigmaHPD;
25  timeSigmaSiPM_ = iTimeSigmaSiPM;
26 
27  activeBXs_ = iActiveBXs;
28 
29  nMaxItersMin_ = iNMaxItersMin;
30  nMaxItersNNLS_ = iNMaxItersNNLS;
31 
32  deltaChiSqThresh_ = iDeltaChiSqThresh;
33  nnlsThresh_ = iNnlsThresh;
34 
35  bxOffsetConf_ = -(*std::min_element(activeBXs_.begin(), activeBXs_.end()));
36  bxSizeConf_ = activeBXs_.size();
37 
38 }
float chiSqSwitch_
Definition: MahiFit.h:184
HcalTimeSlew::BiasSetting slewFlavor_
Definition: MahiFit.h:187
bool applyTimeSlew_
Definition: MahiFit.h:186
float meanTime_
Definition: MahiFit.h:190
int bxOffsetConf_
Definition: MahiFit.h:203
float nnlsThresh_
Definition: MahiFit.h:200
unsigned int bxSizeConf_
Definition: MahiFit.h:202
int nMaxItersMin_
Definition: MahiFit.h:196
std::vector< int > activeBXs_
Definition: MahiFit.h:194
bool dynamicPed_
Definition: MahiFit.h:182
float timeSigmaSiPM_
Definition: MahiFit.h:192
float ts4Thresh_
Definition: MahiFit.h:183
float deltaChiSqThresh_
Definition: MahiFit.h:199
float timeSigmaHPD_
Definition: MahiFit.h:191
int nMaxItersNNLS_
Definition: MahiFit.h:197
void MahiFit::setPulseShapeTemplate ( const HcalPulseShapes::Shape ps,
const HcalTimeSlew hcalTimeSlewDelay 
)

Definition at line 475 of file MahiFit.cc.

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

475  {
476 
477  if (!(&ps == currentPulseShape_ ))
478  {
479 
480  hcalTimeSlewDelay_ = hcalTimeSlewDelay;
481  tsDelay1GeV_= hcalTimeSlewDelay->delay(1.0, slewFlavor_);
482 
484  currentPulseShape_ = &ps;
485  }
486 }
HcalTimeSlew::BiasSetting slewFlavor_
Definition: MahiFit.h:187
double delay(double 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:14
void resetPulseShapeTemplate(const HcalPulseShapes::Shape &ps)
Definition: MahiFit.cc:488
double tsDelay1GeV_
Definition: MahiFit.h:188
const HcalTimeSlew * hcalTimeSlewDelay_
Definition: MahiFit.h:148
const HcalPulseShapes::Shape * currentPulseShape_
Definition: MahiFit.h:147
void MahiFit::solveSubmatrix ( PulseMatrix mat,
PulseVector invec,
PulseVector outvec,
unsigned  nP 
) const
private

Definition at line 594 of file MahiFit.cc.

References Exception, and groupFilesInBlocks::temp.

Referenced by nnls().

594  {
595  using namespace Eigen;
596  switch( nP ) { // pulse matrix is always square.
597  case 10:
598  {
599  Matrix<double,10,10> temp = mat;
600  outvec.head<10>() = temp.ldlt().solve(invec.head<10>());
601  }
602  break;
603  case 9:
604  {
605  Matrix<double,9,9> temp = mat.topLeftCorner<9,9>();
606  outvec.head<9>() = temp.ldlt().solve(invec.head<9>());
607  }
608  break;
609  case 8:
610  {
611  Matrix<double,8,8> temp = mat.topLeftCorner<8,8>();
612  outvec.head<8>() = temp.ldlt().solve(invec.head<8>());
613  }
614  break;
615  case 7:
616  {
617  Matrix<double,7,7> temp = mat.topLeftCorner<7,7>();
618  outvec.head<7>() = temp.ldlt().solve(invec.head<7>());
619  }
620  break;
621  case 6:
622  {
623  Matrix<double,6,6> temp = mat.topLeftCorner<6,6>();
624  outvec.head<6>() = temp.ldlt().solve(invec.head<6>());
625  }
626  break;
627  case 5:
628  {
629  Matrix<double,5,5> temp = mat.topLeftCorner<5,5>();
630  outvec.head<5>() = temp.ldlt().solve(invec.head<5>());
631  }
632  break;
633  case 4:
634  {
635  Matrix<double,4,4> temp = mat.topLeftCorner<4,4>();
636  outvec.head<4>() = temp.ldlt().solve(invec.head<4>());
637  }
638  break;
639  case 3:
640  {
641  Matrix<double,3,3> temp = mat.topLeftCorner<3,3>();
642  outvec.head<3>() = temp.ldlt().solve(invec.head<3>());
643  }
644  break;
645  case 2:
646  {
647  Matrix<double,2,2> temp = mat.topLeftCorner<2,2>();
648  outvec.head<2>() = temp.ldlt().solve(invec.head<2>());
649  }
650  break;
651  case 1:
652  {
653  Matrix<double,1,1> temp = mat.topLeftCorner<1,1>();
654  outvec.head<1>() = temp.ldlt().solve(invec.head<1>());
655  }
656  break;
657  default:
658  throw cms::Exception("HcalMahiWeirdState")
659  << "Weird number of pulses encountered in Mahi, module is configured incorrectly!";
660  }
661 }
void MahiFit::updateCov ( ) const
private

Definition at line 308 of file MahiFit.cc.

References MahiNnlsWorkspace::ampVec, MahiNnlsWorkspace::bxOffset, MahiNnlsWorkspace::bxs, MahiNnlsWorkspace::covDecomp, MahiNnlsWorkspace::fullTSOffset, MahiNnlsWorkspace::invCovMat, nnlsWork_, MahiNnlsWorkspace::noiseTerms, MahiNnlsWorkspace::nPulseTot, PFRecoTauDiscriminationByIsolation_cfi::offset, MahiNnlsWorkspace::pedConstraint, pedestalBX_, MahiNnlsWorkspace::pulseCovArray, and MahiNnlsWorkspace::tsSize.

Referenced by minimize().

308  {
309 
311 
312  nnlsWork_.invCovMat = nnlsWork_.noiseTerms.asDiagonal();
314 
315  for (unsigned int iBX=0; iBX<nnlsWork_.nPulseTot; ++iBX) {
316  if (nnlsWork_.ampVec.coeff(iBX)==0) continue;
317 
318  int offset=nnlsWork_.bxs.coeff(iBX);
319 
320  if (offset==pedestalBX_) continue;
321  else {
322  nnlsWork_.invCovMat += nnlsWork_.ampVec.coeff(iBX)*nnlsWork_.ampVec.coeff(iBX)
324  }
325  }
326 
328 }
unsigned int nPulseTot
Definition: MahiFit.h:19
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
PulseVector ampVec
Definition: MahiFit.h:67
SampleDecompLLT covDecomp
Definition: MahiFit.h:77
SampleMatrix pedConstraint
Definition: MahiFit.h:39
SampleMatrix invCovMat
Definition: MahiFit.h:33
unsigned int fullTSOffset
Definition: MahiFit.h:22
std::array< FullSampleMatrix, MaxPVSize > pulseCovArray
Definition: MahiFit.h:43
SampleVector noiseTerms
Definition: MahiFit.h:36
unsigned int tsSize
Definition: MahiFit.h:20
BXVector bxs
Definition: MahiFit.h:27
static int pedestalBX_
Definition: MahiFit.h:175
void MahiFit::updatePulseShape ( double  itQ,
FullSampleVector pulseShape,
FullSampleVector pulseDeriv,
FullSampleMatrix pulseCov 
) const
private

Definition at line 254 of file MahiFit.cc.

References applyTimeSlew_, HcalTimeSlew::delay(), delta, MahiNnlsWorkspace::dt, MahiNnlsWorkspace::fullTSOffset, hcalTimeSlewDelay_, meanTime_, nnlsWork_, psfPtr_, MahiNnlsWorkspace::pulseM, MahiNnlsWorkspace::pulseN, MahiNnlsWorkspace::pulseP, slewFlavor_, cscNeutronWriter_cfi::t0, tmp, tsDelay1GeV_, MahiNnlsWorkspace::tsSize, and geometryCSVtoXML::xx.

Referenced by doFit().

255  {
256 
257  float t0=meanTime_;
258 
259  if(applyTimeSlew_) {
260  if(itQ<=1.0) t0+=tsDelay1GeV_;
261  else t0+=hcalTimeSlewDelay_->delay(itQ,slewFlavor_);
262  }
263 
264  nnlsWork_.pulseN.fill(0);
265  nnlsWork_.pulseM.fill(0);
266  nnlsWork_.pulseP.fill(0);
267 
268  const double xx[4]={t0, 1.0, 0.0, 3};
269  const double xxm[4]={-nnlsWork_.dt+t0, 1.0, 0.0, 3};
270  const double xxp[4]={ nnlsWork_.dt+t0, 1.0, 0.0, 3};
271 
272  (*pfunctor_)(&xx[0]);
273  psfPtr_->getPulseShape(nnlsWork_.pulseN);
274 
275  (*pfunctor_)(&xxm[0]);
276  psfPtr_->getPulseShape(nnlsWork_.pulseM);
277 
278  (*pfunctor_)(&xxp[0]);
279  psfPtr_->getPulseShape(nnlsWork_.pulseP);
280 
281  //in the 2018+ case where the sample of interest (SOI) is in TS3, add an extra offset to align
282  //with previous SOI=TS4 case assumed by psfPtr_->getPulseShape()
283  int delta =nnlsWork_. tsOffset == 3 ? 1 : 0;
284 
285  for (unsigned int iTS=nnlsWork_.fullTSOffset; iTS<nnlsWork_.fullTSOffset + nnlsWork_.tsSize; ++iTS) {
286 
287  pulseShape.coeffRef(iTS) = nnlsWork_.pulseN[iTS-nnlsWork_.fullTSOffset+delta];
288  pulseDeriv.coeffRef(iTS) = 0.5*(nnlsWork_.pulseM[iTS-nnlsWork_.fullTSOffset+delta]+nnlsWork_.pulseP[iTS-nnlsWork_.fullTSOffset+delta])/(2*nnlsWork_.dt);
289 
292  }
293 
294  for (unsigned int iTS=nnlsWork_.fullTSOffset; iTS<nnlsWork_.fullTSOffset+nnlsWork_.tsSize; ++iTS) {
295  for (unsigned int jTS=nnlsWork_.fullTSOffset; jTS<iTS+1; ++jTS) {
296 
299 
300  pulseCov(iTS,jTS) += tmp;
301  pulseCov(jTS,iTS) += tmp;
302 
303  }
304  }
305 
306 }
dbl * delta
Definition: mlp_gen.cc:36
HcalTimeSlew::BiasSetting slewFlavor_
Definition: MahiFit.h:187
MahiNnlsWorkspace nnlsWork_
Definition: MahiFit.h:169
double delay(double 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:14
bool applyTimeSlew_
Definition: MahiFit.h:186
float meanTime_
Definition: MahiFit.h:190
std::array< double, MaxSVSize > pulseM
Definition: MahiFit.h:53
double tsDelay1GeV_
Definition: MahiFit.h:188
std::array< double, MaxSVSize > pulseN
Definition: MahiFit.h:52
std::unique_ptr< FitterFuncs::PulseShapeFunctor > psfPtr_
Definition: MahiFit.h:207
unsigned int fullTSOffset
Definition: MahiFit.h:22
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
const HcalTimeSlew * hcalTimeSlewDelay_
Definition: MahiFit.h:148
unsigned int tsSize
Definition: MahiFit.h:20
std::array< double, MaxSVSize > pulseP
Definition: MahiFit.h:54

Member Data Documentation

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

Definition at line 194 of file MahiFit.h.

Referenced by doFit(), and setParameters().

bool MahiFit::applyTimeSlew_
private

Definition at line 186 of file MahiFit.h.

Referenced by setParameters(), and updatePulseShape().

int MahiFit::bxOffsetConf_
private

Definition at line 203 of file MahiFit.h.

Referenced by doFit(), and setParameters().

unsigned int MahiFit::bxSizeConf_
private

Definition at line 202 of file MahiFit.h.

Referenced by doFit(), and setParameters().

float MahiFit::chiSqSwitch_
private

Definition at line 184 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

int MahiFit::cntsetPulseShape_
private

Definition at line 206 of file MahiFit.h.

Referenced by resetPulseShapeTemplate().

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

Definition at line 147 of file MahiFit.h.

Referenced by setPulseShapeTemplate().

float MahiFit::deltaChiSqThresh_
private

Definition at line 199 of file MahiFit.h.

Referenced by minimize(), and setParameters().

bool MahiFit::dynamicPed_
private

Definition at line 182 of file MahiFit.h.

Referenced by doFit(), and setParameters().

const unsigned int MahiFit::fullTSofInterest_
private

Definition at line 173 of file MahiFit.h.

Referenced by phase1Apply().

const unsigned int MahiFit::fullTSSize_
private

Definition at line 172 of file MahiFit.h.

const HcalTimeSlew* MahiFit::hcalTimeSlewDelay_ = 0

Definition at line 148 of file MahiFit.h.

Referenced by setPulseShapeTemplate(), and updatePulseShape().

float MahiFit::meanTime_
private

Definition at line 190 of file MahiFit.h.

Referenced by setParameters(), and updatePulseShape().

int MahiFit::nMaxItersMin_
private

Definition at line 196 of file MahiFit.h.

Referenced by minimize(), and setParameters().

int MahiFit::nMaxItersNNLS_
private

Definition at line 197 of file MahiFit.h.

Referenced by nnls(), and setParameters().

float MahiFit::nnlsThresh_
private

Definition at line 200 of file MahiFit.h.

Referenced by nnls(), and setParameters().

MahiNnlsWorkspace MahiFit::nnlsWork_
mutableprivate
int MahiFit::pedestalBX_ = 100
staticprivate

Definition at line 175 of file MahiFit.h.

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

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

Definition at line 208 of file MahiFit.h.

Referenced by resetPulseShapeTemplate().

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

Definition at line 207 of file MahiFit.h.

Referenced by resetPulseShapeTemplate(), and updatePulseShape().

HcalTimeSlew::BiasSetting MahiFit::slewFlavor_
private

Definition at line 187 of file MahiFit.h.

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

float MahiFit::timeLimit_ = 12.5
staticprivate

Definition at line 179 of file MahiFit.h.

Referenced by calculateArrivalTime().

float MahiFit::timeSigmaHPD_
private

Definition at line 191 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

float MahiFit::timeSigmaSiPM_
private

Definition at line 192 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

float MahiFit::ts4Thresh_
private

Definition at line 183 of file MahiFit.h.

Referenced by phase1Apply(), and setParameters().

double MahiFit::tsDelay1GeV_ =0
private

Definition at line 188 of file MahiFit.h.

Referenced by setPulseShapeTemplate(), and updatePulseShape().