GaussTest.cpp 9.38 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;
13 14 15 16 17

    const double expectedR11 = 1.978;
    const double expectedR22 = 0.7998;
    const double expectedR33 = 2.498;
    const double expectedR44 = 0.6212;
18
    const double expectedR55 = 1.537;
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
    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;

37
    Attributes::setUpperCaseString(dist.itsAttr[Attrib::Distribution::TYPE], "GAUSS");
38 39 40 41 42 43
    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

    dist.setDistType();
    dist.checkIfEmitted();
    size_t numParticles = 1000000;
57
    dist.totalNumberParticles_m = numParticles;
frey_m's avatar
frey_m committed
58
    dist.create(numParticles, Physics::m_p, Physics::z_p);
59 60 61 62 63 64 65 66 67 68

    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));
69 70
    double R43 = (gsl_stats_covariance(&(dist.yDist_m[0]), 1, &(dist.pyDist_m[0]), 1, dist.yDist_m.size()) * 1e3 /
                  (expectedR33 * expectedR44));
71 72 73 74 75 76 77 78 79
    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));

80 81 82 83 84 85 86 87 88 89 90 91 92
    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);
93 94 95
}

TEST(GaussTest, FullSigmaTest2) {
96
    OpalTestUtilities::SilenceTest silencer;
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114

    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;

115
    Attributes::setUpperCaseString(dist.itsAttr[Attrib::Distribution::TYPE], "GAUSS");
116 117 118 119 120 121 122 123 124 125 126 127 128
    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);
129 130 131 132 133 134
    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);
135
    Attributes::setBool(dist.itsAttr[Attrib::Distribution::EMITTED], false);
136
    Attributes::setBool(dist.itsAttr[Attrib::Distribution::WRITETOFILE], true);
137 138 139 140 141

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

    size_t numParticles = 1000000;
142
    dist.totalNumberParticles_m = numParticles;
frey_m's avatar
frey_m committed
143
    dist.create(numParticles, Physics::m_p, Physics::z_p);
144 145 146 147 148 149 150 151 152 153

    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));
154 155
    double R43 = (gsl_stats_covariance(&(dist.yDist_m[0]), 1, &(dist.pyDist_m[0]), 1, dist.yDist_m.size()) * 1e3 /
                  (expectedR33 * expectedR44));
156 157 158 159 160 161 162 163 164
    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));

165 166 167 168 169 170 171 172 173 174 175 176 177
    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);
178
}