CMS 3D CMS Logo

List of all members | Public Types | Public Member Functions | Protected Member Functions | Protected Attributes
PulseChiSqSNNLS Class Reference

#include <PulseChiSqSNNLS.h>

Public Types

typedef BXVector::Index Index
 

Public Member Functions

const BXVectorBXs () const
 
double ChiSq () const
 
void disableErrorCalculation ()
 
bool DoFit (const SampleVector &samples, const SampleMatrix &samplecov, const BXVector &bxs, const FullSampleVector &fullpulse, const FullSampleMatrix &fullpulsecov, const SampleGainVector &gains=-1 *SampleGainVector::Ones(), const SampleGainVector &badSamples=SampleGainVector::Zero())
 
const PulseVectorErrors () const
 
const SampleMatrixinvcov () const
 
 PulseChiSqSNNLS ()
 
const SamplePulseMatrixpulsemat () const
 
void setMaxIters (int n)
 
void setMaxIterWarnings (bool b)
 
const PulseVectorX () const
 
 ~PulseChiSqSNNLS ()
 

Protected Member Functions

double ComputeApproxUncertainty (unsigned int ipulse)
 
double ComputeChiSq ()
 
bool Minimize (const SampleMatrix &samplecov, const FullSampleMatrix &fullpulsecov)
 
bool NNLS ()
 
void NNLSConstrainParameter (Index minratioidx)
 
void NNLSUnconstrainParameter (Index idxp)
 
bool OnePulseMinimize ()
 
bool updateCov (const SampleMatrix &samplecov, const FullSampleMatrix &fullpulsecov)
 

Protected Attributes

PulseVector _ampvec
 
PulseVector _ampvecmin
 
BXVector _bxs
 
BXVector _bxsmin
 
double _chisq
 
bool _computeErrors
 
SampleDecompLLT _covdecomp
 
SampleMatrix _covdecompLinv
 
PulseVector _errvec
 
SampleMatrix _invcov
 
int _maxiters
 
bool _maxiterwarnings
 
unsigned int _nP
 
unsigned int _npulsetot
 
PulseDecompLDLT _pulsedecomp
 
SamplePulseMatrix _pulsemat
 
SampleVector _sampvec
 
PulseMatrix _topleft_work
 
PulseVector ampvecpermtest
 
PulseMatrix aTamat
 
PulseVector aTbvec
 
SamplePulseMatrix invcovp
 
PulseVector updatework
 

Detailed Description

Definition at line 10 of file PulseChiSqSNNLS.h.

Member Typedef Documentation

◆ Index

typedef BXVector::Index PulseChiSqSNNLS::Index

Definition at line 12 of file PulseChiSqSNNLS.h.

Constructor & Destructor Documentation

◆ PulseChiSqSNNLS()

PulseChiSqSNNLS::PulseChiSqSNNLS ( )

Definition at line 55 of file PulseChiSqSNNLS.cc.

◆ ~PulseChiSqSNNLS()

PulseChiSqSNNLS::~PulseChiSqSNNLS ( )

Definition at line 57 of file PulseChiSqSNNLS.cc.

57 {}

Member Function Documentation

◆ BXs()

const BXVector& PulseChiSqSNNLS::BXs ( void  ) const
inline

Definition at line 30 of file PulseChiSqSNNLS.h.

References _bxsmin.

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

30 { return _bxsmin; }

◆ ChiSq()

double PulseChiSqSNNLS::ChiSq ( ) const
inline

Definition at line 32 of file PulseChiSqSNNLS.h.

References _chisq.

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

32 { return _chisq; }

◆ ComputeApproxUncertainty()

double PulseChiSqSNNLS::ComputeApproxUncertainty ( unsigned int  ipulse)
protected

Definition at line 289 of file PulseChiSqSNNLS.cc.

References _covdecomp, and _pulsemat.

Referenced by DoFit().

