Cleaned up vector algebra. implemented ThreeVector tests. Using unique

pointers in ParticleFactory

Author: greenwald

include/Constants.h

@@ -62,7 +62,7 @@ constexpr auto ThreeAxis_X = ThreeVector<double>({1, 0, 0});
 constexpr auto ThreeAxis_Y = ThreeVector<double>({0, 1, 0});
 
 /// Z axis (ThreeVector)
-constexpr auto ThreeAxis_Z = ThreeVector<double>({0, 1, 0});
+constexpr auto ThreeAxis_Z = ThreeVector<double>({0, 0, 1});
 
 /// Standard 3D coordinate system
 constexpr auto ThreeAxes = CoordinateSystem<double, 3> {ThreeAxis_X, ThreeAxis_Y, ThreeAxis_Z};

include/FourVector.h

@@ -39,24 +39,18 @@ template <typename T>
 class FourVector : public Vector<T, 4>
 {
 public:
-    /// Default constructor
-    constexpr FourVector() noexcept : Vector<T, 4> () {}
-
     /// intializer_list constructor
     constexpr FourVector(const std::array<T, 4>& l) noexcept : Vector<T, 4>(l) {}
 
-    // /// initializer with individual elements
-    // constexpr FourVector(const T& E, const T& P0, const T& P1, const T& P2) : FourVector({E, P0, P1, P2}) {}
-
     /// energy + ThreeVector constructor
     /// \param E 0th component
     /// \param P #ThreeVector component
     constexpr FourVector(const T& E, const ThreeVector<T>& P) noexcept : Vector<T, 4>( {E, P[0], P[1], P[2]}) {}
 
-    // /// Vector<T, 4> constructor
-    // constexpr FourVector(const Vector<T, 4>&& V) : Vector<T, 4>(V) {}
+    /// Default constructor
+    FourVector() = default;
 
-    /// using assignment with rhs = brace-enclosed list
+    /// Use std::array's assignment operators
     using Vector<T, 4>::operator=;
 
     /// \return inner (dot) product for 4-vectors

include/ParticleFactory.h

@@ -65,21 +65,21 @@ class ParticleFactory
 
     /// Create a FinalStateParticle from a PDG code
     /// \param PDG PDG code of particle to create
-    /// \return shared pointer to new final state particle
-    std::shared_ptr<FinalStateParticle> createFinalStateParticle(int PDG);
+    /// \return unique pointer to new final state particle
+    std::unique_ptr<FinalStateParticle> createFinalStateParticle(int PDG);
 
     /// Create an InitialStateParticle from a PDG code and a MassShape
     /// \param PDG PDG code of particle to create
     /// \param radialSize radial size of particle to create [GeV^-1]
-    /// \return pointer to new InitialStateParticle object
-    std::shared_ptr<InitialStateParticle> createInitialStateParticle(int PDG, double radialSize);
+    /// \return unique pointer to new InitialStateParticle object
+    std::unique_ptr<InitialStateParticle> createInitialStateParticle(int PDG, double radialSize);
 
     /// Create a Resonance from a PDG code and a MassShape
     /// \param PDG PDG code of particle to create
     /// \param radialSize Radial size of particle to create [GeV^-1]
     /// \param massShape Pointer to MassShape object describing resonance
-    /// \return pointer to new Resonance object
-    std::shared_ptr<Resonance> createResonance(int PDG, double radialSize, std::unique_ptr<MassShape>&& massShape);
+    /// \return unique pointer to new Resonance object
+    std::unique_ptr<Resonance> createResonance(int PDG, double radialSize, std::unique_ptr<MassShape>&& massShape);
 
     /// \name Particle table access
     /// @{

include/Vector.h

@@ -29,6 +29,7 @@
 #include <cassert>
 #include <cmath>
 #include <numeric>
+#include <ostream>
 #include <stdexcept>
 #include <string>
 
@@ -42,19 +43,19 @@ template <typename T, size_t N>
 class Vector : public std::array<T, N>
 {
 public:
-    /// Default Constructor
-    constexpr Vector() noexcept : std::array<T, N>() {}
-
     /// Constructor
-    constexpr Vector(const std::array<T, N>& list) noexcept : std::array<T, N>(list) {}
+    constexpr Vector(const std::array<T, N>& v) noexcept : std::array<T, N>(v) {}
+
+    /// Default constructor
+    Vector() = default;
+
+    /// Use std::array's assignment operators
+    using std::array<T, N>::operator=;
 
     /// inner (dot) product of #Vector's
     virtual constexpr T operator*(const Vector<T, N>& B) const
     { return std::inner_product(this->begin(), this->end(), B.begin(), T(0)); }
 
-    /// using assignment with rhs = brace-enclosed list
-    using std::array<T, N>::operator=;
-
     /// unary minus
     constexpr virtual Vector<T, N> operator-() const
     { return T(-1) * *(this); }
@@ -73,6 +74,10 @@ std::string to_string(const Vector<T, N>& V)
     return s;
 }
 
+template <typename T, size_t N>
+std::ostream& operator<<(std::ostream& os, const Vector<T, N>& V)
+{ os << to_string(V); return os; }
+
 /// addition assignment
 template <typename T, size_t N>
 Vector<T, N>& operator+=(Vector<T, N>& A, const Vector<T, N>& B)

programs/D3piTest.cxx

@@ -24,42 +24,42 @@ int main( int argc, char** argv)
     yap::ParticleFactory factory((::getenv("YAPDIR") ? (std::string)::getenv("YAPDIR") : ".") + "/evt.pdl");
 
     // initial state particle
-    auto D = factory.createInitialStateParticle(factory.pdgCode("D+"), radialSize);
+    std::shared_ptr<yap::InitialStateParticle> D = factory.createInitialStateParticle(factory.pdgCode("D+"), radialSize);
 
     // final state particles
-    auto piPlus = factory.createFinalStateParticle(211);
-    auto piMinus = factory.createFinalStateParticle(-211);
+    std::shared_ptr<yap::FinalStateParticle> piPlus = factory.createFinalStateParticle(211);
+    std::shared_ptr<yap::FinalStateParticle> piMinus = factory.createFinalStateParticle(-211);
 
     // set final state
     D->setFinalStateParticles({piPlus, piMinus, piPlus});
 
     // rho
-    auto rho = std::make_shared<yap::Resonance>(factory.quantumNumbers("rho0"), 0.775, "rho", radialSize, std::make_unique<yap::BreitWigner>());
+    std::shared_ptr<yap::Resonance> rho = std::make_shared<yap::Resonance>(factory.quantumNumbers("rho0"), 0.775, "rho", radialSize, std::make_unique<yap::BreitWigner>());
     static_cast<yap::BreitWigner&>(rho->massShape()).width()->setValue(0.149);
     rho->addChannels(piPlus, piMinus, max2L);
 
     // f_2(1270)
-    auto f_2 = std::make_shared<yap::Resonance>(factory.quantumNumbers("f_2"), 1.275, "f_2", radialSize, std::make_unique<yap::BreitWigner>());
+    std::shared_ptr<yap::Resonance> f_2 = std::make_shared<yap::Resonance>(factory.quantumNumbers("f_2"), 1.275, "f_2", radialSize, std::make_unique<yap::BreitWigner>());
     static_cast<yap::BreitWigner&>(f_2->massShape()).width()->setValue(0.185);
     f_2->addChannels(piPlus, piMinus, max2L);
 
     // f_0(980)
-    auto f_0_980 = std::make_shared<yap::Resonance>(factory.quantumNumbers("f_0"), 0.980, "f_0_980", radialSize, std::make_unique<yap::BreitWigner>());
+    std::shared_ptr<yap::Resonance> f_0_980 = std::make_shared<yap::Resonance>(factory.quantumNumbers("f_0"), 0.980, "f_0_980", radialSize, std::make_unique<yap::BreitWigner>());
     static_cast<yap::BreitWigner&>(f_0_980->massShape()).width()->setValue(0.329);
     f_0_980->addChannels(piPlus, piMinus, max2L);
 
     // f_0(1370)
-    auto f_0_1370 = std::make_shared<yap::Resonance>(factory.quantumNumbers("f_0"), 1.350, "f_0_1370", radialSize, std::make_unique<yap::BreitWigner>());
+    std::shared_ptr<yap::Resonance> f_0_1370 = std::make_shared<yap::Resonance>(factory.quantumNumbers("f_0"), 1.350, "f_0_1370", radialSize, std::make_unique<yap::BreitWigner>());
     static_cast<yap::BreitWigner&>(f_0_1370->massShape()).width()->setValue(0.250);
     f_0_1370->addChannels(piPlus, piMinus, max2L);
 
     // f_0(1500)
-    auto f_0_1500 = std::make_shared<yap::Resonance>(factory.quantumNumbers("f_0"), 1.507, "f_0_1500", radialSize, std::make_unique<yap::BreitWigner>());
+    std::shared_ptr<yap::Resonance> f_0_1500 = std::make_shared<yap::Resonance>(factory.quantumNumbers("f_0"), 1.507, "f_0_1500", radialSize, std::make_unique<yap::BreitWigner>());
     static_cast<yap::BreitWigner&>(f_0_1500->massShape()).width()->setValue(0.109);
     f_0_1500->addChannels(piPlus, piMinus, max2L);
 
     // sigma a.k.a. f_0(500)
-    auto sigma = std::make_shared<yap::Resonance>(factory.quantumNumbers("f_0"), 0.800, "sigma", radialSize, std::make_unique<yap::BreitWigner>());
+    std::shared_ptr<yap::Resonance> sigma = std::make_shared<yap::Resonance>(factory.quantumNumbers("f_0"), 0.800, "sigma", radialSize, std::make_unique<yap::BreitWigner>());
     static_cast<yap::BreitWigner&>(sigma->massShape()).width()->setValue(0.800);
     sigma->addChannels(piPlus, piMinus, max2L);
 

programs/D4piTest.cxx

@@ -17,6 +17,7 @@
 
 #include <assert.h>
 #include <iostream>
+#include <memory>
 #include <string>
 
 //#include <callgrind.h>
@@ -35,23 +36,23 @@ int main( int argc, char** argv)
 
     // initial state particle
     double radialSize = 1.;
-    auto D = factory.createInitialStateParticle(421, radialSize);
+    std::shared_ptr<yap::InitialStateParticle> D = factory.createInitialStateParticle(421, radialSize);
 
     // final state particles
-    auto piPlus = factory.createFinalStateParticle(211);
-    auto piMinus = factory.createFinalStateParticle(-211);
+    std::shared_ptr<yap::FinalStateParticle> piPlus = factory.createFinalStateParticle(211);
+    std::shared_ptr<yap::FinalStateParticle> piMinus = factory.createFinalStateParticle(-211);
 
     // Set final-state particles
     D->setFinalStateParticles({piPlus, piMinus, piPlus, piMinus});
 
     // rho rho
-    auto rho = factory.createResonance(113, radialSize, std::make_unique<yap::BreitWigner>());
+    std::shared_ptr<yap::Resonance> rho = factory.createResonance(113, radialSize, std::make_unique<yap::BreitWigner>());
     rho->addChannels(piPlus, piMinus, max2L);
 
     D->addChannels(rho, rho, max2L);
 
     // omega omega
-    auto omega = factory.createResonance(223, radialSize, std::make_unique<yap::BreitWigner>());
+    std::shared_ptr<yap::Resonance> omega = factory.createResonance(223, radialSize, std::make_unique<yap::BreitWigner>());
     omega->addChannels(piPlus, piMinus, max2L);
 
     D->addChannels(omega, omega, max2L);
@@ -60,10 +61,10 @@ int main( int argc, char** argv)
     D->addChannels(rho, omega, max2L);
 
     // a_1 channels
-    auto sigma = factory.createResonance(9000221, radialSize, std::make_unique<yap::BreitWigner>());
+    std::shared_ptr<yap::Resonance> sigma = factory.createResonance(9000221, radialSize, std::make_unique<yap::BreitWigner>());
     sigma->addChannels(piPlus, piMinus, max2L);
 
-    auto a_1 = factory.createResonance(20213, radialSize, std::make_unique<yap::BreitWigner>());
+    std::shared_ptr<yap::Resonance> a_1 = factory.createResonance(20213, radialSize, std::make_unique<yap::BreitWigner>());
     a_1->addChannels(sigma, piPlus, max2L);
 
     a_1->addChannels(rho, piPlus, max2L);

programs/DKKpiTest.cxx

@@ -24,12 +24,12 @@ int main( int argc, char** argv)
     yap::ParticleFactory factory((::getenv("YAPDIR") ? (std::string)::getenv("YAPDIR") : ".") + "/evt.pdl");
 
     // create final state particles
-    auto kPlus  = factory.createFinalStateParticle(+321);
-    auto kMinus = factory.createFinalStateParticle(-321);
-    auto piPlus = factory.createFinalStateParticle(+211);
+    std::shared_ptr<yap::FinalStateParticle> kPlus  = factory.createFinalStateParticle(+321);
+    std::shared_ptr<yap::FinalStateParticle> kMinus = factory.createFinalStateParticle(-321);
+    std::shared_ptr<yap::FinalStateParticle> piPlus = factory.createFinalStateParticle(+211);
 
     // create initial state particle and set final state
-    auto D = factory.createInitialStateParticle(factory.pdgCode("D+"), radialSize);
+    std::shared_ptr<yap::InitialStateParticle> D = factory.createInitialStateParticle(factory.pdgCode("D+"), radialSize);
     D->setFinalStateParticles({kPlus, kMinus, piPlus});
 
     // create a phi

src/ParticleFactory.cxx

@@ -4,6 +4,7 @@
 #include "FinalStateParticle.h"
 #include "InitialStateParticle.h"
 #include "logging.h"
+#include "make_unique.h"
 #include "Resonance.h"
 #include "SpinAmplitude.h"
 
@@ -42,32 +43,32 @@ ParticleFactory::ParticleFactory(const std::string pdlFile)
 }
 
 //-------------------------
