Skip to content

Commit b655919

Browse files
greenwaldgreenwald
committed
VectorAlgebra complete
1 parent 4982d20 commit b655919

11 files changed

+312
-195
lines changed

‎.gitignore

+2-1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
build/
22
docs/html
33
docs/latex
4-
4+
*~
5+
\#*\#

‎examples/bat_gen/.gitignore

+7
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
*.root
2+
*.pdf
3+
*.o
4+
*.d
5+
runBatGen
6+
log.txt
7+
logs

‎include/Constants.h

+33-11
Original file line numberDiff line numberDiff line change
@@ -21,16 +21,16 @@
2121
#ifndef yap_Constants_h
2222
#define yap_Constants_h
2323

24+
#include "FourVector.h"
2425
#include "SquareMatrix.h"
2526
#include "ThreeVector.h"
2627

2728
#include <complex>
2829

2930
namespace yap {
3031

31-
/* template <typename T> */
32-
/* typedef FourVector; */
33-
32+
/// \name Complex constants
33+
/// @{
3434
/// complex zero
3535
extern const std::complex<double> Complex_0;
3636

@@ -39,30 +39,52 @@ extern const std::complex<double> Complex_1;
3939

4040
/// complex i
4141
extern const std::complex<double> Complex_i;
42+
/// @}
43+
44+
/// \name real constants
45+
/// @{
4246

4347
/// pi (11 digits)
4448
extern const double PI;
4549

4650
/// convert deg to rad by multiplying by; rad to deg by dividing by
4751
extern const double DEGTORAD;
4852

49-
/// X axis
53+
/// @}
54+
55+
/// \name #ThreeVector constants
56+
/// @{
57+
58+
/// X axis (ThreeVector)
5059
extern const ThreeVector<double> Axis_X;
5160

52-
/// Y axis
61+
/// Y axis (ThreeVector)
5362
extern const ThreeVector<double> Axis_Y;
5463

55-
/// Z axis
64+
/// Z axis (ThreeVector)
5665
extern const ThreeVector<double> Axis_Z;
5766

5867
/// 0 as ThreeVector;
59-
extern const ThreeVector<double> Vect3_0;
68+
extern const ThreeVector<double> ThreeVector_0;
69+
70+
/// @}
71+
72+
/// \name #FourVector constants
73+
/// @{
74+
75+
/// X axis (FourVector)
76+
extern const FourVector<double> FourAxis_X;
77+
78+
/// Y axis (FourVector)
79+
extern const FourVector<double> FourAxis_Y;
80+
81+
/// Z axis (FourVector)
82+
extern const FourVector<double> FourAxis_Z;
6083

61-
/* /// 0 as FourVector; */
62-
/* extern const FourVector<double> Vect4_0; */
84+
/// 0 as FourVector;
85+
extern const FourVector<double> FourVector_0;
6386

64-
/// 3x3 Unit Matrix
65-
extern const SquareMatrix<double, 3> Unit3x3;
87+
/// @}
6688

6789
}
6890

‎include/FourVector.h

+23-26
Original file line numberDiff line numberDiff line change
@@ -21,55 +21,52 @@
2121
#ifndef yap_FourVector_h
2222
#define yap_FourVector_h
2323

24+
#include "SquareMatrix.h"
2425
#include "ThreeVector.h"
25-
#include "ThreeVectorRotation.h"
2626

2727
#include <algorithm>
28-
#include <array>
28+
#include <type_traits>
2929

