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