Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:30:07

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