3030
namespace yap {
3131

32-
33-
/// \class FourVector
32+
/// \typedef FourVector
3433
/// \author Johannes Rauch, Daniel Greenwald
3534
/// \ingroup VectorAlgebra
3635
template <typename T>
37-
class FourVector : public NVector<T, 4>
38-
{
39-
public:
40-
41-
/// Constructor
42-
FourVector(std::initializer_list<T> list) : NVector<T, 4>(list) {}
43-
44-
/// constructor
45-
FourVector(const T& E = 0, const NVector<T, 3> P = Vect3_0) : FourVector( {E, P[0], P[1], P[2]}) {}
36+
using FourVector = typename std::enable_if<std::is_arithmetic<T>::value, NVector<T, 4> >::type;
37+
// using FourVector = NVector<T, 4>;
4638

47-
/// inner (dot) product of FourVectors
48-
T operator*(const NVector<T, 4>& B) const override
49-
{ return (*this)[0] * B[0] - std::inner_product(this->begin() + 1, this->end(), B.begin() + 1, 0); }
39+
/// \return #FourVector
40+
/// \param E 0th component
41+
/// \param P #ThreeVector component
42+
template <typename T>
43+
FourVector<T> fourVector(const T& E, const ThreeVector<T>& P)
44+
{ return FourVector<T>{E, P[0], P[1], P[2]}; }
5045

51-
};
46+
/// \return inner (dot) product for 4-vectors
47+
template <typename T>
48+
T operator*(const FourVector<T>& A, const FourVector<T>& B)
49+
{ return A.front() * B.front() - std::inner_product(A.begin() + 1, A.end(), B.begin() + 1, 0); }
5250

51+
/// unary minus for 4-vector
52+
template <typename T>
53+
FourVector<T> operator-(const FourVector<T>& A)
54+
{ FourVector<T> res = A; std::transform(res.begin() + 1, res.end(), res.begin() + 1, [](const T& t){return -t;}); return res; }
5355

5456
/// \return Spatial #ThreeVector inside #FourVector
5557
template <typename T>
5658
ThreeVector<T> vect(const FourVector<T>& V)
57-
{ return {V[1], V[2], V[3]}; }
59+
{ return ThreeVector<T>{V[1], V[2], V[3]}; }
5860

5961
/// \return boost vector of this #FourVector
6062
template <typename T>
6163
ThreeVector<T> boost(const FourVector<T>& V)
62-
{ return (V[0] != 0) ? (T(1) / V[0]) * vect(V) : Vect3_0; }
63-
64-
/// inner (dot) product of #FourVector's
65-
template <typename T>
66-
T operator*(const FourVector<T>& A, const FourVector<T>& B )
67-
{ return A.front() * B.front() - std::inner_product(A.begin() + 1, A.end(), B.begin() + 1, 0); }
64+
{ return (V[0] != 0) ? (T(1) / V[0]) * vect(V) : ThreeVector<T>{0, 0, 0}; }
6865

6966
/// apply a three-rotation to a FourVector (rotating only the spatial components)
7067
template <typename T>
71-
FourVector<T>& operator*(const ThreeVectorRotation<T>& R, const FourVector<T>& V)
72-
{ return FourVector<T>(V[0], R * vect(V)); }
68+
FourVector<T>& operator*(const ThreeMatrix<T>& R, const FourVector<T>& V)
69+
{ return fourVector<T>(V[0], R * vect(V)); }
7370

7471
}
7572
#endif

‎include/LorentzTransformation.h