-std::shared_ptr<FinalStateParticle> ParticleFactory::createFinalStateParticle(int PDG)
+std::unique_ptr<FinalStateParticle> ParticleFactory::createFinalStateParticle(int PDG)
 {
-    const ParticleTableEntry& p = particleTableEntry(PDG);
+    const auto& p = particleTableEntry(PDG);
     DEBUG("make FinalStateParticle " << p.Name_ << " with quantum numbers " << p);
-    return std::make_shared<FinalStateParticle>(p, p.Mass_, p.Name_);
+    return std::make_unique<FinalStateParticle>(p, p.Mass_, p.Name_);
 }
 
 //-------------------------
-std::shared_ptr<InitialStateParticle> ParticleFactory::createInitialStateParticle(int PDG, double radialSize)
+std::unique_ptr<InitialStateParticle> ParticleFactory::createInitialStateParticle(int PDG, double radialSize)
 {
-    const ParticleTableEntry& p = particleTableEntry(PDG);
+    const auto& p = particleTableEntry(PDG);
 
     if (p.twoJ() != 0)
         LOG(ERROR) << "InitialStateParticle has spin != 0. ";
 
     DEBUG("make InitialStateParticle " << p.Name_ << " with quantum numbers " << p);
-    return std::make_shared<InitialStateParticle>(p, p.Mass_, p.Name_, radialSize);
+    return std::make_unique<InitialStateParticle>(p, p.Mass_, p.Name_, radialSize);
 }
 
 //-------------------------
