Merge branch '609-transversely-shifted-distribution-in-temporal-monitors-with-fixes' into 'master'

Resolve "Transversely shifted distribution in temporal monitors"  with additional fixes

Closes #609

See merge request !449

src/Classic/AbsBeamline/Monitor.cpp

@@ -74,15 +74,17 @@ bool Monitor::apply(const size_t &i, const double &t, Vector_t &/*E*/, Vector_t
     const Vector_t &R = RefPartBunch_m->R[i];
     const Vector_t &P = RefPartBunch_m->P[i];
     const double &dt = RefPartBunch_m->dt[i];
-    const double recpgamma = Physics::c * dt / Util::getGamma(P);
-    const double middle = 0.5 * getElementLength();
+    const Vector_t singleStep  = Physics::c * dt * Util::getBeta(P);
     if (online_m && type_m == SPATIAL) {
-        if (dt * R(2) < dt * middle && // multiply with dt to allow back tracking
-            dt * (R(2) + P(2) * recpgamma) > dt * middle) {
-            double frac = (middle - R(2)) / (P(2) * recpgamma);
-
-            lossDs_m->addParticle(R + frac * recpgamma * P,
-                                  P, RefPartBunch_m->ID[i], t + frac * dt, 0);
+        if (dt * R(2) < 0.0 &&
+            dt * (R(2) + singleStep(2)) > 0.0) {
+            double frac = R(2) / singleStep(2);
+
+            lossDs_m->addParticle(R + frac * singleStep,
+                                  P,
+                                  RefPartBunch_m->ID[i],
+                                  t + frac * dt,
+                                  0);
         }
     }
 
@@ -96,15 +98,15 @@ bool Monitor::applyToReferenceParticle(const Vector_t &R,
                                        Vector_t &) {
     if (!OpalData::getInstance()->isInPrepState()) {
         const double dt = RefPartBunch_m->getdT();
-        const double recpgamma = Physics::c * dt / Util::getGamma(P);
-        const double middle = 0.5 * getElementLength();
+        const double cdt = Physics::c * dt;
+        const Vector_t singleStep = cdt * Util::getBeta(P);
 
-        if (dt * R(2) < dt * middle && // multiply with dt to allow back tracking
-            dt * (R(2) + P(2) * recpgamma) > dt * middle) {
-            double frac = (middle - R(2)) / (P(2) * recpgamma);
+        if (dt * R(2) < 0.0 &&
+            dt * (R(2) + singleStep(2)) > 0.0) {
+            double frac = -R(2) / singleStep(2);
             double time = t + frac * dt;
-            Vector_t dR = frac * P * recpgamma;
-            double ds = euclidean_norm(dR);
+            Vector_t dR = frac * singleStep;
+            double ds = euclidean_norm(dR + 0.5 * singleStep);
             lossDs_m->addReferenceParticle(csTrafoGlobal2Local_m.transformFrom(R + dR),
                                            csTrafoGlobal2Local_m.rotateFrom(P),
                                            time,
@@ -115,8 +117,8 @@ bool Monitor::applyToReferenceParticle(const Vector_t &R,
                 const unsigned int localNum = RefPartBunch_m->getLocalNum();
 
                 for (unsigned int i = 0; i < localNum; ++ i) {
-                    const double recpgamma = Physics::c * dt / Util::getGamma(RefPartBunch_m->P[i]);
-                    Vector_t shift = frac * recpgamma * RefPartBunch_m->P[i] - Vector_t(0, 0, middle);
+                    Vector_t shift = ((frac - 0.5) * cdt * Util::getBeta(RefPartBunch_m->P[i])
+                                      - singleStep);
                     lossDs_m->addParticle(RefPartBunch_m->R[i] + shift,
                                           RefPartBunch_m->P[i],
                                           RefPartBunch_m->ID[i],

src/Classic/AbsBeamline/Monitor.h

@@ -107,6 +107,8 @@ public:
     static void writeStatistics();
 
     virtual int getRequiredNumberOfTimeSteps() const override;
+
+    virtual bool isInside(const Vector_t &r) const override;
 private:
 
     // Not implemented.
@@ -133,4 +135,12 @@ int Monitor::getRequiredNumberOfTimeSteps() const
     return 1;
 }
 
+inline
+bool Monitor::isInside(const Vector_t &r) const
+{
+    const double length = getElementLength();
+    return std::abs(r(2)) <= 0.5 * length && isInsideTransverse(r);
+}
+
+
 #endif // CLASSIC_Monitor_HH
\ No newline at end of file

src/Classic/Utilities/Util.h

@@ -28,6 +28,11 @@ namespace Util {
         return std::sqrt(dot(p, p) + 1.0);
     }
 
+    inline
+    Vector_t getBeta(Vector_t p) {
+        return p / getGamma(p);
+    }
+
     inline
     double getKineticEnergy(Vector_t p, double mass) {
         return (getGamma(p) - 1.0) * mass;
@@ -45,7 +50,7 @@ namespace Util {
     double convertMomentumeVToBetaGamma(double p, double mass) {
         return p / mass;
     }
-  
+
     inline
     std::string getTimeString(double time, unsigned int precision = 3) {
         std::string timeUnit(" [ps]");
@@ -170,7 +175,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/Elements/OpalBeamline.cpp

@@ -291,9 +291,6 @@ void OpalBeamline::compute3DLattice() {
         if (element->getType() == ElementBase::SOURCE) {
             beginThis3D(2) -= thisLength;
         }
-        if (element->getType() == ElementBase::MONITOR) {
-            beginThis3D(2) -= 0.5 * thisLength;
-        }
 
         Vector_t endThis3D;
         if (element->getType() == ElementBase::SBEND ||