+90
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
/* YAP - Yet another PWA toolkit
2+
Copyright 2015, Technische Universitaet Muenchen,
3+
Authors: Daniel Greenwald, Johannes Rauch
4+
5+
This program is free software: you can redistribute it and/or modify
6+
it under the terms of the GNU General Public License as published by
7+
the Free Software Foundation, either version 3 of the License, or
8+
(at your option) any later version.
9+
10+
This program is distributed in the hope that it will be useful,
11+
but WITHOUT ANY WARRANTY; without even the implied warranty of
12+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+
GNU General Public License for more details.
14+
15+
You should have received a copy of the GNU General Public License
16+
along with this program. If not, see <http://www.gnu.org/licenses/>.
17+
*/
18+
19+
/// \file
20+
21+
#ifndef yap_LorentzTransformation_h
22+
#define yap_LorentzTransformation_h
23+
24+
#include "Constants.h"
25+
#include "SquareMatrix.h"
26+
#include "ThreeVector.h"
27+
28+
namespace yap {
29+
30+
/// \return a 4D Lorentz-transformation matrix for a pure rotation
31+
/// \param R #ThreeMatrix defining the rotation
32+
template <typename T>
33+
FourMatrix<T> lorentzTransformation(const ThreeMatrix<T>& R)
34+
{
35+
FourMatrix<T> L = unitMatrix<T, 4>();
36+
for (unsigned i = 0; i < R.size(); ++i)
37+
std::copy(R[i].begin(), R[i].end(), L[i].begin());
38+
return L;
39+
}
40+
41+
/// \return a 4D Lorentz-transformation matrix for a pure boost
42+
/// \param V #FourVector defining boost
43+
template <typename T>
44+
FourMatrix<T> lorentzTransformation(const FourVector<T>& V)
45+
{
46+
FourVector<T> B = (T(1) / V[0]) * V;
47+
T gamma = T(1) / sqrt(B * B);
48+
B[0] = -(gamma + 1);
49+
50+
FourMatrix<T> L = unitMatrix<T, 4>() + outer(B, B) * (gamma / (gamma + 1));
51+
L[0][0] -= gamma * gamma + 1;
52+
return L;
53+
}
54+
55+
/// \return a 4D Lorentz-transformation matrix for a pure boost
56+
/// \param V #ThreeVector defining boost
57+
template <typename T>
58+
FourMatrix<T> lorentzTransformation(const ThreeVector<T>& V)
59+
{ return lorentzTransformation<T>(fourVector<T>(T(1), V)); }
60+
61+
/// \return a 4D Lorentz-transformation matrix for a rotation followed by a boost
62+
/// \param R #ThreeMatrix defining rotation
63+
/// \param V #FourVector defining boost
64+
template <typename T>
65+
FourMatrix<T> lorentzTransformation(const ThreeMatrix<T>& R, const FourVector<T>& V)
66+
{ return lorentzTransformation<T>(V) * lorentzTransformation<T>(R); }
67+
68+
/// \return a 4D Lorentz-transformation matrix for a rotation followed by a boost
69+
/// \param R #ThreeMatrix defining rotation
70+
/// \param V #ThreeVector defining boost
71+
template <typename T>
72+
FourMatrix<T> lorentzTransformation(const ThreeMatrix<T>& R, const ThreeVector<T>& V)
73+
{ return lorentzTransformation<T>(V) * lorentzTransformation<T>(R); }
74+
75+
/// \return a 4D Lorentz-transformation matrix for a boost followed by a rotation
76+
/// \param R #ThreeMatrix defining rotation
77+
/// \param V #FourVector defining boost
78+
template <typename T>
79+
FourMatrix<T> lorentzTransformation(const FourVector<T>& V, const ThreeMatrix<T> R)
80+
{ return lorentzTransformation<T>(R) * lorentzTransformation<T>(V); }
81+
82+
/// \return a 4D Lorentz-transformation matrix for a boost followed by a rotation
83+
/// \param R #ThreeMatrix defining rotation
84+
/// \param V #ThreeVector defining boost
85+
template <typename T>
86+
FourMatrix<T> lorentzTransformation(const ThreeVector<T>& V, const ThreeMatrix<T>& R)
87+
{ return lorentzTransformation<T>(R) * lorentzTransformation<T>(V); }
88+
89+
}
90+
#endif

‎include/NVector.h

+53-48
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,9 @@
2323

2424
#include <algorithm>
2525
#include <array>
26+
#include <cassert>
2627
#include <numeric>
28+
#include <type_traits>
2729

