From 3d4f7c0f9f8069a7ba410a10dea9413d9cf5ec6b Mon Sep 17 00:00:00 2001
From: kraus <christof.kraus@psi.ch>
Date: Fri, 9 Oct 2020 14:02:03 +0200
Subject: [PATCH] fix position of particles, spos and position of reference
particle in temporal monitor; fix positioning of monitor in 3D space
---
src/Classic/AbsBeamline/Monitor.cpp | 36 +++++++++++++++--------------
src/Classic/AbsBeamline/Monitor.h | 10 ++++++++
src/Classic/Utilities/Util.h | 9 ++++++--
src/Elements/OpalBeamline.cpp | 3 ---
4 files changed, 36 insertions(+), 22 deletions(-)
diff --git a/src/Classic/AbsBeamline/Monitor.cpp b/src/Classic/AbsBeamline/Monitor.cpp
index 7d9b11f61..71aac8a96 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 = 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,8 +117,8 @@ 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]);
- Vector_t shift = frac * recpgamma * RefPartBunch_m->P[i] - Vector_t(0, 0, middle);
+ 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],
diff --git a/src/Classic/AbsBeamline/Monitor.h b/src/Classic/AbsBeamline/Monitor.h
index b44009703..54a418b2f 100644
--- a/src/Classic/AbsBeamline/Monitor.h
+++ b/src/Classic/AbsBeamline/Monitor.h
@@ -107,6 +107,8 @@ public:
static void writeStatistics();
virtual int getRequiredNumberOfTimeSteps() const override;
+
+ virtual bool isInside(const Vector_t &r) const override;
private:
// Not implemented.
@@ -133,4 +135,12 @@ int Monitor::getRequiredNumberOfTimeSteps() const
return 1;
}
+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
\ No newline at end of file
diff --git a/src/Classic/Utilities/Util.h b/src/Classic/Utilities/Util.h
index 42479de4b..44816c5a1 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 83637239f..23c27fb73 100644
--- a/src/Elements/OpalBeamline.cpp
+++ b/src/Elements/OpalBeamline.cpp
@@ -291,9 +291,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 ||
--
GitLab