GaussTest.cpp 9.26 KB
Newer Older
1 2 3
#include "gtest/gtest.h"

#include "Distribution/Distribution.h"
4
#include "Attributes/Attributes.h"
5 6
#include "Physics/Physics.h"

7 8
#include "opal_test_utilities/SilenceTest.h"

9 10
#include "gsl/gsl_statistics_double.h"

11
TEST(GaussTest, FullSigmaTest1) {
12
    OpalTestUtilities::SilenceTest silencer(true);
13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

    const double expectedR11 = 1.978;
    const double expectedR22 = 0.7998;
    const double expectedR33 = 2.498;
    const double expectedR44 = 0.6212;
    const double expectedR55 =  1.537;
    const double expectedR66 = 0.9457;

    const double expectedR21 = -0.40993;
    const double expectedR43 = 0.77208;
    const double expectedR65 = 0.12051;
    const double expectedR51 = 0.14935;
    const double expectedR52 = 0.59095;
    const double expectedR61 = 0.72795;
    const double expectedR62 = -0.3550;

    std::vector<double> expectedR({expectedR21, 0, 0,           expectedR51, expectedR61, \
                /*                           */ 0, 0,           expectedR52, expectedR62, \
                /*                              */ expectedR43, 0,           0, \
                /*                                            */0,           0,                                         \
                /*                                                         */expectedR65});

    Distribution dist;

    Attributes::setString(dist.itsAttr[Attrib::Distribution::TYPE], "GAUSS");
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAX], expectedR11 * 1e-3);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAPX], expectedR22);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAY], expectedR33 * 1e-3);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAPY], expectedR44);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAZ], expectedR55 * 1e-3);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAPZ], expectedR66);
44 45 46 47 48 49
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFX], 5.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFY], 5.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFLONG], 5.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFPX], 5.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFPY], 5.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFPZ], 5.0);
50 51
    Attributes::setRealArray(dist.itsAttr[Attrib::Distribution::R], expectedR);
    Attributes::setBool(dist.itsAttr[Attrib::Distribution::EMITTED], false);
52
    Attributes::setBool(dist.itsAttr[Attrib::Distribution::WRITETOFILE], true);
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67

    dist.setDistType();
    dist.checkIfEmitted();
    size_t numParticles = 1000000;
    dist.create(numParticles, Physics::m_p);

    double R11 = sqrt(gsl_stats_variance(&(dist.xDist_m[0]), 1, dist.xDist_m.size())) * 1e3;
    double R22 = sqrt(gsl_stats_variance(&(dist.pxDist_m[0]), 1, dist.pxDist_m.size()));
    double R33 = sqrt(gsl_stats_variance(&(dist.yDist_m[0]), 1, dist.yDist_m.size())) * 1e3;
    double R44 = sqrt(gsl_stats_variance(&(dist.pyDist_m[0]), 1, dist.pyDist_m.size()));
    double R55 = sqrt(gsl_stats_variance(&(dist.tOrZDist_m[0]), 1, dist.tOrZDist_m.size())) * 1e3;
    double R66 = sqrt(gsl_stats_variance(&(dist.pzDist_m[0]), 1, dist.pzDist_m.size()));

    double R21 = (gsl_stats_covariance(&(dist.xDist_m[0]), 1, &(dist.pxDist_m[0]), 1, dist.xDist_m.size()) * 1e3 /
                  (expectedR11 * expectedR22));
68 69
    double R43 = (gsl_stats_covariance(&(dist.yDist_m[0]), 1, &(dist.pyDist_m[0]), 1, dist.yDist_m.size()) * 1e3 /
                  (expectedR33 * expectedR44));
70 71 72 73 74 75 76 77 78
    double R51 = (gsl_stats_covariance(&(dist.xDist_m[0]), 1, &(dist.tOrZDist_m[0]), 1, dist.xDist_m.size()) * 1e6 /
                  (expectedR11 * expectedR55));
    double R52 = (gsl_stats_covariance(&(dist.pxDist_m[0]), 1, &(dist.tOrZDist_m[0]), 1, dist.pxDist_m.size()) * 1e3 /
                  (expectedR22 * expectedR55));
    double R61 = (gsl_stats_covariance(&(dist.xDist_m[0]), 1, &(dist.pzDist_m[0]), 1, dist.xDist_m.size()) * 1e3 /
                  (expectedR11 * expectedR66));
    double R62 = (gsl_stats_covariance(&(dist.pxDist_m[0]), 1, &(dist.pzDist_m[0]), 1, dist.pxDist_m.size()) /
                  (expectedR22 * expectedR66));

79 80 81 82 83 84 85 86 87 88 89 90 91
    EXPECT_NEAR(R11 / expectedR11, 1.0, 0.002);
    EXPECT_NEAR(R22 / expectedR22, 1.0, 0.002);
    EXPECT_NEAR(R33 / expectedR33, 1.0, 0.002);
    EXPECT_NEAR(R44 / expectedR44, 1.0, 0.002);
    EXPECT_NEAR(R55 / expectedR55, 1.0, 0.002);
    EXPECT_NEAR(R66 / expectedR66, 1.0, 0.002);

    EXPECT_NEAR(R21 / expectedR21, 1.0, 0.02);
    EXPECT_NEAR(R43 / expectedR43, 1.0, 0.02);
    EXPECT_NEAR(R51 / expectedR51, 1.0, 0.05);
    EXPECT_NEAR(R52 / expectedR52, 1.0, 0.02);
    EXPECT_NEAR(R61 / expectedR61, 1.0, 0.02);
    EXPECT_NEAR(R62 / expectedR62, 1.0, 0.02);
