refactor code: force more structure to check routine; streamline and avoid...

refactor code: force more structure to check routine; streamline and avoid reduce by adding reference particle type to PartBunchBase

src/Algorithms/ParallelCyclotronTracker.cpp

@@ -90,7 +90,7 @@ using Physics::pi;
 using Physics::q_e;
 
 
-const double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
+constexpr double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
 
 Vector_t const ParallelCyclotronTracker::xaxis = Vector_t(1.0, 0.0, 0.0);
 Vector_t const ParallelCyclotronTracker::yaxis = Vector_t(0.0, 1.0, 0.0);
@@ -952,7 +952,7 @@ void ParallelCyclotronTracker::visitSBend(const SBend &bend) {
  * @param sep
  */
 void ParallelCyclotronTracker::visitSeparator(const Separator &sep) {
-    *gmsg << "In Seapator L= " << sep.getElementLength() << " however the element is missing " << endl;
+    *gmsg << "In Separator L= " << sep.getElementLength() << " however the element is missing " << endl;
     myElements.push_back(dynamic_cast<Separator *>(sep.clone()));
 }
 
@@ -2936,8 +2936,7 @@ std::tuple<double, double, double> ParallelCyclotronTracker::initializeTracking_
 
     double harm       = getHarmonicNumber();
     double dt         = itsBunch_m->getdT() * 1.0e9 * harm; // time step size (s --> ns)
-    double t          = itsBunch_m->getT() * 1.0e9;               // current time   (s --> ns)
-
+    double t          = itsBunch_m->getT()  * 1.0e9;        // current time   (s --> ns)
 
     double oldReferenceTheta      = referenceTheta * Physics::deg2rad; // init here, reset each step
     setup_m.deltaTheta            = pi / (setup_m.stepsPerTurn);    // half of the average angle per step

src/Classic/AbsBeamline/CCollimator.cpp

@@ -36,11 +36,7 @@ void CCollimator::accept(BeamlineVisitor &visitor) const {
     visitor.visitCCollimator(*this);
 }
 
-// rectangle collimators in cyclotron cylindrical coordinates
-// without particlematterinteraction, the particle hitting collimator is deleted directly
-bool CCollimator::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {
-
-    bool flagNeedUpdate = false;
+bool CCollimator::doPreCheck(PartBunchBase<double, 3> *bunch) {
     Vector_t rmin, rmax;
     bunch->get_bounds(rmin, rmax);
 
@@ -54,26 +50,41 @@ bool CCollimator::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumbe
         double rbunch_max = std::hypot(xmax, ymax);
 
         if( rbunch_max > rmin_m && rbunch_min < rmax_m ){ // check similar to z
-            size_t tempnum = bunch->getLocalNum();
-            int pflag = 0;
-            // now check each particle in bunch
-            for (unsigned int i = 0; i < tempnum; ++i) {
-                if (bunch->PType[i] == ParticleType::REGULAR && bunch->R[i](2) < zend_m && bunch->R[i](2) > zstart_m ) {
-                    // only now careful check in r
-                    pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
-                    /// bunch->Bin[i] != -1 makes sure the particle is not stored in more than one collimator
-                    if ((pflag != 0) && (bunch->Bin[i] != -1)) {
-                        if (!parmatint_m)
-                            lossDs_m->addParticle(bunch->R[i], bunch->P[i], bunch->ID[i]);
-                        bunch->Bin[i] = -1;
-                        flagNeedUpdate = true;
-                    }
-                }
+            return true;
+        }
+    }
+    return false;
+}
+
+// rectangle collimators in cyclotron cylindrical coordinates
+// when there is no particlematterinteraction, the particle hitting collimator is deleted directly
+bool CCollimator::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {
+
+    bool flagNeedUpdate = false;
+    size_t tempnum = bunch->getLocalNum();
+    int pflag = 0;
+    // now check each particle in bunch
+    for (unsigned int i = 0; i < tempnum; ++i) {
+        if (bunch->PType[i] == ParticleType::REGULAR && bunch->R[i](2) < zend_m && bunch->R[i](2) > zstart_m ) {
+            // only now careful check in r
+            pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
+            /// bunch->Bin[i] != -1 makes sure the particle is not stored in more than one collimator
+            if ((pflag != 0) && (bunch->Bin[i] != -1)) {
+                if (!parmatint_m)
+                    lossDs_m->addParticle(bunch->R[i], bunch->P[i], bunch->ID[i]);
+                bunch->Bin[i] = -1;
+                flagNeedUpdate = true;
             }
         }
     }
+    return flagNeedUpdate;
+}
+
+bool CCollimator::doFinaliseCheck(PartBunchBase<double, 3> *bunch, bool flagNeedUpdate) {
     reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
     if (flagNeedUpdate && parmatint_m) {
+        Vector_t rmin, rmax;
+        bunch->get_bounds(rmin, rmax);
         std::pair<Vector_t, double> boundingSphere;
         boundingSphere.first = 0.5 * (rmax + rmin);
         boundingSphere.second = euclidean_norm(rmax - boundingSphere.first);
@@ -144,13 +155,10 @@ ElementBase::ElementType CCollimator::getType() const {
 }
 
 void CCollimator::doSetGeom() {
-    // calculate maximum and mininum r from these
-    // current implementation not perfect minimum does not need to lie on corner, or in middle
-    rmin_m = std::hypot(xstart_m, ystart_m);
+    // calculate maximum r from these
     rmax_m = -std::numeric_limits<double>::max();
     for (int i=0; i<4; i++) {
         double r = std::hypot(geom_m[i].x, geom_m[i].y);
-        rmin_m = std::min(rmin_m, r);
         rmax_m = std::max(rmax_m, r);
     }
     // some debug printout
@@ -158,19 +166,4 @@ void CCollimator::doSetGeom() {
     //     *gmsg << "point " << i << " ( " << geom_m[i].x << ", " << geom_m[i].y << ")" << endl;
     // }
     // *gmsg << "rmin " << rmin_m << " rmax " << rmax_m << endl;
-}
-
-// bool CCollimator::checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax) {
-
-//     double r1 = std::hypot(rmax(0), rmax(1));
-//     bool isDead = false;
-//     if (rmax(2) >= zstart_m && rmin(2) <= zend_m) {
-//         if ( r1 > rstart_m - 10.0 && r1 < rend_m + 10.0 ){
-//             if (r(2) < zend_m && r(2) > zstart_m ) {
-//                 int pflag = checkPoint(r(0), r(1));
-//                 isDead = (pflag != 0);
-//             }
-//         }
-//     }
-//     return isDead;
-// }
\ No newline at end of file
+}
\ No newline at end of file

src/Classic/AbsBeamline/CCollimator.h

@@ -57,6 +57,10 @@ private:
     virtual void doSetGeom() override;
     /// Virtual hook for finalise
     virtual void doFinalise() override;
+    /// Virtual hook for preCheck
+    virtual bool doPreCheck(PartBunchBase<double, 3>*) override;
+    /// Virtual hook for finaliseCheck
+    virtual bool doFinaliseCheck(PartBunchBase<double, 3> *bunch, bool flagNeedUpdate) override;
 
     bool informed_m = false; ///< Flag if error information already printed
 
@@ -65,7 +69,6 @@ private:
     double zend_m;
     double width_m;
     ///@}
-    double rmin_m; ///< minimum extend in r
     double rmax_m; ///< maximum extend in r
 
     ParticleMatterInteractionHandler *parmatint_m = nullptr;

src/Classic/AbsBeamline/PluginElement.cpp

@@ -95,6 +95,17 @@ void PluginElement::setDimensions(double xstart, double xend, double ystart, dou
     B_m = xstart_m - xend_m;
     R_m = std::sqrt(A_m*A_m+B_m*B_m);
     C_m = ystart_m*xend_m - xstart_m*yend_m;
+
+    // element equation: A*X + B*Y + C = 0
+    // point closest to origin https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
+    double x_close = 0.0;
+    if (R_m > 0.0)
+        x_close = - A_m * C_m / (R_m * R_m);
+
+    if (x_close > std::min(xstart_m, xend_m) && x_close < std::max(xstart_m, xend_m) )
+        rmin_m  = std::abs(C_m) / std::hypot(A_m,B_m);
+    else
+        rmin_m = rstart_m;
 }
 
 void PluginElement::setGeom(const double dist) {
@@ -126,6 +137,49 @@ void PluginElement::setGeom(const double dist) {
     doSetGeom();
 }
 
+void PluginElement::changeWidth(PartBunchBase<double, 3> *bunch, const double tstep) {
+
+    Vector_t meanP(0.0, 0.0, 0.0);
+    for(unsigned int i = 0; i < bunch->getLocalNum(); ++i) {
+        for(int d = 0; d < 3; ++d) {
+            meanP(d) += bunch->P[i](d);
+        }
+    }
+    reduce(meanP, meanP, OpAddAssign());
+    meanP = meanP / Vector_t(bunch->getTotalNum());
+
+    double stangle = calculateIncidentAngle(meanP(0), meanP(1));
+    constexpr double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
+    double lstep   = euclidean_norm(meanP) / Util::getGamma(meanP) * c_mmtns * tstep; // [mm]
+    double sWidth  = lstep / sqrt( 1 + 1/stangle/stangle );
+    setGeom(sWidth);
+}
+
+double PluginElement::calculateIncidentAngle(double xp, double yp) const {
+    double k1, k2, tangle = 0.0;
+    if ( B_m == 0.0 && xp == 0.0) {
+        // width is 0.0, keep non-zero
+        tangle = 0.1;
+    } else if ( B_m == 0.0 ){
+        k1 = yp / xp;
+        if (k1 == 0.0)
+            tangle = 1.0e12;
+        else
+            tangle = std::abs(1 / k1);
+    } else if ( xp == 0.0 ) {
+        k2 = - A_m/B_m;
+        if ( k2 == 0.0 )
+            tangle = 1.0e12;
+        else
+            tangle = std::abs(1 / k2);
+    } else {
+        k1 = yp / xp;
+        k2 = - A_m / B_m;
+        tangle = std::abs(( k1-k2 ) / (1 + k1*k2));
+    }
+    return tangle;
+}
+
 double PluginElement::getXStart() const {
     return xstart_m;
 }
@@ -144,7 +198,15 @@ double PluginElement::getYEnd() const {
 
 bool PluginElement::check(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) {
     bool flag = false;
-    flag = doCheck(bunch, turnnumber, t, tstep); // virtual hook
+    // check if bunch close
+    bool bunchClose = preCheck(bunch);
+
+    if (bunchClose == true) {
+        flag = doCheck(bunch, turnnumber, t, tstep); // virtual hook
+    }
+    // finalise, can have reduce
+    flag = finaliseCheck(bunch, flag);
+
     return flag;
 }
 
@@ -156,8 +218,8 @@ void PluginElement::getDimensions(double &zBegin, double &zEnd) const {
 int PluginElement::checkPoint(const double &x, const double &y) const {
     int    cn = 0;
     for(int i = 0; i < 4; i++) {
-        if(((geom_m[i].y <= y) && (geom_m[i+1].y > y))
-           || ((geom_m[i].y > y) && (geom_m[i+1].y <= y))) {
+        if((   (geom_m[i].y <= y) && (geom_m[i+1].y >  y))
+           || ((geom_m[i].y >  y) && (geom_m[i+1].y <= y))) {
 
             float vt = (float)(y - geom_m[i].y) / (geom_m[i+1].y - geom_m[i].y);
             if(x < geom_m[i].x + vt * (geom_m[i+1].x - geom_m[i].x))

src/Classic/AbsBeamline/PluginElement.h

@@ -69,14 +69,26 @@ public:
 protected:
     /// Sets geometry geom_m with element width dist
     void setGeom(const double dist);
+    /// Change probe width depending on step size and angle of particle
+    void changeWidth(PartBunchBase<double, 3> *bunch, const double tstep);
+    /// Calculate angle of particle/bunch wrt to element
+    double calculateIncidentAngle(double xp, double yp) const;
 
 private:
+    /// Check if bunch is close to element
+    bool preCheck(PartBunchBase<double, 3> *bunch) {return doPreCheck(bunch);}
+    /// Finalise call after check
+    bool finaliseCheck(PartBunchBase<double, 3> *bunch, bool flagNeedUpdate) {return doFinaliseCheck(bunch, flagNeedUpdate);}
     /// Pure virtual hook for initialise
-    virtual void doInitialise(PartBunchBase<double, 3> *bunch) {};
+    virtual void doInitialise(PartBunchBase<double, 3> *bunch) {}
     /// Pure virtual hook for check
     virtual bool doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) = 0;
     /// Virtual hook for setGeom
     virtual void doSetGeom() {};
+    /// Virtual hook for preCheck
+    virtual bool doPreCheck(PartBunchBase<double, 3>*) {return true;}
+    /// Virtual hook for finaliseCheck
+    virtual bool doFinaliseCheck(PartBunchBase<double, 3> *, bool flagNeedUpdate) {return flagNeedUpdate;}
     /// Virtual hook for finalise
     virtual void doFinalise() {};
     /// Virtual hook for goOffline
@@ -94,6 +106,7 @@ protected:
     double rstart_m;
     double rend_m;
     ///@}
+    double rmin_m;    ///< radius closest to the origin
     Point  geom_m[5]; ///< actual geometry positions with adaptive width such that each particle hits element once per turn
     double A_m, B_m, R_m, C_m; ///< Geometric lengths used in calculations
 

src/Classic/AbsBeamline/Probe.cpp

@@ -30,7 +30,7 @@ void Probe::accept(BeamlineVisitor &visitor) const {
 void Probe::doInitialise(PartBunchBase<double, 3> *bunch) {
     *gmsg << "* Initialize probe" << endl;
     bool singlemode = (bunch->getTotalNum() == 1) ? true : false;
-    peakfinder_m = std::unique_ptr<PeakFinder> (new PeakFinder(getOutputFN(), rstart_m, rend_m, step_m, singlemode));
+    peakfinder_m = std::unique_ptr<PeakFinder> (new PeakFinder(getOutputFN(), rmin_m, rend_m, step_m, singlemode));
 }
 
 void Probe::doGoOffline() {
@@ -48,7 +48,7 @@ double Probe::getStep() const {
     return step_m;
 }
 
-bool Probe::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) {
+bool Probe::doPreCheck(PartBunchBase<double, 3> *bunch) {
     Vector_t rmin, rmax;
     bunch->get_bounds(rmin, rmax);
     // interested in absolute minimum and maximum
@@ -59,90 +59,49 @@ bool Probe::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const
     double rbunch_min = std::hypot(xmin, ymin);
     double rbunch_max = std::hypot(xmax, ymax);
 
-    if( rbunch_max > rstart_m - 10.0 && rbunch_min < rend_m + 10.0 ) {
-        size_t tempnum = bunch->getLocalNum();
-
-        Vector_t meanP(0.0, 0.0, 0.0);
-        for(unsigned int i = 0; i < bunch->getLocalNum(); ++i) {
-            for(int d = 0; d < 3; ++d) {
-                meanP(d) += bunch->P[i](d);
-            }
-        }
-        reduce(meanP, meanP, OpAddAssign());
-        meanP = meanP / Vector_t(bunch->getTotalNum());
-
-        // change probe width depending on step size and angle of particle
-        double sk1, sk2, stangle = 0.0;
-        if ( B_m == 0.0 ){
-            sk1 = meanP(1)/meanP(0);
-            if (sk1 == 0.0)
-                stangle = 1.0e12;
-            else
-                stangle = std::abs(1/sk1);
-        } else if (meanP(0) == 0.0 ) {
-            sk2 = - A_m/B_m;
-            if ( sk2 == 0.0 )
-                stangle = 1.0e12;
-            else
-                stangle = std::abs(1/sk2);
-        } else {
-            sk1 = meanP(1)/meanP(0);
-            sk2 = - A_m/B_m;
-            stangle = std::abs(( sk1-sk2 )/(1 + sk1*sk2));
-        }
-        double lstep = (sqrt(1.0-1.0/(1.0+dot(meanP, meanP))) * Physics::c) * tstep*1.0e-6; // [mm]
-        double Swidth = lstep / sqrt( 1 + 1/stangle/stangle );
-        setGeom(Swidth);
-        Vector_t probepoint;
-
-        for(unsigned int i = 0; i < tempnum; ++i) {
-            int pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
-            if(pflag == 0) continue;
-            // calculate closest point at probe -> better to use momentum direction
-            // probe: y = -A/B * x - C/B
-            // perpendicular line through R: y = B/A * x + R(1) - B/A * R(0)
-            // probepoint(0) = (B_m*B_m*bunch->R[i](0) - A_m*B_m*bunch->R[i](1) - A_m*C_m) / (R_m*R_m);
-            // probepoint(1) = (A_m*A_m*bunch->R[i](1) - A_m*B_m*bunch->R[i](0) - B_m*C_m) / (R_m*R_m);
-            // probepoint(2) = bunch->R[i](2);
-            // calculate time correction for probepoint
-            // dist1 > 0, right hand, dt > 0; dist1 < 0, left hand, dt < 0
-            double dist1 = (A_m*bunch->R[i](0)+B_m*bunch->R[i](1)+C_m)/R_m/1000.0; // [m]
-            double k1, k2, tangle = 0.0;
-            if ( B_m == 0.0 ){
-                k1 = bunch->P[i](1)/bunch->P[i](0);
-                if (k1 == 0.0)
-                    tangle = 1.0e12;
-                else
-                    tangle = std::abs(1/k1);
-            } else if ( bunch->P[i](0) == 0.0 ) {
-                k2 = -A_m/B_m;
-                if (k2 == 0.0)
-                    tangle = 1.0e12;
-                else
-                    tangle = std::abs(1/k2);
-            } else {
-                k1 = bunch->P[i](1)/bunch->P[i](0);
-                k2 = -A_m/B_m;
-                tangle = std::abs(( k1-k2 )/(1 + k1*k2));
-            }
-            double dist2 = dist1 * sqrt( 1+1/tangle/tangle );
-            double dt = dist2/(sqrt(1.0-1.0/(1.0 + dot(bunch->P[i], bunch->P[i]))) * Physics::c)*1.0e9;
-
-            probepoint = bunch->R[i] + dist2 * 1000.0 * bunch->P[i] / euclidean_norm(bunch->P[i]);
-
-            // peak finder uses millimetre not metre
-            peakfinder_m->addParticle(probepoint);
-
-            /*FIXME Issue #45: mm --> m (when OPAL-Cycl uses metre insteas of millimetre,
-             * this can be removed.
-             */
-            probepoint *= 0.001;
-
-            lossDs_m->addParticle(probepoint, bunch->P[i], bunch->ID[i], t+dt, turnnumber);
-        }
-
-        peakfinder_m->evaluate(turnnumber);
+    if( rbunch_max > rmin_m - 10.0 && rbunch_min < rend_m + 10.0 ) {
+        return true;
     }
+    return false;
+}
+
+bool Probe::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) {
+    changeWidth(bunch, tstep);
+
+    size_t tempnum = bunch->getLocalNum();
+
+    Vector_t probepoint;
+
+    for(unsigned int i = 0; i < tempnum; ++i) {
+        int pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
+        if(pflag == 0) continue;
+        // calculate closest point at probe -> better to use momentum direction
+        // probe: y = -A/B * x - C/B or A*X + B*Y + C = 0
+        // perpendicular line through R: y = B/A * x + R(1) - B/A * R(0)
+        // probepoint(0) = (B_m*B_m*bunch->R[i](0) - A_m*B_m*bunch->R[i](1) - A_m*C_m) / (R_m*R_m);
+        // probepoint(1) = (A_m*A_m*bunch->R[i](1) - A_m*B_m*bunch->R[i](0) - B_m*C_m) / (R_m*R_m);
+        // probepoint(2) = bunch->R[i](2);
+        // calculate time correction for probepoint
+        // dist1 > 0, right hand, dt > 0; dist1 < 0, left hand, dt < 0
+        double dist1 = (A_m * bunch->R[i](0) + B_m * bunch->R[i](1) + C_m) / R_m * 1.0e-3; // [m]
+        double tangle = calculateIncidentAngle(bunch->P[i](0), bunch->P[i](1));
+        double dist2 = dist1 * sqrt(1.0 + 1.0 / tangle / tangle);
+        double dt = dist2 / (sqrt(1.0 - 1.0 / (1.0 + dot(bunch->P[i], bunch->P[i]))) * Physics::c) * 1.0e9;
+
+        probepoint = bunch->R[i] + dist2 * 1000.0 * bunch->P[i] / euclidean_norm(bunch->P[i]);
+
+        // peak finder uses millimetre not metre
+        peakfinder_m->addParticle(probepoint);
+
+        /*FIXME Issue #45: mm --> m (when OPAL-Cycl uses metre instead of millimetre,
+         * this can be removed.
+         */
+        probepoint *= 0.001;
+
+        lossDs_m->addParticle(probepoint, bunch->P[i], bunch->ID[i], t+dt, turnnumber);
+    }
+
+    peakfinder_m->evaluate(turnnumber);
 
     // we do not lose particles in the probe
     return false;

src/Classic/AbsBeamline/Probe.h

@@ -41,6 +41,8 @@ private:
     virtual bool doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) override;
     /// Hook for goOffline
     virtual void doGoOffline() override;
+    /// Virtual hook for preCheck
+    virtual bool doPreCheck(PartBunchBase<double, 3>*) override;
 
     double step_m; ///< Step size of the probe (bin width in histogram file)
     std::unique_ptr<PeakFinder> peakfinder_m; ///< Pointer to Peakfinder instance

src/Classic/AbsBeamline/Septum.cpp

@@ -16,13 +16,11 @@ Septum::Septum():Septum("")
 Septum::Septum(const std::string &name):
     PluginElement(name),
     width_m(0.0) {
-    setDimensions(0.0, 0.0, 0.0, 0.0);
 }
 
 Septum::Septum(const Septum &right):
     PluginElement(right),
     width_m(right.width_m) {
-    setDimensions(right.xstart_m, right.xend_m, right.ystart_m, right.yend_m);
     setGeom(width_m);
 }
 
@@ -52,9 +50,7 @@ void Septum::setWidth(double width) {
     setGeom(width_m);
 }
 
-bool Septum::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {
-
-    bool flag = false;
+bool Septum::doPreCheck(PartBunchBase<double, 3> *bunch) {
     Vector_t rmin;
     Vector_t rmax;
     bunch->get_bounds(rmin, rmax);
@@ -64,22 +60,36 @@ bool Septum::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/,
     double rbunch_max = std::hypot(xmax, ymax);
 
     if(rbunch_max > rstart_m - 100)  {
-        for(unsigned int i = 0; i < bunch->getLocalNum(); ++i) {
-            Vector_t R = bunch->R[i];
-            double slope = (yend_m - ystart_m) / (xend_m - xstart_m);
-            double intcept = ystart_m - slope * xstart_m;
-            double intcept1 = intcept - width_m / 2.0 * sqrt(slope * slope + 1);
-            double intcept2 = intcept + width_m / 2.0 * sqrt(slope * slope + 1);
-
-            double line1 = fabs(slope * R(0) + intcept1);
-            double line2 = fabs(slope * R(0) + intcept2);
-
-            if(fabs(R(1)) > line2 && fabs(R(1)) < line1 && R(0) > xstart_m && R(0) < xend_m && R(1) > ystart_m && R(1) < yend_m) {
-
-                bunch->lossDs_m->addParticle(R, bunch->P[i], bunch->ID[i]);
-                bunch->Bin[i] = -1;
-                flag = true;
-            }
+        return true;
+    }
+    return false;
+}
+
+bool Septum::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {
+
+    bool flag = false;
+    const double slope = (yend_m - ystart_m) / (xend_m - xstart_m);
+    const double halfLength = width_m / 2.0 * std::hypot(slope, 1);
+    const double intcept = ystart_m - slope * xstart_m;
+    const double intcept1 = intcept - halfLength;
+    const double intcept2 = intcept + halfLength;
+
+    for(unsigned int i = 0; i < bunch->getLocalNum(); ++i) {
+        const Vector_t& R = bunch->R[i];
+
+        double line1 = fabs(slope * R(0) + intcept1);
+        double line2 = fabs(slope * R(0) + intcept2);
+
+        if(fabs(R(1)) > line2 &&
+           fabs(R(1)) < line1 &&
+           R(0) > xstart_m    &&
+           R(0) < xend_m      &&
+           R(1) > ystart_m    &&
+           R(1) < yend_m) {
+
+            bunch->lossDs_m->addParticle(R, bunch->P[i], bunch->ID[i]);
+            bunch->Bin[i] = -1;
+            flag = true;
         }
     }
     return flag;

src/Classic/AbsBeamline/Septum.h

@@ -38,6 +38,8 @@ private:
     virtual void doInitialise(PartBunchBase<double, 3> *bunch) override;
     /// Record hits when bunch particles pass
     virtual bool doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) override;
+    /// Virtual hook for preCheck
+    virtual bool doPreCheck(PartBunchBase<double, 3>*) override;
 
     ///@{ input geometry positions
     double width_m;

src/Classic/AbsBeamline/Stripper.cpp

@@ -68,11 +68,7 @@ bool  Stripper::getStop () const {
     return stop_m;
 }
 
-//change the stripped particles to outcome particles
-bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) {
-
-    bool flagNeedUpdate = false;
-    bool flagresetMQ = false;
+bool Stripper::doPreCheck(PartBunchBase<double, 3> *bunch) {
     Vector_t rmin, rmax, strippoint;
     bunch->get_bounds(rmin, rmax);
     // interested in absolute maximum
@@ -80,118 +76,82 @@ bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, co
     double ymax = std::max(std::abs(rmin(1)), std::abs(rmax(1)));
     double rbunch_max = std::hypot(xmax, ymax);
 
-    if(rbunch_max > rstart_m - 10.0 ){
+    if(rbunch_max > rmin_m - 10.0) {
+        return true;
+    }
+    return false;
+}
 
-        size_t count = 0;
-        size_t tempnum = bunch->getLocalNum();
-        int pflag = 0;
+//change the stripped particles to outcome particles
+bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) {
+    changeWidth(bunch, tstep);
 
-        Vector_t meanP(0.0, 0.0, 0.0);
-        for(unsigned int i = 0; i < bunch->getLocalNum(); ++i) {
-          for(int d = 0; d < 3; ++d) {
-              meanP(d) += bunch->P[i](d);
-          }
-        }
-        reduce(meanP, meanP, OpAddAssign());
-        meanP = meanP / Vector_t(bunch->getTotalNum());
-
-        double sk1, sk2, stangle = 0.0;
-        if ( B_m == 0.0 ){
-            sk1 = meanP(1)/meanP(0);
-            if(sk1 == 0.0)
-                stangle =1.0e12;
-            else
-                stangle = std::abs(1/sk1);
-        }else if (meanP(0) == 0.0 ){
-            sk2 = - A_m/B_m;
-            if