2 // ********************************************************************
3 // * License and Disclaimer *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
27 // $Id: G4VTwistSurface.icc 66356 2012-12-18 09:02:32Z gcosmo $
30 // --------------------------------------------------------------------
31 // G4VTwistSurface class inline methods
34 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
37 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
38 // from original version in Jupiter-2.5.02 application.
39 // --------------------------------------------------------------------
41 //=====================================================================
42 //* DistanceToPlaneWithV ----------------------------------------------
45 G4double G4VTwistSurface::DistanceToPlaneWithV(const G4ThreeVector &p,
46 const G4ThreeVector &v,
47 const G4ThreeVector &x0,
48 const G4ThreeVector &n0,
52 G4double t = kInfinity;
53 if (q) { t = (n0 * (x0 - p)) / q; }
58 //=====================================================================
59 //* DistanceToPlane ---------------------------------------------------
62 G4double G4VTwistSurface::DistanceToPlane(const G4ThreeVector &p,
63 const G4ThreeVector &x0,
64 const G4ThreeVector &n0,
68 // Calculate distance to plane in local coordinate,
69 // then return distance and global intersection points.
71 // p - location of flying particle
72 // x0 - reference point of surface
73 // xx - a foot of perpendicular line from p to the plane
74 // t - distance from xx to p
75 // n - a unit normal of this plane from plane to p.
84 // t = n * (p - x0) / std::fabs(n)
87 G4ThreeVector n = n0.unit();
93 //=====================================================================
94 //* DistanceToPlane ---------------------------------------------------
97 G4double G4VTwistSurface::DistanceToPlane(const G4ThreeVector &p,
98 const G4ThreeVector &x0,
99 const G4ThreeVector &t1,
100 const G4ThreeVector &t2,
105 // Calculate distance to plane in local coordinate,
106 // then return distance and global intersection points.
107 // t1 - 1st. vector lying on the plane
108 // t2 - 2nd. vector lying on the plane
110 n = (t1.cross(t2)).unit();
111 return DistanceToPlane(p, x0, n, xx);
114 //=====================================================================
115 //* DistanceToLine ----------------------------------------------------
118 G4double G4VTwistSurface::DistanceToLine(const G4ThreeVector &p,
119 const G4ThreeVector &x0,
120 const G4ThreeVector &d,
124 // Calculate distance to line,
125 // then return distance and global intersection points.
127 // p - location of flying particle
128 // x0 - reference point of line
129 // d - direction vector of line
130 // xx - a foot of perpendicular line from p to the plane
131 // t - distance from xx to p
135 // distance^2 = |(xx - p)|^2
139 // (d/dt)distance^2 = (d/dt)|((x0 + t*d) - p)|^2
140 // = 2*t*|d|^2 + 2*d*(x0 - p)
141 // = 0 // smallest distance
143 // t = - d*(x0 - p) / |d|^2
147 G4ThreeVector dir = d.unit();
148 t = - dir * (x0 - p); // |dir|^2 = 1.
151 G4ThreeVector dist = xx - p;
155 //=====================================================================
156 //* IsAxis0 -----------------------------------------------------------
159 G4bool G4VTwistSurface::IsAxis0(G4int areacode) const
161 if (areacode & sAxis0) return true;
165 //=====================================================================
166 //* IsAxis1 -----------------------------------------------------------
169 G4bool G4VTwistSurface::IsAxis1(G4int areacode) const
171 if (areacode & sAxis1) return true;
175 //=====================================================================
176 //* IsOutside ---------------------------------------------------------
179 G4bool G4VTwistSurface::IsOutside(G4int areacode) const
181 if (areacode & sInside) return false;
185 //=====================================================================
186 //* IsInside ----------------------------------------------------------
189 G4bool G4VTwistSurface::IsInside(G4int areacode, G4bool testbitmode) const
191 if (areacode & sInside) {
195 if (!((areacode & sBoundary) || (areacode & sCorner))) return true;
201 //=====================================================================
202 //* IsBoundary --------------------------------------------------------
205 G4bool G4VTwistSurface::IsBoundary(G4int areacode, G4bool testbitmode) const
207 if ((areacode & sBoundary) == sBoundary) {
211 if ((areacode & sInside) == sInside) return true;
217 //=====================================================================
218 //* IsCorner ----------------------------------------------------------
221 G4bool G4VTwistSurface::IsCorner(G4int areacode, G4bool testbitmode) const
223 if ((areacode & sCorner) == sCorner) {
227 if ((areacode & sInside) == sInside) return true;
233 //=====================================================================
234 //* GetAxisType -------------------------------------------------------
237 G4int G4VTwistSurface::GetAxisType(G4int areacode, G4int whichaxis) const
239 G4int axiscode = areacode & sAxisMask & whichaxis;
241 if (axiscode == (sAxisX & sAxis0) ||
242 axiscode == (sAxisX & sAxis1)) {
244 } else if (axiscode == (sAxisY & sAxis0) ||
245 axiscode == (sAxisY & sAxis1)) {
247 } else if (axiscode == (sAxisZ & sAxis0) ||
248 axiscode == (sAxisZ & sAxis1)) {
250 } else if (axiscode == (sAxisRho & sAxis0) ||
251 axiscode == (sAxisRho & sAxis1)) {
253 } else if (axiscode == (sAxisPhi & sAxis0) ||
254 axiscode == (sAxisPhi & sAxis1)) {
257 std::ostringstream message;
258 message << "Configuration not supported." << G4endl
259 << " areacode = " << areacode;
260 G4Exception("G4VTwistSurface::GetAxisType()","GeomSolids0001",
261 FatalException, message);
266 //=====================================================================
267 //* ComputeGlobalPoint ------------------------------------------------
270 G4ThreeVector G4VTwistSurface::ComputeGlobalPoint(const G4ThreeVector &lp) const
272 return fRot * G4ThreeVector(lp) + fTrans;
275 //=====================================================================
276 //* ComputeGlobalPoint ------------------------------------------------
279 G4ThreeVector G4VTwistSurface::ComputeLocalPoint(const G4ThreeVector &gp) const
281 return fRot.inverse() * ( G4ThreeVector(gp) - fTrans ) ;
284 //=====================================================================
285 //* ComputeGlobalDirection --------------------------------------------
288 G4ThreeVector G4VTwistSurface::ComputeGlobalDirection(const G4ThreeVector &lp) const
290 return fRot * G4ThreeVector(lp);
293 //=====================================================================
294 //* ComputeLocalDirection ---------------------------------------------
297 G4ThreeVector G4VTwistSurface::ComputeLocalDirection(const G4ThreeVector &gp) const
299 return fRot.inverse() * G4ThreeVector(gp);
302 //=====================================================================
303 //* SetNeighbours -----------------------------------------------------
306 void G4VTwistSurface::SetNeighbours(G4VTwistSurface* axis0min, G4VTwistSurface* axis1min,
307 G4VTwistSurface* axis0max, G4VTwistSurface* axis1max)
309 fNeighbours[0] = axis0min;
310 fNeighbours[1] = axis1min;
311 fNeighbours[2] = axis0max;
312 fNeighbours[3] = axis1max;
315 //=====================================================================
316 //* GetNeighbours -----------------------------------------------------
319 G4int G4VTwistSurface::GetNeighbours(G4int areacode, G4VTwistSurface** surfaces)
322 int sAxis0Min = sAxis0 & sAxisMin ;
323 int sAxis1Min = sAxis1 & sAxisMin ;
324 int sAxis0Max = sAxis0 & sAxisMax ;
325 int sAxis1Max = sAxis1 & sAxisMax ;
329 if ( (areacode & sAxis0Min ) == sAxis0Min ) {
330 surfaces[i] = fNeighbours[0] ;
334 if ( ( areacode & sAxis1Min ) == sAxis1Min ) {
335 surfaces[i] = fNeighbours[1] ;
337 if ( i == 2 ) return i ;
340 if ( ( areacode & sAxis0Max ) == sAxis0Max ) {
341 surfaces[i] = fNeighbours[2] ;
343 if ( i == 2 ) return i ;
346 if ( ( areacode & sAxis1Max ) == sAxis1Max ) {
347 surfaces[i] = fNeighbours[3] ;
349 if ( i == 2 ) return i ;
355 //=====================================================================
356 //* GetCorner ---------------------------------------------------------
359 G4ThreeVector G4VTwistSurface::GetCorner(G4int areacode) const
361 if (!(areacode & sCorner)){
362 std::ostringstream message;
363 message << "Area code must represent corner." << G4endl
364 << " areacode = " << areacode;
365 G4Exception("G4VTwistSurface::GetCorner()","GeomSolids0002",
366 FatalException, message);
369 if ((areacode & sC0Min1Min) == sC0Min1Min) {
371 } else if ((areacode & sC0Max1Min) == sC0Max1Min) {
373 } else if ((areacode & sC0Max1Max) == sC0Max1Max) {
375 } else if ((areacode & sC0Min1Max) == sC0Min1Max) {
378 std::ostringstream message;
379 message << "Configuration not supported." << G4endl
380 << " areacode = " << areacode;
381 G4Exception("G4VTwistSurface::GetCorner()", "GeomSolids0001",
382 FatalException, message);