Resolve "Formula momentum conversion"

src/Algorithms/CavityAutophaser.cpp

@@ -72,7 +72,7 @@ double CavityAutophaser::getPhaseAtMaxEnergy(const Vector_t &R,
                    << std::left << std::setw(68) << std::setfill('*') << ss.str()
                    << std::setfill(' ') << endl);
     if (!apVeto) {
-        double initialEnergy = Util::getEnergy(P, itsReference_m.getM()) * 1e-6;
+        double initialEnergy = Util::getKineticEnergy(P, itsReference_m.getM()) * 1e-6;
         double AstraPhase    = 0.0;
         double designEnergy  = element->getDesignEnergy();
 
@@ -191,7 +191,7 @@ double CavityAutophaser::guessCavityPhase(double t) {
         return orig_phi;
     }
 
-    Phimax = element->getAutoPhaseEstimate(getEnergyMeV(refP),
+    Phimax = element->getAutoPhaseEstimate(Util::getKineticEnergy(refP, itsReference_m.getM()) * 1e-6,
                                            t,
                                            itsReference_m.getQ(),
                                            itsReference_m.getM() * 1e-6);
@@ -288,7 +288,7 @@ double CavityAutophaser::track(Vector_t /*R*/,
                                                             out);
     rfc->setPhasem(initialPhase);
 
-    double finalKineticEnergy = Util::getEnergy(Vector_t(0.0, 0.0, pe.first), itsReference_m.getM() * 1e-6);
+    double finalKineticEnergy = Util::getKineticEnergy(Vector_t(0.0, 0.0, pe.first), itsReference_m.getM() * 1e-6);
 
     return finalKineticEnergy;
 }
\ No newline at end of file

src/Algorithms/CavityAutophaser.h

@@ -28,8 +28,6 @@ private:
                  const double dt,
                  const double phase,
                  std::ofstream *out = NULL) const;
-    double getEnergyMeV(const Vector_t &P);
-    double getMomentum(double kineticEnergyMeV);
 
     const PartData &itsReference_m;
     std::shared_ptr<Component> itsCavity_m;
@@ -39,14 +37,4 @@ private:
 
 };
 
-inline
-double CavityAutophaser::getEnergyMeV(const Vector_t &P) {
-    return itsReference_m.getM() * 1e-6 * (sqrt(dot(P,P) + 1) - 1);
-}
-
-inline
-double CavityAutophaser::getMomentum(double kineticEnergyMeV) {
-    return sqrt(std::pow(kineticEnergyMeV / (itsReference_m.getM() * 1e-6) + 1, 2) - 1);
-}
-
 #endif
\ No newline at end of file

src/Algorithms/ParallelTTracker.cpp