289  {
290  //compute approximate uncertainties
291  //(using 1/second derivative since full Hessian is not meaningful in
292  //presence of positive amplitude boundaries.)
293 
294  return 1. / _covdecomp.matrixL().solve(_pulsemat.col(ipulse)).norm();
295 }
SamplePulseMatrix _pulsemat
SampleDecompLLT _covdecomp

◆ ComputeChiSq()

double PulseChiSqSNNLS::ComputeChiSq ( )
protected

Definition at line 282 of file PulseChiSqSNNLS.cc.

References _ampvec, _covdecomp, _pulsemat, and _sampvec.

Referenced by DoFit(), and Minimize().

282  {
283  // SampleVector resvec = _pulsemat*_ampvec - _sampvec;
284  // return resvec.transpose()*_covdecomp.solve(resvec);
285 
286  return _covdecomp.matrixL().solve(_pulsemat * _ampvec - _sampvec).squaredNorm();
287 }
SamplePulseMatrix _pulsemat
SampleDecompLLT _covdecomp
SampleVector _sampvec
PulseVector _ampvec

◆ disableErrorCalculation()

void PulseChiSqSNNLS::disableErrorCalculation ( )
inline

◆ DoFit()

bool PulseChiSqSNNLS::DoFit ( const SampleVector samples,
const SampleMatrix samplecov,
const BXVector bxs,
const FullSampleVector fullpulse,
const FullSampleMatrix fullpulsecov,
const SampleGainVector gains = -1 * SampleGainVector::Ones(),
const SampleGainVector badSamples = SampleGainVector::Zero() 
)

Definition at line 59 of file PulseChiSqSNNLS.cc.

References _ampvec, _ampvecmin, _bxs, _bxsmin, _chisq, _computeErrors, _errvec, _nP, _npulsetot, _pulsemat, _sampvec, funct::abs(), aTamat, nano_mu_digi_cff::bx, ComputeApproxUncertainty(), ComputeChiSq(), ecalph2::gains, mps_fire::i, ALPAKA_ACCELERATOR_NAMESPACE::pixelClustering::pixelStatus::mask, WZElectronSkims53X_cff::max, Minimize(), NNLSUnconstrainParameter(), HLT_IsoTrack_cff::offset, Hcal_Conditions_forGlobalTag_cff::pedestal, BXVector< T >::resize(), EgammaValidation_cff::samples, mathSSE::sqrt(), mps_update::status, and edm::swap().

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

65  {
66  int npulse = bxs.rows();
67 
68  _sampvec = samples;
69  _bxs = bxs;
70  _pulsemat.resize(Eigen::NoChange, npulse);
71 
72  //construct dynamic pedestals if applicable
73  int ngains = gains.maxCoeff() + 1;
74  int nPedestals = 0;
75  for (int gainidx = 0; gainidx < ngains; ++gainidx) {
76  SampleGainVector mask = gainidx * SampleGainVector::Ones();
77  SampleVector pedestal = (gains.array() == mask.array()).cast<SampleVector::value_type>();
78  if (pedestal.maxCoeff() > 0.) {
79  ++nPedestals;
80  _bxs.resize(npulse + nPedestals);
81  _bxs[npulse + nPedestals - 1] = 100 + gainidx; //bx values >=100 indicate dynamic pedestals
82  _pulsemat.resize(Eigen::NoChange, npulse + nPedestals);
83  _pulsemat.col(npulse + nPedestals - 1) = pedestal;
84  }
85  }
86 
87  //construct negative step functions for saturated or potentially slew-rate-limited samples
88  for (int isample = 0; isample < SampleVector::RowsAtCompileTime; ++isample) {
89  if (badSamples.coeff(isample) > 0) {
90  SampleVector step = SampleVector::Zero();
91  //step correction has negative sign for saturated or slew-limited samples which have been forced to zero
92  step[isample] = -1.;
93 
94  ++nPedestals;
95  _bxs.resize(npulse + nPedestals);
96  _bxs[npulse + nPedestals - 1] =
97  -100 - isample; //bx values <=-100 indicate step corrections for saturated or slew-limited samples
98  _pulsemat.resize(Eigen::NoChange, npulse + nPedestals);
99  _pulsemat.col(npulse + nPedestals - 1) = step;
100  }
101  }
102 
103  _npulsetot = npulse + nPedestals;
104 
105  _ampvec = PulseVector::Zero(_npulsetot);
106  _errvec = PulseVector::Zero(_npulsetot);
107  _nP = 0;
108  _chisq = 0.;
109 
110  if (_bxs.rows() == 1 && std::abs(_bxs.coeff(0)) < 100) {
111  _ampvec.coeffRef(0) = _sampvec.coeff(_bxs.coeff(0) + 5);
112  }
113 
114  aTamat.resize(_npulsetot, _npulsetot);
115 
116  //initialize pulse template matrix
117  for (int ipulse = 0; ipulse < npulse; ++ipulse) {
118  int bx = _bxs.coeff(ipulse);
119  int offset = 7 - 3 - bx;
120  _pulsemat.col(ipulse) = fullpulse.segment<SampleVector::RowsAtCompileTime>(offset);
121  }
122 
123  //unconstrain pedestals already for first iteration since they should always be non-zero
124  if (nPedestals > 0) {
125  for (int i = 0; i < _bxs.rows(); ++i) {
126  int bx = _bxs.coeff(i);
127  if (bx >= 100) {
129  }
130  }
131  }
132 
133  //do the actual fit
134  bool status = Minimize(samplecov, fullpulsecov);
136  _bxsmin = _bxs;
137 
138  if (!status)
139  return status;
140 
141  if (!_computeErrors)
142  return status;
143 
144  //compute MINOS-like uncertainties for in-time amplitude
145  bool foundintime = false;
146  unsigned int ipulseintime = 0;
147  for (unsigned int ipulse = 0; ipulse < _npulsetot; ++ipulse) {
148  if (_bxs.coeff(ipulse) == 0) {
149  ipulseintime = ipulse;
150  foundintime = true;
151  break;
152  }
153  }
154  if (!foundintime)
155  return status;
156 
157  const unsigned int ipulseintimemin = ipulseintime;
158 
159  double approxerr = ComputeApproxUncertainty(ipulseintime);
160  double chisq0 = _chisq;
161  double x0 = _ampvecmin[ipulseintime];
162 
163  //move in time pulse first to active set if necessary
164  if (ipulseintime < _nP) {
165  _pulsemat.col(_nP - 1).swap(_pulsemat.col(ipulseintime));
166  std::swap(_ampvec.coeffRef(_nP - 1), _ampvec.coeffRef(ipulseintime));
167  std::swap(_bxs.coeffRef(_nP - 1), _bxs.coeffRef(ipulseintime));
168  ipulseintime = _nP - 1;
169  --_nP;
170  }
171 
172  SampleVector pulseintime = _pulsemat.col(ipulseintime);
173  _pulsemat.col(ipulseintime).setZero();
174 
175  //two point interpolation for upper uncertainty when amplitude is away from boundary
176  double xplus100 = x0 + approxerr;
177  _ampvec.coeffRef(ipulseintime) = xplus100;
178  _sampvec = samples - _ampvec.coeff(ipulseintime) * pulseintime;
179  status &= Minimize(samplecov, fullpulsecov);
180  if (!status)
181  return status;
182  double chisqplus100 = ComputeChiSq();
183 
184  double sigmaplus = std::abs(xplus100 - x0) / sqrt(chisqplus100 - chisq0);
185 
186  //if amplitude is sufficiently far from the boundary, compute also the lower uncertainty and average them
187  if ((x0 / sigmaplus) > 0.5) {
188  for (unsigned int ipulse = 0; ipulse < _npulsetot; ++ipulse) {
189  if (_bxs.coeff(ipulse) == 0) {
190  ipulseintime = ipulse;
191  break;
192  }
193  }
194  double xminus100 = std::max(0., x0 - approxerr);
195  _ampvec.coeffRef(ipulseintime) = xminus100;
196  _sampvec = samples - _ampvec.coeff(ipulseintime) * pulseintime;
197  status &= Minimize(samplecov, fullpulsecov);
198  if (!status)
199  return status;
200  double chisqminus100 = ComputeChiSq();
201 
202  double sigmaminus = std::abs(xminus100 - x0) / sqrt(chisqminus100 - chisq0);
203  _errvec[ipulseintimemin] = 0.5 * (sigmaplus + sigmaminus);
204 
205  } else {
206  _errvec[ipulseintimemin] = sigmaplus;
207  }
208 
209  _chisq = chisq0;
210 
211  return status;
212 }
PulseMatrix aTamat
void NNLSUnconstrainParameter(Index idxp)
unsigned int _npulsetot
unsigned int _nP
PulseVector _ampvecmin
void swap(Association< C > &lhs, Association< C > &rhs)
Definition: Association.h:112
SamplePulseMatrix _pulsemat
double ComputeApproxUncertainty(unsigned int ipulse)
T sqrt(T t)
Definition: SSEVec.h:23
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool Minimize(const SampleMatrix &samplecov, const FullSampleMatrix &fullpulsecov)
constexpr float gains[NGAINS]
Definition: EcalConstants.h:20
SampleVector _sampvec
void resize(int bx, unsigned size)
Eigen::Matrix< double, SampleVectorSize, 1 > SampleVector
Eigen::Matrix< char, SampleVectorSize, 1 > SampleGainVector
PulseVector _errvec
step
Definition: StallMonitor.cc:83
PulseVector _ampvec