92 93 94
}

TEST(GaussTest, FullSigmaTest2) {
95
    OpalTestUtilities::SilenceTest silencer(true);
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127

    const double expectedR11 = 1.978;
    const double expectedR22 = 0.7998;
    const double expectedR33 = 2.498;
    const double expectedR44 = 0.6212;
    const double expectedR55 =  1.537;
    const double expectedR66 = 0.9457;

    const double expectedR21 = -0.40993;
    const double expectedR43 = 0.77208;
    const double expectedR65 = 0.12051;
    const double expectedR51 = 0.14935;
    const double expectedR52 = 0.59095;
    const double expectedR61 = 0.72795;
    const double expectedR62 = -0.3550;

    Distribution dist;

    Attributes::setString(dist.itsAttr[Attrib::Distribution::TYPE], "GAUSS");
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAX], expectedR11 * 1e-3);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAPX], expectedR22);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAY], expectedR33 * 1e-3);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAPY], expectedR44);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAZ], expectedR55 * 1e-3);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::SIGMAPZ], expectedR66);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CORRX], expectedR21);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CORRY], expectedR43);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CORRZ], expectedR65);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::R51], expectedR51);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::R61], expectedR61);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::R52], expectedR52);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::R62], expectedR62);
128 129 130 131 132 133
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFX], 5.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFY], 5.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFLONG], 5.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFPX], 5.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFPY], 5.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::CUTOFFPZ], 5.0);
134
    Attributes::setBool(dist.itsAttr[Attrib::Distribution::EMITTED], false);
135
    Attributes::setBool(dist.itsAttr[Attrib::Distribution::WRITETOFILE], true);
136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151

    dist.setDistType();
    dist.checkIfEmitted();

    size_t numParticles = 1000000;
    dist.create(numParticles, Physics::m_p);

    double R11 = sqrt(gsl_stats_variance(&(dist.xDist_m[0]), 1, dist.xDist_m.size())) * 1e3;
    double R22 = sqrt(gsl_stats_variance(&(dist.pxDist_m[0]), 1, dist.pxDist_m.size()));
    double R33 = sqrt(gsl_stats_variance(&(dist.yDist_m[0]), 1, dist.yDist_m.size())) * 1e3;
    double R44 = sqrt(gsl_stats_variance(&(dist.pyDist_m[0]), 1, dist.pyDist_m.size()));
    double R55 = sqrt(gsl_stats_variance(&(dist.tOrZDist_m[0]), 1, dist.tOrZDist_m.size())) * 1e3;
    double R66 = sqrt(gsl_stats_variance(&(dist.pzDist_m[0]), 1, dist.pzDist_m.size()));

    double R21 = (gsl_stats_covariance(&(dist.xDist_m[0]), 1, &(dist.pxDist_m[0]), 1, dist.xDist_m.size()) * 1e3 /
                  (expectedR11 * expectedR22));
152 153
    double R43 = (gsl_stats_covariance(&(dist.yDist_m[0]), 1, &(dist.pyDist_m[0]), 1, dist.yDist_m.size()) * 1e3 /
                  (expectedR33 * expectedR44));
154 155 156 157 158 159 160 161 162
    double R51 = (gsl_stats_covariance(&(dist.xDist_m[0]), 1, &(dist.tOrZDist_m[0]), 1, dist.xDist_m.size()) * 1e6 /
                  (expectedR11 * expectedR55));
    double R52 = (gsl_stats_covariance(&(dist.pxDist_m[0]), 1, &(dist.tOrZDist_m[0]), 1, dist.pxDist_m.size()) * 1e3 /
                  (expectedR22 * expectedR55));
    double R61 = (gsl_stats_covariance(&(dist.xDist_m[0]), 1, &(dist.pzDist_m[0]), 1, dist.xDist_m.size()) * 1e3 /
                  (expectedR11 * expectedR66));
    double R62 = (gsl_stats_covariance(&(dist.pxDist_m[0]), 1, &(dist.pzDist_m[0]), 1, dist.pxDist_m.size()) /
                  (expectedR22 * expectedR66));

163 164 165 166 167 168 169 170 171 172 173 174 175
    EXPECT_NEAR(R11 / expectedR11, 1.0, 0.002);
    EXPECT_NEAR(R22 / expectedR22, 1.0, 0.002);
    EXPECT_NEAR(R33 / expectedR33, 1.0, 0.002);
    EXPECT_NEAR(R44 / expectedR44, 1.0, 0.002);
    EXPECT_NEAR(R55 / expectedR55, 1.0, 0.002);
    EXPECT_NEAR(R66 / expectedR66, 1.0, 0.002);

    EXPECT_NEAR(R21 / expectedR21, 1.0, 0.02);
    EXPECT_NEAR(R43 / expectedR43, 1.0, 0.02);
    EXPECT_NEAR(R51 / expectedR51, 1.0, 0.05);
    EXPECT_NEAR(R52 / expectedR52, 1.0, 0.02);
    EXPECT_NEAR(R61 / expectedR61, 1.0, 0.02);
    EXPECT_NEAR(R62 / expectedR62, 1.0, 0.02);
176
}