diff --git a/src/Classic/AbsBeamline/Monitor.cpp b/src/Classic/AbsBeamline/Monitor.cpp index 442bf9ecaea260fd9a325c54eb05105dbf048c85..71aac8a96336007e93bd490f653396a141cd3f6a 100644 --- a/src/Classic/AbsBeamline/Monitor.cpp +++ b/src/Classic/AbsBeamline/Monitor.cpp @@ -74,15 +74,17 @@ bool Monitor::apply(const size_t &i, const double &t, Vector_t &/*E*/, Vector_t const Vector_t &R = RefPartBunch_m->R[i]; const Vector_t &P = RefPartBunch_m->P[i]; const double &dt = RefPartBunch_m->dt[i]; - const double recpgamma = Physics::c * dt / Util::getGamma(P); - const double middle = 0.5 * getElementLength(); + const Vector_t singleStep = Physics::c * dt * Util::getBeta(P); if (online_m && type_m == SPATIAL) { - if (dt * R(2) < dt * middle && // multiply with dt to allow back tracking - dt * (R(2) + P(2) * recpgamma) > dt * middle) { - double frac = (middle - R(2)) / (P(2) * recpgamma); - - lossDs_m->addParticle(R + frac * recpgamma * P, - P, RefPartBunch_m->ID[i], t + frac * dt, 0); + if (dt * R(2) < 0.0 && + dt * (R(2) + singleStep(2)) > 0.0) { + double frac = R(2) / singleStep(2); + + lossDs_m->addParticle(R + frac * singleStep, + P, + RefPartBunch_m->ID[i], + t + frac * dt, + 0); } } @@ -96,15 +98,15 @@ bool Monitor::applyToReferenceParticle(const Vector_t &R, Vector_t &) { if (!OpalData::getInstance()->isInPrepState()) { const double dt = RefPartBunch_m->getdT(); - const double recpgamma = Physics::c * dt / Util::getGamma(P); - const double middle = 0.5 * getElementLength(); + const double cdt = Physics::c * dt; + const Vector_t singleStep = cdt * Util::getBeta(P); - if (dt * R(2) < dt * middle && // multiply with dt to allow back tracking - dt * (R(2) + P(2) * recpgamma) > dt * middle) { - double frac = (middle - R(2)) / (P(2) * recpgamma); + if (dt * R(2) < 0.0 && + dt * (R(2) + singleStep(2)) > 0.0) { + double frac = -R(2) / singleStep(2); double time = t + frac * dt; - Vector_t dR = (0.5 + frac) * P * recpgamma; - double ds = euclidean_norm(dR); + Vector_t dR = frac * singleStep; + double ds = euclidean_norm(dR + 0.5 * singleStep); lossDs_m->addReferenceParticle(csTrafoGlobal2Local_m.transformFrom(R + dR), csTrafoGlobal2Local_m.rotateFrom(P), time, @@ -115,10 +117,13 @@ bool Monitor::applyToReferenceParticle(const Vector_t &R, const unsigned int localNum = RefPartBunch_m->getLocalNum(); for (unsigned int i = 0; i < localNum; ++ i) { - const double recpgamma = Physics::c * dt / Util::getGamma(RefPartBunch_m->P[i]); - lossDs_m->addParticle(RefPartBunch_m->R[i] + frac * RefPartBunch_m->P[i] * recpgamma - halfLength_s, - RefPartBunch_m->P[i], RefPartBunch_m->ID[i], - time, 0); + Vector_t shift = ((frac - 0.5) * cdt * Util::getBeta(RefPartBunch_m->P[i]) + - singleStep); + lossDs_m->addParticle(RefPartBunch_m->R[i] + shift, + RefPartBunch_m->P[i], + RefPartBunch_m->ID[i], + time, + 0); } OpalData::OPENMODE openMode; if (numPassages_m > 0) { diff --git a/src/Classic/AbsBeamline/Monitor.h b/src/Classic/AbsBeamline/Monitor.h index 16cffa89dbe8576dd668516467519139f059f2db..81231d1c6dc97c3fb5dee4c92f682c3cf2989fb8 100644 --- a/src/Classic/AbsBeamline/Monitor.h +++ b/src/Classic/AbsBeamline/Monitor.h @@ -105,6 +105,8 @@ public: void setType(Type type); static void writeStatistics(); + + virtual bool isInside(const Vector_t &r) const override; private: // Not implemented. @@ -125,4 +127,11 @@ void Monitor::setType(Monitor::Type type) { type_m = type; } -#endif // CLASSIC_Monitor_HH \ No newline at end of file +inline +bool Monitor::isInside(const Vector_t &r) const +{ + const double length = getElementLength(); + return std::abs(r(2)) <= 0.5 * length && isInsideTransverse(r); +} + +#endif // CLASSIC_Monitor_HH diff --git a/src/Classic/Utilities/Util.h b/src/Classic/Utilities/Util.h index 42479de4b024055ced8e0e5cc28a80a702a933e2..44816c5a14fb74d097b343aeeea91d2952c61930 100644 --- a/src/Classic/Utilities/Util.h +++ b/src/Classic/Utilities/Util.h @@ -28,6 +28,11 @@ namespace Util { return std::sqrt(dot(p, p) + 1.0); } + inline + Vector_t getBeta(Vector_t p) { + return p / getGamma(p); + } + inline double getKineticEnergy(Vector_t p, double mass) { return (getGamma(p) - 1.0) * mass; @@ -45,7 +50,7 @@ namespace Util { double convertMomentumeVToBetaGamma(double p, double mass) { return p / mass; } - + inline std::string getTimeString(double time, unsigned int precision = 3) { std::string timeUnit(" [ps]"); @@ -170,7 +175,7 @@ namespace Util { } Vector_t getTaitBryantAngles(Quaternion rotation, const std::string &elementName = ""); - + std::string toUpper(const std::string &str); std::string combineFilePath(std::initializer_list<std::string>); diff --git a/src/Elements/OpalBeamline.cpp b/src/Elements/OpalBeamline.cpp index 054b74f2d82ed73e26a8fdf7bfb3d0eec7626b59..d57d1cd6774f8eab9a414219204dcf18d341701a 100644 --- a/src/Elements/OpalBeamline.cpp +++ b/src/Elements/OpalBeamline.cpp @@ -322,9 +322,6 @@ void OpalBeamline::compute3DLattice() { if (element->getType() == ElementBase::SOURCE) { beginThis3D(2) -= thisLength; } - if (element->getType() == ElementBase::MONITOR) { - beginThis3D(2) -= 0.5 * thisLength; - } Vector_t endThis3D; if (element->getType() == ElementBase::SBEND ||