WIP: Resolve "Add phase-space halo parameter to stat and smb file"

Closes #379 (closed)

src/Algorithms/MultiBunchHandler.cpp

@@ -401,7 +401,7 @@ bool MultiBunchHandler::calcBunchBeamParameters(PartBunchBase<double, 3>* beam,
      * <x^3>, <y^3>, <z^3>, <x^4>, <y^4>, <z^4>
      */
     const unsigned int dim = PartBunchBase<double, 3>::Dimension;
-    std::vector<double> local(7 * dim + 1);
+    std::vector<double> local(14 * dim + 1);
 
     beaminfo_t& binfo = getBunchInfo(bunchNr);
 
@@ -432,7 +432,8 @@ bool MultiBunchHandler::calcBunchBeamParameters(PartBunchBase<double, 3>* beam,
             local[i + 2 * dim + 1] += r2;
 
             // <p_x^2>, <p_y^2>, <p_z^2>
-            local[i + 3 * dim + 1] += p * p;
+            double p2 = p * p;
+            local[i + 3 * dim + 1] += p2;
 
             // <xp_x>, <y_py>, <zp_z>
             local[i + 4 * dim + 1] += r * p;
@@ -440,8 +441,30 @@ bool MultiBunchHandler::calcBunchBeamParameters(PartBunchBase<double, 3>* beam,
             // <x^3>, <y^3>, <z^3>
             local[i + 5 * dim + 1] += r2 * r;
 
+            // <p_x^3>, <p_y^3>, <p_z^3>
+            local[i + 6 * dim + 1] += p2 * p;
+
+            // <x^2p_x>, <y^2p_y>, <z^2p_z>
+            local[i + 7 * dim + 1] += r2 * p;
+
+            // <xp_x^2>, <yp_y^2>, <zp_z^2>
+            local[i + 8 * dim + 1] += r * p2;
+
             // <x^4>, <y^4>, <z^4>
-            local[i + 6 * dim + 1] += r2 * r2;
+            local[i + 9 * dim + 1] += r2 * r2;
+
+            // <p_x^4>, <p_y^4>, <p_z^4>
+            local[i + 10 * dim + 1] += p2 * p2;
+
+            // <x^2p_x^2>, <y^2p_y^2>, <z^2p_z^2>
+            local[i + 11 * dim + 1] += r2 * p2;
+
+            // <x^3p_x>, <y^3p_y>, <z^3p_z>
+            local[i + 12 * dim + 1] += r2 * r * p;
+
+            // <xp_x^3>, <yp_y^3>, <zp_z^3>
+            local[i + 13 * dim + 1] += r * p2 * p;
+
         }
     }
 
@@ -475,13 +498,20 @@ bool MultiBunchHandler::calcBunchBeamParameters(PartBunchBase<double, 3>* beam,
 
     for (unsigned int i = 0; i < dim; ++i) {
 
-        double w   = local[i + 1] * invN;
-        double pw  = local[i + dim + 1] * invN;
-        double w2  = local[i + 2 * dim + 1] * invN;
-        double pw2 = local[i + 3 * dim + 1] * invN;
-        double wpw = local[i + 4 * dim + 1] * invN;
-        double w3  = local[i + 5 * dim + 1] * invN;
-        double w4  = local[i + 6 * dim + 1] * invN;
+        double w     = local[i + 1] * invN;
+        double pw    = local[i +      dim + 1] * invN;
+        double w2    = local[i +  2 * dim + 1] * invN;
+        double pw2   = local[i +  3 * dim + 1] * invN;
+        double wpw   = local[i +  4 * dim + 1] * invN;
+        double w3    = local[i +  5 * dim + 1] * invN;
+        double pw3   = local[i +  6 * dim + 1] * invN;
+        double w2pw  = local[i +  7 * dim + 1] * invN;
+        double wpw2  = local[i +  8 * dim + 1] * invN;
+        double w4    = local[i +  9 * dim + 1] * invN;
+        double pw4   = local[i + 10 * dim + 1] * invN;
+        double w2pw2 = local[i + 11 * dim + 1] * invN;
+        double w3pw  = local[i + 12 * dim + 1] * invN;
+        double wpw3  = local[i + 13 * dim + 1] * invN;
 
         // <x>, <y>, <z>
         binfo.mean[i] = w;
@@ -502,12 +532,52 @@ bool MultiBunchHandler::calcBunchBeamParameters(PartBunchBase<double, 3>* beam,
         binfo.emit[i] = (binfo.rrms[i] * binfo.prms[i] - wpw * wpw);
         binfo.emit[i] =  std::sqrt(std::max(binfo.emit[i], 0.0));
 
+        /* The halo parameters are computed according to the formulas of
+         * Allen, C. K. and Wangler, T. P., Phys. Rev. ST Accel. Beams 5 (2002) 124202
+         *
+         * Note: We compute all quantities with central moments!
+         */
+
         // central: <w^4> - 4 * <w> * <w^3> + 6 * <w>^2 * <w^2> - 3 * <w>^4
-        double tmp = w4
-                   - 4.0 * w * w3
-                   + 6.0 * w * w * w2
-                   - 3.0 * w * w * w * w;
-        binfo.halo[i] = tmp / ( binfo.rrms[i] * binfo.rrms[i] );
+        w4 = w4
+           - 4.0 * w * w3
+           + 6.0 * w * w * w2
+           - 3.0 * w * w * w * w;
+
+        binfo.halo[2 * i] = w4 / ( binfo.rrms[i] * binfo.rrms[i] ) - Options::haloShift;
+
+        pw4 = pw4
+            - 4.0 * pw * pw3
+            + 6.0 * pw * pw * pw2
+            - 3.0 * pw * pw * pw * pw;
+
+        w2pw2 = w2pw2
+              - 3.0 * w * w * pw * pw
+              - 2.0 * w2pw * pw
+              - 2.0 * wpw2 * w
+              + 4.0 * wpw * w * pw
+              + w2 * pw * pw
+              + pw2 * w * w;
+
+        wpw3 = wpw3
+             - 3.0 * wpw2 * pw
+             + 3.0 * wpw * pw * pw
+             - 3.0 * pw * pw * pw * w
+             + 3.0 * pw2 * pw * w
+             - pw3 * w;
+
+        w3pw = w3pw
+             - 3.0 * w2pw * w
+             + 3.0 * wpw * w * w
+             - 3.0 * w * w * w * pw
+             + 3.0 * w2 * w * pw
+             - w3 * pw;
+
+        double I2 = binfo.rrms[i] * binfo.prms[i] - wpw * wpw;
+
+        double I4 = w4 * pw4 + 3.0 * w2pw2 * w2pw2 - 4.0 * wpw3 * w3pw;
+
+        binfo.halo[2 * i + 1] = 0.5 * std::sqrt(3.0 * I4) / I2 - Options::haloShift;
 
         // central: sqrt(<w^2> - <w>^2) (w = x, y, z)
         binfo.rrms[i] = std::sqrt(binfo.rrms[i]);