NoiseGenerators/interface/CorrelatedNoisifier.h

0001 #ifndef SimAlgos_CorrelatedNoisifier_h
0002 #define SimAlgos_CorrelatedNoisifier_h
0003
0004 /**
0005    \class CorrelatedNoisifier
0006
0007    \brief adds noise to the given frame.
0008
0009 Takes input correlation matrix C and from it creates
0010 an upper-triangular matrix H such that H*Htranspose = C.
0011
0012 Algorithm taken from
0013 http://cg.scs.carleton.ca/~luc/chapter_eleven.pdf, p 564-567
0014 Uses a Cholesky decomposition
0015
0016 Any input array f is "noisified" with f += H*r
0017 where r is an array of random numbers
0018
0019 The above matrix multiplication is expedited in the
0020 trivial cases of a purely diagonal or identity correlation matrix.
0021 */
0022
0023 #include "DataFormats/Math/interface/Error.h"
0024 #include <vector>
0025
0026 namespace CLHEP {
0027   class HepRandomEngine;
0028 }
0029
0030 template <class M>
0031 class CorrelatedNoisifier {
0032 public:
0033   typedef std::vector<double> VecDou;
0034
0035   CorrelatedNoisifier(const M &symCorMat);  // correlation matrix
0036
0037   CorrelatedNoisifier(int *dummy, const M &cholDecMat);  // decomposition matrix
0038
0039   virtual ~CorrelatedNoisifier();
0040
0041   void resetCorrelationMatrix(const M &symCorMat);
0042
0043   void resetCholDecompMatrix(const M &cholMat);
0044
0045   template <class T>
0046   void noisify(T &frame,  // applies random noise to frame
0047                CLHEP::HepRandomEngine *,
0048                const VecDou *rangau = nullptr) const;  // use these
0049
0050   const M &cholMat() const;  // return decomposition
0051
0052   const VecDou &vecgau() const;
0053
0054 private:
0055   static const double k_precision;  // precision to which 0 & 1 are compared
0056
0057   void init(const M &symCorMat);
0058
0059   void initChol();
0060
0061   bool computeDecomposition(const M &symCorMat);
0062
0063   bool checkDecomposition(const M &symCorMat, M &HHtDiff) const;
0064
0065   void checkOffDiagonal(const M &symCorMat);
0066
0067   mutable VecDou m_vecgau;
0068
0069   bool m_isDiagonal;
0070   bool m_isIdentity;
0071
0072   M m_H;
0073 };
0074
0075 #endif