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

typedef BXVector::Index PulseChiSqSNNLS::Index

Definition at line 13 of file PulseChiSqSNNLS.h.

Constructor & Destructor Documentation

PulseChiSqSNNLS::PulseChiSqSNNLS ( )

Definition at line 75 of file PulseChiSqSNNLS.cc.

75  :
76  _chisq(0.),
77  _computeErrors(true),
78  _maxiters(50),
79  _maxiterwarnings(true)
80 {
81 
82 }
PulseChiSqSNNLS::~PulseChiSqSNNLS ( )

Definition at line 84 of file PulseChiSqSNNLS.cc.

84  {
85 
86 }

Member Function Documentation

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

Definition at line 26 of file PulseChiSqSNNLS.h.

References _bxsmin.

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

26 { return _bxsmin; }
double PulseChiSqSNNLS::ChiSq ( ) const
inline

Definition at line 28 of file PulseChiSqSNNLS.h.

References _chisq.

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

28 { return _chisq; }
double PulseChiSqSNNLS::ComputeApproxUncertainty ( unsigned int  ipulse)
protected

Definition at line 314 of file PulseChiSqSNNLS.cc.

References _covdecomp, and _pulsemat.

Referenced by DoFit(), and setMaxIterWarnings().

314  {
315  //compute approximate uncertainties
316  //(using 1/second derivative since full Hessian is not meaningful in
317  //presence of positive amplitude boundaries.)
318 
319  return 1./_covdecomp.matrixL().solve(_pulsemat.col(ipulse)).norm();
320 
321 }
SamplePulseMatrix _pulsemat
SampleDecompLLT _covdecomp
double PulseChiSqSNNLS::ComputeChiSq ( )
protected

Definition at line 305 of file PulseChiSqSNNLS.cc.

References _ampvec, _covdecomp, _pulsemat, and _sampvec.

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

305  {
306 
307 // SampleVector resvec = _pulsemat*_ampvec - _sampvec;
308 // return resvec.transpose()*_covdecomp.solve(resvec);
309 
310  return _covdecomp.matrixL().solve(_pulsemat*_ampvec - _sampvec).squaredNorm();
311 
312 }
SamplePulseMatrix _pulsemat
SampleDecompLLT _covdecomp
SampleVector _sampvec
PulseVector _ampvec
void PulseChiSqSNNLS::disableErrorCalculation ( )
inline
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 88 of file PulseChiSqSNNLS.cc.

References _ampvec, _ampvecmin, _bxs, _bxsmin, _chisq, _computeErrors, _errvec, _nP, _npulsetot, _pulsemat, _sampvec, funct::abs(), aTamat, ComputeApproxUncertainty(), ComputeChiSq(), mps_fire::i, RecoTauDiscriminantConfiguration::mask, SiStripPI::max, Minimize(), ngains, NNLSUnconstrainParameter(), PFRecoTauDiscriminationByIsolation_cfi::offset, muonCSCDigis_cfi::pedestal, BXVector< T >::resize(), hcalTTPDigis_cfi::samples, mathSSE::sqrt(), mps_update::status, and std::swap().

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

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

Definition at line 25 of file PulseChiSqSNNLS.h.

References _errvec.

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

25 { return _errvec; }
PulseVector _errvec
const SampleMatrix& PulseChiSqSNNLS::invcov ( ) const
inline

Definition at line 22 of file PulseChiSqSNNLS.h.

References _invcov.

22 { return _invcov; }
SampleMatrix _invcov
bool PulseChiSqSNNLS::Minimize ( const SampleMatrix samplecov,
const FullSampleMatrix fullpulsecov 
)
protected

Definition at line 235 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(), and setMaxIterWarnings().

235  {
236 
237  const unsigned int npulse = _bxs.rows();
238 
239  int iter = 0;
240  bool status = false;
241  while (true) {
242 
243  if (iter>=_maxiters) {
244  if (_maxiterwarnings) {
245  LogDebug("PulseChiSqSNNLS::Minimize") << "Max Iterations reached at iter " << iter;
246  }
247  break;
248  }
249 
250  status = updateCov(samplecov,fullpulsecov);
251  if (!status) break;
252  if (npulse>1) {
253  status = NNLS();
254  }
255  else {
256  //special case for one pulse fit (performance optimized)
257  status = OnePulseMinimize();
258  }
259  if (!status) break;
260 
261  double chisqnow = ComputeChiSq();
262  double deltachisq = chisqnow-_chisq;
263 
264  _chisq = chisqnow;
265  if (std::abs(deltachisq)<1e-3) {
266  break;
267  }
268  ++iter;
269  }
270 
271  return status;
272 
273 }
#define LogDebug(id)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool updateCov(const SampleMatrix &samplecov, const FullSampleMatrix &fullpulsecov)
bool PulseChiSqSNNLS::NNLS ( )
protected

