PluginElement.cpp

Go to the documentation of this file.

//
// Class PluginElement
//   Abstract Interface for (Cyclotron) Plugin Elements (CCollimator, Probe, Stripper, Septum)
//   Implementation via Non-Virtual Interface Template Method
//
// Copyright (c) 2018-2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "AbsBeamline/PluginElement.h"
 
#include "AbsBeamline/BeamlineVisitor.h"
#include "PartBunch/PartBunch.h"
#include "Physics/Physics.h"
#include "Physics/Units.h"
#include "Structure/LossDataSink.h"
#include "Utilities/Options.h"
#include "Utilities/Util.h"
 
PluginElement::PluginElement() : PluginElement("") {}
 
PluginElement::PluginElement(const std::string& name) : Component(name) {
    setDimensions(0.0, 0.0, 0.0, 0.0);
}
 
PluginElement::PluginElement(const PluginElement& right) : Component(right) {
    setDimensions(right.xstart_m, right.xend_m, right.ystart_m, right.yend_m);
}
 
PluginElement::~PluginElement() {
    if (online_m) goOffline();
}
 
void PluginElement::initialise(PartBunch_t* bunch, double&, double&) { initialise(bunch); }
 
void PluginElement::initialise(PartBunch_t* bunch) {
    RefPartBunch_m = bunch;
    lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getOutputFN(), !Options::asciidump));
    // virtual hook
    doInitialise(bunch);
    goOnline(-1e6);
}
 
void PluginElement::finalise() {
    doFinalise();
    if (online_m) goOffline();
}
 
void PluginElement::goOffline() {
    if (online_m && lossDs_m) lossDs_m->save();
    lossDs_m.reset(nullptr);
    doGoOffline();
    online_m = false;
}
 
bool PluginElement::bends() const { return false; }
 
bool PluginElement::apply(const std::shared_ptr<ParticleContainer_t>& /*pc*/) { return false; }
 
bool PluginElement::apply(
        const size_t& /*i*/, const double&, Vector_t<double, 3>&, Vector_t<double, 3>&) {
    return false;
}
 
bool PluginElement::apply(
        const Vector_t<double, 3>& /*R*/, const Vector_t<double, 3>& /*P*/, const double& /*t*/,
        Vector_t<double, 3>& /*E*/, Vector_t<double, 3>& /*B*/) {
    *gmsg << "passed R, P, t, E, B arguments not used in PluginElement::apply" << endl;
    return false;
}
 
bool PluginElement::applyToReferenceParticle(
        const Vector_t<double, 3>&, const Vector_t<double, 3>&, const double&, Vector_t<double, 3>&,
        Vector_t<double, 3>&) {
    return false;
}
 
void PluginElement::setDimensions(double xstart, double xend, double ystart, double yend) {
    xstart_m = xstart;
    ystart_m = ystart;
    xend_m   = xend;
    yend_m   = yend;
    rstart_m = std::hypot(xstart, ystart);
    rend_m   = std::hypot(xend, yend);
    // start position is the one with lowest radius
    if (rstart_m > rend_m) {
        std::swap(xstart_m, xend_m);
        std::swap(ystart_m, yend_m);
        std::swap(rstart_m, rend_m);
    }
    A_m = yend_m - ystart_m;
    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) {
    double slope;
    if (xend_m == xstart_m)
        slope = 1.0e12;
    else
        slope = (yend_m - ystart_m) / (xend_m - xstart_m);
 
    double coeff2   = std::sqrt(1 + slope * slope);
    double coeff1   = slope / coeff2;
    double halfdist = dist / 2.0;
    geom_m[0].x     = xstart_m - halfdist * coeff1;
    geom_m[0].y     = ystart_m + halfdist / coeff2;
 
    geom_m[1].x = xstart_m + halfdist * coeff1;
    geom_m[1].y = ystart_m - halfdist / coeff2;
 
    geom_m[2].x = xend_m + halfdist * coeff1;
    geom_m[2].y = yend_m - halfdist / coeff2;
 
    geom_m[3].x = xend_m - halfdist * coeff1;
    geom_m[3].y = yend_m + halfdist / coeff2;
 
    geom_m[4].x = geom_m[0].x;
    geom_m[4].y = geom_m[0].y;
 
    doSetGeom();
}
 
void PluginElement::changeWidth(
        PartBunch_t* bunch, int i, const double tstep, const double tangle) {
    constexpr double c_mtns = Physics::c / Units::s2ns;  // m/s --> m/ns
 
    const double tmp = std::sqrt(dot(bunch->P(i), bunch->P(i)));
    double lstep     = tmp / Util::getGamma(bunch->P(i)) * c_mtns * tstep;  // [m]
    double sWidth    = lstep / std::sqrt(1 + 1 / tangle / tangle);
    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; }
 
double PluginElement::getXEnd() const { return xend_m; }
 
double PluginElement::getYStart() const { return ystart_m; }
 
double PluginElement::getYEnd() const { return yend_m; }
 
bool PluginElement::check(
        PartBunch_t* bunch, const int turnnumber, const double t, const double tstep) {
    bool flag = false;
    // 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;
}
 
void PluginElement::getFieldExtend(double& zBegin, double& zEnd) const {
    zBegin = -0.005;
    zEnd   = 0.005;
}
 
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))) {
            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)) ++cn;
        }
    }
    return (cn & 1);  // 0 if even (out), and 1 if odd (in)
}
 
void PluginElement::save() {
    OpalData::OpenMode openMode;
    if (numPassages_m > 0) {
        openMode = OpalData::OpenMode::APPEND;
    } else {
        openMode = OpalData::getInstance()->getOpenMode();
    }
    lossDs_m->save(1, openMode);
    numPassages_m++;
}