File indexing completed on 2024-04-06 12:30:39
0001 #include "CLHEP/Random/JamesRandom.h"
0002 #include "CommonTools/Statistics/interface/AutocorrelationAnalyzer.h"
0003 #include "SimGeneral/NoiseGenerators/interface/CorrelatedNoisifier.h"
0004 #include "SimGeneral/NoiseGenerators/interface/CorrelatedNoisifier.icc"
0005
0006 namespace CLHEP {
0007 class HepRandomEngine;
0008 }
0009
0010 typedef math::ErrorD<10>::type MyMat;
0011
0012 template class CorrelatedNoisifier<MyMat>;
0013
0014 template void CorrelatedNoisifier<MyMat>::noisify(std::vector<double> &,
0015 CLHEP::HepRandomEngine *,
0016 const std::vector<double> *) const;
0017
0018 int main() {
0019 math::ErrorD<10>::type input;
0020 for (int k = 0; k < 10; k++) {
0021 for (int kk = k; kk < 10; kk++) {
0022 input(k, kk) = kk == k ? 1
0023 : kk == k + 1 ? 0.67
0024 : kk == k + 2 ? 0.53
0025 : kk == k + 3 ? 0.44
0026 : kk == k + 4 ? 0.39
0027 : kk == k + 5 ? 0.36
0028 : kk == k + 6 ? 0.38
0029 : kk == k + 7 ? 0.35
0030 : kk == k + 8 ? 0.36
0031 : kk == k + 9 ? 0.32
0032 : 0.;
0033 }
0034 }
0035 CLHEP::HepJamesRandom engine;
0036
0037 typedef math::ErrorD<10>::type MatType;
0038 typedef CorrelatedNoisifier<MatType> Noisifier;
0039
0040 MatType chol;
0041
0042 for (unsigned int itry(0); itry != 2; ++itry) {
0043
0044 Noisifier noisifier(0 == itry ? Noisifier(input) : Noisifier(nullptr, chol));
0045
0046 if (0 == itry)
0047 chol = noisifier.cholMat();
0048
0049 AutocorrelationAnalyzer analyzer(10);
0050
0051 const unsigned int nTotal = 1000000;
0052 for (unsigned int i = 0; i < nTotal; ++i) {
0053 std::vector<double> samples(10);
0054 noisifier.noisify(samples, &engine);
0055 analyzer.analyze(samples);
0056 }
0057
0058 double big(-999);
0059 math::ErrorD<10>::type ratdif;
0060 for (int k = 0; k < 10; k++) {
0061 for (int kk = k; kk < 10; kk++) {
0062 const double diff(input(k, kk) - analyzer.correlation(k, kk));
0063 ratdif(k, kk) = diff / input(k, kk);
0064 if (fabs(ratdif(k, kk)) > big)
0065 big = fabs(ratdif(k, kk));
0066 }
0067 }
0068
0069 std::cout << "In " << nTotal << " trials, the biggest fractional deviation\n"
0070 << "of observed correlations from input correlations =" << big << std::endl;
0071
0072 std::cout << std::endl
0073 << "Initial correlations:"
0074 << "\n"
0075 << input << std::endl;
0076 ;
0077
0078 std::cout << analyzer << std::endl;
0079
0080 std::cout << ratdif << std::endl;
0081 std::cout << std::endl << "\nSQUARE of initial matrix:\n" << input * input << std::endl;
0082 }
0083 }
