code improvements and cleanup

src/Classic/AbsBeamline/Degrader.cpp

@@ -51,22 +51,13 @@ Degrader::Degrader(const Degrader &right):
 {}
 
 Degrader::Degrader(const std::string &name):
-    Component(name),
-    filename_m(""),
-    PosX_m(0),
-    PosY_m(0),
-    PosZ_m(0),
-    MomentumX_m(0),
-    MomentumY_m(0),
-    MomentumZ_m(0),
-    time_m(0),
-    id_m(0)
+    Component(name)
 {}
 
 
 Degrader::~Degrader() {
 
-  if (online_m)
+    if(online_m)
         goOffline();
 }
 
@@ -128,13 +119,8 @@ bool Degrader::applyToReferenceParticle(const Vector_t &R,
 }
 
 void Degrader::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) {
-    RefPartBunch_m = bunch;
+    initialise(bunch);
     endField = startField + getElementLength();
-
-    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));
 }
 
 void Degrader::initialise(PartBunchBase<double, 3> *bunch) {
@@ -154,14 +140,16 @@ void Degrader::finalise()
 void Degrader::goOnline(const double &) {
     Inform msg("Degrader::goOnline ");
 
-    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()));
+    int maximumSize = (int)(1.1 * RefPartBunch_m->getLocalNum());
+
+    PosX_m.reserve(maximumSize);
+    PosY_m.reserve(maximumSize);
+    PosZ_m.reserve(maximumSize);
+    MomentumX_m.reserve(maximumSize);
+    MomentumY_m.reserve(maximumSize);
+    MomentumZ_m.reserve(maximumSize);
+    time_m.reserve(maximumSize);
+    id_m.reserve(maximumSize);
     online_m = true;
 }
 

src/Classic/Algorithms/Tracker.cpp

@@ -81,10 +81,8 @@ typedef FTps<double, 6> Series;
 
 Tracker::Tracker(const Beamline &beamline, const PartData &reference,
                  bool backBeam, bool backTrack):
