Commit c5bb34a4 authored by snuverink_j's avatar snuverink_j
Browse files

Resolve "Reduce code duplication in CCollimator, Probe, Septum and Stripper with common base class"

parent 0d91ba47
......@@ -44,6 +44,7 @@
#include "AbsBeamline/MultipoleTStraight.h"
#include "AbsBeamline/MultipoleTCurvedConstRadius.h"
#include "AbsBeamline/MultipoleTCurvedVarRadius.h"
#include "AbsBeamline/PluginElement.h"
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RFCavity.h"
......@@ -835,16 +836,16 @@ void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
Probe *elptr = dynamic_cast<Probe *>(prob.clone());
myElements.push_back(elptr);
double xstart = elptr->getXstart();
double xstart = elptr->getXStart();
*gmsg << "XStart= " << xstart << " [mm]" << endl;
double xend = elptr->getXend();
double xend = elptr->getXEnd();
*gmsg << "XEnd= " << xend << " [mm]" << endl;
double ystart = elptr->getYstart();
double ystart = elptr->getYStart();
*gmsg << "YStart= " << ystart << " [mm]" << endl;
double yend = elptr->getYend();
double yend = elptr->getYEnd();
*gmsg << "YEnd= " << yend << " [mm]" << endl;
// initialise, do nothing
......@@ -1043,16 +1044,16 @@ void ParallelCyclotronTracker::visitSeptum(const Septum &sept) {
Septum *elptr = dynamic_cast<Septum *>(sept.clone());
myElements.push_back(elptr);
double xstart = elptr->getXstart();
double xstart = elptr->getXStart();
*gmsg << "XStart = " << xstart << " [mm]" << endl;
double xend = elptr->getXend();
double xend = elptr->getXEnd();
*gmsg << "XEnd = " << xend << " [mm]" << endl;
double ystart = elptr->getYstart();
double ystart = elptr->getYStart();
*gmsg << "YStart = " << ystart << " [mm]" << endl;
double yend = elptr->getYend();
double yend = elptr->getYEnd();
*gmsg << "YEnd = " << yend << " [mm]" << endl;
double width = elptr->getWidth();
......@@ -1135,16 +1136,16 @@ void ParallelCyclotronTracker::visitStripper(const Stripper &stripper) {
Stripper *elptr = dynamic_cast<Stripper *>(stripper.clone());
myElements.push_back(elptr);
double xstart = elptr->getXstart();
double xstart = elptr->getXStart();
*gmsg << "XStart= " << xstart << " [mm]" << endl;
double xend = elptr->getXend();
double xend = elptr->getXEnd();
*gmsg << "XEnd= " << xend << " [mm]" << endl;
double ystart = elptr->getYstart();
double ystart = elptr->getYStart();
*gmsg << "YStart= " << ystart << " [mm]" << endl;
double yend = elptr->getYend();
double yend = elptr->getYEnd();
*gmsg << "YEnd= " << yend << " [mm]" << endl;
double opcharge = elptr->getOPCharge();
......@@ -2256,35 +2257,25 @@ void ParallelCyclotronTracker::borisExternalFields(double h) {
void ParallelCyclotronTracker::applyPluginElements(const double dt) {
// 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::SEPTUM) {
(static_cast<Septum *>(((*sindex)->second).second))->checkSeptum(itsBunch_m);
}
if(((*sindex)->first) == ElementBase::PROBE) {
(static_cast<Probe *>(((*sindex)->second).second))->checkProbe(itsBunch_m,
turnnumber_m,
itsBunch_m->getT() * 1e9 /*[ns]*/,
dt);
}
if(((*sindex)->first) == ElementBase::STRIPPER) {
bool flag_stripper = (static_cast<Stripper *>(((*sindex)->second).second))
-> checkStripper(itsBunch_m, turnnumber_m, itsBunch_m->getT() * 1e9 /*[ns]*/, dt);
if(flag_stripper) {
ElementBase::ElementType type = ((*sindex)->first);
if( type == ElementBase::CCOLLIMATOR ||
type == ElementBase::PROBE ||
type == ElementBase::SEPTUM ||
type == ElementBase::STRIPPER) {
PluginElement* element = static_cast<PluginElement *>(((*sindex)->second).second);
bool flag = element->check(itsBunch_m,
turnnumber_m,
itsBunch_m->getT() * 1e9 /*[ns]*/,
dt);
if (type == ElementBase::STRIPPER && flag == true) {
itsBunch_m->updateNumTotal();
*gmsg << "* Total number of particles after stripping = " << itsBunch_m->getTotalNum() << endl;
}
}
if(((*sindex)->first) == ElementBase::CCOLLIMATOR) {
CCollimator * collim;
collim = static_cast<CCollimator *>(((*sindex)->second).second);
collim->checkCollimator(itsBunch_m, turnnumber_m, itsBunch_m->getT() * 1e9 /*[ns]*/, dt);
}
}
itsBunch_m->R *= Vector_t(0.001);
......
// ------------------------------------------------------------------------
// $RCSfile: CCollimator.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: CCollimator
// Defines the abstract interface for a beam collimator for cyclotrons.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
//
// ------------------------------------------------------------------------
#include "AbsBeamline/CCollimator.h"
#include "AbsBeamline/BeamlineVisitor.h"
......@@ -25,93 +5,40 @@
#include "Fields/Fieldmap.h"
#include "Solvers/ParticleMatterInteractionHandler.hh"
#include "Structure/LossDataSink.h"
#include "Utilities/Options.h"
#include <cmath>
#include <fstream>
#include <memory>
extern Inform *gmsg;
// Class Collimator
// ------------------------------------------------------------------------
CCollimator::CCollimator():CCollimator("")
{}
CCollimator::CCollimator():
Component(),
filename_m(""),
informed_m(false),
losses_m(0),
lossDs_m(nullptr),
parmatint_m(NULL) {
CCollimator::CCollimator(const std::string &name):
PluginElement(name) {
setDimensions(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
setGeom(0.0);
}
CCollimator::CCollimator(const CCollimator &right):
Component(right),
filename_m(right.filename_m),
informed_m(right.informed_m),
losses_m(0),
lossDs_m(nullptr),
parmatint_m(NULL)
{
PluginElement(right),
informed_m(right.informed_m) {
setDimensions(right.xstart_m, right.xend_m,
right.ystart_m, right.yend_m,
right.zstart_m, right.zend_m,
right.width_m);
setGeom();
}
CCollimator::CCollimator(const std::string &name):
Component(name),
filename_m(""),
informed_m(false),
losses_m(0),
lossDs_m(nullptr),
parmatint_m(NULL) {
setDimensions(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
setGeom();
}
CCollimator::~CCollimator() {
if (online_m)
goOffline();
setGeom(width_m);
}
CCollimator::~CCollimator() {}
void CCollimator::accept(BeamlineVisitor &visitor) const {
visitor.visitCCollimator(*this);
}
bool CCollimator::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
return false;
}
bool CCollimator::applyToReferenceParticle(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B) {
return false;
}
bool CCollimator::checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax) {
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 ){
if (r(2) < zend_m && r(2) > zstart_m ) {
int pflag = checkPoint(r(0), r(1));
isDead = (pflag != 0);
}
}
}
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::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {
bool flagNeedUpdate = false;
Vector_t rmin, rmax;
......@@ -155,36 +82,12 @@ bool CCollimator::checkCollimator(PartBunchBase<double, 3> *bunch, const int /*t
return flagNeedUpdate;
}
void CCollimator::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) {
endField = startField + getElementLength();
initialise(bunch);
}
void CCollimator::initialise(PartBunchBase<double, 3> *bunch) {
RefPartBunch_m = bunch;
void CCollimator::doInitialise(PartBunchBase<double, 3> *bunch) {
parmatint_m = getParticleMatterInteraction();
// if (!parmatint_m) {
if (filename_m == std::string(""))
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
else
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
// }
goOnline(-1e6);
}
void CCollimator::finalise()
{
if (online_m)
goOffline();
*gmsg << "* Finalize cyclotron collimator " << getName() << endl;
}
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()));
......@@ -196,8 +99,12 @@ void CCollimator::goOnline(const double &) {
online_m = true;
}
void CCollimator::doFinalise() {
*gmsg << "* Finalize cyclotron collimator " << getName() << endl;
}
void CCollimator::print() {
if (RefPartBunch_m == NULL) {
if (RefPartBunch_m == nullptr) {
if (!informed_m) {
std::string errormsg = Fieldmap::typeset_msg("BUNCH SIZE NOT SET", "warning");
ERRORMSG(errormsg << endl);
......@@ -215,48 +122,16 @@ void CCollimator::print() {
<< "R start " << xend_m << " (mm) R rend " << yend_m << " (mm)" << endl;
}
void CCollimator::goOffline() {
if (online_m && lossDs_m)
lossDs_m->save();
lossDs_m.reset(0);
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;
setDimensions(xstart, xend, ystart, yend);
zstart_m = zstart;
zend_m = zend;
width_m = width;
rstart_m = std::hypot(xstart, ystart);
rend_m = std::hypot(xend, yend);
// start position is the one with lowest radius
if (rstart_m > rend_m) {
std::swap(xstart_m, xend_m);
std::swap(ystart_m, yend_m);
std::swap(rstart_m, rend_m);
}
// zstart and zend are independent from x, y
if (zstart_m > zend_m){
std::swap(zstart_m, zend_m);
}
}
bool CCollimator::bends() const {
return false;
}
void CCollimator::setOutputFN(std::string fn) {
filename_m = fn;
}
std::string CCollimator::getOutputFN() {
if (filename_m == std::string(""))
return getName();
else
return filename_m.substr(0, filename_m.rfind("."));
setGeom(width_m);
}
void CCollimator::getDimensions(double &zBegin, double &zEnd) const {
......@@ -268,36 +143,7 @@ ElementBase::ElementType CCollimator::getType() const {
return CCOLLIMATOR;
}
std::string CCollimator::getCollimatorShape() {
return "CCollimator";
}
void CCollimator::setGeom() {
double slope;
if (xend_m == xstart_m)
slope = 1.0e12;
else
slope = (yend_m - ystart_m) / (xend_m - xstart_m);
double coeff2 = sqrt(1 + slope * slope);
double coeff1 = slope / coeff2;
double halfdist = width_m / 2.0;
geom_m[0].x = xstart_m - halfdist * coeff1;
geom_m[0].y = ystart_m + halfdist / coeff2;
geom_m[1].x = xstart_m + halfdist * coeff1;
geom_m[1].y = ystart_m - halfdist / coeff2;
geom_m[2].x = xend_m + halfdist * coeff1;
geom_m[2].y = yend_m - halfdist / coeff2;
geom_m[3].x = xend_m - halfdist * coeff1;
geom_m[3].y = yend_m + halfdist / coeff2;
geom_m[4].x = geom_m[0].x;
geom_m[4].y = geom_m[0].y;
void CCollimator::doSetGeom() {
// calculate maximum and mininum r from these
// current implementation not perfect minimum does not need to lie on corner, or in middle
rmin_m = std::hypot(xstart_m, ystart_m);
......@@ -314,18 +160,17 @@ void CCollimator::setGeom() {
// *gmsg << "rmin " << rmin_m << " rmax " << rmax_m << endl;
}
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))) {
float vt = (float)(y - geom_m[i].y) / (geom_m[i+1].y - geom_m[i].y);
if (x < geom_m[i].x + vt * (geom_m[i+1].x - geom_m[i].x))
++cn;
}
}
return (cn & 1); // 0 if even (out), and 1 if odd (in)
}
\ No newline at end of file
// bool CCollimator::checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax) {
// 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 ){
// if (r(2) < zend_m && r(2) > zstart_m ) {
// int pflag = checkPoint(r(0), r(1));
// isDead = (pflag != 0);
// }
// }
// }
// return isDead;
// }
\ No newline at end of file
#ifndef CLASSIC_CCollimator_HH
#define CLASSIC_CCollimator_HH
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: CCollimator
// Defines the abstract interface for a beam CCollimator.
// *** MISSING *** CCollimator interface is still incomplete.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
//
// ------------------------------------------------------------------------
#include "AbsBeamline/Component.h"
#include "AbsBeamline/PluginElement.h"
#include <memory>
#include <string>
class BeamlineVisitor;
class LossDataSink;
class ParticleMatterInteractionHandler;
// Class CCollimator
// ------------------------------------------------------------------------
/// Abstract collimator.
/// Interface for cyclotron collimator.
// Class CCollimator defines the abstract interface for a collimator.
class CCollimator: public Component {
class CCollimator: public PluginElement {
public:
/// Constructor with given name.
explicit CCollimator(const std::string &name);
CCollimator();
CCollimator(const CCollimator &rhs);
virtual ~CCollimator();
void operator=(const CCollimator &) = delete;
virtual ~CCollimator();
/// Apply visitor to CCollimator.
virtual void accept(BeamlineVisitor &) const;
///@{ Override implementation of PluginElement
virtual void goOnline(const double &kineticEnergy) override;
virtual ElementBase::ElementType getType() const override;
virtual void getDimensions(double &zBegin, double &zEnd) const override;
///@}
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool applyToReferenceParticle(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B);
virtual bool checkCollimator(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep);
virtual bool checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax);
virtual void initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField);
virtual void initialise(PartBunchBase<double, 3> *bunch);
virtual void finalise();
virtual bool bends() const;
virtual void goOnline(const double &kineticEnergy);
virtual void goOffline();
virtual ElementBase::ElementType getType() const;
virtual void getDimensions(double &zBegin, double &zEnd) const;
/// unused check method
// bool checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax);
/// Some debug print
void print();
std::string getCollimatorShape();
void setOutputFN(std::string fn);
std::string getOutputFN();
unsigned int getLosses() const;
// --------Cyclotron collimator
/// Set dimensions and consistency checks
void setDimensions(double xstart, double xend,
double ystart, double yend,
double zstart, double zend,
double width);
/// unhide PluginElement::setDimensions(double xstart, double xend, double ystart, double yend)
using PluginElement::setDimensions;
double getXStart() ;
double getYStart() ;
///@{ Member variable access
double getZStart() ;
double getXEnd() ;
double getYEnd() ;
double getZEnd() ;
double getWidth() ;
int checkPoint(const double & x, const double & y );
///@}
private:
/// Initialise particle matter interaction
virtual void doInitialise(PartBunchBase<double, 3> *bunch) override;
/// Record hits when bunch particles pass
virtual bool doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) override;
/// Calculate extend in r
virtual void doSetGeom() override;
/// Virtual hook for finalise
virtual void doFinalise() override;
std::string filename_m; /**< The name of the outputfile*/
bool informed_m;
bool informed_m = false; ///< Flag if error information already printed
///@{ 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;
///@}
/// 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;
std::unique_ptr<LossDataSink> lossDs_m;
ParticleMatterInteractionHandler *parmatint_m;
ParticleMatterInteractionHandler *parmatint_m = nullptr;
};
inline
unsigned int CCollimator::getLosses() const {
return losses_m;
}
inline
double CCollimator::getXStart() {
return xstart_m;
}
inline
double CCollimator::getXEnd() {
return xend_m;
}
inline
double CCollimator::getYStart() {
return ystart_m;
}
inline
double CCollimator::getYEnd() {
return yend_m;