◆ Errors()

const PulseVector& PulseChiSqSNNLS::Errors ( ) const
inline

Definition at line 29 of file PulseChiSqSNNLS.h.

References _errvec.

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

29 { return _errvec; }
PulseVector _errvec

◆ invcov()

const SampleMatrix& PulseChiSqSNNLS::invcov ( ) const
inline

Definition at line 26 of file PulseChiSqSNNLS.h.

References _invcov.

26 { return _invcov; }
SampleMatrix _invcov

◆ Minimize()

bool PulseChiSqSNNLS::Minimize ( const SampleMatrix samplecov,
const FullSampleMatrix fullpulsecov 
)
protected

Definition at line 214 of file PulseChiSqSNNLS.cc.

References _bxs, _chisq, _maxiters, _maxiterwarnings, funct::abs(), ComputeChiSq(), MillePedeFileConverter_cfg::e, LogDebug, NNLS(), OnePulseMinimize(), mps_update::status, and updateCov().

Referenced by DoFit().

214  {
215  const unsigned int npulse = _bxs.rows();
216 
217  int iter = 0;
218  bool status = false;
219  while (true) {
220  if (iter >= _maxiters) {
221  if (_maxiterwarnings) {
222  LogDebug("PulseChiSqSNNLS::Minimize") << "Max Iterations reached at iter " << iter;
223  }
224  break;
225  }
226 
227  status = updateCov(samplecov, fullpulsecov);
228  if (!status)
229  break;
230  if (npulse > 1) {
231  status = NNLS();
232  } else {
233  //special case for one pulse fit (performance optimized)
235  }
236  if (!status)
237  break;
238 
239  double chisqnow = ComputeChiSq();
240  double deltachisq = chisqnow - _chisq;
241 
242  _chisq = chisqnow;
243  if (std::abs(deltachisq) < 1e-3) {
244  break;
245  }
246  ++iter;
247  }
248 
249  return status;
250 }
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool updateCov(const SampleMatrix &samplecov, const FullSampleMatrix &fullpulsecov)
#define LogDebug(id)

◆ NNLS()

bool PulseChiSqSNNLS::NNLS ( )
protected

