Commit 3695de09 authored by snuverink_j's avatar snuverink_j
Browse files

fix #206: require start position to have lower radius

parent b52d004c
......@@ -39,32 +39,25 @@ CCollimator::CCollimator():
Component(),
filename_m(""),
informed_m(false),
xstart_m(0.0),
xend_m(0.0),
ystart_m(0.0),
yend_m(0.0),
width_m(0.0),
losses_m(0),
lossDs_m(nullptr),
parmatint_m(NULL)
{}
parmatint_m(NULL) {
setDimensions(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
}
CCollimator::CCollimator(const CCollimator &right):
Component(right),
filename_m(right.filename_m),
informed_m(right.informed_m),
xstart_m(right.xstart_m),
xend_m(right.xend_m),
ystart_m(right.ystart_m),
yend_m(right.yend_m),
zstart_m(right.zstart_m),
zend_m(right.zend_m),
width_m(right.width_m),
losses_m(0),
lossDs_m(nullptr),
parmatint_m(NULL)
{
setDimensions(right.xstart_m, right.xend_m,
right.ystart_m, right.yend_m,
right.zstart_m, right.zend_m,
right.width_m);
setGeom();
}
......@@ -73,17 +66,11 @@ CCollimator::CCollimator(const std::string &name):
Component(name),
filename_m(""),
informed_m(false),
xstart_m(0.0),
xend_m(0.0),
ystart_m(0.0),
yend_m(0.0),
zstart_m(0.0),
zend_m(0.0),
width_m(0.0),
losses_m(0),
lossDs_m(nullptr),
parmatint_m(NULL)
{}
parmatint_m(NULL) {
setDimensions(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
}
CCollimator::~CCollimator() {
......@@ -107,12 +94,10 @@ bool CCollimator::applyToReferenceParticle(const Vector_t &R, const Vector_t &P,
bool CCollimator::checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax) {
double r_start = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
bool isDead = false;
if (rmax(2) >= zstart_m && rmin(2) <= zend_m) {
if ( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
if ( r1 > rstart_m - 10.0 && r1 < rend_m + 10.0 ){
if (r(2) < zend_m && r(2) > zstart_m ) {
int pflag = checkPoint(r(0), r(1));
isDead = (pflag != 0);
......@@ -131,27 +116,25 @@ bool CCollimator::checkCollimator(PartBunchBase<double, 3> *bunch, const int tur
Vector_t rmin, rmax;
bunch->get_bounds(rmin, rmax);
double r_start = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
std::pair<Vector_t, double> boundingSphere;
boundingSphere.first = 0.5 * (rmax + rmin);
boundingSphere.second = euclidean_norm(rmax - boundingSphere.first);
if (rmax(2) >= zstart_m && rmin(2) <= zend_m) {
// if ( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
if ( r1 > r_start - 100.0 && r1 < r_end + 100.0 ){
// if ( r1 > rstart_m - 10.0 && r1 < rend_m + 10.0 ){
if ( r1 > rstart_m - 100.0 && r1 < rend_m + 100.0 ){
size_t tempnum = bunch->getLocalNum();
int pflag = 0;
for (unsigned int i = 0; i < tempnum; ++i) {
if (bunch->PType[i] == ParticleType::REGULAR && bunch->R[i](2) < zend_m && bunch->R[i](2) > zstart_m ) {
pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
/// bunch->Bin[i] != -1 makes sure the partcile is not stored in more than one collimator
/// bunch->Bin[i] != -1 makes sure the partcile is not stored in more than one collimator
if ((pflag != 0) && (bunch->Bin[i] != -1)) {
if (!parmatint_m)
lossDs_m->addParticle(bunch->R[i], bunch->P[i], bunch->ID[i]);
bunch->Bin[i] = -1;
flagNeedUpdate = true;
if (!parmatint_m)
lossDs_m->addParticle(bunch->R[i], bunch->P[i], bunch->ID[i]);
bunch->Bin[i] = -1;
flagNeedUpdate = true;
}
}
}
......@@ -231,6 +214,25 @@ void CCollimator::goOffline() {
online_m = false;
}
void CCollimator::setDimensions(double xstart, double xend, double ystart, double yend, double zstart, double zend, double width) {
xstart_m = xstart;
ystart_m = ystart;
xend_m = xend;
yend_m = yend;
zstart_m = zstart;
zend_m = zend;
width_m = width;
rstart_m = std::sqrt(xstart*xstart + ystart * ystart);
rend_m = std::sqrt(xend * xend + yend * yend);
// start position is the one with lowest radius
if (rstart_m > rend_m) {
std::swap(xstart_m, xend_m);
std::swap(ystart_m, yend_m);
std::swap(zstart_m, zend_m);
std::swap(rstart_m, rend_m);
}
}
bool CCollimator::bends() const {
return false;
}
......
......@@ -83,31 +83,13 @@ public:
unsigned int getLosses() const;
// void setXsize(double a) ;
// void setYsize(double b) ;
// void setXpos(double x0) ;
// void setYpos(double y0) ;
// double getXsize(double a) ;
// double getYsize(double b) ;
// double getXpos() ;
// double getYpos() ;
// --------Cyclotron collimator
void setXStart(double xstart) ;
void setYStart(double ystart) ;
void setZStart(double zstart) ;
void setXEnd(double xend) ;
void setYEnd(double yend) ;
void setZEnd(double zend) ;
void setWidth(double width) ;
/// Set dimensions and consistency checks
void setDimensions(double xstart, double xend,
double ystart, double yend,
double zstart, double zend,
double width);
double getXStart() ;
double getYStart() ;
......@@ -128,14 +110,17 @@ private:
bool informed_m;
//parameters for CCollimator
///@{ input geometry positions
double xstart_m;
double xend_m;
double ystart_m;
double yend_m;
double zstart_m;
double zend_m;
double rstart_m;
double rend_m;
double width_m;
///@}
Point geom_m[5];
void setGeom();
......@@ -152,41 +137,6 @@ unsigned int CCollimator::getLosses() const {
return losses_m;
}
inline
void CCollimator::setXStart(double xstart) {
xstart_m = xstart;
}
inline
void CCollimator::setXEnd(double xend) {
xend_m = xend;
}
inline
void CCollimator::setYStart(double ystart) {
ystart_m = ystart;
}
inline
void CCollimator::setYEnd(double yend) {
yend_m = yend;
}
inline
void CCollimator::setZStart(double zstart) {
zstart_m = zstart;
}
inline
void CCollimator::setZEnd(double zend) {
zend_m = zend;
}
inline
void CCollimator::setWidth(double width) {
width_m = width;
}
inline
double CCollimator::getXStart() {
return xstart_m;
......
......@@ -28,7 +28,6 @@
#include "Utilities/Options.h"
#include <iostream>
#include <fstream>
using Physics::pi;
extern Inform *gmsg;
......@@ -39,16 +38,8 @@ Probe::Probe():
Component(),
filename_m(""),
position_m(0.0),
xstart_m(0.0),
xend_m(0.0),
ystart_m(0.0),
yend_m(0.0),
width_m(0.0),
step_m(0){
A_m = yend_m - ystart_m;
B_m = xstart_m - xend_m;
R_m = sqrt(A_m*A_m+B_m*B_m);
C_m = ystart_m*xend_m - xstart_m*yend_m;
step_m(0.0){
setDimensions(0.0, 0.0, 0.0, 0.0, 0.0);
}
......@@ -56,16 +47,8 @@ Probe::Probe(const Probe &right):
Component(right),
filename_m(right.filename_m),
position_m(right.position_m),
xstart_m(right.xstart_m),
xend_m(right.xend_m),
ystart_m(right.ystart_m),
yend_m(right.yend_m),
width_m(right.width_m),
step_m(right.step_m){
A_m = yend_m - ystart_m;
B_m = xstart_m - xend_m;
R_m = sqrt(A_m*A_m+B_m*B_m);
C_m = ystart_m*xend_m - xstart_m*yend_m;
setDimensions(right.xstart_m, right.xend_m, right.ystart_m, right.yend_m, right.width_m);
}
......@@ -73,16 +56,8 @@ Probe::Probe(const std::string &name):
Component(name),
filename_m(""),
position_m(0.0),
xstart_m(0.0),
xend_m(0.0),
ystart_m(0.0),
yend_m(0.0),
width_m(0.0),
step_m(0){
A_m = yend_m - ystart_m;
B_m = xstart_m - xend_m;
R_m = sqrt(A_m*A_m+B_m*B_m);
C_m = ystart_m*xend_m - xstart_m*yend_m;
step_m(0.0){
setDimensions(0.0, 0.0, 0.0, 0.0, 0.0);
}
Probe::~Probe() {}
......@@ -100,9 +75,6 @@ void Probe::initialise(PartBunchBase<double, 3> *bunch, double &startField, doub
void Probe::initialise(PartBunchBase<double, 3> *bunch) {
*gmsg << "* Initialize probe" << endl;
double r_start = std::sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
double r_end = std::sqrt(xend_m * xend_m + yend_m * yend_m);
std::string name;
if (filename_m == std::string(""))
name = getName();
......@@ -110,7 +82,7 @@ void Probe::initialise(PartBunchBase<double, 3> *bunch) {
name = filename_m.substr(0, filename_m.rfind("."));
bool singlemode = (bunch->getTotalNum() == 1) ? true : false;
peakfinder_m = std::unique_ptr<PeakFinder> (new PeakFinder (name, r_start, r_end, step_m, singlemode));
peakfinder_m = std::unique_ptr<PeakFinder> (new PeakFinder (name, rstart_m, rend_m, step_m, singlemode));
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(name, !Options::asciidump));
}
......@@ -131,24 +103,25 @@ void Probe::goOffline() {
lossDs_m->save();
}
void Probe::setXstart(double xstart) {
void Probe::setDimensions(double xstart, double xend, double ystart, double yend, double width) {
xstart_m = xstart;
}
void Probe::setXend(double xend) {
xend_m = xend;
}
void Probe::setYstart(double ystart) {
ystart_m = ystart;
}
void Probe::setYend(double yend) {
yend_m = yend;
}
xend_m = xend;
yend_m = yend;
width_m = width;
rstart_m = std::sqrt(xstart*xstart + ystart * ystart);
rend_m = std::sqrt(xend * xend + yend * yend);
// start position is the one with lowest radius
if (rstart_m > rend_m) {
std::swap(xstart_m, xend_m);
std::swap(ystart_m, yend_m);
std::swap(rstart_m, rend_m);
}
void Probe::setWidth(double width) {
width_m = width;
A_m = yend_m - ystart_m;
B_m = xstart_m - xend_m;
R_m = std::sqrt(A_m*A_m+B_m*B_m);
C_m = ystart_m*xend_m - xstart_m*yend_m;
}
void Probe::setStep(double step) {
......@@ -211,12 +184,10 @@ bool Probe::checkProbe(PartBunchBase<double, 3> *bunch, const int turnnumber, c
bool flagprobed = false;
Vector_t rmin, rmax, probepoint;
bunch->get_bounds(rmin, rmax);
double r_start = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
double r2 = sqrt(rmin(0) * rmin(0) + rmin(1) * rmin(1));
if( r1 > r_start - 10.0 && r2 < r_end + 10.0 ) {
if( r1 > rstart_m - 10.0 && r2 < rend_m + 10.0 ) {
size_t tempnum = bunch->getLocalNum();
int pflag = 0;
......@@ -311,7 +282,7 @@ int Probe::checkPoint(const double &x, const double &y) {
void Probe::getDimensions(double &zBegin, double &zEnd) const {
zBegin = position_m - 0.005;
zEnd = position_m + 0.005;
zEnd = position_m + 0.005;
}
ElementBase::ElementType Probe::getType() const {
......
......@@ -59,13 +59,9 @@ public:
virtual void goOffline();
void setXstart(double xstart);
void setXend(double xend);
/// Set dimensions and consistency checks
void setDimensions(double xstart, double xend, double ystart, double yend, double width);
void setYstart(double ystart);
void setYend(double yend);
void setWidth(double width);
void setStep(double step);
virtual double getXstart() const;
......@@ -90,6 +86,8 @@ private:
double xend_m;
double ystart_m;
double yend_m;
double rstart_m;
double rend_m;
///@}
double width_m; ///< bin width, not used
Point geom_m[5]; ///< actual geometry positions with adaptive width such that each particle hits probe once per turn
......
......@@ -23,9 +23,6 @@
#include "Algorithms/PartBunchBase.h"
#include "Physics/Physics.h"
#include "Structure/LossDataSink.h" // OPAL file
#include <iostream>
#include <fstream>
using Physics::pi;
extern Inform *gmsg;
......@@ -35,37 +32,24 @@ extern Inform *gmsg;
Septum::Septum():
Component(),
filename_m(""),
position_m(0.0),
xstart_m(0.0),
xend_m(0.0),
ystart_m(0.0),
yend_m(0.0),
width_m(0.0)
{}
position_m(0.0) {
setDimensions(0.0, 0.0, 0.0, 0.0, 0.0);
}
Septum::Septum(const Septum &right):
Component(right),
filename_m(right.filename_m),
position_m(right.position_m),
xstart_m(right.xstart_m),
xend_m(right.xend_m),
ystart_m(right.ystart_m),
yend_m(right.yend_m),
width_m(right.width_m)
{}
position_m(right.position_m) {
setDimensions(right.xstart_m, right.xend_m, right.ystart_m, right.yend_m, right.width_m);}
Septum::Septum(const std::string &name):
Component(name),
filename_m(""),
position_m(0.0),
xstart_m(0.0),
xend_m(0.0),
ystart_m(0.0),
yend_m(0.0),
width_m(0.0)
{}
position_m(0.0) {
setDimensions(0.0, 0.0, 0.0, 0.0, 0.0);
}
Septum::~Septum() {
......@@ -77,15 +61,14 @@ void Septum::accept(BeamlineVisitor &visitor) const {
}
void Septum::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) {
RefPartBunch_m = bunch;
position_m = startField;
startField -= 0.005;
endField = position_m + 0.005;
initialise(bunch);
}
void Septum::initialise(PartBunchBase<double, 3> *bunch) {
RefPartBunch_m = bunch;
*gmsg << "Septum initialise" << endl;
}
......@@ -100,27 +83,22 @@ void Septum::goOffline() {
online_m = false;
}
void Septum::setXstart(double xstart) {
void Septum::setDimensions(double xstart, double xend, double ystart, double yend, double width) {
xstart_m = xstart;
}
void Septum::setXend(double xend) {
xend_m = xend;
}
void Septum::setYstart(double ystart) {
ystart_m = ystart;
xend_m = xend;
yend_m = yend;
width_m = width;
rstart_m = std::sqrt(xstart*xstart + ystart * ystart);
rend_m = std::sqrt(xend * xend + yend * yend);
// start position is the one with lowest radius
if (rstart_m > rend_m) {
std::swap(xstart_m, xend_m);
std::swap(ystart_m, yend_m);
std::swap(rstart_m, rend_m);
}
}
void Septum::setYend(double yend) {
yend_m = yend;
}
void Septum::setWidth(double width) {
width_m = width;
}
double Septum::getXstart() const {
return xstart_m;
......@@ -152,7 +130,7 @@ bool Septum::checkSeptum(PartBunchBase<double, 3> *bunch) {
Vector_t rmax;
bunch->get_bounds(rmin, rmax);
double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
if(r1 > sqrt(xstart_m * xstart_m + ystart_m * ystart_m) - 100) {
if(r1 > rstart_m - 100) {
for(unsigned int i = 0; i < bunch->getLocalNum(); ++i) {
Vector_t R = bunch->R[i];
double slope = (yend_m - ystart_m) / (xend_m - xstart_m);
......@@ -163,16 +141,12 @@ bool Septum::checkSeptum(PartBunchBase<double, 3> *bunch) {
double line1 = fabs(slope * R(0) + intcept1);
double line2 = fabs(slope * R(0) + intcept2);
if(fabs(R(1)) > line2 && fabs(R(1)) < line1 && R(0) > xstart_m && R(0) < xend_m && R(1) > ystart_m && R(1) < yend_m) {
bunch->lossDs_m->addParticle(R, bunch->P[i], bunch->ID[i]);
bunch->Bin[i] = -1;
flag = true;
}
}
}
return flag;
......
......@@ -52,15 +52,8 @@ public:
virtual void goOffline();
void setXstart(double xstart);
void setXend(double xend);
void setYstart(double ystart);
void setYend(double yend);
void setWidth(double width);
/// Set dimensions and consistency checks
void setDimensions(double xstart, double xend, double ystart, double yend, double width);
virtual double getXstart() const;
......@@ -81,11 +74,16 @@ public:
private:
std::string filename_m; /**< The name of the inputfile*/
double position_m;
///@{ input geometry positions
double xstart_m;
double xend_m;
double ystart_m;
double yend_m;
double rstart_m;
double rend_m;
double width_m;
///@}
// Not implemented.
void operator=(const Septum &);
};
......
......@@ -27,10 +27,6 @@
#include <iostream>
#include <fstream>
using Physics::pi;
using Physics::q_e;
extern Inform *gmsg;
using namespace std;
......@@ -42,20 +38,12 @@ Stripper::Stripper():
Component(),
filename_m(""),
position_m(0.0),
xstart_m(0.0),
xend_m(0.0),
ystart_m(0.0),
yend_m(0.0),
width_m(0.0),
opcharge_m(0.0),
opmass_m(0.0),
opyield_m(1.0),
stop_m(true),
step_m(0) {
A_m = yend_m - ystart_m;
B_m = xstart_m - xend_m;
R_m = sqrt(A_m*A_m+B_m*B_m);
C_m = ystart_m*xend_m - xstart_m*yend_m;
setDimensions(0.0, 0.0, 0.0, 0.0, 0.0);
}
......@@ -63,20 +51,12 @@ Stripper::Stripper(const Stripper &right):
Component(