File indexing completed on 2024-04-06 12:30:39
0001 #include "FWCore/Utilities/interface/Exception.h"
0002 #include "SimGeneral/NoiseGenerators/interface/CorrelatedNoisifier.h"
0003
0004 #include "CLHEP/Random/RandGaussQ.h"
0005
0006 template <class M>
0007 const double CorrelatedNoisifier<M>::k_precision(1.e-7);
0008
0009 template <class M>
0010 CorrelatedNoisifier<M>::~CorrelatedNoisifier() {}
0011
0012 template <class M>
0013 CorrelatedNoisifier<M>::CorrelatedNoisifier(const M &symCorMat) : m_H(ROOT::Math::SMatrixIdentity()) {
0014 init(symCorMat);
0015 }
0016
0017 template <class M>
0018 CorrelatedNoisifier<M>::CorrelatedNoisifier(int *dummy, const M &cholDecMat) : m_H(cholDecMat) {
0019 initChol();
0020 }
0021
0022 template <class M>
0023 void CorrelatedNoisifier<M>::initChol() {
0024 M symCorMat;
0025 for (unsigned int i = 0; i < m_H.kRows; ++i) {
0026 for (unsigned int j = 0; j < m_H.kRows; ++j) {
0027 double sum(0);
0028 for (unsigned int k(0); k != m_H.kRows; ++k) {
0029 const double hik(k > i ? 0 : m_H(i, k));
0030 const double htkj(k > j ? 0 : m_H(k, j));
0031 sum += hik * htkj;
0032 }
0033 symCorMat(i, j) = sum;
0034 }
0035 }
0036 checkOffDiagonal(symCorMat);
0037 }
0038
0039 template <class M>
0040 void CorrelatedNoisifier<M>::resetCorrelationMatrix(const M &symCorMat) {
0041 m_H = ROOT::Math::SMatrixIdentity();
0042 init(symCorMat);
0043 }
0044
0045 template <class M>
0046 void CorrelatedNoisifier<M>::resetCholDecompMatrix(const M &cholDecMat) {
0047 m_H = cholDecMat;
0048 initChol();
0049 }
0050
0051 template <class M>
0052 void CorrelatedNoisifier<M>::init(const M &symCorMat) {
0053 checkOffDiagonal(symCorMat);
0054
0055 const bool okDecomp(m_isIdentity || computeDecomposition(symCorMat));
0056
0057 if (!okDecomp)
0058 throw cms::Exception("CorrelatedNoisifier") << "Bad covariance matrix. " << symCorMat;
0059
0060 M HHtDiff;
0061 const bool check(checkDecomposition(symCorMat, HHtDiff));
0062 if (!check)
0063 throw cms::Exception("CorrelatedNoisifier") << "Decomposition failed, difference = " << HHtDiff;
0064 }
0065
0066 template <class M>
0067 void CorrelatedNoisifier<M>::checkOffDiagonal(const M &symCorMat) {
0068 m_isDiagonal = true;
0069 m_isIdentity = true;
0070
0071 for (unsigned int i = 0; i < m_H.kRows; i++) {
0072 if (symCorMat(i, i) < 0.)
0073 throw cms::Exception("CorrelatedNoisifier") << "Bad correlation matrix. Negative diagonal";
0074 if (m_isIdentity && fabs(symCorMat(i, i) - 1.) > k_precision)
0075 m_isIdentity = false;
0076 for (unsigned int j = 0; j < i; ++j) {
0077 if (fabs(symCorMat(i, j)) > k_precision) {
0078 m_isDiagonal = false;
0079 m_isIdentity = false;
0080 return;
0081 }
0082 }
0083 }
0084 }
0085
0086 template <class M>
0087 bool CorrelatedNoisifier<M>::computeDecomposition(const M &symCorMat) {
0088 const double s00 = std::sqrt(symCorMat(0, 0));
0089
0090 for (unsigned int i = 0; i < m_H.kRows; ++i) {
0091 m_H(i, 0) = symCorMat(i, 0) / s00;
0092 if (0 < i) {
0093 if (!m_isDiagonal) {
0094 for (unsigned int j = 1; j < i; ++j) {
0095 double sum(symCorMat(i, j));
0096 for (unsigned int k = 0; k < j; ++k)
0097 sum -= m_H(i, k) * m_H(j, k);
0098 m_H(i, j) = sum / m_H(j, j);
0099 }
0100 }
0101
0102 double hii = symCorMat(i, i);
0103 for (unsigned int j = 0; j < i; ++j) {
0104 const double hij(m_H(i, j));
0105 hii -= hij * hij;
0106 }
0107
0108 if (hii <= k_precision)
0109 return false;
0110
0111 m_H(i, i) = sqrt(hii);
0112 }
0113 }
0114 return true;
0115 }
0116
0117 template <class M>
0118 bool CorrelatedNoisifier<M>::checkDecomposition(const M &symCorMat, M &HHtDiff) const {
0119 bool equal(true);
0120 for (unsigned int i = 0; i < m_H.kRows; ++i) {
0121 for (unsigned int j = 0; j < m_H.kRows; ++j) {
0122 double sum(0);
0123 for (unsigned int k(0); k != m_H.kRows; ++k) {
0124 const double hik(k > i ? 0 : m_H(i, k));
0125 const double htkj(k > j ? 0 : m_H(k, j));
0126 sum += hik * htkj;
0127 }
0128 const double diff(symCorMat(i, j) - sum);
0129 HHtDiff(i, j) = diff;
0130 equal = equal && (fabs(diff) < k_precision);
0131 }
0132 }
0133 return equal;
0134 }
0135
0136 template <class M>
0137 template <class T>
0138 void CorrelatedNoisifier<M>::noisify(T &frame,
0139 CLHEP::HepRandomEngine *engine,
0140 const typename CorrelatedNoisifier<M>::VecDou *rangau) const {
0141 if ((unsigned int)frame.size() != (unsigned int)m_H.kRows)
0142 throw cms::Exception("Configuration") << "CorrelatedNoisifier::noisify(frame) was passed wrong size frame "
0143 "object";
0144
0145 if (nullptr != rangau) {
0146 assert(m_H.kRows == rangau->size());
0147 m_vecgau = *rangau;
0148 } else {
0149 if (m_H.kRows != m_vecgau.size())
0150 m_vecgau = VecDou(m_H.kRows, (double)0.0);
0151 CLHEP::RandGaussQ::shootArray(engine, m_H.kRows, &m_vecgau.front());
0152 }
0153
0154 for (unsigned int i(0); i < m_H.kRows; ++i) {
0155 frame[i] += (m_isIdentity ? m_vecgau[i] : m_H(i, i) * m_vecgau[i]);
0156 if (!m_isDiagonal) {
0157 for (unsigned int j = 0; j < i; ++j)
0158 frame[i] += m_H(j, i) * m_vecgau[j];
0159 }
0160 }
0161 }
0162
0163 template <class M>
0164 const M &CorrelatedNoisifier<M>::cholMat() const {
0165 return m_H;
0166 }
0167
0168 template <class M>
0169 const typename CorrelatedNoisifier<M>::VecDou &CorrelatedNoisifier<M>::vecgau() const {
0170 return m_vecgau;
0171 }