Skip to content

GitLab

S
src
Commit 73edf13f authored by snuverink_j's avatar snuverink_j
Browse files
Options

Merge branch '307-pluginelements-use-unit-meter-internally-to-avoid-conversions' into 'master' 

Resolve "PluginElements: use unit meter internally to avoid conversions"

Closes #539 and #307

See merge request !368
parents 59d78b5b 0d0bf6e2
Hide whitespace changes
Inline Side-by-side
Showing with 267 additions and 218 deletions
src/Algorithms/ParallelCyclotronTracker.cpp
View file @ 73edf13f
......@@ -621,35 +621,35 @@ void ParallelCyclotronTracker::visitCCollimator(const CCollimator &coll) {
myElements.push_back(elptr);
double xstart = elptr->getXStart();
*gmsg << "* Xstart = " << xstart << " [mm]" << endl;
*gmsg << "* Xstart = " << xstart << " [m]" << endl;
double xend = elptr->getXEnd();
*gmsg << "* Xend = " << xend << " [mm]" << endl;
*gmsg << "* Xend = " << xend << " [m]" << endl;
double ystart = elptr->getYStart();
*gmsg << "* Ystart = " << ystart << " [mm]" << endl;
*gmsg << "* Ystart = " << ystart << " [m]" << endl;
double yend = elptr->getYEnd();
*gmsg << "* Yend = " << yend << " [mm]" << endl;
*gmsg << "* Yend = " << yend << " [m]" << endl;
double zstart = elptr->getZStart();
*gmsg << "* Zstart = " << zstart << " [mm]" << endl;
*gmsg << "* Zstart = " << zstart << " [m]" << endl;
double zend = elptr->getZEnd();
*gmsg << "* Zend = " << zend << " [mm]" << endl;
*gmsg << "* Zend = " << zend << " [m]" << endl;
double width = elptr->getWidth();
*gmsg << "* Width = " << width << " [mm]" << endl;
*gmsg << "* Width = " << width << " [m]" << endl;
elptr->initialise(itsBunch_m);
double BcParameter[8] = {};
BcParameter[0] = 0.001 * xstart ;
BcParameter[1] = 0.001 * xend;
BcParameter[2] = 0.001 * ystart ;
BcParameter[3] = 0.001 * yend;
BcParameter[4] = 0.001 * width ;
BcParameter[0] = xstart;
BcParameter[1] = xend;
BcParameter[2] = ystart;
BcParameter[3] = yend;
BcParameter[4] = width;
buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
}
......@@ -835,27 +835,27 @@ void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
*gmsg << "* Name = " << elptr->getName() << endl;
double xstart = elptr->getXStart();
*gmsg << "* XStart = " << xstart << " [mm]" << endl;
*gmsg << "* XStart = " << xstart << " [m]" << endl;
double xend = elptr->getXEnd();
*gmsg << "* XEnd = " << xend << " [mm]" << endl;
*gmsg << "* XEnd = " << xend << " [m]" << endl;
double ystart = elptr->getYStart();
*gmsg << "* YStart = " << ystart << " [mm]" << endl;
*gmsg << "* YStart = " << ystart << " [m]" << endl;
double yend = elptr->getYEnd();
*gmsg << "* YEnd = " << yend << " [mm]" << endl;
*gmsg << "* YEnd = " << yend << " [m]" << endl;
// initialise, do nothing
elptr->initialise(itsBunch_m);
double BcParameter[8] = {};
BcParameter[0] = 0.001 * xstart ;
BcParameter[1] = 0.001 * xend;
BcParameter[2] = 0.001 * ystart ;
BcParameter[3] = 0.001 * yend;
BcParameter[4] = 0.001 ; // width
BcParameter[0] = xstart;
BcParameter[1] = xend;
BcParameter[2] = ystart;
BcParameter[3] = yend;
BcParameter[4] = 1 ; // width
// store probe parameters in the list
buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
......@@ -1052,19 +1052,19 @@ void ParallelCyclotronTracker::visitSeptum(const Septum &sept) {
myElements.push_back(elptr);
double xstart = elptr->getXStart();
*gmsg << "* XStart = " << xstart << " [mm]" << endl;
*gmsg << "* XStart = " << xstart << " [m]" << endl;
double xend = elptr->getXEnd();
*gmsg << "* XEnd = " << xend << " [mm]" << endl;
*gmsg << "* XEnd = " << xend << " [m]" << endl;
double ystart = elptr->getYStart();
*gmsg << "* YStart = " << ystart << " [mm]" << endl;
*gmsg << "* YStart = " << ystart << " [m]" << endl;
double yend = elptr->getYEnd();
*gmsg << "* YEnd = " << yend << " [mm]" << endl;
*gmsg << "* YEnd = " << yend << " [m]" << endl;
double width = elptr->getWidth();
*gmsg << "* Width = " << width << " [mm]" << endl;
*gmsg << "* Width = " << width << " [m]" << endl;
// initialise, do nothing
......@@ -1072,11 +1072,11 @@ void ParallelCyclotronTracker::visitSeptum(const Septum &sept) {
double BcParameter[8] = {};
BcParameter[0] = 0.001 * xstart ;
BcParameter[1] = 0.001 * xend;
BcParameter[2] = 0.001 * ystart ;
BcParameter[3] = 0.001 * yend;
BcParameter[4] = 0.001 * width ;
BcParameter[0] = xstart;
BcParameter[1] = xend;
BcParameter[2] = ystart;
BcParameter[3] = yend;
BcParameter[4] = width;
// store septum parameters in the list
buildupFieldList(BcParameter, ElementBase::SEPTUM, elptr);
......@@ -1119,18 +1119,18 @@ void ParallelCyclotronTracker::visitStripper(const Stripper &stripper) {
myElements.push_back(elptr);
*gmsg << "* Name = " << elptr->getName() << endl;
double xstart = elptr->getXStart();
*gmsg << "* XStart = " << xstart << " [mm]" << endl;
*gmsg << "* XStart = " << xstart << " [m]" << endl;
double xend = elptr->getXEnd();
*gmsg << "* XEnd = " << xend << " [mm]" << endl;
*gmsg << "* XEnd = " << xend << " [m]" << endl;
double ystart = elptr->getYStart();
*gmsg << "* YStart = " << ystart << " [mm]" << endl;
*gmsg << "* YStart = " << ystart << " [m]" << endl;
double yend = elptr->getYEnd();
*gmsg << "* YEnd = " << yend << " [mm]" << endl;
*gmsg << "* YEnd = " << yend << " [m]" << endl;
double opcharge = elptr->getOPCharge();
*gmsg << "* Charge of outcoming particle = +e * " << opcharge << endl;
......@@ -1140,16 +1140,16 @@ void ParallelCyclotronTracker::visitStripper(const Stripper &stripper) {
bool stop = elptr->getStop();
*gmsg << std::boolalpha << "* Particles stripped will be deleted after interaction -> " << stop << endl;
elptr->initialise(itsBunch_m);
double BcParameter[8] = {};
BcParameter[0] = 0.001 * xstart ;
BcParameter[1] = 0.001 * xend;
BcParameter[2] = 0.001 * ystart ;
BcParameter[3] = 0.001 * yend;
BcParameter[4] = 0.001; // width
BcParameter[0] = xstart;
BcParameter[1] = xend;
BcParameter[2] = ystart;
BcParameter[3] = yend;
BcParameter[4] = 1; // width
BcParameter[5] = opcharge;
BcParameter[6] = opmass;
......@@ -2120,8 +2120,6 @@ void ParallelCyclotronTracker::borisExternalFields(double h) {
bool ParallelCyclotronTracker::applyPluginElements(const double dt) {
IpplTimings::startTimer(PluginElemTimer_m);
// Plugin Elements are all defined in mm, change beam to mm before applying
itsBunch_m->R *= Vector_t(1000.0);
for(beamline_list::iterator sindex = ++(FieldDimensions.begin()); sindex != FieldDimensions.end(); ++sindex) {
if(((*sindex)->first) == ElementBase::BEAMSTRIPPING) {
......@@ -2146,7 +2144,6 @@ bool ParallelCyclotronTracker::applyPluginElements(const double dt) {
}
}
itsBunch_m->R *= Vector_t(0.001);
IpplTimings::stopTimer(PluginElemTimer_m);
return flag;
}
......@@ -3061,7 +3058,7 @@ void ParallelCyclotronTracker::seoMode_m(double& t, const double dt, bool& /*fin
// 2 particles: Trigger SEO mode
// (Switch off cavity and calculate betatron oscillation tuning)
double r_tuning[2], z_tuning[2] ;
double r_tuning[2], z_tuning[2];
IpplTimings::startTimer(IntegrationTimer_m);
for(size_t i = 0; i < (itsBunch_m->getLocalNum()); i++) {
......
src/Classic/AbsBeamline/BeamStripping.cpp
View file @ 73edf13f
......@@ -177,10 +177,10 @@ bool BeamStripping::checkBeamStripping(PartBunchBase<double, 3> *bunch, Cyclotro
boundingSphere.first = 0.5 * (rmax + rmin);
boundingSphere.second = euclidean_norm(rmax - boundingSphere.first);
maxr_m = 1000 * cycl->getMaxR();
minr_m = 1000 * cycl->getMinR();
maxz_m = 1000 * cycl->getMaxZ();
minz_m = 1000 * cycl->getMinZ();
maxr_m = cycl->getMaxR();
minr_m = cycl->getMinR();
maxz_m = cycl->getMaxZ();
minz_m = cycl->getMinZ();
size_t tempnum = bunch->getLocalNum();
for (unsigned int i = 0; i < tempnum; ++i) {
......@@ -256,8 +256,7 @@ std::string BeamStripping::getBeamStrippingShape() {
int BeamStripping::checkPoint(const double &x, const double &y, const double &z) {
int cn;
double rpos = std::sqrt(x * x + y * y);
double zpos = z;
if (zpos >= maxz_m || zpos <= minz_m || rpos >= maxr_m || rpos <= minr_m)
if (z >= maxz_m || z <= minz_m || rpos >= maxr_m || rpos <= minr_m)
cn = 0;
else
cn = 1;
......
src/Classic/AbsBeamline/CCollimator.cpp
View file @ 73edf13f
//
// Class CCollimator
// Interface for cyclotron collimator
//
// 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/CCollimator.h"
#include "AbsBeamline/BeamlineVisitor.h"
......@@ -99,14 +116,6 @@ void CCollimator::doInitialise(PartBunchBase<double, 3> */*bunch*/) {
void CCollimator::goOnline(const double &) {
print();
// PosX_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
// PosY_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
// PosZ_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
// MomentumX_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
// MomentumY_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
// MomentumZ_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
// time_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
// id_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
online_m = true;
}
......@@ -128,9 +137,6 @@ void CCollimator::print() {
}
return;
}
*gmsg << "* CCollimator angle start " << xstart_m << " (Deg) angle end " << ystart_m << " (Deg) "
<< "R start " << xend_m << " (mm) R rend " << yend_m << " (mm)" << endl;
}
void CCollimator::setDimensions(double xstart, double xend, double ystart, double yend, double zstart, double zend, double width) {
......
src/Classic/AbsBeamline/CCollimator.h
View file @ 73edf13f
//
// Class CCollimator
// Interface for cyclotron collimator
//
// 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/>.
//
#ifndef CLASSIC_CCollimator_HH
#define CLASSIC_CCollimator_HH
......@@ -5,11 +22,6 @@
class ParticleMatterInteractionHandler;
// Class CCollimator
// ------------------------------------------------------------------------
/// Interface for cyclotron collimator.
// Class CCollimator defines the abstract interface for a collimator.
class CCollimator: public PluginElement {
public:
......
src/Classic/AbsBeamline/PluginElement.cpp
View file @ 73edf13f
//
// 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"
......@@ -11,15 +29,13 @@ PluginElement::PluginElement():PluginElement("")
PluginElement::PluginElement(const std::string &name):
Component(name),
filename_m(""),
position_m(0.0) {
filename_m("") {
setDimensions(0.0, 0.0, 0.0, 0.0);
}
PluginElement::PluginElement(const PluginElement &right):
Component(right),
filename_m(right.filename_m),
position_m(right.position_m) {
filename_m(right.filename_m) {
setDimensions(right.xstart_m, right.xend_m, right.ystart_m, right.yend_m);
}
......@@ -138,9 +154,9 @@ void PluginElement::setGeom(const double dist) {
void PluginElement::changeWidth(PartBunchBase<double, 3> *bunch, int i, const double tstep, const double tangle) {
constexpr double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns
double lstep = euclidean_norm(bunch->P[i]) / Util::getGamma(bunch->P[i]) * c_mmtns * tstep; // [mm]
double sWidth = lstep / std::sqrt( 1 + 1/tangle/tangle );
constexpr double c_mtns = Physics::c * 1.0e-9; // m/s --> m/ns
double lstep = euclidean_norm(bunch->P[i]) / Util::getGamma(bunch->P[i]) * c_mtns * tstep; // [m]
double sWidth = lstep / std::sqrt( 1 + 1/tangle/tangle );
setGeom(sWidth);
}
......@@ -156,7 +172,7 @@ double PluginElement::calculateIncidentAngle(double xp, double yp) const {
else
tangle = std::abs(1 / k1);
} else if ( xp == 0.0 ) {
k2 = - A_m/B_m;
k2 = - A_m / B_m;
if ( k2 == 0.0 )
tangle = 1.0e12;
else
......@@ -200,15 +216,15 @@ bool PluginElement::check(PartBunchBase<double, 3> *bunch, const int turnnumber,
}
void PluginElement::getDimensions(double &zBegin, double &zEnd) const {
zBegin = position_m - 0.005;
zEnd = position_m + 0.005;
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))) {
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))
......
src/Classic/AbsBeamline/PluginElement.h
View file @ 73edf13f
//
// 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/>.
//
#ifndef CLASSIC_PluginElement_HH
#define CLASSIC_PluginElement_HH
......@@ -9,12 +27,6 @@ template <class T, unsigned Dim>
class PartBunchBase;
class LossDataSink;
// Class PluginElement
// ------------------------------------------------------------------------
/// Abstract Interface for (Cyclotron) Plugin Elements (Probe, Stripper, Septum etc.)
/// Implementation via Non-Virtual Interface Template Method
class PluginElement: public Component {
public:
......@@ -99,7 +111,6 @@ private:
protected:
/* Members */
std::string filename_m; /**< The name of the outputfile*/
double position_m;
///@{ input geometry positions
double xstart_m;
double xend_m;
......
src/Classic/AbsBeamline/Probe.cpp
View file @ 73edf13f
//
// Class Probe
// Interface for a probe
//
// Copyright (c) 2016-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/Probe.h"
#include "AbsBeamline/BeamlineVisitor.h"
......@@ -59,7 +76,7 @@ bool Probe::doPreCheck(PartBunchBase<double, 3> *bunch) {
double rbunch_min = std::hypot(xmin, ymin);
double rbunch_max = std::hypot(xmax, ymax);
if( rbunch_max > rmin_m - 10.0 && rbunch_min < rend_m + 10.0 ) {
if (rbunch_max > rmin_m - 0.01 && rbunch_min < rend_m + 0.01 ) {
return true;
}
return false;
......@@ -69,11 +86,11 @@ bool Probe::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const
Vector_t probepoint;
size_t tempnum = bunch->getLocalNum();
for(unsigned int i = 0; i < tempnum; ++i) {
for (unsigned int i = 0; i < tempnum; ++i) {
double tangle = calculateIncidentAngle(bunch->P[i](0), bunch->P[i](1));
changeWidth(bunch, i, tstep, tangle);
int pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
if(pflag == 0) continue;
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)
......@@ -82,19 +99,14 @@ bool Probe::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const
// 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 dist1 = (A_m * bunch->R[i](0) + B_m * bunch->R[i](1) + C_m) / R_m; // [m]
double dist2 = dist1 * std::sqrt(1.0 + 1.0 / tangle / tangle);
double dt = dist2 / (std::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]);
probepoint = bunch->R[i] + dist2 * 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;
peakfinder_m->addParticle(probepoint * 1e3);
lossDs_m->addParticle(probepoint, bunch->P[i], bunch->ID[i], t+dt,
turnnumber, bunch->bunchNum[i]);
......
src/Classic/AbsBeamline/Probe.h
View file @ 73edf13f
//
// Class Probe
// Interface for a probe
//
// Copyright (c) 2016-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/>.
//
#ifndef CLASSIC_Probe_HH
#define CLASSIC_Probe_HH
......@@ -8,11 +25,6 @@
class PeakFinder;
// Class Probe
// ------------------------------------------------------------------------
/// Interface for probe.
// Class Probe defines the abstract interface for a probe.
class Probe: public PluginElement {
public:
......
src/Classic/AbsBeamline/Septum.cpp
View file @ 73edf13f
//
// Class Septum
// Interface for a septum magnet
//
// Copyright (c) 2016-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/Septum.h"
#include "AbsBeamline/BeamlineVisitor.h"
......@@ -7,9 +24,6 @@
extern Inform *gmsg;
// Class Septum
// ------------------------------------------------------------------------
Septum::Septum():Septum("")
{}
......@@ -31,9 +45,8 @@ void Septum::accept(BeamlineVisitor &visitor) const {
}
void Septum::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) {
position_m = startField;
endField = startField + 0.005;
startField -= 0.005;
endField = position_m + 0.005;
initialise(bunch);
}
......@@ -59,7 +72,7 @@ bool Septum::doPreCheck(PartBunchBase<double, 3> *bunch) {
double ymax = std::max(std::abs(rmin(1)), std::abs(rmax(1)));
double rbunch_max = std::hypot(xmax, ymax);
if(rbunch_max > rstart_m - 100) {
if (rbunch_max > rstart_m - 0.1) {
return true;
}
return false;
......@@ -74,20 +87,20 @@ bool Septum::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/,
const double intcept1 = intcept - halfLength;
const double intcept2 = intcept + halfLength;
for(unsigned int i = 0; i < bunch->getLocalNum(); ++i) {
for (unsigned int i = 0; i < bunch->getLocalNum(); ++i) {
const Vector_t& R = bunch->R[i];
double line1 = std::abs(slope * R(0) + intcept1);
double line2 = std::abs(slope * R(0) + intcept2);
if(std::abs(R(1)) > line2 &&
if (std::abs(R(1)) > line2 &&
std::abs(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]);
lossDs_m->addParticle(R, bunch->P[i], bunch->ID[i]);
bunch->Bin[i] = -1;
flag = true;
}
......
src/Classic/AbsBeamline/Septum.h
View file @ 73edf13f
//
// Class Septum
// Interface for a septum magnet
//
// Copyright (c) 2016-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/>.
//
#ifndef CLASSIC_Septum_HH
#define CLASSIC_Septum_HH
#include "AbsBeamline/PluginElement.h"
// Class Septum
// ------------------------------------------------------------------------
/// Interface for septum magnet.
// Class Septum defines the abstract interface for a septum magnet.
class Septum: public PluginElement {
public:
......
src/Classic/AbsBeamline/Stripper.cpp
View file @ 73edf13f
......@@ -3,7 +3,7 @@
// The Stripper element defines the interface for a stripping foil
//
// Copyright (c) 2011, Jianjun Yang, Paul Scherrer Institut, Villigen PSI, Switzerland
// Copyright (c) 2017-2019, Paul Scherrer Institut, Villigen PSI, Switzerland
// Copyright (c) 2017-2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// Implemented as part of the PhD thesis
......@@ -101,7 +101,7 @@ bool Stripper::doPreCheck(PartBunchBase<double, 3> *bunch) {
double ymax = std::max(std::abs(rmin(1)), std::abs(rmax(1)));
double rbunch_max = std::hypot(xmax, ymax);
if (rbunch_max > rmin_m - 10.0) {
if (rbunch_max > rmin_m - 1e-2) {
return true;
}
return false;
......@@ -125,13 +125,12 @@ bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, co
if (pflag == 0) continue;
// 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 dist1 = (A_m * bunch->R[i](0) + B_m * bunch->R[i](1) + C_m) / R_m; // [m]
double dist2 = dist1 * std::sqrt(1.0 + 1.0 / tangle / tangle);
double dt = dist2 / (std::sqrt(1.0 - 1.0 / (1.0 + dot(bunch->P[i], bunch->P[i]))) * Physics::c) * 1.0e9; // [ns]
strippoint(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);
strippoint(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);
strippoint(2) = bunch->R[i](2);
strippoint = strippoint*0.001; // mm->m
lossDs_m->addParticle(strippoint, bunch->P[i], bunch->ID[i], t+dt,
turnnumber, bunch->bunchNum[i]);
......