Commit 3d4f7c0f authored by kraus's avatar kraus

fix position of particles, spos and position of reference particle in temporal...

fix position of particles, spos and position of reference particle in temporal monitor; fix positioning of monitor in 3D space
parent 131ad823
......@@ -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],
......
......@@ -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
......@@ -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>);
......
......@@ -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 ||
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment