autophase algorithm doesn't like negative amplitudes

src/Algorithms/CavityAutophaser.cpp

@@ -35,21 +35,28 @@ double CavityAutophaser::getPhaseAtMaxEnergy(const Vector_t &R,
     }
 
     initialP_m = Vector_t(0, 0, euclidean_norm(P));
-    double tErr  = (initialR_m(2) - R(2)) * sqrt(dot(P,P) + 1.0) / (P(2) * Physics::c);
-    double initialEnergy = Util::getEnergy(P, itsReference_m.getM()) * 1e-6;
-    double newPhase = 0.0, AstraPhase = 0.0;
-    double initialPhase = guessCavityPhase(t + tErr);
-    double optimizedPhase = 0.0;
-    double finalEnergy = 0.0;
 
-    RFCavity *element = static_cast<RFCavity *>(itsCavity_m.get());
-    double amplitude     = element->getAmplitudem();
-    double designEnergy  = element->getDesignEnergy();
-    double originalPhase = element->getPhasem();
-    bool   apVeto        = element->getAutophaseVeto();
-    double basePhase = std::fmod(element->getFrequencym() * (t + tErr), Physics::two_pi);
+    RFCavity *element     = static_cast<RFCavity *>(itsCavity_m.get());
+    bool apVeto           = element->getAutophaseVeto();
+    double originalPhase  = element->getPhasem();
+    double tErr           = (initialR_m(2) - R(2)) * sqrt(dot(P,P) + 1.0) / (P(2) * Physics::c);
+    double optimizedPhase = 0.0;
+    double finalEnergy    = 0.0;
+    double newPhase       = 0.0;
+    double amplitude      = element->getAmplitudem();
+    double basePhase      = std::fmod(element->getFrequencym() * (t + tErr), Physics::two_pi);
 
     if (!apVeto) {
+        double initialEnergy = Util::getEnergy(P, itsReference_m.getM()) * 1e-6;
+        double AstraPhase    = 0.0;
+        double initialPhase  = guessCavityPhase(t + tErr);
+        double designEnergy  = element->getDesignEnergy();
+
+        if (amplitude < 0.0) {
+            amplitude = -amplitude;
+            element->setAmplitudem(amplitude);
+        }
+
         if (amplitude == 0.0 && designEnergy <= 0.0) {
             throw OpalException("CavityAutophaser::getPhaseAtMaxEnergy()",
                                 "neither amplitude or design energy given to cavity " + element->getName());

src/Classic/AbsBeamline/RFCavity.cpp

@@ -604,6 +604,35 @@ ElementBase::ElementType RFCavity::getType() const {
     return RFCAVITY;
 }
 
+double RFCavity::getAutoPhaseEstimateFallback(double E0, double t0, double q, double mass) {
+    const double dt = 1e-13;
+    const double dphi = pi / 18;
+    const double p0 = Util::getP(E0, mass);
+
+    double phi = 0.0;
+    setPhasem(phi);
+    std::pair<double, double> ret = trackOnAxisParticle(E0 / mass, t0, dt, q, mass);
+    double phimax = ret.first;
+    double Emax = ret.second;
+    phi += dphi;
+
+    for (unsigned int i = 1; i < 36; ++ i, phi += dphi) {
+        setPhasem(phi);
+        ret = trackOnAxisParticle(p0, t0, dt, q, mass);
+        if (ret.second > Emax) {
+            Emax = ret.second;
+            phimax = ret.second;
+        }
+    }
+
+    const int prevPrecision = Ippl::Info->precision(8);
+    INFOMSG(level2 << "estimated phase= " << phimax << " rad = "
+            << phimax * Physics::rad2deg << " deg \n"
+            << "Ekin= " << Emax << " MeV" << setprecision(prevPrecision) << endl);
+
+    return phimax;
+}
+
 double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const double &q, const double &mass) {
     vector<double> t, E, t2, E2;
     std::vector< double > F;
@@ -695,8 +724,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
             E2[i] += q * scale_m * getdE(i, t2, dz, phi + dphi, frequency_m, F);
             double a = E[i], b = E2[i];
             if (std::isnan(a) || std::isnan(b)) {
-                throw GeneralClassicException("RFCavity::getAutoPhaseEstimate",
-                                              "solution is nan");
+                return getAutoPhaseEstimateFallback(E0, t0, q, mass);
             }
             t[i] = t[i - 1] + getdT(i, E, dz, mass);
             t2[i] = t2[i - 1] + getdT(i, E2, dz, mass);

src/Classic/AbsBeamline/RFCavity.h

@@ -97,12 +97,13 @@ public:
     virtual bool getAutophaseVeto() const;
 
     virtual double getAutoPhaseEstimate(const double & E0, const double & t0, const double & q, const double & m);
+    virtual double getAutoPhaseEstimateFallback(double E0, double t0, double q, double m);
 
     virtual std::pair<double, double> trackOnAxisParticle(const double & p0,
-                                                  const double & t0,
-                                                  const double & dt,
-                                                  const double & q,
-                                                  const double & mass);
+                                                          const double & t0,
+                                                          const double & dt,
+                                                          const double & q,
+                                                          const double & mass);
 
     virtual void addKR(int i, double t, Vector_t &K);
 
@@ -503,4 +504,4 @@ std::string RFCavity::getFrequencyModelName() {
     return frequency_name_m;
 }
 
-#endif // CLASSIC_RFCavity_HH
+#endif // CLASSIC_RFCavity_HH
\ No newline at end of file

src/Classic/Utilities/Util.h

@@ -26,6 +26,12 @@ namespace Util {
         return (getGamma(p) - 1.0) * mass;
     }
 
+    inline
+    double getP(double E, double mass) {
+        double gamma = E / mass + 1;
+        return sqrt(std::pow(gamma, 2.0) - 1.0);
+    }
+
     inline
     std::string getTimeString(double time, unsigned int precision = 3) {
         std::string timeUnit(" [ps]");
@@ -127,12 +133,12 @@ namespace Util {
         std::string chargeUnit(" [fC]");
 
         charge *= 1e15;
-        
+
         if (std::abs(charge) > 1000.0) {
             charge /= 1000.0;
             chargeUnit = std::string(" [pC]");
         }
-        
+
         if (std::abs(charge) > 1000.0) {
             charge /= 1000.0;
             chargeUnit = std::string(" [nC]");