diff --git a/src/Algorithms/ParallelCyclotronTracker.cpp b/src/Algorithms/ParallelCyclotronTracker.cpp
index 908d7dabddaa0bfe230e38b893afeabd5c6a2765..cb556ae79eb28e441d6877dd7af1f768806c96eb 100644
--- a/src/Algorithms/ParallelCyclotronTracker.cpp
+++ b/src/Algorithms/ParallelCyclotronTracker.cpp
@@ -429,19 +429,15 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
 
         // If the user wants to save the restarted run in local frame,
         // make sure the previous h5 file was local too
-        if (Options::psDumpLocalFrame) {
-
-	    if (!previousH5Local) {
-
+      if (Options::psDumpLocalFrame != Options::GLOBAL) {
+        if (!previousH5Local) {
                 throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                     "You are trying a local restart from a global h5 file!");
-	    }
+        }
             // Else, if the user wants to save the restarted run in global frame,
             // make sure the previous h5 file was global too
-        } else {
-
-	    if (previousH5Local) {
-
+      } else {
+        if (previousH5Local) {
                 throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                     "You are trying a global restart from a local h5 file!");
             }
@@ -451,9 +447,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
         referencePhi *= Physics::deg2rad;
         referencePsi *= Physics::deg2rad;
         referencePtot = bega;
-
         if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
-
             throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                 "PHIINIT is out of [-180, 180)!");
         }
@@ -3515,7 +3509,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
             itsBunch->P[i](0) += referencePr;
             itsBunch->P[i](1) += referencePt;
             itsBunch->P[i](2) += referencePz;
-	}
+        }
 
         // Out of the three coordinates of meanR (R, Theta, Z) only the angle
         // changes the momentum vector...
@@ -3545,7 +3539,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
         // Do a local frame restart (we have already checked that the old h5 file was saved in local
         // frame as well).
         // Cave: Multi-bunch must not be done in the local frame! (TODO: Is this still true? -DW)
