Commit 955a6f6d authored by snuverink_j's avatar snuverink_j
Browse files

fix checkCollimator, fix calculation of smallest and largest r, fixes #253

parent 1c7b2f3d
......@@ -70,6 +70,7 @@ CCollimator::CCollimator(const std::string &name):
lossDs_m(nullptr),
parmatint_m(NULL) {
setDimensions(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
setGeom();
}
......@@ -94,7 +95,7 @@ bool CCollimator::applyToReferenceParticle(const Vector_t &R, const Vector_t &P,
bool CCollimator::checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax) {
double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
double r1 = std::hypot(rmax(0), rmax(1));
bool isDead = false;
if (rmax(2) >= zstart_m && rmin(2) <= zend_m) {
if ( r1 > rstart_m - 10.0 && r1 < rend_m + 10.0 ){
......@@ -107,33 +108,33 @@ bool CCollimator::checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax) {
return isDead;
}
// rectangle collimators in cyclotron cylindrical coordinates
// without particlematterinteraction, the particle hitting collimator is deleted directly
bool CCollimator::checkCollimator(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) {
bool CCollimator::checkCollimator(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {
bool flagNeedUpdate = false;
Vector_t rmin, rmax;
bunch->get_bounds(rmin, rmax);
std::pair<Vector_t, double> boundingSphere;
boundingSphere.first = 0.5 * (rmax + rmin);
boundingSphere.second = euclidean_norm(rmax - boundingSphere.first);
if (rmax(2) >= zstart_m && rmin(2) <= zend_m) {
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 rbunch_min = sqrt(rmin(0) * rmin(0) + rmin(1) * rmin(1));
double rbunch_max = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
// if ( rbunch_max > rstart_m - 100.0 && rbunch_max < rend_m + 100.0 ){ // old check, only checks outer bound
if( rbunch_max > r_start && rbunch_min < r_end ){ // check similar to z
// interested in absolute minimum and maximum
double xmin = std::min(std::abs(rmin(0)), std::abs(rmax(0)));
double xmax = std::max(std::abs(rmin(0)), std::abs(rmax(0)));
double ymin = std::min(std::abs(rmin(1)), std::abs(rmax(1)));
double ymax = std::max(std::abs(rmin(1)), std::abs(rmax(1)));
double rbunch_min = std::hypot(xmin, ymin);
double rbunch_max = std::hypot(xmax, ymax);
if( rbunch_max > rmin_m && rbunch_min < rmax_m ){ // check similar to z
size_t tempnum = bunch->getLocalNum();
int pflag = 0;
// now check each particle in bunch
for (unsigned int i = 0; i < tempnum; ++i) {
if (bunch->PType[i] == ParticleType::REGULAR && bunch->R[i](2) < zend_m && bunch->R[i](2) > zstart_m ) {
// only now careful check in r
pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
/// bunch->Bin[i] != -1 makes sure the particle is not stored in more than one collimator
if ((pflag != 0) && (bunch->Bin[i] != -1)) {
if ((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;
......@@ -145,6 +146,9 @@ bool CCollimator::checkCollimator(PartBunchBase<double, 3> *bunch, const int tur
}
reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
if (flagNeedUpdate && parmatint_m) {
std::pair<Vector_t, double> boundingSphere;
boundingSphere.first = 0.5 * (rmax + rmin);
boundingSphere.second = euclidean_norm(rmax - boundingSphere.first);
parmatint_m->apply(bunch, boundingSphere);
}
return flagNeedUpdate;
......@@ -225,15 +229,18 @@ void CCollimator::setDimensions(double xstart, double xend, double ystart, doubl
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);
rstart_m = std::hypot(xstart, ystart);
rend_m = std::hypot(xend, yend);
// start position is the one with lowest radius
if (rstart_m > rend_m) {
std::swap(xstart_m, xend_m);
std::swap(ystart_m, yend_m);
std::swap(zstart_m, zend_m);
std::swap(rstart_m, rend_m);
}
// zstart and zend are independent from x, y
if (zstart_m > zend_m){
std::swap(zstart_m, zend_m);
}
}
bool CCollimator::bends() const {
......@@ -290,11 +297,20 @@ void CCollimator::setGeom() {
geom_m[4].x = geom_m[0].x;
geom_m[4].y = geom_m[0].y;
if (zstart_m > zend_m){
double tempz = zstart_m;
zstart_m = zend_m;
zend_m = tempz;
// calculate maximum and mininum r from these
// current implementation not perfect minimum does not need to lie on corner, on in middle
rmin_m = std::hypot(xstart_m, ystart_m);
rmax_m = -std::numeric_limits<double>::max();
for (int i=0; i<4; i++) {
double r = std::hypot(geom_m[i].x, geom_m[i].y);
rmin_m = std::min(rmin_m, r);
rmax_m = std::max(rmax_m, r);
}
// some debug printout
// for (int i=0; i<4; i++) {
// *gmsg << "point " << i << " ( " << geom_m[i].x << ", " << geom_m[i].y << ")" << endl;
// }
// *gmsg << "rmin " << rmin_m << " rmax " << rmax_m << endl;
}
......@@ -302,8 +318,8 @@ int CCollimator::checkPoint(const double &x, const double &y) {
int cn = 0;
for (int i = 0; i < 4; i++) {
if (((geom_m[i].y <= y) && (geom_m[i+1].y > y))
|| ((geom_m[i].y > y) && (geom_m[i+1].y <= y))) {
if (((geom_m[i].y <= y) && (geom_m[i+1].y > y)) ||
((geom_m[i].y > y) && (geom_m[i+1].y <= y))) {
float vt = (float)(y - geom_m[i].y) / (geom_m[i+1].y - geom_m[i].y);
if (x < geom_m[i].x + vt * (geom_m[i+1].x - geom_m[i].x))
......
......@@ -121,8 +121,11 @@ private:
double rend_m;
double width_m;
///@}
/// 4 end points in (x,y) (5th point == 1st)
Point geom_m[5];
double rmin_m; ///< minimum extend in r
double rmax_m; ///< maximum extend in r
/// Sets geom_m and maximal radii
void setGeom();
unsigned int losses_m;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment