Skip to content

Commit de7d215

Browse files
ewanwmewanwm
authored
Fix meter to ev^-1 conversion factor (#131)
* fix value of meter -> eV^-1 conversion factor, remove extra factor of 2pi * update test values in the brager test * also fix values in python barger test
1 parent 182fa05 commit de7d215

File tree

5 files changed

+59
-52
lines changed

5 files changed

+59
-52
lines changed

‎nuTens/propagator/propagator.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -81,7 +81,7 @@ class Propagator
8181
_energies = newEnergies;
8282
_weightMatrix = Tensor::ones({_energies.getBatchDim(), _nGenerations, _nGenerations}, dtypes::kComplexFloat)
8383
.requiresGrad(false);
84-
_weightArgDenom = Tensor::scale(Tensor::scale(_energies, 2.0), std::complex<float>(1.0) / (std::complex<float>(-1.0J) * _baseline));
84+
_weightArgDenom = Tensor::scale(Tensor::scale(_energies, 2.0), std::complex<float>(1.0) / (std::complex<float>(-1.0J) * _baseline * 2.0f * (float)M_PI));
8585

8686
if (_matterSolver)
8787
{

‎nuTens/propagator/units.hpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,8 @@
44
/// @brief Defines some suuuper basic units to convert to eV
55

66

7-
// 1 eV^-1 = 0.197 e-6 m
8-
// => 1m = 5.07614213198 e^6 eV^-1
7+
// 1 eV^-1 = 1.239841984 e-6 m
8+
// => 1m = 0.806554394 e^6 eV^-1
99

1010
namespace nuTens
1111
{
@@ -17,7 +17,7 @@ namespace units
1717
static constexpr double MeV = 1e6; // eV
1818
static constexpr double GeV = 1e9; // eV
1919

20-
static constexpr double m = 5.07614213198e6; // eV^-1
20+
static constexpr double m = 0.806554394e6; // eV^-1
2121
static constexpr double cm = 1e-2 * m; // eV^-1
2222
static constexpr double km = 1e3 * m; // eV^-1
2323

‎tests/barger-propagator.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -160,7 +160,7 @@ class TwoFlavourBarger
160160

161161
// now get the actual probabilities
162162
float sin2Gamma = std::sin(2.0 * gamma);
163-
float sinPhi = std::sin(dM2 * _baseline / (4.0 * energy));
163+
float sinPhi = std::sin(dM2 * 2.0 * M_PI * _baseline / (4.0 * energy));
164164

165165
float offAxis = sin2Gamma * sin2Gamma * sinPhi * sinPhi;
166166
float onAxis = 1.0 - offAxis;
@@ -195,6 +195,6 @@ class TwoFlavourBarger
195195
bool _antiNeutrino;
196196
};
197197

198-
} // testing
198+
} // namespace testing
199199

200200
} // namespace nuTens

‎tests/python/test_two-flavour-barger.py‎

Lines changed: 21 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -80,47 +80,49 @@ class TestVacuumOscProbs(unittest.TestCase):
8080
"""test matter propagations for some fixed param values
8181
"""
8282

