PulseChiSqSNNLS.h

Go to the documentation of this file.

#ifndef PulseChiSqSNNLS_h
#define PulseChiSqSNNLS_h

#define EIGEN_NO_DEBUG  // kill throws in eigen code
#include "RecoLocalCalo/EcalRecAlgos/interface/EigenMatrixTypes.h"

#include <set>
#include <array>

class PulseChiSqSNNLS {
public:
  typedef BXVector::Index Index;

  PulseChiSqSNNLS();
  ~PulseChiSqSNNLS();

  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 SamplePulseMatrix &pulsemat() const { return _pulsemat; }
  const SampleMatrix &invcov() const { return _invcov; }

  const PulseVector &X() const { return _ampvecmin; }
  const PulseVector &Errors() const { return _errvec; }
  const BXVector &BXs() const { return _bxsmin; }

  double ChiSq() const { return _chisq; }
  void disableErrorCalculation() { _computeErrors = false; }
  void setMaxIters(int n) { _maxiters = n; }
  void setMaxIterWarnings(bool b) { _maxiterwarnings = b; }

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

  SampleVector _sampvec;
  SampleMatrix _invcov;
  SamplePulseMatrix _pulsemat;
  PulseVector _ampvec;
  PulseVector _errvec;
  PulseVector _ampvecmin;

  SampleDecompLLT _covdecomp;
  SampleMatrix _covdecompLinv;
  PulseMatrix _topleft_work;
  PulseDecompLDLT _pulsedecomp;

  BXVector _bxs;
  BXVector _bxsmin;
  unsigned int _npulsetot;
  unsigned int _nP;

  SamplePulseMatrix invcovp;
  PulseMatrix aTamat;
  PulseVector aTbvec;
  PulseVector updatework;

  PulseVector ampvecpermtest;

  double _chisq;
  bool _computeErrors;
  int _maxiters;
  bool _maxiterwarnings;
};

#endif