13 inline const Vector2D DiscSurface::localPolarToCartesian(
16 lpolar[
eLOC_R] * sin(lpolar[eLOC_PHI]));
19 inline const Vector2D DiscSurface::localCartesianToPolar(
23 atan2(lcart[eLOC_Y], lcart[
eLOC_X]));
37 const double x = direction(0);
38 const double y = direction(1);
39 const double z = direction(2);
42 const double cos_theta =
z;
43 const double sin_theta = sqrt(x * x + y * y);
44 const double inv_sin_theta = 1. / sin_theta;
45 const double cos_phi = x * inv_sin_theta;
46 const double sin_phi = y * inv_sin_theta;
48 const auto rframe = referenceFrame(gctx, position, direction);
51 double lrad = pars[
eLOC_0];
52 double lphi = pars[
eLOC_1];
53 double lcos_phi = cos(lphi);
54 double lsin_phi = sin(lphi);
56 jacobian.block<3, 1>(0,
eLOC_0) =
57 lcos_phi * rframe.block<3, 1>(0, 0) + lsin_phi * rframe.block<3, 1>(0, 1);
58 jacobian.block<3, 1>(0,
eLOC_1) =
59 lrad * (lcos_phi * rframe.block<3, 1>(0, 1) -
60 lsin_phi * rframe.block<3, 1>(0, 0));
64 jacobian(4,
ePHI) = (-sin_theta) * sin_phi;
65 jacobian(4,
eTHETA) = cos_theta * cos_phi;
66 jacobian(5,
ePHI) = sin_theta * cos_phi;
67 jacobian(5,
eTHETA) = cos_theta * sin_phi;
68 jacobian(6,
eTHETA) = (-sin_theta);
69 jacobian(7,
eQOP) = 1;
78 const double x = direction(0);
79 const double y = direction(1);
80 const double z = direction(2);
82 const double cosTheta =
z;
83 const double sinTheta = sqrt(x * x + y * y);
84 const double invSinTheta = 1. / sinTheta;
85 const double cosPhi = x * invSinTheta;
86 const double sinPhi = y * invSinTheta;
89 referenceFrame(gctx, position, direction).transpose();
92 const double lr =
perp(pos_loc);
93 const double lphi =
phi(pos_loc);
94 const double lcphi = cos(lphi);
95 const double lsphi = sin(lphi);
97 auto lx = rframeT.block<1, 3>(0, 0);
98 auto ly = rframeT.block<1, 3>(1, 0);
99 jacobian.block<1, 3>(0, 0) = lcphi * lx + lsphi * ly;
100 jacobian.block<1, 3>(1, 0) = (lcphi * ly - lsphi * lx) / lr;
104 jacobian(
ePHI, 4) = -sinPhi * invSinTheta;
105 jacobian(
ePHI, 5) = cosPhi * invSinTheta;
106 jacobian(
eTHETA, 4) = cosPhi * cosTheta;
107 jacobian(
eTHETA, 5) = sinPhi * cosTheta;
108 jacobian(
eTHETA, 6) = -sinTheta;
109 jacobian(
eQOP, 7) = 1;
116 const Vector3D& direction,
const BoundaryCheck& bcheck)
const {
123 if (intersection.status != Intersection::Status::unreachable and bcheck and
124 m_bounds !=
nullptr) {
126 const auto& tMatrix = gctxTransform.matrix();
127 const Vector3D vecLocal(intersection.position - tMatrix.block<3, 1>(0, 3));
129 tMatrix.block<3, 2>(0, 0).transpose() * vecLocal;
130 if (bcheck.type() == BoundaryCheck::Type::eAbsolute and
131 m_bounds->coversFullAzimuth()) {
135 intersection.status = Intersection::Status::missed;
137 }
else if (not insideBounds(localCartesianToPolar(lcartesian), bcheck)) {
138 intersection.status = Intersection::Status::missed;
147 const auto& tMatrix =
transform(gctx).matrix();
148 return Vector3D(tMatrix(0, 2), tMatrix(1, 2), tMatrix(2, 2));
153 if (bValue ==
binR) {
154 double r = m_bounds->binningValueR();
155 double phi = m_bounds->binningValuePhi();
156 return Vector3D(r * cos(phi), r * sin(phi), center(gctx).
z());
161 inline double DiscSurface::binningPositionValue(
const GeometryContext& gctx,
164 if (bValue ==
binR) {
167 return GeometryObject::binningPositionValue(gctx, bValue);