BinomialTest.cpp 9.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
#include "gtest/gtest.h"

#include "Distribution/Distribution.h"
#include "Attributes/Attributes.h"
#include "Physics/Physics.h"

#include "opal_test_utilities/SilenceTest.h"

#include "gsl/gsl_statistics_double.h"

TEST(BinomialTest, FullSigmaTest1) {
    OpalTestUtilities::SilenceTest silencer(true);

    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 37
    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], "BINOMIAL");
38 39 40
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::MX], 999999999.9);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::MY], 999999999.9);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::MZ], 999999999.9);
41 42 43 44 45 46
    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);
47 48 49 50 51 52 53
    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);
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
    Attributes::setBool(dist.itsAttr[Attrib::Distribution::EMITTED], false);

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

81 82 83 84 85 86 87 88 89 90 91 92 93
    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.01);
    EXPECT_NEAR(R43 / expectedR43, 1.0, 0.01);
    EXPECT_NEAR(R51 / expectedR51, 1.0, 0.01);
    EXPECT_NEAR(R52 / expectedR52, 1.0, 0.01);
    EXPECT_NEAR(R61 / expectedR61, 1.0, 0.01);
    EXPECT_NEAR(R62 / expectedR62, 1.0, 0.01);
94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
}

TEST(BinomialTest, FullSigmaTest2) {
    OpalTestUtilities::SilenceTest silencer(true);

    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;

116 117 118 119
    Attributes::setString(dist.itsAttr[Attrib::Distribution::TYPE], "BINOMIAL");
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::MX], 1.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::MY], 1.0);
    Attributes::setReal(dist.itsAttr[Attrib::Distribution::MZ], 1.0);
120 121 122 123 124 125 126 127 128 129 130 131 132 133
    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);
    Attributes::setBool(dist.itsAttr[Attrib::Distribution::EMITTED], false);
134
    Attributes::setBool(dist.itsAttr[Attrib::Distribution::WRITETOFILE], true);
135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150

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

162 163 164 165 166 167 168 169 170 171 172 173 174
    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.01);
    EXPECT_NEAR(R43 / expectedR43, 1.0, 0.01);
    EXPECT_NEAR(R51 / expectedR51, 1.0, 0.01);
    EXPECT_NEAR(R52 / expectedR52, 1.0, 0.01);
    EXPECT_NEAR(R61 / expectedR61, 1.0, 0.01);
    EXPECT_NEAR(R62 / expectedR62, 1.0, 0.01);
175
}