2830
namespace yap {
2931

@@ -35,62 +37,65 @@ template <typename T, size_t N>
3537
class NVector : public std::array<T, N>
3638
{
3739
public:
38-
3940
/// Default Constructor
4041
NVector() : std::array<T, N>() {}
4142

4243
/// Constructor
43-
NVector(std::initializer_list<T> list) : std::array<T, N>()
44-
{ std::copy(list.begin(), list.begin() + N, this->begin()); }
45-
46-
/// Destructor
47-
virtual ~NVector() {}
48-
49-
/// addition assignment
50-
NVector<T, N>& operator+=(const NVector<T, N>& B)
51-
{ std::transform(this->begin(), this->end(), B.begin(), this->begin(), std::operator+); return *this; }
52-
53-
/// subtraction assignment
54-
NVector<T, N>& operator-=(const NVector<T, N>& B)
55-
{ std::transform(this->begin(), this->end(), B.begin(), this->begin(), std::operator-); return *this; }
56-
57-
/// (assignment) multiplication by a single element
58-
NVector<T, N>& operator*=(const T& B)
59-
{ std::transform(this->begin(), this->end(), this->begin(), [&](const T & v) {return B * v;}); return *this; }
60-
61-
/// addition
62-
NVector<T, N> operator+(const NVector<T, N>& B) const
63-
{
64-
NVector<T, N> res;
65-
std::transform(this->begin(), this-> end(), B.begin(), res.begin(), std::operator+);
66-
return res;
67-
}
68-
69-
/// subtraction
70-
NVector<T, N> operator-(const NVector<T, N>& B) const
71-
{
72-
NVector<T, N> res;
73-
std::transform(this->begin(), this->end(), B.begin(), res.begin(), std::operator-);
74-
return res;
75-
}
76-
77-
/// inner (dot) product of #NVector's
78-
virtual T operator*(const NVector<T, N>& B) const
79-
{ return std::inner_product(this->begin(), this->end(), B.begin(), 0); }
80-
81-
/// multiplication: #NVector<T> * T
82-
NVector<T, N> operator*(const T& c) const
83-
{
84-
NVector<T, N> res;
85-
std::transform(this->begin(), this->end(), res.begin(), [&](const T & t) {return t * c;});
86-
return res;
87-
}
44+
NVector(std::initializer_list<T> list) : NVector<T, N>()
45+
{
46+
assert(list.size() == N);
47+
std::copy(list.begin(), list.begin() + N, this->begin());
48+
}
8849
};
8950

51+
/// addition assignment
52+
template <typename T, size_t N>
53+
NVector<T, N>& operator+=(NVector<T, N>& A, const NVector<T, N>& B)
54+
{ std::transform(A.begin(), A.end(), B.begin(), A.begin(), [](const T& a, const T& b){return a + b;}); return A; }
55+
56+
/// subtraction assignment
57+
template <typename T, size_t N>
58+
NVector<T, N>& operator-=(NVector<T, N>& A, const NVector<T, N>& B)
59+
{ std::transform(A.begin(), A.end(), B.begin(), A.begin(), [](const T& a, const T& b){return a - b;}); return A; }
60+
61+
/// (assignment) multiplication by a single element
62+
template <typename T, size_t N>
63+
typename std::enable_if<std::is_arithmetic<T>::value, NVector<T, N>& >::type
64+
operator*=(NVector<T, N>& A, const T& B)
65+
{ std::transform(A.begin(), A.end(), A.begin(), [&](const T & v) {return B * v;}); return A; }
66+
67+
/// addition
68+
template <typename T, size_t N>
69+
NVector<T, N> operator+(const NVector<T, N>& A, const NVector<T, N>& B)
70+
{ NVector<T, N> res = A; res += B; return res; }
71+
72+
/// subtraction
73+
template <typename T, size_t N>
74+
NVector<T, N> operator-(const NVector<T, N>& A, const NVector<T, N>& B)
75+
{ NVector<T, N> res = A; A -= B; return res; }
76+
77+
/// inner (dot) product of #NVector's
78+
template <typename T, size_t N>
79+
typename std::enable_if<std::is_arithmetic<T>::value, T>::type
80+
operator*(const NVector<T, N>& A, const NVector<T, N>& B)
81+
{ return std::inner_product(A.begin(), A.end(), B.begin(), T(0)); }
82+
83+
/// multiplication: #NVector<T> * T
84+
template <typename T, size_t N>
85+
typename std::enable_if<std::is_arithmetic<T>::value, NVector<T, N> >::type
86+
operator*(const NVector<T, N>& A, const T& c)
87+
{ NVector<T, N> res = A; res *= c; return res; }
88+
9089
/// multiplication: T * #NVector<T>
9190
template <typename T, size_t N>
92-
NVector<T, N> operator*(const T& c, const NVector<T, N>& V)
93-
{ return V * c; }
91+
typename std::enable_if<std::is_arithmetic<T>::value, NVector<T, N> >::type
92+
operator*(const T& c, const NVector<T, N>& A)
93+
{ return operator*<T, N>(A, c); }
94+
95+
/// unary minus
96+
template <typename T, size_t N>
97+
NVector<T, N> operator-(const NVector<T, N>& A)
98+
{ NVector<T, N> res = A; std::transform(res.begin(), res.end(), res.begin(), [](const T& t){return -t;}); return res; }
9499

95100
}
96101
#endif

0 commit comments

Comments
 (0)