-    AbstractTracker(beamline, reference, backBeam, backTrack),
-    itsBeamline_m(beamline),
-    itsBunch_m(nullptr)
-{ }
+    Tracker(beamline, nullptr, reference, backBeam, backTrack)
+{}
 
 
 Tracker::Tracker(const Beamline &beamline,

src/Classic/Solvers/BeamStrippingPhysics.cpp

@@ -635,12 +635,6 @@ bool BeamStrippingPhysics::stillActive() {
     return locPartsInMat_m != 0;
 }
 
-bool BeamStrippingPhysics::stillAlive(PartBunchBase<double, 3> */*bunch*/) {
-    bool beamstrippingAlive = true;
-    return beamstrippingAlive;
-}
-
-
 /*
     Cross sections parameters for interaction with air
     -- [1] -> Nitrogen

src/Classic/Solvers/BeamStrippingPhysics.hh

@@ -52,15 +52,14 @@ public:
                        const std::pair<Vector_t, double> &boundingSphere);
 
     virtual const std::string getType() const;
-    void print(Inform& msg);
-    bool stillActive();
-    bool stillAlive(PartBunchBase<double, 3> *bunch);
-
-    inline double getTime() {return T_m;}
-    std::string getName() {return element_ref_m->getName();}
-    size_t getParticlesInMat() {return locPartsInMat_m;}
-    unsigned getRediffused() {return rediffusedStat_m;}
-    unsigned int getNumEntered() {return bunchToMatStat_m;}
+    virtual void print(Inform& msg);
+    virtual bool stillActive();
+
+    virtual inline double getTime() {return T_m;}
+    virtual std::string getName() {return element_ref_m->getName();}
+    virtual size_t getParticlesInMat() {return locPartsInMat_m;}
+    virtual unsigned getRediffused() {return rediffusedStat_m;}
+    virtual unsigned int getNumEntered() {return bunchToMatStat_m;}
     inline void doPhysics(PartBunchBase<double, 3> *bunch);
 
 private:

src/Classic/Solvers/CollimatorPhysics.cpp

@@ -303,7 +303,7 @@ bool CollimatorPhysics::computeEnergyLoss(Vector_t &P,
     double gamma = Util::getGamma(P);
     const double gammaSqr = std::pow(gamma, 2);
     const double betaSqr = 1.0 - 1.0 / gammaSqr;
-    double beta = sqrt(betaSqr);
+    double beta = std::sqrt(betaSqr);
     double Ekin = (gamma - 1) * massProton_keV;
     double dEdx = 0.0;
 
@@ -331,12 +331,12 @@ bool CollimatorPhysics::computeEnergyLoss(Vector_t &P,
     }
 
     if (includeFluctuations) {
-        double sigma_E = sqrt(K * massElectron_keV * rho_m * (Z_m / A_m) * deltas * m2cm);
+        double sigma_E = std::sqrt(K * massElectron_keV * rho_m * (Z_m / A_m) * deltas * m2cm);
         Ekin += gsl_ran_gaussian(rGen_m, sigma_E);
     }
 
     gamma = Ekin / massProton_keV + 1.0;
-    beta = sqrt(1.0 - 1.0 / std::pow(gamma, 2));
+    beta = std::sqrt(1.0 - 1.0 / std::pow(gamma, 2));
     P = gamma * beta * P / euclidean_norm(P);
 
     bool stopped = (Ekin < 10 || dEdx > 0);
@@ -366,10 +366,10 @@ void  CollimatorPhysics::applyRotation(Vector_t &P,
 
 void CollimatorPhysics::applyRandomRotation(Vector_t &P, double theta0) {
 
-    double thetaru = 2.5 / sqrt(gsl_rng_uniform(rGen_m)) * 2.0 * theta0;
+    double thetaru = 2.5 / std::sqrt(gsl_rng_uniform(rGen_m)) * 2.0 * theta0;
     double phiru = Physics::two_pi * gsl_rng_uniform(rGen_m);
 
-    double normPtrans = sqrt(P(0) * P(0) + P(1) * P(1));
+    double normPtrans = std::sqrt(P(0) * P(0) + P(1) * P(1));
     double Theta = std::atan(normPtrans / std::abs(P(2)));
     double CosT = cos(Theta);
     double SinT = sin(Theta);
@@ -399,10 +399,10 @@ void  CollimatorPhysics::computeCoulombScattering(Vector_t &R,
     constexpr double sqrtThreeInv = 0.57735026918962576451; // sqrt(1.0 / 3.0)
     const double normP = euclidean_norm(P);
     const double gammaSqr = std::pow(normP, 2) + 1.0;
-    const double beta = sqrt(1.0 - 1.0 / gammaSqr);
+    const double beta = std::sqrt(1.0 - 1.0 / gammaSqr);
     const double deltas = dt * beta * Physics::c;
     const double theta0 = (13.6e6 / (beta * normP * massProton_eV) *
-                           chargeProton * sqrt(deltas / X0_m) *
+                           chargeProton * std::sqrt(deltas / X0_m) *
                            (1.0 + 0.038 * log(deltas / X0_m)));
 
     double phi = Physics::two_pi * gsl_rng_uniform(rGen_m);
@@ -505,7 +505,7 @@ void CollimatorPhysics::copyFromBunch(PartBunchBase<double, 3> *bunch,
             hitTester_m->checkHit(bunch->R[i]))
         {
             // adjust the time step for those particles that enter the material
-            // such that it corresponds to the time needed to reach the curren
+            // such that it corresponds to the time needed to reach the current
             // location form the edge of the material. Only use this time step
             // for the computation of the interaction with the material, not for
             // the integration of the particles. This will ensure that the momenta
@@ -576,31 +576,6 @@ bool CollimatorPhysics::stillActive() {
     return totalPartsInMat_m != 0;
 }
 
-bool CollimatorPhysics::stillAlive(PartBunchBase<double, 3> *bunch) {
-
-    bool degraderAlive = true;
-
-    //free GPU memory in case element is degrader, it is empty and bunch has moved past it
-    if (collshape_m == ElementBase::DEGRADER && totalPartsInMat_m == 0) {
-        Degrader *deg = static_cast<Degrader *>(element_ref_m);
-
-        //get the size of the degrader
-        double zBegin, zEnd;
-        deg->getDimensions(zBegin, zEnd);
-
-        //get the average Z position of the bunch
-        Vector_t bunchOrigin = bunch->get_origin();
-
-        //if bunch has moved past degrader free GPU memory
-        if (bunchOrigin[2] > zBegin) {
-            degraderAlive = false;
-        }
-    }
-
-    return degraderAlive;
-
-}
-
 namespace {
     bool myCompF(PART x, PART y) {
       return x.label > y.label;

src/Classic/Solvers/CollimatorPhysics.hh

@@ -75,15 +75,14 @@ public:
 
     virtual const std::string getType() const;
 
-    void print(Inform& os);
-    bool stillActive();
-    bool stillAlive(PartBunchBase<double, 3> *bunch);
-
-    double getTime();
-    std::string getName();
-    size_t getParticlesInMat();
-    unsigned getRediffused();
-    unsigned int getNumEntered();
+    virtual void print(Inform& os);
+    virtual bool stillActive();
+
+    virtual double getTime();
+    virtual std::string getName();
+    virtual size_t getParticlesInMat();
+    virtual unsigned getRediffused();
+    virtual unsigned int getNumEntered();
     void computeInteraction();
 
     virtual bool computeEnergyLoss(Vector_t &P,

src/Classic/Solvers/ParticleMatterInteractionHandler.hh

@@ -19,7 +19,6 @@ public:
     virtual const std::string getType() const = 0;
     virtual void print(Inform& os) = 0;
     virtual bool stillActive() = 0;
-    virtual bool stillAlive(PartBunchBase<double, 3> *bunch) = 0;
     virtual double getTime() = 0;
     virtual std::string getName() = 0;
     virtual size_t getParticlesInMat() = 0;

src/Classic/Utilities/ClassicField.h

 #ifndef CLASSIC_FIELD_H
 #define CLASSIC_FIELD_H
 
-#include <vector>
 #include <list>
 #include <memory>
 #include "AbsBeamline/Component.h"
-#include "Algorithms/Quaternion.h"
-#include "Algorithms/CoordinateSystemTrafo.h"
 
 class ClassicField {
 public:
@@ -34,10 +31,6 @@ public:
 
     ElementBase::BoundingBox getBoundingBoxInLabCoords() const;
 
-    CoordinateSystemTrafo getCoordTransformationTo() const ;
-    void setCoordTransformationTo(const CoordinateSystemTrafo &trafo);
-    bool isPositioned() const;
-    void fixPosition();
     unsigned int order_m;
 private:
     std::shared_ptr<Component> element_m;
@@ -80,26 +73,6 @@ inline void ClassicField::setEnd(const double & z) {
     end_m = z;
 }
 
-inline
-CoordinateSystemTrafo ClassicField::getCoordTransformationTo() const {
-    return element_m->getCSTrafoGlobal2Local();
-}
-
-inline
-void ClassicField::setCoordTransformationTo(const CoordinateSystemTrafo &trafo) {
-    element_m->setCSTrafoGlobal2Local(trafo);
-}
-
-inline
-bool ClassicField::isPositioned() const {
-    return element_m->isPositioned();
-}
-
-inline
-void ClassicField::fixPosition() {
-    element_m->fixPosition();
-}
-
 inline
 ElementBase::BoundingBox ClassicField::getBoundingBoxInLabCoords() const {
     return element_m->getBoundingBoxInLabCoords();

src/Elements/OpalBeamline.cpp

@@ -53,7 +53,7 @@ std::set<std::shared_ptr<Component>> OpalBeamline::getElements(const Vector_t &x
     const FieldList::iterator end = elements_m.end();
     for (; it != end; ++ it) {
         std::shared_ptr<Component> element = (*it).getElement();
-        Vector_t r = (*it).getCoordTransformationTo().transformTo(x);
+        Vector_t r = element->getCSTrafoGlobal2Local().transformTo(x);
 
         if (element->isInside(r)) {
             elementSet.insert(element);
@@ -212,11 +212,11 @@ void OpalBeamline::compute3DLattice() {
 
         FieldList::iterator it = elements_m.begin();
         for (; it != end; ++ it) {
-            if ((*it).isPositioned()) {
+            std::shared_ptr<Component> element = (*it).getElement();
+            if (element->isPositioned()) {
                 continue;
             }
             (*it).order_m = minOrder;
-            std::shared_ptr<Component> element = (*it).getElement();
 
             if (element->getType() != ElementBase::SBEND &&
                 element->getType() != ElementBase::RBEND &&
@@ -266,7 +266,7 @@ void OpalBeamline::compute3DLattice() {
             CoordinateSystemTrafo fromEndLastToEndThis(endThis3D,
                                                        rotationAboutAxis.conjugate());
 
-            (*it).setCoordTransformationTo(fromEndLastToBeginThis * currentCoordTrafo);
+            element->setCSTrafoGlobal2Local(fromEndLastToBeginThis * currentCoordTrafo);
 
             currentCoordTrafo = (fromEndLastToEndThis * currentCoordTrafo);
 
@@ -279,11 +279,11 @@ void OpalBeamline::compute3DLattice() {
 
     FieldList::iterator it = elements_m.begin();
     for (; it != end; ++ it) {
-        if ((*it).isPositioned()) continue;
+        std::shared_ptr<Component> element = (*it).getElement();
+        if (element->isPositioned()) continue;
 
         (*it).order_m = order ++;
 
-        std::shared_ptr<Component> element = (*it).getElement();
         double beginThisPathLength = element->getElementPosition();
         double thisLength = element->getElementLength();
         Vector_t beginThis3D(0, 0, beginThisPathLength - endPriorPathLength);
@@ -344,10 +344,10 @@ void OpalBeamline::compute3DLattice() {
 
             CoordinateSystemTrafo fromLastToThis(beginThis3D, rotationAboutZ);
 
-            (*it).setCoordTransformationTo(fromLastToThis * currentCoordTrafo);
+            element->setCSTrafoGlobal2Local(fromLastToThis * currentCoordTrafo);
         }
 
-        (*it).fixPosition();
+        element->fixPosition();
     }
 }
 
@@ -376,8 +376,8 @@ void OpalBeamline::save3DLattice() {
     MeshGenerator mesh;
     for (; it != end; ++ it) {
         std::shared_ptr<Component> element = (*it).getElement();
-        CoordinateSystemTrafo toBegin = element->getEdgeToBegin() * (*it).getCoordTransformationTo();
-        CoordinateSystemTrafo toEnd = element->getEdgeToEnd() * (*it).getCoordTransformationTo();
+        CoordinateSystemTrafo toBegin = element->getEdgeToBegin() * element->getCSTrafoGlobal2Local();
+        CoordinateSystemTrafo toEnd = element->getEdgeToEnd() * element->getCSTrafoGlobal2Local();
         Vector_t entry3D = toBegin.getOrigin();
         Vector_t exit3D = toEnd.getOrigin();
 
@@ -388,7 +388,7 @@ void OpalBeamline::save3DLattice() {
 
             Bend2D * bendElement = static_cast<Bend2D*>(element.get());
             std::vector<Vector_t> designPath = bendElement->getDesignPath();
-            toEnd = bendElement->getBeginToEnd_local() * (*it).getCoordTransformationTo();
+            toEnd = bendElement->getBeginToEnd_local() * element->getCSTrafoGlobal2Local();
             exit3D = toEnd.getOrigin();
 
             unsigned int size = designPath.size();
@@ -402,7 +402,7 @@ void OpalBeamline::save3DLattice() {
                 << std::setw(18) << std::setprecision(10) << entry3D(1)
                 << "\n";
 
-            Vector_t position = (*it).getCoordTransformationTo().transformFrom(designPath.front());
+            Vector_t position = element->getCSTrafoGlobal2Local().transformFrom(designPath.front());
             pos << std::setw(30) << std::left << std::string("\"BEGIN: ") + element->getName() + std::string("\"")
                 << std::setw(18) << std::setprecision(10) << position(2)
                 << std::setw(18) << std::setprecision(10) << position(0)
@@ -411,7 +411,7 @@ void OpalBeamline::save3DLattice() {
 
             for (unsigned int i = frequency; i + 1 < size; i += frequency) {
 
-                position = (*it).getCoordTransformationTo().transformFrom(designPath[i]);
+                position = element->getCSTrafoGlobal2Local().transformFrom(designPath[i]);
                 pos << std::setw(30) << std::left << std::string("\"MID: ") + element->getName() + std::string("\"")
                     << std::setw(18) << std::setprecision(10) << position(2)
                     << std::setw(18) << std::setprecision(10) << position(0)
@@ -419,7 +419,7 @@ void OpalBeamline::save3DLattice() {
                     << std::endl;
             }
 
-            position = (*it).getCoordTransformationTo().transformFrom(designPath.back());
+            position = element->getCSTrafoGlobal2Local().transformFrom(designPath.back());
             pos << std::setw(30) << std::left << std::string("\"END: ") + element->getName() + std::string("\"")
                 << std::setw(18) << std::setprecision(10) << position(2)
                 << std::setw(18) << std::setprecision(10) << position(0)
@@ -534,15 +534,16 @@ void OpalBeamline::save3DInput() {
     std::ofstream pos(fname);
 
     for (; it != end; ++ it) {
-        std::string element = (*it).getElement()->getName();
-        const boost::regex replacePSI("(" + element + "\\s*:[^\\n]*)PSI\\s*=[^,;]*,?", boost::regex::icase);
+        std::shared_ptr<Component> element = (*it).getElement();
+        std::string elementName = element->getName();
+        const boost::regex replacePSI("(" + elementName + "\\s*:[^\\n]*)PSI\\s*=[^,;]*,?", boost::regex::icase);
         input = boost::regex_replace(input, replacePSI, "\\1\\2");
 
-        const boost::regex replaceELEMEDGE("(" + element + "\\s*:[^\\n]*)ELEMEDGE\\s*=[^,;]*(.)", boost::regex::icase);
+        const boost::regex replaceELEMEDGE("(" + elementName + "\\s*:[^\\n]*)ELEMEDGE\\s*=[^,;]*(.)", boost::regex::icase);
 
-        CoordinateSystemTrafo cst = (*it).getCoordTransformationTo();
+        CoordinateSystemTrafo cst = element->getCSTrafoGlobal2Local();
         Vector_t origin = cst.getOrigin();
-        Vector_t orient = Util::getTaitBryantAngles(cst.getRotation().conjugate(), element);
+        Vector_t orient = Util::getTaitBryantAngles(cst.getRotation().conjugate(), elementName);
         for (unsigned int d = 0; d < 3; ++ d)
             orient(d) *= Physics::rad2deg;
 
@@ -562,20 +563,20 @@ void OpalBeamline::save3DInput() {
 
         input = boost::regex_replace(input, replaceELEMEDGE, position);
 
-        if ((*it).getElement()->getType() == ElementBase::RBEND ||
-            (*it).getElement()->getType() == ElementBase::SBEND) {
-            const Bend2D* dipole = static_cast<const Bend2D*>((*it).getElement().get());
+        if (element->getType() == ElementBase::RBEND ||
+            element->getType() == ElementBase::SBEND) {
+            const Bend2D* dipole = static_cast<const Bend2D*>(element.get());
             double angle = dipole->getBendAngle();
             double E1 = dipole->getEntranceAngle();
             double E2 = dipole->getExitAngle();
 
-            const boost::regex angleR("(" + element + "\\s*:[^\\n]*ANGLE\\s*=)[^,;]*(.)");
+            const boost::regex angleR("(" + elementName + "\\s*:[^\\n]*ANGLE\\s*=)[^,;]*(.)");
             const std::string angleF("\\1 " + round2string(angle * 180 / Physics::pi, 6) + " / 180 * PI\\2");
-            const boost::regex E1R("(" + element + "\\s*:[^\\n]*E1\\s*=)[^,;]*(.)");
+            const boost::regex E1R("(" + elementName + "\\s*:[^\\n]*E1\\s*=)[^,;]*(.)");
             const std::string E1F("\\1 " + round2string(E1 * 180 / Physics::pi, 6) + " / 180 * PI\\2");
-            const boost::regex E2R("(" + element + "\\s*:[^\\n]*E2\\s*=)[^,;]*(.)");
+            const boost::regex E2R("(" + elementName + "\\s*:[^\\n]*E2\\s*=)[^,;]*(.)");
             const std::string E2F("\\1 " + round2string(E2 * 180 / Physics::pi, 6) + " / 180 * PI\\2");
-            const boost::regex noRotation("(" + element + "\\s*:[^\\n]*),\\s*ROTATION\\s*=[^,;]*(.)");
+            const boost::regex noRotation("(" + elementName + "\\s*:[^\\n]*),\\s*ROTATION\\s*=[^,;]*(.)");
             const std::string noRotationFormat("\\1\\2  ");
 
             input = boost::regex_replace(input, angleR, angleF);