EcalRecAlgos/interface/PulseChiSqSNNLS.h

0001 #ifndef PulseChiSqSNNLS_h
0002 #define PulseChiSqSNNLS_h
0003
0004 #define EIGEN_NO_DEBUG  // kill throws in eigen code
0005 #include "RecoLocalCalo/EcalRecAlgos/interface/EigenMatrixTypes.h"
0006
0007 #include <set>
0008 #include <array>
0009
0010 class PulseChiSqSNNLS {
0011 public:
0012   typedef BXVector::Index Index;
0013
0014   PulseChiSqSNNLS();
0015   ~PulseChiSqSNNLS();
0016
0017   bool DoFit(const SampleVector &samples,
0018              const SampleMatrix &samplecov,
0019              const BXVector &bxs,
0020              const FullSampleVector &fullpulse,
0021              const FullSampleMatrix &fullpulsecov,
0022              const SampleGainVector &gains = -1 * SampleGainVector::Ones(),
0023              const SampleGainVector &badSamples = SampleGainVector::Zero());
0024
0025   const SamplePulseMatrix &pulsemat() const { return _pulsemat; }
0026   const SampleMatrix &invcov() const { return _invcov; }
0027
0028   const PulseVector &X() const { return _ampvecmin; }
0029   const PulseVector &Errors() const { return _errvec; }
0030   const BXVector &BXs() const { return _bxsmin; }
0031
0032   double ChiSq() const { return _chisq; }
0033   void disableErrorCalculation() { _computeErrors = false; }
0034   void setMaxIters(int n) { _maxiters = n; }
0035   void setMaxIterWarnings(bool b) { _maxiterwarnings = b; }
0036
0037 protected:
0038   bool Minimize(const SampleMatrix &samplecov, const FullSampleMatrix &fullpulsecov);
0039   bool NNLS();
0040   void NNLSUnconstrainParameter(Index idxp);
0041   void NNLSConstrainParameter(Index minratioidx);
0042   bool OnePulseMinimize();
0043   bool updateCov(const SampleMatrix &samplecov, const FullSampleMatrix &fullpulsecov);
0044   double ComputeChiSq();
0045   double ComputeApproxUncertainty(unsigned int ipulse);
0046
0047   SampleVector _sampvec;
0048   SampleMatrix _invcov;
0049   SamplePulseMatrix _pulsemat;
0050   PulseVector _ampvec;
0051   PulseVector _errvec;
0052   PulseVector _ampvecmin;
0053
0054   SampleDecompLLT _covdecomp;
0055   SampleMatrix _covdecompLinv;
0056   PulseMatrix _topleft_work;
0057   PulseDecompLDLT _pulsedecomp;
0058
0059   BXVector _bxs;
0060   BXVector _bxsmin;
0061   unsigned int _npulsetot;
0062   unsigned int _nP;
0063
0064   SamplePulseMatrix invcovp;
0065   PulseMatrix aTamat;
0066   PulseVector aTbvec;
0067   PulseVector updatework;
0068
0069   PulseVector ampvecpermtest;
0070
0071   double _chisq;
0072   bool _computeErrors;
0073   int _maxiters;
0074   bool _maxiterwarnings;
0075 };
0076
0077 #endif