herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
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 
18 using namespace Herwig::Kinematics;
19 using 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  */
33 BOOST_AUTO_TEST_SUITE(utilitiesKinematicsTest)
34 
35 /*
36  * Boost unit tests
37  *
38  */
39 BOOST_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 
51 BOOST_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 
66 BOOST_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 
79 BOOST_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 
100 BOOST_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 
122 BOOST_AUTO_TEST_SUITE_END()
123 
124 #endif /* HERWIG_Utilities_Test_Kinematics_H */
void generateAngles(double &ct, double &az)
This just generates angles.
Definition: Kinematics.h:103
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.
This is a namespace which provides some useful methods for kinematics computation, as the two body decays.
Definition: Kinematics.h:33
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
Axis unitDirection(const double cosTheta, const double phi)
It returns the unit 3-vector with the given cosTheta and phi.
Definition: Kinematics.h:54