Verified Commit 9672a572 authored by Guilherme Amadio's avatar Guilherme Amadio
Browse files

Replace private dx/dy/dz with public x()/y()/z() in Hep3Vector

This is in preparation to change the data layout of Hep3Vector
by replacing the private dx/dy/dz members with a vector, data[3],
to allow us to remove the inefficient switch from operator[]().
parent d0fa71ed
......@@ -128,9 +128,7 @@ inline Hep3Vector::Hep3Vector(const Hep3Vector & p)
inline Hep3Vector::~Hep3Vector() {}
inline Hep3Vector & Hep3Vector::operator = (const Hep3Vector & p) {
dx = p.dx;
dy = p.dy;
dz = p.dz;
set(p.x(), p.y(), p.z());
return *this;
}
......@@ -140,15 +138,15 @@ inline Hep3Vector & Hep3Vector::operator = (const Hep3Vector & p) {
// r, theta, phi
inline double Hep3Vector::mag2() const { return dx*dx + dy*dy + dz*dz; }
inline double Hep3Vector::mag2() const { return x()*x() + y()*y() + z()*z(); }
inline double Hep3Vector::mag() const { return std::sqrt(mag2()); }
inline double Hep3Vector::r() const { return mag(); }
inline double Hep3Vector::theta() const {
return dx == 0.0 && dy == 0.0 && dz == 0.0 ? 0.0 : std::atan2(perp(),dz);
return x() == 0.0 && y() == 0.0 && z() == 0.0 ? 0.0 : std::atan2(perp(),z());
}
inline double Hep3Vector::phi() const {
return dx == 0.0 && dy == 0.0 ? 0.0 : std::atan2(dy,dx);
return x() == 0.0 && y() == 0.0 ? 0.0 : std::atan2(y(),x());
}
inline double Hep3Vector::getR() const { return mag(); }
......@@ -158,12 +156,12 @@ inline double Hep3Vector::angle() const { return theta(); }
inline double Hep3Vector::cosTheta() const {
double ptot = mag();
return ptot == 0.0 ? 1.0 : dz/ptot;
return ptot == 0.0 ? 1.0 : z()/ptot;
}
inline double Hep3Vector::cos2Theta() const {
double ptot2 = mag2();
return ptot2 == 0.0 ? 1.0 : dz*dz/ptot2;
return ptot2 == 0.0 ? 1.0 : z()*z()/ptot2;
}
inline void Hep3Vector::setR(double r1) { setMag(r1); }
......@@ -184,7 +182,7 @@ inline void Hep3Vector::setPhi(double ph) {
// perp, eta,
inline double Hep3Vector::perp2() const { return dx*dx + dy*dy; }
inline double Hep3Vector::perp2() const { return x()*x() + y()*y(); }
inline double Hep3Vector::perp() const { return std::sqrt(perp2()); }
inline double Hep3Vector::rho() const { return perp(); }
inline double Hep3Vector::eta() const { return pseudoRapidity();}
......@@ -236,7 +234,7 @@ inline Hep3Vector& Hep3Vector::operator -= (const Hep3Vector & p) {
}
inline Hep3Vector Hep3Vector::operator - () const {
return Hep3Vector(-dx, -dy, -dz);
return Hep3Vector(-x(), -y(), -z());
}
inline Hep3Vector& Hep3Vector::operator *= (double a) {
......@@ -255,11 +253,11 @@ inline double Hep3Vector::diff2(const Hep3Vector & p) const {
}
inline double Hep3Vector::dot(const Hep3Vector & p) const {
return dx*p.x() + dy*p.y() + dz*p.z();
return x()*p.x() + y()*p.y() + z()*p.z();
}
inline Hep3Vector Hep3Vector::cross(const Hep3Vector & p) const {
return Hep3Vector(dy*p.z()-p.y()*dz, dz*p.x()-p.z()*dx, dx*p.y()-p.x()*dy);
return Hep3Vector(y()*p.z()-p.y()*z(), z()*p.x()-p.z()*x(), x()*p.y()-p.x()*y());
}
inline double Hep3Vector::perp2(const Hep3Vector & p) const {
......@@ -273,10 +271,10 @@ inline double Hep3Vector::perp(const Hep3Vector & p) const {
}
inline Hep3Vector Hep3Vector::perpPart () const {
return Hep3Vector (dx, dy, 0);
return Hep3Vector (x(), y(), 0);
}
inline Hep3Vector Hep3Vector::project () const {
return Hep3Vector (0, 0, dz);
return Hep3Vector (0, 0, z());
}
inline Hep3Vector Hep3Vector::perpPart (const Hep3Vector & v2) const {
......@@ -306,13 +304,13 @@ inline Hep3Vector Hep3Vector::unit() const {
}
inline Hep3Vector Hep3Vector::orthogonal() const {
double xx = dx < 0.0 ? -dx : dx;
double yy = dy < 0.0 ? -dy : dy;
double zz = dz < 0.0 ? -dz : dz;
double xx = x() < 0.0 ? -x() : x();
double yy = y() < 0.0 ? -y() : y();
double zz = z() < 0.0 ? -z() : z();
if (xx < yy) {
return xx < zz ? Hep3Vector(0,dz,-dy) : Hep3Vector(dy,-dx,0);
return xx < zz ? Hep3Vector(0,z(),-y()) : Hep3Vector(y(),-x(),0);
}else{
return yy < zz ? Hep3Vector(-dz,0,dx) : Hep3Vector(dy,-dx,0);
return yy < zz ? Hep3Vector(-z(),0,x()) : Hep3Vector(y(),-x(),0);
}
}
......
......@@ -44,10 +44,10 @@ void Hep3Vector::setSpherical (
"Spherical coordinates set with theta not in [0, PI]"));
// No special return needed if warning is ignored.
}
dz = r1 * std::cos(theta1);
double rho1 ( r1*std::sin(theta1));
dy = rho1 * std::sin (phi1);
dx = rho1 * std::cos (phi1);
setZ(r1 * std::cos(theta1));
setY(rho1 * std::sin (phi1));
setX(rho1 * std::cos (phi1));
return;
} /* setSpherical (r, theta1, phi) */
......@@ -60,9 +60,9 @@ void Hep3Vector::setCylindrical (
"Cylindrical coordinates supplied with negative Rho"));
// No special return needed if warning is ignored.
}
dz = z1;
dy = rho1 * std::sin (phi1);
dx = rho1 * std::cos (phi1);
setZ(z1);
setY(rho1 * std::sin (phi1));
setX(rho1 * std::cos (phi1));
return;
} /* setCylindrical (r, phi, z) */
......@@ -74,7 +74,7 @@ void Hep3Vector::setRhoPhiTheta (
ZMthrowC (ZMxpvZeroVector(
"Attempt set vector components rho, phi, theta with zero rho -- "
"zero vector is returned, ignoring theta and phi"));
dx = 0; dy = 0; dz = 0;
set(0.0, 0.0, 0.0);
return;
}
if ( (theta1 == 0) || (theta1 == CLHEP::pi) ) {
......@@ -87,9 +87,9 @@ void Hep3Vector::setRhoPhiTheta (
"Rho, phi, theta set with theta not in [0, PI]"));
// No special return needed if warning is ignored.
}
dz = rho1 / std::tan (theta1);
dy = rho1 * std::sin (phi1);
dx = rho1 * std::cos (phi1);
setZ(rho1 / std::tan (theta1));
setY(rho1 * std::sin (phi1));
setX(rho1 * std::cos (phi1));
return;
} /* setCyl (rho, phi, theta) */
......@@ -101,13 +101,13 @@ void Hep3Vector::setRhoPhiEta (
ZMthrowC (ZMxpvZeroVector(
"Attempt set vector components rho, phi, eta with zero rho -- "
"zero vector is returned, ignoring eta and phi"));
dx = 0; dy = 0; dz = 0;
set(0.0, 0.0, 0.0);
return;
}
double theta1 (2 * std::atan ( std::exp (-eta1) ));
dz = rho1 / std::tan (theta1);
dy = rho1 * std::sin (phi1);
dx = rho1 * std::cos (phi1);
setZ(rho1 / std::tan (theta1));
setY(rho1 * std::sin (phi1));
setX(rho1 * std::cos (phi1));
return;
} /* setCyl (rho, phi, eta) */
......@@ -119,17 +119,17 @@ void Hep3Vector::setRhoPhiEta (
//************
int Hep3Vector::compare (const Hep3Vector & v) const {
if ( dz > v.dz ) {
if ( z() > v.z() ) {
return 1;
} else if ( dz < v.dz ) {
} else if ( z() < v.z() ) {
return -1;
} else if ( dy > v.dy ) {
} else if ( y() > v.y() ) {
return 1;
} else if ( dy < v.dy ) {
} else if ( y() < v.y() ) {
return -1;
} else if ( dx > v.dx ) {
} else if ( x() > v.x() ) {
return 1;
} else if ( dx < v.dx ) {
} else if ( x() < v.x() ) {
return -1;
} else {
return 0;
......@@ -205,9 +205,9 @@ bool Hep3Vector::isParallel (const Hep3Vector & v,
// At this point we know v1v2 can be squared.
Hep3Vector v1Xv2 ( cross(v) );
if ( (std::fabs (v1Xv2.dx) > TOOBIG) ||
(std::fabs (v1Xv2.dy) > TOOBIG) ||
(std::fabs (v1Xv2.dz) > TOOBIG) ) {
if ( (std::fabs (v1Xv2.x()) > TOOBIG) ||
(std::fabs (v1Xv2.y()) > TOOBIG) ||
(std::fabs (v1Xv2.z()) > TOOBIG) ) {
return false;
}
......@@ -256,9 +256,9 @@ bool Hep3Vector::isOrthogonal (const Hep3Vector & v,
// At this point we know v1v2 can be squared.
Hep3Vector eps_v1Xv2 ( cross(epsilon*v) );
if ( (std::fabs (eps_v1Xv2.dx) > TOOBIG) ||
(std::fabs (eps_v1Xv2.dy) > TOOBIG) ||
(std::fabs (eps_v1Xv2.dz) > TOOBIG) ) {
if ( (std::fabs (eps_v1Xv2.x()) > TOOBIG) ||
(std::fabs (eps_v1Xv2.y()) > TOOBIG) ||
(std::fabs (eps_v1Xv2.z()) > TOOBIG) ) {
return true;
}
......
......@@ -49,18 +49,18 @@ double Hep3Vector::gamma() const {
}
double Hep3Vector::rapidity() const {
if (std::fabs(dz) == 1) {
if (std::fabs(z()) == 1) {
ZMthrowC (ZMxpvTachyonic(
"Rapidity in Z direction taken for Hep3Vector with |Z| = 1 -- \n"
"the log should return infinity"));
}
if (std::fabs(dz) > 1) {
if (std::fabs(z()) > 1) {
ZMthrowA (ZMxpvTachyonic(
"Rapidity in Z direction taken for Hep3Vector with |Z| > 1 -- \n"
"the log would return a NAN" ));
}
// Want inverse std::tanh(dz):
return (.5 * std::log((1+dz)/(1-dz)) );
// Want inverse std::tanh(z()):
return (.5 * std::log((1+z())/(1-z())) );
}
double Hep3Vector::coLinearRapidity() const {
......
......@@ -40,27 +40,24 @@ Hep3Vector & Hep3Vector::rotate (const Hep3Vector & axis,
double rz;
{ double ocdux = ocd * ux;
rx = dx * ( cd + ocdux * ux ) +
dy * ( ocdux * uy - sd * uz ) +
dz * ( ocdux * uz + sd * uy ) ;
rx = x() * ( cd + ocdux * ux ) +
y() * ( ocdux * uy - sd * uz ) +
z() * ( ocdux * uz + sd * uy ) ;
}
{ double ocduy = ocd * uy;
ry = dy * ( cd + ocduy * uy ) +
dz * ( ocduy * uz - sd * ux ) +
dx * ( ocduy * ux + sd * uz ) ;
ry = y() * ( cd + ocduy * uy ) +
z() * ( ocduy * uz - sd * ux ) +
x() * ( ocduy * ux + sd * uz ) ;
}
{ double ocduz = ocd * uz;
rz = dz * ( cd + ocduz * uz ) +
dx * ( ocduz * ux - sd * uy ) +
dy * ( ocduz * uy + sd * ux ) ;
rz = z() * ( cd + ocduz * uz ) +
x() * ( ocduz * ux - sd * uy ) +
y() * ( ocduz * uy + sd * ux ) ;
}
dx = rx;
dy = ry;
dz = rz;
set(rx, ry, rz);
return *this;
} /* rotate */
......@@ -82,22 +79,19 @@ Hep3Vector & Hep3Vector::rotate (double phi1,
double sinTheta = std::sin( theta1 ), cosTheta1 = std::cos( theta1 );
double sinPsi = std::sin( psi1 ), cosPsi = std::cos( psi1 );
rx = (cosPsi * cosPhi - cosTheta1 * sinPsi * sinPhi) * dx +
(cosPsi * sinPhi + cosTheta1 * sinPsi * cosPhi) * dy +
(sinPsi * sinTheta) * dz ;
ry = (- sinPsi * cosPhi - cosTheta1 * cosPsi * sinPhi) * dx +
(- sinPsi * sinPhi + cosTheta1 * cosPsi * cosPhi) * dy +
(cosPsi * sinTheta) * dz ;
rx = (cosPsi * cosPhi - cosTheta1 * sinPsi * sinPhi) * x() +
(cosPsi * sinPhi + cosTheta1 * sinPsi * cosPhi) * y() +
(sinPsi * sinTheta) * z() ;
rz = (sinTheta * sinPhi) * dx +
(- sinTheta * cosPhi) * dy +
(cosTheta1) * dz ;
ry = (- sinPsi * cosPhi - cosTheta1 * cosPsi * sinPhi) * x() +
(- sinPsi * sinPhi + cosTheta1 * cosPsi * cosPhi) * y() +
(cosPsi * sinTheta) * z() ;
dx = rx;
dy = ry;
dz = rz;
rz = (sinTheta * sinPhi) * x() +
(- sinTheta * cosPhi) * y() +
(cosTheta1) * z() ;
set(rx, ry, rz);
return *this;
} /* rotate */
......
......@@ -41,15 +41,17 @@ Hep3Vector & Hep3Vector::rotateUz(const Hep3Vector& NewUzVector) {
double u3 = NewUzVector.z();
double up = u1*u1 + u2*u2;
if (up>0) {
up = std::sqrt(up);
double px = dx, py = dy, pz = dz;
dx = (u1*u3*px - u2*py)/up + u1*pz;
dy = (u2*u3*px + u1*py)/up + u2*pz;
dz = -up*px + u3*pz;
}
else if (u3 < 0.) { dx = -dx; dz = -dz; } // phi=0 teta=pi
else {};
if (up > 0) {
up = std::sqrt(up);
double px = (u1 * u3 * x() - u2 * y()) / up + u1 * z();
double py = (u2 * u3 * x() + u1 * y()) / up + u2 * z();
double pz = -up * x() + u3 * z();
set(px, py, pz);
} else if (u3 < 0.) {
setX(-x());
setZ(-z());
} // phi=0 teta=pi
return *this;
}
......@@ -88,30 +90,30 @@ const Hep3Vector HepZHat(0.0, 0.0, 1.0);
Hep3Vector & Hep3Vector::rotateX (double phi1) {
double sinphi = std::sin(phi1);
double cosphi = std::cos(phi1);
double ty;
ty = dy * cosphi - dz * sinphi;
dz = dz * cosphi + dy * sinphi;
dy = ty;
double ty = y() * cosphi - z() * sinphi;
double tz = z() * cosphi + y() * sinphi;
setY(ty);
setZ(tz);
return *this;
} /* rotateX */
Hep3Vector & Hep3Vector::rotateY (double phi1) {
double sinphi = std::sin(phi1);
double cosphi = std::cos(phi1);
double tz;
tz = dz * cosphi - dx * sinphi;
dx = dx * cosphi + dz * sinphi;
dz = tz;
double tx = x() * cosphi + z() * sinphi;
double tz = z() * cosphi - x() * sinphi;
setX(tx);
setZ(tz);
return *this;
} /* rotateY */
Hep3Vector & Hep3Vector::rotateZ (double phi1) {
double sinphi = std::sin(phi1);
double cosphi = std::cos(phi1);
double tx;
tx = dx * cosphi - dy * sinphi;
dy = dy * cosphi + dx * sinphi;
dx = tx;
double tx = x() * cosphi - y() * sinphi;
double ty = y() * cosphi + x() * sinphi;
setX(tx);
setY(ty);
return *this;
} /* rotateZ */
......@@ -181,15 +183,15 @@ double Hep3Vector::cos2Theta(const Hep3Vector & q) const {
void Hep3Vector::setEta (double eta1) {
double phi1 = 0;
double r1;
if ( (dx == 0) && (dy == 0) ) {
if (dz == 0) {
if ( (x() == 0) && (y() == 0) ) {
if (z() == 0) {
ZMthrowC (ZMxpvZeroVector(
"Attempt to set eta of zero vector -- vector is unchanged"));
return;
}
ZMthrowC (ZMxpvZeroVector(
"Attempt to set eta of vector along Z axis -- will use phi = 0"));
r1 = std::fabs(dz);
r1 = std::fabs(z());
} else {
r1 = getR();
phi1 = getPhi();
......@@ -197,10 +199,10 @@ void Hep3Vector::setEta (double eta1) {
double tanHalfTheta = std::exp ( -eta1 );
double cosTheta1 =
(1 - tanHalfTheta*tanHalfTheta) / (1 + tanHalfTheta*tanHalfTheta);
dz = r1 * cosTheta1;
double rho1 = r1*std::sqrt(1 - cosTheta1*cosTheta1);
dy = rho1 * std::sin (phi1);
dx = rho1 * std::cos (phi1);
setZ(r1 * cosTheta1);
setY(rho1 * std::sin (phi1));
setX(rho1 * std::cos (phi1));
return;
}
......@@ -208,25 +210,25 @@ void Hep3Vector::setCylTheta (double theta1) {
// In cylindrical coords, set theta while keeping rho and phi fixed
if ( (dx == 0) && (dy == 0) ) {
if (dz == 0) {
if ( (x() == 0) && (y() == 0) ) {
if (z() == 0) {
ZMthrowC (ZMxpvZeroVector(
"Attempt to set cylTheta of zero vector -- vector is unchanged"));
return;
}
if (theta1 == 0) {
dz = std::fabs(dz);
setZ(std::fabs(z()));
return;
}
if (theta1 == CLHEP::pi) {
dz = -std::fabs(dz);
setZ(-std::fabs(z()));
return;
}
ZMthrowC (ZMxpvZeroVector(
"Attempt set cylindrical theta of vector along Z axis "
"to a non-trivial value, while keeping rho fixed -- "
"will return zero vector"));
dz = 0;
setZ(0.0);
return;
}
if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) {
......@@ -240,12 +242,12 @@ void Hep3Vector::setCylTheta (double theta1) {
ZMthrowC (ZMxpvInfiniteVector(
"Attempt to set cylindrical theta to 0 or PI "
"while keeping rho fixed -- infinite Z will be computed"));
dz = (theta1==0) ? 1.0E72 : -1.0E72;
setZ((theta1==0) ? 1.0E72 : -1.0E72);
return;
}
dz = rho1 / std::tan (theta1);
dy = rho1 * std::sin (phi1);
dx = rho1 * std::cos (phi1);
setZ(rho1 / std::tan (theta1));
setY(rho1 * std::sin (phi1));
setX(rho1 * std::cos (phi1));
} /* setCylTheta */
......@@ -260,46 +262,43 @@ void Hep3Vector::setCylEta (double eta1) {
//-| ZMthrows to say eta rather than theta. Besides, we assumedly
//-| need not check for theta of 0 or PI.
if ( (dx == 0) && (dy == 0) ) {
if (dz == 0) {
if ( (x() == 0) && (y() == 0) ) {
if (z() == 0) {
ZMthrowC (ZMxpvZeroVector(
"Attempt to set cylEta of zero vector -- vector is unchanged"));
return;
}
if (theta1 == 0) {
dz = std::fabs(dz);
setZ(std::fabs(z()));
return;
}
if (theta1 == CLHEP::pi) {
dz = -std::fabs(dz);
setZ(-std::fabs(z()));
return;
}
ZMthrowC (ZMxpvZeroVector(
"Attempt set cylindrical eta of vector along Z axis "
"to a non-trivial value, while keeping rho fixed -- "
"will return zero vector"));
dz = 0;
setZ(0.0);
return;
}
double phi1 (getPhi());
double rho1 = getRho();
dz = rho1 / std::tan (theta1);
dy = rho1 * std::sin (phi1);
dx = rho1 * std::cos (phi1);
setZ(rho1 / std::tan (theta1));
setY(rho1 * std::sin (phi1));
setX(rho1 * std::cos (phi1));
} /* setCylEta */
Hep3Vector operator/ ( const Hep3Vector & v1, double c ) {
Hep3Vector operator/ ( const Hep3Vector & v1, double c ) {
if (c == 0) {
ZMthrowA ( ZMxpvInfiniteVector (
"Attempt to divide vector by 0 -- "
"will produce infinities and/or NANs"));
}
double oneOverC = 1.0/c;
return Hep3Vector ( v1.x() * oneOverC,
v1.y() * oneOverC,
v1.z() * oneOverC );
return v1 * (1.0/c);
} /* v / c */
Hep3Vector & Hep3Vector::operator/= (double c) {
......@@ -308,10 +307,7 @@ Hep3Vector & Hep3Vector::operator/= (double c) {
"Attempt to do vector /= 0 -- "
"division by zero would produce infinite or NAN components"));
}
double oneOverC = 1.0/c;
dx *= oneOverC;
dy *= oneOverC;
dz *= oneOverC;
*this *= 1.0/c;
return *this;
}
......
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