Commit 3d0b60e0 authored by kraus's avatar kraus Committed by 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 993ba266
......@@ -75,15 +75,17 @@ bool Monitor::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B)
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);
}
}
......@@ -97,15 +99,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,
......@@ -116,8 +118,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],
......
......@@ -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
......@@ -25,6 +25,11 @@ namespace Util {
return sqrt(dot(p, p) + 1.0);
}
inline
Vector_t getBeta(Vector_t p) {
return p / getGamma(p);
}
inline
double getEnergy(Vector_t p, double mass) {
return (getGamma(p) - 1.0) * mass;
......@@ -36,6 +41,11 @@ namespace Util {
return sqrt(std::pow(gamma, 2.0) - 1.0);
}
inline
double convertMomentumeVToBetaGamma(double p, double mass) {
return p / mass;
}
inline
std::string getTimeString(double time, unsigned int precision = 3) {
std::string timeUnit(" [ps]");
......@@ -215,4 +225,4 @@ std::string Util::toStringWithThousandSep(T value, char sep) {
return ret.str();
}
#endif
\ No newline at end of file
#endif
......@@ -288,9 +288,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