Definition at line 297 of file PulseChiSqSNNLS.cc.

References _ampvec, _bxs, _covdecomp, _nP, _pulsemat, _sampvec, ampvecpermtest, aTamat, aTbvec, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), MillePedeFileConverter_cfg::e, eigen_solve_submatrix(), mps_fire::i, invcovp, LogDebug, WZElectronSkims53X_cff::max, SiStripPI::min, NNLSConstrainParameter(), NNLSUnconstrainParameter(), DiMuonV_cfg::threshold, updatework, and cms::cuda::wmax.

Referenced by Minimize().

297  {
298  //Fast NNLS (fnnls) algorithm as per http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.9203&rep=rep1&type=pdf
299 
300  const unsigned int npulse = _bxs.rows();
301  constexpr unsigned int nsamples = SampleVector::RowsAtCompileTime;
302 
303  invcovp = _covdecomp.matrixL().solve(_pulsemat);
304  aTamat.noalias() = invcovp.transpose().lazyProduct(invcovp);
305  aTbvec.noalias() = invcovp.transpose().lazyProduct(_covdecomp.matrixL().solve(_sampvec));
306 
307  int iter = 0;
308  Index idxwmax = 0;
309  double wmax = 0.0;
310  double threshold = 1e-11;
311  while (true) {
312  //can only perform this step if solution is guaranteed viable
313  if (iter > 0 || _nP == 0) {
314  if (_nP == std::min(npulse, nsamples))
315  break;
316 
317  const unsigned int nActive = npulse - _nP;
318 
320  Index idxwmaxprev = idxwmax;
321  double wmaxprev = wmax;
322  wmax = updatework.tail(nActive).maxCoeff(&idxwmax);
323 
324  //convergence
325  if (wmax < threshold || (idxwmax == idxwmaxprev && wmax == wmaxprev))
326  break;
327 
328  //worst case protection
329  if (iter >= 500) {
330  LogDebug("PulseChiSqSNNLS::NNLS()") << "Max Iterations reached at iter " << iter;
331  break;
332  }
333 
334  //unconstrain parameter
335  Index idxp = _nP + idxwmax;
337  }
338 
339  while (true) {
340  //printf("iter in, idxsP = %i\n",int(_idxsP.size()));
341 
342  if (_nP == 0)
343  break;
344 
346 
347  //solve for unconstrained parameters
348  //need to have specialized function to call optimized versions
349  // of matrix solver... this is truly amazing...
351 
352  //check solution
353  bool positive = true;
354  for (unsigned int i = 0; i < _nP; ++i)
355  positive &= (ampvecpermtest(i) > 0);
356  if (positive) {
357  _ampvec.head(_nP) = ampvecpermtest.head(_nP);
358  break;
359  }
360 
361  //update parameter vector
362  Index minratioidx = 0;
363 
364  // no realizable optimization here (because it autovectorizes!)
365  double minratio = std::numeric_limits<double>::max();
366  for (unsigned int ipulse = 0; ipulse < _nP; ++ipulse) {
367  if (ampvecpermtest.coeff(ipulse) <= 0.) {
368  const double c_ampvec = _ampvec.coeff(ipulse);
369  const double ratio = c_ampvec / (c_ampvec - ampvecpermtest.coeff(ipulse));
370  if (ratio < minratio) {
371  minratio = ratio;
372  minratioidx = ipulse;
373  }
374  }
375  }
376 
377  _ampvec.head(_nP) += minratio * (ampvecpermtest.head(_nP) - _ampvec.head(_nP));
378 
379  //avoid numerical problems with later ==0. check
380  _ampvec.coeffRef(minratioidx) = 0.;
381 
382  //printf("removing index %i, orig idx %i\n",int(minratioidx),int(_bxs.coeff(minratioidx)));
383  NNLSConstrainParameter(minratioidx);
384  }
385  ++iter;
386 
387  //adaptive convergence threshold to avoid infinite loops but still
388  //ensure best value is used
389  if (iter % 16 == 0) {
390  threshold *= 2;
391  }
392  }
393 
394  return true;
395 }
PulseMatrix aTamat
void NNLSUnconstrainParameter(Index idxp)
unsigned int _nP
PulseVector aTbvec
PulseVector ampvecpermtest
SamplePulseMatrix _pulsemat
SampleDecompLLT _covdecomp
PulseVector updatework
SampleVector _sampvec
void eigen_solve_submatrix(PulseMatrix &mat, PulseVector &invec, PulseVector &outvec, unsigned NP)
BXVector::Index Index
void NNLSConstrainParameter(Index minratioidx)
SamplePulseMatrix invcovp
PulseVector _ampvec
#define LogDebug(id)
__host__ __device__ V V wmax

