diff --git a/src/Algorithms/ParallelCyclotronTracker.cpp b/src/Algorithms/ParallelCyclotronTracker.cpp
index 19fb8f4e08c520a91ed2d4ab060436abb87dc13e..015d6e4ce3134107ae0dc09eaf615e18bd6806ae 100644
--- a/src/Algorithms/ParallelCyclotronTracker.cpp
+++ b/src/Algorithms/ParallelCyclotronTracker.cpp
@@ -41,6 +41,7 @@
#include "AbsBeamline/MultipoleTCurvedVarRadius.h"
#include "AbsBeamline/MultipoleTStraight.h"
#include "AbsBeamline/Offset.h"
+#include "AbsBeamline/OutputPlane.h"
#include "AbsBeamline/PluginElement.h"
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
@@ -670,6 +671,37 @@ void ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius(const MultipoleTCu
myElements.push_back(dynamic_cast<MultipoleTCurvedVarRadius*>(multTvcurv.clone()));
}
+
+/**
+ *
+ *
+ * @param plane
+ */
+void ParallelCyclotronTracker::visitOutputPlane(const OutputPlane& plane) {
+ if (opalRing_m == nullptr) {
+ throw OpalException("ParallelCylcotronTracker::visitOutputPlane",
+ "Attempt to place an OutputPlane when Ring not defined");
+ }
+
+ OutputPlane* elptr = dynamic_cast<OutputPlane*>(plane.clone());
+ elptr->setGlobalFieldMap(opalRing_m);
+ myElements.push_back(elptr);
+ elptr->initialise(itsBunch_m);
+
+ double BcParameter[8] = {};
+
+ Vector_t centre = elptr->getCentre();
+ Vector_t norm = elptr->getNormal();
+ double width = elptr->getHorizontalExtent();
+ BcParameter[0] = centre[0]-width*norm[1];
+ BcParameter[1] = centre[0]+width*norm[1];
+ BcParameter[2] = centre[1]-width*norm[0];
+ BcParameter[3] = centre[1]+width*norm[0];
+ BcParameter[4] = Units::mm2m; // thickness, not used
+ buildupFieldList(BcParameter, ElementType::OUTPUTPLANE, elptr);
+}
+
+
/**
*
*
@@ -1173,6 +1205,7 @@ void ParallelCyclotronTracker::execute() {
CavityCrossData ccd = {cav, cav->getSinAzimuth(), cav->getCosAzimuth(), cav->getPerpenDistance()};
cavCrossDatas_m.push_back(ccd);
} else if ( type == ElementType::CCOLLIMATOR ||
+ type == ElementType::OUTPUTPLANE ||
type == ElementType::PROBE ||
type == ElementType::SEPTUM ||
type == ElementType::STRIPPER) {
diff --git a/src/Algorithms/ParallelCyclotronTracker.h b/src/Algorithms/ParallelCyclotronTracker.h
index a338d90310ca7589b776227565e8b07cd1a68d2b..cd1a6580ddc4dfc79496edecaf3f2b9f8c01c9d0 100644
--- a/src/Algorithms/ParallelCyclotronTracker.h
+++ b/src/Algorithms/ParallelCyclotronTracker.h
@@ -143,6 +143,9 @@ public:
/// Apply the algorithm to a offset (placement).
virtual void visitOffset(const Offset&);
+ /// Apply the algorithm to a outputplane.
+ virtual void visitOutputPlane(const OutputPlane&);
+
/// Apply the algorithm to a probe.
virtual void visitProbe(const Probe&);
diff --git a/src/Classic/AbsBeamline/BeamlineVisitor.h b/src/Classic/AbsBeamline/BeamlineVisitor.h
index d56131f916a97ce7375c4c4f612dfa8def75b4b2..2eaed611f25d3f29ae2f6211ca007532b14b59cd 100644
--- a/src/Classic/AbsBeamline/BeamlineVisitor.h
+++ b/src/Classic/AbsBeamline/BeamlineVisitor.h
@@ -55,6 +55,7 @@ class MultipoleTStraight;
class MultipoleTCurvedConstRadius;
class MultipoleTCurvedVarRadius;
class Offset;
+class OutputPlane;
class Probe;
class RBend;
class RBend3D;
@@ -136,7 +137,10 @@ public:
/// Apply the algorithm to an offset (placement).
virtual void visitOffset(const Offset &) = 0;
-
+
+ /// Apply the algorithm to an outputplane.
+ virtual void visitOutputPlane(const OutputPlane &) = 0;
+
/// Apply the algorithm to a probe.
virtual void visitProbe(const Probe &) = 0;
diff --git a/src/Classic/AbsBeamline/CMakeLists.txt b/src/Classic/AbsBeamline/CMakeLists.txt
index 4f30c20d974799cfaf3a9c935f20b3ffa87032a6..46509618ef7f494a1b4788840ad35832be06af7d 100644
--- a/src/Classic/AbsBeamline/CMakeLists.txt
+++ b/src/Classic/AbsBeamline/CMakeLists.txt
@@ -23,6 +23,7 @@ set (_SRCS
MultipoleTCurvedConstRadius.cpp
MultipoleTCurvedVarRadius.cpp
Offset.cpp
+ OutputPlane.cpp
PluginElement.cpp
Probe.cpp
RBend.cpp
@@ -65,6 +66,7 @@ set (HDRS
MultipoleTCurvedConstRadius.h
MultipoleTCurvedVarRadius.h
Offset.h
+ OutputPlane.h
PluginElement.h
Probe.h
RBend.h
diff --git a/src/Classic/AbsBeamline/ElementBase.cpp b/src/Classic/AbsBeamline/ElementBase.cpp
index 08e40dcf069f1b2d4e404111806d0f2cec6e3f00..869d23628c5e76d5f562b5e2836d64e0de75fa68 100644
--- a/src/Classic/AbsBeamline/ElementBase.cpp
+++ b/src/Classic/AbsBeamline/ElementBase.cpp
@@ -85,6 +85,7 @@ const std::map<ElementType, std::string> ElementBase::elementTypeToString_s = {
{ElementType::MULTIPOLE, "Multipole"},
{ElementType::MULTIPOLET, "MultipoleT"},
{ElementType::OFFSET, "Offset"},
+ {ElementType::OUTPUTPLANE, "OutputPlane"},
{ElementType::PROBE, "Probe"},
{ElementType::RBEND, "RBend"},
{ElementType::RBEND3D, "RBend3D"},
diff --git a/src/Classic/AbsBeamline/ElementBase.h b/src/Classic/AbsBeamline/ElementBase.h
index 498bf8470e52d9b9e596fdc20f35e9b00b90fd15..10f1774305cf6f1e459ebbdb9db6e7e80c8dae12 100644
--- a/src/Classic/AbsBeamline/ElementBase.h
+++ b/src/Classic/AbsBeamline/ElementBase.h
@@ -100,6 +100,7 @@ enum class ElementType: unsigned short {
MULTIPOLE,
MULTIPOLET,
OFFSET,
+ OUTPUTPLANE,
PROBE,
RBEND,
RBEND3D,
diff --git a/src/Classic/AbsBeamline/MultipoleT.cpp b/src/Classic/AbsBeamline/MultipoleT.cpp
index 0d47f805bf4c997056a3dfb14ac41e893cea25e5..81f9390760dd63eb797b3c1fe087688bb91aafd3 100644
--- a/src/Classic/AbsBeamline/MultipoleT.cpp
+++ b/src/Classic/AbsBeamline/MultipoleT.cpp
@@ -1,6 +1,7 @@
/*
* Copyright (c) 2017, Titus Dascalu
* Copyright (c) 2018, Martin Duy Tat
+ * Copyright (c) 2019-2023, Chris Rogers
* All rights reserved.
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
diff --git a/src/Classic/AbsBeamline/OutputPlane.cpp b/src/Classic/AbsBeamline/OutputPlane.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..90e8d056a42a2876369ff2a4502e995321db3d01
--- /dev/null
+++ b/src/Classic/AbsBeamline/OutputPlane.cpp
@@ -0,0 +1,322 @@
+// Copyright (c) 2023, Chris Rogers
+// 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/OutputPlane.h"
+
+#include "AbsBeamline/BeamlineVisitor.h"
+#include "Algorithms/PartBunchBase.h"
+#include "Physics/Physics.h"
+#include "Structure/LossDataSink.h"
+#include "Physics/Units.h"
+#include "Utilities/GeneralClassicException.h"
+
+extern Inform *gmsg;
+
+OutputPlane::OutputPlane() : OutputPlane("")
+{}
+
+OutputPlane::OutputPlane(const std::string& name):
+ PluginElement(name),
+ maxIterations_m(10),
+ algorithm_m(algorithm::INTERPOLATION),
+ verbose_m(0) {
+}
+
+OutputPlane::OutputPlane(const OutputPlane& right):
+ PluginElement(right),
+ maxIterations_m(10) {
+
+ field_m = right.field_m;
+ normal_m = right.normal_m;
+ centre_m = right.centre_m;
+ radialExtent_m = right.radialExtent_m;
+ verticalExtent_m = right.verticalExtent_m;
+ horizontalExtent_m = right.horizontalExtent_m;
+ maxIterations_m = right.maxIterations_m;
+ tolerance_m = right.tolerance_m;
+ nullfield_m = right.nullfield_m;
+ geom_m = right.geom_m;
+ recentre_m = right.recentre_m;
+ algorithm_m = right.algorithm_m;
+ verbose_m = right.verbose_m;
+}
+
+OutputPlane::~OutputPlane()
+{}
+
+ElementBase* OutputPlane::clone() const {
+ ElementBase* element = dynamic_cast<ElementBase*>(new OutputPlane(*this));
+ return element;
+}
+
+void OutputPlane::accept(BeamlineVisitor& visitor) const {
+ visitor.visitOutputPlane(*this);
+}
+
+void OutputPlane::doInitialise(PartBunchBase<double, 3>* /*bunch*/) {
+ //void OutputPlane::doInitialise() {
+ *gmsg << "* Initialize OutputPlane at " << centre_m << " with normal " << normal_m << "\n*";
+ if (radialExtent_m >= 0) {
+ *gmsg << " Radial extent " << radialExtent_m;
+ }
+ if (horizontalExtent_m >= 0) {
+ *gmsg << " Horizontal extent " << horizontalExtent_m;
+ }
+ if (verticalExtent_m >= 0) {
+ *gmsg << " Vertical extent " << verticalExtent_m;
+ }
+ if (field_m != nullptr) {
+ *gmsg << " Using field map " << field_m->getName();
+ } else {
+ *gmsg << " Using empty field map";
+ }
+ if (recentre_m >= 0) {
+ *gmsg << " Recentre event " << recentre_m;
+ }
+ *gmsg << endl;
+}
+
+void OutputPlane::doGoOffline() {
+ *gmsg << "* OutputPlane goes offline " << getName() << endl;
+}
+
+void OutputPlane::RK4Step(
+ const double& tStep,
+ const double& chargeToMass,
+ const double& t,
+ Vector_t& R1,
+ Vector_t& P1) const {
+ double thalfStep = tStep / 2.;
+ double tPlusHalf = t + thalfStep;
+ double tPlusStep = t + tStep;
+ // f = dy/dt
+ // f1 = f(x,t)
+ double deriv1[6];
+
+ getDerivatives(R1, P1, t, chargeToMass, deriv1);
+ // f2 = f(x+dt*f1/2, t+dt/2 )
+ double deriv2[6];
+ Vector_t R2(R1[0] + thalfStep * deriv1[0], R1[1] + thalfStep * deriv1[1], R1[2] + thalfStep * deriv1[2]);
+ Vector_t P2(P1[0] + thalfStep * deriv1[3], P1[1] + thalfStep * deriv1[4], P1[2] + thalfStep * deriv1[5]);
+ getDerivatives(R2, P2, tPlusHalf, chargeToMass, deriv2);
+
+ // f3 = f( x+dt*f2/2, t+dt/2 )
+ double deriv3[6];
+ Vector_t R3(R1[0] + thalfStep * deriv2[0], R1[1] + thalfStep * deriv2[1], R1[2] + thalfStep * deriv2[2]);
+ Vector_t P3(P1[0] + thalfStep * deriv2[3], P1[1] + thalfStep * deriv2[4], P1[2] + thalfStep * deriv2[5]);
+ getDerivatives(R3, P3, tPlusHalf, chargeToMass, deriv3);
+
+ // f4 = f(x+dt*f3, t+dt ).
+ double deriv4[6];
+ Vector_t R4(R1[0] + tStep * deriv3[0], R1[1] + tStep * deriv3[1], R1[2] + tStep * deriv3[2]);
+ Vector_t P4(P1[0] + tStep * deriv3[3], P1[1] + tStep * deriv3[4], P1[2] + tStep * deriv3[5]);
+ getDerivatives(R4, P4, tPlusStep, chargeToMass, deriv4);
+
+ // Return x(t+dt) computed from fourth-order R-K.
+ R1[0] += (deriv1[0] + deriv4[0] + 2. * (deriv2[0] + deriv3[0])) * tStep / 6.;
+ R1[1] += (deriv1[1] + deriv4[1] + 2. * (deriv2[1] + deriv3[1])) * tStep / 6.;
+ R1[2] += (deriv1[2] + deriv4[2] + 2. * (deriv2[2] + deriv3[2])) * tStep / 6.;
+ P1[0] += (deriv1[3] + deriv4[3] + 2. * (deriv2[3] + deriv3[3])) * tStep / 6.;
+ P1[1] += (deriv1[4] + deriv4[4] + 2. * (deriv2[4] + deriv3[4])) * tStep / 6.;
+ P1[2] += (deriv1[5] + deriv4[5] + 2. * (deriv2[5] + deriv3[5])) * tStep / 6.;
+}
+
+void OutputPlane::getDerivatives(const Vector_t& R,
+ const Vector_t& P,
+ const double& t,
+ const double& chargeToMass,
+ double* yp) const {
+
+ double gamma = std::sqrt(1 + (P[0] * P[0] + P[1] * P[1] + P[2] * P[2]));
+
+ Vector_t beta = P / gamma;
+ double betax = beta[0];
+ double betay = beta[1];
+ double betaz = beta[2];
+
+ yp[0] = Physics::c * betax;
+ yp[1] = Physics::c * betay;
+ yp[2] = Physics::c * betaz;
+
+ if (field_m == nullptr) {
+ throw GeneralClassicException("OutputPlane::getDerivatives",
+ "Field was null");
+ }
+ Vector_t externalB, externalE;
+ field_m->apply(R, P, t, externalE, externalB);
+ //double kiloGaussToTesla = 0.1;
+ externalB *= Units::kG2T;
+ externalE *= Units::kV2V / Units::mm2m / Physics::c;
+
+ yp[3] = chargeToMass * (externalB(2) * betay - externalB(1) * betaz+externalE(0));
+ yp[4] = chargeToMass * (externalB(0) * betaz - externalB(2) * betax + externalE(1));
+ yp[5] = chargeToMass * (externalB(1) * betax - externalB(0) * betay + externalE(2));
+}
+
+bool OutputPlane::checkOne(const int index, const double tstep, double chargeToMass,
+ double& t, Vector_t& R, Vector_t& P) {
+
+ // distance from particle to the output plane
+ // time units are ns
+ Vector_t delta = R-centre_m;
+ double distance = dot(normal_m, delta);
+
+ // maximum step, assuming no curvature, rough guess for relativistic beta
+ double betaEstimate = euclidean_norm(P);
+ if (betaEstimate > 1) {
+ betaEstimate = 1.0;
+ }
+ double sStep = tstep*betaEstimate*Physics::c;
+ if (verbose_m > 3) {
+ *gmsg << "* First check crossing of plane " << getName() << " at " << centre_m << " with normal " << normal_m << endl;
+ *gmsg << " Particle " << index << " with R " << R << " P " << P << " t0 " << t << " tstep " << tstep << endl;
+ *gmsg << " Distance prestep " << distance << " compared to s step length " << sStep << endl;
+ }
+ if (std::abs(distance) > sStep) {
+ // we can't cross the plane
+ return false;
+ }
+
+ Vector_t rTest(R);
+ Vector_t pTest(P);
+ RK4Step(tstep, chargeToMass, t, rTest, pTest);
+ double distanceTest = dot(normal_m, (rTest-centre_m));
+ if (verbose_m > 2) {
+ *gmsg << "* Second check crossing of plane " << getName() << " at " << centre_m << " with normal " << normal_m << endl;
+ *gmsg << " Particle " << index << " with R " << R << " P " << P << " tstep " << tstep << endl;
+ *gmsg << " After RK4 " << rTest << " " << pTest << endl;
+ *gmsg << " Step distance " << distanceTest << " compared to " << distance << endl;
+ }
+
+ if (distance != 0 && distanceTest/distance > 0) {
+ // step does not cross the plane - give up
+ // note that particle could cross and cross back the plane in a single
+ // time-step; in this case the particle is not registered
+ return false;
+ }
+
+ if (algorithm_m == algorithm::RK4STEP) {
+ rk4Test(tstep, chargeToMass, t, R, P);
+ } else if (algorithm_m == algorithm::INTERPOLATION) {
+ interpolation(t, R, P);
+ }
+
+ delta = R-centre_m;
+ if (verbose_m > 1) {
+ *gmsg << "* Track estimate RK4? " << bool(algorithm_m == algorithm::RK4STEP)
+ << " INTERPOLATION? " << bool(algorithm_m == algorithm::INTERPOLATION)
+ << " R " << R << " P " << P << " t " << t << endl
+ << " delta " << delta << endl;
+ }
+ if (horizontalExtent_m > 0 && delta[0] * delta[0] + delta[1] * delta[1] > horizontalExtent_m*horizontalExtent_m) {
+ // out of horizontal extent. Again, defined in global coordinates
+ return false;
+ }
+ if (verticalExtent_m > 0 && abs(delta[2]) > verticalExtent_m) {
+ // out of vertical extent; note that this is defined in global coordinates
+ // not in the coordinates of the plane
+ return false;
+ }
+ if (radialExtent_m > 0 &&
+ delta[0] * delta[0] + delta[1] * delta[1] + delta[2] * delta[2] > radialExtent_m * radialExtent_m) {
+ // out of radial extent
+ return false;
+ }
+ if (recentre_m == index) {
+ recentre(R, P);
+ *gmsg << "* Recentred output plane to " << centre_m
+ << " with normal " << normal_m << " by event " << index << endl;
+ }
+ if (verbose_m > 0) {
+ *gmsg << "* Found track" << endl;
+ }
+ return true;
+}
+
+void OutputPlane::rk4Test(double tstep, double chargeToMass,
+ double& t, Vector_t& R, Vector_t& P) {
+ double preStepDistance = 0.0;
+ double postStepDistance = dot(normal_m, (R-centre_m));
+ size_t iteration = 0;
+ while (std::abs(postStepDistance) > tolerance_m && iteration < maxIterations_m) {
+ preStepDistance = postStepDistance;
+ RK4Step(tstep, chargeToMass, t, R, P); // this updates R and P
+ if (verbose_m > 2) {
+ *gmsg << " RK4 iteration " << iteration << " step distance " << preStepDistance << " R " << R << " P " << P
+ << " centre " << centre_m << " d: " << R-centre_m << " t: " << t << " dt: " << tstep << endl;
+ Vector_t externalB, externalE;
+ field_m->apply(R, P, t, externalE, externalB);
+ *gmsg << " B " << externalB << " [kG] E " << externalE << " [MV/m] " << endl;
+ }
+ t += tstep;
+ postStepDistance = dot(normal_m, (R - centre_m));
+ if (postStepDistance/preStepDistance < 0) { // straddling the plane
+ // we stepped too far; step in opposite direction
+ // step length in ratio of distance to plane
+ tstep *= -abs(postStepDistance)/abs(postStepDistance - preStepDistance);
+ } else {
+ // we didn't step far enough; step in same direction
+ tstep *= abs(postStepDistance)/abs(postStepDistance - preStepDistance);
+ }
+ iteration++;
+ }
+}
+
+void OutputPlane::interpolation(double& t, Vector_t& R, Vector_t& P) {
+ // trajectory R = R0 + V dt
+ // plane = (X-X0).n
+ // intersection time t0 = (X0-R0).N / (V.N)
+ // intersection position = R0 + V dt
+ // relativistic gamma:
+ double gamma = std::sqrt(1 + P[0] * P[0] + P[1] * P[1] + P[2] * P[2]);
+ Vector_t velocity = P/gamma*Physics::c; // m/ns
+ double dt = dot((centre_m-R), normal_m) / dot(velocity, normal_m);
+ R += velocity*dt;
+ t += dt;
+}
+
+
+bool OutputPlane::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber,
+ const double t, const double tstep) {
+ size_t tempnum = bunch->getLocalNum();
+ for(unsigned int i = 0; i < tempnum; ++i) {
+ if (verbose_m > 2) {
+ *gmsg << "OutputPlane checking at time " << t
+ << " turn number " << turnnumber << " track id " << i << endl;
+ }
+ Vector_t R(bunch->R[i]);
+ Vector_t P(bunch->P[i]);
+ double t0(t);
+ double chargeToMass = bunch->Q[i] / Physics::q_e * Physics::c * Physics::c
+ / bunch->M[i]/Units::GeV2eV; // electron charge cancels
+ bool crossing = checkOne(i, tstep, chargeToMass, t0, R, P);
+ if (crossing && lossDs_m) {
+ nHits_m += 1;
+ lossDs_m->addParticle(OpalParticle(bunch->ID[i], R, P,
+ t0, bunch->Q[i], bunch->M[i]),
+ std::make_pair(turnnumber, bunch->bunchNum[i]));
+ }
+ }
+ return false;
+}
+
+void OutputPlane::recentre(Vector_t R, Vector_t P) {
+ setCentre(R);
+ setNormal(P);
+ recentre_m = -1; // don't recentre again
+}
+
+ElementType OutputPlane::getType() const {
+ return ElementType::OUTPUTPLANE;
+}
diff --git a/src/Classic/AbsBeamline/OutputPlane.h b/src/Classic/AbsBeamline/OutputPlane.h
new file mode 100644
index 0000000000000000000000000000000000000000..28bc5a87d9c2768d232fc45e8b219a82e5e0f8f3
--- /dev/null
+++ b/src/Classic/AbsBeamline/OutputPlane.h
@@ -0,0 +1,357 @@
+// Copyright (c) 2023, Chris Rogers
+// 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/>.
+//
+
+#ifndef CLASSIC_OutputPlane_HH
+#define CLASSIC_OutputPlane_HH
+
+#include "AbsBeamline/PluginElement.h"
+#include "BeamlineGeometry/StraightGeometry.h"
+#include "Fields/NullField.h"
+#include "Steppers/Steppers.h"
+#include "Steppers/RK4.h"
+
+#include <string>
+
+class Component;
+
+/** Class OutputPlane
+ * ------------------------------------------------------------------------
+ * Class OutputPlane defines an output plane, which enables output of tracking
+ * data at an arbitrary spatial plane. The crossing of the plane in a given time
+ * step is detected by (1) checking if the step is within c*dt of the plane and
+ * then (2) checking if an RK4 step crosses the plane.
+ *
+ * If the track crosses, then we estimate the position of the track either by
+ * doing a linear interpolation between the initial position and the position
+ * of the test step; OR we do an iterative RK4 step, reducing the RK4 step size
+ * each time until we find crossing position to some tolerance or the maximum
+ * number of iterations is reached.
+ *
+ * Additional features - we can recentre the plane on a particular track (for
+ * example, one might seek to make a plane perpendicular to a reference
+ * trajectory); we can limit the extent of the plane so only particles crossing
+ * within a certain distance are registered.
+ *
+ * BUG: if the "Global" Opal tracking is not RK4, then there is a small
+ * possibility that the "Global" track and the OutputPlane step will be not
+ * perfectly aligned; tracks could register crossing twice or not crossing at
+ * all. An improved algorithm would be to use the "global" stepper or make some
+ * hack to consider e.g. two steps.
+ */
+
+class OutputPlane: public PluginElement {
+
+public:
+ enum class algorithm {RK4STEP, INTERPOLATION};
+
+ /** Constructor with given name. */
+ explicit OutputPlane(const std::string& name);
+
+ /** Default constructor (initialise to empty */
+ OutputPlane();
+
+ /** Copy constructor */
+ OutputPlane(const OutputPlane&);
+
+ /** Destructor */
+ virtual ~OutputPlane();
+
+ /** Inheritable copy operation */
+ ElementBase* clone() const override;
+
+ /** Apply visitor to output plane. */
+ virtual void accept(BeamlineVisitor&) const override;
+
+ /** Get the field map */
+ inline Component* getGlobalFieldMap() const;
+ /** Set the field map */
+ inline void setGlobalFieldMap(Component* field);
+
+ /** Get the normal to the plane */
+ inline Vector_t getNormal() const;
+ /** Set the normal to the plane */
+ inline void setNormal(Vector_t normal);
+
+ /** Get the centre of the plane */
+ inline Vector_t getCentre() const;
+ /** Set the centre of the plane */
+ inline void setCentre(Vector_t centre);
+
+ /** Get the tolerance, used when finding intersection with the plane */
+ inline double getTolerance() const;
+ /** Set the tolerance, used when finding intersection with the plane */
+ inline void setTolerance(double tolerance);
+
+ /** Get the maximum allowed iteratiosn when finding intersection */
+ inline size_t getMaxIterations() const;
+ /** Set the maximum allowed iteratiosn when finding intersection */
+ inline void setMaxIterations(size_t max);
+
+ /** Return the index of the recentring particle*/
+ inline int getRecentre() const;
+ /** Set the index of the recentring particle*/
+ inline void setRecentre(int willRecentre);
+
+ /** Get the radial extent from the plane centre */
+ inline double getRadialExtent() const;
+ /** Set the radial extent from the centre */
+ inline void setRadialExtent(double r);
+
+ /** Get the vertical extent from the plane centre
+ */
+ inline double getVerticalExtent() const;
+ /** Set the vertical extent from the plane centre */
+ inline void setVerticalExtent(double z);
+
+ /** Get the horizontal extent from the plane centre
+ */
+ inline double getHorizontalExtent() const;
+ /** Set the horizontal extent from the plane centre */
+ inline void setHorizontalExtent(double width);
+
+ /** Get the algorithm used to interpolate from step to the plane */
+ inline algorithm getAlgorithm() const;
+ /** Set the algorithm used to interpolate from step to the plane */
+ inline void setAlgorithm(algorithm alg);
+
+ /** Get the horizontal extent from the plane centre
+ */
+ inline int getVerboseLevel() const;
+ /** Set the verbose level
+ *
+ * set to 0 for silent running; 1 when hit detected;
+ * 2 when second check passes; 3 when first check passes;
+ * 4 to output every step
+ **/
+ inline void setVerboseLevel(int verbose);
+
+ /** Check for plane crossing. If crossed fill t, R, P with the intercept.
+ *
+ * @li index: Particle number of the particle crossing
+ * @li tstep: Time step used by the global tracking
+ * @li chargeToMass: charge to mass ratio of the particle, in units of
+ * (positron charge)/GeV/c^2
+ * @li t: time of the initial particle [s]. Unchanged if no crossing is
+ * found, else filled with the estimated time of the particle at the
+ * plane crossing point.
+ * @li R: position three-vector of the initial particle [m]. Unchanged if
+ * no crossing is found, else filled with the estimated position of
+ * the particle at the plane crossing point.
+ * @li P: (gamma beta) three-vector of the initial particle [].
+ * Unchanged if no crossing is found, else filled with the estimated
+ * position of the particle at the plane crossing point.
+ *
+ * The routine checks for crossing by performing a full RK4 step using the
+ * user-supplied time step. The distance to the plane is calculated. If the
+ * distance changes sign, the particle is assumed to cross; if the distance
+ * does not change sign, the particle is assumed to stay on the same side.
+ * Note that particles that step across the boundary and back again will
+ * not be registered as a crossing.
+ *
+ * The intercept is found by iteratively doing RK4 steps. Successive step
+ * sizes are found by linear interpolation based on the distance from the
+ * plane and the time step. I could use the Steppers/RK4.h routines but
+ * these automatically update the global particle information and I don't
+ * want to do that, I just want a local stepping.
+ *
+ * If recentre_m is positive then centre_m is set to the position of
+ * the n^th crossing particle and normal_m is set to be parallel to the
+ * momentum vector. recentre_m is set to -1. The idea is to generate
+ * output planes in the coordinate system of some trajectory (e.g. for
+ * studying beam transport/transfer matrix). Normally this will be the
+ * reference trajectory but it is possible to make it another trajectory if
+ * required.
+ *
+ * If radialExtent_m is > 0, crossings more than radialExtent_m from the
+ * centre are ignored. If horizontalExtent_m is > 0, crossings more than
+ * horizontalExtent_m from the centre, in global x-y plane, are ignored. If
+ * verticalExtent_m is > 0, crossings more than verticalExtent_m in
+ * direction parallel to global z axis, are ignored.
+ *
+ * @returns True if the particle crosses the plane; else false.
+ */
+ bool checkOne(const int index, const double tstep, double chargeToMass,
+ double& t, Vector_t& R, Vector_t& P);
+
+ /** Returns empty field */
+ NullField& getField() {return nullfield_m;}
+ /** Returns empty field */
+ const NullField& getField() const {return nullfield_m;}
+
+ /** Returns empty geometry */
+ StraightGeometry& getGeometry() {return geom_m;}
+ /** Returns empty geometry */
+ const StraightGeometry& getGeometry() const {return geom_m;}
+ /** Make an RK4Step
+ * @li tstep: time step [s]
+ * @li chargeToMass: chargeToMass ratio [m^2/s^2/GV]
+ * @li t: t [s]
+ * @li R: R [m]
+ * @li P: beta gamma [] - note dimensionless
+ */
+ void RK4Step(const double& tstep,
+ const double& chargeToMass,
+ const double& t,
+ Vector_t& R,
+ Vector_t& P) const;
+
+ void operator=(const OutputPlane&) = delete;
+
+ //ElementBase::ElementType getType() const;
+ ElementType getType() const;
+ void recentre(Vector_t R, Vector_t P);
+
+
+private:
+ /// Initialise peakfinder file
+ virtual void doInitialise(PartBunchBase<double, 3>* /*bunch*/);
+
+ /// Record probe hits when bunch particles pass
+ inline bool doPreCheck(PartBunchBase<double, 3> *bunch) override;
+
+ /// Record probe hits when bunch particles pass
+ virtual bool doCheck(PartBunchBase<double, 3> *bunch,
+ const int turnnumber,
+ const double t,
+ const double tstep) override;
+ /// Hook for goOffline
+ virtual void doGoOffline() override;
+ void getDerivatives(const Vector_t& R,
+ const Vector_t& P,
+ const double& t,
+ const double& chargeToMass,
+ double* yp) const;
+ double distanceToPlane(Vector_t point) const;
+ void rk4Test(double tstep, double chargeToMass,
+ double& t, Vector_t& R, Vector_t& P);
+ void interpolation(double& t, Vector_t& R, Vector_t& P);
+ bool getFieldsAtPoint(const Vector_t& R, const Vector_t& P, const double& t, Vector_t& Efield, Vector_t& Bfield);
+
+ Component* field_m = NULL; // field map - this is a borrowed pointer
+ Vector_t normal_m; // normal to the output plane
+ Vector_t centre_m; // centre of the output plane
+ double radialExtent_m = -1.0; // maximum radial extent of the plane (circular)
+ double verticalExtent_m = -1.0; // maximum vertical extent of the plane (rectangular)
+ double horizontalExtent_m = -1.0; // maximum horizontal extent of the plane (rectangular)
+ size_t maxIterations_m = 10; // maximum number of iterations when finding intercept
+ double tolerance_m = 1e-9; // tolerance on distance from plane when finding intercept
+ NullField nullfield_m; // dummy variable for inheritance
+ StraightGeometry geom_m; // dummy variable for inheritance
+ algorithm algorithm_m; // sets interpolation or RK4
+ int recentre_m = -1; // particle index for recentring
+ int nHits_m = 0; // counter for number of hits on the plane
+ int verbose_m = 0; // verbosity
+
+ typedef std::function<bool(const double&,
+ const size_t&,
+ Vector_t&,
+ Vector_t&)> function_t;
+ std::unique_ptr<Stepper<function_t> > stepper_m;
+};
+
+Component* OutputPlane::getGlobalFieldMap() const {
+ return field_m;
+}
+
+void OutputPlane::setGlobalFieldMap(Component* field) {
+ field_m = field;
+}
+
+Vector_t OutputPlane::getNormal() const {
+ return normal_m;
+}
+
+void OutputPlane::setNormal(Vector_t normal) {
+ normal_m = normal;
+ normal_m /= euclidean_norm(normal);
+}
+
+Vector_t OutputPlane::getCentre() const {
+ return centre_m;
+}
+
+void OutputPlane::setCentre(Vector_t centre) {
+ centre_m = centre;
+}
+
+double OutputPlane::getTolerance() const {
+ return tolerance_m;
+}
+
+void OutputPlane::setTolerance(double tolerance) {
+ tolerance_m = tolerance;
+}
+
+size_t OutputPlane::getMaxIterations() const {
+ return maxIterations_m;
+}
+
+void OutputPlane::setMaxIterations(size_t max) {
+ maxIterations_m = max;
+}
+
+double OutputPlane::getHorizontalExtent() const {
+ return horizontalExtent_m;
+}
+
+void OutputPlane::setHorizontalExtent(double width) {
+ horizontalExtent_m = width;
+}
+
+double OutputPlane::getVerticalExtent() const {
+ return verticalExtent_m;
+}
+
+void OutputPlane::setVerticalExtent(double width) {
+ verticalExtent_m = width;
+}
+
+double OutputPlane::getRadialExtent() const {
+ return radialExtent_m;
+}
+
+void OutputPlane::setRadialExtent(double radius) {
+ radialExtent_m = radius;
+}
+
+int OutputPlane::getRecentre() const {
+ return recentre_m;
+}
+
+void OutputPlane::setRecentre(int recentre) {
+ recentre_m = recentre;
+}
+
+OutputPlane::algorithm OutputPlane::getAlgorithm() const {
+ return algorithm_m;
+}
+
+void OutputPlane::setAlgorithm(OutputPlane::algorithm alg) {
+ algorithm_m = alg;
+}
+
+bool OutputPlane::doPreCheck(PartBunchBase<double, 3>* /*bunch*/) {
+ return true;
+}
+
+int OutputPlane::getVerboseLevel() const {
+ return verbose_m;
+}
+
+void OutputPlane::setVerboseLevel(int verbose) {
+ verbose_m = verbose;
+}
+
+
+#endif // CLASSIC_OutputPlane_HH
diff --git a/src/Classic/AbsBeamline/SpecificElementVisitor.h b/src/Classic/AbsBeamline/SpecificElementVisitor.h
index 1d01a43c95709c78cea7cd4b0d084b2513de9dc2..b9dd50d4c8c7798d3a4b911de4d12a0720445f76 100644
--- a/src/Classic/AbsBeamline/SpecificElementVisitor.h
+++ b/src/Classic/AbsBeamline/SpecificElementVisitor.h
@@ -36,6 +36,7 @@
#include "AbsBeamline/MultipoleTStraight.h"
#include "AbsBeamline/MultipoleTCurvedConstRadius.h"
#include "AbsBeamline/MultipoleTCurvedVarRadius.h"
+#include "AbsBeamline/OutputPlane.h"
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/Ring.h"
@@ -136,6 +137,9 @@ public:
/// Apply the algorithm to an offset (placement).
virtual void visitOffset(const Offset &);
+ /// Apply the algorithm to an output plane.
+ virtual void visitOutputPlane(const OutputPlane &);
+
/// Apply the algorithm to a probe.
virtual void visitProbe(const Probe &prob);
@@ -314,6 +318,11 @@ void SpecificElementVisitor<ELEM>::visitOffset(const Offset &element) {
CastsTrait<ELEM, Offset>::apply(allElementsOfTypeE, element);
}
+template<class ELEM>
+void SpecificElementVisitor<ELEM>::visitOutputPlane(const OutputPlane &element) {
+ CastsTrait<ELEM, OutputPlane>::apply(allElementsOfTypeE, element);
+}
+
template<class ELEM>
void SpecificElementVisitor<ELEM>::visitProbe(const Probe &element) {
CastsTrait<ELEM, Probe>::apply(allElementsOfTypeE, element);
diff --git a/src/Classic/Algorithms/DefaultVisitor.cpp b/src/Classic/Algorithms/DefaultVisitor.cpp
index ddc6f79207833c0274a9b846f5a7f371d1990296..f6c36aa6ec5f4d240602d76c77c7344a4fc23160 100644
--- a/src/Classic/Algorithms/DefaultVisitor.cpp
+++ b/src/Classic/Algorithms/DefaultVisitor.cpp
@@ -37,6 +37,7 @@
#include "AbsBeamline/MultipoleTCurvedVarRadius.h"
#include "AbsBeamline/MultipoleTStraight.h"
#include "AbsBeamline/Offset.h"
+#include "AbsBeamline/OutputPlane.h"
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RBend3D.h"
@@ -144,6 +145,10 @@ void DefaultVisitor::visitOffset(const Offset& off) {
applyDefault(off);
}
+void DefaultVisitor::visitOutputPlane(const OutputPlane& out) {
+ applyDefault(out);
+}
+
void DefaultVisitor::visitProbe(const Probe &probe) {
applyDefault(probe);
}
diff --git a/src/Classic/Algorithms/DefaultVisitor.h b/src/Classic/Algorithms/DefaultVisitor.h
index 94fbe608526eefb7a33bfc89a065411308ac56ed..0f8e82d32be4e3b13416c991fd11b900431c7fbf 100644
--- a/src/Classic/Algorithms/DefaultVisitor.h
+++ b/src/Classic/Algorithms/DefaultVisitor.h
@@ -96,6 +96,9 @@ public:
/// Apply the algorithm to an offset (placement).
virtual void visitOffset(const Offset &);
+ /// Apply the algorithm to an output plane.
+ virtual void visitOutputPlane(const OutputPlane &);
+
/// Apply the algorithm to a probe.
virtual void visitProbe(const Probe &prob);
diff --git a/src/Elements/CMakeLists.txt b/src/Elements/CMakeLists.txt
index 4874230ae30c4ea125ada726fd5d14e37b4e7c75..f83f37b73c7d87fc278148c6b2553f09fe456676 100644
--- a/src/Elements/CMakeLists.txt
+++ b/src/Elements/CMakeLists.txt
@@ -25,6 +25,7 @@ set (_SRCS
OpalOffset/OpalLocalCartesianOffset.cpp
OpalOffset/OpalGlobalCylindricalOffset.cpp
OpalOffset/OpalGlobalCartesianOffset.cpp
+ OpalOutputPlane.cpp
OpalProbe.cpp
OpalPepperPot.cpp
OpalPolynomialTimeDependence.cpp
diff --git a/src/Elements/OpalOutputPlane.cpp b/src/Elements/OpalOutputPlane.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..da90c7fbb041a3a330adb54ae41226dd602ea303
--- /dev/null
+++ b/src/Elements/OpalOutputPlane.cpp
@@ -0,0 +1,174 @@
+// Copyright (c) 2023, Chris Rogers
+// 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 "Elements/OpalOutputPlane.h"
+#include "Classic/AbsBeamline/OutputPlane.h"
+#include "AbstractObjects/Attribute.h"
+#include "Attributes/Attributes.h"
+#include "Utilities/OpalException.h"
+
+const std::string docstring =
+ std::string("The \"OUTPUTPLANE\" element writes out position at which")+
+ std::string("trajectories cross a given plane.");
+
+OpalOutputPlane::OpalOutputPlane():
+ OpalElement(SIZE, "OUTPUTPLANE", docstring.c_str()) {
+
+ itsAttr[CENTRE] = Attributes::makeRealArray
+ ("CENTRE", "3-vector position of the plane centre [m]");
+ itsAttr[NORMAL] = Attributes::makeRealArray
+ ("NORMAL", "3-vector normal to the plane");
+ itsAttr[REFERENCE_ALIGNMENT_PARTICLE] = Attributes::makeReal
+ ("REFERENCE_ALIGNMENT_PARTICLE",
+ "Set to a particle number (usually 0, the reference particle). "
+ "The first time that the particle crosses the reference plane, then "
+ "the plane will be moved to centre on that particle and point in S of the particle.", -1);
+ itsAttr[TOLERANCE] = Attributes::makeReal
+ ("TOLERANCE", "Tolerance on position of track intercept [m]", 1e-6);
+ itsAttr[WIDTH] = Attributes::makeReal
+ ("WIDTH", "Full width of the output plane [m], defined in the lab coordinate system (*not* the output plane coordinate system).");
+ itsAttr[HEIGHT] = Attributes::makeReal
+ ("HEIGHT", "Full height of the output plane [m], defined in the lab coordinate system (*not* the output plane coordinate system)");
+ itsAttr[RADIUS] = Attributes::makeReal
+ ("RADIUS", "Maximum distance from centre of plane for crossings [m].");
+ itsAttr[ALGORITHM] = Attributes::makePredefinedString
+ ("ALGORITHM",
+ "The algorithm used to step from the track point to the plane", {"INTERPOLATION", "RK4"}, "RK4");
+ itsAttr[XSTART] = Attributes::makeReal
+ ("XSTART", "Define a plane with horizontal extent [m] from XSTART to XEND and vertical extent from YSTART to YEND", 0.0);
+ itsAttr[YSTART] = Attributes::makeReal
+ ("YSTART", "Define a plane with horizontal extent [m] from XSTART to XEND and vertical extent from YSTART to YEND", 1.0);
+ itsAttr[XEND] = Attributes::makeReal
+ ("XEND", "Define a plane with horizontal extent [m] from XSTART to XEND and vertical extent from YSTART to YEND", 0.0);
+ itsAttr[YEND] = Attributes::makeReal
+ ("YEND", "Define a plane with horizontal extent [m] from XSTART to XEND and vertical extent from YSTART to YEND", 0.0);
+
+ itsAttr[PLACEMENT_STYLE] = Attributes::makePredefinedString("PLACEMENT_STYLE",
+ "Set to PROBE to define the plane using XSTART, XEND, YSTART, YEND or CENTRE_NORMAL to define the plane using centre and normal", {"CENTRE_NORMAL", "PROBE"}, "PROBE");
+ itsAttr[VERBOSE] = Attributes::makeReal("VERBOSE",
+ "Set to 0 for minimal output up to 4 to output diagnostics on every track step. Output is sent to OPAL", 0);
+ registerOwnership();
+ setElement(new OutputPlane("OUTPUTPLANE"));
+}
+
+
+OpalOutputPlane::OpalOutputPlane(const std::string &name, OpalOutputPlane *parent):
+ OpalElement(name, parent) {
+ OutputPlane* plane = dynamic_cast<OutputPlane*>(parent->getElement());
+ setElement(new OutputPlane(*plane));
+}
+
+
+OpalOutputPlane::~OpalOutputPlane() {
+}
+
+
+OpalOutputPlane *OpalOutputPlane::clone(const std::string &name) {
+ return new OpalOutputPlane(name, this);
+}
+
+void OpalOutputPlane::update() {
+ OpalElement::update();
+ std::string placementStyle = Attributes::getString(itsAttr[PLACEMENT_STYLE]);
+ OutputPlane *output = dynamic_cast<OutputPlane *>(getElement());
+ Vector_t centre;
+ Vector_t normal;
+ double tolerance = Attributes::getReal(itsAttr[TOLERANCE]);
+ output->setTolerance(tolerance);
+
+ if (placementStyle == "CENTRE_NORMAL") {
+ std::vector<double> centreVec = Attributes::getRealArray(itsAttr[CENTRE]);
+ if (centreVec.size() != 3) {
+ throw OpalException("OpalOutputPlane::update()",
+ "OutputPlane centre should be a 3-vector");
+ }
+ centre = Vector_t(centreVec[0], centreVec[1], centreVec[2]);
+ std::vector<double> normalVec = Attributes::getRealArray(itsAttr[NORMAL]);
+ if (normalVec.size() != 3) {
+ throw OpalException("OpalOutputPlane::update()",
+ "OutputPlane normal should be a 3-vector");
+ }
+ normal = Vector_t(normalVec[0], normalVec[1], normalVec[2]);
+ if (euclidean_norm(normal) < tolerance) {
+ throw OpalException("OpalOutputPlane::update()",
+ "OutputPlane normal was not defined or almost zero");
+ }
+ } else if (placementStyle == "PROBE") {
+ centre = Vector_t(
+ Attributes::getReal(itsAttr[XSTART])+Attributes::getReal(itsAttr[XEND]),
+ Attributes::getReal(itsAttr[YSTART])+Attributes::getReal(itsAttr[YEND]),
+ 0.0
+ );
+ centre *= 0.5; // average
+ normal = Vector_t(
+ Attributes::getReal(itsAttr[YSTART])-Attributes::getReal(itsAttr[YEND]),
+ Attributes::getReal(itsAttr[XEND])-Attributes::getReal(itsAttr[XSTART]),
+ 0.0
+ );
+ double width = std::sqrt(normal[0]*normal[0]+normal[1]*normal[1]);
+ if (width < tolerance) {
+ throw OpalException("OpalOutputPlane::update()",
+ "OutputPlane had very small width or size was not defined");
+ }
+ output->setHorizontalExtent(width/2.0);
+ }
+ output->setCentre(centre);
+ output->setNormal(normal);
+
+
+ if (itsAttr[WIDTH]) {
+ double width = Attributes::getReal(itsAttr[WIDTH]);
+ if (width < 0) {
+ throw OpalException("OpalOutputPlane::update()",
+ "OutputPlane had negative width");
+ }
+ output->setHorizontalExtent(width/2.0);
+ }
+ if (itsAttr[HEIGHT]) {
+ double height = Attributes::getReal(itsAttr[HEIGHT]);
+ if (height < 0) {
+ throw OpalException("OpalOutputPlane::update()",
+ "OutputPlane had negative height");
+ }
+ output->setVerticalExtent(height/2.0);
+ }
+ if (itsAttr[RADIUS]) {
+ double radius = Attributes::getReal(itsAttr[RADIUS]);
+ if (radius < 0) {
+ throw OpalException("OpalOutputPlane::update()",
+ "OutputPlane had negative radius");
+ }
+ output->setRadialExtent(radius);
+ }
+
+ if (itsAttr[REFERENCE_ALIGNMENT_PARTICLE]) {
+ double ref = Attributes::getReal(itsAttr[REFERENCE_ALIGNMENT_PARTICLE]);
+ int refAlign = std::floor(ref + 0.5);
+ if (refAlign >= 0) { // disabled if negative
+ output->setRecentre(refAlign);
+ }
+ } else {
+ output->setRecentre(-1);
+ }
+ std::string algorithm = Attributes::getString(itsAttr[ALGORITHM]);
+ if (algorithm == "RK4") {
+ output->setAlgorithm(OutputPlane::algorithm::RK4STEP);
+ } else if (algorithm == "INTERPOLATION") {
+ output->setAlgorithm(OutputPlane::algorithm::INTERPOLATION);
+ }
+ output->setVerboseLevel(Attributes::getReal(itsAttr[VERBOSE]));
+ output->setOutputFN(Attributes::getString(itsAttr[OUTFN]));
+ // Transmit "unknown" attributes.
+ OpalElement::updateUnknown(output);
+}
\ No newline at end of file
diff --git a/src/Elements/OpalOutputPlane.h b/src/Elements/OpalOutputPlane.h
new file mode 100644
index 0000000000000000000000000000000000000000..15b6bf6a337e5c0f7537a78f3aa74421becfe546
--- /dev/null
+++ b/src/Elements/OpalOutputPlane.h
@@ -0,0 +1,75 @@
+// Copyright (c) 2023, Chris Rogers
+// 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/>.
+//
+
+#ifndef OPAL_OpalOutputPlane_H
+#define OPAL_OpalOutputPlane_H
+
+#include "Elements/OpalElement.h"
+
+// Class OpalOutputPlane
+// ------------------------------------------------------------------------
+/// Interface for output plane element.
+
+
+class OpalOutputPlane: public OpalElement {
+
+public:
+
+ /// The attributes of class OpalOutputPlane.
+ enum {
+ CENTRE = COMMON, // plane centre position (3-vector)
+ NORMAL, // plane normal (3-vector)
+ WIDTH, // maximum allowed horizontal displacement from centre, in global x-y plane
+ HEIGHT, // maximum allowed vertical displacement from centre, in global vertical direction
+ RADIUS, // maximum allowed radial displacement from centre
+
+ XSTART, // "PROBE"-style placement - horizontal start position
+ XEND, // "PROBE"-style placement - horizontal end position
+ YSTART, // "PROBE"-style placement - vertical start position
+ YEND, // "PROBE"-style placement - vertical end position
+
+ PLACEMENT_STYLE, // "PROBE" or "CENTRE_NORMAL"
+
+ ALGORITHM, // algorithm used to find crossing-point
+ TOLERANCE, // tolerance on position estimate
+ REFERENCE_ALIGNMENT_PARTICLE, // if true, centres on reference particle
+ VERBOSE, // set to 0 - 4 to get more verbose output
+
+ SIZE
+ };
+ /// Exemplar constructor.
+ OpalOutputPlane();
+
+ virtual ~OpalOutputPlane();
+
+ /// Make clone.
+ virtual OpalOutputPlane *clone(const std::string &name);
+
+ /// Fill in all registered attributes.
+ //virtual void fillRegisteredAttributes(const ElementBase &, ValueFlag);
+
+ /// Update the embedded CLASSIC septum.
+ virtual void update();
+
+private:
+
+ // Not implemented.
+ OpalOutputPlane(const OpalOutputPlane &);
+ void operator=(const OpalOutputPlane &);
+
+ // Clone constructor.
+ OpalOutputPlane(const std::string &name, OpalOutputPlane *parent);
+};
+
+#endif // OPAL_OpalOutputPlane_H
diff --git a/src/OpalConfigure/Configure.cpp b/src/OpalConfigure/Configure.cpp
index f1d19cb055ab2d57bffb5ed317a0dbdaf6887433..74bc38b353509c1edb91e72f3618586ca66e33c8 100644
--- a/src/OpalConfigure/Configure.cpp
+++ b/src/OpalConfigure/Configure.cpp
@@ -89,6 +89,7 @@
#include "Elements/OpalOffset/OpalLocalCylindricalOffset.h"
#include "Elements/OpalOffset/OpalGlobalCartesianOffset.h"
#include "Elements/OpalOffset/OpalGlobalCylindricalOffset.h"
+#include "Elements/OpalOutputPlane.h"
#include "Elements/OpalPepperPot.h"
#include "Elements/OpalPolynomialTimeDependence.h"
#include "Elements/OpalProbe.h"
@@ -212,6 +213,7 @@ namespace {
// opal->create(new OpalOffset::OpalLocalCylindricalOffset());
// opal->create(new OpalOffset::OpalGlobalCartesianOffset());
// opal->create(new OpalOffset::OpalGlobalCylindricalOffset());
+ opal->create(new OpalOutputPlane());
opal->create(new OpalPepperPot());
opal->create(new OpalPolynomialTimeDependence());
opal->create(new OpalProbe());
diff --git a/src/PyOpal/PyElements/CMakeLists.txt b/src/PyOpal/PyElements/CMakeLists.txt
index 5796f9e418abd600ac7b49350b0819e35b1b00dc..a8888e3de0f295f30d31ab2fcca08b3d2ef4d77e 100644
--- a/src/PyOpal/PyElements/CMakeLists.txt
+++ b/src/PyOpal/PyElements/CMakeLists.txt
@@ -38,3 +38,6 @@ python_module(variable_rf_cavity_fringe_field "${VARIABLERFCAVITYFRINGEFIELD_SRC
set (POLYNOMIALTIMEDEPENDENCE_SRCS PyPolynomialTimeDependence.cpp $<TARGET_OBJECTS:libPyOPALPyCore>)
python_module(polynomial_time_dependence "${POLYNOMIALTIMEDEPENDENCE_SRCS}" "elements")
+
+set (OUTPUTPLANE_SRCS PyOutputPlane.cpp $<TARGET_OBJECTS:libPyOPALPyCore>)
+python_module(output_plane "${OUTPUTPLANE_SRCS}" "elements")
diff --git a/src/PyOpal/PyElements/PyOutputPlane.cpp b/src/PyOpal/PyElements/PyOutputPlane.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ff963db4ac3a67fcac98906139040c2ab544dadf
--- /dev/null
+++ b/src/PyOpal/PyElements/PyOutputPlane.cpp
@@ -0,0 +1,58 @@
+// Copyright (c) 2023, Chris Rogers
+// 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 "PyOpal/PyCore/ExceptionTranslation.h"
+#include "PyOpal/PyCore/PyOpalObject.h"
+#include "PyOpal/PyCore/Globals.h"
+
+#include "Elements/OpalOutputPlane.h"
+
+namespace PyOpal {
+namespace PyOutputPlane {
+
+const char* module_docstring =
+ "output_plane contains the OutputPlane class";
+
+template <>
+std::vector<PyOpalObjectNS::AttributeDef> PyOpalObjectNS::PyOpalObject<OpalOutputPlane>::attributes = {
+ {"CENTRE", "centre", "", PyOpalObjectNS::FLOAT_LIST},
+ {"NORMAL", "normal", "", PyOpalObjectNS::FLOAT_LIST},
+ {"XSTART", "x_start", "", PyOpalObjectNS::DOUBLE},
+ {"XEND", "x_end", "", PyOpalObjectNS::DOUBLE},
+ {"YSTART", "y_start", "", PyOpalObjectNS::DOUBLE},
+ {"YEND", "y_end", "", PyOpalObjectNS::DOUBLE},
+ {"PLACEMENT_STYLE", "placement_style", "", PyOpalObjectNS::PREDEFINED_STRING},
+ {"ALGORITHM", "algorithm", "", PyOpalObjectNS::PREDEFINED_STRING},
+ {"TOLERANCE", "tolerance", "", PyOpalObjectNS::DOUBLE},
+ {"REFERENCE_ALIGNMENT_PARTICLE", "reference_alignment_particle", "", PyOpalObjectNS::INT},
+ {"OUTFN", "output_filename", "", PyOpalObjectNS::STRING}, // OUTFN comes from OpalElement (yes, all elements can have a filename!)
+ {"VERBOSE", "verbose_level", "", PyOpalObjectNS::INT},
+ {"WIDTH", "width", "", PyOpalObjectNS::DOUBLE},
+ {"HEIGHT", "height", "", PyOpalObjectNS::DOUBLE},
+ {"RADIUS", "radius", "", PyOpalObjectNS::DOUBLE}
+};
+
+template <>
+std::string PyOpalObjectNS::PyOpalObject<OpalOutputPlane>::classDocstring =
+"OutputPlane is used to generate output data based on particle tracks crossing a plane.";
+
+BOOST_PYTHON_MODULE(output_plane) {
+ PyOpal::Globals::Initialise();
+ ExceptionTranslation::registerExceptions();
+ PyOpalObjectNS::PyOpalObject<OpalOutputPlane> element;
+ auto elementClass = element.make_element_class("OutputPlane");
+}
+
+}
+}
diff --git a/src/ValueDefinitions/StringConstant.cpp b/src/ValueDefinitions/StringConstant.cpp
index 9ec2897ff0d8be0db3548fde6c4e7e1ba70f9a2e..ab1a9a6c26fef115051633b4bf844012d7bfb65a 100644
--- a/src/ValueDefinitions/StringConstant.cpp
+++ b/src/ValueDefinitions/StringConstant.cpp
@@ -125,6 +125,14 @@ StringConstant::StringConstant():
CREATE_STRINGCONSTANT("AC");
CREATE_STRINGCONSTANT("DC");
+ // OutputPlane / PLACEMENT_STYLE
+ CREATE_STRINGCONSTANT("CENTRE_NORMAL");
+ CREATE_STRINGCONSTANT("PROBE");
+
+ // OutputPlane / ALGORITHM
+ CREATE_STRINGCONSTANT("INTERPOLATION");
+ // CREATE_STRINGCONSTANT("RK4");
+
// ParticleMatterInteraction / TYPE
CREATE_STRINGCONSTANT("SCATTERING");
CREATE_STRINGCONSTANT("BEAMSTRIPPING");
diff --git a/tests/classic_src/AbsBeamline/CMakeLists.txt b/tests/classic_src/AbsBeamline/CMakeLists.txt
index 06c6230bdef5073b6ea3a51798e4f98f61a9b013..35ec3b6c5bd452b5ffeee12d194426e5c14f90ed 100644
--- a/tests/classic_src/AbsBeamline/CMakeLists.txt
+++ b/tests/classic_src/AbsBeamline/CMakeLists.txt
@@ -4,6 +4,7 @@ set (_SRCS
DipoleFieldTest.cpp
MultipoleTTest.cpp
OffsetTest.cpp
+ OutputPlaneTest.cpp
PolynomialTest.cpp
RingTest.cpp
SBend3DTest.cpp
diff --git a/tests/classic_src/AbsBeamline/OutputPlaneTest.cpp b/tests/classic_src/AbsBeamline/OutputPlaneTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a0a6b80e389ca94795995146dbca97dc5dad33fe
--- /dev/null
+++ b/tests/classic_src/AbsBeamline/OutputPlaneTest.cpp
@@ -0,0 +1,352 @@
+// Copyright (c) 2023, Chris Rogers
+// 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 <cmath>
+#include <fstream>
+#include <sstream>
+
+#include "gtest/gtest.h"
+#include "opal_test_utilities/SilenceTest.h"
+
+#include "Classic/Physics/Physics.h"
+#include "Classic/AbsBeamline/OutputPlane.h"
+#include "Classic/AbsBeamline/Component.h"
+
+
+class MockDipole : public Component {
+public:
+ MockDipole(const std::string /*&name*/) : bfield_m(0, 1, 0), efield_m() {}
+ NullField &getField() {return nullfield_m;}
+ const NullField &getField() const {return nullfield_m;}
+ void accept(BeamlineVisitor& /*&visitor*/) const {}
+ void initialise(PartBunchBase<double, 3>* /*bunch*/, double& /*start*/, double& /*end*/) {}
+ void finalise() {}
+ bool bends() const {return true;}
+ void getDimensions(double& /*zBegin*/, double& /*zEnd*/) const {}
+ StraightGeometry& getGeometry() {return geom_m;}
+ const StraightGeometry& getGeometry() const {return geom_m;}
+ ElementBase* clone() const {return NULL;}
+
+ void setField(Vector_t bfield, Vector_t efield) {bfield_m = bfield; efield_m = efield;}
+
+ bool apply(const Vector_t& /*R*/, const Vector_t& /*P*/, const double& /*t*/, Vector_t& E, Vector_t& B) override {
+ B = bfield_m;
+ E = efield_m;
+ return false;
+ }
+private:
+ NullField nullfield_m;
+ StraightGeometry geom_m;
+ Vector_t bfield_m;
+ Vector_t efield_m;
+};
+
+// Set also - lots of testing done in pyopal level
+class OutputPlaneTest : public ::testing::Test {
+public:
+ OutputPlaneTest() : dipole_m("") {
+ }
+
+ void SetUp( ) {
+ output_m.reset(new OutputPlane());
+ output_m->setRecentre(-1);
+
+ }
+
+ void TearDown( ) {
+ }
+
+ ~OutputPlaneTest() {
+ output_m.release();
+ }
+
+ std::unique_ptr<OutputPlane> output_m;
+ MockDipole dipole_m;
+ const double q2m = 1/Physics::c/Physics::c/Physics::m_p;
+
+ // OpalTestUtilities::SilenceTest silencer_m;
+};
+
+TEST_F(OutputPlaneTest, TestSetGet) {
+ // test basic set and get methods
+ output_m->setGlobalFieldMap(&dipole_m);
+ EXPECT_EQ(output_m->getGlobalFieldMap(), &dipole_m);
+
+ output_m->setNormal(Vector_t(1.0, 2.0, 3.0));
+ Vector_t ref = Vector_t(1.0, 2.0, 3.0)/std::sqrt(14);
+ EXPECT_NEAR(output_m->getNormal()[0], ref[0], 1e-9);
+ EXPECT_NEAR(output_m->getNormal()[1], ref[1], 1e-9);
+ EXPECT_NEAR(output_m->getNormal()[2], ref[2], 1e-9);
+
+ output_m->setCentre(Vector_t(4.0, 5.0, 6.0));
+ EXPECT_EQ(output_m->getCentre(), Vector_t(4.0, 5.0, 6.0));
+
+ output_m->setTolerance(1e-6);
+ EXPECT_EQ(output_m->getTolerance(), 1e-6);
+}
+
+
+TEST_F(OutputPlaneTest, TestCheckOne_FieldOff) {
+ //checkOne(const double tstep, const double chargeToMass,
+ // double& t, Vector_t& R, Vector_t& P)
+ output_m->setGlobalFieldMap(&dipole_m);
+ output_m->setCentre(Vector_t(0.0, 0.0, 1.0));
+ output_m->setNormal(Vector_t(1.0, 0.0, 1.0));
+ output_m->setTolerance(1e-9);
+
+ Vector_t R, P;
+ bool crossing;
+ double t = 0.0;
+
+ R = Vector_t(1.0, 0.0, 1.0);
+ P = Vector_t(0.1, 0.0, 0.0);
+ // c = 0.3 m/ns and particle is 1.0 m away; so should fail without doing RK4
+ // step at all
+ crossing = output_m->checkOne(0, 1.0e-9, 0.0, t, R, P);
+ EXPECT_EQ(R, Vector_t(1.0, 0.0, 1.0));
+ EXPECT_EQ(P, Vector_t(0.1, 0.0, 0.0));
+ EXPECT_FALSE(crossing);
+ // betagamma = 0.1 and particle is ~0.07 m away; so should fail after trying
+ // a single RK4 step (which goes about 0.03 m)
+ R = Vector_t(1-0.05, 0.0, -0.05);
+ crossing = output_m->checkOne(0, 1.0e-9, 0.0, t, R, P);
+ EXPECT_EQ(R, Vector_t(1.0-0.05, 0.0, -0.05));
+ EXPECT_EQ(P, Vector_t(0.1, 0.0, 0.0));
+ EXPECT_FALSE(crossing);
+ // betagamma = 0.1 and particle is 0.014 m away; so should succeed after a
+ // couple of RK4 step (which goes about 0.03 m)
+ R = Vector_t(1.0-0.01, 0.0, -0.01);
+ crossing = output_m->checkOne(0, 1.0, 0.0, t, R, P);
+ EXPECT_NEAR(R[0], 1.01, 1e-5);
+ EXPECT_NEAR(R[1], 0.0, 1e-5);
+ EXPECT_NEAR(R[2], -0.01, 1e-5);
+ EXPECT_EQ(P, Vector_t(0.1, 0.0, 0.0));
+ EXPECT_TRUE(crossing);
+ // particle is past the plane but close; should fail after trying a single
+ // RK4 step
+ R = Vector_t(1.0+0.01, 0.0, +0.01);
+ crossing = output_m->checkOne(0, 1.0, 0.0, t, R, P);
+ EXPECT_NEAR(R[0], 1.0+0.01, 1e-5);
+ EXPECT_NEAR(R[1], 0.0, 1e-5);
+ EXPECT_NEAR(R[2], +0.01, 1e-5);
+ EXPECT_EQ(P, Vector_t(0.1, 0.0, 0.0));
+ EXPECT_FALSE(crossing);
+ // particle is past the plane but travelling backwards; should pass
+ R = Vector_t(1.0+0.01, 0.0, +0.01);
+ P = Vector_t(-0.1, 0.0, 0.0);
+ crossing = output_m->checkOne(0, 1.0, 0.0, t, R, P);
+ EXPECT_NEAR(R[0], 1.0-0.01, 1e-5);
+ EXPECT_NEAR(R[1], 0.0, 1e-5);
+ EXPECT_NEAR(R[2], +0.01, 1e-5);
+ EXPECT_EQ(P, Vector_t(-0.1, 0.0, 0.0));
+ EXPECT_TRUE(crossing);
+}
+
+double magnitude(Vector_t v) {
+ return std::sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+}
+
+
+// ps: trying to figure out OPAL units, I put some print statements in ParallelCyclotronTracker
+//
+// itsBunch_m->Q[i]/itsBunch_m->M[i]*Physics::c*Physics::c/Physics::q_e/1e3 = 9.57883e+13
+// std::cerr << "CHARGE TO MASS " << i << " Qi " << itsBunch_m->Q[i] << " Mi "
+// << itsBunch_m->M[i] << " c " << Physics::c << " qe " << Physics::q_e << std::endl;
+// CHARGE TO MASS 0 Qi 1.60218e-19 Mi 0.938272 c 2.99792e+08 qe 1.60218e-19
+
+TEST_F(OutputPlaneTest, TestRK4StepNoField) {
+ Vector_t pVecRef(-1.0, -0.0, -0.4); // moving at beta gamma = 1.414
+ Vector_t rVecRef(2., 3., 4.);
+ Vector_t rVec = rVecRef;
+ Vector_t pVec = pVecRef;
+ double time = 2.0e-9;
+ double tstep = 1.0e-9;
+ double gammabeta = magnitude(pVecRef);
+ double beta = gammabeta/std::sqrt(1+gammabeta*gammabeta);
+ double gamma = gammabeta/beta;
+
+
+ // field not set
+ EXPECT_THROW(output_m->RK4Step(tstep, q2m, time, rVec, pVec), GeneralClassicException);
+
+ // zero field
+ Vector_t speed = pVecRef/gamma*Physics::c; // (beta gamma)* c / gamma = velocity
+ rVecRef = speed*tstep+rVec;
+ std::shared_ptr<MockDipole> testDipole;
+ testDipole.reset(new MockDipole("mock"));
+ testDipole->setField(Vector_t(), Vector_t());
+
+ output_m->setGlobalFieldMap(testDipole.get());
+ output_m->RK4Step(tstep, q2m, time, rVec, pVec);
+ EXPECT_EQ(pVec, pVecRef);
+ for (size_t i = 0; i < 3; ++i) {
+ EXPECT_NEAR(rVec[i], rVecRef[i], 1e-6);
+ }
+}
+
+TEST_F(OutputPlaneTest, TestRK4StepDipole) {
+ // dipole in z direction
+ // B = 2.0 Tesla
+ // p = 1.0 GeV/c
+ // rho = 1e9/2/c
+ double time = 0.0;
+ double tstep = 1.0e-2;
+
+ std::shared_ptr<MockDipole> testDipole;
+ testDipole.reset(new MockDipole("mock"));
+ testDipole->setField(Vector_t(0.0, 0.0, -2.0), Vector_t());
+ Vector_t rVecRef = Vector_t(0.0, 0.0, 0.0);
+ Vector_t pVecRef = Vector_t(0.0, 1.0/Physics::m_p, 1.0);
+ Vector_t rVec = rVecRef;
+ Vector_t pVec = pVecRef;
+ output_m->setGlobalFieldMap(testDipole.get());
+ output_m->RK4Step(tstep, q2m, time, rVec, pVec);
+ EXPECT_NEAR(pVecRef[2], pVec[2], 1e-9);
+ EXPECT_NEAR(pVecRef[0]*pVecRef[0]+pVecRef[1]*pVecRef[1],
+ pVec[0]*pVec[0]+pVec[1]*pVec[1], 1e-8);
+}
+
+TEST_F(OutputPlaneTest, TestRK4StepEField) {
+ // dipole in z direction
+ // B = 2.0 Tesla
+ // p = 1.0 GeV/c
+ // rho = 1e9/2/c
+ double time = 0.0;
+ double tstep = 1.0e-9;
+
+ std::shared_ptr<MockDipole> testDipole;
+ testDipole.reset(new MockDipole("mock"));
+ testDipole->setField(Vector_t(), Vector_t(0.0, 0.0, 2.0));
+ Vector_t rVecRef = Vector_t(1e9/2/Physics::c, 0.0, 0.0);
+ Vector_t pVecRef = Vector_t(0.0, 1.0/0.938272, 0.0);
+ Vector_t rVec = rVecRef;
+ Vector_t pVec = pVecRef;
+ output_m->setGlobalFieldMap(testDipole.get());
+ output_m->RK4Step(tstep, q2m, time, rVec, pVec);
+ EXPECT_NEAR(pVec[0], pVecRef[0], 1e-9);
+ EXPECT_NEAR(pVec[1], pVecRef[1], 1e-9);
+}
+
+TEST_F(OutputPlaneTest, TestCheckOne_FieldOn) {
+ //checkOne(const double tstep, const double chargeToMass,
+ // double& t, Vector_t& R, Vector_t& P)
+ output_m->setGlobalFieldMap(&dipole_m);
+ output_m->setCentre(Vector_t(0.0, 0.0, 0.0));
+ output_m->setNormal(Vector_t(1.0, 0.0, 0.0));
+ output_m->setTolerance(1e-9);
+ output_m->setAlgorithm(OutputPlane::algorithm::RK4STEP);
+
+ Vector_t R, P;
+ bool crossing;
+ double q2m=1.0*Physics::q_e; // proton mass is 0.938
+ // mass excluded since RK4 solvers works in terms of momentum
+ double t = 0.0;
+
+ // betagamma = 0.1 and particle is 0.02 m away; so should succeed after a
+ // couple of RK4 step (which goes about 0.03 m)
+ R = Vector_t(-0.02, 0.0, 0.0);
+ P = Vector_t(0.1, 0.1, 0.0);
+ crossing = output_m->checkOne(0, 1.0, q2m, t, R, P);
+
+ EXPECT_NEAR(R[0], 0.0, 1e-5);
+ EXPECT_NEAR(R[1], 0.02, 1e-5);
+ EXPECT_NEAR(R[2], 0.0, 1e-5);
+ EXPECT_NEAR(P[0], 0.1, 1e-5);
+ EXPECT_NEAR(P[1], 0.1, 1e-5);
+ EXPECT_NEAR(P[2], 0.0, 1e-5);
+
+ EXPECT_TRUE(crossing);
+}
+
+TEST_F(OutputPlaneTest, TestCentre) {
+ // set centre
+ output_m->setCentre(Vector_t(2.0, 3.0, 4.0));
+ output_m->setNormal(Vector_t(1.0, 0.0, 0.0));
+ output_m->setTolerance(1e-9);
+
+ Vector_t R, P, normP, R1, R2, a, b, normaltoplane, unitnormaltoplane;
+ R = Vector_t(-0.02, 5.0, 7.1);
+ P = Vector_t(0.5, 0.8, 0.3);
+
+ //re centre plane
+ output_m->recentre(R, P);
+
+ double Ptot = std::sqrt(P[0]*P[0]+P[1]*P[1]+P[2]*P[2]);
+ normP = Vector_t(P[0]/Ptot, P[1]/Ptot, P[2]/Ptot);
+
+ // check plane is recentred properly
+ EXPECT_EQ(output_m->getCentre(), R);
+ EXPECT_EQ(output_m->getNormal(), normP);
+ }
+
+TEST_F(OutputPlaneTest, TestExtent) {
+ output_m->setCentre(Vector_t(2.0, 3.0, 4.0));
+ output_m->setNormal(Vector_t(1.0, 0.0, 0.0));
+ output_m->setTolerance(1e-9);
+ double halfwidth, halfheight;
+
+ halfheight = 10;
+ halfwidth = 10;
+
+ double horizextent = output_m->getHorizontalExtent();
+ // calling getHorizontalExtent before settingthe horizontal extent returns a value of -1
+ EXPECT_EQ(-1, horizextent);
+
+ output_m->setHorizontalExtent(halfwidth);
+ output_m->setVerticalExtent(halfheight);
+ double vertextent = output_m->getVerticalExtent();
+ horizextent = output_m->getHorizontalExtent();
+
+ EXPECT_EQ( horizextent, halfwidth);
+ EXPECT_EQ( vertextent, halfheight);
+
+}
+
+
+TEST_F(OutputPlaneTest, checkOneTest){
+ OutputPlane* myOutputPlane = new OutputPlane("TestPlane");
+ MockDipole* myDipole = new MockDipole("TestDipole");
+ // check that checkOne returns true when the particle is within one time step
+ // of the OutputPlane and false when it is not within one timestep
+ myOutputPlane->setGlobalFieldMap(myDipole);
+
+ Vector_t R, P;
+ Vector_t R1, P1;
+
+ double tStep = 1.0e-9; // s
+
+ Vector_t normal = {1.0, 0.0, 1.0};
+ Vector_t centre = {0.0, 0.0, 1.0};
+ myOutputPlane->setCentre(centre);
+ myOutputPlane->setNormal(normal);
+ myOutputPlane->setTolerance(1e-9);
+
+ R = Vector_t(0.1, 0.0, 0.1);
+ P = Vector_t(0.1, 0.0, 0.0);
+ double t = 0.0;
+ bool cross;
+ cross = myOutputPlane->checkOne(0, tStep, 0.0, t, R, P);
+ EXPECT_FALSE(cross);
+
+ // opposite case where it is near to the plane
+ // distance from plane
+ R1 = Vector_t(1.0-0.01, 0.0, -0.01);
+ P1 = Vector_t(0.1, 0.0, 0.0);
+
+ bool cross1;
+ cross1 = myOutputPlane->checkOne(0, tStep, 0.0, t, R1, P1);
+ EXPECT_TRUE(cross1);
+ delete myDipole;
+ delete myOutputPlane;
+}
diff --git a/tests/opal_src/Elements/CMakeLists.txt b/tests/opal_src/Elements/CMakeLists.txt
index 127c036de671623fead20c9a9749d53ab976dbde..18c7fe3b06231a2346ecafcea0d56e9965772cc1 100644
--- a/tests/opal_src/Elements/CMakeLists.txt
+++ b/tests/opal_src/Elements/CMakeLists.txt
@@ -2,6 +2,7 @@ set (_SRCS
OpalAsymmetricEngeTest.cpp
OpalEngeTest.cpp
OpalOffsetTest.cpp
+ OpalOutputPlaneTest.cpp
OpalPolynomialTimeDependenceTest.cpp
OpalScalingFFAMagnetTest.cpp
OpalSplineTimeDependenceTest.cpp
diff --git a/tests/opal_src/Elements/OpalOutputPlaneTest.cpp b/tests/opal_src/Elements/OpalOutputPlaneTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..400e5f201212feac377f3d8697902ace9beef475
--- /dev/null
+++ b/tests/opal_src/Elements/OpalOutputPlaneTest.cpp
@@ -0,0 +1,85 @@
+// Copyright (c) 2023, Chris Rogers
+// 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 "gtest/gtest.h"
+
+#include "Attributes/Attributes.h"
+#include "AbsBeamline/Offset.h"
+#include "Elements/OpalOutputPlane.h"
+#include "AbsBeamline/OutputPlane.h"
+
+TEST(OpalOutputPlaneTest, TestCtor) {
+ OpalOutputPlane plane1;
+ EXPECT_EQ((&plane1)->getOpalName(), "OUTPUTPLANE");
+}
+
+
+void setReal(OpalOutputPlane& plane, std::string attName, double value) {
+ Attribute* att = plane.findAttribute(attName);
+ ASSERT_NE(att, nullptr) << "Failed with attribute " << attName;
+ Attributes::setReal(*att, value);
+}
+
+TEST(OpalOutputPlaneTest, TestProbeMode) {
+ OpalOutputPlane plane;
+ Attribute* att = plane.findAttribute("PLACEMENT_STYLE");
+ Attributes::setPredefinedString(*att, "PROBE");
+ setReal(plane, "XSTART", 1.0);
+ setReal(plane, "XEND", 3.0);
+ setReal(plane, "YSTART", -1.0);
+ setReal(plane, "YEND", -3.0);
+ plane.update();
+ OutputPlane *output = dynamic_cast<OutputPlane *>(plane.getElement());
+ double roottwo = std::sqrt(2);
+ EXPECT_NEAR(output->getCentre()[0], 2.0, 1e-12);
+ EXPECT_NEAR(output->getCentre()[1], -2.0, 1e-12);
+ EXPECT_NEAR(output->getCentre()[2], 0.0, 1e-12);
+ EXPECT_NEAR(output->getHorizontalExtent(), roottwo, 1e-12); // that is half width
+ EXPECT_NEAR(output->getNormal()[0], 1/roottwo, 1e-12);
+ EXPECT_NEAR(output->getNormal()[1], 1/roottwo, 1e-12);
+
+ // swap and check normal is reversed
+ setReal(plane, "XSTART", 3.0);
+ setReal(plane, "XEND", 1.0);
+ setReal(plane, "YSTART", -3.0);
+ setReal(plane, "YEND", -1.0);
+ plane.update();
+ EXPECT_NEAR(output->getNormal()[0], -1/roottwo, 1e-12);
+ EXPECT_NEAR(output->getNormal()[1], -1/roottwo, 1e-12);
+}
+
+TEST(OpalOutputPlaneTest, TestPlacement) {
+ OpalOutputPlane plane;
+ Attribute* att0 = plane.findAttribute("PLACEMENT_STYLE");
+ Attributes::setPredefinedString(*att0, "CENTRE_NORMAL");
+ Attribute* att1 = plane.findAttribute("CENTRE");
+ Attributes::setRealArray(*att1, {1.0, 2.0, 3.0});
+ Attribute* att2 = plane.findAttribute("NORMAL");
+ Attributes::setRealArray(*att2, {-4.0, 4.0, 2.0});
+ setReal(plane, "WIDTH", 4.0);
+ setReal(plane, "HEIGHT", 5.0);
+ setReal(plane, "RADIUS", 6.0);
+ plane.update();
+ OutputPlane *output = dynamic_cast<OutputPlane *>(plane.getElement());
+ EXPECT_NEAR(output->getCentre()[0], 1.0, 1e-12);
+ EXPECT_NEAR(output->getCentre()[1], 2.0, 1e-12);
+ EXPECT_NEAR(output->getCentre()[2], 3.0, 1e-12);
+ EXPECT_NEAR(output->getNormal()[0], -2/3.0, 1e-12);
+ EXPECT_NEAR(output->getNormal()[1], 2/3.0, 1e-12);
+ EXPECT_NEAR(output->getNormal()[2], 1.0/3.0, 1e-12);
+ EXPECT_NEAR(output->getHorizontalExtent(), 2.0, 1e-12); // that is half width
+ EXPECT_NEAR(output->getVerticalExtent(), 2.5, 1e-12); // that is half height
+ EXPECT_NEAR(output->getRadialExtent(), 6.0, 1e-12); // that is half height
+}
diff --git a/tests/opal_src/PyOpal/PyElements/test_output_plane.py b/tests/opal_src/PyOpal/PyElements/test_output_plane.py
new file mode 100644
index 0000000000000000000000000000000000000000..f4beddea88c06eb0b59a8f4cbc0c3909e9970954
--- /dev/null
+++ b/tests/opal_src/PyOpal/PyElements/test_output_plane.py
@@ -0,0 +1,265 @@
+# Copyright (c) 2023, Chris Rogers, STFC Rutherford Appleton Laboratory, Didcot, UK
+#
+# 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/>.
+
+import os
+import unittest
+import math
+import pyopal.objects.minimal_runner
+
+import pyopal.objects.field
+import pyopal.objects.option
+import pyopal.elements.output_plane
+import pyopal.elements.multipolet
+import pyopal.elements.global_cartesian_offset
+import pyopal.elements.polynomial_time_dependence
+import pyopal.elements.variable_rf_cavity
+
+
+class OutputPlaneRunner(pyopal.objects.minimal_runner.MinimalRunner):
+ """Class to set up a probe and run it"""
+ def __init__(self, number_of_steps=100, beta_rel=0.001, bz=0.0, algorithm="RK4", verbose=0):
+ """Set up the test"""
+ super().__init__()
+ self.verbose = verbose
+ self.probe = None
+ self.options = None
+ self.r0 = 0.9
+ self.yprobe = 0.1 # that should be in [m]
+ self.algorithm = algorithm
+
+ self.bz = bz
+ # Tracking is set up so the beam will go 0.11 metres
+ self.steps_per_turn = number_of_steps-1
+ self.max_steps = number_of_steps
+ self.time_per_turn = 0.11/(beta_rel*3e8) # [seconds] beam goes 0.11 metres
+ self.momentum = beta_rel/(1-beta_rel**2)*self.mass # [GeV/c]
+ self.distribution_str = """1
+0.0 0.0 0.0 0.0 0.0 0.0
+"""
+ self.setup_plane()
+
+ def make_options(self):
+ """Don't want to force user to have hdf2py - so dump to ASCII"""
+ self.options = pyopal.objects.option.Option()
+ self.options.ascii_dump = True
+ self.options.enable_hdf5 = False
+ self.options.execute()
+
+ def setup_plane(self):
+ """Make a probe, just offset from x axis so should catch first step"""
+ self.output_plane = pyopal.elements.output_plane.OutputPlane()
+ self.output_plane.set_attributes(
+ x_start = 0.1,
+ x_end = self.r0*2,
+ y_start = self.yprobe,
+ y_end = self.yprobe,
+ placement_style = "PROBE",
+ output_filename = "test_plane",
+ algorithm = self.algorithm,
+ verbose_level = self.verbose
+ )
+ self.output_plane.set_opal_name("my_plane")
+
+ def make_dipole(self):
+ dipole = pyopal.elements.multipolet.MultipoleT()
+ dipole.set_attributes(
+ t_p = [self.bz, self.bz*100],
+ left_fringe = 0.01,
+ right_fringe = 0.01,
+ length = 0.5,
+ horizontal_aperture = 100,
+ vertical_aperture = 100,
+ bounding_box_length = 1.0
+ )
+ return [dipole]
+
+ @classmethod
+ def make_time_dependence(cls, name, pol0, pol1):
+ """Make a time dependence"""
+ ptd = pyopal.elements.polynomial_time_dependence.PolynomialTimeDependence()
+ ptd.p0 = pol0
+ ptd.p1 = pol1
+ ptd.set_opal_name(name)
+ ptd.update()
+ return ptd
+
+ def make_rf(self):
+ rf = pyopal.elements.variable_rf_cavity.VariableRFCavity()
+ phase = self.make_time_dependence("phase", math.pi/2, 0.0)
+ # nb: 4 GeV protons; significant energy change in 0.1 m -> O(GV/m)
+ voltage = self.make_time_dependence("voltage", self.bz*1000.0, 0.0)
+ # time step is ~ ns so frequency ~ GHz is appropriate
+ frequency = self.make_time_dependence("frequency", 1e3, 0.0)
+ rf.set_attributes(
+ phase_model = "phase",
+ amplitude_model = "voltage",
+ frequency_model = "frequency",
+ height = 10.0,
+ width = 10.0,
+ length = 100.0,
+ )
+ offset = pyopal.elements.global_cartesian_offset.GlobalCartesianOffset()
+ offset.set_attributes(
+ end_position_x = 0.9,
+ end_position_y = -0.5,
+ end_normal_x = 0.0,
+ end_normal_y = 1.0,
+ )
+ return [offset, rf]
+
+ def make_element_iterable(self):
+ self.make_options()
+ return [self.output_plane]+self.make_dipole()+self.make_rf()
+
+class TestOutputPlane(unittest.TestCase):
+
+ def run_one(self, n_steps, beta_rel, bz, algorithm):
+ plane_runner = OutputPlaneRunner(n_steps, beta_rel, bz, algorithm, self.verbose)
+ plane_runner.execute_fork() # forked process for memory safety
+ return self.read_lines(plane_runner)
+
+ def read_lines(self, plane_runner, name = "test_plane"):
+ plane_file = os.path.join(plane_runner.tmp_dir, name+".loss")
+ try:
+ with open(plane_file) as a_file:
+ lines = [line for line in a_file.readlines()]
+ except FileNotFoundError:
+ print(f"File {plane_file} not found - ls returns:")
+ print(" ", os.listdir(plane_runner.tmp_dir))
+ raise
+ return lines
+
+ def test_no_field_rk4(self):
+ """Simple test - track through a probe and check we get an output"""
+ lines = self.run_one(1000, 0.001, 0.0, "INTERPOLATION")
+ self.assertEqual(len(lines), 2) # header and one track object
+ words = lines[1].split()
+ self.assertAlmostEqual(float(words[0]), 0.9)
+ self.assertAlmostEqual(float(words[1]), 0.1) # output is in [m]
+ self.assertAlmostEqual(float(words[2]), 0.0)
+ self.assertAlmostEqual(float(words[3]), 0.0)
+ self.assertAlmostEqual(float(words[4]), 0.001) # beta_gamma
+ self.assertAlmostEqual(float(words[5]), 0.0)
+
+ def test_field_rk4(self):
+ """
+ Track through a probe and check we get an output; check that RK4 is
+ consistent with linear interpolation in the limit of many interpolation
+ steps. Note fields are set up just to be "reasonably" complicated i.e. I
+ want to check fields in several different directions, time varying, etc
+ all works okay.
+ """
+ int_lines = self.run_one(1000, 0.9, 2.0, "INTERPOLATION")
+ rk4_lines = self.run_one(3, 0.9, 2.0, "RK4")
+ if self.verbose:
+ print(rk4_lines)
+ print("Interpolation\n ", int_lines[1])
+ print("RK4\n ", rk4_lines[1])
+ self.assertEqual(len(int_lines), 2) # header and one track object
+ self.assertEqual(len(rk4_lines), 2) # header and one track object
+ rk4_words = rk4_lines[1].split()
+ int_words = int_lines[1].split()
+ self.assertAlmostEqual(float(rk4_words[3]), float(int_words[3]), 4)
+ self.assertAlmostEqual(float(rk4_words[4]), float(int_words[4]), 4) # beta_gamma
+ self.assertAlmostEqual(float(rk4_words[5]), float(int_words[5]), 4)
+ self.assertAlmostEqual(float(rk4_words[0]), float(int_words[0]), 9)
+ self.assertAlmostEqual(float(rk4_words[1]), float(int_words[1]), 9) # output is in [m]
+ self.assertAlmostEqual(float(rk4_words[2]), float(int_words[2]), 9)
+ self.assertAlmostEqual(float(rk4_words[9])*1e9, float(int_words[9])*1e9, 3) # ns
+
+ def test_reference_alignment_particle(self):
+ """Check plane realignment works"""
+ plane_runner = OutputPlaneRunner(100, 0.001, 0.0, "INTERPOLATION", self.verbose)
+ plane_runner.output_plane.reference_alignment_particle = 1
+ plane_runner.output_plane.tolerance = 1e-10
+ # momentum units are eV(!)
+ dpx = 0.0001*0.938272e9 # dxds = 0.1
+ plane_runner.distribution_str = """3
+0.0 0.0 0.0 0.0 0.0 0.0
+0.0 """+str(dpx)+""" 0.0 0.0 0.0 0.0
+-0.01 0.0 0.0 0.0 0.0 0.0
+"""
+ plane_runner.execute_fork() # forked process for memory safety
+ lines = self.read_lines(plane_runner)
+ if self.verbose:
+ print("Reference plane test")
+ for line in lines:
+ print(line[:-1])
+ self.assertEqual(len(lines), 4)
+ dxds = 0.1 # the alignment particle angle
+ dx = 0.02 # alignment particle offset
+ test_particle = lines[3].split()
+ self.assertAlmostEqual(float(test_particle[0]), 0.89)
+ self.assertAlmostEqual(float(test_particle[1])-0.1, dx*dxds)
+
+ def test_rectangular_aperture(self):
+ """Check we cut approporiately with rectangular aperture"""
+ width = 0.01
+ plane_runner = OutputPlaneRunner(100, 0.001, 0.0, "INTERPOLATION", self.verbose)
+ plane_runner.output_plane.placement_style = "CENTRE-NORMAL"
+ plane_runner.output_plane.centre = [0.9, 0.1, 0.0]
+ plane_runner.output_plane.normal = [1.0, 1.0, 0.0]
+ plane_runner.output_plane.height = 0.1
+ plane_runner.output_plane.width = width*2
+ plane_runner.distribution_str = f"""9
+0.0 0.0 0.0 0.0 0.0 0.0
+{(width-0.001)*0.5**0.5} 0.0 0.0 0.0 0.0 0.0
+{(width+0.001)*0.5**0.5} 0.0 0.0 0.0 0.0 0.0
+{(-width-0.001)*0.5**0.5} 0.0 0.0 0.0 0.0 0.0
+{(-width+0.001)*0.5**0.5} 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.049 0.0
+0.0 0.0 0.0 0.0 0.051 0.0
+0.0 0.0 0.0 0.0 -0.049 0.0
+0.0 0.0 0.0 0.0 -0.051 0.0
+"""
+
+ plane_runner.execute_fork() # forked process for memory safety
+ lines = self.read_lines(plane_runner)
+ if self.verbose:
+ print("Reference plane test")
+ for line in lines:
+ print(line[:-1])
+ self.assertEqual(len(lines), 6)
+
+ def test_radial_aperture(self):
+ """Check we cut appropriately with radial aperture"""
+ plane_runner = OutputPlaneRunner(100, 0.001, 0.0, "INTERPOLATION", self.verbose)
+ plane_runner.output_plane.placement_style = "CENTRE-NORMAL"
+ plane_runner.output_plane.output_filename = "" # defaults to "my_plane"
+ plane_runner.output_plane.centre = [0.9, 0.1, 0.0]
+ plane_runner.output_plane.normal = [1.0, 1.0, 0.0]
+ plane_runner.output_plane.radius = 0.1
+ width = 0.1
+ plane_runner.distribution_str = f"""9
+0.0 0.0 0.0 0.0 0.0 0.0
+{(width-0.001)*0.5**0.5} 0.0 0.0 0.0 0.0 0.0
+{(width+0.001)*0.5**0.5} 0.0 0.0 0.0 0.0 0.0
+{(-width-0.001)*0.5**0.5} 0.0 0.0 0.0 0.0 0.0
+{(-width+0.001)*0.5**0.5} 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 {width-0.001} 0.0
+0.0 0.0 0.0 0.0 {width+0.001} 0.0
+0.0 0.0 0.0 0.0 {-width-0.001} 0.0
+0.0 0.0 0.0 0.0 {-width+0.001} 0.0
+"""
+
+ plane_runner.execute_fork() # forked process for memory safety
+ lines = self.read_lines(plane_runner, "my_plane")
+ if self.verbose:
+ print("Reference plane test")
+ for line in lines:
+ print(line[:-1])
+ self.assertEqual(len(lines), 5)
+
+ verbose = 0
+
+if __name__ == "__main__":
+ unittest.main()
\ No newline at end of file