Definition at line 323 of file PulseChiSqSNNLS.cc.

References _ampvec, _bxs, _covdecomp, _nP, _pulsemat, _sampvec, ampvecpermtest, aTamat, aTbvec, constexpr, MillePedeFileConverter_cfg::e, eigen_solve_submatrix(), mps_fire::i, invcovp, LogDebug, SiStripPI::max, min(), NNLSConstrainParameter(), NNLSUnconstrainParameter(), particleFlowDisplacedVertex_cfi::ratio, electronIdCutBased_cfi::threshold, and updatework.

Referenced by Minimize(), and setMaxIterWarnings().

323  {
324 
325  //Fast NNLS (fnnls) algorithm as per http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.9203&rep=rep1&type=pdf
326 
327  const unsigned int npulse = _bxs.rows();
328  constexpr unsigned int nsamples = SampleVector::RowsAtCompileTime;
329 
330  invcovp = _covdecomp.matrixL().solve(_pulsemat);
331  aTamat.noalias() = invcovp.transpose().lazyProduct(invcovp);
332  aTbvec.noalias() = invcovp.transpose().lazyProduct(_covdecomp.matrixL().solve(_sampvec));
333 
334  int iter = 0;
335  Index idxwmax = 0;
336  double wmax = 0.0;
337  double threshold = 1e-11;
338  while (true) {
339  //can only perform this step if solution is guaranteed viable
340  if (iter>0 || _nP==0) {
341  if ( _nP==std::min(npulse,nsamples) ) break;
342 
343  const unsigned int nActive = npulse - _nP;
344 
346  Index idxwmaxprev = idxwmax;
347  double wmaxprev = wmax;
348  wmax = updatework.tail(nActive).maxCoeff(&idxwmax);
349 
350  //convergence
351  if (wmax<threshold || (idxwmax==idxwmaxprev && wmax==wmaxprev)) break;
352 
353  //worst case protection
354  if (iter>=500) {
355  LogDebug("PulseChiSqSNNLS::NNLS()") << "Max Iterations reached at iter " << iter;
356  break;
357  }
358 
359  //unconstrain parameter
360  Index idxp = _nP + idxwmax;
362  }
363 
364 
365  while (true) {
366  //printf("iter in, idxsP = %i\n",int(_idxsP.size()));
367 
368  if (_nP==0) break;
369 
371 
372  //solve for unconstrained parameters
373  //need to have specialized function to call optimized versions
374  // of matrix solver... this is truly amazing...
376 
377  //check solution
378  bool positive = true;
379  for (unsigned int i = 0; i < _nP; ++i)
380  positive &= (ampvecpermtest(i) > 0);
381  if (positive) {
382  _ampvec.head(_nP) = ampvecpermtest.head(_nP);
383  break;
384  }
385 
386  //update parameter vector
387  Index minratioidx=0;
388 
389  // no realizable optimization here (because it autovectorizes!)
390  double minratio = std::numeric_limits<double>::max();
391  for (unsigned int ipulse=0; ipulse<_nP; ++ipulse) {
392  if (ampvecpermtest.coeff(ipulse)<=0.) {
393  const double c_ampvec = _ampvec.coeff(ipulse);
394  const double ratio = c_ampvec/(c_ampvec-ampvecpermtest.coeff(ipulse));
395  if (ratio<minratio) {
396  minratio = ratio;
397  minratioidx = ipulse;
398  }
399  }
400  }
401 
402  _ampvec.head(_nP) += minratio*(ampvecpermtest.head(_nP)- _ampvec.head(_nP));
403 
404  //avoid numerical problems with later ==0. check
405  _ampvec.coeffRef(minratioidx) = 0.;
406 
407  //printf("removing index %i, orig idx %i\n",int(minratioidx),int(_bxs.coeff(minratioidx)));
408  NNLSConstrainParameter(minratioidx);
409  }
410  ++iter;
411 
412  //adaptive convergence threshold to avoid infinite loops but still
413  //ensure best value is used
414  if (iter % 16 == 0) {
415  threshold *= 2;
416  }
417  }
418 
419  return true;
420 
421 
422 }
#define LogDebug(id)
PulseMatrix aTamat
void NNLSUnconstrainParameter(Index idxp)
unsigned int _nP
PulseVector aTbvec
PulseVector ampvecpermtest
#define constexpr
SamplePulseMatrix _pulsemat
SampleDecompLLT _covdecomp
T min(T a, T b)
Definition: MathUtil.h:58
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
void PulseChiSqSNNLS::NNLSConstrainParameter ( Index  minratioidx)
protected