@@ -1162,7 +1162,7 @@ void ParallelTTracker::findStartPosition(const BorisPusher &pusher) {
     selectDT();
 
     if ((back_track && itsOpalBeamline_m.containsSource()) ||
-        Util::getEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
+        Util::getKineticEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
         double gamma = 0.1 / itsBunch_m->getM() + 1.0;
         double beta = sqrt(1.0 - 1.0 / std::pow(gamma, 2));
         itsBunch_m->RefPartP_m = itsBunch_m->toLabTrafo_m.rotateTo(beta * gamma * Vector_t(0, 0, 1));

src/Classic/AbsBeamline/RFCavity.cpp

@@ -494,22 +494,22 @@ ElementBase::ElementType RFCavity::getType() const {
 double RFCavity::getAutoPhaseEstimateFallback(double E0, double t0, double q, double mass) {
 
     const double dt = 1e-13;
-    const double p0 = Util::getP(E0, mass);
+    const double p0 = Util::getBetaGamma(E0, mass);
     const double origPhase =getPhasem();
     double dphi = Physics::pi / 18;
 
     double phi = 0.0;
     setPhasem(phi);
-    std::pair<double, double> ret = trackOnAxisParticle(E0 / mass, t0, dt, q, mass);
+    std::pair<double, double> ret = trackOnAxisParticle(p0, t0, dt, q, mass);
     double phimax = 0.0;
-    double Emax = Util::getEnergy(Vector_t(0.0, 0.0, ret.first), mass);
+    double Emax = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
     phi += dphi;
 
     for (unsigned int j = 0; j < 2; ++ j) {
         for (unsigned int i = 0; i < 36; ++ i, phi += dphi) {
             setPhasem(phi);
             ret = trackOnAxisParticle(p0, t0, dt, q, mass);
-            double Ekin = Util::getEnergy(Vector_t(0.0, 0.0, ret.first), mass);
+            double Ekin = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
             if (Ekin > Emax) {
                 Emax = Ekin;
                 phimax = phi;
@@ -635,7 +635,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
         }
 
         double cosine_part = 0.0, sine_part = 0.0;
-        double p0 = std::sqrt((E0 / mass + 1) * (E0 / mass + 1) - 1);
+        double p0 = Util::getBetaGamma(E0, mass);
         cosine_part += scale_m * std::cos(frequency_m * t0) * F[0];
         sine_part += scale_m * std::sin(frequency_m * t0) * F[0];
 
@@ -675,11 +675,11 @@ std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
     const double zend = length_m + startField_m;
 
     Vector_t z(0.0, 0.0, zbegin);
-    double dz = 0.5 * p(2) / std::sqrt(1.0 + dot(p, p)) * cdt;
+    double dz = 0.5 * p(2) / Util::getGamma(p) * cdt;
     Vector_t Ef(0.0), Bf(0.0);
 
     if (out) *out << std::setw(18) << z[2]
-                  << std::setw(18) << Util::getEnergy(p, mass)
+                  << std::setw(18) << Util::getKineticEnergy(p, mass)
                   << std::endl;
     while (z(2) + dz < zend && z(2) + dz > zbegin) {
         z /= cdt;
@@ -700,7 +700,7 @@ std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
         t += dt;
 
         if (out) *out << std::setw(18) << z[2]
-                      << std::setw(18) << Util::getEnergy(p, mass)
+                      << std::setw(18) << Util::getKineticEnergy(p, mass)
                       << std::endl;
     }
 

src/Classic/Utilities/Util.h

@@ -6,6 +6,7 @@
 
 #include <string>
 #include <cstring>
+#include <limits>
 #include <sstream>
 #include <type_traits>
 #include <functional>
@@ -28,16 +29,23 @@ namespace Util {
     }
 
     inline
-    double getEnergy(Vector_t p, double mass) {
+    double getKineticEnergy(Vector_t p, double mass) {
         return (getGamma(p) - 1.0) * mass;
     }
 
     inline
-    double getP(double E, double mass) {
-        double gamma = E / mass + 1;
-        return std::sqrt(std::pow(gamma, 2.0) - 1.0);
+    double getBetaGamma(double Ekin, double mass) {
+        double value = std::sqrt(std::pow(Ekin / mass + 1.0, 2) - 1.0);
+        if (value < std::numeric_limits<double>::epsilon())
+            value = std::sqrt(2 * Ekin / mass);
+        return value;
     }
 
+    inline
+    double convertMomentumeVToBetaGamma(double p, double mass) {
+        return p / mass;
+    }
+  
     inline
     std::string getTimeString(double time, unsigned int precision = 3) {
         std::string timeUnit(" [ps]");
@@ -162,7 +170,7 @@ namespace Util {
     }
 
     Vector_t getTaitBryantAngles(Quaternion rotation, const std::string &elementName = "");
-
+  
     std::string toUpper(const std::string &str);
 
     std::string combineFilePath(std::initializer_list<std::string>);

src/Distribution/Distribution.cpp

@@ -17,7 +17,6 @@
 #include <string>
 #include <vector>
 #include <numeric>
-#include <limits>
 
 // IPPL
 #include "DataSource/DataConnect.h"
@@ -885,14 +884,6 @@ void Distribution::chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUn
 
 }
 
-double Distribution::converteVToBetaGamma(double valueIneV, double massIneV) {
-    double value = std::copysign(std::sqrt(std::pow(std::abs(valueIneV) / massIneV + 1.0, 2) - 1.0), valueIneV);
-    if (std::abs(value) < std::numeric_limits<double>::epsilon())
-        value = std::copysign(std::sqrt(2 * std::abs(valueIneV) / massIneV), valueIneV);
-
-    return value;
-}
-
 void Distribution::createDistributionBinomial(size_t numberOfParticles, double massIneV) {
 
     setDistParametersBinomial(massIneV);
@@ -1079,9 +1070,9 @@ void Distribution::createDistributionFromFile(size_t /*numberOfParticles*/, doub
                 saveProcessor = 0;
 
             if (inputMoUnits_m == InputMomentumUnitsT::EV) {
-                P(0) = converteVToBetaGamma(P(0), massIneV);
-                P(1) = converteVToBetaGamma(P(1), massIneV);
-                P(2) = converteVToBetaGamma(P(2), massIneV);
+                P(0) = Util::convertMomentumeVToBetaGamma(P(0), massIneV);
+                P(1) = Util::convertMomentumeVToBetaGamma(P(1), massIneV);
+                P(2) = Util::convertMomentumeVToBetaGamma(P(2), massIneV);
             }
 
             pmean_m += P;
@@ -1480,7 +1471,7 @@ void Distribution::createOpalT(PartBunchBase<double, 3> *beam,
     // This is PC from BEAM
     double deltaP = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETP]);
     if (inputMoUnits_m == InputMomentumUnitsT::EV) {
-        deltaP = converteVToBetaGamma(deltaP, beam->getM());
+        deltaP = Util::convertMomentumeVToBetaGamma(deltaP, beam->getM());
     }
 
     avrgpz_m = beam->getP()/beam->getM() + deltaP;
@@ -3577,9 +3568,9 @@ void Distribution::setSigmaP_m(double massIneV) {
 
     // Check what input units we are using for momentum.
     if (inputMoUnits_m == InputMomentumUnitsT::EV) {
-        sigmaP_m[0] = converteVToBetaGamma(sigmaP_m[0], massIneV);
-        sigmaP_m[1] = converteVToBetaGamma(sigmaP_m[1], massIneV);
-        sigmaP_m[2] = converteVToBetaGamma(sigmaP_m[2], massIneV);
+        sigmaP_m[0] = Util::convertMomentumeVToBetaGamma(sigmaP_m[0], massIneV);
+        sigmaP_m[1] = Util::convertMomentumeVToBetaGamma(sigmaP_m[1], massIneV);
+        sigmaP_m[2] = Util::convertMomentumeVToBetaGamma(sigmaP_m[2], massIneV);
     }
 }
 
@@ -3912,14 +3903,14 @@ void Distribution::setupEmissionModel(PartBunchBase<double, 3> *beam) {
 void Distribution::setupEmissionModelAstra(PartBunchBase<double, 3> *beam) {
 
     double wThermal = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]));
-    pTotThermal_m = converteVToBetaGamma(wThermal, beam->getM());
+    pTotThermal_m = Util::getBetaGamma(wThermal, beam->getM());
     pmean_m = Vector_t(0.0, 0.0, 0.5 * pTotThermal_m);
 }
 
 void Distribution::setupEmissionModelNone(PartBunchBase<double, 3> *beam) {
 
     double wThermal = std::abs(Attributes::getReal(itsAttr[Attrib::Distribution::EKIN]));
-    pTotThermal_m = converteVToBetaGamma(wThermal, beam->getM());
+    pTotThermal_m = Util::getBetaGamma(wThermal, beam->getM());
     double avgPz = std::accumulate(pzDist_m.begin(), pzDist_m.end(), 0.0);
     size_t numParticles = pzDist_m.size();
     reduce(avgPz, avgPz, OpAddAssign());
@@ -4073,9 +4064,9 @@ void Distribution::shiftDistCoordinates(double massIneV) {
 
         // Check input momentum units.
         if (inputMoUnits_m == InputMomentumUnitsT::EV) {
-            deltaPx = converteVToBetaGamma(deltaPx, massIneV);
-            deltaPy = converteVToBetaGamma(deltaPy, massIneV);
-            deltaPz = converteVToBetaGamma(deltaPz, massIneV);
+            deltaPx = Util::convertMomentumeVToBetaGamma(deltaPx, massIneV);
+            deltaPy = Util::convertMomentumeVToBetaGamma(deltaPy, massIneV);
+            deltaPz = Util::convertMomentumeVToBetaGamma(deltaPz, massIneV);
         }
 
         size_t endIdx = startIdx + particlesPerDist_m[i];
@@ -4356,8 +4347,8 @@ void Distribution::adjustPhaseSpace(double massIneV) {
     double deltaPy = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETPY]);
     // Check input momentum units.
     if (inputMoUnits_m == InputMomentumUnitsT::EV) {
-        deltaPx = converteVToBetaGamma(deltaPx, massIneV);
-        deltaPy = converteVToBetaGamma(deltaPy, massIneV);
+        deltaPx = Util::convertMomentumeVToBetaGamma(deltaPx, massIneV);
+        deltaPy = Util::convertMomentumeVToBetaGamma(deltaPy, massIneV);
     }
 
     double avrg[6];

src/Distribution/Distribution.h

@@ -291,7 +291,6 @@ private:
     void checkIfEmitted();
     void checkParticleNumber(size_t &numberOfParticles);
     void chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits);
-    double converteVToBetaGamma(double valueIneV, double massIneV);
     size_t getNumberOfParticlesInFile(std::ifstream &inputFile);
 
     class BinomialBehaviorSplitter {