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 // G4VTwistSurface implementation 26 // G4VTwistSurface implementation 27 // 27 // 28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu 28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp), created. 29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), 29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4 30 // from original version in Jupi 30 // from original version in Jupiter-2.5.02 application. 31 // ------------------------------------------- 31 // -------------------------------------------------------------------- 32 32 33 #include <iomanip> 33 #include <iomanip> 34 34 35 #include "G4VTwistSurface.hh" 35 #include "G4VTwistSurface.hh" 36 #include "G4GeometryTolerance.hh" 36 #include "G4GeometryTolerance.hh" 37 37 38 const G4int G4VTwistSurface::sOutside 38 const G4int G4VTwistSurface::sOutside = 0x00000000; 39 const G4int G4VTwistSurface::sInside 39 const G4int G4VTwistSurface::sInside = 0x10000000; 40 const G4int G4VTwistSurface::sBoundary 40 const G4int G4VTwistSurface::sBoundary = 0x20000000; 41 const G4int G4VTwistSurface::sCorner 41 const G4int G4VTwistSurface::sCorner = 0x40000000; 42 const G4int G4VTwistSurface::sC0Min1Min 42 const G4int G4VTwistSurface::sC0Min1Min = 0x40000101; 43 const G4int G4VTwistSurface::sC0Max1Min 43 const G4int G4VTwistSurface::sC0Max1Min = 0x40000201; 44 const G4int G4VTwistSurface::sC0Max1Max 44 const G4int G4VTwistSurface::sC0Max1Max = 0x40000202; 45 const G4int G4VTwistSurface::sC0Min1Max 45 const G4int G4VTwistSurface::sC0Min1Max = 0x40000102; 46 const G4int G4VTwistSurface::sAxisMin 46 const G4int G4VTwistSurface::sAxisMin = 0x00000101; 47 const G4int G4VTwistSurface::sAxisMax 47 const G4int G4VTwistSurface::sAxisMax = 0x00000202; 48 const G4int G4VTwistSurface::sAxisX 48 const G4int G4VTwistSurface::sAxisX = 0x00000404; 49 const G4int G4VTwistSurface::sAxisY 49 const G4int G4VTwistSurface::sAxisY = 0x00000808; 50 const G4int G4VTwistSurface::sAxisZ 50 const G4int G4VTwistSurface::sAxisZ = 0x00000C0C; 51 const G4int G4VTwistSurface::sAxisRho 51 const G4int G4VTwistSurface::sAxisRho = 0x00001010; 52 const G4int G4VTwistSurface::sAxisPhi 52 const G4int G4VTwistSurface::sAxisPhi = 0x00001414; 53 53 54 // mask 54 // mask 55 const G4int G4VTwistSurface::sAxis0 55 const G4int G4VTwistSurface::sAxis0 = 0x0000FF00; 56 const G4int G4VTwistSurface::sAxis1 56 const G4int G4VTwistSurface::sAxis1 = 0x000000FF; 57 const G4int G4VTwistSurface::sSizeMask 57 const G4int G4VTwistSurface::sSizeMask = 0x00000303; 58 const G4int G4VTwistSurface::sAxisMask 58 const G4int G4VTwistSurface::sAxisMask = 0x0000FCFC; 59 const G4int G4VTwistSurface::sAreaMask 59 const G4int G4VTwistSurface::sAreaMask = 0XF0000000; 60 60 61 //============================================ 61 //===================================================================== 62 //* constructors ----------------------------- 62 //* constructors ------------------------------------------------------ 63 63 64 G4VTwistSurface::G4VTwistSurface(const G4Strin 64 G4VTwistSurface::G4VTwistSurface(const G4String &name) 65 : fIsValidNorm(false), fName(name) 65 : fIsValidNorm(false), fName(name) 66 { 66 { 67 67 68 fAxis[0] = kUndefined; 68 fAxis[0] = kUndefined; 69 fAxis[1] = kUndefined; 69 fAxis[1] = kUndefined; 70 fAxisMin[0] = kInfinity; 70 fAxisMin[0] = kInfinity; 71 fAxisMin[1] = kInfinity; 71 fAxisMin[1] = kInfinity; 72 fAxisMax[0] = kInfinity; 72 fAxisMax[0] = kInfinity; 73 fAxisMax[1] = kInfinity; 73 fAxisMax[1] = kInfinity; 74 fHandedness = 1; 74 fHandedness = 1; 75 75 76 for (auto i=0; i<4; ++i) 76 for (auto i=0; i<4; ++i) 77 { 77 { 78 fCorners[i].set(kInfinity, kInfinity, kI 78 fCorners[i].set(kInfinity, kInfinity, kInfinity); 79 fNeighbours[i] = nullptr; 79 fNeighbours[i] = nullptr; 80 } 80 } 81 81 82 fCurrentNormal.p.set(kInfinity, kInfinity, 82 fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity); 83 83 84 fAmIOnLeftSide.me.set(kInfinity, kInfinity, 84 fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity); 85 fAmIOnLeftSide.vec.set(kInfinity, kInfinity 85 fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity); 86 kCarTolerance = G4GeometryTolerance::GetIns 86 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 87 } 87 } 88 88 89 G4VTwistSurface::G4VTwistSurface(const G4Strin 89 G4VTwistSurface::G4VTwistSurface(const G4String& name, 90 const G4RotationMatrix& 90 const G4RotationMatrix& rot, 91 const G4ThreeVector& 91 const G4ThreeVector& tlate, 92 G4int 92 G4int handedness, 93 const EAxis 93 const EAxis axis0 , 94 const EAxis 94 const EAxis axis1 , 95 G4double 95 G4double axis0min, 96 G4double 96 G4double axis1min, 97 G4double 97 G4double axis0max, 98 G4double 98 G4double axis1max ) 99 : fIsValidNorm(false), fName(name) 99 : fIsValidNorm(false), fName(name) 100 { 100 { 101 fAxis[0] = axis0; 101 fAxis[0] = axis0; 102 fAxis[1] = axis1; 102 fAxis[1] = axis1; 103 fAxisMin[0] = axis0min; 103 fAxisMin[0] = axis0min; 104 fAxisMin[1] = axis1min; 104 fAxisMin[1] = axis1min; 105 fAxisMax[0] = axis0max; 105 fAxisMax[0] = axis0max; 106 fAxisMax[1] = axis1max; 106 fAxisMax[1] = axis1max; 107 fHandedness = handedness; 107 fHandedness = handedness; 108 fRot = rot; 108 fRot = rot; 109 fTrans = tlate; 109 fTrans = tlate; 110 110 111 for (auto i=0; i<4; ++i) 111 for (auto i=0; i<4; ++i) 112 { 112 { 113 fCorners[i].set(kInfinity, kInfinity, kI 113 fCorners[i].set(kInfinity, kInfinity, kInfinity); 114 fNeighbours[i] = nullptr; 114 fNeighbours[i] = nullptr; 115 } 115 } 116 116 117 fCurrentNormal.p.set(kInfinity, kInfinity, 117 fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity); 118 118 119 fAmIOnLeftSide.me.set(kInfinity, kInfinity, 119 fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity); 120 fAmIOnLeftSide.vec.set(kInfinity, kInfinity 120 fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity); 121 kCarTolerance = G4GeometryTolerance::GetIns 121 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 122 } 122 } 123 123 124 //============================================ 124 //===================================================================== 125 //* Fake default constructor ----------------- 125 //* Fake default constructor ------------------------------------------ 126 126 127 G4VTwistSurface::G4VTwistSurface( __void__& ) 127 G4VTwistSurface::G4VTwistSurface( __void__& ) 128 : fHandedness(0), fIsValidNorm(false), kCarT 128 : fHandedness(0), fIsValidNorm(false), kCarTolerance(0.), 129 fName("") 129 fName("") 130 { 130 { 131 fAxis[0] = fAxis[1] = kXAxis; 131 fAxis[0] = fAxis[1] = kXAxis; 132 fAxisMin[0] = fAxisMin[1] = 0.; 132 fAxisMin[0] = fAxisMin[1] = 0.; 133 fAxisMax[0] = fAxisMax[1] = 0.; 133 fAxisMax[0] = fAxisMax[1] = 0.; 134 fNeighbours[0] = fNeighbours[1] = fNeighbou 134 fNeighbours[0] = fNeighbours[1] = fNeighbours[2] = fNeighbours[3] = nullptr; 135 } 135 } 136 136 137 //============================================ 137 //===================================================================== >> 138 //* destructor -------------------------------------------------------- >> 139 >> 140 G4VTwistSurface::~G4VTwistSurface() >> 141 { >> 142 } >> 143 >> 144 //===================================================================== 138 //* AmIOnLeftSide ---------------------------- 145 //* AmIOnLeftSide ----------------------------------------------------- 139 146 140 G4int G4VTwistSurface::AmIOnLeftSide(const G4T 147 G4int G4VTwistSurface::AmIOnLeftSide(const G4ThreeVector& me, 141 const G4T 148 const G4ThreeVector& vec, 142 G4b 149 G4bool withtol) 143 { 150 { 144 // AmIOnLeftSide returns phi-location of "m 151 // AmIOnLeftSide returns phi-location of "me" 145 // (phi relation between me and vec project 152 // (phi relation between me and vec projected on z=0 plane). 146 // If "me" is on -ve-phi-side of "vec", it 153 // If "me" is on -ve-phi-side of "vec", it returns 1. 147 // On the other hand, if "me" is on +ve-phi 154 // On the other hand, if "me" is on +ve-phi-side of "vec", 148 // it returns -1. 155 // it returns -1. 149 // (The return value represents z-coordinat 156 // (The return value represents z-coordinate of normal vector 150 // of me.cross(vec).) 157 // of me.cross(vec).) 151 // If me is on boundary of vec, return 0. 158 // If me is on boundary of vec, return 0. 152 159 153 const G4double kAngTolerance 160 const G4double kAngTolerance 154 = G4GeometryTolerance::GetInstance()->Get 161 = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 155 162 156 G4RotationMatrix unitrot; 163 G4RotationMatrix unitrot; 157 const G4RotationMatrix rottol = unitrot. 164 const G4RotationMatrix rottol = unitrot.rotateZ(0.5*kAngTolerance); 158 const G4RotationMatrix invrottol = unitrot. 165 const G4RotationMatrix invrottol = unitrot.rotateZ(-1.*kAngTolerance); 159 166 160 if (fAmIOnLeftSide.me == me 167 if (fAmIOnLeftSide.me == me 161 && fAmIOnLeftSide.vec == vec 168 && fAmIOnLeftSide.vec == vec 162 && fAmIOnLeftSide.withTol == withtol) 169 && fAmIOnLeftSide.withTol == withtol) 163 { 170 { 164 return fAmIOnLeftSide.amIOnLeftSide; 171 return fAmIOnLeftSide.amIOnLeftSide; 165 } 172 } 166 173 167 fAmIOnLeftSide.me = me; 174 fAmIOnLeftSide.me = me; 168 fAmIOnLeftSide.vec = vec; 175 fAmIOnLeftSide.vec = vec; 169 fAmIOnLeftSide.withTol = withtol; 176 fAmIOnLeftSide.withTol = withtol; 170 177 171 G4ThreeVector met = (G4ThreeVector(me.x() 178 G4ThreeVector met = (G4ThreeVector(me.x(), me.y(), 0.)).unit(); 172 G4ThreeVector vect = (G4ThreeVector(vec.x( 179 G4ThreeVector vect = (G4ThreeVector(vec.x(), vec.y(), 0.)).unit(); 173 180 174 G4ThreeVector ivect = invrottol * vect; 181 G4ThreeVector ivect = invrottol * vect; 175 G4ThreeVector rvect = rottol * vect; 182 G4ThreeVector rvect = rottol * vect; 176 183 177 G4double metcrossvect = met.x() * vect.y() 184 G4double metcrossvect = met.x() * vect.y() - met.y() * vect.x(); 178 185 179 if (withtol) 186 if (withtol) 180 { 187 { 181 if (met.x() * ivect.y() - met.y() * ivec 188 if (met.x() * ivect.y() - met.y() * ivect.x() > 0 && 182 metcrossvect >= 0) { 189 metcrossvect >= 0) { 183 fAmIOnLeftSide.amIOnLeftSide = 1; 190 fAmIOnLeftSide.amIOnLeftSide = 1; 184 } else if (met.x() * rvect.y() - met.y() 191 } else if (met.x() * rvect.y() - met.y() * rvect.x() < 0 && 185 metcrossvect <= 0) { 192 metcrossvect <= 0) { 186 fAmIOnLeftSide.amIOnLeftSide = -1; 193 fAmIOnLeftSide.amIOnLeftSide = -1; 187 } else { 194 } else { 188 fAmIOnLeftSide.amIOnLeftSide = 0; 195 fAmIOnLeftSide.amIOnLeftSide = 0; 189 } 196 } 190 } 197 } 191 else 198 else 192 { 199 { 193 if (metcrossvect > 0) { 200 if (metcrossvect > 0) { 194 fAmIOnLeftSide.amIOnLeftSide = 1; 201 fAmIOnLeftSide.amIOnLeftSide = 1; 195 } else if (metcrossvect < 0 ) { 202 } else if (metcrossvect < 0 ) { 196 fAmIOnLeftSide.amIOnLeftSide = -1; 203 fAmIOnLeftSide.amIOnLeftSide = -1; 197 } else { 204 } else { 198 fAmIOnLeftSide.amIOnLeftSide = 0; 205 fAmIOnLeftSide.amIOnLeftSide = 0; 199 } 206 } 200 } 207 } 201 208 202 #ifdef G4TWISTDEBUG 209 #ifdef G4TWISTDEBUG 203 G4cout << " === G4VTwistSurface::Am 210 G4cout << " === G4VTwistSurface::AmIOnLeftSide() ==============" 204 << G4endl; 211 << G4endl; 205 G4cout << " Name , returncode 212 G4cout << " Name , returncode : " << fName << " " 206 << fAmIOnLeftSide.amIOn 213 << fAmIOnLeftSide.amIOnLeftSide << G4endl; 207 G4cout << " me, vec : " << s 214 G4cout << " me, vec : " << std::setprecision(14) << me 208 << " 215 << " " << vec << G4endl; 209 G4cout << " met, vect : " << m 216 G4cout << " met, vect : " << met << " " << vect << G4endl; 210 G4cout << " ivec, rvec : " << i 217 G4cout << " ivec, rvec : " << ivect << " " << rvect << G4endl; 211 G4cout << " met x vect : " << m 218 G4cout << " met x vect : " << metcrossvect << G4endl; 212 G4cout << " met x ivec : " << m 219 G4cout << " met x ivec : " << met.cross(ivect) << G4endl; 213 G4cout << " met x rvec : " << m 220 G4cout << " met x rvec : " << met.cross(rvect) << G4endl; 214 G4cout << " ======================= 221 G4cout << " ==============================================" 215 << G4endl; 222 << G4endl; 216 #endif 223 #endif 217 224 218 return fAmIOnLeftSide.amIOnLeftSide; 225 return fAmIOnLeftSide.amIOnLeftSide; 219 } 226 } 220 227 221 //============================================ 228 //===================================================================== 222 //* DistanceToBoundary ----------------------- 229 //* DistanceToBoundary ------------------------------------------------ 223 230 224 G4double G4VTwistSurface::DistanceToBoundary(G 231 G4double G4VTwistSurface::DistanceToBoundary(G4int areacode, 225 G 232 G4ThreeVector& xx, 226 const G 233 const G4ThreeVector& p) 227 { 234 { 228 // DistanceToBoundary 235 // DistanceToBoundary 229 // 236 // 230 // return distance to nearest boundary from 237 // return distance to nearest boundary from arbitrary point p 231 // in local coodinate. 238 // in local coodinate. 232 // Argument areacode must be one of them: 239 // Argument areacode must be one of them: 233 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 240 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 234 // sAxis1 & sAxisMin, sAxis1 & sAxisMax. 241 // sAxis1 & sAxisMin, sAxis1 & sAxisMax. 235 // 242 // 236 243 237 G4ThreeVector d; // direction vector of 244 G4ThreeVector d; // direction vector of the boundary 238 G4ThreeVector x0; // reference point of t 245 G4ThreeVector x0; // reference point of the boundary 239 G4double dist = kInfinity; 246 G4double dist = kInfinity; 240 G4int boundarytype; 247 G4int boundarytype; 241 248 242 if (IsAxis0(areacode) && IsAxis1(areacode)) 249 if (IsAxis0(areacode) && IsAxis1(areacode)) 243 { 250 { 244 std::ostringstream message; 251 std::ostringstream message; 245 message << "Point is in the corner area. 252 message << "Point is in the corner area." << G4endl 246 << " Point is in the corn 253 << " Point is in the corner area. This function returns" 247 << G4endl 254 << G4endl 248 << " a direction vector o 255 << " a direction vector of a boundary line." << G4endl 249 << " areacode = " << area 256 << " areacode = " << areacode; 250 G4Exception("G4VTwistSurface::DistanceTo 257 G4Exception("G4VTwistSurface::DistanceToBoundary()", "GeomSolids0003", 251 FatalException, message); 258 FatalException, message); 252 } 259 } 253 else if (IsAxis0(areacode) || IsAxis1(areac 260 else if (IsAxis0(areacode) || IsAxis1(areacode)) 254 { 261 { 255 GetBoundaryParameters(areacode, d, x0, b 262 GetBoundaryParameters(areacode, d, x0, boundarytype); 256 if (boundarytype == sAxisPhi) 263 if (boundarytype == sAxisPhi) 257 { 264 { 258 G4double t = x0.getRho() / p.getRho() 265 G4double t = x0.getRho() / p.getRho(); 259 xx.set(t*p.x(), t*p.y(), x0.z()); 266 xx.set(t*p.x(), t*p.y(), x0.z()); 260 dist = (xx - p).mag(); 267 dist = (xx - p).mag(); 261 } 268 } 262 else 269 else 263 { 270 { 264 // linear boundary 271 // linear boundary 265 // sAxisX, sAxisY, sAxisZ, sAxisRho 272 // sAxisX, sAxisY, sAxisZ, sAxisRho 266 dist = DistanceToLine(p, x0, d, xx); 273 dist = DistanceToLine(p, x0, d, xx); 267 } 274 } 268 } 275 } 269 else 276 else 270 { 277 { 271 std::ostringstream message; 278 std::ostringstream message; 272 message << "Bad areacode of boundary." < 279 message << "Bad areacode of boundary." << G4endl 273 << " areacode = " << area 280 << " areacode = " << areacode; 274 G4Exception("G4VTwistSurface::DistanceTo 281 G4Exception("G4VTwistSurface::DistanceToBoundary()", "GeomSolids0003", 275 FatalException, message); 282 FatalException, message); 276 } 283 } 277 return dist; 284 return dist; 278 } 285 } 279 286 280 //============================================ 287 //===================================================================== 281 //* DistanceToIn ----------------------------- 288 //* DistanceToIn ------------------------------------------------------ 282 289 283 G4double G4VTwistSurface::DistanceToIn(const G 290 G4double G4VTwistSurface::DistanceToIn(const G4ThreeVector& gp, 284 const G 291 const G4ThreeVector& gv, 285 G 292 G4ThreeVector& gxxbest) 286 { 293 { 287 #ifdef G4TWISTDEBUG 294 #ifdef G4TWISTDEBUG 288 G4cout << " ~~~~ G4VTwistSurface::DistanceT 295 G4cout << " ~~~~ G4VTwistSurface::DistanceToIn(p,v) - Start ~~~~" << G4endl; 289 G4cout << " Name : " << fName << G4end 296 G4cout << " Name : " << fName << G4endl; 290 G4cout << " gp : " << gp << G4endl; 297 G4cout << " gp : " << gp << G4endl; 291 G4cout << " gv : " << gv << G4endl; 298 G4cout << " gv : " << gv << G4endl; 292 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 299 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; 293 #endif 300 #endif 294 301 295 G4ThreeVector gxx[G4VSURFACENXX]; 302 G4ThreeVector gxx[G4VSURFACENXX]; 296 G4double distance[G4VSURFACENXX] ; 303 G4double distance[G4VSURFACENXX] ; 297 G4int areacode[G4VSURFACENXX] ; 304 G4int areacode[G4VSURFACENXX] ; 298 G4bool isvalid[G4VSURFACENXX] ; 305 G4bool isvalid[G4VSURFACENXX] ; 299 306 300 for (G4int i = 0 ; i<G4VSURFACENXX ; ++i ) 307 for (G4int i = 0 ; i<G4VSURFACENXX ; ++i ) 301 { 308 { 302 distance[i] = kInfinity ; 309 distance[i] = kInfinity ; 303 areacode[i] = sOutside ; 310 areacode[i] = sOutside ; 304 isvalid[i] = false ; 311 isvalid[i] = false ; 305 } 312 } 306 313 307 G4double bestdistance = kInfinity; 314 G4double bestdistance = kInfinity; 308 #ifdef G4TWISTDEBUG 315 #ifdef G4TWISTDEBUG 309 G4int besti = -1; 316 G4int besti = -1; 310 #endif 317 #endif 311 G4ThreeVector bestgxx(kInfinity, kInfinity, 318 G4ThreeVector bestgxx(kInfinity, kInfinity, kInfinity); 312 319 313 G4int nxx = DistanceToSurface(gp, g 320 G4int nxx = DistanceToSurface(gp, gv, gxx, distance, areacode, 314 isval 321 isvalid, kValidateWithTol); 315 322 316 for (G4int i=0; i<nxx; ++i) 323 for (G4int i=0; i<nxx; ++i) 317 { 324 { 318 325 319 // skip this intersection if: 326 // skip this intersection if: 320 // - invalid intersection 327 // - invalid intersection 321 // - particle goes outword the surface 328 // - particle goes outword the surface 322 329 323 if (!isvalid[i]) 330 if (!isvalid[i]) 324 { 331 { 325 // xx[i] is sOutside or distance[i] < 332 // xx[i] is sOutside or distance[i] < 0 326 continue; 333 continue; 327 } 334 } 328 335 329 G4ThreeVector normal = GetNormal(gxx[i], 336 G4ThreeVector normal = GetNormal(gxx[i], true); 330 337 331 if ((normal * gv) >= 0) 338 if ((normal * gv) >= 0) 332 { 339 { 333 340 334 #ifdef G4TWISTDEBUG 341 #ifdef G4TWISTDEBUG 335 G4cout << " G4VTwistSurface::Distan 342 G4cout << " G4VTwistSurface::DistanceToIn(p,v): " 336 << "particle goes outword the 343 << "particle goes outword the surface." << G4endl; 337 #endif 344 #endif 338 continue; 345 continue; 339 } 346 } 340 347 341 // 348 // 342 // accept this intersection if the inter 349 // accept this intersection if the intersection is inside. 343 // 350 // 344 351 345 if (IsInside(areacode[i])) 352 if (IsInside(areacode[i])) 346 { 353 { 347 if (distance[i] < bestdistance) 354 if (distance[i] < bestdistance) 348 { 355 { 349 bestdistance = distance[i]; 356 bestdistance = distance[i]; 350 bestgxx = gxx[i]; 357 bestgxx = gxx[i]; 351 #ifdef G4TWISTDEBUG 358 #ifdef G4TWISTDEBUG 352 besti = i; 359 besti = i; 353 G4cout << " G4VTwistSurface::Dis 360 G4cout << " G4VTwistSurface::DistanceToIn(p,v): " 354 << " areacode sInside name, 361 << " areacode sInside name, distance = " 355 << fName << " "<< bestdist 362 << fName << " "<< bestdistance << G4endl; 356 #endif 363 #endif 357 } 364 } 358 365 359 // 366 // 360 // else, the intersection is on boundary 367 // else, the intersection is on boundary or corner. 361 // 368 // 362 369 363 } 370 } 364 else 371 else 365 { 372 { 366 G4VTwistSurface* neighbours[2]; 373 G4VTwistSurface* neighbours[2]; 367 G4bool isaccepted[2] = {false, f 374 G4bool isaccepted[2] = {false, false}; 368 G4int nneighbours = GetNeighb 375 G4int nneighbours = GetNeighbours(areacode[i], neighbours); 369 376 370 for (G4int j=0; j<nneighbours; ++j) 377 for (G4int j=0; j<nneighbours; ++j) 371 { 378 { 372 // if on corner, nneighbours = 2. 379 // if on corner, nneighbours = 2. 373 // if on boundary, nneighbours = 1 380 // if on boundary, nneighbours = 1. 374 381 375 G4ThreeVector tmpgxx[G4VSURFACENXX 382 G4ThreeVector tmpgxx[G4VSURFACENXX]; 376 G4double tmpdist[G4VSURFACENX 383 G4double tmpdist[G4VSURFACENXX] ; 377 G4int tmpareacode[G4VSURFA 384 G4int tmpareacode[G4VSURFACENXX] ; 378 G4bool tmpisvalid[G4VSURFAC 385 G4bool tmpisvalid[G4VSURFACENXX] ; 379 386 380 for (G4int l = 0 ; l<G4VSURFACENXX 387 for (G4int l = 0 ; l<G4VSURFACENXX ; ++l ) 381 { 388 { 382 tmpdist[l] = kInfinity ; 389 tmpdist[l] = kInfinity ; 383 tmpareacode[l] = sOutside ; 390 tmpareacode[l] = sOutside ; 384 tmpisvalid[l] = false ; 391 tmpisvalid[l] = false ; 385 } 392 } 386 393 387 G4int tmpnxx = neighbours[j]->Dist 394 G4int tmpnxx = neighbours[j]->DistanceToSurface( 388 gp, 395 gp, gv, tmpgxx, tmpdist, 389 tmpa 396 tmpareacode, tmpisvalid, 390 kVal 397 kValidateWithTol); 391 G4ThreeVector neighbournormal; 398 G4ThreeVector neighbournormal; 392 399 393 for (G4int k=0; k< tmpnxx; ++k) 400 for (G4int k=0; k< tmpnxx; ++k) 394 { 401 { 395 // 402 // 396 // if tmpxx[k] is valid && sIns 403 // if tmpxx[k] is valid && sInside, the final winner must 397 // be neighbour surface. return 404 // be neighbour surface. return kInfinity. 398 // else , choose tmpxx on same 405 // else , choose tmpxx on same boundary of xx, then check normal 399 // 406 // 400 407 401 if (IsInside(tmpareacode[k])) 408 if (IsInside(tmpareacode[k])) 402 { 409 { 403 #ifdef G4TWISTDEBUG 410 #ifdef G4TWISTDEBUG 404 G4cout << " G4VTwistSurfac 411 G4cout << " G4VTwistSurface:DistanceToIn(p,v): " 405 << " intersection "<< 412 << " intersection "<< tmpgxx[k] << G4endl 406 << " is inside of n 413 << " is inside of neighbour surface of " << fName 407 << " . returning kInf 414 << " . returning kInfinity." << G4endl; 408 G4cout << "~~ G4VTwistSurfac 415 G4cout << "~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~" 409 << G4endl; 416 << G4endl; 410 G4cout << " No intersec 417 G4cout << " No intersections " << G4endl; 411 G4cout << " Name : " << 418 G4cout << " Name : " << fName << G4endl; 412 G4cout << "~~~~~~~~~~~~~~~~~ 419 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" 413 << G4endl; 420 << G4endl; 414 #endif 421 #endif 415 if (tmpisvalid[k]) return k 422 if (tmpisvalid[k]) return kInfinity; 416 continue; 423 continue; 417 424 418 // 425 // 419 // if tmpxx[k] is valid && sIns 426 // if tmpxx[k] is valid && sInside, the final winner must 420 // be neighbour surface. return 427 // be neighbour surface. return . 421 // 428 // 422 429 423 } 430 } 424 else if (IsSameBoundary(this,ar 431 else if (IsSameBoundary(this,areacode[i], 425 neighbo 432 neighbours[j], tmpareacode[k])) 426 { 433 { 427 // tmpxx[k] is same boundary 434 // tmpxx[k] is same boundary (or corner) of xx. 428 435 429 neighbournormal = neighbours 436 neighbournormal = neighbours[j]->GetNormal(tmpgxx[k], true); 430 if (neighbournormal * gv < 0 437 if (neighbournormal * gv < 0) isaccepted[j] = true; 431 } 438 } 432 } 439 } 433 440 434 // if nneighbours = 1, chabge isac 441 // if nneighbours = 1, chabge isaccepted[1] before 435 // exiting neighboursurface loop. 442 // exiting neighboursurface loop. 436 443 437 if (nneighbours == 1) isaccepted[1 444 if (nneighbours == 1) isaccepted[1] = true; 438 445 439 } // neighboursurface loop end 446 } // neighboursurface loop end 440 447 441 // now, we can accept xx intersection 448 // now, we can accept xx intersection 442 449 443 if (isaccepted[0] && isaccepted[1]) << 450 if (isaccepted[0] == true && isaccepted[1] == true) 444 { 451 { 445 if (distance[i] < bestdistance) 452 if (distance[i] < bestdistance) 446 { 453 { 447 bestdistance = distance[i]; 454 bestdistance = distance[i]; 448 gxxbest = gxx[i]; 455 gxxbest = gxx[i]; 449 #ifdef G4TWISTDEBUG 456 #ifdef G4TWISTDEBUG 450 besti = i; 457 besti = i; 451 G4cout << " G4VTwistSurface:: 458 G4cout << " G4VTwistSurface::DistanceToIn(p,v): " 452 << " areacode sBoundary 459 << " areacode sBoundary & sBoundary distance = " 453 << fName << " " << dist 460 << fName << " " << distance[i] << G4endl; 454 #endif 461 #endif 455 } 462 } 456 } 463 } 457 } // else end 464 } // else end 458 } // intersection loop end 465 } // intersection loop end 459 466 460 gxxbest = bestgxx; 467 gxxbest = bestgxx; 461 468 462 #ifdef G4TWISTDEBUG 469 #ifdef G4TWISTDEBUG 463 if (besti < 0) 470 if (besti < 0) 464 { 471 { 465 G4cout << "~~~ G4VTwistSurface::Distance 472 G4cout << "~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~" << G4endl; 466 G4cout << " No intersections " << G 473 G4cout << " No intersections " << G4endl; 467 G4cout << " Name : " << fName << G4 474 G4cout << " Name : " << fName << G4endl; 468 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 475 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; 469 } 476 } 470 else 477 else 471 { 478 { 472 G4cout << "~~~ G4VTwistSurface::Distance 479 G4cout << "~~~ G4VTwistSurface::DistanceToIn(p,v) : return ~~~" << G4endl; 473 G4cout << " Name, i : " << fName < 480 G4cout << " Name, i : " << fName << " , " << besti << G4endl; 474 G4cout << " gxx[i] : " << gxxbest 481 G4cout << " gxx[i] : " << gxxbest << G4endl; 475 G4cout << " bestdist : " << bestdis 482 G4cout << " bestdist : " << bestdistance << G4endl; 476 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 483 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; 477 } 484 } 478 485 479 #endif 486 #endif 480 487 481 return bestdistance; 488 return bestdistance; 482 } 489 } 483 490 484 //============================================ 491 //===================================================================== 485 //* DistanceToOut(p, v) ---------------------- 492 //* DistanceToOut(p, v) ----------------------------------------------- 486 493 487 G4double G4VTwistSurface::DistanceToOut(const 494 G4double G4VTwistSurface::DistanceToOut(const G4ThreeVector& gp, 488 const 495 const G4ThreeVector& gv, 489 496 G4ThreeVector& gxxbest) 490 { 497 { 491 #ifdef G4TWISTDEBUG 498 #ifdef G4TWISTDEBUG 492 G4cout << "~~~~~ G4VTwistSurface::DistanceT 499 G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) - Start ~~~~" << G4endl; 493 G4cout << " Name : " << fName << G4end 500 G4cout << " Name : " << fName << G4endl; 494 G4cout << " gp : " << gp << G4endl; 501 G4cout << " gp : " << gp << G4endl; 495 G4cout << " gv : " << gv << G4endl; 502 G4cout << " gv : " << gv << G4endl; 496 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 503 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; 497 #endif 504 #endif 498 505 499 G4ThreeVector gxx[G4VSURFACENXX]; 506 G4ThreeVector gxx[G4VSURFACENXX]; 500 G4double distance[G4VSURFACENXX]; 507 G4double distance[G4VSURFACENXX]; 501 G4int areacode[G4VSURFACENXX]; 508 G4int areacode[G4VSURFACENXX]; 502 G4bool isvalid[G4VSURFACENXX]; 509 G4bool isvalid[G4VSURFACENXX]; 503 510 504 for ( G4int i = 0 ; i<G4VSURFACENXX ; ++i ) 511 for ( G4int i = 0 ; i<G4VSURFACENXX ; ++i ) 505 { 512 { 506 distance[i] = kInfinity ; 513 distance[i] = kInfinity ; 507 areacode[i] = sOutside ; 514 areacode[i] = sOutside ; 508 isvalid[i] = false ; 515 isvalid[i] = false ; 509 } 516 } 510 517 511 G4int nxx; 518 G4int nxx; 512 G4double bestdistance = kInfinity; 519 G4double bestdistance = kInfinity; 513 520 514 nxx = DistanceToSurface(gp, gv, gxx, distan 521 nxx = DistanceToSurface(gp, gv, gxx, distance, areacode, 515 isvalid, kValidateW 522 isvalid, kValidateWithTol); 516 523 517 for (G4int i=0; i<nxx; ++i) 524 for (G4int i=0; i<nxx; ++i) 518 { 525 { 519 if (!(isvalid[i])) 526 if (!(isvalid[i])) 520 { 527 { 521 continue; 528 continue; 522 } 529 } 523 530 524 G4ThreeVector normal = GetNormal(gxx[i], 531 G4ThreeVector normal = GetNormal(gxx[i], true); 525 if (normal * gv <= 0) 532 if (normal * gv <= 0) 526 { 533 { 527 // particle goes toword inside of sol 534 // particle goes toword inside of solid, return kInfinity 528 #ifdef G4TWISTDEBUG 535 #ifdef G4TWISTDEBUG 529 G4cout << " G4VTwistSurface::Dista 536 G4cout << " G4VTwistSurface::DistanceToOut(p,v): normal*gv < 0 " 530 << fName << " " << normal 537 << fName << " " << normal 531 << G4endl; 538 << G4endl; 532 #endif 539 #endif 533 } 540 } 534 else 541 else 535 { 542 { 536 // gxx[i] is accepted. 543 // gxx[i] is accepted. 537 if (distance[i] < bestdistance) 544 if (distance[i] < bestdistance) 538 { 545 { 539 bestdistance = distance[i]; 546 bestdistance = distance[i]; 540 gxxbest = gxx[i]; 547 gxxbest = gxx[i]; 541 } 548 } 542 } 549 } 543 } 550 } 544 551 545 #ifdef G4TWISTDEBUG 552 #ifdef G4TWISTDEBUG 546 if (besti < 0) 553 if (besti < 0) 547 { 554 { 548 G4cout << "~~ G4VTwistSurface::DistanceT 555 G4cout << "~~ G4VTwistSurface::DistanceToOut(p,v) - return ~~" << G4endl; 549 G4cout << " No intersections " << 556 G4cout << " No intersections " << G4endl; 550 G4cout << " Name : " << fName < 557 G4cout << " Name : " << fName << G4endl; 551 G4cout << " bestdist : " << bestdis 558 G4cout << " bestdist : " << bestdistance << G4endl; 552 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 559 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; 553 } 560 } 554 else 561 else 555 { 562 { 556 G4cout << "~~ G4VTwistSurface::DistanceT 563 G4cout << "~~ G4VTwistSurface::DistanceToOut(p,v) : return ~~" << G4endl; 557 G4cout << " Name, i : " << fName < 564 G4cout << " Name, i : " << fName << " , " << i << G4endl; 558 G4cout << " gxx[i] : " << gxxbest 565 G4cout << " gxx[i] : " << gxxbest << G4endl; 559 G4cout << " bestdist : " << bestdis 566 G4cout << " bestdist : " << bestdistance << G4endl; 560 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 567 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; 561 } 568 } 562 #endif 569 #endif 563 570 564 return bestdistance; 571 return bestdistance; 565 } 572 } 566 573 567 //============================================ 574 //===================================================================== 568 //* DistanceTo(p) ---------------------------- 575 //* DistanceTo(p) ----------------------------------------------------- 569 576 570 G4double G4VTwistSurface::DistanceTo(const G4T 577 G4double G4VTwistSurface::DistanceTo(const G4ThreeVector& gp, 571 G4T 578 G4ThreeVector& gxxbest) 572 { 579 { 573 #ifdef G4TWISTDEBUG 580 #ifdef G4TWISTDEBUG 574 G4cout << "~~~~~ G4VTwistSurface::DistanceT 581 G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - Start ~~~~~~~~~" << G4endl; 575 G4cout << " Name : " << fName << G4end 582 G4cout << " Name : " << fName << G4endl; 576 G4cout << " gp : " << gp << G4endl; 583 G4cout << " gp : " << gp << G4endl; 577 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 584 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; 578 #endif 585 #endif 579 586 580 587 581 G4ThreeVector gxx[G4VSURFACENXX]; 588 G4ThreeVector gxx[G4VSURFACENXX]; 582 G4double distance[G4VSURFACENXX] ; 589 G4double distance[G4VSURFACENXX] ; 583 G4int areacode[G4VSURFACENXX] ; 590 G4int areacode[G4VSURFACENXX] ; 584 591 585 for (G4int i = 0 ; i<G4VSURFACENXX ; ++i ) 592 for (G4int i = 0 ; i<G4VSURFACENXX ; ++i ) 586 { 593 { 587 distance[i] = kInfinity ; 594 distance[i] = kInfinity ; 588 areacode[i] = sOutside ; 595 areacode[i] = sOutside ; 589 } 596 } 590 597 591 DistanceToSurface(gp, gxx, distance, areaco 598 DistanceToSurface(gp, gxx, distance, areacode); 592 gxxbest = gxx[0]; 599 gxxbest = gxx[0]; 593 600 594 #ifdef G4TWISTDEBUG 601 #ifdef G4TWISTDEBUG 595 G4cout << "~~~~~ G4VTwistSurface::DistanceT 602 G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - return ~~~~~~~~" << G4endl; 596 G4cout << " Name : " << fName << G 603 G4cout << " Name : " << fName << G4endl; 597 G4cout << " gxx : " << gxxbest << 604 G4cout << " gxx : " << gxxbest << G4endl; 598 G4cout << " bestdist : " << distance[0 605 G4cout << " bestdist : " << distance[0] << G4endl; 599 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 606 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl; 600 #endif 607 #endif 601 608 602 return distance[0]; 609 return distance[0]; 603 } 610 } 604 611 605 //============================================ 612 //===================================================================== 606 //* IsSameBoundary --------------------------- 613 //* IsSameBoundary ---------------------------------------------------- 607 614 608 G4bool 615 G4bool 609 G4VTwistSurface::IsSameBoundary(G4VTwistSurfac 616 G4VTwistSurface::IsSameBoundary(G4VTwistSurface* surf1, G4int areacode1, 610 G4VTwistSurfac 617 G4VTwistSurface* surf2, G4int areacode2 ) const 611 { 618 { 612 // 619 // 613 // IsSameBoundary 620 // IsSameBoundary 614 // 621 // 615 // checking tool whether two boundaries on 622 // checking tool whether two boundaries on different surfaces are same or not. 616 // 623 // 617 624 618 G4bool testbitmode = true; 625 G4bool testbitmode = true; 619 G4bool iscorner[2] = {IsCorner(areacode1, t 626 G4bool iscorner[2] = {IsCorner(areacode1, testbitmode), 620 IsCorner(areacode2, t 627 IsCorner(areacode2, testbitmode)}; 621 628 622 if (iscorner[0] && iscorner[1]) 629 if (iscorner[0] && iscorner[1]) 623 { 630 { 624 // on corner 631 // on corner 625 G4ThreeVector corner1 = 632 G4ThreeVector corner1 = 626 surf1->ComputeGlobalPoint(surf1->Ge 633 surf1->ComputeGlobalPoint(surf1->GetCorner(areacode1)); 627 G4ThreeVector corner2 = 634 G4ThreeVector corner2 = 628 surf2->ComputeGlobalPoint(surf2->Ge 635 surf2->ComputeGlobalPoint(surf2->GetCorner(areacode2)); 629 636 630 return (corner1 - corner2).mag() < kCarT << 637 if ((corner1 - corner2).mag() < kCarTolerance) { return true; } >> 638 else { return false; } 631 } 639 } 632 else if ((IsBoundary(areacode1, testbitmode 640 else if ((IsBoundary(areacode1, testbitmode) && (!iscorner[0])) && 633 (IsBoundary(areacode2, testbitmode 641 (IsBoundary(areacode2, testbitmode) && (!iscorner[1]))) 634 { 642 { 635 // on boundary 643 // on boundary 636 G4ThreeVector d1, d2, ld1, ld2; 644 G4ThreeVector d1, d2, ld1, ld2; 637 G4ThreeVector x01, x02, lx01, lx02; 645 G4ThreeVector x01, x02, lx01, lx02; 638 G4int type1, type2; 646 G4int type1, type2; 639 surf1->GetBoundaryParameters(areacode1, 647 surf1->GetBoundaryParameters(areacode1, ld1, lx01, type1); 640 surf2->GetBoundaryParameters(areacode2, 648 surf2->GetBoundaryParameters(areacode2, ld2, lx02, type2); 641 649 642 x01 = surf1->ComputeGlobalPoint(lx01); 650 x01 = surf1->ComputeGlobalPoint(lx01); 643 x02 = surf2->ComputeGlobalPoint(lx02); 651 x02 = surf2->ComputeGlobalPoint(lx02); 644 d1 = surf1->ComputeGlobalDirection(ld1) 652 d1 = surf1->ComputeGlobalDirection(ld1); 645 d2 = surf2->ComputeGlobalDirection(ld2) 653 d2 = surf2->ComputeGlobalDirection(ld2); 646 654 647 return (x01 - x02).mag() < kCarTolerance << 655 if ((x01 - x02).mag() < kCarTolerance && 648 && (d1 - d2).mag() < kCarTolerance << 656 (d1 - d2).mag() < kCarTolerance) { return true; } >> 657 else { return false; } 649 } 658 } 650 else 659 else 651 { 660 { 652 return false; 661 return false; 653 } 662 } 654 } 663 } 655 664 656 //============================================ 665 //===================================================================== 657 //* GetBoundaryParameters -------------------- 666 //* GetBoundaryParameters --------------------------------------------- 658 667 659 void G4VTwistSurface::GetBoundaryParameters(co 668 void G4VTwistSurface::GetBoundaryParameters(const G4int& areacode, 660 G 669 G4ThreeVector& d, 661 G 670 G4ThreeVector& x0, 662 G 671 G4int& boundarytype) const 663 { 672 { 664 // areacode must be one of them: 673 // areacode must be one of them: 665 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 674 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 666 // sAxis1 & sAxisMin, sAxis1 & sAxisMax. 675 // sAxis1 & sAxisMin, sAxis1 & sAxisMax. 667 676 668 for (const auto & boundary : fBoundaries) << 677 for (G4int i=0; i<4; ++i) 669 { 678 { 670 if (boundary.GetBoundaryParameters(areac << 679 if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0, >> 680 boundarytype)) 671 { 681 { 672 return; 682 return; 673 } 683 } 674 } 684 } 675 685 676 std::ostringstream message; 686 std::ostringstream message; 677 message << "Not registered boundary." << G4 687 message << "Not registered boundary." << G4endl 678 << " Boundary at areacode " 688 << " Boundary at areacode " << std::hex << areacode 679 << std::dec << G4endl 689 << std::dec << G4endl 680 << " is not registered."; 690 << " is not registered."; 681 G4Exception("G4VTwistSurface::GetBoundaryPa 691 G4Exception("G4VTwistSurface::GetBoundaryParameters()", "GeomSolids0002", 682 FatalException, message); 692 FatalException, message); 683 } 693 } 684 694 685 //============================================ 695 //===================================================================== 686 //* GetBoundaryAtPZ -------------------------- 696 //* GetBoundaryAtPZ --------------------------------------------------- 687 697 688 G4ThreeVector G4VTwistSurface::GetBoundaryAtPZ 698 G4ThreeVector G4VTwistSurface::GetBoundaryAtPZ(G4int areacode, 689 699 const G4ThreeVector& p) const 690 { 700 { 691 // areacode must be one of them: 701 // areacode must be one of them: 692 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 702 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 693 // sAxis1 & sAxisMin, sAxis1 & sAxisMax. 703 // sAxis1 & sAxisMin, sAxis1 & sAxisMax. 694 704 695 if (((areacode & sAxis0) != 0) && ((areacod << 705 if (areacode & sAxis0 && areacode & sAxis1) 696 { 706 { 697 std::ostringstream message; 707 std::ostringstream message; 698 message << "Point is in the corner area." 708 message << "Point is in the corner area." << G4endl 699 << " This function returns 709 << " This function returns " 700 << "a direction vector of a bound 710 << "a direction vector of a boundary line." << G4endl 701 << " areacode = " << areac 711 << " areacode = " << areacode; 702 G4Exception("G4VTwistSurface::GetBoundary 712 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "GeomSolids0003", 703 FatalException, message); 713 FatalException, message); 704 } 714 } 705 715 706 G4ThreeVector d; 716 G4ThreeVector d; 707 G4ThreeVector x0; 717 G4ThreeVector x0; 708 G4int boundarytype = 0; << 718 G4int boundarytype; 709 G4bool found = false; 719 G4bool found = false; 710 720 711 for (const auto & boundary : fBoundaries) << 721 for (G4int i=0; i<4; ++i) 712 { 722 { 713 if (boundary.GetBoundaryParameters(areac << 723 if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0, >> 724 boundarytype)) 714 { 725 { 715 found = true; 726 found = true; 716 continue; 727 continue; 717 } 728 } 718 } 729 } 719 730 720 if (!found) 731 if (!found) 721 { 732 { 722 std::ostringstream message; 733 std::ostringstream message; 723 message << "Not registered boundary." << 734 message << "Not registered boundary." << G4endl 724 << " Boundary at areacode 735 << " Boundary at areacode " << areacode << G4endl 725 << " is not registered."; 736 << " is not registered."; 726 G4Exception("G4VTwistSurface::GetBoundary 737 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "GeomSolids0002", 727 FatalException, message); 738 FatalException, message); 728 } 739 } 729 740 730 if (((boundarytype & sAxisPhi) == sAxisPhi) 741 if (((boundarytype & sAxisPhi) == sAxisPhi) || 731 ((boundarytype & sAxisRho) == sAxisRho) 742 ((boundarytype & sAxisRho) == sAxisRho)) 732 { 743 { 733 std::ostringstream message; 744 std::ostringstream message; 734 message << "Not a z-depended line boundar 745 message << "Not a z-depended line boundary." << G4endl 735 << " Boundary at areacode 746 << " Boundary at areacode " << areacode << G4endl 736 << " is not a z-depended l 747 << " is not a z-depended line."; 737 G4Exception("G4VTwistSurface::GetBoundary 748 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "GeomSolids0002", 738 FatalException, message); 749 FatalException, message); 739 } 750 } 740 return ((p.z() - x0.z()) / d.z()) * d + x0; 751 return ((p.z() - x0.z()) / d.z()) * d + x0; 741 } 752 } 742 753 743 //============================================ 754 //===================================================================== 744 //* SetCorner -------------------------------- 755 //* SetCorner --------------------------------------------------------- 745 756 746 void G4VTwistSurface::SetCorner(G4int areacode 757 void G4VTwistSurface::SetCorner(G4int areacode, 747 G4double x, G4 758 G4double x, G4double y, G4double z) 748 { 759 { 749 if ((areacode & sCorner) != sCorner) 760 if ((areacode & sCorner) != sCorner) 750 { 761 { 751 std::ostringstream message; 762 std::ostringstream message; 752 message << "Area code must represents cor 763 message << "Area code must represents corner." << G4endl 753 << " areacode " << areacod 764 << " areacode " << areacode; 754 G4Exception("G4VTwistSurface::SetCorner() 765 G4Exception("G4VTwistSurface::SetCorner()", "GeomSolids0002", 755 FatalException, message); 766 FatalException, message); 756 } 767 } 757 768 758 if ((areacode & sC0Min1Min) == sC0Min1Min) 769 if ((areacode & sC0Min1Min) == sC0Min1Min) { 759 fCorners[0].set(x, y, z); 770 fCorners[0].set(x, y, z); 760 } else if ((areacode & sC0Max1Min) == sC0Ma 771 } else if ((areacode & sC0Max1Min) == sC0Max1Min) { 761 fCorners[1].set(x, y, z); 772 fCorners[1].set(x, y, z); 762 } else if ((areacode & sC0Max1Max) == sC0Ma 773 } else if ((areacode & sC0Max1Max) == sC0Max1Max) { 763 fCorners[2].set(x, y, z); 774 fCorners[2].set(x, y, z); 764 } else if ((areacode & sC0Min1Max) == sC0Mi 775 } else if ((areacode & sC0Min1Max) == sC0Min1Max) { 765 fCorners[3].set(x, y, z); 776 fCorners[3].set(x, y, z); 766 } 777 } 767 } 778 } 768 779 769 //============================================ 780 //===================================================================== 770 //* SetBoundaryAxis -------------------------- 781 //* SetBoundaryAxis --------------------------------------------------- 771 782 772 void G4VTwistSurface::GetBoundaryAxis(G4int ar 783 void G4VTwistSurface::GetBoundaryAxis(G4int areacode, EAxis axis[]) const 773 { 784 { 774 if ((areacode & sBoundary) != sBoundary) { 785 if ((areacode & sBoundary) != sBoundary) { 775 G4Exception("G4VTwistSurface::GetBoundary 786 G4Exception("G4VTwistSurface::GetBoundaryAxis()", "GeomSolids0003", 776 FatalException, "Not located 787 FatalException, "Not located on a boundary!"); 777 } 788 } 778 for (G4int i=0; i<2; ++i) 789 for (G4int i=0; i<2; ++i) 779 { 790 { 780 G4int whichaxis = 0 ; 791 G4int whichaxis = 0 ; 781 if (i == 0) { 792 if (i == 0) { 782 whichaxis = sAxis0; 793 whichaxis = sAxis0; 783 } else if (i == 1) { 794 } else if (i == 1) { 784 whichaxis = sAxis1; 795 whichaxis = sAxis1; 785 } 796 } 786 797 787 // extracted axiscode of whichaxis 798 // extracted axiscode of whichaxis 788 G4int axiscode = whichaxis & sAxisMask & 799 G4int axiscode = whichaxis & sAxisMask & areacode ; 789 if (axiscode != 0) { << 800 if (axiscode) { 790 if (axiscode == (whichaxis & sAxisX)) 801 if (axiscode == (whichaxis & sAxisX)) { 791 axis[i] = kXAxis; 802 axis[i] = kXAxis; 792 } else if (axiscode == (whichaxis & s 803 } else if (axiscode == (whichaxis & sAxisY)) { 793 axis[i] = kYAxis; 804 axis[i] = kYAxis; 794 } else if (axiscode == (whichaxis & s 805 } else if (axiscode == (whichaxis & sAxisZ)) { 795 axis[i] = kZAxis; 806 axis[i] = kZAxis; 796 } else if (axiscode == (whichaxis & s 807 } else if (axiscode == (whichaxis & sAxisRho)) { 797 axis[i] = kRho; 808 axis[i] = kRho; 798 } else if (axiscode == (whichaxis & s 809 } else if (axiscode == (whichaxis & sAxisPhi)) { 799 axis[i] = kPhi; 810 axis[i] = kPhi; 800 } else { 811 } else { 801 std::ostringstream message; 812 std::ostringstream message; 802 message << "Not supported areacode. 813 message << "Not supported areacode." << G4endl 803 << " areacode " << a 814 << " areacode " << areacode; 804 G4Exception("G4VTwistSurface::GetBo 815 G4Exception("G4VTwistSurface::GetBoundaryAxis()", "GeomSolids0001", 805 FatalException, message 816 FatalException, message); 806 } 817 } 807 } 818 } 808 } 819 } 809 } 820 } 810 821 811 //============================================ 822 //===================================================================== 812 //* SetBoundaryLimit ------------------------- 823 //* SetBoundaryLimit -------------------------------------------------- 813 824 814 void G4VTwistSurface::GetBoundaryLimit(G4int a 825 void G4VTwistSurface::GetBoundaryLimit(G4int areacode, G4double limit[]) const 815 { 826 { 816 if ((areacode & sCorner) != 0) { << 827 if (areacode & sCorner) { 817 if ((areacode & sC0Min1Min) != 0) { << 828 if (areacode & sC0Min1Min) { 818 limit[0] = fAxisMin[0]; 829 limit[0] = fAxisMin[0]; 819 limit[1] = fAxisMin[1]; 830 limit[1] = fAxisMin[1]; 820 } else if ((areacode & sC0Max1Min) != 0) << 831 } else if (areacode & sC0Max1Min) { 821 limit[0] = fAxisMax[0]; 832 limit[0] = fAxisMax[0]; 822 limit[1] = fAxisMin[1]; 833 limit[1] = fAxisMin[1]; 823 } else if ((areacode & sC0Max1Max) != 0) << 834 } else if (areacode & sC0Max1Max) { 824 limit[0] = fAxisMax[0]; 835 limit[0] = fAxisMax[0]; 825 limit[1] = fAxisMax[1]; 836 limit[1] = fAxisMax[1]; 826 } else if ((areacode & sC0Min1Max) != 0) << 837 } else if (areacode & sC0Min1Max) { 827 limit[0] = fAxisMin[0]; 838 limit[0] = fAxisMin[0]; 828 limit[1] = fAxisMax[1]; 839 limit[1] = fAxisMax[1]; 829 } 840 } 830 } else if ((areacode & sBoundary) != 0) { << 841 } else if (areacode & sBoundary) { 831 if ((areacode & (sAxis0 | sAxisMin)) != << 842 if (areacode & (sAxis0 | sAxisMin)) { 832 limit[0] = fAxisMin[0]; 843 limit[0] = fAxisMin[0]; 833 } else if ((areacode & (sAxis1 | sAxisMi << 844 } else if (areacode & (sAxis1 | sAxisMin)) { 834 limit[0] = fAxisMin[1]; 845 limit[0] = fAxisMin[1]; 835 } else if ((areacode & (sAxis0 | sAxisMa << 846 } else if (areacode & (sAxis0 | sAxisMax)) { 836 limit[0] = fAxisMax[0]; 847 limit[0] = fAxisMax[0]; 837 } else if ((areacode & (sAxis1 | sAxisMa << 848 } else if (areacode & (sAxis1 | sAxisMax)) { 838 limit[0] = fAxisMax[1]; 849 limit[0] = fAxisMax[1]; 839 } 850 } 840 } else { 851 } else { 841 std::ostringstream message; 852 std::ostringstream message; 842 message << "Not located on a boundary!" < 853 message << "Not located on a boundary!" << G4endl 843 << " areacode " << areac 854 << " areacode " << areacode; 844 G4Exception("G4VTwistSurface::GetBoundary 855 G4Exception("G4VTwistSurface::GetBoundaryLimit()", "GeomSolids1002", 845 JustWarning, message); 856 JustWarning, message); 846 } 857 } 847 } 858 } 848 859 849 //============================================ 860 //===================================================================== 850 //* SetBoundary ------------------------------ 861 //* SetBoundary ------------------------------------------------------- 851 862 852 void G4VTwistSurface::SetBoundary(const G4int& 863 void G4VTwistSurface::SetBoundary(const G4int& axiscode, 853 const G4Thre 864 const G4ThreeVector& direction, 854 const G4Thre 865 const G4ThreeVector& x0, 855 const G4int& 866 const G4int& boundarytype) 856 { 867 { 857 G4int code = (~sAxisMask) & axiscode; 868 G4int code = (~sAxisMask) & axiscode; 858 if ((code == (sAxis0 & sAxisMin)) || 869 if ((code == (sAxis0 & sAxisMin)) || 859 (code == (sAxis0 & sAxisMax)) || 870 (code == (sAxis0 & sAxisMax)) || 860 (code == (sAxis1 & sAxisMin)) || 871 (code == (sAxis1 & sAxisMin)) || 861 (code == (sAxis1 & sAxisMax))) 872 (code == (sAxis1 & sAxisMax))) 862 { 873 { 863 G4bool done = false; 874 G4bool done = false; 864 for (auto & boundary : fBoundaries) << 875 for (auto i=0; i<4; ++i) 865 { 876 { 866 if (boundary.IsEmpty()) << 877 if (fBoundaries[i].IsEmpty()) 867 { 878 { 868 boundary.SetFields(axiscode, direc << 879 fBoundaries[i].SetFields(axiscode, direction, >> 880 x0, boundarytype); 869 done = true; 881 done = true; 870 break; 882 break; 871 } 883 } 872 } 884 } 873 885 874 if (!done) 886 if (!done) 875 { 887 { 876 G4Exception("G4VTwistSurface::SetBoun 888 G4Exception("G4VTwistSurface::SetBoundary()", "GeomSolids0003", 877 FatalException, "Number 889 FatalException, "Number of boundary exceeding 4!"); 878 } 890 } 879 } 891 } 880 else 892 else 881 { 893 { 882 std::ostringstream message; 894 std::ostringstream message; 883 message << "Invalid axis-code." << G4end 895 message << "Invalid axis-code." << G4endl 884 << " axiscode = " 896 << " axiscode = " 885 << std::hex << axiscode << std:: 897 << std::hex << axiscode << std::dec; 886 G4Exception("G4VTwistSurface::SetBoundar 898 G4Exception("G4VTwistSurface::SetBoundary()", "GeomSolids0003", 887 FatalException, message); 899 FatalException, message); 888 } 900 } 889 } 901 } 890 902 891 //============================================ 903 //===================================================================== 892 //* GetFace ---------------------------------- 904 //* GetFace ----------------------------------------------------------- 893 905 894 G4int G4VTwistSurface::GetFace( G4int i, G4int 906 G4int G4VTwistSurface::GetFace( G4int i, G4int j, G4int k, 895 G4int n, G4int 907 G4int n, G4int iside ) 896 { 908 { 897 // this is the face mapping function 909 // this is the face mapping function 898 // (i,j) -> face number 910 // (i,j) -> face number 899 911 900 if ( iside == 0 ) { 912 if ( iside == 0 ) { 901 return i * ( k - 1 ) + j ; 913 return i * ( k - 1 ) + j ; 902 } 914 } 903 915 904 else if ( iside == 1 ) { 916 else if ( iside == 1 ) { 905 return (k-1)*(k-1) + i*(k-1) + j ; 917 return (k-1)*(k-1) + i*(k-1) + j ; 906 } 918 } 907 919 908 else if ( iside == 2 ) { 920 else if ( iside == 2 ) { 909 return 2*(k-1)*(k-1) + i*(k-1) + j ; 921 return 2*(k-1)*(k-1) + i*(k-1) + j ; 910 } 922 } 911 923 912 else if ( iside == 3 ) { 924 else if ( iside == 3 ) { 913 return 2*(k-1)*(k-1) + (n-1)*(k-1) + i*(k- 925 return 2*(k-1)*(k-1) + (n-1)*(k-1) + i*(k-1) + j ; 914 } 926 } 915 927 916 else if ( iside == 4 ) { 928 else if ( iside == 4 ) { 917 return 2*(k-1)*(k-1) + 2*(n-1)*(k-1) + i*( 929 return 2*(k-1)*(k-1) + 2*(n-1)*(k-1) + i*(k-1) + j ; 918 } 930 } 919 931 920 else if ( iside == 5 ) { 932 else if ( iside == 5 ) { 921 return 2*(k-1)*(k-1) + 3*(n-1)*(k-1) + i*( 933 return 2*(k-1)*(k-1) + 3*(n-1)*(k-1) + i*(k-1) + j ; 922 } 934 } 923 935 924 else { 936 else { 925 std::ostringstream message; 937 std::ostringstream message; 926 message << "Not correct side number: " 938 message << "Not correct side number: " 927 << GetName() << G4endl 939 << GetName() << G4endl 928 << "iside is " << iside << " but s 940 << "iside is " << iside << " but should be " 929 << "0,1,2,3,4 or 5" << "."; 941 << "0,1,2,3,4 or 5" << "."; 930 G4Exception("G4TwistSurface::G4GetFace()", 942 G4Exception("G4TwistSurface::G4GetFace()", "GeomSolids0002", 931 FatalException, message); 943 FatalException, message); 932 } 944 } 933 945 934 return -1 ; // wrong face 946 return -1 ; // wrong face 935 } 947 } 936 948 937 //============================================ 949 //===================================================================== 938 //* GetNode ---------------------------------- 950 //* GetNode ----------------------------------------------------------- 939 951 940 G4int G4VTwistSurface::GetNode( G4int i, G4int 952 G4int G4VTwistSurface::GetNode( G4int i, G4int j, G4int k, 941 G4int n, G4int 953 G4int n, G4int iside ) 942 { 954 { 943 // this is the node mapping function 955 // this is the node mapping function 944 // (i,j) -> node number 956 // (i,j) -> node number 945 // Depends on the side iside and the used me 957 // Depends on the side iside and the used meshing of the surface 946 958 947 if ( iside == 0 ) 959 if ( iside == 0 ) 948 { 960 { 949 // lower endcap is kxk squared. 961 // lower endcap is kxk squared. 950 // n = k 962 // n = k 951 return i * k + j ; 963 return i * k + j ; 952 } 964 } 953 965 954 if ( iside == 1 ) 966 if ( iside == 1 ) 955 { 967 { 956 // upper endcap is kxk squared. Shift by k 968 // upper endcap is kxk squared. Shift by k*k 957 // n = k 969 // n = k 958 return k*k + i*k + j ; 970 return k*k + i*k + j ; 959 } 971 } 960 972 961 else if ( iside == 2 ) 973 else if ( iside == 2 ) 962 { 974 { 963 // front side. 975 // front side. 964 if ( i == 0 ) { return j ; } 976 if ( i == 0 ) { return j ; } 965 else if ( i == n-1 ) { return k*k + j ; } 977 else if ( i == n-1 ) { return k*k + j ; } 966 else { return 2*k*k + 4*(i 978 else { return 2*k*k + 4*(i-1)*(k-1) + j ; } 967 } 979 } 968 980 969 else if ( iside == 3 ) 981 else if ( iside == 3 ) 970 { 982 { 971 // right side 983 // right side 972 if ( i == 0 ) { return (j+1)* 984 if ( i == 0 ) { return (j+1)*k - 1 ; } 973 else if ( i == n-1 ) { return k*k + (j+1)* 985 else if ( i == n-1 ) { return k*k + (j+1)*k - 1 ; } 974 else 986 else 975 { 987 { 976 return 2*k*k + 4*(i-1)*(k-1) + (k-1) + j 988 return 2*k*k + 4*(i-1)*(k-1) + (k-1) + j ; 977 } 989 } 978 } 990 } 979 else if ( iside == 4 ) 991 else if ( iside == 4 ) 980 { 992 { 981 // back side 993 // back side 982 if ( i == 0 ) { return k*k - 1 - 994 if ( i == 0 ) { return k*k - 1 - j ; } // reversed order 983 else if ( i == n-1 ) { return 2*k*k - 1 - 995 else if ( i == n-1 ) { return 2*k*k - 1 - j ; } // reversed order 984 else 996 else 985 { 997 { 986 return 2*k*k + 4*(i-1)*(k-1) + 2*(k-1) + 998 return 2*k*k + 4*(i-1)*(k-1) + 2*(k-1) + j ; // normal order 987 } 999 } 988 } 1000 } 989 else if ( iside == 5 ) 1001 else if ( iside == 5 ) 990 { 1002 { 991 // left side 1003 // left side 992 if ( i == 0 ) { return k*k - (j+1 1004 if ( i == 0 ) { return k*k - (j+1)*k ; } // reversed order 993 else if ( i == n-1) { return 2*k*k - (j+1 1005 else if ( i == n-1) { return 2*k*k - (j+1)*k ; } // reverded order 994 else 1006 else 995 { 1007 { 996 if ( j == k-1 ) { return 2*k*k + 4*(i-1) 1008 if ( j == k-1 ) { return 2*k*k + 4*(i-1)*(k-1) ; } // special case 997 else 1009 else 998 { 1010 { 999 return 2*k*k + 4*(i-1)*(k-1) + 3*(k-1) 1011 return 2*k*k + 4*(i-1)*(k-1) + 3*(k-1) + j ; // normal order 1000 } 1012 } 1001 } 1013 } 1002 } 1014 } 1003 else 1015 else 1004 { 1016 { 1005 std::ostringstream message; 1017 std::ostringstream message; 1006 message << "Not correct side number: " 1018 message << "Not correct side number: " 1007 << GetName() << G4endl 1019 << GetName() << G4endl 1008 << "iside is " << iside << " but 1020 << "iside is " << iside << " but should be " 1009 << "0,1,2,3,4 or 5" << "."; 1021 << "0,1,2,3,4 or 5" << "."; 1010 G4Exception("G4TwistSurface::G4GetNode()" 1022 G4Exception("G4TwistSurface::G4GetNode()", "GeomSolids0002", 1011 FatalException, message); 1023 FatalException, message); 1012 } 1024 } 1013 return -1 ; // wrong node 1025 return -1 ; // wrong node 1014 } 1026 } 1015 1027 1016 //=========================================== 1028 //===================================================================== 1017 //* GetEdgeVisiblility ---------------------- 1029 //* GetEdgeVisiblility ------------------------------------------------ 1018 1030 1019 G4int G4VTwistSurface::GetEdgeVisibility( G4i 1031 G4int G4VTwistSurface::GetEdgeVisibility( G4int i, G4int j, G4int k, G4int n, 1020 G4i 1032 G4int number, G4int orientation ) 1021 { 1033 { 1022 // clockwise filling -> positive or 1034 // clockwise filling -> positive orientation 1023 // counter clockwise filling -> negative or 1035 // counter clockwise filling -> negative orientation 1024 1036 1025 // 1037 // 1026 // d C c 1038 // d C c 1027 // +------+ 1039 // +------+ 1028 // | | 1040 // | | 1029 // | | 1041 // | | 1030 // | | 1042 // | | 1031 // D | |B 1043 // D | |B 1032 // | | 1044 // | | 1033 // | | 1045 // | | 1034 // | | 1046 // | | 1035 // +------+ 1047 // +------+ 1036 // a A b 1048 // a A b 1037 // 1049 // 1038 // a = +--+ A = ---+ 1050 // a = +--+ A = ---+ 1039 // b = --++ B = --+- 1051 // b = --++ B = --+- 1040 // c = -++- C = -+-- 1052 // c = -++- C = -+-- 1041 // d = ++-- D = +--- 1053 // d = ++-- D = +--- 1042 1054 1043 1055 1044 // check first invisible faces 1056 // check first invisible faces 1045 1057 1046 if ( ( i>0 && i<n-2 ) && ( j>0 && j<k-2 ) ) 1058 if ( ( i>0 && i<n-2 ) && ( j>0 && j<k-2 ) ) 1047 { 1059 { 1048 return -1 ; // always invisible, signs: 1060 return -1 ; // always invisible, signs: ---- 1049 } 1061 } 1050 1062 1051 // change first the vertex number (depends 1063 // change first the vertex number (depends on the orientation) 1052 // 0,1,2,3 -> 3,2,1,0 1064 // 0,1,2,3 -> 3,2,1,0 1053 if ( orientation < 0 ) { number = ( 3 - num 1065 if ( orientation < 0 ) { number = ( 3 - number ) ; } 1054 1066 1055 // check true edges 1067 // check true edges 1056 if ( ( j>=1 && j<=k-3 ) ) 1068 if ( ( j>=1 && j<=k-3 ) ) 1057 { 1069 { 1058 if ( i == 0 ) { // signs (A): -- 1070 if ( i == 0 ) { // signs (A): ---+ 1059 return ( number == 3 ) ? 1 : -1 ; 1071 return ( number == 3 ) ? 1 : -1 ; 1060 } 1072 } 1061 1073 1062 else if ( i == n-2 ) { // signs (C): -+ 1074 else if ( i == n-2 ) { // signs (C): -+-- 1063 return ( number == 1 ) ? 1 : -1 ; 1075 return ( number == 1 ) ? 1 : -1 ; 1064 } 1076 } 1065 1077 1066 else 1078 else 1067 { 1079 { 1068 std::ostringstream message; 1080 std::ostringstream message; 1069 message << "Not correct face number: " 1081 message << "Not correct face number: " << GetName() << " !"; 1070 G4Exception("G4TwistSurface::G4GetEdgeV 1082 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", 1071 "GeomSolids0003", FatalExce 1083 "GeomSolids0003", FatalException, message); 1072 } 1084 } 1073 } 1085 } 1074 1086 1075 if ( ( i>=1 && i<=n-3 ) ) 1087 if ( ( i>=1 && i<=n-3 ) ) 1076 { 1088 { 1077 if ( j == 0 ) { // signs (D): +- 1089 if ( j == 0 ) { // signs (D): +--- 1078 return ( number == 0 ) ? 1 : -1 ; 1090 return ( number == 0 ) ? 1 : -1 ; 1079 } 1091 } 1080 1092 1081 else if ( j == k-2 ) { // signs (B): -- 1093 else if ( j == k-2 ) { // signs (B): --+- 1082 return ( number == 2 ) ? 1 : -1 ; 1094 return ( number == 2 ) ? 1 : -1 ; 1083 } 1095 } 1084 1096 1085 else 1097 else 1086 { 1098 { 1087 std::ostringstream message; 1099 std::ostringstream message; 1088 message << "Not correct face number: " 1100 message << "Not correct face number: " << GetName() << " !"; 1089 G4Exception("G4TwistSurface::G4GetEdgeV 1101 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", 1090 "GeomSolids0003", FatalExce 1102 "GeomSolids0003", FatalException, message); 1091 } 1103 } 1092 } 1104 } 1093 1105 1094 // now the corners 1106 // now the corners 1095 if ( i == 0 && j == 0 ) { // signs 1107 if ( i == 0 && j == 0 ) { // signs (a) : +--+ 1096 return ( number == 0 || number == 3 ) ? 1 1108 return ( number == 0 || number == 3 ) ? 1 : -1 ; 1097 } 1109 } 1098 else if ( i == 0 && j == k-2 ) { // signs 1110 else if ( i == 0 && j == k-2 ) { // signs (b) : --++ 1099 return ( number == 2 || number == 3 ) ? 1 1111 return ( number == 2 || number == 3 ) ? 1 : -1 ; 1100 } 1112 } 1101 else if ( i == n-2 && j == k-2 ) { // signs 1113 else if ( i == n-2 && j == k-2 ) { // signs (c) : -++- 1102 return ( number == 1 || number == 2 ) ? 1 1114 return ( number == 1 || number == 2 ) ? 1 : -1 ; 1103 } 1115 } 1104 else if ( i == n-2 && j == 0 ) { // signs 1116 else if ( i == n-2 && j == 0 ) { // signs (d) : ++-- 1105 return ( number == 0 || number == 1 ) ? 1 1117 return ( number == 0 || number == 1 ) ? 1 : -1 ; 1106 } 1118 } 1107 else 1119 else 1108 { 1120 { 1109 std::ostringstream message; 1121 std::ostringstream message; 1110 message << "Not correct face number: " << 1122 message << "Not correct face number: " << GetName() << " !"; 1111 G4Exception("G4TwistSurface::G4GetEdgeVis 1123 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", 1112 "GeomSolids0003", FatalExcept 1124 "GeomSolids0003", FatalException, message); 1113 } 1125 } 1114 1126 1115 std::ostringstream message; 1127 std::ostringstream message; 1116 message << "Not correct face number: " << G 1128 message << "Not correct face number: " << GetName() << " !"; 1117 G4Exception("G4TwistSurface::G4GetEdgeVisib 1129 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "GeomSolids0003", 1118 FatalException, message); 1130 FatalException, message); 1119 1131 1120 return 0 ; 1132 return 0 ; 1121 } 1133 } 1122 1134 1123 1135 1124 //=========================================== 1136 //===================================================================== 1125 //* DebugPrint ------------------------------ 1137 //* DebugPrint -------------------------------------------------------- 1126 1138 1127 void G4VTwistSurface::DebugPrint() const 1139 void G4VTwistSurface::DebugPrint() const 1128 { 1140 { 1129 G4ThreeVector A = fRot * GetCorner(sC0Min1 1141 G4ThreeVector A = fRot * GetCorner(sC0Min1Min) + fTrans; 1130 G4ThreeVector B = fRot * GetCorner(sC0Max1 1142 G4ThreeVector B = fRot * GetCorner(sC0Max1Min) + fTrans; 1131 G4ThreeVector C = fRot * GetCorner(sC0Max1 1143 G4ThreeVector C = fRot * GetCorner(sC0Max1Max) + fTrans; 1132 G4ThreeVector D = fRot * GetCorner(sC0Min1 1144 G4ThreeVector D = fRot * GetCorner(sC0Min1Max) + fTrans; 1133 1145 1134 G4cout << "/* G4VTwistSurface::DebugPrint( 1146 G4cout << "/* G4VTwistSurface::DebugPrint():--------------------------" 1135 << G4endl; 1147 << G4endl; 1136 G4cout << "/* Name = " << fName << G4endl; 1148 G4cout << "/* Name = " << fName << G4endl; 1137 G4cout << "/* Axis = " << std::hex << fAxi 1149 G4cout << "/* Axis = " << std::hex << fAxis[0] << " " 1138 << std::hex << fAxis[1] 1150 << std::hex << fAxis[1] 1139 << " (0,1,2,3,5 = kXAxis,kYAxis,kZA 1151 << " (0,1,2,3,5 = kXAxis,kYAxis,kZAxis,kRho,kPhi)" 1140 << std::dec << G4endl; 1152 << std::dec << G4endl; 1141 G4cout << "/* BoundaryLimit(in local) fAxi 1153 G4cout << "/* BoundaryLimit(in local) fAxis0(min, max) = ("<<fAxisMin[0] 1142 << ", " << fAxisMax[0] << ")" << G4 1154 << ", " << fAxisMax[0] << ")" << G4endl; 1143 G4cout << "/* BoundaryLimit(in local) fAxi 1155 G4cout << "/* BoundaryLimit(in local) fAxis1(min, max) = ("<<fAxisMin[1] 1144 << ", " << fAxisMax[1] << ")" << G4 1156 << ", " << fAxisMax[1] << ")" << G4endl; 1145 G4cout << "/* Cornar point sC0Min1Min = " 1157 G4cout << "/* Cornar point sC0Min1Min = " << A << G4endl; 1146 G4cout << "/* Cornar point sC0Max1Min = " 1158 G4cout << "/* Cornar point sC0Max1Min = " << B << G4endl; 1147 G4cout << "/* Cornar point sC0Max1Max = " 1159 G4cout << "/* Cornar point sC0Max1Max = " << C << G4endl; 1148 G4cout << "/* Cornar point sC0Min1Max = " 1160 G4cout << "/* Cornar point sC0Min1Max = " << D << G4endl; 1149 G4cout << "/*----------------------------- 1161 G4cout << "/*---------------------------------------------------------" 1150 << G4endl; 1162 << G4endl; 1151 } 1163 } 1152 1164 1153 //=========================================== 1165 //===================================================================== 1154 // G4VTwistSurface::CurrentStatus class 1166 // G4VTwistSurface::CurrentStatus class 1155 //=========================================== 1167 //===================================================================== 1156 1168 1157 //=========================================== 1169 //===================================================================== 1158 //* CurrentStatus::CurrentStatus ------------ 1170 //* CurrentStatus::CurrentStatus -------------------------------------- 1159 1171 1160 G4VTwistSurface::CurrentStatus::CurrentStatus 1172 G4VTwistSurface::CurrentStatus::CurrentStatus() 1161 { 1173 { 1162 for (size_t i=0; i<G4VSURFACENXX; ++i) 1174 for (size_t i=0; i<G4VSURFACENXX; ++i) 1163 { 1175 { 1164 fDistance[i] = kInfinity; 1176 fDistance[i] = kInfinity; 1165 fAreacode[i] = sOutside; 1177 fAreacode[i] = sOutside; 1166 fIsValid[i] = false; 1178 fIsValid[i] = false; 1167 fXX[i].set(kInfinity, kInfinity, kInfinit 1179 fXX[i].set(kInfinity, kInfinity, kInfinity); 1168 } 1180 } 1169 fNXX = 0; 1181 fNXX = 0; 1170 fLastp.set(kInfinity, kInfinity, kInfinity) 1182 fLastp.set(kInfinity, kInfinity, kInfinity); 1171 fLastv.set(kInfinity, kInfinity, kInfinity) 1183 fLastv.set(kInfinity, kInfinity, kInfinity); 1172 fLastValidate = kUninitialized; 1184 fLastValidate = kUninitialized; 1173 fDone = false; 1185 fDone = false; 1174 } 1186 } 1175 1187 1176 //=========================================== 1188 //===================================================================== 1177 //* CurrentStatus::~CurrentStatus ----------- 1189 //* CurrentStatus::~CurrentStatus ------------------------------------- 1178 1190 1179 G4VTwistSurface::CurrentStatus::~CurrentStatu 1191 G4VTwistSurface::CurrentStatus::~CurrentStatus() 1180 = default; << 1192 { >> 1193 } 1181 1194 1182 //=========================================== 1195 //===================================================================== 1183 //* CurrentStatus::SetCurrentStatus --------- 1196 //* CurrentStatus::SetCurrentStatus ----------------------------------- 1184 1197 1185 void 1198 void 1186 G4VTwistSurface::CurrentStatus::SetCurrentSta 1199 G4VTwistSurface::CurrentStatus::SetCurrentStatus(G4int i, 1187 1200 G4ThreeVector& xx, 1188 1201 G4double& dist, 1189 1202 G4int& areacode, 1190 1203 G4bool& isvalid, 1191 1204 G4int nxx, 1192 1205 EValidate validate, 1193 co 1206 const G4ThreeVector* p, 1194 co 1207 const G4ThreeVector* v) 1195 { 1208 { 1196 fDistance[i] = dist; 1209 fDistance[i] = dist; 1197 fAreacode[i] = areacode; 1210 fAreacode[i] = areacode; 1198 fIsValid[i] = isvalid; 1211 fIsValid[i] = isvalid; 1199 fXX[i] = xx; 1212 fXX[i] = xx; 1200 fNXX = nxx; 1213 fNXX = nxx; 1201 fLastValidate = validate; 1214 fLastValidate = validate; 1202 if (p != nullptr) 1215 if (p != nullptr) 1203 { 1216 { 1204 fLastp = *p; 1217 fLastp = *p; 1205 } 1218 } 1206 else 1219 else 1207 { 1220 { 1208 G4Exception("G4VTwistSurface::CurrentStat 1221 G4Exception("G4VTwistSurface::CurrentStatus::SetCurrentStatus()", 1209 "GeomSolids0003", FatalExcept 1222 "GeomSolids0003", FatalException, "SetCurrentStatus: p = 0!"); 1210 } 1223 } 1211 if (v != nullptr) 1224 if (v != nullptr) 1212 { 1225 { 1213 fLastv = *v; 1226 fLastv = *v; 1214 } 1227 } 1215 else 1228 else 1216 { 1229 { 1217 fLastv.set(kInfinity, kInfinity, kInfinit 1230 fLastv.set(kInfinity, kInfinity, kInfinity); 1218 } 1231 } 1219 fDone = true; 1232 fDone = true; 1220 } 1233 } 1221 1234 1222 //=========================================== 1235 //===================================================================== 1223 //* CurrentStatus::ResetfDone --------------- 1236 //* CurrentStatus::ResetfDone ----------------------------------------- 1224 1237 1225 void 1238 void 1226 G4VTwistSurface::CurrentStatus::ResetfDone(EV 1239 G4VTwistSurface::CurrentStatus::ResetfDone(EValidate validate, 1227 const G4 1240 const G4ThreeVector* p, 1228 const G4 1241 const G4ThreeVector* v) 1229 1242 1230 { 1243 { 1231 if (validate == fLastValidate && p != nullp 1244 if (validate == fLastValidate && p != nullptr && *p == fLastp) 1232 { 1245 { 1233 if (v == nullptr || (*v == fLastv)) retu 1246 if (v == nullptr || (*v == fLastv)) return; 1234 } 1247 } 1235 G4ThreeVector xx(kInfinity, kInfinity, kInf 1248 G4ThreeVector xx(kInfinity, kInfinity, kInfinity); 1236 for (size_t i=0; i<G4VSURFACENXX; ++i) 1249 for (size_t i=0; i<G4VSURFACENXX; ++i) 1237 { 1250 { 1238 fDistance[i] = kInfinity; 1251 fDistance[i] = kInfinity; 1239 fAreacode[i] = sOutside; 1252 fAreacode[i] = sOutside; 1240 fIsValid[i] = false; 1253 fIsValid[i] = false; 1241 fXX[i] = xx; // bug in old code ( was f 1254 fXX[i] = xx; // bug in old code ( was fXX[i] = xx[i] ) 1242 } 1255 } 1243 fNXX = 0; 1256 fNXX = 0; 1244 fLastp.set(kInfinity, kInfinity, kInfinity) 1257 fLastp.set(kInfinity, kInfinity, kInfinity); 1245 fLastv.set(kInfinity, kInfinity, kInfinity) 1258 fLastv.set(kInfinity, kInfinity, kInfinity); 1246 fLastValidate = kUninitialized; 1259 fLastValidate = kUninitialized; 1247 fDone = false; 1260 fDone = false; 1248 } 1261 } 1249 1262 1250 //=========================================== 1263 //===================================================================== 1251 //* CurrentStatus::DebugPrint --------------- 1264 //* CurrentStatus::DebugPrint ----------------------------------------- 1252 1265 1253 void 1266 void 1254 G4VTwistSurface::CurrentStatus::DebugPrint() 1267 G4VTwistSurface::CurrentStatus::DebugPrint() const 1255 { 1268 { 1256 G4cout << "CurrentStatus::Dist0,1= " << fDi 1269 G4cout << "CurrentStatus::Dist0,1= " << fDistance[0] 1257 << " " << fDistance[1] << " areacode 1270 << " " << fDistance[1] << " areacode = " << fAreacode[0] 1258 << " " << fAreacode[1] << G4endl; 1271 << " " << fAreacode[1] << G4endl; 1259 } 1272 } 1260 1273 1261 //=========================================== 1274 //===================================================================== 1262 // G4VTwistSurface::Boundary class 1275 // G4VTwistSurface::Boundary class 1263 //=========================================== 1276 //===================================================================== 1264 1277 1265 //=========================================== 1278 //===================================================================== >> 1279 //* Boundary::Boundary ------------------------------------------------ >> 1280 >> 1281 G4VTwistSurface::Boundary::Boundary() >> 1282 : fBoundaryAcode(-1), fBoundaryType(0) >> 1283 { >> 1284 } >> 1285 >> 1286 //===================================================================== >> 1287 //* Boundary::~Boundary ----------------------------------------------- >> 1288 >> 1289 G4VTwistSurface::Boundary::~Boundary() >> 1290 { >> 1291 } >> 1292 >> 1293 //===================================================================== 1266 //* Boundary::SetFields --------------------- 1294 //* Boundary::SetFields ----------------------------------------------- 1267 1295 1268 void 1296 void 1269 G4VTwistSurface::Boundary::SetFields(const G4 1297 G4VTwistSurface::Boundary::SetFields(const G4int& areacode, 1270 const G4Three 1298 const G4ThreeVector& d, 1271 const G4Three 1299 const G4ThreeVector& x0, 1272 const G4int& 1300 const G4int& boundarytype) 1273 { 1301 { 1274 fBoundaryAcode = areacode; 1302 fBoundaryAcode = areacode; 1275 fBoundaryDirection = d; 1303 fBoundaryDirection = d; 1276 fBoundaryX0 = x0; 1304 fBoundaryX0 = x0; 1277 fBoundaryType = boundarytype; 1305 fBoundaryType = boundarytype; 1278 } 1306 } 1279 1307 1280 //=========================================== 1308 //===================================================================== 1281 //* Boundary::IsEmpty ----------------------- 1309 //* Boundary::IsEmpty ------------------------------------------------- 1282 1310 1283 G4bool G4VTwistSurface::Boundary::IsEmpty() c 1311 G4bool G4VTwistSurface::Boundary::IsEmpty() const 1284 { 1312 { 1285 return fBoundaryAcode == -1; << 1313 if (fBoundaryAcode == -1) return true; >> 1314 return false; 1286 } 1315 } 1287 1316 1288 //=========================================== 1317 //===================================================================== 1289 //* Boundary::GetBoundaryParameters --------- 1318 //* Boundary::GetBoundaryParameters ----------------------------------- 1290 1319 1291 G4bool 1320 G4bool 1292 G4VTwistSurface::Boundary::GetBoundaryParamet 1321 G4VTwistSurface::Boundary::GetBoundaryParameters(const G4int& areacode, 1293 1322 G4ThreeVector& d, 1294 1323 G4ThreeVector& x0, 1295 1324 G4int& boundarytype) const 1296 { 1325 { 1297 // areacode must be one of them: 1326 // areacode must be one of them: 1298 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 1327 // sAxis0 & sAxisMin, sAxis0 & sAxisMax, 1299 // sAxis1 & sAxisMin, sAxis1 & sAxisMax 1328 // sAxis1 & sAxisMin, sAxis1 & sAxisMax 1300 // 1329 // 1301 if (((areacode & sAxis0) != 0) && ((areacod << 1330 if ((areacode & sAxis0) && (areacode & sAxis1)) 1302 { 1331 { 1303 std::ostringstream message; 1332 std::ostringstream message; 1304 message << "Located in the corner area." 1333 message << "Located in the corner area." << G4endl 1305 << " This function returns 1334 << " This function returns a direction vector of " 1306 << "a boundary line." << G4endl 1335 << "a boundary line." << G4endl 1307 << " areacode = " << areac 1336 << " areacode = " << areacode; 1308 G4Exception("G4VTwistSurface::Boundary::G 1337 G4Exception("G4VTwistSurface::Boundary::GetBoundaryParameters()", 1309 "GeomSolids0003", FatalExcept 1338 "GeomSolids0003", FatalException, message); 1310 } 1339 } 1311 if ((areacode & sSizeMask) != (fBoundaryAco 1340 if ((areacode & sSizeMask) != (fBoundaryAcode & sSizeMask)) 1312 { 1341 { 1313 return false; 1342 return false; 1314 } 1343 } 1315 d = fBoundaryDirection; 1344 d = fBoundaryDirection; 1316 x0 = fBoundaryX0; 1345 x0 = fBoundaryX0; 1317 boundarytype = fBoundaryType; 1346 boundarytype = fBoundaryType; 1318 return true; 1347 return true; 1319 } 1348 } 1320 1349