Commit 52dce53c authored by snuverink_j's avatar snuverink_j
Browse files

improve calculation of point where particles crosses probe

parent 16bb99c4
......@@ -175,9 +175,8 @@ void Probe::setGeom(const double dist) {
}
bool Probe::checkProbe(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) {
bool flagprobed = false;
Vector_t rmin, rmax, probepoint;
Vector_t rmin, rmax;
bunch->get_bounds(rmin, rmax);
// interested in absolute minimum and maximum
double xmin = std::min(std::abs(rmin(0)), std::abs(rmax(0)));
......@@ -222,41 +221,46 @@ bool Probe::checkProbe(PartBunchBase<double, 3> *bunch, const int turnnumber, co
double lstep = (sqrt(1.0-1.0/(1.0+dot(meanP, meanP))) * Physics::c) * tstep*1.0e-6; // [mm]
double Swidth = lstep / sqrt( 1 + 1/stangle/stangle );
setGeom(Swidth);
Vector_t probepoint;
for(unsigned int i = 0; i < tempnum; ++i) {
pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
if(pflag != 0) {
// dist1 > 0, right hand, dt > 0; dist1 < 0, left hand, dt < 0
double dist1 = (A_m*bunch->R[i](0)+B_m*bunch->R[i](1)+C_m)/R_m/1000.0;
double k1, k2, tangle = 0.0;
if ( B_m == 0.0 ){
k1 = bunch->P[i](1)/bunch->P[i](0);
if (k1 == 0.0)
tangle = 1.0e12;
else
tangle = std::abs(1/k1);
} else if ( bunch->P[i](0) == 0.0 ) {
k2 = -A_m/B_m;
if (k2 == 0.0)
tangle = 1.0e12;
else
tangle = std::abs(1/k2);
} else {
k1 = bunch->P[i](1)/bunch->P[i](0);
k2 = -A_m/B_m;
tangle = std::abs(( k1-k2 )/(1 + k1*k2));
}
double dist2 = dist1 * sqrt( 1+1/tangle/tangle );
double dt = dist2/(sqrt(1.0-1.0/(1.0 + dot(bunch->P[i], bunch->P[i]))) * Physics::c)*1.0e9;
probepoint(0) = (B_m*B_m*bunch->R[i](0) - A_m*B_m*bunch->R[i](1)-A_m*C_m)/(R_m*R_m);
probepoint(1) = (A_m*A_m*bunch->R[i](1) - A_m*B_m*bunch->R[i](0)-B_m*C_m)/(R_m*R_m);
probepoint(2) = bunch->R[i](2);
lossDs_m->addParticle(probepoint, bunch->P[i], bunch->ID[i], t+dt, turnnumber);
flagprobed = true;
peakfinder_m->addParticle(bunch->R[i]);
if(pflag == 0) continue;
// calculate closest point at probe -> better to use momentum direction
// probe: y = -A/B * x - C/B
// perpendicular line through R: y = B/A * x + R(1) - B/A * R(0)
// probepoint(0) = (B_m*B_m*bunch->R[i](0) - A_m*B_m*bunch->R[i](1) - A_m*C_m) / (R_m*R_m);
// probepoint(1) = (A_m*A_m*bunch->R[i](1) - A_m*B_m*bunch->R[i](0) - B_m*C_m) / (R_m*R_m);
// probepoint(2) = bunch->R[i](2);
// calculate time correction for probepoint
// dist1 > 0, right hand, dt > 0; dist1 < 0, left hand, dt < 0
double dist1 = (A_m*bunch->R[i](0)+B_m*bunch->R[i](1)+C_m)/R_m/1000.0; // [m]
double k1, k2, tangle = 0.0;
if ( B_m == 0.0 ){
k1 = bunch->P[i](1)/bunch->P[i](0);
if (k1 == 0.0)
tangle = 1.0e12;
else
tangle = std::abs(1/k1);
} else if ( bunch->P[i](0) == 0.0 ) {
k2 = -A_m/B_m;
if (k2 == 0.0)
tangle = 1.0e12;
else
tangle = std::abs(1/k2);
} else {
k1 = bunch->P[i](1)/bunch->P[i](0);
k2 = -A_m/B_m;
tangle = std::abs(( k1-k2 )/(1 + k1*k2));
}
double dist2 = dist1 * sqrt( 1+1/tangle/tangle );
double dt = dist2/(sqrt(1.0-1.0/(1.0 + dot(bunch->P[i], bunch->P[i]))) * Physics::c)*1.0e9;
probepoint = bunch->R[i] + dist2 * 1000.0 * bunch->P[i] / euclidean_norm(bunch->P[i]);
lossDs_m->addParticle(probepoint, bunch->P[i], bunch->ID[i], t+dt, turnnumber);
peakfinder_m->addParticle(probepoint);
flagprobed = true;
}
}
......@@ -266,7 +270,6 @@ bool Probe::checkProbe(PartBunchBase<double, 3> *bunch, const int turnnumber, co
int Probe::checkPoint(const double &x, const double &y) {
int cn = 0;
for(int i = 0; i < 4; i++) {
if(((geom_m[i].y <= y) && (geom_m[i+1].y > y))
|| ((geom_m[i].y > y) && (geom_m[i+1].y <= y))) {
......@@ -277,7 +280,6 @@ int Probe::checkPoint(const double &x, const double &y) {
}
}
return (cn & 1); // 0 if even (out), and 1 if odd (in)
}
void Probe::getDimensions(double &zBegin, double &zEnd) const {
......
......@@ -63,7 +63,7 @@ public:
void setDimensions(double xstart, double xend, double ystart, double yend);
void setStep(double step);
///@{ Member variable access
virtual double getXstart() const;
virtual double getXend() const;
......@@ -71,7 +71,8 @@ public:
virtual double getYend() const;
virtual double getStep() const;
///@}
/// Record probe hits when bunch particles pass
bool checkProbe(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep);
virtual ElementBase::ElementType getType() const;
......@@ -91,13 +92,12 @@ private:
Point geom_m[5]; ///< actual geometry positions with adaptive width such that each particle hits probe once per turn
double step_m; ///< Step size of the probe (bin width in histogram file)
double A_m, B_m,R_m, C_m;
void setGeom(const double dist);
int checkPoint( const double & x, const double & y );
std::unique_ptr<PeakFinder> peakfinder_m;
double A_m, B_m, R_m, C_m; ///< Geometric lengths used in calculations
void setGeom(const double dist); ///< Sets geometry geom_m with probe width dist
int checkPoint( const double & x, const double & y ); ///< Checks if coordinate is within probe
std::unique_ptr<LossDataSink> lossDs_m;
std::unique_ptr<PeakFinder> peakfinder_m; ///< Pointer to Peakfinder instance
std::unique_ptr<LossDataSink> lossDs_m; ///< Pointer to Loss instance
// Not implemented.
void operator=(const Probe &);
......
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