◆ NNLSConstrainParameter()

void PulseChiSqSNNLS::NNLSConstrainParameter ( Index  minratioidx)
protected

Definition at line 407 of file PulseChiSqSNNLS.cc.

References _ampvec, _bxs, _nP, _pulsemat, aTamat, aTbvec, and edm::swap().

Referenced by NNLS().

407  {
408  aTamat.col(_nP - 1).swap(aTamat.col(minratioidx));
409  aTamat.row(_nP - 1).swap(aTamat.row(minratioidx));
410  _pulsemat.col(_nP - 1).swap(_pulsemat.col(minratioidx));
411  std::swap(aTbvec.coeffRef(_nP - 1), aTbvec.coeffRef(minratioidx));
412  std::swap(_ampvec.coeffRef(_nP - 1), _ampvec.coeffRef(minratioidx));
413  std::swap(_bxs.coeffRef(_nP - 1), _bxs.coeffRef(minratioidx));
414  --_nP;
415 }
PulseMatrix aTamat
unsigned int _nP
PulseVector aTbvec
void swap(Association< C > &lhs, Association< C > &rhs)
Definition: Association.h:112
SamplePulseMatrix _pulsemat
PulseVector _ampvec

◆ NNLSUnconstrainParameter()

void PulseChiSqSNNLS::NNLSUnconstrainParameter ( Index  idxp)
protected

Definition at line 397 of file PulseChiSqSNNLS.cc.

References _ampvec, _bxs, _nP, _pulsemat, aTamat, aTbvec, and edm::swap().

Referenced by DoFit(), and NNLS().

397  {
398  aTamat.col(_nP).swap(aTamat.col(idxp));
399  aTamat.row(_nP).swap(aTamat.row(idxp));
400  _pulsemat.col(_nP).swap(_pulsemat.col(idxp));
401  std::swap(aTbvec.coeffRef(_nP), aTbvec.coeffRef(idxp));
402  std::swap(_ampvec.coeffRef(_nP), _ampvec.coeffRef(idxp));
403  std::swap(_bxs.coeffRef(_nP), _bxs.coeffRef(idxp));
404  ++_nP;
405 }
PulseMatrix aTamat
unsigned int _nP
PulseVector aTbvec
void swap(Association< C > &lhs, Association< C > &rhs)
Definition: Association.h:112
SamplePulseMatrix _pulsemat
PulseVector _ampvec

◆ OnePulseMinimize()

bool PulseChiSqSNNLS::OnePulseMinimize ( )
protected

Definition at line 417 of file PulseChiSqSNNLS.cc.

References _ampvec, _covdecomp, _pulsemat, _sampvec, invcovp, and WZElectronSkims53X_cff::max.

Referenced by Minimize().

