diff --git a/src/Classic/AbsBeamline/Degrader.cpp b/src/Classic/AbsBeamline/Degrader.cpp index f9f8bcc20e357d1174d0d54be1d989dafcdb35d8..13930d0bcb082581531d3ec566fd941b9d922cb8 100644 --- a/src/Classic/AbsBeamline/Degrader.cpp +++ b/src/Classic/AbsBeamline/Degrader.cpp @@ -51,23 +51,14 @@ 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) - goOffline(); + if(online_m) + goOffline(); } @@ -128,13 +119,8 @@ bool Degrader::applyToReferenceParticle(const Vector_t &R, } void Degrader::initialise(PartBunchBase *bunch, double &startField, double &endField) { - RefPartBunch_m = bunch; + initialise(bunch); endField = startField + getElementLength(); - - if (filename_m == std::string("")) - lossDs_m = std::unique_ptr(new LossDataSink(getName(), !Options::asciidump)); - else - lossDs_m = std::unique_ptr(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump)); } void Degrader::initialise(PartBunchBase *bunch) { @@ -148,20 +134,22 @@ void Degrader::initialise(PartBunchBase *bunch) { void Degrader::finalise() { - *gmsg << "* Finalize Degrader" << endl; + *gmsg << "* Finalize Degrader" << endl; } 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; } diff --git a/src/Classic/Algorithms/Tracker.cpp b/src/Classic/Algorithms/Tracker.cpp index 2e90c97a521fdd317fe0ce80e320f671864a323d..76c0623f09f3a5dbb3060716d78f7d8c1557df69 100644 --- a/src/Classic/Algorithms/Tracker.cpp +++ b/src/Classic/Algorithms/Tracker.cpp @@ -81,10 +81,8 @@ typedef FTps 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, diff --git a/src/Classic/Solvers/BeamStrippingPhysics.cpp b/src/Classic/Solvers/BeamStrippingPhysics.cpp index 823207a13b287bf84fafc4285595cf4887045ed6..bb4cafcb6f11d84cd35af762d3bb7ad8d92e0d76 100644 --- a/src/Classic/Solvers/BeamStrippingPhysics.cpp +++ b/src/Classic/Solvers/BeamStrippingPhysics.cpp @@ -635,12 +635,6 @@ bool BeamStrippingPhysics::stillActive() { return locPartsInMat_m != 0; } -bool BeamStrippingPhysics::stillAlive(PartBunchBase */*bunch*/) { - bool beamstrippingAlive = true; - return beamstrippingAlive; -} - - /* Cross sections parameters for interaction with air -- [1] -> Nitrogen diff --git a/src/Classic/Solvers/BeamStrippingPhysics.hh b/src/Classic/Solvers/BeamStrippingPhysics.hh index 931ea0b3109e8909d14ac75e9354020d450b7196..8aa88ae3b72cc416dd63bb961fbaa830c633822a 100644 --- a/src/Classic/Solvers/BeamStrippingPhysics.hh +++ b/src/Classic/Solvers/BeamStrippingPhysics.hh @@ -52,15 +52,14 @@ public: const std::pair &boundingSphere); virtual const std::string getType() const; - void print(Inform& msg); - bool stillActive(); - bool stillAlive(PartBunchBase *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 *bunch); private: diff --git a/src/Classic/Solvers/CollimatorPhysics.cpp b/src/Classic/Solvers/CollimatorPhysics.cpp index 8c5812a5e718a32dfff77b90021f38f0ea797bc4..2041184624321fa72be55361fcbc15128dfcc28b 100644 --- a/src/Classic/Solvers/CollimatorPhysics.cpp +++ b/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 *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 *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(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; diff --git a/src/Classic/Solvers/CollimatorPhysics.hh b/src/Classic/Solvers/CollimatorPhysics.hh index d0690b62877464db53375068f0a27dc6f7d938f5..bccf439b776605a1c964ecf58966581edca980cf 100644 --- a/src/Classic/Solvers/CollimatorPhysics.hh +++ b/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 *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, diff --git a/src/Classic/Solvers/ParticleMatterInteractionHandler.hh b/src/Classic/Solvers/ParticleMatterInteractionHandler.hh index f952bdb51946951158d6c5d8e122878fb27deeb7..af4fc75ad42382060a326fa9746be44156b08d85 100644 --- a/src/Classic/Solvers/ParticleMatterInteractionHandler.hh +++ b/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 *bunch) = 0; virtual double getTime() = 0; virtual std::string getName() = 0; virtual size_t getParticlesInMat() = 0; diff --git a/src/Classic/Utilities/ClassicField.h b/src/Classic/Utilities/ClassicField.h index 766780492036d840888bc34ae27bf6d9299490a3..97cc566a2c6c05bf5da4b6bb1780d4ed44ff49b9 100644 --- a/src/Classic/Utilities/ClassicField.h +++ b/src/Classic/Utilities/ClassicField.h @@ -1,12 +1,9 @@ #ifndef CLASSIC_FIELD_H #define CLASSIC_FIELD_H -#include #include #include #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 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(); diff --git a/src/Elements/OpalBeamline.cpp b/src/Elements/OpalBeamline.cpp index 3a9cddc3cda176241cbec2142107ae2e8006cc06..83637239f980e737125ca1a9c0576ac0d93df284 100644 --- a/src/Elements/OpalBeamline.cpp +++ b/src/Elements/OpalBeamline.cpp @@ -53,7 +53,7 @@ std::set> OpalBeamline::getElements(const Vector_t &x const FieldList::iterator end = elements_m.end(); for (; it != end; ++ it) { std::shared_ptr 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 element = (*it).getElement(); + if (element->isPositioned()) { continue; } (*it).order_m = minOrder; - std::shared_ptr 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 element = (*it).getElement(); + if (element->isPositioned()) continue; (*it).order_m = order ++; - std::shared_ptr 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 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(element.get()); std::vector 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 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((*it).getElement().get()); + if (element->getType() == ElementBase::RBEND || + element->getType() == ElementBase::SBEND) { + const Bend2D* dipole = static_cast(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);