-std::shared_ptr<Resonance> ParticleFactory::createResonance(int PDG, double radialSize, std::unique_ptr<MassShape>&& massShape)
+std::unique_ptr<Resonance> ParticleFactory::createResonance(int PDG, double radialSize, std::unique_ptr<MassShape>&& massShape)
 {
-    const ParticleTableEntry& p = particleTableEntry(PDG);
+    const auto& p = particleTableEntry(PDG);
     DEBUG("make Resonance " << p.Name_ << " with quantum numbers " << p);
     massShape->setParameters(p);
-    return std::make_shared<Resonance>(p, p.Mass_, p.Name_, radialSize, std::move(massShape));
+    return std::make_unique<Resonance>(p, p.Mass_, p.Name_, radialSize, std::move(massShape));
 }
 
 //-------------------------

test/test_Vector.cxx

@@ -2,6 +2,7 @@
 
 #include <Constants.h>
 #include <FourVector.h>
+#include <ThreeVector.h>
 #include <logging.h>
 
 #include <cmath>
@@ -12,6 +13,62 @@ TEST_CASE( "Vector" )
     // disable logs in text
     yap::disableLogs(el::Level::Global);
 
+    SECTION( "ThreeVector" ) {
+
+        auto v1 = yap::ThreeVector<double>({1, 2, 3});
+        auto v2 = yap::ThreeVector<double>({4, 5, 6});
+
+        SECTION( "addition-assignment" ) {
+            v1 += v2;
+            REQUIRE( v1 == yap::ThreeVector<double>({5, 7, 9}) );
+        }
+
+        SECTION( "subtraction-assignment" ) {
+            v1 -= v2;
+            REQUIRE( v1 == yap::ThreeVector<double>({ -3, -3, -3}) );
+        }
+
+        SECTION( "multiplication-assignment" ) {
+            v1 *= 3.;
+            REQUIRE( v1 == yap::ThreeVector<double>({3, 6, 9}) );
+        }
+
+        SECTION( "arithmetic operations" ) {
+
+            // +
+            REQUIRE( v1 + v2 == yap::ThreeVector<double>({5, 7, 9}) );
+
+            // -
+            REQUIRE( v1 - v2 == yap::ThreeVector<double>({ -3, -3, -3}) );
+
+            // inner product
+            REQUIRE( v1 * v2 == 32 );
+
+            // norm
+            REQUIRE( norm(v1) == 14 );
+            REQUIRE( norm(v2) == 77 );
+
+            // abs
+            REQUIRE( abs(v1) == sqrt(14) );
+            REQUIRE( abs(v2) == sqrt(77) );
+
+            // cross product
+            auto x = yap::ThreeVector<double>({1, 0, 0});
+            auto y = yap::ThreeVector<double>({0, 1, 0});
+            auto z = yap::ThreeVector<double>({0, 0, 1});
+            REQUIRE( cross(x, y) == z );
+
+        }
+
+        SECTION( "constants" ) {
+
+            // axes:
+            REQUIRE( cross(yap::ThreeAxis_X, yap::ThreeAxis_Y) == yap::ThreeAxis_Z );
+            REQUIRE( isRightHanded(yap::ThreeAxes) );
+        }
+
+    }
+
     SECTION( "FourVector" ) {
 
         auto v1 = yap::FourVector<double>({4, 3, 2, 1});