Definition at line 435 of file PulseChiSqSNNLS.cc.

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

Referenced by NNLS(), and setMaxIterWarnings().

435  {
436  aTamat.col(_nP-1).swap(aTamat.col(minratioidx));
437  aTamat.row(_nP-1).swap(aTamat.row(minratioidx));
438  _pulsemat.col(_nP-1).swap(_pulsemat.col(minratioidx));
439  std::swap(aTbvec.coeffRef(_nP-1),aTbvec.coeffRef(minratioidx));
440  std::swap(_ampvec.coeffRef(_nP-1),_ampvec.coeffRef(minratioidx));
441  std::swap(_bxs.coeffRef(_nP-1),_bxs.coeffRef(minratioidx));
442  --_nP;
443 
444 }
PulseMatrix aTamat
unsigned int _nP
PulseVector aTbvec
SamplePulseMatrix _pulsemat
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
PulseVector _ampvec
void PulseChiSqSNNLS::NNLSUnconstrainParameter ( Index  idxp)
protected

Definition at line 424 of file PulseChiSqSNNLS.cc.

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

Referenced by DoFit(), NNLS(), and setMaxIterWarnings().

424  {
425 
426  aTamat.col(_nP).swap(aTamat.col(idxp));
427  aTamat.row(_nP).swap(aTamat.row(idxp));
428  _pulsemat.col(_nP).swap(_pulsemat.col(idxp));
429  std::swap(aTbvec.coeffRef(_nP),aTbvec.coeffRef(idxp));
430  std::swap(_ampvec.coeffRef(_nP),_ampvec.coeffRef(idxp));
431  std::swap(_bxs.coeffRef(_nP),_bxs.coeffRef(idxp));
432  ++_nP;
433 }
PulseMatrix aTamat
unsigned int _nP
PulseVector aTbvec
SamplePulseMatrix _pulsemat
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
PulseVector _ampvec
bool PulseChiSqSNNLS::OnePulseMinimize ( )
protected

Definition at line 446 of file PulseChiSqSNNLS.cc.

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

Referenced by Minimize(), and setMaxIterWarnings().

446  {
447 
448  //Fast NNLS (fnnls) algorithm as per http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.9203&rep=rep1&type=pdf
449 
450 // const unsigned int npulse = 1;
451 
452  invcovp = _covdecomp.matrixL().solve(_pulsemat);
453 // aTamat = invcovp.transpose()*invcovp;
454 // aTbvec = invcovp.transpose()*_covdecomp.matrixL().solve(_sampvec);
455 
456  SingleMatrix aTamatval = invcovp.transpose()*invcovp;
457  SingleVector aTbvecval = invcovp.transpose()*_covdecomp.matrixL().solve(_sampvec);
458  _ampvec.coeffRef(0) = std::max(0.,aTbvecval.coeff(0)/aTamatval.coeff(0));
459 
460  return true;
461 
462 }
SamplePulseMatrix _pulsemat
Eigen::Matrix< double, 1, 1 > SingleVector
SampleDecompLLT _covdecomp
SampleVector _sampvec
Eigen::Matrix< double, 1, 1 > SingleMatrix
SamplePulseMatrix invcovp
PulseVector _ampvec
const SamplePulseMatrix& PulseChiSqSNNLS::pulsemat ( ) const
inline

Definition at line 21 of file PulseChiSqSNNLS.h.

References _pulsemat.

21 { return _pulsemat; }
SamplePulseMatrix _pulsemat
void PulseChiSqSNNLS::setMaxIters ( int  n)
inline
void PulseChiSqSNNLS::setMaxIterWarnings ( bool  b)
inline
bool PulseChiSqSNNLS::updateCov ( const SampleMatrix samplecov,
const FullSampleMatrix fullpulsecov 
)
protected

Definition at line 275 of file PulseChiSqSNNLS.cc.

References _ampvec, _bxs, _covdecomp, _invcov, funct::abs(), SiStripPI::max, PFRecoTauDiscriminationByIsolation_cfi::offset, and mps_update::status.

Referenced by Minimize(), and setMaxIterWarnings().

