diff --git a/src/Algorithms/BoostMatrix.h b/src/Algorithms/BoostMatrix.h index 64d2ba178142402eb224305606492aaac8b8c2b5..e6c02859a633aa4792d8d675c5427e277fcad457 100644 --- a/src/Algorithms/BoostMatrix.h +++ b/src/Algorithms/BoostMatrix.h @@ -23,20 +23,22 @@ typedef boost::numeric::ublas::matrix<double> matrix_t; template<class T> -T prod_boost_vector(boost::numeric::ublas::matrix<double> rotation, const T& vector) { - boost::numeric::ublas::vector<double> boostVector(3); - boostVector(0) = vector(0); - boostVector(1) = vector(1); - boostVector(2) = vector(2); - - boostVector = boost::numeric::ublas::prod(rotation, boostVector); - - T prodVector; - prodVector(0) = boostVector(0); - prodVector(1) = boostVector(1); - prodVector(2) = boostVector(2); - - return prodVector; +T prod_boost_vector(const boost::numeric::ublas::matrix<double>& rotation, const T& vector) { + // Ensure that 'vector' has the correct size (3x1 or 1x3) + assert(vector.size() == 3); + + // Create a temporary result vector + T result; + + // Perform matrix-vector multiplication directly without copying + for (std::size_t i = 0; i < 3; ++i) { + result(i) = 0; + for (std::size_t j = 0; j < 3; ++j) { + result(i) += rotation(i, j) * vector(j); + } + } + + return result; } #endif diff --git a/src/Algorithms/ParallelCyclotronTracker.cpp b/src/Algorithms/ParallelCyclotronTracker.cpp index 06496a758e6a7b634c9715a15f19cdd5150f045c..19fb8f4e08c520a91ed2d4ab060436abb87dc13e 100644 --- a/src/Algorithms/ParallelCyclotronTracker.cpp +++ b/src/Algorithms/ParallelCyclotronTracker.cpp @@ -60,7 +60,6 @@ #include "AbstractObjects/Element.h" #include "AbstractObjects/OpalData.h" -#include "Algorithms/BoostMatrix.h" #include "Algorithms/Ctunes.h" #include "Algorithms/PolynomialTimeDependence.h" @@ -135,6 +134,7 @@ ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline& beamline, , itsStepper_mp(nullptr) , mode_m(TrackingMode::UNDEFINED) , stepper_m(timeintegrator) + , rotation_m(3,3) { itsBeamline = dynamic_cast<Beamline*>(beamline.clone()); itsDataSink = &ds; @@ -1664,19 +1664,18 @@ void ParallelCyclotronTracker::globalToLocal(ParticleAttrib<Vector_t>& particleV IpplTimings::startTimer(TransformTimer_m); particleVectors -= translationToGlobal; - matrix_t rotation(3,3); - rotation(0,0) = std::cos(phi); - rotation(0,1) = std::sin(phi); - rotation(0,2) = 0; - rotation(1,0) = -std::sin(phi); - rotation(1,1) = std::cos(phi); - rotation(1,2) = 0; - rotation(2,0) = 0; - rotation(2,1) = 0; - rotation(2,2) = 1; + rotation_m(0,0) = std::cos(phi); + rotation_m(0,1) = std::sin(phi); + rotation_m(0,2) = 0; + rotation_m(1,0) = -std::sin(phi); + rotation_m(1,1) = std::cos(phi); + rotation_m(1,2) = 0; + rotation_m(2,0) = 0; + rotation_m(2,1) = 0; + rotation_m(2,2) = 1; for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { - particleVectors[i] = prod_boost_vector(rotation, particleVectors[i]); + particleVectors[i] = prod_boost_vector(rotation_m, particleVectors[i]); } IpplTimings::stopTimer(TransformTimer_m); } @@ -1685,19 +1684,18 @@ void ParallelCyclotronTracker::localToGlobal(ParticleAttrib<Vector_t>& particleV double phi, Vector_t const translationToGlobal) { IpplTimings::startTimer(TransformTimer_m); - matrix_t rotation(3,3); - rotation(0,0) = std::cos(phi); - rotation(0,1) = -std::sin(phi); - rotation(0,2) = 0; - rotation(1,0) = std::sin(phi); - rotation(1,1) = std::cos(phi); - rotation(1,2) = 0; - rotation(2,0) = 0; - rotation(2,1) = 0; - rotation(2,2) = 1; + rotation_m(0,0) = std::cos(phi); + rotation_m(0,1) = -std::sin(phi); + rotation_m(0,2) = 0; + rotation_m(1,0) = std::sin(phi); + rotation_m(1,1) = std::cos(phi); + rotation_m(1,2) = 0; + rotation_m(2,0) = 0; + rotation_m(2,1) = 0; + rotation_m(2,2) = 1; for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { - particleVectors[i] = prod_boost_vector(rotation, particleVectors[i]); + particleVectors[i] = prod_boost_vector(rotation_m, particleVectors[i]); } particleVectors += translationToGlobal; @@ -1859,73 +1857,69 @@ inline void ParallelCyclotronTracker::normalizeVector(Vector_t& vector) { inline void ParallelCyclotronTracker::rotateAroundZ(ParticleAttrib<Vector_t>& particleVectors, double const phi) { // Clockwise rotation of particles 'particleVectors' by 'phi' around Z axis - matrix_t rotation(3,3); - rotation(0,0) = std::cos(phi); - rotation(0,1) = std::sin(phi); - rotation(0,2) = 0; - rotation(1,0) = -std::sin(phi); - rotation(1,1) = std::cos(phi); - rotation(1,2) = 0; - rotation(2,0) = 0; - rotation(2,1) = 0; - rotation(2,2) = 1; + rotation_m(0,0) = std::cos(phi); + rotation_m(0,1) = std::sin(phi); + rotation_m(0,2) = 0; + rotation_m(1,0) = -std::sin(phi); + rotation_m(1,1) = std::cos(phi); + rotation_m(1,2) = 0; + rotation_m(2,0) = 0; + rotation_m(2,1) = 0; + rotation_m(2,2) = 1; for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { - particleVectors[i] = prod_boost_vector(rotation, particleVectors[i]); + particleVectors[i] = prod_boost_vector(rotation_m, particleVectors[i]); } } inline void ParallelCyclotronTracker::rotateAroundZ(Vector_t& myVector, double const phi) { // Clockwise rotation of single vector 'myVector' by 'phi' around Z axis - matrix_t rotation(3,3); - rotation(0,0) = std::cos(phi); - rotation(0,1) = std::sin(phi); - rotation(0,2) = 0; - rotation(1,0) = -std::sin(phi); - rotation(1,1) = std::cos(phi); - rotation(1,2) = 0; - rotation(2,0) = 0; - rotation(2,1) = 0; - rotation(2,2) = 1; + rotation_m(0,0) = std::cos(phi); + rotation_m(0,1) = std::sin(phi); + rotation_m(0,2) = 0; + rotation_m(1,0) = -std::sin(phi); + rotation_m(1,1) = std::cos(phi); + rotation_m(1,2) = 0; + rotation_m(2,0) = 0; + rotation_m(2,1) = 0; + rotation_m(2,2) = 1; - myVector = prod_boost_vector(rotation, myVector); + myVector = prod_boost_vector(rotation_m, myVector); } inline void ParallelCyclotronTracker::rotateAroundX(ParticleAttrib<Vector_t>& particleVectors, double const psi) { // Clockwise rotation of particles 'particleVectors' by 'psi' around X axis - matrix_t rotation(3,3); - rotation(0,0) = 1; - rotation(0,1) = 0; - rotation(0,2) = 0; - rotation(1,0) = 0; - rotation(1,1) = std::cos(psi); - rotation(1,2) = std::sin(psi); - rotation(2,0) = 0; - rotation(2,1) = -std::sin(psi); - rotation(2,2) = std::cos(psi); + rotation_m(0,0) = 1; + rotation_m(0,1) = 0; + rotation_m(0,2) = 0; + rotation_m(1,0) = 0; + rotation_m(1,1) = std::cos(psi); + rotation_m(1,2) = std::sin(psi); + rotation_m(2,0) = 0; + rotation_m(2,1) = -std::sin(psi); + rotation_m(2,2) = std::cos(psi); for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { - particleVectors[i] = prod_boost_vector(rotation, particleVectors[i]); + particleVectors[i] = prod_boost_vector(rotation_m, particleVectors[i]); } } inline void ParallelCyclotronTracker::rotateAroundX(Vector_t& myVector, double const psi) { // Clockwise rotation of single vector 'myVector' by 'psi' around X axis - matrix_t rotation(3,3); - rotation(0,0) = 1; - rotation(0,1) = 0; - rotation(0,2) = 0; - rotation(1,0) = 0; - rotation(1,1) = std::cos(psi); - rotation(1,2) = std::sin(psi); - rotation(2,0) = 0; - rotation(2,1) = -std::sin(psi); - rotation(2,2) = std::cos(psi); - - myVector = prod_boost_vector(rotation, myVector); + rotation_m(0,0) = 1; + rotation_m(0,1) = 0; + rotation_m(0,2) = 0; + rotation_m(1,0) = 0; + rotation_m(1,1) = std::cos(psi); + rotation_m(1,2) = std::sin(psi); + rotation_m(2,0) = 0; + rotation_m(2,1) = -std::sin(psi); + rotation_m(2,2) = std::cos(psi); + + myVector = prod_boost_vector(rotation_m, myVector); } inline void ParallelCyclotronTracker::getQuaternionTwoVectors(Vector_t u, Vector_t v, Quaternion_t& quaternion) { diff --git a/src/Algorithms/ParallelCyclotronTracker.h b/src/Algorithms/ParallelCyclotronTracker.h index a612ec2d433456ea18e5822db3c8d1c3bc478017..a338d90310ca7589b776227565e8b07cd1a68d2b 100644 --- a/src/Algorithms/ParallelCyclotronTracker.h +++ b/src/Algorithms/ParallelCyclotronTracker.h @@ -29,6 +29,7 @@ #define OPAL_ParallelCyclotronTracker_HH #include "AbsBeamline/ElementBase.h" +#include "Algorithms/BoostMatrix.h" #include "Algorithms/MultiBunchHandler.h" #include "Algorithms/Tracker.h" #include "Steppers/Steppers.h" @@ -472,6 +473,9 @@ private: Steppers::TimeIntegrator stepper_m; + // for change of coordinates + matrix_t rotation_m; + /// Check if turn done bool isTurnDone();