Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // G4TwistedTubs implementation << 27 // 26 // 28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu << 27 // $Id: G4TwistedTubs.cc,v 1.24 2007/05/18 07:39:56 gcosmo Exp $ 29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), << 28 // GEANT4 tag $Name: geant4-09-00 $ 30 // from original version in Jupi << 29 // >> 30 // >> 31 // -------------------------------------------------------------------- >> 32 // GEANT 4 class source file >> 33 // >> 34 // >> 35 // G4TwistTubsSide.cc >> 36 // >> 37 // Author: >> 38 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp) >> 39 // >> 40 // History: >> 41 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4 >> 42 // from original version in Jupiter-2.5.02 application. 31 // ------------------------------------------- 43 // -------------------------------------------------------------------- 32 44 33 #include "G4TwistedTubs.hh" 45 #include "G4TwistedTubs.hh" 34 46 35 #include "G4GeomTools.hh" << 36 #include "G4PhysicalConstants.hh" << 37 #include "G4SystemOfUnits.hh" << 38 #include "G4GeometryTolerance.hh" 47 #include "G4GeometryTolerance.hh" 39 #include "G4VoxelLimits.hh" 48 #include "G4VoxelLimits.hh" 40 #include "G4AffineTransform.hh" 49 #include "G4AffineTransform.hh" 41 #include "G4BoundingEnvelope.hh" << 50 #include "G4SolidExtentList.hh" 42 #include "G4ClippablePolygon.hh" 51 #include "G4ClippablePolygon.hh" 43 #include "G4VPVParameterisation.hh" 52 #include "G4VPVParameterisation.hh" 44 #include "meshdefs.hh" 53 #include "meshdefs.hh" 45 54 46 #include "G4VGraphicsScene.hh" 55 #include "G4VGraphicsScene.hh" 47 #include "G4Polyhedron.hh" 56 #include "G4Polyhedron.hh" 48 #include "G4VisExtent.hh" 57 #include "G4VisExtent.hh" >> 58 #include "G4NURBS.hh" >> 59 #include "G4NURBStube.hh" >> 60 #include "G4NURBScylinder.hh" >> 61 #include "G4NURBStubesector.hh" 49 62 50 #include "Randomize.hh" 63 #include "Randomize.hh" 51 64 52 #include "G4AutoLock.hh" << 53 << 54 namespace << 55 { << 56 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE << 57 } << 58 << 59 << 60 //============================================ 65 //===================================================================== 61 //* constructors ----------------------------- 66 //* constructors ------------------------------------------------------ 62 67 63 G4TwistedTubs::G4TwistedTubs(const G4String& p << 68 G4TwistedTubs::G4TwistedTubs(const G4String &pname, 64 G4double t 69 G4double twistedangle, 65 G4double e 70 G4double endinnerrad, 66 G4double e 71 G4double endouterrad, 67 G4double h 72 G4double halfzlen, 68 G4double d 73 G4double dphi) 69 : G4VSolid(pname), fDPhi(dphi), 74 : G4VSolid(pname), fDPhi(dphi), 70 fLowerEndcap(nullptr), fUpperEndcap(nullp << 75 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0), 71 fFormerTwisted(nullptr), fInnerHype(nullp << 76 fFormerTwisted(0), fInnerHype(0), fOuterHype(0), >> 77 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0) 72 { 78 { 73 if (endinnerrad < DBL_MIN) 79 if (endinnerrad < DBL_MIN) 74 { 80 { 75 G4Exception("G4TwistedTubs::G4TwistedTub << 81 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup", 76 FatalErrorInArgument, "Inval << 82 FatalException, "Invalid end-inner-radius!"); 77 } 83 } 78 84 79 G4double sinhalftwist = std::sin(0.5 * twis 85 G4double sinhalftwist = std::sin(0.5 * twistedangle); 80 86 81 G4double endinnerradX = endinnerrad * sinha 87 G4double endinnerradX = endinnerrad * sinhalftwist; 82 G4double innerrad = std::sqrt( endinner 88 G4double innerrad = std::sqrt( endinnerrad * endinnerrad 83 - endinnerrad 89 - endinnerradX * endinnerradX ); 84 90 85 G4double endouterradX = endouterrad * sinha 91 G4double endouterradX = endouterrad * sinhalftwist; 86 G4double outerrad = std::sqrt( endouter 92 G4double outerrad = std::sqrt( endouterrad * endouterrad 87 - endouterrad 93 - endouterradX * endouterradX ); 88 94 89 // temporary treatment!! 95 // temporary treatment!! 90 SetFields(twistedangle, innerrad, outerrad, 96 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen); 91 CreateSurfaces(); 97 CreateSurfaces(); 92 } 98 } 93 99 94 G4TwistedTubs::G4TwistedTubs(const G4String& p << 100 G4TwistedTubs::G4TwistedTubs(const G4String &pname, 95 G4double t 101 G4double twistedangle, 96 G4double e 102 G4double endinnerrad, 97 G4double e 103 G4double endouterrad, 98 G4double h 104 G4double halfzlen, 99 G4int n 105 G4int nseg, 100 G4double t 106 G4double totphi) 101 : G4VSolid(pname), 107 : G4VSolid(pname), 102 fLowerEndcap(nullptr), fUpperEndcap(nullp << 108 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0), 103 fFormerTwisted(nullptr), fInnerHype(nullp << 109 fFormerTwisted(0), fInnerHype(0), fOuterHype(0), >> 110 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0) 104 { 111 { 105 112 106 if (nseg == 0) << 113 if (!nseg) 107 { 114 { 108 std::ostringstream message; << 115 G4cerr << "ERROR - G4TwistedTubs::G4TwistedTubs()" << G4endl 109 message << "Invalid number of segments." << 116 << " Invalid nseg. nseg = " << nseg << G4endl; 110 << " nseg = " << nseg; << 111 G4Exception("G4TwistedTubs::G4TwistedTub << 112 FatalErrorInArgument, messag << 113 } 117 } 114 if (totphi == DBL_MIN || endinnerrad < DBL_ 118 if (totphi == DBL_MIN || endinnerrad < DBL_MIN) 115 { 119 { 116 G4Exception("G4TwistedTubs::G4TwistedTub << 120 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup", 117 FatalErrorInArgument, "Invalid << 121 FatalException, "Invalid total-phi or end-inner-radius!"); 118 } << 122 } 119 123 120 G4double sinhalftwist = std::sin(0.5 * twis 124 G4double sinhalftwist = std::sin(0.5 * twistedangle); 121 125 122 G4double endinnerradX = endinnerrad * sinha 126 G4double endinnerradX = endinnerrad * sinhalftwist; 123 G4double innerrad = std::sqrt( endinner 127 G4double innerrad = std::sqrt( endinnerrad * endinnerrad 124 - endinnerrad 128 - endinnerradX * endinnerradX ); 125 129 126 G4double endouterradX = endouterrad * sinha 130 G4double endouterradX = endouterrad * sinhalftwist; 127 G4double outerrad = std::sqrt( endouter 131 G4double outerrad = std::sqrt( endouterrad * endouterrad 128 - endouterrad 132 - endouterradX * endouterradX ); 129 133 130 // temporary treatment!! 134 // temporary treatment!! 131 fDPhi = totphi / nseg; 135 fDPhi = totphi / nseg; 132 SetFields(twistedangle, innerrad, outerrad, 136 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen); 133 CreateSurfaces(); 137 CreateSurfaces(); 134 } 138 } 135 139 136 G4TwistedTubs::G4TwistedTubs(const G4String& p << 140 G4TwistedTubs::G4TwistedTubs(const G4String &pname, 137 G4double t 141 G4double twistedangle, 138 G4double i 142 G4double innerrad, 139 G4double o 143 G4double outerrad, 140 G4double n 144 G4double negativeEndz, 141 G4double p 145 G4double positiveEndz, 142 G4double d 146 G4double dphi) 143 : G4VSolid(pname), fDPhi(dphi), 147 : G4VSolid(pname), fDPhi(dphi), 144 fLowerEndcap(nullptr), fUpperEndcap(nullp << 148 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0), 145 fFormerTwisted(nullptr), fInnerHype(nullp << 149 fFormerTwisted(0), fInnerHype(0), fOuterHype(0), >> 150 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0) 146 { 151 { 147 if (innerrad < DBL_MIN) 152 if (innerrad < DBL_MIN) 148 { 153 { 149 G4Exception("G4TwistedTubs::G4TwistedTub << 154 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup", 150 FatalErrorInArgument, "Inval << 155 FatalException, "Invalid end-inner-radius!"); 151 } 156 } 152 157 153 SetFields(twistedangle, innerrad, outerrad, 158 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz); 154 CreateSurfaces(); 159 CreateSurfaces(); 155 } 160 } 156 161 157 G4TwistedTubs::G4TwistedTubs(const G4String& p << 162 G4TwistedTubs::G4TwistedTubs(const G4String &pname, 158 G4double t 163 G4double twistedangle, 159 G4double i 164 G4double innerrad, 160 G4double o 165 G4double outerrad, 161 G4double n 166 G4double negativeEndz, 162 G4double p 167 G4double positiveEndz, 163 G4int n 168 G4int nseg, 164 G4double t 169 G4double totphi) 165 : G4VSolid(pname), 170 : G4VSolid(pname), 166 fLowerEndcap(nullptr), fUpperEndcap(nullp << 171 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0), 167 fFormerTwisted(nullptr), fInnerHype(nullp << 172 fFormerTwisted(0), fInnerHype(0), fOuterHype(0), >> 173 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0) 168 { 174 { 169 if (nseg == 0) << 175 if (!nseg) 170 { 176 { 171 std::ostringstream message; << 177 G4cerr << "ERROR - G4TwistedTubs::G4TwistedTubs()" << G4endl 172 message << "Invalid number of segments." << 178 << " Invalid nseg. nseg = " << nseg << G4endl; 173 << " nseg = " << nseg; << 174 G4Exception("G4TwistedTubs::G4TwistedTub << 175 FatalErrorInArgument, messag << 176 } 179 } 177 if (totphi == DBL_MIN || innerrad < DBL_MIN 180 if (totphi == DBL_MIN || innerrad < DBL_MIN) 178 { 181 { 179 G4Exception("G4TwistedTubs::G4TwistedTub << 182 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup", 180 FatalErrorInArgument, "Invalid << 183 FatalException, "Invalid total-phi or end-inner-radius!"); 181 } 184 } 182 185 183 fDPhi = totphi / nseg; 186 fDPhi = totphi / nseg; 184 SetFields(twistedangle, innerrad, outerrad, 187 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz); 185 CreateSurfaces(); 188 CreateSurfaces(); 186 } 189 } 187 190 188 //============================================ 191 //===================================================================== 189 //* Fake default constructor ----------------- 192 //* Fake default constructor ------------------------------------------ 190 193 191 G4TwistedTubs::G4TwistedTubs( __void__& a ) 194 G4TwistedTubs::G4TwistedTubs( __void__& a ) 192 : G4VSolid(a), << 195 : G4VSolid(a), fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0), 193 fLowerEndcap(nullptr), fUpperEndcap(nullpt << 196 fFormerTwisted(0), fInnerHype(0), fOuterHype(0), fCubicVolume(0.), 194 fLatterTwisted(nullptr), fFormerTwisted(nu << 197 fSurfaceArea(0.), fpPolyhedron(0) 195 fInnerHype(nullptr), fOuterHype(nullptr) << 196 { 198 { 197 } 199 } 198 200 199 //============================================ 201 //===================================================================== 200 //* destructor ------------------------------- 202 //* destructor -------------------------------------------------------- 201 203 202 G4TwistedTubs::~G4TwistedTubs() 204 G4TwistedTubs::~G4TwistedTubs() 203 { 205 { 204 delete fLowerEndcap; << 206 if (fLowerEndcap) { delete fLowerEndcap; } 205 delete fUpperEndcap; << 207 if (fUpperEndcap) { delete fUpperEndcap; } 206 delete fLatterTwisted; << 208 if (fLatterTwisted) { delete fLatterTwisted; } 207 delete fFormerTwisted; << 209 if (fFormerTwisted) { delete fFormerTwisted; } 208 delete fInnerHype; << 210 if (fInnerHype) { delete fInnerHype; } 209 delete fOuterHype; << 211 if (fOuterHype) { delete fOuterHype; } 210 delete fpPolyhedron; fpPolyhedron = nullptr << 212 if (fpPolyhedron) { delete fpPolyhedron; } 211 } << 212 << 213 //============================================ << 214 //* Copy constructor ------------------------- << 215 << 216 G4TwistedTubs::G4TwistedTubs(const G4TwistedTu << 217 : G4VSolid(rhs), fPhiTwist(rhs.fPhiTwist), << 218 fInnerRadius(rhs.fInnerRadius), fOuterRadi << 219 fDPhi(rhs.fDPhi), fZHalfLength(rhs.fZHalfL << 220 fInnerStereo(rhs.fInnerStereo), fOuterSter << 221 fTanInnerStereo(rhs.fTanInnerStereo), fTan << 222 fKappa(rhs.fKappa), fInnerRadius2(rhs.fInn << 223 fOuterRadius2(rhs.fOuterRadius2), fTanInne << 224 fTanOuterStereo2(rhs.fTanOuterStereo2), << 225 fLowerEndcap(nullptr), fUpperEndcap(nullpt << 226 fInnerHype(nullptr), fOuterHype(nullptr), << 227 fCubicVolume(rhs.fCubicVolume), fSurfaceAr << 228 { << 229 for (auto i=0; i<2; ++i) << 230 { << 231 fEndZ[i] = rhs.fEndZ[i]; << 232 fEndInnerRadius[i] = rhs.fEndInnerRadius[i << 233 fEndOuterRadius[i] = rhs.fEndOuterRadius[i << 234 fEndPhi[i] = rhs.fEndPhi[i]; << 235 fEndZ2[i] = rhs.fEndZ2[i]; << 236 } << 237 CreateSurfaces(); << 238 } << 239 << 240 << 241 //============================================ << 242 //* Assignment operator ---------------------- << 243 << 244 G4TwistedTubs& G4TwistedTubs::operator = (cons << 245 { << 246 // Check assignment to self << 247 // << 248 if (this == &rhs) { return *this; } << 249 << 250 // Copy base class data << 251 // << 252 G4VSolid::operator=(rhs); << 253 << 254 // Copy data << 255 // << 256 fPhiTwist= rhs.fPhiTwist; << 257 fInnerRadius= rhs.fInnerRadius; fOuterRadiu << 258 fDPhi= rhs.fDPhi; fZHalfLength= rhs.fZHalfL << 259 fInnerStereo= rhs.fInnerStereo; fOuterStere << 260 fTanInnerStereo= rhs.fTanInnerStereo; fTanO << 261 fKappa= rhs.fKappa; fInnerRadius2= rhs.fInn << 262 fOuterRadius2= rhs.fOuterRadius2; fTanInner << 263 fTanOuterStereo2= rhs.fTanOuterStereo2; << 264 fLowerEndcap= fUpperEndcap= fLatterTwisted= << 265 fInnerHype= fOuterHype= nullptr; << 266 fCubicVolume= rhs.fCubicVolume; fSurfaceAre << 267 << 268 for (auto i=0; i<2; ++i) << 269 { << 270 fEndZ[i] = rhs.fEndZ[i]; << 271 fEndInnerRadius[i] = rhs.fEndInnerRadius[ << 272 fEndOuterRadius[i] = rhs.fEndOuterRadius[ << 273 fEndPhi[i] = rhs.fEndPhi[i]; << 274 fEndZ2[i] = rhs.fEndZ2[i]; << 275 } << 276 << 277 CreateSurfaces(); << 278 fRebuildPolyhedron = false; << 279 delete fpPolyhedron; fpPolyhedron = nullptr << 280 << 281 return *this; << 282 } 213 } 283 214 284 //============================================ 215 //===================================================================== 285 //* ComputeDimensions ------------------------ 216 //* ComputeDimensions ------------------------------------------------- 286 217 287 void G4TwistedTubs::ComputeDimensions(G4VPVPar 218 void G4TwistedTubs::ComputeDimensions(G4VPVParameterisation* /* p */ , 288 const G4 219 const G4int /* n */ , 289 const G4 220 const G4VPhysicalVolume* /* pRep */ ) 290 { 221 { 291 G4Exception("G4TwistedTubs::ComputeDimension 222 G4Exception("G4TwistedTubs::ComputeDimensions()", 292 "GeomSolids0001", FatalException << 223 "NotSupported", FatalException, 293 "G4TwistedTubs does not support 224 "G4TwistedTubs does not support Parameterisation."); 294 } 225 } 295 226 >> 227 296 //============================================ 228 //===================================================================== 297 //* BoundingLimits --------------------------- << 229 //* CalculateExtent --------------------------------------------------- 298 230 299 void G4TwistedTubs::BoundingLimits(G4ThreeVect << 231 G4bool G4TwistedTubs::CalculateExtent( const EAxis axis, 300 G4ThreeVect << 232 const G4VoxelLimits &voxelLimit, >> 233 const G4AffineTransform &transform, >> 234 G4double &min, >> 235 G4double &max ) const 301 { 236 { 302 // Find bounding tube << 303 G4double rmin = GetInnerRadius(); << 304 G4double rmax = GetEndOuterRadius(); << 305 237 306 G4double zmin = std::min(GetEndZ(0), GetEndZ << 238 G4SolidExtentList extentList( axis, voxelLimit ); 307 G4double zmax = std::max(GetEndZ(0), GetEndZ << 239 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? >> 240 fEndOuterRadius[0] : fEndOuterRadius[1]); >> 241 G4double maxEndInnerRad = (fEndInnerRadius[0] > fEndInnerRadius[1] ? >> 242 fEndInnerRadius[0] : fEndInnerRadius[1]); >> 243 G4double maxphi = (std::fabs(fEndPhi[0]) > std::fabs(fEndPhi[1]) ? >> 244 std::fabs(fEndPhi[0]) : std::fabs(fEndPhi[1])); >> 245 >> 246 // >> 247 // Choose phi size of our segment(s) based on constants as >> 248 // defined in meshdefs.hh >> 249 // >> 250 // G4int numPhi = kMaxMeshSections; >> 251 G4double sigPhi = 2*maxphi + fDPhi; >> 252 G4double rFudge = 1.0/std::cos(0.5*sigPhi); >> 253 G4double fudgeEndOuterRad = rFudge * maxEndOuterRad; >> 254 >> 255 // >> 256 // We work around in phi building polygons along the way. >> 257 // As a reasonable compromise between accuracy and >> 258 // complexity (=cpu time), the following facets are chosen: >> 259 // >> 260 // 1. If fOuterRadius/maxEndOuterRad > 0.95, approximate >> 261 // the outer surface as a cylinder, and use one >> 262 // rectangular polygon (0-1) to build its mesh. >> 263 // >> 264 // Otherwise, use two trapazoidal polygons that >> 265 // meet at z = 0 (0-4-1) >> 266 // >> 267 // 2. If there is no inner surface, then use one >> 268 // polygon for each entire endcap. (0) and (1) >> 269 // >> 270 // Otherwise, use a trapazoidal polygon for each >> 271 // phi segment of each endcap. (0-2) and (1-3) >> 272 // >> 273 // 3. For the inner surface, if fInnerRadius/maxEndInnerRad > 0.95, >> 274 // approximate the inner surface as a cylinder of >> 275 // radius fInnerRadius and use one rectangular polygon >> 276 // to build each phi segment of its mesh. (2-3) >> 277 // >> 278 // Otherwise, use one rectangular polygon centered >> 279 // at z = 0 (5-6) and two connecting trapazoidal polygons >> 280 // for each phi segment (2-5) and (3-6). >> 281 // 308 282 309 G4double dphi = 0.5*GetDPhi(); << 283 G4bool splitOuter = (fOuterRadius/maxEndOuterRad < 0.95); 310 G4double sphi = std::min(GetEndPhi(0), GetEn << 284 G4bool splitInner = (fInnerRadius/maxEndInnerRad < 0.95); 311 G4double ephi = std::max(GetEndPhi(0), GetEn << 285 312 G4double totalphi = ephi - sphi; << 286 // >> 287 // Vertex assignments (v and w arrays) >> 288 // [0] and [1] are mandatory >> 289 // the rest are optional >> 290 // >> 291 // + - >> 292 // [0]------[4]------[1] <--- outer radius >> 293 // | | >> 294 // | | >> 295 // [2]---[5]---[6]---[3] <--- inner radius >> 296 // >> 297 >> 298 G4ClippablePolygon endPoly1, endPoly2; >> 299 >> 300 G4double phimax = maxphi + 0.5*fDPhi; >> 301 G4double phimin = - phimax; >> 302 >> 303 G4ThreeVector v0, v1, v2, v3, v4, v5, v6; // -ve phi verticies for polygon >> 304 G4ThreeVector w0, w1, w2, w3, w4, w5, w6; // +ve phi verticies for polygon >> 305 >> 306 // >> 307 // decide verticies of -ve phi boundary >> 308 // >> 309 >> 310 G4double cosPhi = std::cos(phimin); >> 311 G4double sinPhi = std::sin(phimin); >> 312 >> 313 // Outer hyperbolic surface >> 314 >> 315 v0 = transform.TransformPoint( >> 316 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, >> 317 + fZHalfLength)); >> 318 v1 = transform.TransformPoint( >> 319 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, >> 320 - fZHalfLength)); >> 321 if (splitOuter) >> 322 { >> 323 v4 = transform.TransformPoint( >> 324 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, 0)); >> 325 } >> 326 >> 327 // Inner hyperbolic surface 313 328 314 // Find bounding box << 329 G4double zInnerSplit = 0.; 315 if (dphi <= 0 || totalphi >= CLHEP::twopi) << 330 if (splitInner) 316 { 331 { 317 pMin.set(-rmax,-rmax, zmin); << 332 v2 = transform.TransformPoint( 318 pMax.set( rmax, rmax, zmax); << 333 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, >> 334 + fZHalfLength)); >> 335 v3 = transform.TransformPoint( >> 336 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, >> 337 - fZHalfLength)); >> 338 >> 339 // Find intersection of tangential line of inner >> 340 // surface at z = fZHalfLength and line r=fInnerRadius. >> 341 G4double dr = fZHalfLength * fTanInnerStereo2; >> 342 G4double dz = maxEndInnerRad; >> 343 zInnerSplit = fZHalfLength + (fInnerRadius - maxEndInnerRad) * dz / dr; >> 344 >> 345 // Build associated vertices >> 346 v5 = transform.TransformPoint( >> 347 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, >> 348 + zInnerSplit)); >> 349 v6 = transform.TransformPoint( >> 350 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, >> 351 - zInnerSplit)); 319 } 352 } 320 else 353 else 321 { 354 { 322 G4TwoVector vmin,vmax; << 355 v2 = transform.TransformPoint( 323 G4GeomTools::DiskExtent(rmin, rmax, sphi, << 356 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, 324 pMin.set(vmin.x(), vmin.y(), zmin); << 357 + fZHalfLength)); 325 pMax.set(vmax.x(), vmax.y(), zmax); << 358 v3 = transform.TransformPoint( >> 359 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, >> 360 - fZHalfLength)); 326 } 361 } 327 362 328 // Check correctness of the bounding box << 329 // 363 // 330 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 364 // decide vertices of +ve phi boundary >> 365 // >> 366 >> 367 cosPhi = std::cos(phimax); >> 368 sinPhi = std::sin(phimax); >> 369 >> 370 // Outer hyperbolic surface >> 371 >> 372 w0 = transform.TransformPoint( >> 373 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, >> 374 + fZHalfLength)); >> 375 w1 = transform.TransformPoint( >> 376 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, >> 377 - fZHalfLength)); >> 378 if (splitOuter) >> 379 { >> 380 G4double r = rFudge*fOuterRadius; >> 381 >> 382 w4 = transform.TransformPoint(G4ThreeVector( r*cosPhi, r*sinPhi, 0 )); >> 383 >> 384 AddPolyToExtent( v0, v4, w4, w0, voxelLimit, axis, extentList ); >> 385 AddPolyToExtent( v4, v1, w1, w4, voxelLimit, axis, extentList ); >> 386 } >> 387 else >> 388 { >> 389 AddPolyToExtent( v0, v1, w1, w0, voxelLimit, axis, extentList ); >> 390 } >> 391 >> 392 // Inner hyperbolic surface >> 393 >> 394 if (splitInner) 331 { 395 { 332 std::ostringstream message; << 396 w2 = transform.TransformPoint( 333 message << "Bad bounding box (min >= max) << 397 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, 334 << GetName() << " !" << 398 + fZHalfLength)); 335 << "\npMin = " << pMin << 399 w3 = transform.TransformPoint( 336 << "\npMax = " << pMax; << 400 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi, 337 G4Exception("G4TwistedTubs::BoundingLimits << 401 - fZHalfLength)); 338 JustWarning, message); << 402 339 DumpInfo(); << 403 w5 = transform.TransformPoint( >> 404 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, >> 405 + zInnerSplit)); >> 406 w6 = transform.TransformPoint( >> 407 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, >> 408 - zInnerSplit)); >> 409 >> 410 AddPolyToExtent( v3, v6, w6, w3, voxelLimit, axis, extentList ); >> 411 AddPolyToExtent( v6, v5, w5, w6, voxelLimit, axis, extentList ); >> 412 AddPolyToExtent( v5, v2, w2, w5, voxelLimit, axis, extentList ); >> 413 340 } 414 } >> 415 else >> 416 { >> 417 w2 = transform.TransformPoint( >> 418 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, >> 419 + fZHalfLength)); >> 420 w3 = transform.TransformPoint( >> 421 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi, >> 422 - fZHalfLength)); >> 423 >> 424 AddPolyToExtent( v3, v2, w2, w3, voxelLimit, axis, extentList ); >> 425 } >> 426 >> 427 // >> 428 // Endplate segments >> 429 // >> 430 AddPolyToExtent( v1, v3, w3, w1, voxelLimit, axis, extentList ); >> 431 AddPolyToExtent( v2, v0, w0, w2, voxelLimit, axis, extentList ); >> 432 >> 433 // >> 434 // Return min/max value >> 435 // >> 436 return extentList.GetExtent( min, max ); 341 } 437 } 342 438 >> 439 343 //============================================ 440 //===================================================================== 344 //* CalculateExtent -------------------------- << 441 //* AddPolyToExtent --------------------------------------------------- 345 442 346 G4bool << 443 void G4TwistedTubs::AddPolyToExtent( const G4ThreeVector &v0, 347 G4TwistedTubs::CalculateExtent(const EAxis pAx << 444 const G4ThreeVector &v1, 348 const G4VoxelLi << 445 const G4ThreeVector &w1, 349 const G4AffineT << 446 const G4ThreeVector &w0, 350 G4double& << 447 const G4VoxelLimits &voxelLimit, 351 { << 448 const EAxis axis, 352 G4ThreeVector bmin, bmax; << 449 G4SolidExtentList &extentList ) 353 << 450 { 354 // Get bounding box << 451 // Utility function for CalculateExtent 355 BoundingLimits(bmin,bmax); << 452 // 356 << 453 G4ClippablePolygon phiPoly; 357 // Find extent << 454 358 G4BoundingEnvelope bbox(bmin,bmax); << 455 phiPoly.AddVertexInOrder( v0 ); 359 return bbox.CalculateExtent(pAxis,pVoxelLimi << 456 phiPoly.AddVertexInOrder( v1 ); >> 457 phiPoly.AddVertexInOrder( w1 ); >> 458 phiPoly.AddVertexInOrder( w0 ); >> 459 >> 460 if (phiPoly.PartialClip( voxelLimit, axis )) >> 461 { >> 462 phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() ); >> 463 extentList.AddSurface( phiPoly ); >> 464 } 360 } 465 } 361 466 362 467 363 //============================================ 468 //===================================================================== 364 //* Inside ----------------------------------- 469 //* Inside ------------------------------------------------------------ 365 470 366 EInside G4TwistedTubs::Inside(const G4ThreeVec 471 EInside G4TwistedTubs::Inside(const G4ThreeVector& p) const 367 { 472 { 368 473 369 const G4double halftol << 474 static const G4double halftol 370 = 0.5 * G4GeometryTolerance::GetInstance( 475 = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance(); 371 // static G4int timerid = -1; 476 // static G4int timerid = -1; 372 // G4Timer timer(timerid, "G4TwistedTubs", 477 // G4Timer timer(timerid, "G4TwistedTubs", "Inside"); 373 // timer.Start(); 478 // timer.Start(); >> 479 >> 480 374 481 >> 482 G4ThreeVector *tmpp; >> 483 EInside *tmpinside; >> 484 if (fLastInside.p == p) >> 485 { >> 486 return fLastInside.inside; >> 487 } >> 488 else >> 489 { >> 490 tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p)); >> 491 tmpinside = const_cast<EInside*>(&(fLastInside.inside)); >> 492 tmpp->set(p.x(), p.y(), p.z()); >> 493 } 375 494 376 EInside outerhypearea = ((G4TwistTubsHypeS 495 EInside outerhypearea = ((G4TwistTubsHypeSide *)fOuterHype)->Inside(p); 377 G4double innerhyperho = ((G4TwistTubsHypeS 496 G4double innerhyperho = ((G4TwistTubsHypeSide *)fInnerHype)->GetRhoAtPZ(p); 378 G4double distanceToOut = p.getRho() - inner 497 G4double distanceToOut = p.getRho() - innerhyperho; // +ve: inside 379 EInside tmpinside; << 498 380 if ((outerhypearea == kOutside) || (distanc 499 if ((outerhypearea == kOutside) || (distanceToOut < -halftol)) 381 { 500 { 382 tmpinside = kOutside; << 501 *tmpinside = kOutside; 383 } 502 } 384 else if (outerhypearea == kSurface) 503 else if (outerhypearea == kSurface) 385 { 504 { 386 tmpinside = kSurface; << 505 *tmpinside = kSurface; 387 } 506 } 388 else 507 else 389 { 508 { 390 if (distanceToOut <= halftol) 509 if (distanceToOut <= halftol) 391 { 510 { 392 tmpinside = kSurface; << 511 *tmpinside = kSurface; 393 } 512 } 394 else 513 else 395 { 514 { 396 tmpinside = kInside; << 515 *tmpinside = kInside; 397 } 516 } 398 } 517 } 399 518 400 return tmpinside; << 519 return fLastInside.inside; 401 } 520 } 402 521 403 //============================================ 522 //===================================================================== 404 //* SurfaceNormal ---------------------------- 523 //* SurfaceNormal ----------------------------------------------------- 405 524 406 G4ThreeVector G4TwistedTubs::SurfaceNormal(con 525 G4ThreeVector G4TwistedTubs::SurfaceNormal(const G4ThreeVector& p) const 407 { 526 { 408 // 527 // 409 // return the normal unit vector to the Hyp 528 // return the normal unit vector to the Hyperbolical Surface at a point 410 // p on (or nearly on) the surface 529 // p on (or nearly on) the surface 411 // 530 // 412 // Which of the three or four surfaces are 531 // Which of the three or four surfaces are we closest to? 413 // 532 // 414 533 >> 534 if (fLastNormal.p == p) >> 535 { >> 536 return fLastNormal.vec; >> 537 } >> 538 G4ThreeVector *tmpp = >> 539 const_cast<G4ThreeVector*>(&(fLastNormal.p)); >> 540 G4ThreeVector *tmpnormal = >> 541 const_cast<G4ThreeVector*>(&(fLastNormal.vec)); >> 542 G4VTwistSurface **tmpsurface = >> 543 const_cast<G4VTwistSurface**>(fLastNormal.surface); >> 544 tmpp->set(p.x(), p.y(), p.z()); 415 545 416 G4double distance = kInfinity; 546 G4double distance = kInfinity; 417 547 418 G4VTwistSurface *surfaces[6]; 548 G4VTwistSurface *surfaces[6]; 419 surfaces[0] = fLatterTwisted; 549 surfaces[0] = fLatterTwisted; 420 surfaces[1] = fFormerTwisted; 550 surfaces[1] = fFormerTwisted; 421 surfaces[2] = fInnerHype; 551 surfaces[2] = fInnerHype; 422 surfaces[3] = fOuterHype; 552 surfaces[3] = fOuterHype; 423 surfaces[4] = fLowerEndcap; 553 surfaces[4] = fLowerEndcap; 424 surfaces[5] = fUpperEndcap; 554 surfaces[5] = fUpperEndcap; 425 555 426 G4ThreeVector xx; 556 G4ThreeVector xx; 427 G4ThreeVector bestxx; 557 G4ThreeVector bestxx; >> 558 G4int i; 428 G4int besti = -1; 559 G4int besti = -1; 429 for (auto i=0; i<6; ++i) << 560 for (i=0; i< 6; i++) 430 { 561 { 431 G4double tmpdistance = surfaces[i]->Dist 562 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx); 432 if (tmpdistance < distance) 563 if (tmpdistance < distance) 433 { 564 { 434 distance = tmpdistance; 565 distance = tmpdistance; 435 bestxx = xx; 566 bestxx = xx; 436 besti = i; 567 besti = i; 437 } 568 } 438 } 569 } 439 570 440 return surfaces[besti]->GetNormal(bestxx, tr << 571 tmpsurface[0] = surfaces[besti]; >> 572 *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true); >> 573 >> 574 return fLastNormal.vec; 441 } 575 } 442 576 443 //============================================ 577 //===================================================================== 444 //* DistanceToIn (p, v) ---------------------- 578 //* DistanceToIn (p, v) ----------------------------------------------- 445 579 446 G4double G4TwistedTubs::DistanceToIn (const G4 580 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p, 447 const G4 581 const G4ThreeVector& v ) const 448 { 582 { 449 583 450 // DistanceToIn (p, v): 584 // DistanceToIn (p, v): 451 // Calculate distance to surface of shape f 585 // Calculate distance to surface of shape from `outside' 452 // along with the v, allowing for tolerance 586 // along with the v, allowing for tolerance. 453 // The function returns kInfinity if no int 587 // The function returns kInfinity if no intersection or 454 // just grazing within tolerance. 588 // just grazing within tolerance. 455 589 456 // 590 // >> 591 // checking last value >> 592 // >> 593 >> 594 G4ThreeVector *tmpp; >> 595 G4ThreeVector *tmpv; >> 596 G4double *tmpdist; >> 597 if ((fLastDistanceToInWithV.p == p) && (fLastDistanceToInWithV.vec == v)) >> 598 { >> 599 return fLastDistanceToIn.value; >> 600 } >> 601 else >> 602 { >> 603 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p)); >> 604 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec)); >> 605 tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value)); >> 606 tmpp->set(p.x(), p.y(), p.z()); >> 607 tmpv->set(v.x(), v.y(), v.z()); >> 608 } >> 609 >> 610 // 457 // Calculate DistanceToIn(p,v) 611 // Calculate DistanceToIn(p,v) 458 // 612 // 459 613 460 EInside currentside = Inside(p); 614 EInside currentside = Inside(p); 461 615 462 if (currentside == kInside) 616 if (currentside == kInside) 463 { 617 { 464 } 618 } 465 else 619 else 466 { 620 { 467 if (currentside == kSurface) 621 if (currentside == kSurface) 468 { 622 { 469 // particle is just on a boundary. 623 // particle is just on a boundary. 470 // If the particle is entering to the v 624 // If the particle is entering to the volume, return 0. 471 // 625 // 472 G4ThreeVector normal = SurfaceNormal(p) 626 G4ThreeVector normal = SurfaceNormal(p); 473 if (normal*v < 0) 627 if (normal*v < 0) 474 { 628 { 475 return 0; << 629 *tmpdist = 0; >> 630 return fLastDistanceToInWithV.value; 476 } 631 } 477 } 632 } 478 } 633 } 479 634 480 // now, we can take smallest positive dista 635 // now, we can take smallest positive distance. 481 636 482 // Initialize 637 // Initialize 483 // 638 // 484 G4double distance = kInfinity; 639 G4double distance = kInfinity; 485 640 486 // find intersections and choose nearest on 641 // find intersections and choose nearest one. 487 // 642 // 488 G4VTwistSurface* surfaces[6]; << 643 G4VTwistSurface *surfaces[6]; 489 surfaces[0] = fLowerEndcap; 644 surfaces[0] = fLowerEndcap; 490 surfaces[1] = fUpperEndcap; 645 surfaces[1] = fUpperEndcap; 491 surfaces[2] = fLatterTwisted; 646 surfaces[2] = fLatterTwisted; 492 surfaces[3] = fFormerTwisted; 647 surfaces[3] = fFormerTwisted; 493 surfaces[4] = fInnerHype; 648 surfaces[4] = fInnerHype; 494 surfaces[5] = fOuterHype; 649 surfaces[5] = fOuterHype; 495 650 496 G4ThreeVector xx; 651 G4ThreeVector xx; 497 G4ThreeVector bestxx; 652 G4ThreeVector bestxx; 498 for (const auto & surface : surfaces) << 653 G4int i; >> 654 G4int besti = -1; >> 655 for (i=0; i< 6; i++) 499 { 656 { 500 G4double tmpdistance = surface->Distance << 657 G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx); 501 if (tmpdistance < distance) 658 if (tmpdistance < distance) 502 { 659 { 503 distance = tmpdistance; 660 distance = tmpdistance; 504 bestxx = xx; 661 bestxx = xx; >> 662 besti = i; 505 } 663 } 506 } 664 } 507 return distance; << 665 *tmpdist = distance; >> 666 >> 667 return fLastDistanceToInWithV.value; 508 } 668 } 509 669 510 //============================================ 670 //===================================================================== 511 //* DistanceToIn (p) ------------------------- 671 //* DistanceToIn (p) -------------------------------------------------- 512 672 513 G4double G4TwistedTubs::DistanceToIn (const G4 673 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p) const 514 { 674 { 515 // DistanceToIn(p): 675 // DistanceToIn(p): 516 // Calculate distance to surface of shape f 676 // Calculate distance to surface of shape from `outside', 517 // allowing for tolerance 677 // allowing for tolerance >> 678 >> 679 // >> 680 // checking last value >> 681 // >> 682 >> 683 G4ThreeVector *tmpp; >> 684 G4double *tmpdist; >> 685 if (fLastDistanceToIn.p == p) >> 686 { >> 687 return fLastDistanceToIn.value; >> 688 } >> 689 else >> 690 { >> 691 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p)); >> 692 tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value)); >> 693 tmpp->set(p.x(), p.y(), p.z()); >> 694 } 518 695 519 // 696 // 520 // Calculate DistanceToIn(p) 697 // Calculate DistanceToIn(p) 521 // 698 // 522 699 523 EInside currentside = Inside(p); 700 EInside currentside = Inside(p); 524 701 525 switch (currentside) 702 switch (currentside) 526 { 703 { 527 case (kInside) : 704 case (kInside) : 528 {} 705 {} 529 case (kSurface) : 706 case (kSurface) : 530 { 707 { 531 return 0; << 708 *tmpdist = 0.; >> 709 return fLastDistanceToIn.value; 532 } 710 } 533 case (kOutside) : 711 case (kOutside) : 534 { 712 { 535 // Initialize 713 // Initialize 536 G4double distance = kInfinity; << 714 G4double distance = kInfinity; 537 715 538 // find intersections and choose near 716 // find intersections and choose nearest one. 539 G4VTwistSurface *surfaces[6]; 717 G4VTwistSurface *surfaces[6]; 540 surfaces[0] = fLowerEndcap; 718 surfaces[0] = fLowerEndcap; 541 surfaces[1] = fUpperEndcap; 719 surfaces[1] = fUpperEndcap; 542 surfaces[2] = fLatterTwisted; 720 surfaces[2] = fLatterTwisted; 543 surfaces[3] = fFormerTwisted; 721 surfaces[3] = fFormerTwisted; 544 surfaces[4] = fInnerHype; 722 surfaces[4] = fInnerHype; 545 surfaces[5] = fOuterHype; 723 surfaces[5] = fOuterHype; 546 724 >> 725 G4int i; >> 726 G4int besti = -1; 547 G4ThreeVector xx; 727 G4ThreeVector xx; 548 G4ThreeVector bestxx; 728 G4ThreeVector bestxx; 549 for (const auto & surface : surfaces) << 729 for (i=0; i< 6; i++) 550 { 730 { 551 G4double tmpdistance = surface->Di << 731 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx); 552 if (tmpdistance < distance) 732 if (tmpdistance < distance) 553 { 733 { 554 distance = tmpdistance; 734 distance = tmpdistance; 555 bestxx = xx; 735 bestxx = xx; >> 736 besti = i; 556 } 737 } 557 } 738 } 558 return distance; << 739 *tmpdist = distance; >> 740 return fLastDistanceToIn.value; 559 } 741 } 560 default : 742 default : 561 { 743 { 562 G4Exception("G4TwistedTubs::DistanceT << 744 G4Exception("G4TwistedTubs::DistanceToIn(p)", "InvalidCondition", 563 FatalException, "Unknown 745 FatalException, "Unknown point location!"); 564 } 746 } 565 } // switch end 747 } // switch end 566 748 567 return kInfinity; 749 return kInfinity; 568 } 750 } 569 751 570 //============================================ 752 //===================================================================== 571 //* DistanceToOut (p, v) --------------------- 753 //* DistanceToOut (p, v) ---------------------------------------------- 572 754 573 G4double G4TwistedTubs::DistanceToOut( const G 755 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p, 574 const G 756 const G4ThreeVector& v, 575 const G 757 const G4bool calcNorm, 576 G4bool 758 G4bool *validNorm, 577 G4Three 759 G4ThreeVector *norm ) const 578 { 760 { 579 // DistanceToOut (p, v): 761 // DistanceToOut (p, v): 580 // Calculate distance to surface of shape f 762 // Calculate distance to surface of shape from `inside' 581 // along with the v, allowing for tolerance 763 // along with the v, allowing for tolerance. 582 // The function returns kInfinity if no int 764 // The function returns kInfinity if no intersection or 583 // just grazing within tolerance. 765 // just grazing within tolerance. 584 766 585 // 767 // >> 768 // checking last value >> 769 // >> 770 >> 771 G4ThreeVector *tmpp; >> 772 G4ThreeVector *tmpv; >> 773 G4double *tmpdist; >> 774 if ((fLastDistanceToOutWithV.p == p) && (fLastDistanceToOutWithV.vec == v) ) >> 775 { >> 776 return fLastDistanceToOutWithV.value; >> 777 } >> 778 else >> 779 { >> 780 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p)); >> 781 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec)); >> 782 tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value)); >> 783 tmpp->set(p.x(), p.y(), p.z()); >> 784 tmpv->set(v.x(), v.y(), v.z()); >> 785 } >> 786 >> 787 // 586 // Calculate DistanceToOut(p,v) 788 // Calculate DistanceToOut(p,v) 587 // 789 // 588 790 589 EInside currentside = Inside(p); 791 EInside currentside = Inside(p); >> 792 590 if (currentside == kOutside) 793 if (currentside == kOutside) 591 { 794 { 592 } 795 } 593 else 796 else 594 { 797 { 595 if (currentside == kSurface) 798 if (currentside == kSurface) 596 { 799 { 597 // particle is just on a boundary. 800 // particle is just on a boundary. 598 // If the particle is exiting from the 801 // If the particle is exiting from the volume, return 0. 599 // 802 // 600 G4ThreeVector normal = SurfaceNormal(p) 803 G4ThreeVector normal = SurfaceNormal(p); >> 804 G4VTwistSurface *blockedsurface = fLastNormal.surface[0]; 601 if (normal*v > 0) 805 if (normal*v > 0) 602 { 806 { 603 if (calcNorm) 807 if (calcNorm) 604 { 808 { 605 *norm = normal; << 809 *norm = (blockedsurface->GetNormal(p, true)); 606 *validNorm = true; << 810 *validNorm = blockedsurface->IsValidNorm(); 607 } 811 } 608 return 0; << 812 *tmpdist = 0.; >> 813 return fLastDistanceToOutWithV.value; 609 } 814 } 610 } 815 } 611 } 816 } 612 817 613 // now, we can take smallest positive dista 818 // now, we can take smallest positive distance. 614 819 615 // Initialize 820 // Initialize 616 // 821 // 617 G4double distance = kInfinity; << 822 G4double distance = kInfinity; 618 823 619 // find intersections and choose nearest on 824 // find intersections and choose nearest one. 620 // 825 // 621 G4VTwistSurface* surfaces[6]; << 826 G4VTwistSurface *surfaces[6]; 622 surfaces[0] = fLatterTwisted; 827 surfaces[0] = fLatterTwisted; 623 surfaces[1] = fFormerTwisted; 828 surfaces[1] = fFormerTwisted; 624 surfaces[2] = fInnerHype; 829 surfaces[2] = fInnerHype; 625 surfaces[3] = fOuterHype; 830 surfaces[3] = fOuterHype; 626 surfaces[4] = fLowerEndcap; 831 surfaces[4] = fLowerEndcap; 627 surfaces[5] = fUpperEndcap; 832 surfaces[5] = fUpperEndcap; 628 833 >> 834 G4int i; 629 G4int besti = -1; 835 G4int besti = -1; 630 G4ThreeVector xx; 836 G4ThreeVector xx; 631 G4ThreeVector bestxx; 837 G4ThreeVector bestxx; 632 for (auto i=0; i<6; ++i) << 838 for (i=0; i< 6; i++) 633 { 839 { 634 G4double tmpdistance = surfaces[i]->Dist 840 G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx); 635 if (tmpdistance < distance) 841 if (tmpdistance < distance) 636 { 842 { 637 distance = tmpdistance; 843 distance = tmpdistance; 638 bestxx = xx; 844 bestxx = xx; 639 besti = i; 845 besti = i; 640 } 846 } 641 } 847 } 642 848 643 if (calcNorm) 849 if (calcNorm) 644 { 850 { 645 if (besti != -1) 851 if (besti != -1) 646 { 852 { 647 *norm = (surfaces[besti]->GetNormal(p 853 *norm = (surfaces[besti]->GetNormal(p, true)); 648 *validNorm = surfaces[besti]->IsValid 854 *validNorm = surfaces[besti]->IsValidNorm(); 649 } 855 } 650 } 856 } 651 857 652 return distance; << 858 *tmpdist = distance; >> 859 >> 860 return fLastDistanceToOutWithV.value; 653 } 861 } 654 862 655 863 656 //============================================ 864 //===================================================================== 657 //* DistanceToOut (p) ------------------------ 865 //* DistanceToOut (p) ---------------------------------------------- 658 866 659 G4double G4TwistedTubs::DistanceToOut( const G 867 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p ) const 660 { 868 { 661 // DistanceToOut(p): 869 // DistanceToOut(p): 662 // Calculate distance to surface of shape f 870 // Calculate distance to surface of shape from `inside', 663 // allowing for tolerance 871 // allowing for tolerance 664 872 665 // 873 // >> 874 // checking last value >> 875 // >> 876 >> 877 G4ThreeVector *tmpp; >> 878 G4double *tmpdist; >> 879 if (fLastDistanceToOut.p == p) >> 880 { >> 881 return fLastDistanceToOut.value; >> 882 } >> 883 else >> 884 { >> 885 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p)); >> 886 tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value)); >> 887 tmpp->set(p.x(), p.y(), p.z()); >> 888 } >> 889 >> 890 // 666 // Calculate DistanceToOut(p) 891 // Calculate DistanceToOut(p) 667 // 892 // 668 893 669 EInside currentside = Inside(p); 894 EInside currentside = Inside(p); 670 895 671 switch (currentside) 896 switch (currentside) 672 { 897 { 673 case (kOutside) : 898 case (kOutside) : 674 { 899 { 675 } 900 } 676 case (kSurface) : 901 case (kSurface) : 677 { 902 { 678 return 0; << 903 *tmpdist = 0.; >> 904 return fLastDistanceToOut.value; 679 } 905 } 680 case (kInside) : 906 case (kInside) : 681 { 907 { 682 // Initialize 908 // Initialize 683 G4double distance = kInfinity; 909 G4double distance = kInfinity; 684 910 685 // find intersections and choose near 911 // find intersections and choose nearest one. 686 G4VTwistSurface* surfaces[6]; << 912 G4VTwistSurface *surfaces[6]; 687 surfaces[0] = fLatterTwisted; 913 surfaces[0] = fLatterTwisted; 688 surfaces[1] = fFormerTwisted; 914 surfaces[1] = fFormerTwisted; 689 surfaces[2] = fInnerHype; 915 surfaces[2] = fInnerHype; 690 surfaces[3] = fOuterHype; 916 surfaces[3] = fOuterHype; 691 surfaces[4] = fLowerEndcap; 917 surfaces[4] = fLowerEndcap; 692 surfaces[5] = fUpperEndcap; 918 surfaces[5] = fUpperEndcap; 693 919 >> 920 G4int i; >> 921 G4int besti = -1; 694 G4ThreeVector xx; 922 G4ThreeVector xx; 695 G4ThreeVector bestxx; 923 G4ThreeVector bestxx; 696 for (const auto & surface : surfaces) << 924 for (i=0; i< 6; i++) 697 { 925 { 698 G4double tmpdistance = surface->Di << 926 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx); 699 if (tmpdistance < distance) 927 if (tmpdistance < distance) 700 { 928 { 701 distance = tmpdistance; 929 distance = tmpdistance; 702 bestxx = xx; 930 bestxx = xx; >> 931 besti = i; 703 } 932 } 704 } 933 } 705 return distance; << 934 *tmpdist = distance; >> 935 >> 936 return fLastDistanceToOut.value; 706 } 937 } 707 default : 938 default : 708 { 939 { 709 G4Exception("G4TwistedTubs::DistanceT << 940 G4Exception("G4TwistedTubs::DistanceToOut(p)", "InvalidCondition", 710 FatalException, "Unknown 941 FatalException, "Unknown point location!"); 711 } 942 } 712 } // switch end 943 } // switch end 713 944 714 return 0.; << 945 return 0; 715 } 946 } 716 947 717 //============================================ 948 //===================================================================== 718 //* StreamInfo ------------------------------- 949 //* StreamInfo -------------------------------------------------------- 719 950 720 std::ostream& G4TwistedTubs::StreamInfo(std::o 951 std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const 721 { 952 { 722 // 953 // 723 // Stream object contents to an output strea 954 // Stream object contents to an output stream 724 // 955 // 725 G4long oldprc = os.precision(16); << 726 os << "------------------------------------- 956 os << "-----------------------------------------------------------\n" 727 << " *** Dump for solid - " << GetName 957 << " *** Dump for solid - " << GetName() << " ***\n" 728 << " ================================= 958 << " ===================================================\n" 729 << " Solid type: G4TwistedTubs\n" 959 << " Solid type: G4TwistedTubs\n" 730 << " Parameters: \n" 960 << " Parameters: \n" 731 << " -ve end Z : " << fEn 961 << " -ve end Z : " << fEndZ[0]/mm << " mm \n" 732 << " +ve end Z : " << fEn 962 << " +ve end Z : " << fEndZ[1]/mm << " mm \n" 733 << " inner end radius(-ve z): " << fEn 963 << " inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n" 734 << " inner end radius(+ve z): " << fEn 964 << " inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n" 735 << " outer end radius(-ve z): " << fEn 965 << " outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n" 736 << " outer end radius(+ve z): " << fEn 966 << " outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n" 737 << " inner radius (z=0) : " << fIn 967 << " inner radius (z=0) : " << fInnerRadius/mm << " mm \n" 738 << " outer radius (z=0) : " << fOu 968 << " outer radius (z=0) : " << fOuterRadius/mm << " mm \n" 739 << " twisted angle : " << fPh 969 << " twisted angle : " << fPhiTwist/degree << " degrees \n" 740 << " inner stereo angle : " << fIn 970 << " inner stereo angle : " << fInnerStereo/degree << " degrees \n" 741 << " outer stereo angle : " << fOu 971 << " outer stereo angle : " << fOuterStereo/degree << " degrees \n" 742 << " phi-width of a piece : " << fDP 972 << " phi-width of a piece : " << fDPhi/degree << " degrees \n" 743 << "------------------------------------- 973 << "-----------------------------------------------------------\n"; 744 os.precision(oldprc); << 745 974 746 return os; 975 return os; 747 } 976 } 748 977 749 978 750 //============================================ 979 //===================================================================== 751 //* DiscribeYourselfTo ----------------------- 980 //* DiscribeYourselfTo ------------------------------------------------ 752 981 753 void G4TwistedTubs::DescribeYourselfTo (G4VGra 982 void G4TwistedTubs::DescribeYourselfTo (G4VGraphicsScene& scene) const 754 { 983 { 755 scene.AddSolid (*this); 984 scene.AddSolid (*this); 756 } 985 } 757 986 758 //============================================ 987 //===================================================================== 759 //* GetExtent -------------------------------- 988 //* GetExtent --------------------------------------------------------- 760 989 761 G4VisExtent G4TwistedTubs::GetExtent() const << 990 G4VisExtent G4TwistedTubs::GetExtent() const 762 { 991 { 763 // Define the sides of the box into which th 992 // Define the sides of the box into which the G4Tubs instance would fit. 764 // << 993 765 G4ThreeVector pmin,pmax; << 994 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1); 766 BoundingLimits(pmin,pmax); << 995 return G4VisExtent( -maxEndOuterRad, maxEndOuterRad, 767 return { pmin.x(),pmax.x(), << 996 -maxEndOuterRad, maxEndOuterRad, 768 pmin.y(),pmax.y(), << 997 -fZHalfLength, fZHalfLength ); 769 pmin.z(),pmax.z() }; << 770 } 998 } 771 999 772 //============================================ 1000 //===================================================================== 773 //* CreatePolyhedron ------------------------- 1001 //* CreatePolyhedron -------------------------------------------------- 774 1002 775 G4Polyhedron* G4TwistedTubs::CreatePolyhedron 1003 G4Polyhedron* G4TwistedTubs::CreatePolyhedron () const 776 { 1004 { 777 // number of meshes 1005 // number of meshes 778 // 1006 // 779 G4double absPhiTwist = std::abs(fPhiTwist); << 1007 G4double dA = std::max(fDPhi,fPhiTwist); 780 G4double dA = std::max(fDPhi,absPhiTwist); << 1008 const G4int m = 781 const G4int k = << 782 G4int(G4Polyhedron::GetNumberOfRotationSte 1009 G4int(G4Polyhedron::GetNumberOfRotationSteps() * dA / twopi) + 2; 783 const G4int n = 1010 const G4int n = 784 G4int(G4Polyhedron::GetNumberOfRotationSte << 1011 G4int(G4Polyhedron::GetNumberOfRotationSteps() * fPhiTwist / twopi) + 2; 785 1012 786 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ; << 1013 const G4int nnodes = 4*(m-1)*(n-2) + 2*m*m ; 787 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1) << 1014 const G4int nfaces = 4*(m-1)*(n-1) + 2*(m-1)*(m-1) ; 788 1015 789 auto ph = new G4Polyhedron; << 1016 G4Polyhedron *ph=new G4Polyhedron; 790 typedef G4double G4double3[3]; 1017 typedef G4double G4double3[3]; 791 typedef G4int G4int4[4]; 1018 typedef G4int G4int4[4]; 792 auto xyz = new G4double3[nnodes]; // number << 1019 G4double3* xyz = new G4double3[nnodes]; // number of nodes 793 auto faces = new G4int4[nfaces] ; // number << 1020 G4int4* faces = new G4int4[nfaces] ; // number of faces 794 fLowerEndcap->GetFacets(k,k,xyz,faces,0) ; << 1021 fLowerEndcap->GetFacets(m,m,xyz,faces,0) ; 795 fUpperEndcap->GetFacets(k,k,xyz,faces,1) ; << 1022 fUpperEndcap->GetFacets(m,m,xyz,faces,1) ; 796 fInnerHype->GetFacets(k,n,xyz,faces,2) ; << 1023 fInnerHype->GetFacets(m,n,xyz,faces,2) ; 797 fFormerTwisted->GetFacets(k,n,xyz,faces,3) ; << 1024 fFormerTwisted->GetFacets(m,n,xyz,faces,3) ; 798 fOuterHype->GetFacets(k,n,xyz,faces,4) ; << 1025 fOuterHype->GetFacets(m,n,xyz,faces,4) ; 799 fLatterTwisted->GetFacets(k,n,xyz,faces,5) ; << 1026 fLatterTwisted->GetFacets(m,n,xyz,faces,5) ; 800 1027 801 ph->createPolyhedron(nnodes,nfaces,xyz,faces 1028 ph->createPolyhedron(nnodes,nfaces,xyz,faces); 802 1029 803 delete[] xyz; 1030 delete[] xyz; 804 delete[] faces; 1031 delete[] faces; 805 1032 806 return ph; 1033 return ph; 807 } 1034 } 808 1035 809 //============================================ 1036 //===================================================================== >> 1037 //* CreateNUBS -------------------------------------------------------- >> 1038 >> 1039 G4NURBS* G4TwistedTubs::CreateNURBS () const >> 1040 { >> 1041 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1); >> 1042 G4double maxEndInnerRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1); >> 1043 return new G4NURBStube(maxEndInnerRad, maxEndOuterRad, fZHalfLength); >> 1044 // Tube for now!!! >> 1045 } >> 1046 >> 1047 //===================================================================== 810 //* GetPolyhedron ---------------------------- 1048 //* GetPolyhedron ----------------------------------------------------- 811 1049 812 G4Polyhedron* G4TwistedTubs::GetPolyhedron () 1050 G4Polyhedron* G4TwistedTubs::GetPolyhedron () const 813 { 1051 { 814 if (fpPolyhedron == nullptr || << 1052 if ((!fpPolyhedron) || 815 fRebuildPolyhedron || << 1053 (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != 816 fpPolyhedron->GetNumberOfRotationStepsAt << 1054 fpPolyhedron->GetNumberOfRotationSteps())) 817 fpPolyhedron->GetNumberOfRotationSteps() << 818 { 1055 { 819 G4AutoLock l(&polyhedronMutex); << 820 delete fpPolyhedron; 1056 delete fpPolyhedron; 821 fpPolyhedron = CreatePolyhedron(); 1057 fpPolyhedron = CreatePolyhedron(); 822 fRebuildPolyhedron = false; << 823 l.unlock(); << 824 } 1058 } 825 return fpPolyhedron; 1059 return fpPolyhedron; 826 } 1060 } 827 1061 828 //============================================ 1062 //===================================================================== 829 //* CreateSurfaces --------------------------- 1063 //* CreateSurfaces ---------------------------------------------------- 830 1064 831 void G4TwistedTubs::CreateSurfaces() 1065 void G4TwistedTubs::CreateSurfaces() 832 { 1066 { 833 // create 6 surfaces of TwistedTub 1067 // create 6 surfaces of TwistedTub 834 1068 835 G4ThreeVector x0(0, 0, fEndZ[0]); 1069 G4ThreeVector x0(0, 0, fEndZ[0]); 836 G4ThreeVector n (0, 0, -1); 1070 G4ThreeVector n (0, 0, -1); 837 1071 838 fLowerEndcap = new G4TwistTubsFlatSide("Low 1072 fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap", 839 fEndInnerR 1073 fEndInnerRadius, fEndOuterRadius, 840 fDPhi, fEn 1074 fDPhi, fEndPhi, fEndZ, -1) ; 841 1075 842 fUpperEndcap = new G4TwistTubsFlatSide("Upp 1076 fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap", 843 fEndInnerR 1077 fEndInnerRadius, fEndOuterRadius, 844 fDPhi, fEn 1078 fDPhi, fEndPhi, fEndZ, 1) ; 845 1079 846 G4RotationMatrix rotHalfDPhi; 1080 G4RotationMatrix rotHalfDPhi; 847 rotHalfDPhi.rotateZ(0.5*fDPhi); 1081 rotHalfDPhi.rotateZ(0.5*fDPhi); 848 1082 849 fLatterTwisted = new G4TwistTubsSide("Latte 1083 fLatterTwisted = new G4TwistTubsSide("LatterTwisted", 850 fEndI 1084 fEndInnerRadius, fEndOuterRadius, 851 fDPhi 1085 fDPhi, fEndPhi, fEndZ, 852 fInne 1086 fInnerRadius, fOuterRadius, fKappa, 853 1 ) ; 1087 1 ) ; 854 fFormerTwisted = new G4TwistTubsSide("Forme 1088 fFormerTwisted = new G4TwistTubsSide("FormerTwisted", 855 fEndI 1089 fEndInnerRadius, fEndOuterRadius, 856 fDPhi 1090 fDPhi, fEndPhi, fEndZ, 857 fInne 1091 fInnerRadius, fOuterRadius, fKappa, 858 -1 ) 1092 -1 ) ; 859 1093 860 fInnerHype = new G4TwistTubsHypeSide("Inner 1094 fInnerHype = new G4TwistTubsHypeSide("InnerHype", 861 fEndIn 1095 fEndInnerRadius, fEndOuterRadius, 862 fDPhi, 1096 fDPhi, fEndPhi, fEndZ, 863 fInner 1097 fInnerRadius, fOuterRadius,fKappa, 864 fTanIn 1098 fTanInnerStereo, fTanOuterStereo, -1) ; 865 fOuterHype = new G4TwistTubsHypeSide("Outer 1099 fOuterHype = new G4TwistTubsHypeSide("OuterHype", 866 fEndIn 1100 fEndInnerRadius, fEndOuterRadius, 867 fDPhi, 1101 fDPhi, fEndPhi, fEndZ, 868 fInner 1102 fInnerRadius, fOuterRadius,fKappa, 869 fTanIn 1103 fTanInnerStereo, fTanOuterStereo, 1) ; 870 1104 871 1105 872 // set neighbour surfaces 1106 // set neighbour surfaces 873 // 1107 // 874 fLowerEndcap->SetNeighbours(fInnerHype, fLa 1108 fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted, 875 fOuterHype, fFo 1109 fOuterHype, fFormerTwisted); 876 fUpperEndcap->SetNeighbours(fInnerHype, fLa 1110 fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted, 877 fOuterHype, fFo 1111 fOuterHype, fFormerTwisted); 878 fLatterTwisted->SetNeighbours(fInnerHype, f 1112 fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap, 879 fOuterHype, f 1113 fOuterHype, fUpperEndcap); 880 fFormerTwisted->SetNeighbours(fInnerHype, f 1114 fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap, 881 fOuterHype, f 1115 fOuterHype, fUpperEndcap); 882 fInnerHype->SetNeighbours(fLatterTwisted, f 1116 fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap, 883 fFormerTwisted, f 1117 fFormerTwisted, fUpperEndcap); 884 fOuterHype->SetNeighbours(fLatterTwisted, f 1118 fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap, 885 fFormerTwisted, f 1119 fFormerTwisted, fUpperEndcap); 886 } 1120 } 887 1121 888 1122 889 //============================================ 1123 //===================================================================== 890 //* GetEntityType ---------------------------- 1124 //* GetEntityType ----------------------------------------------------- 891 1125 892 G4GeometryType G4TwistedTubs::GetEntityType() 1126 G4GeometryType G4TwistedTubs::GetEntityType() const 893 { 1127 { 894 return {"G4TwistedTubs"}; << 1128 return G4String("G4TwistedTubs"); 895 } << 896 << 897 //============================================ << 898 //* Clone ------------------------------------ << 899 << 900 G4VSolid* G4TwistedTubs::Clone() const << 901 { << 902 return new G4TwistedTubs(*this); << 903 } 1129 } 904 1130 905 //============================================ 1131 //===================================================================== 906 //* GetCubicVolume --------------------------- 1132 //* GetCubicVolume ---------------------------------------------------- 907 1133 908 G4double G4TwistedTubs::GetCubicVolume() 1134 G4double G4TwistedTubs::GetCubicVolume() 909 { 1135 { 910 if (fCubicVolume == 0.) << 1136 if(fCubicVolume != 0.) {;} 911 { << 1137 else { fCubicVolume = fDPhi*fZHalfLength*(fOuterRadius*fOuterRadius 912 G4double DPhi = GetDPhi(); << 1138 -fInnerRadius*fInnerRadius); } 913 G4double Z0 = GetEndZ(0); << 914 G4double Z1 = GetEndZ(1); << 915 G4double Ain = GetInnerRadius(); << 916 G4double Aout = GetOuterRadius(); << 917 G4double R0in = GetEndInnerRadius(0); << 918 G4double R1in = GetEndInnerRadius(1); << 919 G4double R0out = GetEndOuterRadius(0); << 920 G4double R1out = GetEndOuterRadius(1); << 921 << 922 // V_hyperboloid = pi*h*(2*a*a + R*R)/3 << 923 fCubicVolume = (2.*(Z1 - Z0)*(Aout + Ain)* << 924 + Z1*(R1out + R1in)*(R1out << 925 - Z0*(R0out + R0in)*(R0out << 926 } << 927 return fCubicVolume; 1139 return fCubicVolume; 928 } 1140 } 929 1141 930 //============================================ 1142 //===================================================================== 931 //* GetLateralArea --------------------------- << 932 << 933 G4double << 934 G4TwistedTubs::GetLateralArea(G4double a, G4do << 935 { << 936 if (z == 0) return 0.; << 937 G4double h = std::abs(z); << 938 G4double area = h*a; << 939 if (std::abs(a - r) > kCarTolerance) << 940 { << 941 G4double aa = a*a; << 942 G4double hh = h*h; << 943 G4double rr = r*r; << 944 G4double cc = aa*hh/(rr - aa); << 945 G4double k = std::sqrt(aa + cc)/cc; << 946 G4double kh = k*h; << 947 area = 0.5*a*(h*std::sqrt(1. + kh*kh) + st << 948 } << 949 return GetDPhi()*area; << 950 } << 951 << 952 //============================================ << 953 //* GetPhiCutArea ---------------------------- << 954 << 955 G4double << 956 G4TwistedTubs::GetPhiCutArea(G4double a, G4dou << 957 { << 958 if (GetDPhi() >= CLHEP::twopi || r <= 0 || z << 959 G4double h = std::abs(z); << 960 G4double area = h*a; << 961 if (GetPhiTwist() > kCarTolerance) << 962 { << 963 G4double sinw = std::sin(0.5*GetPhiTwist() << 964 G4double p = sinw*r/h; << 965 G4double q = sinw*r/a; << 966 G4double pp = p*p; << 967 G4double qq = q*q; << 968 G4double pq = p*q; << 969 G4double sqroot = std::sqrt(pp + qq + 1); << 970 area = (pq*sqroot + << 971 0.5*p*(pp + 3.)*std::atanh(q/sqroo << 972 0.5*q*(qq + 3.)*std::atanh(p/sqroo << 973 std::atan(sqroot/(pq)) - CLHEP::ha << 974 } << 975 return area; << 976 } << 977 << 978 //============================================ << 979 //* GetSurfaceArea --------------------------- 1143 //* GetSurfaceArea ---------------------------------------------------- 980 1144 981 G4double G4TwistedTubs::GetSurfaceArea() 1145 G4double G4TwistedTubs::GetSurfaceArea() 982 { 1146 { 983 if (fSurfaceArea == 0.) << 1147 if(fSurfaceArea != 0.) {;} 984 { << 1148 else { fSurfaceArea = G4VSolid::GetSurfaceArea(); } 985 G4double dphi = GetDPhi(); << 986 G4double Ainn = GetInnerRadius(); << 987 G4double Aout = GetOuterRadius(); << 988 G4double Rinn0 = GetEndInnerRadius(0); << 989 G4double Rout0 = GetEndOuterRadius(0); << 990 G4double Rinn1 = GetEndInnerRadius(1); << 991 G4double Rout1 = GetEndOuterRadius(1); << 992 G4double z0 = GetEndZ(0); << 993 G4double z1 = GetEndZ(1); << 994 << 995 G4double base0 = 0.5*dphi*(Rout0*Rout0 - R << 996 G4double inner0 = GetLateralArea(Ainn, Rin << 997 G4double outer0 = GetLateralArea(Aout, Rou << 998 G4double cut0 = << 999 GetPhiCutArea(Aout, Rout0, z0) - GetPhiC << 1000 << 1001 G4double base1 = base0; << 1002 G4double inner1 = inner0; << 1003 G4double outer1 = outer0; << 1004 G4double cut1 = cut0; << 1005 if (std::abs(z0) != std::abs(z1)) << 1006 { << 1007 base1 = 0.5*dphi*(Rout1*Rout1 - Rinn1*R << 1008 inner1 = GetLateralArea(Ainn, Rinn1, z1 << 1009 outer1 = GetLateralArea(Aout, Rout1, z1 << 1010 cut1 = << 1011 GetPhiCutArea(Aout, Rout1, z1) - GetPhi << 1012 } << 1013 fSurfaceArea = base0 + base1 + << 1014 ((z0*z1 < 0) ? << 1015 (inner0 + inner1 + outer0 + outer1 + 2. << 1016 std::abs(inner0 - inner1 + outer0 - out << 1017 } << 1018 return fSurfaceArea; 1149 return fSurfaceArea; 1019 } 1150 } 1020 1151 1021 //=========================================== 1152 //===================================================================== 1022 //* GetPointOnSurface ----------------------- 1153 //* GetPointOnSurface ------------------------------------------------- 1023 1154 1024 G4ThreeVector G4TwistedTubs::GetPointOnSurfac 1155 G4ThreeVector G4TwistedTubs::GetPointOnSurface() const 1025 { 1156 { 1026 1157 1027 G4double z = G4RandFlat::shoot(fEndZ[0],fEn << 1158 G4double z = CLHEP::RandFlat::shoot(fEndZ[0],fEndZ[1]); 1028 G4double phi , phimin, phimax ; 1159 G4double phi , phimin, phimax ; 1029 G4double x , xmin, xmax ; 1160 G4double x , xmin, xmax ; 1030 G4double r , rmin, rmax ; 1161 G4double r , rmin, rmax ; 1031 1162 1032 G4double a1 = fOuterHype->GetSurfaceArea() 1163 G4double a1 = fOuterHype->GetSurfaceArea() ; 1033 G4double a2 = fInnerHype->GetSurfaceArea() 1164 G4double a2 = fInnerHype->GetSurfaceArea() ; 1034 G4double a3 = fLatterTwisted->GetSurfaceAre 1165 G4double a3 = fLatterTwisted->GetSurfaceArea() ; 1035 G4double a4 = fFormerTwisted->GetSurfaceAre 1166 G4double a4 = fFormerTwisted->GetSurfaceArea() ; 1036 G4double a5 = fLowerEndcap->GetSurfaceArea( 1167 G4double a5 = fLowerEndcap->GetSurfaceArea() ; 1037 G4double a6 = fUpperEndcap->GetSurfaceArea( 1168 G4double a6 = fUpperEndcap->GetSurfaceArea() ; 1038 1169 1039 G4double chose = G4RandFlat::shoot(0.,a1 + << 1170 G4double chose = CLHEP::RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ; 1040 1171 1041 if(chose < a1) 1172 if(chose < a1) 1042 { 1173 { 1043 1174 1044 phimin = fOuterHype->GetBoundaryMin(z) ; 1175 phimin = fOuterHype->GetBoundaryMin(z) ; 1045 phimax = fOuterHype->GetBoundaryMax(z) ; 1176 phimax = fOuterHype->GetBoundaryMax(z) ; 1046 phi = G4RandFlat::shoot(phimin,phimax) ; << 1177 phi = CLHEP::RandFlat::shoot(phimin,phimax) ; 1047 1178 1048 return fOuterHype->SurfacePoint(phi,z,tru 1179 return fOuterHype->SurfacePoint(phi,z,true) ; 1049 1180 1050 } 1181 } 1051 else if ( (chose >= a1) && (chose < a1 + a2 1182 else if ( (chose >= a1) && (chose < a1 + a2 ) ) 1052 { 1183 { 1053 1184 1054 phimin = fInnerHype->GetBoundaryMin(z) ; 1185 phimin = fInnerHype->GetBoundaryMin(z) ; 1055 phimax = fInnerHype->GetBoundaryMax(z) ; 1186 phimax = fInnerHype->GetBoundaryMax(z) ; 1056 phi = G4RandFlat::shoot(phimin,phimax) ; << 1187 phi = CLHEP::RandFlat::shoot(phimin,phimax) ; 1057 1188 1058 return fInnerHype->SurfacePoint(phi,z,tru 1189 return fInnerHype->SurfacePoint(phi,z,true) ; 1059 1190 1060 } 1191 } 1061 else if ( (chose >= a1 + a2 ) && (chose < a 1192 else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) ) 1062 { 1193 { 1063 1194 1064 xmin = fLatterTwisted->GetBoundaryMin(z) 1195 xmin = fLatterTwisted->GetBoundaryMin(z) ; 1065 xmax = fLatterTwisted->GetBoundaryMax(z) 1196 xmax = fLatterTwisted->GetBoundaryMax(z) ; 1066 x = G4RandFlat::shoot(xmin,xmax) ; << 1197 x = CLHEP::RandFlat::shoot(xmin,xmax) ; 1067 1198 1068 return fLatterTwisted->SurfacePoint(x,z,t 1199 return fLatterTwisted->SurfacePoint(x,z,true) ; 1069 1200 1070 } 1201 } 1071 else if ( (chose >= a1 + a2 + a3 ) && (cho 1202 else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) ) 1072 { 1203 { 1073 1204 1074 xmin = fFormerTwisted->GetBoundaryMin(z) 1205 xmin = fFormerTwisted->GetBoundaryMin(z) ; 1075 xmax = fFormerTwisted->GetBoundaryMax(z) 1206 xmax = fFormerTwisted->GetBoundaryMax(z) ; 1076 x = G4RandFlat::shoot(xmin,xmax) ; << 1207 x = CLHEP::RandFlat::shoot(xmin,xmax) ; 1077 1208 1078 return fFormerTwisted->SurfacePoint(x,z,t 1209 return fFormerTwisted->SurfacePoint(x,z,true) ; 1079 } << 1210 >> 1211 } 1080 else if( (chose >= a1 + a2 + a3 + a4 )&&(c 1212 else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) ) 1081 { 1213 { >> 1214 1082 rmin = GetEndInnerRadius(0) ; 1215 rmin = GetEndInnerRadius(0) ; 1083 rmax = GetEndOuterRadius(0) ; 1216 rmax = GetEndOuterRadius(0) ; 1084 r = std::sqrt(G4RandFlat::shoot()*(sqr(rm << 1217 r = CLHEP::RandFlat::shoot(rmin,rmax) ; 1085 1218 1086 phimin = fLowerEndcap->GetBoundaryMin(r) 1219 phimin = fLowerEndcap->GetBoundaryMin(r) ; 1087 phimax = fLowerEndcap->GetBoundaryMax(r) 1220 phimax = fLowerEndcap->GetBoundaryMax(r) ; 1088 phi = G4RandFlat::shoot(phimin,phimax) << 1221 phi = CLHEP::RandFlat::shoot(phimin,phimax) ; 1089 1222 1090 return fLowerEndcap->SurfacePoint(phi,r,t 1223 return fLowerEndcap->SurfacePoint(phi,r,true) ; >> 1224 1091 } 1225 } 1092 else 1226 else 1093 { 1227 { 1094 rmin = GetEndInnerRadius(1) ; 1228 rmin = GetEndInnerRadius(1) ; 1095 rmax = GetEndOuterRadius(1) ; 1229 rmax = GetEndOuterRadius(1) ; 1096 r = rmin + (rmax-rmin)*std::sqrt(G4RandFl << 1230 r = CLHEP::RandFlat::shoot(rmin,rmax) ; 1097 1231 1098 phimin = fUpperEndcap->GetBoundaryMin(r) 1232 phimin = fUpperEndcap->GetBoundaryMin(r) ; 1099 phimax = fUpperEndcap->GetBoundaryMax(r) 1233 phimax = fUpperEndcap->GetBoundaryMax(r) ; 1100 phi = G4RandFlat::shoot(phimin,phimax) << 1234 phi = CLHEP::RandFlat::shoot(phimin,phimax) ; 1101 1235 1102 return fUpperEndcap->SurfacePoint(phi,r,t 1236 return fUpperEndcap->SurfacePoint(phi,r,true) ; 1103 } 1237 } 1104 } 1238 } 1105 1239