417  {
418  //Fast NNLS (fnnls) algorithm as per http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.9203&rep=rep1&type=pdf
419 
420  // const unsigned int npulse = 1;
421 
422  invcovp = _covdecomp.matrixL().solve(_pulsemat);
423  // aTamat = invcovp.transpose()*invcovp;
424  // aTbvec = invcovp.transpose()*_covdecomp.matrixL().solve(_sampvec);
425 
426  SingleMatrix aTamatval = invcovp.transpose() * invcovp;
427  SingleVector aTbvecval = invcovp.transpose() * _covdecomp.matrixL().solve(_sampvec);
428  _ampvec.coeffRef(0) = std::max(0., aTbvecval.coeff(0) / aTamatval.coeff(0));
429 
430  return true;
431 }
Eigen::Matrix< double, 1, 1 > SingleVector
SamplePulseMatrix _pulsemat
SampleDecompLLT _covdecomp
SampleVector _sampvec
Eigen::Matrix< double, 1, 1 > SingleMatrix
SamplePulseMatrix invcovp
PulseVector _ampvec

◆ pulsemat()

const SamplePulseMatrix& PulseChiSqSNNLS::pulsemat ( ) const
inline

Definition at line 25 of file PulseChiSqSNNLS.h.

References _pulsemat.

25 { return _pulsemat; }
SamplePulseMatrix _pulsemat

◆ setMaxIters()

void PulseChiSqSNNLS::setMaxIters ( int  n)
inline

◆ setMaxIterWarnings()

void PulseChiSqSNNLS::setMaxIterWarnings ( bool  b)
inline

Definition at line 35 of file PulseChiSqSNNLS.h.

References _maxiterwarnings, and b.

Referenced by EcalUncalibRecHitMultiFitAlgo::EcalUncalibRecHitMultiFitAlgo().

35 { _maxiterwarnings = b; }
double b
Definition: hdecay.h:120

◆ updateCov()

bool PulseChiSqSNNLS::updateCov ( const SampleMatrix samplecov,
const FullSampleMatrix fullpulsecov 
)
protected

Definition at line 252 of file PulseChiSqSNNLS.cc.

References _ampvec, _bxs, _covdecomp, _invcov, funct::abs(), nano_mu_digi_cff::bx, WZElectronSkims53X_cff::max, HLT_IsoTrack_cff::offset, and mps_update::status.

Referenced by Minimize().

252  {
253  const unsigned int nsample = SampleVector::RowsAtCompileTime;
254  const unsigned int npulse = _bxs.rows();
255 
256  _invcov = samplecov; //
257 
258  for (unsigned int ipulse = 0; ipulse < npulse; ++ipulse) {
259  if (_ampvec.coeff(ipulse) == 0.)
260  continue;
261  int bx = _bxs.coeff(ipulse);
262  if (std::abs(bx) >= 100)
263  continue; //no contribution to covariance from pedestal or saturation/slew step correction
264 
265  int firstsamplet = std::max(0, bx + 3);
266  int offset = 7 - 3 - bx;
267 
268  const double ampveccoef = _ampvec.coeff(ipulse);
269  const double ampsq = ampveccoef * ampveccoef;
270 
271  const unsigned int nsamplepulse = nsample - firstsamplet;
272  _invcov.block(firstsamplet, firstsamplet, nsamplepulse, nsamplepulse) +=
273  ampsq * fullpulsecov.block(firstsamplet + offset, firstsamplet + offset, nsamplepulse, nsamplepulse);
274  }
275 
276  _covdecomp.compute(_invcov);
277 
278  bool status = true;
279  return status;
280 }
SampleDecompLLT _covdecomp
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
SampleMatrix _invcov
PulseVector _ampvec

◆ X()

const PulseVector& PulseChiSqSNNLS::X ( ) const
inline

Definition at line 28 of file PulseChiSqSNNLS.h.

References _ampvecmin.

Referenced by svgfig.Curve.Sample::__repr__(), and EcalUncalibRecHitMultiFitAlgo::makeRecHit().

28 { return _ampvecmin; }
PulseVector _ampvecmin

Member Data Documentation

◆ _ampvec

PulseVector PulseChiSqSNNLS::_ampvec
protected

◆ _ampvecmin

PulseVector PulseChiSqSNNLS::_ampvecmin
protected

