// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // // $Id: G4TwistTubsSide.cc,v 1.6 2009/11/11 12:23:37 gcosmo Exp $ // GEANT4 tag $Name: geant4-09-03 $ // // // -------------------------------------------------------------------- // GEANT 4 class source file // // // G4TwistTubsSide.cc // // Author: // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp) // // History: // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4 // from original version in Jupiter-2.5.02 application. // 29-Apr-2004 - O.Link. Bug fixed in GetAreaCode // -------------------------------------------------------------------- #include "G4TwistTubsSide.hh" //===================================================================== //* constructors ------------------------------------------------------ G4TwistTubsSide::G4TwistTubsSide(const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, G4int handedness, const G4double kappa, const EAxis axis0, const EAxis axis1, G4double axis0min, G4double axis1min, G4double axis0max, G4double axis1max) : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1, axis0min, axis1min, axis0max, axis1max), fKappa(kappa) { if (axis0 == kZAxis && axis1 == kXAxis) { G4Exception("G4TwistTubsSide::G4TwistTubsSide()", "InvalidSetup", FatalException, "Should swap axis0 and axis1!"); } fIsValidNorm = false; SetCorners(); SetBoundaries(); } G4TwistTubsSide::G4TwistTubsSide(const G4String &name, G4double EndInnerRadius[2], G4double EndOuterRadius[2], G4double DPhi, G4double EndPhi[2], G4double EndZ[2], G4double InnerRadius, G4double OuterRadius, G4double Kappa, G4int handedness) : G4VTwistSurface(name) { fHandedness = handedness; // +z = +ve, -z = -ve fAxis[0] = kXAxis; // in local coordinate system fAxis[1] = kZAxis; fAxisMin[0] = InnerRadius; // Inner-hype radius at z=0 fAxisMax[0] = OuterRadius; // Outer-hype radius at z=0 fAxisMin[1] = EndZ[0]; fAxisMax[1] = EndZ[1]; fKappa = Kappa; fRot.rotateZ( fHandedness > 0 ? -0.5*DPhi : 0.5*DPhi ); fTrans.set(0, 0, 0); fIsValidNorm = false; SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ; SetBoundaries(); } //===================================================================== //* Fake default constructor ------------------------------------------ G4TwistTubsSide::G4TwistTubsSide( __void__& a ) : G4VTwistSurface(a) { } //===================================================================== //* destructor -------------------------------------------------------- G4TwistTubsSide::~G4TwistTubsSide() { } //===================================================================== //* GetNormal --------------------------------------------------------- G4ThreeVector G4TwistTubsSide::GetNormal(const G4ThreeVector &tmpxx, G4bool isGlobal) { // GetNormal returns a normal vector at a surface (or very close // to surface) point at tmpxx. // If isGlobal=true, it returns the normal in global coordinate. // G4ThreeVector xx; if (isGlobal) { xx = ComputeLocalPoint(tmpxx); if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) { return ComputeGlobalDirection(fCurrentNormal.normal); } } else { xx = tmpxx; if (xx == fCurrentNormal.p) { return fCurrentNormal.normal; } } G4ThreeVector er(1, fKappa * xx.z(), 0); G4ThreeVector ez(0, fKappa * xx.x(), 1); G4ThreeVector normal = fHandedness*(er.cross(ez)); if (isGlobal) { fCurrentNormal.normal = ComputeGlobalDirection(normal.unit()); } else { fCurrentNormal.normal = normal.unit(); } return fCurrentNormal.normal; } //===================================================================== //* DistanceToSurface ------------------------------------------------- G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate) { // Coordinate system: // // The coordinate system is so chosen that the intersection of // the twisted surface with the z=0 plane coincides with the // x-axis. // Rotation matrix from this coordinate system (local system) // to global system is saved in fRot field. // So the (global) particle position and (global) velocity vectors, // p and v, should be rotated fRot.inverse() in order to convert // to local vectors. // // Equation of a twisted surface: // // x(rho(z=0), z) = rho(z=0) // y(rho(z=0), z) = rho(z=0)*K*z // z(rho(z=0), z) = z // with // K = std::tan(fPhiTwist/2)/fZHalfLen // // Equation of a line: // // gxx = p + t*v // with // p = fRot.inverse()*gp // v = fRot.inverse()*gv // // Solution for intersection: // // Required time for crossing is given by solving the // following quadratic equation: // // a*t^2 + b*t + c = 0 // // where // // a = K*v_x*v_z // b = K*(v_x*p_z + v_z*p_x) - v_y // c = K*p_x*p_z - p_y // // Out of the possible two solutions you must choose // the one that gives a positive rho(z=0). // // fCurStatWithV.ResetfDone(validate, &gp, &gv); if (fCurStatWithV.IsDone()) { G4int i; for (i=0; i DBL_MIN) { // single solution distance[0] = - c / b; xx[0] = p + distance[0]*v; gxx[0] = ComputeGlobalPoint(xx[0]); if (validate == kValidateWithTol) { areacode[0] = GetAreaCode(xx[0]); if (!IsOutside(areacode[0])) { if (distance[0] >= 0) isvalid[0] = true; } } else if (validate == kValidateWithoutTol) { areacode[0] = GetAreaCode(xx[0], false); if (IsInside(areacode[0])) { if (distance[0] >= 0) isvalid[0] = true; } } else { // kDontValidate // we must omit x(rho,z) = rho(z=0) < 0 if (xx[0].x() > 0) { areacode[0] = sInside; if (distance[0] >= 0) isvalid[0] = true; } else { distance[0] = kInfinity; fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], isvalid[0], 0, validate, &gp, &gv); return 0; } } fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], isvalid[0], 1, validate, &gp, &gv); return 1; } else { // if a=b=0 , v.y=0 and (v.x=0 && p.x=0) or (v.z=0 && p.z=0) . // if v.x=0 && p.x=0, no intersection unless p is on z-axis // (in that case, v is paralell to surface). // if v.z=0 && p.z=0, no intersection unless p is on x-axis // (in that case, v is paralell to surface). // return distance = infinity. fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], isvalid[0], 0, validate, &gp, &gv); return 0; } } else if (D > DBL_MIN) { // double solutions D = std::sqrt(D); G4double factor = 0.5/a; G4double tmpdist[2] = {kInfinity, kInfinity}; G4ThreeVector tmpxx[2]; G4int tmpareacode[2] = {sOutside, sOutside}; G4bool tmpisvalid[2] = {false, false}; G4int i; for (i=0; i<2; i++) { G4double bminusD = - b - D; // protection against round off error //G4double protection = 1.0e-6; G4double protection = 0; if ( b * D < 0 && std::fabs(bminusD / D) < protection ) { G4double acovbb = (a*c)/(b*b); tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb)); } else { tmpdist[i] = factor * bminusD; } D = -D; tmpxx[i] = p + tmpdist[i]*v; if (validate == kValidateWithTol) { tmpareacode[i] = GetAreaCode(tmpxx[i]); if (!IsOutside(tmpareacode[i])) { if (tmpdist[i] >= 0) tmpisvalid[i] = true; continue; } } else if (validate == kValidateWithoutTol) { tmpareacode[i] = GetAreaCode(tmpxx[i], false); if (IsInside(tmpareacode[i])) { if (tmpdist[i] >= 0) tmpisvalid[i] = true; continue; } } else { // kDontValidate // we must choose x(rho,z) = rho(z=0) > 0 if (tmpxx[i].x() > 0) { tmpareacode[i] = sInside; if (tmpdist[i] >= 0) tmpisvalid[i] = true; continue; } else { tmpdist[i] = kInfinity; continue; } } } if (tmpdist[0] <= tmpdist[1]) { distance[0] = tmpdist[0]; distance[1] = tmpdist[1]; xx[0] = tmpxx[0]; xx[1] = tmpxx[1]; gxx[0] = ComputeGlobalPoint(tmpxx[0]); gxx[1] = ComputeGlobalPoint(tmpxx[1]); areacode[0] = tmpareacode[0]; areacode[1] = tmpareacode[1]; isvalid[0] = tmpisvalid[0]; isvalid[1] = tmpisvalid[1]; } else { distance[0] = tmpdist[1]; distance[1] = tmpdist[0]; xx[0] = tmpxx[1]; xx[1] = tmpxx[0]; gxx[0] = ComputeGlobalPoint(tmpxx[1]); gxx[1] = ComputeGlobalPoint(tmpxx[0]); areacode[0] = tmpareacode[1]; areacode[1] = tmpareacode[0]; isvalid[0] = tmpisvalid[1]; isvalid[1] = tmpisvalid[0]; } fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], isvalid[0], 2, validate, &gp, &gv); fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1], isvalid[1], 2, validate, &gp, &gv); // protection against roundoff error for (G4int k=0; k<2; k++) { if (!isvalid[k]) continue; G4ThreeVector xxonsurface(xx[k].x(), fKappa * std::fabs(xx[k].x()) * xx[k].z() , xx[k].z()); G4double deltaY = (xx[k] - xxonsurface).mag(); if ( deltaY > 0.5*kCarTolerance ) { G4int maxcount = 10; G4int l; G4double lastdeltaY = deltaY; for (l=0; l lastdeltaY) { } gxx[k] = ComputeGlobalPoint(xx[k]); if (deltaY <= 0.5*kCarTolerance) { break; } xxonsurface.set(xx[k].x(), fKappa * std::fabs(xx[k].x()) * xx[k].z(), xx[k].z()); } if (l == maxcount) { G4cerr << "ERROR - G4TwistTubsSide::DistanceToSurface(p,v)" << G4endl << " maxloop count " << maxcount << G4endl; G4Exception("G4TwistTubsFlatSide::DistanceToSurface(p,v)", "InvalidSetup", FatalException, "Exceeded maxloop count!"); } } } return 2; } else { // if D<0, no solution // if D=0, just grazing the surfaces, return kInfinity fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], isvalid[0], 0, validate, &gp, &gv); return 0; } G4Exception("G4TwistTubsSide::DistanceToSurface(p,v)", "InvalidCondition", FatalException, "Illegal operation !"); return 1; } //===================================================================== //* DistanceToSurface ------------------------------------------------- G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector &gp, G4ThreeVector gxx[], G4double distance[], G4int areacode[]) { fCurStat.ResetfDone(kDontValidate, &gp); G4int i = 0; if (fCurStat.IsDone()) { for (i=0; i= 0 ? 1 : -1); // // special case! // If p is on surface, or // p is on z-axis, // return here immediatery. // G4ThreeVector lastgxx[2]; G4double distfromlast[2]; for (i=0; i<2; i++) { lastgxx[i] = fCurStatWithV.GetXX(i); distfromlast[i] = (gp - lastgxx[i]).mag(); } if ((gp - lastgxx[0]).mag() < halftol || (gp - lastgxx[1]).mag() < halftol) { // last winner, or last poststep point is on the surface. xx = p; distance[0] = 0; gxx[0] = gp; G4bool isvalid = true; fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], isvalid, 1, kDontValidate, &gp); return 1; } if (p.getRho() == 0) { // p is on z-axis. Namely, p is on twisted surface (invalid area). // We must return here, however, returning distance to x-minimum // boundary is better than return 0-distance. // G4bool isvalid = true; if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) { distance[0] = DistanceToBoundary(sAxis0 & sAxisMin, xx, p); areacode[0] = sInside; } else { distance[0] = 0; xx.set(0., 0., 0.); } gxx[0] = ComputeGlobalPoint(xx); fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], isvalid, 0, kDontValidate, &gp); return 1; } // // special case end // // set corner points of quadrangle try area ... G4ThreeVector A; // foot of normal from p to boundary of sAxis0 & sAxisMin G4ThreeVector C; // foot of normal from p to boundary of sAxis0 & sAxisMax G4ThreeVector B; // point on boundary sAxis0 & sAxisMax at z = A.z() G4ThreeVector D; // point on boundary sAxis0 & sAxisMin at z = C.z() G4double distToA; // distance from p to A G4double distToC; // distance from p to C distToA = DistanceToBoundary(sAxis0 & sAxisMin, A, p); distToC = DistanceToBoundary(sAxis0 & sAxisMax, C, p); // is p.z between a.z and c.z? // p.z must be bracketed a.z and c.z. if (A.z() > C.z()) { if (p.z() > A.z()) { A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p); } else if (p.z() < C.z()) { C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p); } } else { if (p.z() > C.z()) { C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p); } else if (p.z() < A.z()) { A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p); } } G4ThreeVector d[2]; // direction vectors of boundary G4ThreeVector x0[2]; // foot of normal from line to p G4int btype[2]; // boundary type for (i=0; i<2; i++) { if (i == 0) { GetBoundaryParameters((sAxis0 & sAxisMax), d[i], x0[i], btype[i]); B = x0[i] + ((A.z() - x0[i].z()) / d[i].z()) * d[i]; // x0 + t*d , d is direction unit vector. } else { GetBoundaryParameters((sAxis0 & sAxisMin), d[i], x0[i], btype[i]); D = x0[i] + ((C.z() - x0[i].z()) / d[i].z()) * d[i]; } } // In order to set correct diagonal, swap A and D, C and B if needed. G4ThreeVector pt(p.x(), p.y(), 0.); G4double rc = std::fabs(p.x()); G4ThreeVector surfacevector(rc, rc * fKappa * p.z(), 0.); G4int pside = AmIOnLeftSide(pt, surfacevector); G4double test = (A.z() - C.z()) * parity * pside; if (test == 0) { if (pside == 0) { // p is on surface. xx = p; distance[0] = 0; gxx[0] = gp; G4bool isvalid = true; fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], isvalid, 1, kDontValidate, &gp); return 1; } else { // A.z = C.z(). return distance to line. d[0] = C - A; distance[0] = DistanceToLine(p, A, d[0], xx); areacode[0] = sInside; gxx[0] = ComputeGlobalPoint(xx); G4bool isvalid = true; fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], isvalid, 1, kDontValidate, &gp); return 1; } } else if (test < 0) { // wrong diagonal. vector AC is crossing the surface! // swap A and D, C and B G4ThreeVector tmp; tmp = A; A = D; D = tmp; tmp = C; C = B; B = tmp; } else { // correct diagonal. nothing to do. } // Now, we chose correct diaglnal. // First try. divide quadrangle into double triangle by diagonal and // calculate distance to both surfaces. G4ThreeVector xxacb; // foot of normal from plane ACB to p G4ThreeVector nacb; // normal of plane ACD G4ThreeVector xxcad; // foot of normal from plane CAD to p G4ThreeVector ncad; // normal of plane CAD G4ThreeVector AB(A.x(), A.y(), 0); G4ThreeVector DC(C.x(), C.y(), 0); G4double distToACB = G4VTwistSurface::DistanceToPlane(p, A, C-A, AB, xxacb, nacb) * parity; G4double distToCAD = G4VTwistSurface::DistanceToPlane(p, C, C-A, DC, xxcad, ncad) * parity; // if calculated distance = 0, return if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol) { xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad); areacode[0] = sInside; gxx[0] = ComputeGlobalPoint(xx); distance[0] = 0; G4bool isvalid = true; fCurStat.SetCurrentStatus(0, gxx[0], distance[0] , areacode[0], isvalid, 1, kDontValidate, &gp); return 1; } if (distToACB * distToCAD > 0 && distToACB < 0) { // both distToACB and distToCAD are negative. // divide quadrangle into double triangle by diagonal G4ThreeVector normal; distance[0] = DistanceToPlane(p, A, B, C, D, parity, xx, normal); } else { if (distToACB * distToCAD > 0) { // both distToACB and distToCAD are positive. // Take smaller one. if (distToACB <= distToCAD) { distance[0] = distToACB; xx = xxacb; } else { distance[0] = distToCAD; xx = xxcad; } } else { // distToACB * distToCAD is negative. // take positive one if (distToACB > 0) { distance[0] = distToACB; xx = xxacb; } else { distance[0] = distToCAD; xx = xxcad; } } } areacode[0] = sInside; gxx[0] = ComputeGlobalPoint(xx); G4bool isvalid = true; fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], isvalid, 1, kDontValidate, &gp); return 1; } //===================================================================== //* DistanceToPlane --------------------------------------------------- G4double G4TwistTubsSide::DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &A, const G4ThreeVector &B, const G4ThreeVector &C, const G4ThreeVector &D, const G4int parity, G4ThreeVector &xx, G4ThreeVector &n) { static const G4double halftol = 0.5 * kCarTolerance; G4ThreeVector M = 0.5*(A + B); G4ThreeVector N = 0.5*(C + D); G4ThreeVector xxanm; // foot of normal from p to plane ANM G4ThreeVector nanm; // normal of plane ANM G4ThreeVector xxcmn; // foot of normal from p to plane CMN G4ThreeVector ncmn; // normal of plane CMN G4double distToanm = G4VTwistSurface::DistanceToPlane(p, A, (N - A), (M - A), xxanm, nanm) * parity; G4double distTocmn = G4VTwistSurface::DistanceToPlane(p, C, (M - C), (N - C), xxcmn, ncmn) * parity; // if p is behind of both surfaces, abort. if (distToanm * distTocmn > 0 && distToanm < 0) { G4Exception("G4TwistTubsSide::DistanceToPlane()", "InvalidCondition", FatalException, "Point p is behind the surfaces."); } // if p is on surface, return 0. if (std::fabs(distToanm) <= halftol) { xx = xxanm; n = nanm * parity; return 0; } else if (std::fabs(distTocmn) <= halftol) { xx = xxcmn; n = ncmn * parity; return 0; } if (distToanm <= distTocmn) { if (distToanm > 0) { // both distanses are positive. take smaller one. xx = xxanm; n = nanm * parity; return distToanm; } else { // take -ve distance and call the function recursively. return DistanceToPlane(p, A, M, N, D, parity, xx, n); } } else { if (distTocmn > 0) { // both distanses are positive. take smaller one. xx = xxcmn; n = ncmn * parity; return distTocmn; } else { // take -ve distance and call the function recursively. return DistanceToPlane(p, C, N, M, B, parity, xx, n); } } } //===================================================================== //* GetAreaCode ------------------------------------------------------- G4int G4TwistTubsSide::GetAreaCode(const G4ThreeVector &xx, G4bool withTol) { // We must use the function in local coordinate system. // See the description of DistanceToSurface(p,v). static const G4double ctol = 0.5 * kCarTolerance; G4int areacode = sInside; if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) { G4int xaxis = 0; G4int zaxis = 1; if (withTol) { G4bool isoutside = false; // test boundary of xaxis if (xx.x() < fAxisMin[xaxis] + ctol) { areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary; if (xx.x() <= fAxisMin[xaxis] - ctol) isoutside = true; } else if (xx.x() > fAxisMax[xaxis] - ctol) { areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary; if (xx.x() >= fAxisMax[xaxis] + ctol) isoutside = true; } // test boundary of z-axis if (xx.z() < fAxisMin[zaxis] + ctol) { areacode |= (sAxis1 & (sAxisZ | sAxisMin)); if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner. else areacode |= sBoundary; if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true; } else if (xx.z() > fAxisMax[zaxis] - ctol) { areacode |= (sAxis1 & (sAxisZ | sAxisMax)); if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner. else areacode |= sBoundary; if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true; } // if isoutside = true, clear inside bit. // if not on boundary, add axis information. if (isoutside) { G4int tmpareacode = areacode & (~sInside); areacode = tmpareacode; } else if ((areacode & sBoundary) != sBoundary) { areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ); } } else { // boundary of x-axis if (xx.x() < fAxisMin[xaxis] ) { areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary; } else if (xx.x() > fAxisMax[xaxis]) { areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary; } // boundary of z-axis if (xx.z() < fAxisMin[zaxis]) { areacode |= (sAxis1 & (sAxisZ | sAxisMin)); if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner. else areacode |= sBoundary; } else if (xx.z() > fAxisMax[zaxis]) { areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ; if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner. else areacode |= sBoundary; } if ((areacode & sBoundary) != sBoundary) { areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ); } } return areacode; } else { G4Exception("G4TwistTubsSide::GetAreaCode()", "NotImplemented", FatalException, "Feature NOT implemented !"); } return areacode; } //===================================================================== //* SetCorners( arglist ) ------------------------------------------------- void G4TwistTubsSide::SetCorners( G4double endInnerRad[2], G4double endOuterRad[2], G4double endPhi[2], G4double endZ[2]) { // Set Corner points in local coodinate. if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) { G4int zmin = 0 ; // at -ve z G4int zmax = 1 ; // at +ve z G4double x, y, z; // corner of Axis0min and Axis1min x = endInnerRad[zmin]*std::cos(endPhi[zmin]); y = endInnerRad[zmin]*std::sin(endPhi[zmin]); z = endZ[zmin]; SetCorner(sC0Min1Min, x, y, z); // corner of Axis0max and Axis1min x = endOuterRad[zmin]*std::cos(endPhi[zmin]); y = endOuterRad[zmin]*std::sin(endPhi[zmin]); z = endZ[zmin]; SetCorner(sC0Max1Min, x, y, z); // corner of Axis0max and Axis1max x = endOuterRad[zmax]*std::cos(endPhi[zmax]); y = endOuterRad[zmax]*std::sin(endPhi[zmax]); z = endZ[zmax]; SetCorner(sC0Max1Max, x, y, z); // corner of Axis0min and Axis1max x = endInnerRad[zmax]*std::cos(endPhi[zmax]); y = endInnerRad[zmax]*std::sin(endPhi[zmax]); z = endZ[zmax]; SetCorner(sC0Min1Max, x, y, z); } else { G4cerr << "ERROR - G4TwistTubsFlatSide::SetCorners()" << G4endl << " fAxis[0] = " << fAxis[0] << G4endl << " fAxis[1] = " << fAxis[1] << G4endl; G4Exception("G4TwistTubsSide::SetCorners()", "NotImplemented", FatalException, "Feature NOT implemented !"); } } //===================================================================== //* SetCorners() ------------------------------------------------------ void G4TwistTubsSide::SetCorners() { G4Exception("G4TwistTubsSide::SetCorners()", "NotImplemented", FatalException, "Method NOT implemented !"); } //===================================================================== //* SetBoundaries() --------------------------------------------------- void G4TwistTubsSide::SetBoundaries() { // Set direction-unit vector of boundary-lines in local coodinate. // G4ThreeVector direction; if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) { // sAxis0 & sAxisMin direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min); direction = direction.unit(); SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction, GetCorner(sC0Min1Min), sAxisZ) ; // sAxis0 & sAxisMax direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min); direction = direction.unit(); SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction, GetCorner(sC0Max1Min), sAxisZ); // sAxis1 & sAxisMin direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min); direction = direction.unit(); SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction, GetCorner(sC0Min1Min), sAxisX); // sAxis1 & sAxisMax direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max); direction = direction.unit(); SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction, GetCorner(sC0Min1Max), sAxisX); } else { G4cerr << "ERROR - G4TwistTubsFlatSide::SetBoundaries()" << G4endl << " fAxis[0] = " << fAxis[0] << G4endl << " fAxis[1] = " << fAxis[1] << G4endl; G4Exception("G4TwistTubsSide::SetCorners()", "NotImplemented", FatalException, "Feature NOT implemented !"); } } //===================================================================== //* GetFacets() ------------------------------------------------------- void G4TwistTubsSide::GetFacets( G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside ) { G4double z ; // the two parameters for the surface equation G4double x,xmin,xmax ; G4ThreeVector p ; // a point on the surface, given by (z,u) G4int nnode ; G4int nface ; // calculate the (n-1)*(m-1) vertices G4int i,j ; for ( i = 0 ; i