83-
# theta = 0.24, m1 = 0.04eV, m2 = 0.001eV, E = 1GeV, L = 500km, density = 2
84-
# lv = 4pi * E / dm^2 = 7.8588934e+12
85-
# lm = 2pi / ( sqrt(2) * G * density ) = 2.0588727e+13
86-
# gamma = atan( sin( 2theta ) / (cos( 2theta ) - lv / lm) ) / 2.0
87-
# = atan(0.91389598537 ) / 2 = 0.370219805 rad
88-
# dM2 = dm^2 * sqrt( 1 - 2 * (lv / lm) * cos(2theta) + (lv / lm)^2)
89-
# = 0.00109453
90-
#
91-
# => prob_(alpha != beta) = sin^2(2*gamma) * sin^2((L / E) * dM2/4 )
92-
# = 0.186410
93-
# prob_(alpha == beta) = 1 - 0.186410 = 0.81359
83+
## theta = 0.24, m1 = 0.04eV, m2 = 0.001eV, E = 1GeV, L = 250 km, density = 2
84+
## lv = 4pi * E / dm^2 = 7.8588934e+12
85+
## lm = 2pi / ( sqrt(2) * G * density ) = 4.1177454e+13
86+
## gamma = atan( sin( 2theta ) / (cos( 2theta ) - lv / lm) ) / 2.0
87+
## = atan(0.663342 ) / 2 = 0.292848614 rad
88+
## dM2 = dm^2 * sqrt( 1 - 2 * (lv / lm) * cos(2theta) + (lv / lm)^2)
89+
## = 0.001599 * sqrt ( 1 - 2 * 0.19085428 * 0.8869949 + 0.19085428 ^2)
90+
## = 0.001335765
91+
##
92+
## => prob_(alpha != beta) = sin^2(2*gamma) * sin^2( 1.27 (L / E) * dM2 )
93+
## = sin^2( 2 * 0.292848614 ) * sin^2( 1.27 * 250 / 1 * 0.001599)
94+
## = 0.0517436
95+
## prob_(alpha == beta) = 1 - 0.0517436 = 0.9482564
9496

9597
energy = 1.0 * nt.units.GeV
9698
barger = TwoFlavourBarger()
9799
barger.set_params(m1=0.04 * nt.units.eV, m2=0.001 * nt.units.eV, theta=0.24,
98-
baseline=500.0 * nt.units.km, density=2.0)
100+
baseline=250.0 * nt.units.km, density=2.0)
99101

100102
def test_vacuum_osc_length(self):
101103
self.assertAlmostEqual(self.barger.lv(self.energy), 7.8588934e+12, -6,
102104
f"bad vacuum osc length")
103105

