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 // Implementation for G4Trap class 26 // Implementation for G4Trap class 27 // 27 // 28 // 21.03.95 P.Kent: Modified for `tolerant' ge 28 // 21.03.95 P.Kent: Modified for `tolerant' geometry 29 // 09.09.96 V.Grichine: Final modifications be 29 // 09.09.96 V.Grichine: Final modifications before to commit 30 // 08.12.97 J.Allison: Added "nominal" constru 30 // 08.12.97 J.Allison: Added "nominal" constructor and method SetAllParameters 31 // 28.04.05 V.Grichine: new SurfaceNormal acco 31 // 28.04.05 V.Grichine: new SurfaceNormal according to J.Apostolakis proposal 32 // 18.04.17 E.Tcherniaev: complete revision, s 32 // 18.04.17 E.Tcherniaev: complete revision, speed-up 33 // ------------------------------------------- 33 // -------------------------------------------------------------------- 34 34 35 #include "G4Trap.hh" 35 #include "G4Trap.hh" 36 36 37 #if !defined(G4GEOM_USE_UTRAP) 37 #if !defined(G4GEOM_USE_UTRAP) 38 38 39 #include "globals.hh" 39 #include "globals.hh" 40 #include "G4GeomTools.hh" 40 #include "G4GeomTools.hh" 41 41 42 #include "G4VoxelLimits.hh" 42 #include "G4VoxelLimits.hh" 43 #include "G4AffineTransform.hh" 43 #include "G4AffineTransform.hh" 44 #include "G4BoundingEnvelope.hh" 44 #include "G4BoundingEnvelope.hh" 45 45 46 #include "G4VPVParameterisation.hh" 46 #include "G4VPVParameterisation.hh" 47 47 48 #include "G4QuickRand.hh" 48 #include "G4QuickRand.hh" 49 49 50 #include "G4VGraphicsScene.hh" 50 #include "G4VGraphicsScene.hh" 51 #include "G4Polyhedron.hh" 51 #include "G4Polyhedron.hh" 52 52 53 using namespace CLHEP; 53 using namespace CLHEP; 54 54 55 ////////////////////////////////////////////// 55 ////////////////////////////////////////////////////////////////////////// 56 // 56 // 57 // Constructor - check and set half-widths as 57 // Constructor - check and set half-widths as well as angles: 58 // final check of coplanarity 58 // final check of coplanarity 59 59 60 G4Trap::G4Trap( const G4String& pName, 60 G4Trap::G4Trap( const G4String& pName, 61 G4double pDz, 61 G4double pDz, 62 G4double pTheta, G4doubl 62 G4double pTheta, G4double pPhi, 63 G4double pDy1, G4double 63 G4double pDy1, G4double pDx1, G4double pDx2, 64 G4double pAlp1, 64 G4double pAlp1, 65 G4double pDy2, G4double 65 G4double pDy2, G4double pDx3, G4double pDx4, 66 G4double pAlp2 ) 66 G4double pAlp2 ) 67 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 67 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance) 68 { 68 { 69 fDz = pDz; 69 fDz = pDz; 70 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi 70 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi); 71 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi 71 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi); 72 72 73 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp 73 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalpha1 = std::tan(pAlp1); 74 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp 74 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalpha2 = std::tan(pAlp2); 75 75 76 CheckParameters(); 76 CheckParameters(); 77 MakePlanes(); 77 MakePlanes(); 78 } 78 } 79 79 80 ////////////////////////////////////////////// 80 ////////////////////////////////////////////////////////////////////////// 81 // 81 // 82 // Constructor - Design of trapezoid based on 82 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters, 83 // which are its vertices. Checking of planari 83 // which are its vertices. Checking of planarity with preparation of 84 // fPlanes[] and than calculation of other mem 84 // fPlanes[] and than calculation of other members 85 85 86 G4Trap::G4Trap( const G4String& pName, 86 G4Trap::G4Trap( const G4String& pName, 87 const G4ThreeVector pt[8] ) 87 const G4ThreeVector pt[8] ) 88 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 88 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance) 89 { 89 { 90 // Start with check of centering - the cente 90 // Start with check of centering - the center of gravity trap line 91 // should cross the origin of frame 91 // should cross the origin of frame 92 // 92 // 93 if ( pt[0].z() >= 0 93 if ( pt[0].z() >= 0 94 || pt[0].z() != pt[1].z() 94 || pt[0].z() != pt[1].z() 95 || pt[0].z() != pt[2].z() 95 || pt[0].z() != pt[2].z() 96 || pt[0].z() != pt[3].z() 96 || pt[0].z() != pt[3].z() 97 97 98 || pt[4].z() <= 0 98 || pt[4].z() <= 0 99 || pt[4].z() != pt[5].z() 99 || pt[4].z() != pt[5].z() 100 || pt[4].z() != pt[6].z() 100 || pt[4].z() != pt[6].z() 101 || pt[4].z() != pt[7].z() 101 || pt[4].z() != pt[7].z() 102 102 103 || std::fabs( pt[0].z() + pt[4].z() ) 103 || std::fabs( pt[0].z() + pt[4].z() ) >= kCarTolerance 104 104 105 || pt[0].y() != pt[1].y() 105 || pt[0].y() != pt[1].y() 106 || pt[2].y() != pt[3].y() 106 || pt[2].y() != pt[3].y() 107 || pt[4].y() != pt[5].y() 107 || pt[4].y() != pt[5].y() 108 || pt[6].y() != pt[7].y() 108 || pt[6].y() != pt[7].y() 109 109 110 || std::fabs(pt[0].y()+pt[2].y()+pt[4] 110 || std::fabs(pt[0].y()+pt[2].y()+pt[4].y()+pt[6].y()) >= kCarTolerance 111 || std::fabs(pt[0].x()+pt[1].x()+pt[4] 111 || std::fabs(pt[0].x()+pt[1].x()+pt[4].x()+pt[5].x() + 112 pt[2].x()+pt[3].x()+pt[6] 112 pt[2].x()+pt[3].x()+pt[6].x()+pt[7].x()) >= kCarTolerance ) 113 { 113 { 114 std::ostringstream message; 114 std::ostringstream message; 115 message << "Invalid vertice coordinates fo 115 message << "Invalid vertice coordinates for Solid: " << GetName(); 116 G4Exception("G4Trap::G4Trap()", "GeomSolid 116 G4Exception("G4Trap::G4Trap()", "GeomSolids0002", 117 FatalException, message); 117 FatalException, message); 118 } 118 } 119 119 120 // Set parameters 120 // Set parameters 121 // 121 // 122 fDz = (pt[7]).z(); 122 fDz = (pt[7]).z(); 123 123 124 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5; 124 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5; 125 fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5; 125 fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5; 126 fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5; 126 fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5; 127 fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]). 127 fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy1; 128 128 129 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5; 129 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5; 130 fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5; 130 fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5; 131 fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5; 131 fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5; 132 fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]). 132 fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy2; 133 133 134 fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx 134 fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx3)/fDz; 135 fTthetaSphi = ((pt[4]).y()+fDy2)/fDz; 135 fTthetaSphi = ((pt[4]).y()+fDy2)/fDz; 136 136 137 CheckParameters(); 137 CheckParameters(); 138 MakePlanes(pt); 138 MakePlanes(pt); 139 } 139 } 140 140 141 ////////////////////////////////////////////// 141 ////////////////////////////////////////////////////////////////////////// 142 // 142 // 143 // Constructor for Right Angular Wedge from ST 143 // Constructor for Right Angular Wedge from STEP 144 144 145 G4Trap::G4Trap( const G4String& pName, 145 G4Trap::G4Trap( const G4String& pName, 146 G4double pZ, 146 G4double pZ, 147 G4double pY, 147 G4double pY, 148 G4double pX, G4double pL 148 G4double pX, G4double pLTX ) 149 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 149 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance) 150 { 150 { 151 fDz = 0.5*pZ; fTthetaCphi = 0; fTthetaSphi 151 fDz = 0.5*pZ; fTthetaCphi = 0; fTthetaSphi = 0; 152 fDy1 = 0.5*pY; fDx1 = 0.5*pX; fDx2 = 0.5*pLT 152 fDy1 = 0.5*pY; fDx1 = 0.5*pX; fDx2 = 0.5*pLTX; fTalpha1 = 0.5*(pLTX - pX)/pY; 153 fDy2 = fDy1; fDx3 = fDx1; fDx4 = fDx2; 153 fDy2 = fDy1; fDx3 = fDx1; fDx4 = fDx2; fTalpha2 = fTalpha1; 154 154 155 CheckParameters(); 155 CheckParameters(); 156 MakePlanes(); 156 MakePlanes(); 157 } 157 } 158 158 159 ////////////////////////////////////////////// 159 ////////////////////////////////////////////////////////////////////////// 160 // 160 // 161 // Constructor for G4Trd 161 // Constructor for G4Trd 162 162 163 G4Trap::G4Trap( const G4String& pName, 163 G4Trap::G4Trap( const G4String& pName, 164 G4double pDx1, G4double 164 G4double pDx1, G4double pDx2, 165 G4double pDy1, G4double 165 G4double pDy1, G4double pDy2, 166 G4double pDz ) 166 G4double pDz ) 167 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 167 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance), fTrapType(0) 168 { 168 { 169 fDz = pDz; fTthetaCphi = 0; fTthetaSphi = 169 fDz = pDz; fTthetaCphi = 0; fTthetaSphi = 0; 170 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx1; fTalp 170 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx1; fTalpha1 = 0; 171 fDy2 = pDy2; fDx3 = pDx2; fDx4 = pDx2; fTalp 171 fDy2 = pDy2; fDx3 = pDx2; fDx4 = pDx2; fTalpha2 = 0; 172 172 173 CheckParameters(); 173 CheckParameters(); 174 MakePlanes(); 174 MakePlanes(); 175 } 175 } 176 176 177 ////////////////////////////////////////////// 177 ////////////////////////////////////////////////////////////////////////// 178 // 178 // 179 // Constructor for G4Para 179 // Constructor for G4Para 180 180 181 G4Trap::G4Trap( const G4String& pName, 181 G4Trap::G4Trap( const G4String& pName, 182 G4double pDx, G4double p 182 G4double pDx, G4double pDy, 183 G4double pDz, 183 G4double pDz, 184 G4double pAlpha, 184 G4double pAlpha, 185 G4double pTheta, G4doubl 185 G4double pTheta, G4double pPhi ) 186 : G4CSGSolid(pName), halfCarTolerance(0.5*kC 186 : G4CSGSolid(pName), halfCarTolerance(0.5*kCarTolerance) 187 { 187 { 188 fDz = pDz; 188 fDz = pDz; 189 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi 189 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi); 190 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi 190 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi); 191 191 192 fDy1 = pDy; fDx1 = pDx; fDx2 = pDx; fTalpha1 192 fDy1 = pDy; fDx1 = pDx; fDx2 = pDx; fTalpha1 = std::tan(pAlpha); 193 fDy2 = pDy; fDx3 = pDx; fDx4 = pDx; fTalpha2 193 fDy2 = pDy; fDx3 = pDx; fDx4 = pDx; fTalpha2 = fTalpha1; 194 194 195 CheckParameters(); 195 CheckParameters(); 196 MakePlanes(); 196 MakePlanes(); 197 } 197 } 198 198 199 ////////////////////////////////////////////// 199 ////////////////////////////////////////////////////////////////////////// 200 // 200 // 201 // Nominal constructor for G4Trap whose parame 201 // Nominal constructor for G4Trap whose parameters are to be set by 202 // a G4VParamaterisation later. Check and set 202 // a G4VParamaterisation later. Check and set half-widths as well as 203 // angles: final check of coplanarity 203 // angles: final check of coplanarity 204 204 205 G4Trap::G4Trap( const G4String& pName ) 205 G4Trap::G4Trap( const G4String& pName ) 206 : G4CSGSolid (pName), halfCarTolerance(0.5*k 206 : G4CSGSolid (pName), halfCarTolerance(0.5*kCarTolerance), 207 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), 207 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), 208 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.) 208 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.), 209 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.) 209 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.) 210 { 210 { 211 MakePlanes(); 211 MakePlanes(); 212 } 212 } 213 213 214 ////////////////////////////////////////////// 214 ////////////////////////////////////////////////////////////////////////// 215 // 215 // 216 // Fake default constructor - sets only member 216 // Fake default constructor - sets only member data and allocates memory 217 // for usage restri 217 // for usage restricted to object persistency. 218 // 218 // 219 G4Trap::G4Trap( __void__& a ) 219 G4Trap::G4Trap( __void__& a ) 220 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo 220 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTolerance), 221 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), 221 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), 222 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.) 222 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.), 223 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.) 223 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.) 224 { 224 { 225 MakePlanes(); 225 MakePlanes(); 226 } 226 } 227 227 228 ////////////////////////////////////////////// 228 ////////////////////////////////////////////////////////////////////////// 229 // 229 // 230 // Destructor 230 // Destructor 231 231 232 G4Trap::~G4Trap() = default; 232 G4Trap::~G4Trap() = default; 233 233 234 ////////////////////////////////////////////// 234 ////////////////////////////////////////////////////////////////////////// 235 // 235 // 236 // Copy constructor 236 // Copy constructor 237 237 238 G4Trap::G4Trap(const G4Trap& rhs) 238 G4Trap::G4Trap(const G4Trap& rhs) 239 : G4CSGSolid(rhs), halfCarTolerance(rhs.half 239 : G4CSGSolid(rhs), halfCarTolerance(rhs.halfCarTolerance), 240 fDz(rhs.fDz), fTthetaCphi(rhs.fTthetaCphi) 240 fDz(rhs.fDz), fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi), 241 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.f 241 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1), 242 fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.f 242 fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2) 243 { 243 { 244 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs 244 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; } 245 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs. 245 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs.fAreas[i]; } 246 fTrapType = rhs.fTrapType; 246 fTrapType = rhs.fTrapType; 247 } 247 } 248 248 249 ////////////////////////////////////////////// 249 ////////////////////////////////////////////////////////////////////////// 250 // 250 // 251 // Assignment operator 251 // Assignment operator 252 252 253 G4Trap& G4Trap::operator = (const G4Trap& rhs) 253 G4Trap& G4Trap::operator = (const G4Trap& rhs) 254 { 254 { 255 // Check assignment to self 255 // Check assignment to self 256 // 256 // 257 if (this == &rhs) { return *this; } 257 if (this == &rhs) { return *this; } 258 258 259 // Copy base class data 259 // Copy base class data 260 // 260 // 261 G4CSGSolid::operator=(rhs); 261 G4CSGSolid::operator=(rhs); 262 262 263 // Copy data 263 // Copy data 264 // 264 // 265 halfCarTolerance = rhs.halfCarTolerance; 265 halfCarTolerance = rhs.halfCarTolerance; 266 fDz = rhs.fDz; fTthetaCphi = rhs.fTthetaCphi 266 fDz = rhs.fDz; fTthetaCphi = rhs.fTthetaCphi; fTthetaSphi = rhs.fTthetaSphi; 267 fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs 267 fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1; 268 fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs 268 fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2; 269 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs 269 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs.fPlanes[i]; } 270 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs. 270 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs.fAreas[i]; } 271 fTrapType = rhs.fTrapType; 271 fTrapType = rhs.fTrapType; 272 return *this; 272 return *this; 273 } 273 } 274 274 275 ////////////////////////////////////////////// 275 ////////////////////////////////////////////////////////////////////////// 276 // 276 // 277 // Set all parameters, as for constructor - ch 277 // Set all parameters, as for constructor - check and set half-widths 278 // as well as angles: final check of coplanari 278 // as well as angles: final check of coplanarity 279 279 280 void G4Trap::SetAllParameters ( G4double pDz, 280 void G4Trap::SetAllParameters ( G4double pDz, 281 G4double pThet 281 G4double pTheta, 282 G4double pPhi, 282 G4double pPhi, 283 G4double pDy1, 283 G4double pDy1, 284 G4double pDx1, 284 G4double pDx1, 285 G4double pDx2, 285 G4double pDx2, 286 G4double pAlp1 286 G4double pAlp1, 287 G4double pDy2, 287 G4double pDy2, 288 G4double pDx3, 288 G4double pDx3, 289 G4double pDx4, 289 G4double pDx4, 290 G4double pAlp2 290 G4double pAlp2 ) 291 { 291 { 292 // Reset data of the base class 292 // Reset data of the base class 293 fCubicVolume = 0; 293 fCubicVolume = 0; 294 fSurfaceArea = 0; 294 fSurfaceArea = 0; 295 fRebuildPolyhedron = true; 295 fRebuildPolyhedron = true; 296 296 297 // Set parameters 297 // Set parameters 298 fDz = pDz; 298 fDz = pDz; 299 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi 299 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi); 300 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi 300 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi); 301 301 302 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp 302 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalpha1 = std::tan(pAlp1); 303 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp 303 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalpha2 = std::tan(pAlp2); 304 304 305 CheckParameters(); 305 CheckParameters(); 306 MakePlanes(); 306 MakePlanes(); 307 } 307 } 308 308 309 ////////////////////////////////////////////// 309 ////////////////////////////////////////////////////////////////////////// 310 // 310 // 311 // Check length parameters 311 // Check length parameters 312 312 313 void G4Trap::CheckParameters() 313 void G4Trap::CheckParameters() 314 { 314 { 315 if (fDz<=0 || 315 if (fDz<=0 || 316 fDy1<=0 || fDx1<=0 || fDx2<=0 || 316 fDy1<=0 || fDx1<=0 || fDx2<=0 || 317 fDy2<=0 || fDx3<=0 || fDx4<=0) 317 fDy2<=0 || fDx3<=0 || fDx4<=0) 318 { 318 { 319 std::ostringstream message; 319 std::ostringstream message; 320 message << "Invalid Length Parameters for 320 message << "Invalid Length Parameters for Solid: " << GetName() 321 << "\n X - " <<fDx1<<", "<<fDx2<< 321 << "\n X - " <<fDx1<<", "<<fDx2<<", "<<fDx3<<", "<<fDx4 322 << "\n Y - " <<fDy1<<", "<<fDy2 322 << "\n Y - " <<fDy1<<", "<<fDy2 323 << "\n Z - " <<fDz; 323 << "\n Z - " <<fDz; 324 G4Exception("G4Trap::CheckParameters()", " 324 G4Exception("G4Trap::CheckParameters()", "GeomSolids0002", 325 FatalException, message); 325 FatalException, message); 326 } 326 } 327 } 327 } 328 328 329 ////////////////////////////////////////////// 329 ////////////////////////////////////////////////////////////////////////// 330 // 330 // 331 // Compute vertices and set side planes 331 // Compute vertices and set side planes 332 332 333 void G4Trap::MakePlanes() 333 void G4Trap::MakePlanes() 334 { 334 { 335 G4double DzTthetaCphi = fDz*fTthetaCphi; 335 G4double DzTthetaCphi = fDz*fTthetaCphi; 336 G4double DzTthetaSphi = fDz*fTthetaSphi; 336 G4double DzTthetaSphi = fDz*fTthetaSphi; 337 G4double Dy1Talpha1 = fDy1*fTalpha1; 337 G4double Dy1Talpha1 = fDy1*fTalpha1; 338 G4double Dy2Talpha2 = fDy2*fTalpha2; 338 G4double Dy2Talpha2 = fDy2*fTalpha2; 339 339 340 G4ThreeVector pt[8] = 340 G4ThreeVector pt[8] = 341 { 341 { 342 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx 342 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx1,-DzTthetaSphi-fDy1,-fDz), 343 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx 343 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx1,-DzTthetaSphi-fDy1,-fDz), 344 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx 344 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx2,-DzTthetaSphi+fDy1,-fDz), 345 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx 345 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx2,-DzTthetaSphi+fDy1,-fDz), 346 G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx 346 G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx3, DzTthetaSphi-fDy2, fDz), 347 G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx 347 G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx3, DzTthetaSphi-fDy2, fDz), 348 G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx 348 G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx4, DzTthetaSphi+fDy2, fDz), 349 G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx 349 G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx4, DzTthetaSphi+fDy2, fDz) 350 }; 350 }; 351 351 352 MakePlanes(pt); 352 MakePlanes(pt); 353 } 353 } 354 354 355 ////////////////////////////////////////////// 355 ////////////////////////////////////////////////////////////////////////// 356 // 356 // 357 // Set side planes, check planarity 357 // Set side planes, check planarity 358 358 359 void G4Trap::MakePlanes(const G4ThreeVector pt 359 void G4Trap::MakePlanes(const G4ThreeVector pt[8]) 360 { 360 { 361 constexpr G4int iface[4][4] = { {0,4,5,1}, { 361 constexpr G4int iface[4][4] = { {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3} }; 362 const static G4String side[4] = { "~-Y", "~+ 362 const static G4String side[4] = { "~-Y", "~+Y", "~-X", "~+X" }; 363 363 364 for (G4int i=0; i<4; ++i) 364 for (G4int i=0; i<4; ++i) 365 { 365 { 366 if (MakePlane(pt[iface[i][0]], 366 if (MakePlane(pt[iface[i][0]], 367 pt[iface[i][1]], 367 pt[iface[i][1]], 368 pt[iface[i][2]], 368 pt[iface[i][2]], 369 pt[iface[i][3]], 369 pt[iface[i][3]], 370 fPlanes[i])) continue; 370 fPlanes[i])) continue; 371 371 372 // Non planar side face 372 // Non planar side face 373 G4ThreeVector normal(fPlanes[i].a,fPlanes[ 373 G4ThreeVector normal(fPlanes[i].a,fPlanes[i].b,fPlanes[i].c); 374 G4double dmax = 0; 374 G4double dmax = 0; 375 for (G4int k=0; k<4; ++k) 375 for (G4int k=0; k<4; ++k) 376 { 376 { 377 G4double dist = normal.dot(pt[iface[i][k 377 G4double dist = normal.dot(pt[iface[i][k]]) + fPlanes[i].d; 378 if (std::abs(dist) > std::abs(dmax)) dma 378 if (std::abs(dist) > std::abs(dmax)) dmax = dist; 379 } 379 } 380 std::ostringstream message; 380 std::ostringstream message; 381 message << "Side face " << side[i] << " is 381 message << "Side face " << side[i] << " is not planar for solid: " 382 << GetName() << "\nDiscrepancy: " 382 << GetName() << "\nDiscrepancy: " << dmax/mm << " mm\n"; 383 StreamInfo(message); 383 StreamInfo(message); 384 G4Exception("G4Trap::MakePlanes()", "GeomS 384 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002", 385 FatalException, message); 385 FatalException, message); 386 } 386 } 387 387 388 // Re-compute parameters 388 // Re-compute parameters 389 SetCachedValues(); 389 SetCachedValues(); 390 } 390 } 391 391 392 ////////////////////////////////////////////// 392 ////////////////////////////////////////////////////////////////////////// 393 // 393 // 394 // Calculate the coef's of the plane p1->p2->p 394 // Calculate the coef's of the plane p1->p2->p3->p4->p1 395 // where the ThreeVectors 1-4 are in anti-cloc 395 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed 396 // from infront of the plane (i.e. from normal 396 // from infront of the plane (i.e. from normal direction). 397 // 397 // 398 // Return true if the points are coplanar, fal 398 // Return true if the points are coplanar, false otherwise 399 399 400 G4bool G4Trap::MakePlane( const G4ThreeVector& 400 G4bool G4Trap::MakePlane( const G4ThreeVector& p1, 401 const G4ThreeVector& 401 const G4ThreeVector& p2, 402 const G4ThreeVector& 402 const G4ThreeVector& p3, 403 const G4ThreeVector& 403 const G4ThreeVector& p4, 404 TrapSidePlane& 404 TrapSidePlane& plane ) 405 { 405 { 406 G4ThreeVector normal = ((p4 - p2).cross(p3 - 406 G4ThreeVector normal = ((p4 - p2).cross(p3 - p1)).unit(); 407 if (std::abs(normal.x()) < DBL_EPSILON) norm 407 if (std::abs(normal.x()) < DBL_EPSILON) normal.setX(0); 408 if (std::abs(normal.y()) < DBL_EPSILON) norm 408 if (std::abs(normal.y()) < DBL_EPSILON) normal.setY(0); 409 if (std::abs(normal.z()) < DBL_EPSILON) norm 409 if (std::abs(normal.z()) < DBL_EPSILON) normal.setZ(0); 410 normal = normal.unit(); 410 normal = normal.unit(); 411 411 412 G4ThreeVector centre = (p1 + p2 + p3 + p4)*0 412 G4ThreeVector centre = (p1 + p2 + p3 + p4)*0.25; 413 plane.a = normal.x(); 413 plane.a = normal.x(); 414 plane.b = normal.y(); 414 plane.b = normal.y(); 415 plane.c = normal.z(); 415 plane.c = normal.z(); 416 plane.d = -normal.dot(centre); 416 plane.d = -normal.dot(centre); 417 417 418 // compute distances and check planarity 418 // compute distances and check planarity 419 G4double d1 = std::abs(normal.dot(p1) + plan 419 G4double d1 = std::abs(normal.dot(p1) + plane.d); 420 G4double d2 = std::abs(normal.dot(p2) + plan 420 G4double d2 = std::abs(normal.dot(p2) + plane.d); 421 G4double d3 = std::abs(normal.dot(p3) + plan 421 G4double d3 = std::abs(normal.dot(p3) + plane.d); 422 G4double d4 = std::abs(normal.dot(p4) + plan 422 G4double d4 = std::abs(normal.dot(p4) + plane.d); 423 G4double dmax = std::max(std::max(std::max(d 423 G4double dmax = std::max(std::max(std::max(d1,d2),d3),d4); 424 424 425 return dmax <= 1000 * kCarTolerance; 425 return dmax <= 1000 * kCarTolerance; 426 } 426 } 427 427 428 ////////////////////////////////////////////// 428 ////////////////////////////////////////////////////////////////////////// 429 // 429 // 430 // Recompute parameters using planes 430 // Recompute parameters using planes 431 431 432 void G4Trap::SetCachedValues() 432 void G4Trap::SetCachedValues() 433 { 433 { 434 // Set indeces 434 // Set indeces 435 constexpr G4int iface[6][4] = 435 constexpr G4int iface[6][4] = 436 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2, 436 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} }; 437 437 438 // Get vertices 438 // Get vertices 439 G4ThreeVector pt[8]; 439 G4ThreeVector pt[8]; 440 GetVertices(pt); 440 GetVertices(pt); 441 441 442 // Set face areas 442 // Set face areas 443 for (G4int i=0; i<6; ++i) 443 for (G4int i=0; i<6; ++i) 444 { 444 { 445 fAreas[i] = G4GeomTools::QuadAreaNormal(pt 445 fAreas[i] = G4GeomTools::QuadAreaNormal(pt[iface[i][0]], 446 pt 446 pt[iface[i][1]], 447 pt 447 pt[iface[i][2]], 448 pt 448 pt[iface[i][3]]).mag(); 449 } 449 } 450 for (G4int i=1; i<6; ++i) { fAreas[i] += fAr 450 for (G4int i=1; i<6; ++i) { fAreas[i] += fAreas[i - 1]; } 451 451 452 // Define type of trapezoid 452 // Define type of trapezoid 453 fTrapType = 0; 453 fTrapType = 0; 454 if (fPlanes[0].b == -1 && fPlanes[1].b == 1 454 if (fPlanes[0].b == -1 && fPlanes[1].b == 1 && 455 std::abs(fPlanes[0].a) < DBL_EPSILON && 455 std::abs(fPlanes[0].a) < DBL_EPSILON && 456 std::abs(fPlanes[0].c) < DBL_EPSILON && 456 std::abs(fPlanes[0].c) < DBL_EPSILON && 457 std::abs(fPlanes[1].a) < DBL_EPSILON && 457 std::abs(fPlanes[1].a) < DBL_EPSILON && 458 std::abs(fPlanes[1].c) < DBL_EPSILON) 458 std::abs(fPlanes[1].c) < DBL_EPSILON) 459 { 459 { 460 fTrapType = 1; // YZ section is a rectangl 460 fTrapType = 1; // YZ section is a rectangle ... 461 if (std::abs(fPlanes[2].a + fPlanes[3].a) 461 if (std::abs(fPlanes[2].a + fPlanes[3].a) < DBL_EPSILON && 462 std::abs(fPlanes[2].c - fPlanes[3].c) 462 std::abs(fPlanes[2].c - fPlanes[3].c) < DBL_EPSILON && 463 fPlanes[2].b == 0 && 463 fPlanes[2].b == 0 && 464 fPlanes[3].b == 0) 464 fPlanes[3].b == 0) 465 { 465 { 466 fTrapType = 2; // ... and XZ section is 466 fTrapType = 2; // ... and XZ section is a isosceles trapezoid 467 fPlanes[2].a = -fPlanes[3].a; 467 fPlanes[2].a = -fPlanes[3].a; 468 fPlanes[2].c = fPlanes[3].c; 468 fPlanes[2].c = fPlanes[3].c; 469 } 469 } 470 if (std::abs(fPlanes[2].a + fPlanes[3].a) 470 if (std::abs(fPlanes[2].a + fPlanes[3].a) < DBL_EPSILON && 471 std::abs(fPlanes[2].b - fPlanes[3].b) 471 std::abs(fPlanes[2].b - fPlanes[3].b) < DBL_EPSILON && 472 fPlanes[2].c == 0 && 472 fPlanes[2].c == 0 && 473 fPlanes[3].c == 0) 473 fPlanes[3].c == 0) 474 { 474 { 475 fTrapType = 3; // ... and XY section is 475 fTrapType = 3; // ... and XY section is a isosceles trapezoid 476 fPlanes[2].a = -fPlanes[3].a; 476 fPlanes[2].a = -fPlanes[3].a; 477 fPlanes[2].b = fPlanes[3].b; 477 fPlanes[2].b = fPlanes[3].b; 478 } 478 } 479 } 479 } 480 } 480 } 481 481 482 ////////////////////////////////////////////// 482 ////////////////////////////////////////////////////////////////////////// 483 // 483 // 484 // Get volume 484 // Get volume 485 485 486 G4double G4Trap::GetCubicVolume() 486 G4double G4Trap::GetCubicVolume() 487 { 487 { 488 if (fCubicVolume == 0) 488 if (fCubicVolume == 0) 489 { 489 { 490 G4ThreeVector pt[8]; 490 G4ThreeVector pt[8]; 491 GetVertices(pt); 491 GetVertices(pt); 492 492 493 G4double dz = pt[4].z() - pt[0].z(); 493 G4double dz = pt[4].z() - pt[0].z(); 494 G4double dy1 = pt[2].y() - pt[0].y(); 494 G4double dy1 = pt[2].y() - pt[0].y(); 495 G4double dx1 = pt[1].x() - pt[0].x(); 495 G4double dx1 = pt[1].x() - pt[0].x(); 496 G4double dx2 = pt[3].x() - pt[2].x(); 496 G4double dx2 = pt[3].x() - pt[2].x(); 497 G4double dy2 = pt[6].y() - pt[4].y(); 497 G4double dy2 = pt[6].y() - pt[4].y(); 498 G4double dx3 = pt[5].x() - pt[4].x(); 498 G4double dx3 = pt[5].x() - pt[4].x(); 499 G4double dx4 = pt[7].x() - pt[6].x(); 499 G4double dx4 = pt[7].x() - pt[6].x(); 500 500 501 fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(d 501 fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(dy1 + dy2) + 502 (dx4 + dx3 - dx2 - dx1)*(d 502 (dx4 + dx3 - dx2 - dx1)*(dy2 - dy1)/3)*dz*0.125; 503 } 503 } 504 return fCubicVolume; 504 return fCubicVolume; 505 } 505 } 506 506 507 ////////////////////////////////////////////// 507 ////////////////////////////////////////////////////////////////////////// 508 // 508 // 509 // Get surface area 509 // Get surface area 510 510 511 G4double G4Trap::GetSurfaceArea() 511 G4double G4Trap::GetSurfaceArea() 512 { 512 { 513 if (fSurfaceArea == 0) 513 if (fSurfaceArea == 0) 514 { 514 { 515 G4ThreeVector pt[8]; 515 G4ThreeVector pt[8]; 516 G4int iface [6][4] = 516 G4int iface [6][4] = 517 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2, 517 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} }; 518 518 519 GetVertices(pt); 519 GetVertices(pt); 520 for (const auto & i : iface) 520 for (const auto & i : iface) 521 { 521 { 522 fSurfaceArea += G4GeomTools::QuadAreaNor 522 fSurfaceArea += G4GeomTools::QuadAreaNormal(pt[i[0]], 523 523 pt[i[1]], 524 524 pt[i[2]], 525 525 pt[i[3]]).mag(); 526 } 526 } 527 } 527 } 528 return fSurfaceArea; 528 return fSurfaceArea; 529 } 529 } 530 530 531 ////////////////////////////////////////////// 531 ////////////////////////////////////////////////////////////////////////// 532 // 532 // 533 // Dispatch to parameterisation for replicatio 533 // Dispatch to parameterisation for replication mechanism dimension 534 // computation & modification. 534 // computation & modification. 535 535 536 void G4Trap::ComputeDimensions( G4VPVPar 536 void G4Trap::ComputeDimensions( G4VPVParameterisation* p, 537 const G4int n, 537 const G4int n, 538 const G4VPhysi 538 const G4VPhysicalVolume* pRep ) 539 { 539 { 540 p->ComputeDimensions(*this,n,pRep); 540 p->ComputeDimensions(*this,n,pRep); 541 } 541 } 542 542 543 ////////////////////////////////////////////// 543 ////////////////////////////////////////////////////////////////////////// 544 // 544 // 545 // Get bounding box 545 // Get bounding box 546 546 547 void G4Trap::BoundingLimits(G4ThreeVector& pMi 547 void G4Trap::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const 548 { 548 { 549 G4ThreeVector pt[8]; 549 G4ThreeVector pt[8]; 550 GetVertices(pt); 550 GetVertices(pt); 551 551 552 G4double xmin = kInfinity, xmax = -kInfinity 552 G4double xmin = kInfinity, xmax = -kInfinity; 553 G4double ymin = kInfinity, ymax = -kInfinity 553 G4double ymin = kInfinity, ymax = -kInfinity; 554 for (const auto & i : pt) 554 for (const auto & i : pt) 555 { 555 { 556 G4double x = i.x(); 556 G4double x = i.x(); 557 if (x < xmin) xmin = x; 557 if (x < xmin) xmin = x; 558 if (x > xmax) xmax = x; 558 if (x > xmax) xmax = x; 559 G4double y = i.y(); 559 G4double y = i.y(); 560 if (y < ymin) ymin = y; 560 if (y < ymin) ymin = y; 561 if (y > ymax) ymax = y; 561 if (y > ymax) ymax = y; 562 } 562 } 563 563 564 G4double dz = GetZHalfLength(); 564 G4double dz = GetZHalfLength(); 565 pMin.set(xmin,ymin,-dz); 565 pMin.set(xmin,ymin,-dz); 566 pMax.set(xmax,ymax, dz); 566 pMax.set(xmax,ymax, dz); 567 567 568 // Check correctness of the bounding box 568 // Check correctness of the bounding box 569 // 569 // 570 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 570 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z()) 571 { 571 { 572 std::ostringstream message; 572 std::ostringstream message; 573 message << "Bad bounding box (min >= max) 573 message << "Bad bounding box (min >= max) for solid: " 574 << GetName() << " !" 574 << GetName() << " !" 575 << "\npMin = " << pMin 575 << "\npMin = " << pMin 576 << "\npMax = " << pMax; 576 << "\npMax = " << pMax; 577 G4Exception("G4Trap::BoundingLimits()", "G 577 G4Exception("G4Trap::BoundingLimits()", "GeomMgt0001", 578 JustWarning, message); 578 JustWarning, message); 579 DumpInfo(); 579 DumpInfo(); 580 } 580 } 581 } 581 } 582 582 583 ////////////////////////////////////////////// 583 ////////////////////////////////////////////////////////////////////////// 584 // 584 // 585 // Calculate extent under transform and specif 585 // Calculate extent under transform and specified limit 586 586 587 G4bool G4Trap::CalculateExtent( const EAxis pA 587 G4bool G4Trap::CalculateExtent( const EAxis pAxis, 588 const G4VoxelL 588 const G4VoxelLimits& pVoxelLimit, 589 const G4Affine 589 const G4AffineTransform& pTransform, 590 G4double 590 G4double& pMin, G4double& pMax) const 591 { 591 { 592 G4ThreeVector bmin, bmax; 592 G4ThreeVector bmin, bmax; 593 G4bool exist; 593 G4bool exist; 594 594 595 // Check bounding box (bbox) 595 // Check bounding box (bbox) 596 // 596 // 597 BoundingLimits(bmin,bmax); 597 BoundingLimits(bmin,bmax); 598 G4BoundingEnvelope bbox(bmin,bmax); 598 G4BoundingEnvelope bbox(bmin,bmax); 599 #ifdef G4BBOX_EXTENT 599 #ifdef G4BBOX_EXTENT 600 return bbox.CalculateExtent(pAxis,pVoxelLimi 600 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 601 #endif 601 #endif 602 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 602 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax)) 603 { 603 { 604 return exist = pMin < pMax; 604 return exist = pMin < pMax; 605 } 605 } 606 606 607 // Set bounding envelope (benv) and calculat 607 // Set bounding envelope (benv) and calculate extent 608 // 608 // 609 G4ThreeVector pt[8]; 609 G4ThreeVector pt[8]; 610 GetVertices(pt); 610 GetVertices(pt); 611 611 612 G4ThreeVectorList baseA(4), baseB(4); 612 G4ThreeVectorList baseA(4), baseB(4); 613 baseA[0] = pt[0]; 613 baseA[0] = pt[0]; 614 baseA[1] = pt[1]; 614 baseA[1] = pt[1]; 615 baseA[2] = pt[3]; 615 baseA[2] = pt[3]; 616 baseA[3] = pt[2]; 616 baseA[3] = pt[2]; 617 617 618 baseB[0] = pt[4]; 618 baseB[0] = pt[4]; 619 baseB[1] = pt[5]; 619 baseB[1] = pt[5]; 620 baseB[2] = pt[7]; 620 baseB[2] = pt[7]; 621 baseB[3] = pt[6]; 621 baseB[3] = pt[6]; 622 622 623 std::vector<const G4ThreeVectorList *> polyg 623 std::vector<const G4ThreeVectorList *> polygons(2); 624 polygons[0] = &baseA; 624 polygons[0] = &baseA; 625 polygons[1] = &baseB; 625 polygons[1] = &baseB; 626 626 627 G4BoundingEnvelope benv(bmin,bmax,polygons); 627 G4BoundingEnvelope benv(bmin,bmax,polygons); 628 exist = benv.CalculateExtent(pAxis,pVoxelLim 628 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 629 return exist; 629 return exist; 630 } 630 } 631 631 632 ////////////////////////////////////////////// 632 ////////////////////////////////////////////////////////////////////////// 633 // 633 // 634 // Return whether point is inside/outside/on_s 634 // Return whether point is inside/outside/on_surface 635 635 636 EInside G4Trap::Inside( const G4ThreeVector& p 636 EInside G4Trap::Inside( const G4ThreeVector& p ) const 637 { 637 { 638 switch (fTrapType) 638 switch (fTrapType) 639 { 639 { 640 case 0: // General case 640 case 0: // General case 641 { 641 { 642 G4double dz = std::abs(p.z())-fDz; 642 G4double dz = std::abs(p.z())-fDz; 643 G4double dy1 = fPlanes[0].b*p.y()+fPlane 643 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d; 644 G4double dy2 = fPlanes[1].b*p.y()+fPlane 644 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d; 645 G4double dy = std::max(dz,std::max(dy1,d 645 G4double dy = std::max(dz,std::max(dy1,dy2)); 646 646 647 G4double dx1 = fPlanes[2].a*p.x()+fPlane 647 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y() 648 + fPlanes[2].c*p.z()+fPlane 648 + fPlanes[2].c*p.z()+fPlanes[2].d; 649 G4double dx2 = fPlanes[3].a*p.x()+fPlane 649 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y() 650 + fPlanes[3].c*p.z()+fPlane 650 + fPlanes[3].c*p.z()+fPlanes[3].d; 651 G4double dist = std::max(dy,std::max(dx1 651 G4double dist = std::max(dy,std::max(dx1,dx2)); 652 652 653 return (dist > halfCarTolerance) ? kOuts 653 return (dist > halfCarTolerance) ? kOutside : 654 ((dist > -halfCarTolerance) ? kSurface 654 ((dist > -halfCarTolerance) ? kSurface : kInside); 655 } 655 } 656 case 1: // YZ section is a rectangle 656 case 1: // YZ section is a rectangle 657 { 657 { 658 G4double dz = std::abs(p.z())-fDz; 658 G4double dz = std::abs(p.z())-fDz; 659 G4double dy = std::max(dz,std::abs(p.y() 659 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d); 660 G4double dx1 = fPlanes[2].a*p.x()+fPlane 660 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y() 661 + fPlanes[2].c*p.z()+fPlane 661 + fPlanes[2].c*p.z()+fPlanes[2].d; 662 G4double dx2 = fPlanes[3].a*p.x()+fPlane 662 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y() 663 + fPlanes[3].c*p.z()+fPlane 663 + fPlanes[3].c*p.z()+fPlanes[3].d; 664 G4double dist = std::max(dy,std::max(dx1 664 G4double dist = std::max(dy,std::max(dx1,dx2)); 665 665 666 return (dist > halfCarTolerance) ? kOuts 666 return (dist > halfCarTolerance) ? kOutside : 667 ((dist > -halfCarTolerance) ? kSurface 667 ((dist > -halfCarTolerance) ? kSurface : kInside); 668 } 668 } 669 case 2: // YZ section is a rectangle and 669 case 2: // YZ section is a rectangle and 670 { // XZ section is an isosceles trap 670 { // XZ section is an isosceles trapezoid 671 G4double dz = std::abs(p.z())-fDz; 671 G4double dz = std::abs(p.z())-fDz; 672 G4double dy = std::max(dz,std::abs(p.y() 672 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d); 673 G4double dx = fPlanes[3].a*std::abs(p.x( 673 G4double dx = fPlanes[3].a*std::abs(p.x()) 674 + fPlanes[3].c*p.z()+fPlanes 674 + fPlanes[3].c*p.z()+fPlanes[3].d; 675 G4double dist = std::max(dy,dx); 675 G4double dist = std::max(dy,dx); 676 676 677 return (dist > halfCarTolerance) ? kOuts 677 return (dist > halfCarTolerance) ? kOutside : 678 ((dist > -halfCarTolerance) ? kSurface 678 ((dist > -halfCarTolerance) ? kSurface : kInside); 679 } 679 } 680 case 3: // YZ section is a rectangle and 680 case 3: // YZ section is a rectangle and 681 { // XY section is an isosceles trap 681 { // XY section is an isosceles trapezoid 682 G4double dz = std::abs(p.z())-fDz; 682 G4double dz = std::abs(p.z())-fDz; 683 G4double dy = std::max(dz,std::abs(p.y() 683 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d); 684 G4double dx = fPlanes[3].a*std::abs(p.x( 684 G4double dx = fPlanes[3].a*std::abs(p.x()) 685 + fPlanes[3].b*p.y()+fPlanes 685 + fPlanes[3].b*p.y()+fPlanes[3].d; 686 G4double dist = std::max(dy,dx); 686 G4double dist = std::max(dy,dx); 687 687 688 return (dist > halfCarTolerance) ? kOuts 688 return (dist > halfCarTolerance) ? kOutside : 689 ((dist > -halfCarTolerance) ? kSurface 689 ((dist > -halfCarTolerance) ? kSurface : kInside); 690 } 690 } 691 } 691 } 692 return kOutside; 692 return kOutside; 693 } 693 } 694 694 695 ////////////////////////////////////////////// 695 ////////////////////////////////////////////////////////////////////////// 696 // 696 // 697 // Determine side, and return corresponding no 697 // Determine side, and return corresponding normal 698 698 699 G4ThreeVector G4Trap::SurfaceNormal( const G4T 699 G4ThreeVector G4Trap::SurfaceNormal( const G4ThreeVector& p ) const 700 { 700 { 701 G4double nx = 0, ny = 0, nz = 0; 701 G4double nx = 0, ny = 0, nz = 0; 702 G4double dz = std::abs(p.z()) - fDz; 702 G4double dz = std::abs(p.z()) - fDz; 703 nz = std::copysign(G4double(std::abs(dz) <= 703 nz = std::copysign(G4double(std::abs(dz) <= halfCarTolerance), p.z()); 704 704 705 switch (fTrapType) 705 switch (fTrapType) 706 { 706 { 707 case 0: // General case 707 case 0: // General case 708 { 708 { 709 for (G4int i=0; i<2; ++i) 709 for (G4int i=0; i<2; ++i) 710 { 710 { 711 G4double dy = fPlanes[i].b*p.y() + fPl 711 G4double dy = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d; 712 if (std::abs(dy) > halfCarTolerance) c 712 if (std::abs(dy) > halfCarTolerance) continue; 713 ny = fPlanes[i].b; 713 ny = fPlanes[i].b; 714 nz += fPlanes[i].c; 714 nz += fPlanes[i].c; 715 break; 715 break; 716 } 716 } 717 for (G4int i=2; i<4; ++i) 717 for (G4int i=2; i<4; ++i) 718 { 718 { 719 G4double dx = fPlanes[i].a*p.x() + 719 G4double dx = fPlanes[i].a*p.x() + 720 fPlanes[i].b*p.y() + fPl 720 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d; 721 if (std::abs(dx) > halfCarTolerance) c 721 if (std::abs(dx) > halfCarTolerance) continue; 722 nx = fPlanes[i].a; 722 nx = fPlanes[i].a; 723 ny += fPlanes[i].b; 723 ny += fPlanes[i].b; 724 nz += fPlanes[i].c; 724 nz += fPlanes[i].c; 725 break; 725 break; 726 } 726 } 727 break; 727 break; 728 } 728 } 729 case 1: // YZ section - rectangle 729 case 1: // YZ section - rectangle 730 { 730 { 731 G4double dy = std::abs(p.y()) + fPlanes[ 731 G4double dy = std::abs(p.y()) + fPlanes[1].d; 732 ny = std::copysign(G4double(std::abs(dy) 732 ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y()); 733 for (G4int i=2; i<4; ++i) 733 for (G4int i=2; i<4; ++i) 734 { 734 { 735 G4double dx = fPlanes[i].a*p.x() + 735 G4double dx = fPlanes[i].a*p.x() + 736 fPlanes[i].b*p.y() + fPl 736 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d; 737 if (std::abs(dx) > halfCarTolerance) c 737 if (std::abs(dx) > halfCarTolerance) continue; 738 nx = fPlanes[i].a; 738 nx = fPlanes[i].a; 739 ny += fPlanes[i].b; 739 ny += fPlanes[i].b; 740 nz += fPlanes[i].c; 740 nz += fPlanes[i].c; 741 break; 741 break; 742 } 742 } 743 break; 743 break; 744 } 744 } 745 case 2: // YZ section - rectangle, XZ sect 745 case 2: // YZ section - rectangle, XZ section - isosceles trapezoid 746 { 746 { 747 G4double dy = std::abs(p.y()) + fPlanes[ 747 G4double dy = std::abs(p.y()) + fPlanes[1].d; 748 ny = std::copysign(G4double(std::abs(dy) 748 ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y()); 749 G4double dx = fPlanes[3].a*std::abs(p.x( 749 G4double dx = fPlanes[3].a*std::abs(p.x()) + 750 fPlanes[3].c*p.z() + fPlan 750 fPlanes[3].c*p.z() + fPlanes[3].d; 751 G4double k = std::abs(dx) <= halfCarTole 751 G4double k = std::abs(dx) <= halfCarTolerance; 752 nx = std::copysign(k, p.x())*fPlanes[3] 752 nx = std::copysign(k, p.x())*fPlanes[3].a; 753 nz += k*fPlanes[3].c; 753 nz += k*fPlanes[3].c; 754 break; 754 break; 755 } 755 } 756 case 3: // YZ section - rectangle, XY sect 756 case 3: // YZ section - rectangle, XY section - isosceles trapezoid 757 { 757 { 758 G4double dy = std::abs(p.y()) + fPlanes[ 758 G4double dy = std::abs(p.y()) + fPlanes[1].d; 759 ny = std::copysign(G4double(std::abs(dy) 759 ny = std::copysign(G4double(std::abs(dy) <= halfCarTolerance), p.y()); 760 G4double dx = fPlanes[3].a*std::abs(p.x( 760 G4double dx = fPlanes[3].a*std::abs(p.x()) + 761 fPlanes[3].b*p.y() + fPlan 761 fPlanes[3].b*p.y() + fPlanes[3].d; 762 G4double k = std::abs(dx) <= halfCarTole 762 G4double k = std::abs(dx) <= halfCarTolerance; 763 nx = std::copysign(k, p.x())*fPlanes[3] 763 nx = std::copysign(k, p.x())*fPlanes[3].a; 764 ny += k*fPlanes[3].b; 764 ny += k*fPlanes[3].b; 765 break; 765 break; 766 } 766 } 767 } 767 } 768 768 769 // Return normal 769 // Return normal 770 // 770 // 771 G4double mag2 = nx*nx + ny*ny + nz*nz; 771 G4double mag2 = nx*nx + ny*ny + nz*nz; 772 if (mag2 == 1) return { nx,ny,nz }; 772 if (mag2 == 1) return { nx,ny,nz }; 773 else if (mag2 != 0) return G4ThreeVector(nx, 773 else if (mag2 != 0) return G4ThreeVector(nx,ny,nz).unit(); // edge or corner 774 else 774 else 775 { 775 { 776 // Point is not on the surface 776 // Point is not on the surface 777 // 777 // 778 #ifdef G4CSGDEBUG 778 #ifdef G4CSGDEBUG 779 std::ostringstream message; 779 std::ostringstream message; 780 G4long oldprc = message.precision(16); 780 G4long oldprc = message.precision(16); 781 message << "Point p is not on surface (!?) 781 message << "Point p is not on surface (!?) of solid: " 782 << GetName() << G4endl; 782 << GetName() << G4endl; 783 message << "Position:\n"; 783 message << "Position:\n"; 784 message << " p.x() = " << p.x()/mm << " 784 message << " p.x() = " << p.x()/mm << " mm\n"; 785 message << " p.y() = " << p.y()/mm << " 785 message << " p.y() = " << p.y()/mm << " mm\n"; 786 message << " p.z() = " << p.z()/mm << " 786 message << " p.z() = " << p.z()/mm << " mm"; 787 G4cout.precision(oldprc) ; 787 G4cout.precision(oldprc) ; 788 G4Exception("G4Trap::SurfaceNormal(p)", "G 788 G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002", 789 JustWarning, message ); 789 JustWarning, message ); 790 DumpInfo(); 790 DumpInfo(); 791 #endif 791 #endif 792 return ApproxSurfaceNormal(p); 792 return ApproxSurfaceNormal(p); 793 } 793 } 794 } 794 } 795 795 796 ////////////////////////////////////////////// 796 ////////////////////////////////////////////////////////////////////////// 797 // 797 // 798 // Algorithm for SurfaceNormal() following the 798 // Algorithm for SurfaceNormal() following the original specification 799 // for points not on the surface 799 // for points not on the surface 800 800 801 G4ThreeVector G4Trap::ApproxSurfaceNormal( con 801 G4ThreeVector G4Trap::ApproxSurfaceNormal( const G4ThreeVector& p ) const 802 { 802 { 803 G4double dist = -DBL_MAX; 803 G4double dist = -DBL_MAX; 804 G4int iside = 0; 804 G4int iside = 0; 805 for (G4int i=0; i<4; ++i) 805 for (G4int i=0; i<4; ++i) 806 { 806 { 807 G4double d = fPlanes[i].a*p.x() + 807 G4double d = fPlanes[i].a*p.x() + 808 fPlanes[i].b*p.y() + 808 fPlanes[i].b*p.y() + 809 fPlanes[i].c*p.z() + fPlanes[ 809 fPlanes[i].c*p.z() + fPlanes[i].d; 810 if (d > dist) { dist = d; iside = i; } 810 if (d > dist) { dist = d; iside = i; } 811 } 811 } 812 812 813 G4double distz = std::abs(p.z()) - fDz; 813 G4double distz = std::abs(p.z()) - fDz; 814 if (dist > distz) 814 if (dist > distz) 815 return { fPlanes[iside].a, fPlanes[iside]. 815 return { fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c }; 816 else 816 else 817 return { 0, 0, (G4double)((p.z() < 0) ? -1 817 return { 0, 0, (G4double)((p.z() < 0) ? -1 : 1) }; 818 } 818 } 819 819 820 ////////////////////////////////////////////// 820 ////////////////////////////////////////////////////////////////////////// 821 // 821 // 822 // Calculate distance to shape from outside 822 // Calculate distance to shape from outside 823 // - return kInfinity if no intersection 823 // - return kInfinity if no intersection 824 824 825 G4double G4Trap::DistanceToIn(const G4ThreeVec 825 G4double G4Trap::DistanceToIn(const G4ThreeVector& p, 826 const G4ThreeVec 826 const G4ThreeVector& v ) const 827 { 827 { 828 // Z intersections 828 // Z intersections 829 // 829 // 830 if ((std::abs(p.z()) - fDz) >= -halfCarToler 830 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() >= 0) 831 return kInfinity; 831 return kInfinity; 832 G4double invz = (-v.z() == 0) ? DBL_MAX : -1 832 G4double invz = (-v.z() == 0) ? DBL_MAX : -1./v.z(); 833 G4double dz = (invz < 0) ? fDz : -fDz; 833 G4double dz = (invz < 0) ? fDz : -fDz; 834 G4double tzmin = (p.z() + dz)*invz; 834 G4double tzmin = (p.z() + dz)*invz; 835 G4double tzmax = (p.z() - dz)*invz; 835 G4double tzmax = (p.z() - dz)*invz; 836 836 837 // Y intersections 837 // Y intersections 838 // 838 // 839 G4double tymin = 0, tymax = DBL_MAX; 839 G4double tymin = 0, tymax = DBL_MAX; 840 G4int i = 0; 840 G4int i = 0; 841 for ( ; i<2; ++i) 841 for ( ; i<2; ++i) 842 { 842 { 843 G4double cosa = fPlanes[i].b*v.y() + fPlan 843 G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z(); 844 G4double dist = fPlanes[i].b*p.y() + fPlan 844 G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d; 845 if (dist >= -halfCarTolerance) 845 if (dist >= -halfCarTolerance) 846 { 846 { 847 if (cosa >= 0) return kInfinity; 847 if (cosa >= 0) return kInfinity; 848 G4double tmp = -dist/cosa; 848 G4double tmp = -dist/cosa; 849 if (tymin < tmp) tymin = tmp; 849 if (tymin < tmp) tymin = tmp; 850 } 850 } 851 else if (cosa > 0) 851 else if (cosa > 0) 852 { 852 { 853 G4double tmp = -dist/cosa; 853 G4double tmp = -dist/cosa; 854 if (tymax > tmp) tymax = tmp; 854 if (tymax > tmp) tymax = tmp; 855 } 855 } 856 } 856 } 857 857 858 // Z intersections 858 // Z intersections 859 // 859 // 860 G4double txmin = 0, txmax = DBL_MAX; 860 G4double txmin = 0, txmax = DBL_MAX; 861 for ( ; i<4; ++i) 861 for ( ; i<4; ++i) 862 { 862 { 863 G4double cosa = fPlanes[i].a*v.x()+fPlanes 863 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z(); 864 G4double dist = fPlanes[i].a*p.x()+fPlanes 864 G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].c*p.z() + 865 fPlanes[i].d; 865 fPlanes[i].d; 866 if (dist >= -halfCarTolerance) 866 if (dist >= -halfCarTolerance) 867 { 867 { 868 if (cosa >= 0) return kInfinity; 868 if (cosa >= 0) return kInfinity; 869 G4double tmp = -dist/cosa; 869 G4double tmp = -dist/cosa; 870 if (txmin < tmp) txmin = tmp; 870 if (txmin < tmp) txmin = tmp; 871 } 871 } 872 else if (cosa > 0) 872 else if (cosa > 0) 873 { 873 { 874 G4double tmp = -dist/cosa; 874 G4double tmp = -dist/cosa; 875 if (txmax > tmp) txmax = tmp; 875 if (txmax > tmp) txmax = tmp; 876 } 876 } 877 } 877 } 878 878 879 // Find distance 879 // Find distance 880 // 880 // 881 G4double tmin = std::max(std::max(txmin,tymi 881 G4double tmin = std::max(std::max(txmin,tymin),tzmin); 882 G4double tmax = std::min(std::min(txmax,tyma 882 G4double tmax = std::min(std::min(txmax,tymax),tzmax); 883 883 884 if (tmax <= tmin + halfCarTolerance) return 884 if (tmax <= tmin + halfCarTolerance) return kInfinity; // touch or no hit 885 return (tmin < halfCarTolerance ) ? 0. : tmi 885 return (tmin < halfCarTolerance ) ? 0. : tmin; 886 } 886 } 887 887 888 ////////////////////////////////////////////// 888 ////////////////////////////////////////////////////////////////////////// 889 // 889 // 890 // Calculate exact shortest distance to any bo 890 // Calculate exact shortest distance to any boundary from outside 891 // This is the best fast estimation of the sho 891 // This is the best fast estimation of the shortest distance to trap 892 // - return 0 if point is inside 892 // - return 0 if point is inside 893 893 894 G4double G4Trap::DistanceToIn( const G4ThreeVe 894 G4double G4Trap::DistanceToIn( const G4ThreeVector& p ) const 895 { 895 { 896 switch (fTrapType) 896 switch (fTrapType) 897 { 897 { 898 case 0: // General case 898 case 0: // General case 899 { 899 { 900 G4double dz = std::abs(p.z())-fDz; 900 G4double dz = std::abs(p.z())-fDz; 901 G4double dy1 = fPlanes[0].b*p.y()+fPlane 901 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d; 902 G4double dy2 = fPlanes[1].b*p.y()+fPlane 902 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d; 903 G4double dy = std::max(dz,std::max(dy1,d 903 G4double dy = std::max(dz,std::max(dy1,dy2)); 904 904 905 G4double dx1 = fPlanes[2].a*p.x()+fPlane 905 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y() 906 + fPlanes[2].c*p.z()+fPlane 906 + fPlanes[2].c*p.z()+fPlanes[2].d; 907 G4double dx2 = fPlanes[3].a*p.x()+fPlane 907 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y() 908 + fPlanes[3].c*p.z()+fPlane 908 + fPlanes[3].c*p.z()+fPlanes[3].d; 909 G4double dist = std::max(dy,std::max(dx1 909 G4double dist = std::max(dy,std::max(dx1,dx2)); 910 return (dist > 0) ? dist : 0.; 910 return (dist > 0) ? dist : 0.; 911 } 911 } 912 case 1: // YZ section is a rectangle 912 case 1: // YZ section is a rectangle 913 { 913 { 914 G4double dz = std::abs(p.z())-fDz; 914 G4double dz = std::abs(p.z())-fDz; 915 G4double dy = std::max(dz,std::abs(p.y() 915 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d); 916 G4double dx1 = fPlanes[2].a*p.x()+fPlane 916 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y() 917 + fPlanes[2].c*p.z()+fPlane 917 + fPlanes[2].c*p.z()+fPlanes[2].d; 918 G4double dx2 = fPlanes[3].a*p.x()+fPlane 918 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y() 919 + fPlanes[3].c*p.z()+fPlane 919 + fPlanes[3].c*p.z()+fPlanes[3].d; 920 G4double dist = std::max(dy,std::max(dx1 920 G4double dist = std::max(dy,std::max(dx1,dx2)); 921 return (dist > 0) ? dist : 0.; 921 return (dist > 0) ? dist : 0.; 922 } 922 } 923 case 2: // YZ section is a rectangle and 923 case 2: // YZ section is a rectangle and 924 { // XZ section is an isosceles trap 924 { // XZ section is an isosceles trapezoid 925 G4double dz = std::abs(p.z())-fDz; 925 G4double dz = std::abs(p.z())-fDz; 926 G4double dy = std::max(dz,std::abs(p.y() 926 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d); 927 G4double dx = fPlanes[3].a*std::abs(p.x( 927 G4double dx = fPlanes[3].a*std::abs(p.x()) 928 + fPlanes[3].c*p.z()+fPlanes 928 + fPlanes[3].c*p.z()+fPlanes[3].d; 929 G4double dist = std::max(dy,dx); 929 G4double dist = std::max(dy,dx); 930 return (dist > 0) ? dist : 0.; 930 return (dist > 0) ? dist : 0.; 931 } 931 } 932 case 3: // YZ section is a rectangle and 932 case 3: // YZ section is a rectangle and 933 { // XY section is an isosceles trap 933 { // XY section is an isosceles trapezoid 934 G4double dz = std::abs(p.z())-fDz; 934 G4double dz = std::abs(p.z())-fDz; 935 G4double dy = std::max(dz,std::abs(p.y() 935 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d); 936 G4double dx = fPlanes[3].a*std::abs(p.x( 936 G4double dx = fPlanes[3].a*std::abs(p.x()) 937 + fPlanes[3].b*p.y()+fPlanes 937 + fPlanes[3].b*p.y()+fPlanes[3].d; 938 G4double dist = std::max(dy,dx); 938 G4double dist = std::max(dy,dx); 939 return (dist > 0) ? dist : 0.; 939 return (dist > 0) ? dist : 0.; 940 } 940 } 941 } 941 } 942 return 0.; 942 return 0.; 943 } 943 } 944 944 945 ////////////////////////////////////////////// 945 ////////////////////////////////////////////////////////////////////////// 946 // 946 // 947 // Calculate distance to surface of shape from 947 // Calculate distance to surface of shape from inside and 948 // find normal at exit point, if required 948 // find normal at exit point, if required 949 // - when leaving the surface, return 0 949 // - when leaving the surface, return 0 950 950 951 G4double G4Trap::DistanceToOut(const G4ThreeVe 951 G4double G4Trap::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v, 952 const G4bool ca 952 const G4bool calcNorm, 953 G4bool* v 953 G4bool* validNorm, G4ThreeVector* n) const 954 { 954 { 955 // Z intersections 955 // Z intersections 956 // 956 // 957 if ((std::abs(p.z()) - fDz) >= -halfCarToler 957 if ((std::abs(p.z()) - fDz) >= -halfCarTolerance && p.z()*v.z() > 0) 958 { 958 { 959 if (calcNorm) 959 if (calcNorm) 960 { 960 { 961 *validNorm = true; 961 *validNorm = true; 962 n->set(0, 0, (p.z() < 0) ? -1 : 1); 962 n->set(0, 0, (p.z() < 0) ? -1 : 1); 963 } 963 } 964 return 0; 964 return 0; 965 } 965 } 966 G4double vz = v.z(); 966 G4double vz = v.z(); 967 G4double tmax = (vz == 0) ? DBL_MAX : (std:: 967 G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - p.z())/vz; 968 G4int iside = (vz < 0) ? -4 : -2; // little 968 G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1 969 969 970 // Y intersections 970 // Y intersections 971 // 971 // 972 G4int i = 0; 972 G4int i = 0; 973 for ( ; i<2; ++i) 973 for ( ; i<2; ++i) 974 { 974 { 975 G4double cosa = fPlanes[i].b*v.y() + fPlan 975 G4double cosa = fPlanes[i].b*v.y() + fPlanes[i].c*v.z(); 976 if (cosa > 0) 976 if (cosa > 0) 977 { 977 { 978 G4double dist = fPlanes[i].b*p.y() + fPl 978 G4double dist = fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d; 979 if (dist >= -halfCarTolerance) 979 if (dist >= -halfCarTolerance) 980 { 980 { 981 if (calcNorm) 981 if (calcNorm) 982 { 982 { 983 *validNorm = true; 983 *validNorm = true; 984 n->set(0, fPlanes[i].b, fPlanes[i].c 984 n->set(0, fPlanes[i].b, fPlanes[i].c); 985 } 985 } 986 return 0; 986 return 0; 987 } 987 } 988 G4double tmp = -dist/cosa; 988 G4double tmp = -dist/cosa; 989 if (tmax > tmp) { tmax = tmp; iside = i; 989 if (tmax > tmp) { tmax = tmp; iside = i; } 990 } 990 } 991 } 991 } 992 992 993 // X intersections 993 // X intersections 994 // 994 // 995 for ( ; i<4; ++i) 995 for ( ; i<4; ++i) 996 { 996 { 997 G4double cosa = fPlanes[i].a*v.x()+fPlanes 997 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z(); 998 if (cosa > 0) 998 if (cosa > 0) 999 { 999 { 1000 G4double dist = fPlanes[i].a*p.x() + 1000 G4double dist = fPlanes[i].a*p.x() + 1001 fPlanes[i].b*p.y() + fP 1001 fPlanes[i].b*p.y() + fPlanes[i].c*p.z() + fPlanes[i].d; 1002 if (dist >= -halfCarTolerance) 1002 if (dist >= -halfCarTolerance) 1003 { 1003 { 1004 if (calcNorm) 1004 if (calcNorm) 1005 { 1005 { 1006 *validNorm = true; 1006 *validNorm = true; 1007 n->set(fPlanes[i].a, fPlanes[i].b, 1007 n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c); 1008 } 1008 } 1009 return 0; 1009 return 0; 1010 } 1010 } 1011 G4double tmp = -dist/cosa; 1011 G4double tmp = -dist/cosa; 1012 if (tmax > tmp) { tmax = tmp; iside = i 1012 if (tmax > tmp) { tmax = tmp; iside = i; } 1013 } 1013 } 1014 } 1014 } 1015 1015 1016 // Set normal, if required, and return dist 1016 // Set normal, if required, and return distance 1017 // 1017 // 1018 if (calcNorm) 1018 if (calcNorm) 1019 { 1019 { 1020 *validNorm = true; 1020 *validNorm = true; 1021 if (iside < 0) 1021 if (iside < 0) 1022 n->set(0, 0, iside + 3); // (-4+3)=-1, 1022 n->set(0, 0, iside + 3); // (-4+3)=-1, (-2+3)=+1 1023 else 1023 else 1024 n->set(fPlanes[iside].a, fPlanes[iside] 1024 n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c); 1025 } 1025 } 1026 return tmax; 1026 return tmax; 1027 } 1027 } 1028 1028 1029 ///////////////////////////////////////////// 1029 ////////////////////////////////////////////////////////////////////////// 1030 // 1030 // 1031 // Calculate exact shortest distance to any b 1031 // Calculate exact shortest distance to any boundary from inside 1032 // - Returns 0 is ThreeVector outside 1032 // - Returns 0 is ThreeVector outside 1033 1033 1034 G4double G4Trap::DistanceToOut( const G4Three 1034 G4double G4Trap::DistanceToOut( const G4ThreeVector& p ) const 1035 { 1035 { 1036 #ifdef G4CSGDEBUG 1036 #ifdef G4CSGDEBUG 1037 if( Inside(p) == kOutside ) 1037 if( Inside(p) == kOutside ) 1038 { 1038 { 1039 std::ostringstream message; 1039 std::ostringstream message; 1040 G4long oldprc = message.precision(16); 1040 G4long oldprc = message.precision(16); 1041 message << "Point p is outside (!?) of so 1041 message << "Point p is outside (!?) of solid: " << GetName() << G4endl; 1042 message << "Position:\n"; 1042 message << "Position:\n"; 1043 message << " p.x() = " << p.x()/mm << " 1043 message << " p.x() = " << p.x()/mm << " mm\n"; 1044 message << " p.y() = " << p.y()/mm << " 1044 message << " p.y() = " << p.y()/mm << " mm\n"; 1045 message << " p.z() = " << p.z()/mm << " 1045 message << " p.z() = " << p.z()/mm << " mm"; 1046 G4cout.precision(oldprc); 1046 G4cout.precision(oldprc); 1047 G4Exception("G4Trap::DistanceToOut(p)", " 1047 G4Exception("G4Trap::DistanceToOut(p)", "GeomSolids1002", 1048 JustWarning, message ); 1048 JustWarning, message ); 1049 DumpInfo(); 1049 DumpInfo(); 1050 } 1050 } 1051 #endif 1051 #endif 1052 switch (fTrapType) 1052 switch (fTrapType) 1053 { 1053 { 1054 case 0: // General case 1054 case 0: // General case 1055 { 1055 { 1056 G4double dz = std::abs(p.z())-fDz; 1056 G4double dz = std::abs(p.z())-fDz; 1057 G4double dy1 = fPlanes[0].b*p.y()+fPlan 1057 G4double dy1 = fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d; 1058 G4double dy2 = fPlanes[1].b*p.y()+fPlan 1058 G4double dy2 = fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d; 1059 G4double dy = std::max(dz,std::max(dy1, 1059 G4double dy = std::max(dz,std::max(dy1,dy2)); 1060 1060 1061 G4double dx1 = fPlanes[2].a*p.x()+fPlan 1061 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y() 1062 + fPlanes[2].c*p.z()+fPlan 1062 + fPlanes[2].c*p.z()+fPlanes[2].d; 1063 G4double dx2 = fPlanes[3].a*p.x()+fPlan 1063 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y() 1064 + fPlanes[3].c*p.z()+fPlan 1064 + fPlanes[3].c*p.z()+fPlanes[3].d; 1065 G4double dist = std::max(dy,std::max(dx 1065 G4double dist = std::max(dy,std::max(dx1,dx2)); 1066 return (dist < 0) ? -dist : 0.; 1066 return (dist < 0) ? -dist : 0.; 1067 } 1067 } 1068 case 1: // YZ section is a rectangle 1068 case 1: // YZ section is a rectangle 1069 { 1069 { 1070 G4double dz = std::abs(p.z())-fDz; 1070 G4double dz = std::abs(p.z())-fDz; 1071 G4double dy = std::max(dz,std::abs(p.y( 1071 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d); 1072 G4double dx1 = fPlanes[2].a*p.x()+fPlan 1072 G4double dx1 = fPlanes[2].a*p.x()+fPlanes[2].b*p.y() 1073 + fPlanes[2].c*p.z()+fPlan 1073 + fPlanes[2].c*p.z()+fPlanes[2].d; 1074 G4double dx2 = fPlanes[3].a*p.x()+fPlan 1074 G4double dx2 = fPlanes[3].a*p.x()+fPlanes[3].b*p.y() 1075 + fPlanes[3].c*p.z()+fPlan 1075 + fPlanes[3].c*p.z()+fPlanes[3].d; 1076 G4double dist = std::max(dy,std::max(dx 1076 G4double dist = std::max(dy,std::max(dx1,dx2)); 1077 return (dist < 0) ? -dist : 0.; 1077 return (dist < 0) ? -dist : 0.; 1078 } 1078 } 1079 case 2: // YZ section is a rectangle and 1079 case 2: // YZ section is a rectangle and 1080 { // XZ section is an isosceles tra 1080 { // XZ section is an isosceles trapezoid 1081 G4double dz = std::abs(p.z())-fDz; 1081 G4double dz = std::abs(p.z())-fDz; 1082 G4double dy = std::max(dz,std::abs(p.y( 1082 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d); 1083 G4double dx = fPlanes[3].a*std::abs(p.x 1083 G4double dx = fPlanes[3].a*std::abs(p.x()) 1084 + fPlanes[3].c*p.z()+fPlane 1084 + fPlanes[3].c*p.z()+fPlanes[3].d; 1085 G4double dist = std::max(dy,dx); 1085 G4double dist = std::max(dy,dx); 1086 return (dist < 0) ? -dist : 0.; 1086 return (dist < 0) ? -dist : 0.; 1087 } 1087 } 1088 case 3: // YZ section is a rectangle and 1088 case 3: // YZ section is a rectangle and 1089 { // XY section is an isosceles tra 1089 { // XY section is an isosceles trapezoid 1090 G4double dz = std::abs(p.z())-fDz; 1090 G4double dz = std::abs(p.z())-fDz; 1091 G4double dy = std::max(dz,std::abs(p.y( 1091 G4double dy = std::max(dz,std::abs(p.y())+fPlanes[1].d); 1092 G4double dx = fPlanes[3].a*std::abs(p.x 1092 G4double dx = fPlanes[3].a*std::abs(p.x()) 1093 + fPlanes[3].b*p.y()+fPlane 1093 + fPlanes[3].b*p.y()+fPlanes[3].d; 1094 G4double dist = std::max(dy,dx); 1094 G4double dist = std::max(dy,dx); 1095 return (dist < 0) ? -dist : 0.; 1095 return (dist < 0) ? -dist : 0.; 1096 } 1096 } 1097 } 1097 } 1098 return 0.; 1098 return 0.; 1099 } 1099 } 1100 1100 1101 ///////////////////////////////////////////// 1101 ////////////////////////////////////////////////////////////////////////// 1102 // 1102 // 1103 // GetEntityType 1103 // GetEntityType 1104 1104 1105 G4GeometryType G4Trap::GetEntityType() const 1105 G4GeometryType G4Trap::GetEntityType() const 1106 { 1106 { 1107 return {"G4Trap"}; 1107 return {"G4Trap"}; 1108 } 1108 } 1109 1109 1110 ///////////////////////////////////////////// 1110 ////////////////////////////////////////////////////////////////////////// 1111 // 1111 // 1112 // IsFaceted << 1113 << 1114 G4bool G4Trap::IsFaceted() const << 1115 { << 1116 return true; << 1117 } << 1118 << 1119 ///////////////////////////////////////////// << 1120 // << 1121 // Make a clone of the object 1112 // Make a clone of the object 1122 // 1113 // 1123 G4VSolid* G4Trap::Clone() const 1114 G4VSolid* G4Trap::Clone() const 1124 { 1115 { 1125 return new G4Trap(*this); 1116 return new G4Trap(*this); 1126 } 1117 } 1127 1118 1128 ///////////////////////////////////////////// 1119 ////////////////////////////////////////////////////////////////////////// 1129 // 1120 // 1130 // Stream object contents to an output stream 1121 // Stream object contents to an output stream 1131 1122 1132 std::ostream& G4Trap::StreamInfo( std::ostrea 1123 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const 1133 { 1124 { 1134 G4double phi = GetPhi(); 1125 G4double phi = GetPhi(); 1135 G4double theta = GetTheta(); 1126 G4double theta = GetTheta(); 1136 G4double alpha1 = GetAlpha1(); 1127 G4double alpha1 = GetAlpha1(); 1137 G4double alpha2 = GetAlpha2(); 1128 G4double alpha2 = GetAlpha2(); 1138 1129 1139 G4long oldprc = os.precision(16); 1130 G4long oldprc = os.precision(16); 1140 os << "------------------------------------ 1131 os << "-----------------------------------------------------------\n" 1141 << " *** Dump for solid: " << GetName 1132 << " *** Dump for solid: " << GetName() << " ***\n" 1142 << " ================================ 1133 << " ===================================================\n" 1143 << " Solid type: G4Trap\n" 1134 << " Solid type: G4Trap\n" 1144 << " Parameters:\n" 1135 << " Parameters:\n" 1145 << " half length Z: " << fDz/mm << " 1136 << " half length Z: " << fDz/mm << " mm\n" 1146 << " half length Y, face -Dz: " << fD 1137 << " half length Y, face -Dz: " << fDy1/mm << " mm\n" 1147 << " half length X, face -Dz, side -D 1138 << " half length X, face -Dz, side -Dy1: " << fDx1/mm << " mm\n" 1148 << " half length X, face -Dz, side +D 1139 << " half length X, face -Dz, side +Dy1: " << fDx2/mm << " mm\n" 1149 << " half length Y, face +Dz: " << fD 1140 << " half length Y, face +Dz: " << fDy2/mm << " mm\n" 1150 << " half length X, face +Dz, side -D 1141 << " half length X, face +Dz, side -Dy2: " << fDx3/mm << " mm\n" 1151 << " half length X, face +Dz, side +D 1142 << " half length X, face +Dz, side +Dy2: " << fDx4/mm << " mm\n" 1152 << " theta: " << theta/degree << " de 1143 << " theta: " << theta/degree << " degrees\n" 1153 << " phi: " << phi/degree << " degr 1144 << " phi: " << phi/degree << " degrees\n" 1154 << " alpha, face -Dz: " << alpha1/deg 1145 << " alpha, face -Dz: " << alpha1/degree << " degrees\n" 1155 << " alpha, face +Dz: " << alpha2/deg 1146 << " alpha, face +Dz: " << alpha2/degree << " degrees\n" 1156 << "------------------------------------ 1147 << "-----------------------------------------------------------\n"; 1157 os.precision(oldprc); 1148 os.precision(oldprc); 1158 1149 1159 return os; 1150 return os; 1160 } 1151 } 1161 1152 1162 ///////////////////////////////////////////// 1153 ////////////////////////////////////////////////////////////////////////// 1163 // 1154 // 1164 // Compute vertices from planes 1155 // Compute vertices from planes 1165 1156 1166 void G4Trap::GetVertices(G4ThreeVector pt[8]) 1157 void G4Trap::GetVertices(G4ThreeVector pt[8]) const 1167 { 1158 { 1168 for (G4int i=0; i<8; ++i) 1159 for (G4int i=0; i<8; ++i) 1169 { 1160 { 1170 G4int iy = (i==0 || i==1 || i==4 || i==5) 1161 G4int iy = (i==0 || i==1 || i==4 || i==5) ? 0 : 1; 1171 G4int ix = (i==0 || i==2 || i==4 || i==6) 1162 G4int ix = (i==0 || i==2 || i==4 || i==6) ? 2 : 3; 1172 G4double z = (i < 4) ? -fDz : fDz; 1163 G4double z = (i < 4) ? -fDz : fDz; 1173 G4double y = -(fPlanes[iy].c*z + fPlanes[ 1164 G4double y = -(fPlanes[iy].c*z + fPlanes[iy].d)/fPlanes[iy].b; 1174 G4double x = -(fPlanes[ix].b*y + fPlanes[ 1165 G4double x = -(fPlanes[ix].b*y + fPlanes[ix].c*z 1175 + fPlanes[ix].d)/fPlanes[i 1166 + fPlanes[ix].d)/fPlanes[ix].a; 1176 pt[i].set(x,y,z); 1167 pt[i].set(x,y,z); 1177 } 1168 } 1178 } 1169 } 1179 1170 1180 ///////////////////////////////////////////// 1171 ////////////////////////////////////////////////////////////////////////// 1181 // 1172 // 1182 // Generate random point on the surface 1173 // Generate random point on the surface 1183 1174 1184 G4ThreeVector G4Trap::GetPointOnSurface() con 1175 G4ThreeVector G4Trap::GetPointOnSurface() const 1185 { 1176 { 1186 // Set indeces 1177 // Set indeces 1187 constexpr G4int iface [6][4] = 1178 constexpr G4int iface [6][4] = 1188 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6 1179 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6,4}, {1,5,7,3}, {4,6,7,5} }; 1189 1180 1190 // Set vertices 1181 // Set vertices 1191 G4ThreeVector pt[8]; 1182 G4ThreeVector pt[8]; 1192 GetVertices(pt); 1183 GetVertices(pt); 1193 1184 1194 // Select face 1185 // Select face 1195 // 1186 // 1196 G4double select = fAreas[5]*G4QuickRand(); 1187 G4double select = fAreas[5]*G4QuickRand(); 1197 G4int k = 5; 1188 G4int k = 5; 1198 k -= (G4int)(select <= fAreas[4]); 1189 k -= (G4int)(select <= fAreas[4]); 1199 k -= (G4int)(select <= fAreas[3]); 1190 k -= (G4int)(select <= fAreas[3]); 1200 k -= (G4int)(select <= fAreas[2]); 1191 k -= (G4int)(select <= fAreas[2]); 1201 k -= (G4int)(select <= fAreas[1]); 1192 k -= (G4int)(select <= fAreas[1]); 1202 k -= (G4int)(select <= fAreas[0]); 1193 k -= (G4int)(select <= fAreas[0]); 1203 1194 1204 // Select sub-triangle 1195 // Select sub-triangle 1205 // 1196 // 1206 G4int i0 = iface[k][0]; 1197 G4int i0 = iface[k][0]; 1207 G4int i1 = iface[k][1]; 1198 G4int i1 = iface[k][1]; 1208 G4int i2 = iface[k][2]; 1199 G4int i2 = iface[k][2]; 1209 G4int i3 = iface[k][3]; 1200 G4int i3 = iface[k][3]; 1210 G4double s2 = G4GeomTools::TriangleAreaNorm 1201 G4double s2 = G4GeomTools::TriangleAreaNormal(pt[i2],pt[i1],pt[i3]).mag(); 1211 if (select > fAreas[k] - s2) i0 = i2; 1202 if (select > fAreas[k] - s2) i0 = i2; 1212 1203 1213 // Generate point 1204 // Generate point 1214 // 1205 // 1215 G4double u = G4QuickRand(); 1206 G4double u = G4QuickRand(); 1216 G4double v = G4QuickRand(); 1207 G4double v = G4QuickRand(); 1217 if (u + v > 1.) { u = 1. - u; v = 1. - v; } 1208 if (u + v > 1.) { u = 1. - u; v = 1. - v; } 1218 return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3 1209 return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3]; 1219 } 1210 } 1220 1211 1221 ///////////////////////////////////////////// 1212 ////////////////////////////////////////////////////////////////////////// 1222 // 1213 // 1223 // Methods for visualisation 1214 // Methods for visualisation 1224 1215 1225 void G4Trap::DescribeYourselfTo ( G4VGraphics 1216 void G4Trap::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 1226 { 1217 { 1227 scene.AddSolid (*this); 1218 scene.AddSolid (*this); 1228 } 1219 } 1229 1220 1230 G4Polyhedron* G4Trap::CreatePolyhedron () con 1221 G4Polyhedron* G4Trap::CreatePolyhedron () const 1231 { 1222 { 1232 G4double phi = std::atan2(fTthetaSphi, fTth 1223 G4double phi = std::atan2(fTthetaSphi, fTthetaCphi); 1233 G4double alpha1 = std::atan(fTalpha1); 1224 G4double alpha1 = std::atan(fTalpha1); 1234 G4double alpha2 = std::atan(fTalpha2); 1225 G4double alpha2 = std::atan(fTalpha2); 1235 G4double theta = std::atan(std::sqrt(fTthet 1226 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi 1236 +fTthet 1227 +fTthetaSphi*fTthetaSphi)); 1237 1228 1238 return new G4PolyhedronTrap(fDz, theta, phi 1229 return new G4PolyhedronTrap(fDz, theta, phi, 1239 fDy1, fDx1, fDx 1230 fDy1, fDx1, fDx2, alpha1, 1240 fDy2, fDx3, fDx 1231 fDy2, fDx3, fDx4, alpha2); 1241 } 1232 } 1242 1233 1243 #endif 1234 #endif 1244 1235