-        if((Options::psDumpLocalFrame)) {
+        if((Options::psDumpLocalFrame != Options::GLOBAL)) {
 
             *gmsg << "* Restart in the local frame" << endl;
 
@@ -3786,8 +3780,15 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
     // --------------------------------- Get some Values ---------------------------------------- //
     double const E = itsBunch->get_meanKineticEnergy();
     double const temp_t = itsBunch->getT() * 1e9; // s -> ns
-    Vector_t const meanR = calcMeanR();
-    Vector_t const meanP = calcMeanP();
+    Vector_t meanR;
+    Vector_t meanP;
+    if (Options::psDumpLocalFrame == Options::BUNCH_MEAN) {
+        meanR = calcMeanR();
+        meanP = calcMeanP();
+    } else {
+        meanR = itsBunch->R[0];
+        meanP = itsBunch->P[0];
+    }
     double phi = 0;
     double psi = 0;
     // --------------  Calculate the external fields at the center of the bunch ----------------- //
@@ -3800,9 +3801,8 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
 
     // If we are saving in local frame, bunch and fields at the bunch center have to be rotated
     // TODO: Make decision if we maybe want to always save statistics data in local frame? -DW
-    if(Options::psDumpLocalFrame) {
+    if(Options::psDumpLocalFrame != Options::GLOBAL) {
         // -------------------- ----------- Do Transformations ---------------------------------- //
-
         // Bunch (local) azimuth at meanR w.r.t. y-axis
         phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * pi;
 
@@ -3829,7 +3829,7 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
     //itsBunch->R *= Vector_t(1000.0); // m -> mm
 
     // If we are in local mode, transform back after saving
-    if(Options::psDumpLocalFrame) {
+    if(Options::psDumpLocalFrame != Options::GLOBAL) {
         localToGlobal(itsBunch->R, phi, psi);
         localToGlobal(itsBunch->P, phi, psi);
     }
@@ -3850,8 +3850,15 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
     // --------------------------------- Get some Values ---------------------------------------- //
     double const temp_t = itsBunch->getT() * 1.0e9; // s -> ns
 
-    Vector_t const meanR = calcMeanR();
-    Vector_t const meanP = calcMeanP();
+    Vector_t meanR;
+    Vector_t meanP;
+    if (Options::psDumpLocalFrame == Options::BUNCH_MEAN) {
+        meanR = calcMeanR();
+        meanP = calcMeanP();
+    } else {
+        meanR = itsBunch->R[0];
+        meanP = itsBunch->P[0];
+    }
 
     double const betagamma_temp = sqrt(dot(meanP, meanP));
     double const E = itsBunch->get_meanKineticEnergy();
@@ -3893,7 +3900,7 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
 
     // -------------- If flag DumpLocalFrame is not set, dump bunch in global frame ------------- //
     if (Options::psDumpFreq < std::numeric_limits<int>::max() ){
-        if(!(Options::psDumpLocalFrame) && (Options::psDumpFreq < std::numeric_limits<int>::max())) {
+        if (Options::psDumpLocalFrame == Options::GLOBAL) {
 
             FDext_m[0] = extB_m * 0.1; // kgauss --> T
             FDext_m[1] = extE_m;
diff --git a/src/BasicActions/Option.cpp b/src/BasicActions/Option.cpp
index 2e531a6391f5ee573eabb06d12aa0b83e4b615cb..f8b197f4e6d6642742ef48a46e4bde2ac216e0ce 100644
--- a/src/BasicActions/Option.cpp
+++ b/src/BasicActions/Option.cpp
@@ -23,6 +23,8 @@
 #include "Utilities/ClassicRandom.h"
 #include "Utility/IpplInfo.h"
 
+#include "Utilities/OpalException.h"
+
 #include <ctime>
 #include <iostream>
 #include <limits>
@@ -49,6 +51,7 @@ namespace {
         STATDUMPFREQ,
         PSDUMPEACHTURN,
         PSDUMPLOCALFRAME,
+        PSDUMPFRAME,
         SPTDUMPFREQ,
         REPARTFREQ,
         REBINFREQ,
@@ -123,7 +126,10 @@ Option::Option():
       ("REMOTEPARTDEL", "Artifically delete the remote particle if its distance to the beam mass is larger than REMOTEPARTDEL times of the beam rms size, its default values is 0 (no delete) ",0.0);
 
     itsAttr[PSDUMPLOCALFRAME] = Attributes::makeBool
-                                ("PSDUMPLOCALFRAME", "If true, in local Cartesian frame, otherwise in global Cartesian frame, only aviable for OPAL-cycl, its default value is false");
+                                ("PSDUMPLOCALFRAME", "Deprecated; please use PSDUMPFRAME.");
+
+    itsAttr[PSDUMPFRAME] = Attributes::makeString
+                                ("PSDUMPFRAME", "Controls the frame of phase space dump in stat file and h5 file. If 'GLOBAL' OPAL will dump in the lab (global) Cartesian frame; if 'BUNCH_MEAN' OPAL will dump in the local Cartesian frame of the beam mean; if 'REFERENCE'  OPAL will dump in the local Cartesian frame of the reference particle 0. Only aviable for OPAL-cycl, its default value is 'GLOBAL'");
 
     itsAttr[SPTDUMPFREQ] = Attributes::makeReal
                            ("SPTDUMPFREQ", "The frequency to dump single particle trajectory of particles with ID = 0 & 1, its default value is 1. ");
@@ -211,7 +217,8 @@ Option::Option(const std::string &name, Option *parent):
     Attributes::setReal(itsAttr[PSDUMPFREQ], psDumpFreq);
     Attributes::setReal(itsAttr[STATDUMPFREQ], statDumpFreq);
     Attributes::setBool(itsAttr[PSDUMPEACHTURN], psDumpEachTurn);
-    Attributes::setBool(itsAttr[PSDUMPLOCALFRAME], psDumpLocalFrame);
+    Attributes::setBool(itsAttr[PSDUMPLOCALFRAME], psDumpLocalFrame_m);
+    Attributes::setString(itsAttr[PSDUMPFRAME], psDumpFrame_m);
     Attributes::setReal(itsAttr[SPTDUMPFREQ], sptDumpFreq);
     Attributes::setReal(itsAttr[SCSOLVEFREQ], scSolveFreq);
     Attributes::setReal(itsAttr[MTSSUBSTEPS], mtsSubsteps);
@@ -258,7 +265,6 @@ void Option::execute() {
     verify    = Attributes::getBool(itsAttr[VERIFY]);
     warn      = Attributes::getBool(itsAttr[WARN]);
     psDumpEachTurn =   Attributes::getBool(itsAttr[PSDUMPEACHTURN]);
-    psDumpLocalFrame = Attributes::getBool(itsAttr[PSDUMPLOCALFRAME]);
     scan = Attributes::getBool(itsAttr[SCAN]);
     rhoDump = Attributes::getBool(itsAttr[RHODUMP]);
     ebDump = Attributes::getBool(itsAttr[EBDUMP]);
@@ -270,6 +276,10 @@ void Option::execute() {
     IpplInfo::Info->on(info);
     IpplInfo::Warn->on(warn);
 
+    psDumpLocalFrame_m = Attributes::getBool(itsAttr[PSDUMPLOCALFRAME]);
+    psDumpFrame_m = Attributes::getString(itsAttr[PSDUMPFRAME]);
+    handlePsDumpFrame();
+
     if(itsAttr[ASCIIDUMP]) {
         asciidump = Attributes::getBool(itsAttr[ASCIIDUMP]);
     }
@@ -372,4 +382,31 @@ void Option::execute() {
     if(Attributes::getBool(itsAttr[TELL])) {
         *gmsg << "\nCurrent settings of options:\n" << *this << endl;
     }
-}
\ No newline at end of file
+}
+
+void Option::handlePsDumpFrame() {
+    if (psDumpLocalFrame_m) {
+        if (psDumpFrame_m == "GLOBAL") {
+            psDumpFrame_m == "BUNCH_MEAN";
+        } else {
+            std::string msg = std::string("Either set 'PSDUMPLOCALFRAME' Option")+\
+                              std::string(" or 'PSDUMPFRAME' Option but not both."); 
+            throw OpalException("Option::handlePsDumpFrame", msg);
+        }
+    }
+    if (psDumpFrame_m == "GLOBAL") {
+        // do nothing; leave as defaults (this gets called for every option,
+        // so we have to implicitly take default otherwise we will overwrite
+        // other settings)
+    } else if (psDumpFrame_m == "BUNCH_MEAN") {
+        psDumpLocalFrame = BUNCH_MEAN;
+    } else if (psDumpFrame_m == "REFERENCE") {
+        psDumpLocalFrame = REFERENCE;
+    } else {
+        std::string msg = std::string("Did not recognise PSDUMPFRAME '")+\
+                    psDumpFrame_m+std::string("'. It should be one of 'GLOBAL',")+\
+                    std::string(" 'BUNCH_MEAN' or 'REFERENCE'");
+        throw OpalException("Option::handlePsDumpFrame", msg);
+    }
+}
+
diff --git a/src/BasicActions/Option.h b/src/BasicActions/Option.h
index 637c784702d705c8520aeda4b97391e50f1eba2f..fdd39ae453b413bc29418acac4296be85f1c710e 100644
--- a/src/BasicActions/Option.h
+++ b/src/BasicActions/Option.h
@@ -43,6 +43,7 @@ public:
     virtual void execute();
 
 private:
+    void handlePsDumpFrame();
 
     // Not implemented.
     Option(const Option &);
@@ -50,6 +51,9 @@ private:
 
     // Clone constructor.
     Option(const std::string &name, Option *parent);
+
+    bool psDumpLocalFrame_m = false;
+    std::string psDumpFrame_m = "GLOBAL";
 };
 
 #endif // OPAL_Option_HH
diff --git a/src/Classic/Utilities/OptionTypes.h b/src/Classic/Utilities/OptionTypes.h
index cfde4eb555ed16291f9ab313474d2ed6099004de..c68727eb614087c454d0cdb3b20715d0622c5896 100644
--- a/src/Classic/Utilities/OptionTypes.h
+++ b/src/Classic/Utilities/OptionTypes.h
@@ -1,4 +1,9 @@
 namespace Options {
+    enum DumpFrame {
+        GLOBAL=0,
+        BUNCH_MEAN=1,
+        REFERENCE=2
+    };
 
     enum OPENMODE {
         WRITE,
diff --git a/src/Classic/Utilities/Options.cpp b/src/Classic/Utilities/Options.cpp
index 95b4ff9d8d869a51e5028429b56b68f1d9ac6fdb..639f20e0c5598a68ba97bf215340426bc3b1fe3c 100644
--- a/src/Classic/Utilities/Options.cpp
+++ b/src/Classic/Utilities/Options.cpp
@@ -30,7 +30,7 @@ namespace Options {
     bool verify = false;
     bool warn = true;
     bool psDumpEachTurn = false;
-    bool psDumpLocalFrame = false;
+    DumpFrame psDumpLocalFrame = GLOBAL;
     bool scan = false;
     bool rhoDump = false;
     bool ebDump = false;
@@ -94,4 +94,4 @@ namespace Options {
 
     // opal version of input file
     int version = 10000;
-}
\ No newline at end of file
+}
diff --git a/src/Classic/Utilities/Options.h b/src/Classic/Utilities/Options.h
index 4a97f1dec66b60a0cfa7fa65e4e1ce6558618396..b6cf3ee21145e8e911dbbd89199bfeffaa11f11a 100644
--- a/src/Classic/Utilities/Options.h
+++ b/src/Classic/Utilities/Options.h
@@ -28,7 +28,6 @@
 #include "Utilities/ClassicRandom.h"
 
 namespace Options {
-
     /// Echo flag.
     //  If true, print an input echo.
     extern bool echo;
@@ -114,8 +113,10 @@ namespace Options {
     extern int boundpDestroyFreq;
 
     /// flag to decide in which coordinate frame the phase space will be dumped for OPAL-cycl
-    // if true, in local Cartesian frame, otherwise in global Cartesian frame
-    extern bool psDumpLocalFrame;
+    //  - GLOBAL, in Cartesian frame of the global particle
+    //  - BUNCH_MEAN, in Cartesian frame of the bunch mean
+    //  - REFERENCE, in Cartesian frame of the reference (0) particle
+    extern DumpFrame psDumpLocalFrame;
 
     /// The frequency to solve space charge fields.
     extern int scSolveFreq;
@@ -158,4 +159,4 @@ namespace Options {
     extern int version;
 }
 
-#endif // OPAL_Options_HH
\ No newline at end of file
+#endif // OPAL_Options_HH
diff --git a/src/Structure/H5PartWrapperForPC.cpp b/src/Structure/H5PartWrapperForPC.cpp
index 65f1f05705d339bcf0ac0a69b79932205ba7d2e0..cd06e1c259027efb7a2490faae2b02f867901ef6 100644
--- a/src/Structure/H5PartWrapperForPC.cpp
+++ b/src/Structure/H5PartWrapperForPC.cpp
@@ -310,7 +310,7 @@ void H5PartWrapperForPC::writeStepHeader(PartBunch& bunch, const std::map<std::s
     double mass = 1.0e-9 * bunch.getM();
     double charge = bunch.getCharge();
 
-    h5_int64_t localFrame = Options::psDumpLocalFrame? 1: 0;
+    h5_int64_t localFrame = Options::psDumpLocalFrame;
 
     double sposHead = 0.0;
     double sposRef = 0.0;
@@ -585,4 +585,4 @@ void H5PartWrapperForPC::writeStepData(PartBunch& bunch) {
                                                (h5_float64_t)bunch.get_hr()(1),
                                                (h5_float64_t)bunch.get_hr()(2)));
     }
-}
\ No newline at end of file
+}
diff --git a/tests/classic_src/AbsBeamline/RingTest.cpp b/tests/classic_src/AbsBeamline/RingTest.cpp
index 7543e89ecc7c6a984332ae3502de10f23c34b7ab..8fd3c89ad7b96d015444ebbe207b40cf426ca603 100644
--- a/tests/classic_src/AbsBeamline/RingTest.cpp
+++ b/tests/classic_src/AbsBeamline/RingTest.cpp
@@ -37,15 +37,6 @@
 #include <iostream>
 #include <sstream>
 
-namespace {
-    std::string burnAfterReading(std::ostringstream &ostr) {
-        std::string returnValue = ostr.str();
-        ostr.str("");
-
-        return returnValue;
-    }
-}
-
 // generate a set of weird, but closed, elements
 // reaches theta sum after 16 elements
 class OffsetFactory {
diff --git a/tests/opal_src/Distribution/GaussTest.cpp b/tests/opal_src/Distribution/GaussTest.cpp
index ae89e60b7c536bb3c64bf2d0f187774a53db8cec..136f79e9c5562822d1587d81b2963aeb36136f8f 100644
--- a/tests/opal_src/Distribution/GaussTest.cpp
+++ b/tests/opal_src/Distribution/GaussTest.cpp
@@ -25,15 +25,6 @@ Inform *gmsg;
 #include <cstring>
 #include <sstream>
 
-namespace {
-    std::string burnAfterReading(std::ostringstream &ostr) {
-        std::string returnValue = ostr.str();
-        ostr.str("");
-
-        return returnValue;
-    }
-}
-
 TEST(GaussTest, FullSigmaTest1) {
     OpalTestUtilities::SilenceTest silencer(true);
     char inputFileName[] = "GaussDistributionTest.in";