104106
def test_matter_osc_length(self):
105-
self.assertAlmostEqual(self.barger.lm(), 2.0588727e+13, -6,
107+
self.assertAlmostEqual(self.barger.lm(), 4.1177454e+13, -6,
106108
f"bad matter osc length")
107109

108110
def test_effective_mixing_angle(self):
109-
self.assertAlmostEqual(self.barger.calculate_effective_angle(self.energy),0.370219805, 6,
111+
self.assertAlmostEqual(self.barger.calculate_effective_angle(self.energy), 0.292848614, 6,
110112
f"bad effective mixing angle")
111113

112114
def test_effective_dm2(self):
113-
self.assertAlmostEqual(self.barger.calculate_effective_dm2(self.energy),0.00109453, 6,
115+
self.assertAlmostEqual(self.barger.calculate_effective_dm2(self.energy), 0.001335765, 6,
114116
f"bad effective dm^2")
115117

116118
def test_survuval(self):
117-
self.assertAlmostEqual(self.barger.calculate_prob(self.energy, 1, 1), 0.81359, 5,
119+
self.assertAlmostEqual(self.barger.calculate_prob(self.energy, 1, 1), 0.9482564, 3,
118120
f"b->b vacuum prob not as expected")
119-
self.assertAlmostEqual(self.barger.calculate_prob(self.energy, 0, 0), 0.81359, 5,
121+
self.assertAlmostEqual(self.barger.calculate_prob(self.energy, 0, 0), 0.9482564, 3,
120122
f"a->a vacuum prob not as expected")
121123

122124
def test_oscillation(self):
123-
self.assertAlmostEqual(self.barger.calculate_prob(self.energy, 0, 1), 0.186410, 5,
125+
self.assertAlmostEqual(self.barger.calculate_prob(self.energy, 0, 1), 0.0517436, 3,
124126
f"a->b vacuum prob not as expected")
125-
self.assertAlmostEqual(self.barger.calculate_prob(self.energy, 1, 0), 0.186410, 5,
127+
self.assertAlmostEqual(self.barger.calculate_prob(self.energy, 1, 0), 0.0517436, 3,
126128
f"b->a vacuum prob not as expected")

‎tests/test-barger.cpp‎

Lines changed: 32 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ using namespace nuTens::testing;
1111

1212
TEST(TwoFlavourBargerPropTest, zeroThetaNoOscTest) {
1313

14-
constexpr float baseline = 500.0 * units::km;
14+
constexpr float baseline = 5.0e12;
1515

1616
TwoFlavourBarger bargerProp{};
1717

@@ -36,7 +36,7 @@ TEST(TwoFlavourBargerPropTest, zeroThetaNoOscTest) {
3636

3737
TEST(TwoFlavourBargerPropTest, zeroDmsqNoOscTest) {
3838

39-
constexpr float baseline = 500.0 * units::km;
39+
constexpr float baseline = 5.0e12;
4040

4141
TwoFlavourBarger bargerProp{};
4242

@@ -65,54 +65,59 @@ TEST(TwoFlavourBargerPropTest, fixedValuesTest) {
6565

6666
// now check for fixed parameters values against externally calculated values
6767

68-
// theta = pi/8, m1 = 1, m2 = 2, E = 3, L = 4
69-
// => prob_(alpha != beta) = sin^2(Pi/4) * sin^2(1) = 0.35403670913
70-
// prob_(alpha == beta) = 1 - 0.35403670913 = 0.64596329086
68+
// theta = pi/8, dm^2 = 0.01 eV E = 1 GeV L = 100 km
69+
// => prob_(alpha != beta) = sin^2(2 theta) * sin^2( 1.27 * dm^2 * L / E [ eV^2 km / GeV] )
70+
//
71+
// = sin^2(Pi/4) * sin^2( 1.27 * 0.01 * 100 / 1 ) = 0.4561088222
72+
//
73+
// prob_(alpha == beta) = 1 - 0.4561088222 = 0.5438911778
7174

72-
bargerProp.setParams(/*m1=*/1.0, /*m2=*/2.0, /*theta=*/M_PI / 8.0,
73-
/*baseline=*/4.0);
75+
bargerProp.setParams(/*m1=*/0.0, /*m2=*/0.1, /*theta=*/M_PI / 8.0,
76+
/*baseline=*/100.0 * units::km );
7477

75-
ASSERT_NEAR(bargerProp.calculateProb(3.0, 0, 0), 0.64596329086, 1e-5);
78+
ASSERT_NEAR(bargerProp.calculateProb(1.0 * units::GeV, 0, 0), 0.5438911778, 1e-3);
7679

77-
ASSERT_NEAR(bargerProp.calculateProb(3.0, 1, 1), 0.64596329086, 1e-5);
80+
ASSERT_NEAR(bargerProp.calculateProb(1.0 * units::GeV, 1, 1), 0.5438911778, 1e-3);
7881

79-
ASSERT_NEAR(bargerProp.calculateProb(3.0, 0, 1), 0.35403670913, 1e-5);
82+
ASSERT_NEAR(bargerProp.calculateProb(1.0 * units::GeV, 0, 1), 0.4561088222, 1e-3);
8083

81-
ASSERT_NEAR(bargerProp.calculateProb(3.0, 1, 0), 0.35403670913, 1e-5);
84+
ASSERT_NEAR(bargerProp.calculateProb(1.0 * units::GeV, 1, 0), 0.4561088222, 1e-3);
8285

8386

8487
// ##############################################################
8588
// ## Now test matter propagations for some fixed param values ##
8689
// ##############################################################
8790

88-
// theta = 0.24, m1 = 0.04eV, m2 = 0.001eV, E = 1GeV, L = 500km, density = 2
91+
// theta = 0.24, m1 = 0.04eV, m2 = 0.001eV, E = 1GeV, L = 250 km, density = 2
8992
// lv = 4pi * E / dm^2 = 7.8588934e+12
90-
// lm = 2pi / ( sqrt(2) * G * density ) = 2.0588727e+13
93+
// lm = 2pi / ( sqrt(2) * G * density ) = 4.1177454e+13
9194
// gamma = atan( sin( 2theta ) / (cos( 2theta ) - lv / lm) ) / 2.0
92-
// = atan(0.91389598537 ) / 2 = 0.370219805 rad
95+
// = atan(0.663342 ) / 2 = 0.292848614 rad
9396
// dM2 = dm^2 * sqrt( 1 - 2 * (lv / lm) * cos(2theta) + (lv / lm)^2)
94-
// = 0.00109453
97+
// = 0.001599 * sqrt ( 1 - 2 * 0.19085428 * 0.8869949 + 0.19085428 ^2)
98+
// = 0.001335765
9599
//
96-
// => prob_(alpha != beta) = sin^2(2*gamma) * sin^2((L / E) * dM2/4 )
97-
// = 0.186410
98-
// prob_(alpha == beta) = 1 - 0.186410 = 0.81359
100+
// => prob_(alpha != beta) = sin^2(2*gamma) * sin^2( 1.27 (L / E) * dM2 )
101+
// = sin^2( 2 * 0.292848614 ) * sin^2( 1.27 * 250 / 1 * 0.001599)
102+
// = 0.0517436
103+
// prob_(alpha == beta) = 1 - 0.0517436 = 0.9482564
99104

100105
bargerProp.setParams(/*m1=*/0.04, /*m2=*/0.001, /*theta=*/0.24,
101-
/*baseline=*/500.0 * units::km, /*density=*/2.0);
106+
/*baseline=*/250 * units::km, /*density=*/2.0);
102107

103-
ASSERT_NEAR(bargerProp.lv(1.0 * units::GeV), 7.8588934e+12, 1e6) << "vacuum osc length";
108+
ASSERT_NEAR(bargerProp.lv(1.0e9), 7.8588934e+12, 1e6) << "vacuum osc length";
104109

105-
ASSERT_NEAR(bargerProp.lm(), 2.0588727e+13, 1e6) << "matter osc length";
110+
ASSERT_NEAR(bargerProp.lm(), 4.1177454e+13 , 1e6) << "matter osc length";
106111

107-
ASSERT_NEAR(bargerProp.calculateEffectiveAngle(1.0 * units::GeV), 0.370219805, 0.00001) << "effective mixing angle";
112+
ASSERT_NEAR(bargerProp.calculateEffectiveAngle(1.0e9), 0.292848614, 0.00001) << "effective mixing angle";
108113

109-
ASSERT_NEAR(bargerProp.calculateEffectiveDm2(1.0 * units::GeV), 0.00109453, 0.00001) << "effective m^2 diff";
114+
ASSERT_NEAR(bargerProp.calculateEffectiveDm2(1.0e9), 0.001335765, 0.00001) << "effective m^2 diff";
110115

111-
ASSERT_NEAR(bargerProp.calculateProb(1.0 * units::GeV, 0, 0), 0.81359, 0.00001) << "probability for alpha == beta == 0";
116+
ASSERT_NEAR(bargerProp.calculateProb(1.0e9, 0, 0), 0.9482564, 1e-3) << "probability for alpha == beta == 0";
112117

113-
ASSERT_NEAR(bargerProp.calculateProb(1.0 * units::GeV, 1, 1), 0.81359, 0.00001) << "probability for alpha == beta == 1";
118+
ASSERT_NEAR(bargerProp.calculateProb(1.0e9, 1, 1), 0.9482564, 1e-3) << "probability for alpha == beta == 1";
114119

115-
ASSERT_NEAR(bargerProp.calculateProb(1.0 * units::GeV, 0, 1), 0.186410, 0.00001) << "probability for alpha == 0, beta == 1";
120+
ASSERT_NEAR(bargerProp.calculateProb(1.0e9, 0, 1), 0.0517436, 1e-3) << "probability for alpha == 0, beta == 1";
116121

117-
ASSERT_NEAR(bargerProp.calculateProb(1.0 * units::GeV, 1, 0), 0.186410, 0.00001) << "probability for alpha == 1, beta == 0";
122+
ASSERT_NEAR(bargerProp.calculateProb(1.0e9, 1, 0), 0.0517436, 1e-3) << "probability for alpha == 1, beta == 0";
118123
}

0 commit comments

Comments
 (0)