Definition at line 52 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), and X().

◆ _bxs

BXVector PulseChiSqSNNLS::_bxs
protected

◆ _bxsmin

BXVector PulseChiSqSNNLS::_bxsmin
protected

Definition at line 60 of file PulseChiSqSNNLS.h.

Referenced by BXs(), and DoFit().

◆ _chisq

double PulseChiSqSNNLS::_chisq
protected

Definition at line 71 of file PulseChiSqSNNLS.h.

Referenced by ChiSq(), DoFit(), and Minimize().

◆ _computeErrors

bool PulseChiSqSNNLS::_computeErrors
protected

Definition at line 72 of file PulseChiSqSNNLS.h.

Referenced by disableErrorCalculation(), and DoFit().

◆ _covdecomp

SampleDecompLLT PulseChiSqSNNLS::_covdecomp
protected

◆ _covdecompLinv

SampleMatrix PulseChiSqSNNLS::_covdecompLinv
protected

Definition at line 55 of file PulseChiSqSNNLS.h.

◆ _errvec

PulseVector PulseChiSqSNNLS::_errvec
protected

Definition at line 51 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), and Errors().

◆ _invcov

SampleMatrix PulseChiSqSNNLS::_invcov
protected

Definition at line 48 of file PulseChiSqSNNLS.h.

Referenced by invcov(), and updateCov().

◆ _maxiters

int PulseChiSqSNNLS::_maxiters
protected

Definition at line 73 of file PulseChiSqSNNLS.h.

Referenced by Minimize(), and setMaxIters().

◆ _maxiterwarnings

bool PulseChiSqSNNLS::_maxiterwarnings
protected

Definition at line 74 of file PulseChiSqSNNLS.h.

Referenced by Minimize(), and setMaxIterWarnings().

◆ _nP

unsigned int PulseChiSqSNNLS::_nP
protected

Definition at line 62 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), NNLS(), NNLSConstrainParameter(), and NNLSUnconstrainParameter().

◆ _npulsetot

unsigned int PulseChiSqSNNLS::_npulsetot
protected

Definition at line 61 of file PulseChiSqSNNLS.h.

Referenced by DoFit().

◆ _pulsedecomp

PulseDecompLDLT PulseChiSqSNNLS::_pulsedecomp
protected

Definition at line 57 of file PulseChiSqSNNLS.h.

◆ _pulsemat

SamplePulseMatrix PulseChiSqSNNLS::_pulsemat
protected

◆ _sampvec

SampleVector PulseChiSqSNNLS::_sampvec
protected

Definition at line 47 of file PulseChiSqSNNLS.h.

Referenced by ComputeChiSq(), DoFit(), NNLS(), and OnePulseMinimize().

◆ _topleft_work

PulseMatrix PulseChiSqSNNLS::_topleft_work
protected

Definition at line 56 of file PulseChiSqSNNLS.h.

◆ ampvecpermtest

PulseVector PulseChiSqSNNLS::ampvecpermtest
protected

Definition at line 69 of file PulseChiSqSNNLS.h.

Referenced by NNLS().

◆ aTamat

PulseMatrix PulseChiSqSNNLS::aTamat
protected

Definition at line 65 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), NNLS(), NNLSConstrainParameter(), and NNLSUnconstrainParameter().

◆ aTbvec

PulseVector PulseChiSqSNNLS::aTbvec
protected

Definition at line 66 of file PulseChiSqSNNLS.h.

Referenced by NNLS(), NNLSConstrainParameter(), and NNLSUnconstrainParameter().

◆ invcovp

SamplePulseMatrix PulseChiSqSNNLS::invcovp
protected

Definition at line 64 of file PulseChiSqSNNLS.h.

Referenced by NNLS(), and OnePulseMinimize().

◆ updatework

PulseVector PulseChiSqSNNLS::updatework
protected

Definition at line 67 of file PulseChiSqSNNLS.h.

Referenced by NNLS().