275  {
276 
277  const unsigned int nsample = SampleVector::RowsAtCompileTime;
278  const unsigned int npulse = _bxs.rows();
279 
280  _invcov = samplecov; //
281 
282  for (unsigned int ipulse=0; ipulse<npulse; ++ipulse) {
283  if (_ampvec.coeff(ipulse)==0.) continue;
284  int bx = _bxs.coeff(ipulse);
285  if (std::abs(bx)>=100) continue; //no contribution to covariance from pedestal or saturation/slew step correction
286 
287  int firstsamplet = std::max(0,bx + 3);
288  int offset = 7-3-bx;
289 
290  const double ampveccoef = _ampvec.coeff(ipulse);
291  const double ampsq = ampveccoef*ampveccoef;
292 
293  const unsigned int nsamplepulse = nsample-firstsamplet;
294  _invcov.block(firstsamplet,firstsamplet,nsamplepulse,nsamplepulse) +=
295  ampsq*fullpulsecov.block(firstsamplet+offset,firstsamplet+offset,nsamplepulse,nsamplepulse);
296  }
297 
298  _covdecomp.compute(_invcov);
299 
300  bool status = true;
301  return status;
302 
303 }
SampleDecompLLT _covdecomp
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
SampleMatrix _invcov
PulseVector _ampvec
const PulseVector& PulseChiSqSNNLS::X ( ) const
inline

Definition at line 24 of file PulseChiSqSNNLS.h.

References _ampvecmin.

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

24 { return _ampvecmin; }
PulseVector _ampvecmin

Member Data Documentation

PulseVector PulseChiSqSNNLS::_ampvec
protected
PulseVector PulseChiSqSNNLS::_ampvecmin
protected

Definition at line 50 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), and X().

BXVector PulseChiSqSNNLS::_bxs
protected
BXVector PulseChiSqSNNLS::_bxsmin
protected

Definition at line 58 of file PulseChiSqSNNLS.h.

Referenced by BXs(), and DoFit().

double PulseChiSqSNNLS::_chisq
protected

Definition at line 69 of file PulseChiSqSNNLS.h.

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

bool PulseChiSqSNNLS::_computeErrors
protected

Definition at line 70 of file PulseChiSqSNNLS.h.

Referenced by disableErrorCalculation(), and DoFit().

SampleDecompLLT PulseChiSqSNNLS::_covdecomp
protected
SampleMatrix PulseChiSqSNNLS::_covdecompLinv
protected

Definition at line 53 of file PulseChiSqSNNLS.h.

PulseVector PulseChiSqSNNLS::_errvec
protected

Definition at line 49 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), and Errors().

SampleMatrix PulseChiSqSNNLS::_invcov
protected

Definition at line 46 of file PulseChiSqSNNLS.h.

Referenced by invcov(), and updateCov().

int PulseChiSqSNNLS::_maxiters
protected

Definition at line 71 of file PulseChiSqSNNLS.h.

Referenced by Minimize(), and setMaxIters().

bool PulseChiSqSNNLS::_maxiterwarnings
protected

Definition at line 72 of file PulseChiSqSNNLS.h.

Referenced by Minimize(), and setMaxIterWarnings().

unsigned int PulseChiSqSNNLS::_nP
protected

Definition at line 60 of file PulseChiSqSNNLS.h.

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

unsigned int PulseChiSqSNNLS::_npulsetot
protected

Definition at line 59 of file PulseChiSqSNNLS.h.

Referenced by DoFit().

PulseDecompLDLT PulseChiSqSNNLS::_pulsedecomp
protected

Definition at line 55 of file PulseChiSqSNNLS.h.

SamplePulseMatrix PulseChiSqSNNLS::_pulsemat
protected
SampleVector PulseChiSqSNNLS::_sampvec
protected

Definition at line 45 of file PulseChiSqSNNLS.h.

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

PulseMatrix PulseChiSqSNNLS::_topleft_work
protected

Definition at line 54 of file PulseChiSqSNNLS.h.

PulseVector PulseChiSqSNNLS::ampvecpermtest
protected

Definition at line 67 of file PulseChiSqSNNLS.h.

Referenced by NNLS().

PulseMatrix PulseChiSqSNNLS::aTamat
protected

Definition at line 63 of file PulseChiSqSNNLS.h.

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

PulseVector PulseChiSqSNNLS::aTbvec
protected

Definition at line 64 of file PulseChiSqSNNLS.h.

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

SamplePulseMatrix PulseChiSqSNNLS::invcovp
protected

Definition at line 62 of file PulseChiSqSNNLS.h.

Referenced by NNLS(), and OnePulseMinimize().

PulseVector PulseChiSqSNNLS::updatework
protected

Definition at line 65 of file PulseChiSqSNNLS.h.

Referenced by NNLS().