herwig is hosted by Hepforge, IPPP Durham
Herwig 7.3.0
utilitiesTestsKinematics.h
1// -*- C++ -*-
2//
3// utilitiesTestKinematics.h is a part of Herwig - A multi-purpose Monte Carlo event generator
4// Copyright (C) 2002-2019 The Herwig Collaboration, 2015 Marco A. Harrendorf
5//
6// Herwig is licenced under version 3 of the GPL, see COPYING for details.
7// Please respect the MCnet academic guidelines, see GUIDELINES for details.
8//
9#ifndef HERWIG_Utilities_Test_Kinematics_H
10#define HERWIG_Utilities_Test_Kinematics_H
11
12#include <boost/test/unit_test.hpp>
13
14#include "Herwig/Utilities/Kinematics.h"
15
16#include "ThePEG/Config/Unitsystem.h"
17
18using namespace Herwig::Kinematics;
19using namespace ThePEG::Units;
20
23 {BOOST_TEST_MESSAGE( "setup fixture for utilitiesKinematicsTestTest" ); }
24
25 ~FixKinematics1() { BOOST_TEST_MESSAGE( "teardown fixture for utilitiesKinematicsTest" ); }
26};
27
28/*
29 * Start of boost unit tests for Kinematics.h
30 *
31 * @todo Implement unit test for threeBodyDecay
32 */
33BOOST_AUTO_TEST_SUITE(utilitiesKinematicsTest)
34
35/*
36 * Boost unit tests
37 *
38 */
39BOOST_AUTO_TEST_CASE(generateAnglesTest)
40{
41 double flatMinusPiToPlusPi, flatNullToTwoPi;
42 for(int i = 0; i < 100; ++i) {
43 generateAngles(flatMinusPiToPlusPi, flatNullToTwoPi);
44 BOOST_CHECK( -M_PI <= flatMinusPiToPlusPi );
45 BOOST_CHECK( flatMinusPiToPlusPi <= M_PI);
46 BOOST_CHECK( 0. <= flatNullToTwoPi);
47 BOOST_CHECK(flatNullToTwoPi <= 2*M_PI);
48 }
49}
50
51BOOST_AUTO_TEST_CASE(unitDirectionTest)
52{
53 using namespace Herwig::Kinematics;
54 BOOST_CHECK_EQUAL( unitDirection(1.1, -1), Axis() );
55 BOOST_CHECK_EQUAL( unitDirection(-1.1, -1), Axis() );
56 BOOST_CHECK_EQUAL( unitDirection(1.1, 0), Axis() );
57 BOOST_CHECK_EQUAL( unitDirection(-1.1, -1), Axis() );
58
59 BOOST_CHECK_EQUAL( unitDirection(1, 0), Axis(0, 0, 1) );
60 BOOST_CHECK_EQUAL( unitDirection(1, M_PI/2.), Axis(0, 0, 1) );
61
62 BOOST_CHECK(unitDirection(0, M_PI/2).almostEqual(Axis(0, 1, 0), 0.001) );
63 BOOST_CHECK_EQUAL( unitDirection(0, 0), Axis(1, 0, 0) );
64}
65
66BOOST_AUTO_TEST_CASE(pstarTwoBodyDecayTest)
67{
68 BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(-100*GeV), Energy(60*GeV), Energy(60*GeV))/GeV, Energy(0*GeV)/GeV);
69 BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(100*GeV), Energy(-40*GeV), Energy(40*GeV))/GeV, Energy(0*GeV)/GeV);
70 BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(100*GeV), Energy(-40*GeV), Energy(-40*GeV))/GeV, Energy(0*GeV)/GeV);
71
72 BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(100*GeV), Energy(60*GeV), Energy(60*GeV))/GeV, Energy(0*GeV)/GeV);
73 BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(100*GeV), Energy(50*GeV), Energy(50*GeV))/GeV, Energy(0*GeV)/GeV);
74
75 BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(10*GeV), Energy(6*GeV), Energy(3*GeV))/GeV, Energy(std::sqrt(19*91)/20.*GeV)/GeV);
76 BOOST_CHECK_EQUAL(pstarTwoBodyDecay(Energy(10*GeV), Energy(3*GeV), Energy(6*GeV))/GeV, Energy(std::sqrt(19*91)/20.*GeV)/GeV);
77}
78
79BOOST_AUTO_TEST_CASE(twoBodyDecayTest1)
80{
81 Lorentz5Momentum decayProductOne(GeV);
82 Lorentz5Momentum decayProductTwo(GeV);
83 BOOST_CHECK(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(40*GeV), Energy(40*GeV), Axis(), decayProductOne, decayProductTwo));
84 BOOST_CHECK(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(50*GeV), Energy(50*GeV), Axis(), decayProductOne, decayProductTwo));
85
86 BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(60*GeV), Energy(60*GeV), Axis(), decayProductOne, decayProductTwo)));
87 BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(-100*GeV), Energy(40*GeV), Energy(40*GeV), Axis(), decayProductOne, decayProductTwo)));
88 BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(-40*GeV), Energy(40*GeV), Axis(), decayProductOne, decayProductTwo)));
89 BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(40*GeV), Energy(-40*GeV), Axis(), decayProductOne, decayProductTwo)));
90
91 twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(50*GeV), Energy(50*GeV), Axis(1,0,0), decayProductOne, decayProductTwo);
92 BOOST_CHECK_EQUAL(decayProductOne/GeV, Lorentz5Momentum(50*GeV)/GeV);
93 BOOST_CHECK_EQUAL(decayProductTwo/GeV, Lorentz5Momentum(50*GeV)/GeV);
94
95 twoBodyDecay(Lorentz5Momentum(10*GeV), Energy(6*GeV), Energy(3*GeV), Axis(1,0,0), decayProductOne, decayProductTwo);
96 BOOST_CHECK_EQUAL(decayProductOne/GeV, Lorentz5Momentum(6*GeV, Momentum3(std::sqrt(19*91)/20.*GeV, ThePEG::ZERO, ThePEG::ZERO))/GeV);
97 BOOST_CHECK_EQUAL(decayProductTwo/GeV, Lorentz5Momentum(3*GeV, Momentum3(-(std::sqrt(19*91)/20.*GeV), ThePEG::ZERO, ThePEG::ZERO))/GeV);
98}
99
100BOOST_AUTO_TEST_CASE(twoBodyDecayTest2)
101{
102 Lorentz5Momentum decayProductOne(GeV);
103 Lorentz5Momentum decayProductTwo(GeV);
104 BOOST_CHECK(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(40*GeV), Energy(40*GeV), 1, M_PI/2., decayProductOne, decayProductTwo));
105 BOOST_CHECK(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(50*GeV), Energy(50*GeV), 1, M_PI/2., decayProductOne, decayProductTwo));
106
107 BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(60*GeV), Energy(60*GeV), 1, M_PI/2., decayProductOne, decayProductTwo)));
108 BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(-100*GeV), Energy(40*GeV), Energy(40*GeV), 1, M_PI/2., decayProductOne, decayProductTwo)));
109 BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(-40*GeV), Energy(40*GeV), 1, M_PI/2., decayProductOne, decayProductTwo)));
110 BOOST_CHECK(!(twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(40*GeV), Energy(-40*GeV), 1, M_PI/2., decayProductOne, decayProductTwo)));
111
112 twoBodyDecay(Lorentz5Momentum(100*GeV), Energy(50*GeV), Energy(50*GeV), 1, M_PI/2., decayProductOne, decayProductTwo);
113 BOOST_CHECK_EQUAL(decayProductOne/GeV, Lorentz5Momentum(50*GeV)/GeV);
114 BOOST_CHECK_EQUAL(decayProductTwo/GeV, Lorentz5Momentum(50*GeV)/GeV);
115
116 twoBodyDecay(Lorentz5Momentum(10*GeV), Energy(6*GeV), Energy(3*GeV), 1, M_PI/2., decayProductOne, decayProductTwo);
117 BOOST_CHECK_EQUAL(decayProductOne/GeV, Lorentz5Momentum(6*GeV, Momentum3(ThePEG::ZERO, ThePEG::ZERO, std::sqrt(19*91)/20.*GeV))/GeV);
118 BOOST_CHECK_EQUAL(decayProductTwo/GeV, Lorentz5Momentum(3*GeV, Momentum3(ThePEG::ZERO, ThePEG::ZERO, -(std::sqrt(19*91)/20.*GeV)))/GeV);
119}
120
121
122BOOST_AUTO_TEST_SUITE_END()
123
124#endif /* HERWIG_Utilities_Test_Kinematics_H */
This is a namespace which provides some useful methods for kinematics computation,...
Definition: Kinematics.h:33
Axis unitDirection(const double cosTheta, const double phi)
It returns the unit 3-vector with the given cosTheta and phi.
Definition: Kinematics.h:54
void generateAngles(double &ct, double &az)
This just generates angles.
Definition: Kinematics.h:103
Energy pstarTwoBodyDecay(const Energy M, const Energy m1, const Energy m2)
For the two body decay M -> m1 + m2 it gives the module of the 3-momentum of the decay product in the...
Definition: Kinematics.h:93
bool twoBodyDecay(const Lorentz5Momentum &p, const Energy m1, const Energy m2, const Axis &unitDir1, Lorentz5Momentum &p1, Lorentz5Momentum &p2)
Calculate the momenta for a two body decay The return value indicates success or failure.
Qty< 0, 1, 0 > Energy
ThreeVector< Energy > Momentum3
Lorentz5Vector< Energy > Lorentz5Momentum
ThreeVector< double > Axis
constexpr ZeroUnit ZERO