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 // G4TwistTrapParallelSide implementation << 27 // 26 // 28 // Author: Oliver Link (Oliver.Link@cern.ch) << 27 // $Id: G4TwistTrapParallelSide.cc,v >> 28 // GEANT4 tag $Name: not supported by cvs2svn $ >> 29 // >> 30 // >> 31 // -------------------------------------------------------------------- >> 32 // GEANT 4 class source file >> 33 // >> 34 // >> 35 // G4TwistTrapParallelSide.cc >> 36 // >> 37 // Author: >> 38 // >> 39 // Oliver Link (Oliver.Link@cern.ch) >> 40 // 29 // ------------------------------------------- 41 // -------------------------------------------------------------------- 30 42 31 #include <cmath> 43 #include <cmath> 32 44 33 #include "G4TwistTrapParallelSide.hh" 45 #include "G4TwistTrapParallelSide.hh" 34 #include "G4PhysicalConstants.hh" 46 #include "G4PhysicalConstants.hh" 35 #include "G4JTPolynomialSolver.hh" 47 #include "G4JTPolynomialSolver.hh" 36 48 37 //============================================ 49 //===================================================================== 38 //* constructors ----------------------------- 50 //* constructors ------------------------------------------------------ 39 51 40 G4TwistTrapParallelSide::G4TwistTrapParallelSi << 52 G4TwistTrapParallelSide::G4TwistTrapParallelSide(const G4String &name, 41 G4double PhiTwist, << 53 G4double PhiTwist, // twist angle 42 G4double pDz, << 54 G4double pDz, // half z lenght 43 G4double pTheta, << 55 G4double pTheta, // direction between end planes 44 G4double pPhi, << 56 G4double pPhi, // defined by polar and azimutal angles. 45 G4double pDy1, << 57 G4double pDy1, // half y length at -pDz 46 G4double pDx1, << 58 G4double pDx1, // half x length at -pDz,-pDy 47 G4double pDx2, << 59 G4double pDx2, // half x length at -pDz,+pDy 48 G4double pDy2, << 60 G4double pDy2, // half y length at +pDz 49 G4double pDx3, << 61 G4double pDx3, // half x length at +pDz,-pDy 50 G4double pDx4, << 62 G4double pDx4, // half x length at +pDz,+pDy 51 G4double pAlph, << 63 G4double pAlph, // tilt angle at +pDz 52 G4double AngleSide << 64 G4double AngleSide // parity 53 65 ) 54 : G4VTwistSurface(name) << 66 : G4VTwistSurface(name), fDy(0.) 55 { 67 { 56 68 57 fAxis[0] = kXAxis; // in local coordinate 69 fAxis[0] = kXAxis; // in local coordinate system 58 fAxis[1] = kZAxis; 70 fAxis[1] = kZAxis; 59 fAxisMin[0] = -kInfinity ; // X Axis bounda 71 fAxisMin[0] = -kInfinity ; // X Axis boundary 60 fAxisMax[0] = kInfinity ; // depends on 72 fAxisMax[0] = kInfinity ; // depends on z !! 61 fAxisMin[1] = -pDz ; // Z Axis bounda << 73 fAxisMin[1] = -pDz ; // Z Axis boundary 62 fAxisMax[1] = pDz ; 74 fAxisMax[1] = pDz ; 63 75 64 fDx1 = pDx1 ; 76 fDx1 = pDx1 ; 65 fDx2 = pDx2 ; 77 fDx2 = pDx2 ; 66 fDx3 = pDx3 ; 78 fDx3 = pDx3 ; 67 fDx4 = pDx4 ; 79 fDx4 = pDx4 ; 68 80 69 fDy1 = pDy1 ; 81 fDy1 = pDy1 ; 70 fDy2 = pDy2 ; 82 fDy2 = pDy2 ; 71 83 72 fDz = pDz ; 84 fDz = pDz ; 73 85 74 fAlph = pAlph ; 86 fAlph = pAlph ; 75 fTAlph = std::tan(fAlph) ; 87 fTAlph = std::tan(fAlph) ; 76 88 77 fTheta = pTheta ; 89 fTheta = pTheta ; 78 fPhi = pPhi ; 90 fPhi = pPhi ; 79 91 80 // precalculate frequently used parameters << 92 // precalculate frequently used parameters 81 // << 82 fDx4plus2 = fDx4 + fDx2 ; 93 fDx4plus2 = fDx4 + fDx2 ; 83 fDx4minus2 = fDx4 - fDx2 ; 94 fDx4minus2 = fDx4 - fDx2 ; 84 fDx3plus1 = fDx3 + fDx1 ; 95 fDx3plus1 = fDx3 + fDx1 ; 85 fDx3minus1 = fDx3 - fDx1 ; 96 fDx3minus1 = fDx3 - fDx1 ; 86 fDy2plus1 = fDy2 + fDy1 ; 97 fDy2plus1 = fDy2 + fDy1 ; 87 fDy2minus1 = fDy2 - fDy1 ; 98 fDy2minus1 = fDy2 - fDy1 ; 88 99 89 fa1md1 = 2*fDx2 - 2*fDx1 ; 100 fa1md1 = 2*fDx2 - 2*fDx1 ; 90 fa2md2 = 2*fDx4 - 2*fDx3 ; 101 fa2md2 = 2*fDx4 - 2*fDx3 ; 91 102 92 fPhiTwist = PhiTwist ; // dphi << 103 fPhiTwist = PhiTwist ; // dphi 93 fAngleSide = AngleSide ; // 0,90,180,270 de 104 fAngleSide = AngleSide ; // 0,90,180,270 deg 94 105 95 fdeltaX = 2*fDz*std::tan(fTheta)*std::cos(fP << 106 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ; // dx in surface equation 96 fdeltaY = 2*fDz*std::tan(fTheta)*std::sin(fP << 107 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ; // dy in surface equation 97 108 98 fRot.rotateZ( AngleSide ) ; 109 fRot.rotateZ( AngleSide ) ; 99 110 100 fTrans.set(0, 0, 0); // No Translation 111 fTrans.set(0, 0, 0); // No Translation 101 fIsValidNorm = false; 112 fIsValidNorm = false; 102 113 103 SetCorners() ; 114 SetCorners() ; 104 SetBoundaries() ; 115 SetBoundaries() ; >> 116 105 } 117 } 106 118 >> 119 107 //============================================ 120 //===================================================================== 108 //* Fake default constructor ----------------- 121 //* Fake default constructor ------------------------------------------ 109 122 110 G4TwistTrapParallelSide::G4TwistTrapParallelSi 123 G4TwistTrapParallelSide::G4TwistTrapParallelSide( __void__& a ) 111 : G4VTwistSurface(a) << 124 : G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.), >> 125 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.), >> 126 fAngleSide(0.), fdeltaX(0.), fdeltaY(0.), fDx4plus2(0.), fDx4minus2(0.), >> 127 fDx3plus1(0.), fDx3minus1(0.), fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.), >> 128 fa2md2(0.), fDy(0.) 112 { 129 { 113 } 130 } 114 131 >> 132 115 //============================================ 133 //===================================================================== 116 //* destructor ------------------------------- 134 //* destructor -------------------------------------------------------- 117 135 118 G4TwistTrapParallelSide::~G4TwistTrapParallelS << 136 G4TwistTrapParallelSide::~G4TwistTrapParallelSide() >> 137 { >> 138 } 119 139 120 //============================================ 140 //===================================================================== 121 //* GetNormal -------------------------------- 141 //* GetNormal --------------------------------------------------------- 122 142 123 G4ThreeVector G4TwistTrapParallelSide::GetNorm << 143 G4ThreeVector G4TwistTrapParallelSide::GetNormal(const G4ThreeVector &tmpxx, 124 << 144 G4bool isGlobal) 125 { 145 { 126 // GetNormal returns a normal vector at a s 146 // GetNormal returns a normal vector at a surface (or very close 127 // to surface) point at tmpxx. 147 // to surface) point at tmpxx. 128 // If isGlobal=true, it returns the normal 148 // If isGlobal=true, it returns the normal in global coordinate. 129 // 149 // 130 150 131 G4ThreeVector xx; 151 G4ThreeVector xx; 132 if (isGlobal) << 152 if (isGlobal) { 133 { << 134 xx = ComputeLocalPoint(tmpxx); 153 xx = ComputeLocalPoint(tmpxx); 135 if ((xx - fCurrentNormal.p).mag() < 0.5 << 154 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) { 136 { << 137 return ComputeGlobalDirection(fCurren 155 return ComputeGlobalDirection(fCurrentNormal.normal); 138 } 156 } 139 } << 157 } else { 140 else << 141 { << 142 xx = tmpxx; 158 xx = tmpxx; 143 if (xx == fCurrentNormal.p) << 159 if (xx == fCurrentNormal.p) { 144 { << 145 return fCurrentNormal.normal; 160 return fCurrentNormal.normal; 146 } 161 } 147 } 162 } 148 163 149 G4double phi ; 164 G4double phi ; 150 G4double u ; 165 G4double u ; 151 166 152 GetPhiUAtX(xx,phi,u) ; // phi,u for point 167 GetPhiUAtX(xx,phi,u) ; // phi,u for point xx close to surface 153 168 154 G4ThreeVector normal = NormAng(phi,u) ; / 169 G4ThreeVector normal = NormAng(phi,u) ; // the normal vector at phi,u 155 170 156 #ifdef G4TWISTDEBUG 171 #ifdef G4TWISTDEBUG 157 G4cout << "normal vector = " << normal << 172 G4cout << "normal vector = " << normal << G4endl ; 158 G4cout << "phi = " << phi << " , u = " << u 173 G4cout << "phi = " << phi << " , u = " << u << G4endl ; 159 #endif 174 #endif 160 175 161 // normal = normal/normal.mag() ; 176 // normal = normal/normal.mag() ; 162 177 163 if (isGlobal) << 178 if (isGlobal) { 164 { << 165 fCurrentNormal.normal = ComputeGlobalDir 179 fCurrentNormal.normal = ComputeGlobalDirection(normal.unit()); 166 } << 180 } else { 167 else << 168 { << 169 fCurrentNormal.normal = normal.unit(); 181 fCurrentNormal.normal = normal.unit(); 170 } 182 } 171 return fCurrentNormal.normal; 183 return fCurrentNormal.normal; 172 } 184 } 173 185 174 //============================================ 186 //===================================================================== 175 //* DistanceToSurface ------------------------ 187 //* DistanceToSurface ------------------------------------------------- 176 188 177 G4int G4TwistTrapParallelSide::DistanceToSurfa << 189 G4int G4TwistTrapParallelSide::DistanceToSurface(const G4ThreeVector &gp, 178 << 190 const G4ThreeVector &gv, 179 << 191 G4ThreeVector gxx[], 180 << 192 G4double distance[], 181 << 193 G4int areacode[], 182 << 194 G4bool isvalid[], 183 << 195 EValidate validate) 184 { 196 { >> 197 >> 198 static const G4double ctol = 0.5 * kCarTolerance; 185 static const G4double pihalf = pi/2 ; 199 static const G4double pihalf = pi/2 ; 186 const G4double ctol = 0.5 * kCarTolerance; << 187 200 188 G4bool IsParallel = false ; 201 G4bool IsParallel = false ; 189 G4bool IsConverged = false ; 202 G4bool IsConverged = false ; 190 203 191 G4int nxx = 0 ; // number of physical solut 204 G4int nxx = 0 ; // number of physical solutions 192 205 193 fCurStatWithV.ResetfDone(validate, &gp, &gv) 206 fCurStatWithV.ResetfDone(validate, &gp, &gv); 194 207 195 if (fCurStatWithV.IsDone()) << 208 if (fCurStatWithV.IsDone()) { 196 { << 209 G4int i; 197 for (G4int i=0; i<fCurStatWithV.GetNXX(); << 210 for (i=0; i<fCurStatWithV.GetNXX(); i++) { 198 { << 199 gxx[i] = fCurStatWithV.GetXX(i); 211 gxx[i] = fCurStatWithV.GetXX(i); 200 distance[i] = fCurStatWithV.GetDistance( 212 distance[i] = fCurStatWithV.GetDistance(i); 201 areacode[i] = fCurStatWithV.GetAreacode( 213 areacode[i] = fCurStatWithV.GetAreacode(i); 202 isvalid[i] = fCurStatWithV.IsValid(i); 214 isvalid[i] = fCurStatWithV.IsValid(i); 203 } 215 } 204 return fCurStatWithV.GetNXX(); 216 return fCurStatWithV.GetNXX(); 205 } << 217 } else { 206 else // initialize << 218 207 { << 219 // initialize 208 for (G4int i=0; i<G4VSURFACENXX ; ++i) << 220 G4int i; 209 { << 221 for (i=0; i<G4VSURFACENXX ; i++) { 210 distance[i] = kInfinity; 222 distance[i] = kInfinity; 211 areacode[i] = sOutside; 223 areacode[i] = sOutside; 212 isvalid[i] = false; 224 isvalid[i] = false; 213 gxx[i].set(kInfinity, kInfinity, kInfini 225 gxx[i].set(kInfinity, kInfinity, kInfinity); 214 } 226 } 215 } 227 } 216 228 217 G4ThreeVector p = ComputeLocalPoint(gp); 229 G4ThreeVector p = ComputeLocalPoint(gp); 218 G4ThreeVector v = ComputeLocalDirection(gv); 230 G4ThreeVector v = ComputeLocalDirection(gv); 219 231 220 #ifdef G4TWISTDEBUG 232 #ifdef G4TWISTDEBUG 221 G4cout << "Local point p = " << p << G4endl 233 G4cout << "Local point p = " << p << G4endl ; 222 G4cout << "Local direction v = " << v << G4e 234 G4cout << "Local direction v = " << v << G4endl ; 223 #endif 235 #endif 224 236 225 G4double phi,u ; // parameters 237 G4double phi,u ; // parameters 226 238 227 // temporary variables 239 // temporary variables 228 240 229 G4double tmpdist = kInfinity ; 241 G4double tmpdist = kInfinity ; 230 G4ThreeVector tmpxx; 242 G4ThreeVector tmpxx; 231 G4int tmpareacode = sOutside ; 243 G4int tmpareacode = sOutside ; 232 G4bool tmpisvalid = false ; 244 G4bool tmpisvalid = false ; 233 245 234 std::vector<Intersection> xbuf ; 246 std::vector<Intersection> xbuf ; 235 Intersection xbuftmp ; 247 Intersection xbuftmp ; 236 248 237 // prepare some variables for the intersecti 249 // prepare some variables for the intersection finder 238 250 239 G4double L = 2*fDz ; 251 G4double L = 2*fDz ; 240 252 241 G4double phixz = fPhiTwist * ( p.x() * v.z() 253 G4double phixz = fPhiTwist * ( p.x() * v.z() - p.z() * v.x() ) ; 242 G4double phiyz = fPhiTwist * ( p.y() * v.z() 254 G4double phiyz = fPhiTwist * ( p.y() * v.z() - p.z() * v.y() ) ; 243 255 244 // special case vz = 0 256 // special case vz = 0 245 257 246 if ( v.z() == 0. ) << 258 if ( v.z() == 0. ) { 247 { << 259 248 if ( std::fabs(p.z()) <= L ) // inte << 260 if ( std::fabs(p.z()) <= L ) { // intersection possible in z 249 { << 261 250 phi = p.z() * fPhiTwist / L ; // phi is 262 phi = p.z() * fPhiTwist / L ; // phi is determined by the z-position 251 263 252 u = (2*(fdeltaY*phi*v.x() - fPhiTwist*p. << 264 u = (2*(fdeltaY*phi*v.x() - fPhiTwist*p.y()*v.x() - fdeltaX*phi*v.y() + fPhiTwist*p.x()*v.y()) + (fDy2plus1*fPhiTwist + 2*fDy2minus1*phi)*(v.x()*std::cos(phi) + v.y()*std::sin(phi)))/(2.* fPhiTwist*(v.y()*std::cos(phi) - v.x()*std::sin(phi))) ; 253 + fPhiTwist*p.x()*v.y()) + (fDy2 << 254 + 2*fDy2minus1*phi)*(v.x()*std:: << 255 / (2.* fPhiTwist*(v.y()*std::cos(phi) << 256 265 257 xbuftmp.phi = phi ; 266 xbuftmp.phi = phi ; 258 xbuftmp.u = u ; 267 xbuftmp.u = u ; 259 xbuftmp.areacode = sOutside ; 268 xbuftmp.areacode = sOutside ; 260 xbuftmp.distance = kInfinity ; 269 xbuftmp.distance = kInfinity ; 261 xbuftmp.isvalid = false ; 270 xbuftmp.isvalid = false ; 262 271 263 xbuf.push_back(xbuftmp) ; // store it to << 272 xbuf.push_back(xbuftmp) ; // store it to xbuf >> 273 264 } 274 } 265 else // no in << 275 266 { << 276 else { // no intersection possible >> 277 267 distance[0] = kInfinity; 278 distance[0] = kInfinity; 268 gxx[0].set(kInfinity,kInfinity,kInfinity 279 gxx[0].set(kInfinity,kInfinity,kInfinity); 269 isvalid[0] = false ; 280 isvalid[0] = false ; 270 areacode[0] = sOutside ; 281 areacode[0] = sOutside ; 271 fCurStatWithV.SetCurrentStatus(0, gxx[0] 282 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], 272 areacode[ 283 areacode[0], isvalid[0], 273 0, valida 284 0, validate, &gp, &gv); 274 285 275 return 0; 286 return 0; >> 287 >> 288 276 } // end std::fabs(p.z() <= L 289 } // end std::fabs(p.z() <= L >> 290 277 } // end v.z() == 0 291 } // end v.z() == 0 278 else // general solution for non-zero vz << 292 279 { << 293 >> 294 // general solution for non-zero vz >> 295 >> 296 else { >> 297 280 G4double c[9],srd[8],si[8] ; 298 G4double c[9],srd[8],si[8] ; 281 299 282 c[8] = -3600*(-2*phiyz + fDy2plus1*fPhiTwi 300 c[8] = -3600*(-2*phiyz + fDy2plus1*fPhiTwist*v.z()) ; 283 c[7] = -7200*(phixz - 2*fDz*v.y() + (fdelt 301 c[7] = -7200*(phixz - 2*fDz*v.y() + (fdeltaY + fDy2minus1)*v.z()) ; 284 c[6] = 120*(-52*phiyz - 120*fDz*v.x() + 60 << 302 c[6] = 120*(-52*phiyz - 120*fDz*v.x() + 60*fdeltaX*v.z() + 11*fDy2plus1*fPhiTwist*v.z()) ; 285 + 11*fDy2plus1*fPhiTwist*v.z()) ; << 303 c[5] = 240*(16*phixz - 52*fDz*v.y() + 26*fdeltaY*v.z() + 11*fDy2minus1*v.z()) ; 286 c[5] = 240*(16*phixz - 52*fDz*v.y() + 26*f << 304 c[4] = 12*(127*phiyz + 640*fDz*v.x() - 320*fdeltaX*v.z() + 4*fDy2plus1*fPhiTwist*v.z()) ; 287 + 11*fDy2minus1*v.z()) ; << 305 c[3] = -404*phixz + 3048*fDz*v.y() - 1524*fdeltaY*v.z() + 96*fDy2minus1*v.z() ; 288 c[4] = 12*(127*phiyz + 640*fDz*v.x() - 320 << 289 + 4*fDy2plus1*fPhiTwist*v.z()) ; << 290 c[3] = -404*phixz + 3048*fDz*v.y() - 1524* << 291 + 96*fDy2minus1*v.z() ; << 292 c[2] = -72*phiyz + 404*(-2*fDz*v.x() + fde 306 c[2] = -72*phiyz + 404*(-2*fDz*v.x() + fdeltaX*v.z()) ; 293 c[1] = 12*(phixz - 12*fDz*v.y() + 6*fdelta 307 c[1] = 12*(phixz - 12*fDz*v.y() + 6*fdeltaY*v.z()) ; 294 c[0] = 24*fDz*v.x() - 12*fdeltaX*v.z() ; 308 c[0] = 24*fDz*v.x() - 12*fdeltaX*v.z() ; 295 309 296 310 297 #ifdef G4TWISTDEBUG 311 #ifdef G4TWISTDEBUG 298 G4cout << "coef = " << c[0] << " " 312 G4cout << "coef = " << c[0] << " " 299 << c[1] << " " 313 << c[1] << " " 300 << c[2] << " " 314 << c[2] << " " 301 << c[3] << " " 315 << c[3] << " " 302 << c[4] << " " 316 << c[4] << " " 303 << c[5] << " " 317 << c[5] << " " 304 << c[6] << " " 318 << c[6] << " " 305 << c[7] << " " 319 << c[7] << " " 306 << c[8] << G4endl ; 320 << c[8] << G4endl ; 307 #endif 321 #endif 308 322 309 G4JTPolynomialSolver trapEq ; 323 G4JTPolynomialSolver trapEq ; 310 G4int num = trapEq.FindRoots(c,8,srd,si); 324 G4int num = trapEq.FindRoots(c,8,srd,si); >> 325 311 326 312 for (G4int i = 0 ; i<num ; ++i ) // loop << 327 for (G4int i = 0 ; i<num ; i++ ) { // loop over all mathematical solutions 313 { << 328 if ( si[i]==0.0 ) { // only real solutions 314 if ( si[i]==0.0 ) // only real solutio << 315 { << 316 #ifdef G4TWISTDEBUG 329 #ifdef G4TWISTDEBUG 317 G4cout << "Solution " << i << " : " << 330 G4cout << "Solution " << i << " : " << srd[i] << G4endl ; 318 #endif 331 #endif 319 phi = std::fmod(srd[i] , pihalf) ; 332 phi = std::fmod(srd[i] , pihalf) ; 320 u = (1/std::cos(phi)*(2*phixz + 4*fDz* << 333 321 - 2*fdeltaX*phi*v.z() + (fDy2plus1*f << 334 u = (1/std::cos(phi)*(2*phixz + 4*fDz*phi*v.x() - 2*fdeltaX*phi*v.z() + (fDy2plus1*fPhiTwist + 2*fDy2minus1*phi)*v.z()* std::sin(phi)))/(2.*fPhiTwist*v.z()) ; 322 + 2*fDy2minus1*phi)*v.z()* std::sin( << 323 335 324 xbuftmp.phi = phi ; 336 xbuftmp.phi = phi ; 325 xbuftmp.u = u ; 337 xbuftmp.u = u ; 326 xbuftmp.areacode = sOutside ; 338 xbuftmp.areacode = sOutside ; 327 xbuftmp.distance = kInfinity ; 339 xbuftmp.distance = kInfinity ; 328 xbuftmp.isvalid = false ; 340 xbuftmp.isvalid = false ; 329 341 330 xbuf.push_back(xbuftmp) ; // store it 342 xbuf.push_back(xbuftmp) ; // store it to xbuf 331 343 332 #ifdef G4TWISTDEBUG 344 #ifdef G4TWISTDEBUG 333 G4cout << "solution " << i << " = " << 345 G4cout << "solution " << i << " = " << phi << " , " << u << G4endl ; 334 #endif 346 #endif 335 347 336 } // end if real solution 348 } // end if real solution 337 } // end loop i 349 } // end loop i 338 } // end general case << 350 >> 351 } // end general case >> 352 339 353 340 nxx = (G4int)xbuf.size() ; // save the numb << 354 nxx = xbuf.size() ; // save the number of solutions 341 355 342 G4ThreeVector xxonsurface ; // point 356 G4ThreeVector xxonsurface ; // point on surface 343 G4ThreeVector surfacenormal ; // normal 357 G4ThreeVector surfacenormal ; // normal vector 344 G4double deltaX ; // distance between inter << 358 G4double deltaX ; // distance between intersection point and point on surface 345 G4double theta ; // angle 359 G4double theta ; // angle between track and surfacenormal 346 G4double factor ; // a scal 360 G4double factor ; // a scaling factor 347 G4int maxint = 30 ; // number 361 G4int maxint = 30 ; // number of iterations 348 362 349 363 350 for (auto & k : xbuf) << 364 for ( size_t k = 0 ; k<xbuf.size() ; k++ ) { 351 { << 365 352 #ifdef G4TWISTDEBUG 366 #ifdef G4TWISTDEBUG 353 G4cout << "Solution " << k << " : " 367 G4cout << "Solution " << k << " : " 354 << "reconstructed phiR = " << xbuf[ 368 << "reconstructed phiR = " << xbuf[k].phi 355 << ", uR = " << xbuf[k].u << G4endl 369 << ", uR = " << xbuf[k].u << G4endl ; 356 #endif 370 #endif 357 371 358 phi = k.phi ; // get the stored values fo << 372 phi = xbuf[k].phi ; // get the stored values for phi and u 359 u = k.u ; << 373 u = xbuf[k].u ; 360 374 361 IsConverged = false ; // no convergence 375 IsConverged = false ; // no convergence at the beginning 362 376 363 for ( G4int i = 1 ; i<maxint ; ++i ) << 377 for ( G4int i = 1 ; i<maxint ; i++ ) { 364 { << 378 365 xxonsurface = SurfacePoint(phi,u) ; 379 xxonsurface = SurfacePoint(phi,u) ; 366 surfacenormal = NormAng(phi,u) ; 380 surfacenormal = NormAng(phi,u) ; 367 tmpdist = DistanceToPlaneWithV(p, v, xxo 381 tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 368 deltaX = ( tmpxx - xxonsurface ).mag() ; 382 deltaX = ( tmpxx - xxonsurface ).mag() ; 369 theta = std::fabs(std::acos(v*surfacenor 383 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ; 370 if ( theta < 0.001 ) << 384 if ( theta < 0.001 ) { 371 { << 372 factor = 50 ; 385 factor = 50 ; 373 IsParallel = true ; 386 IsParallel = true ; 374 } 387 } 375 else << 388 else { 376 { << 377 factor = 1 ; 389 factor = 1 ; 378 } 390 } 379 391 380 #ifdef G4TWISTDEBUG 392 #ifdef G4TWISTDEBUG 381 G4cout << "Step i = " << i << ", distanc << 393 G4cout << "Step i = " << i << ", distance = " << tmpdist << ", " << deltaX << G4endl ; 382 << tmpdist << ", " << deltaX << G << 383 G4cout << "X = " << tmpxx << G4endl ; 394 G4cout << "X = " << tmpxx << G4endl ; 384 #endif 395 #endif 385 396 386 GetPhiUAtX(tmpxx, phi, u); // new point << 397 GetPhiUAtX(tmpxx, phi, u) ; // the new point xx is accepted and phi/u replaced 387 398 388 #ifdef G4TWISTDEBUG 399 #ifdef G4TWISTDEBUG 389 G4cout << "approximated phi = " << phi < 400 G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 390 #endif 401 #endif 391 402 392 if ( deltaX <= factor*ctol ) { IsConverg 403 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; } >> 404 393 } // end iterative loop (i) 405 } // end iterative loop (i) 394 406 >> 407 >> 408 // new code 21.09.05 O.Link 395 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 409 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ; 396 410 >> 411 397 #ifdef G4TWISTDEBUG 412 #ifdef G4TWISTDEBUG 398 G4cout << "refined solution " << phi << " 413 G4cout << "refined solution " << phi << " , " << u << G4endl ; 399 G4cout << "distance = " << tmpdist << G4en 414 G4cout << "distance = " << tmpdist << G4endl ; 400 G4cout << "local X = " << tmpxx << G4endl 415 G4cout << "local X = " << tmpxx << G4endl ; 401 #endif 416 #endif 402 417 403 tmpisvalid = false ; // init 418 tmpisvalid = false ; // init 404 419 405 if ( IsConverged ) << 420 if ( IsConverged ) { 406 { << 421 407 if (validate == kValidateWithTol) << 422 if (validate == kValidateWithTol) { 408 { << 409 tmpareacode = GetAreaCode(tmpxx); 423 tmpareacode = GetAreaCode(tmpxx); 410 if (!IsOutside(tmpareacode)) << 424 if (!IsOutside(tmpareacode)) { 411 { << 412 if (tmpdist >= 0) tmpisvalid = true; 425 if (tmpdist >= 0) tmpisvalid = true; 413 } 426 } 414 } << 427 } else if (validate == kValidateWithoutTol) { 415 else if (validate == kValidateWithoutTol << 416 { << 417 tmpareacode = GetAreaCode(tmpxx, false 428 tmpareacode = GetAreaCode(tmpxx, false); 418 if (IsInside(tmpareacode)) << 429 if (IsInside(tmpareacode)) { 419 { << 420 if (tmpdist >= 0) tmpisvalid = true; 430 if (tmpdist >= 0) tmpisvalid = true; 421 } 431 } 422 } << 432 } else { // kDontValidate 423 else // kDontValidate << 424 { << 425 G4Exception("G4TwistTrapParallelSide:: 433 G4Exception("G4TwistTrapParallelSide::DistanceToSurface()", 426 "GeomSolids0001", FatalExc 434 "GeomSolids0001", FatalException, 427 "Feature NOT implemented ! 435 "Feature NOT implemented !"); 428 } 436 } >> 437 429 } 438 } 430 else << 439 else { 431 { << 432 tmpdist = kInfinity; // no convergen 440 tmpdist = kInfinity; // no convergence after 10 steps 433 tmpisvalid = false ; // solution is 441 tmpisvalid = false ; // solution is not vaild 434 } 442 } 435 443 >> 444 436 // store the found values 445 // store the found values 437 k.xx = tmpxx ; << 446 xbuf[k].xx = tmpxx ; 438 k.distance = tmpdist ; << 447 xbuf[k].distance = tmpdist ; 439 k.areacode = tmpareacode ; << 448 xbuf[k].areacode = tmpareacode ; 440 k.isvalid = tmpisvalid ; << 449 xbuf[k].isvalid = tmpisvalid ; >> 450 >> 451 441 } // end loop over physical solutions (vari 452 } // end loop over physical solutions (variable k) 442 453 >> 454 443 std::sort(xbuf.begin() , xbuf.end(), Distanc 455 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting 444 456 445 #ifdef G4TWISTDEBUG 457 #ifdef G4TWISTDEBUG 446 G4cout << G4endl << "list xbuf after sorting 458 G4cout << G4endl << "list xbuf after sorting : " << G4endl ; 447 G4cout << G4endl << G4endl ; 459 G4cout << G4endl << G4endl ; 448 #endif 460 #endif 449 461 >> 462 450 // erase identical intersection (within kCar 463 // erase identical intersection (within kCarTolerance) 451 xbuf.erase(std::unique(xbuf.begin(),xbuf.end << 464 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) , xbuf.end() ) ; 452 465 453 466 454 // add guesses 467 // add guesses 455 468 456 auto nxxtmp = (G4int)xbuf.size() ; << 469 G4int nxxtmp = xbuf.size() ; >> 470 >> 471 if ( nxxtmp<2 || IsParallel ) { 457 472 458 if ( nxxtmp<2 || IsParallel ) << 459 { << 460 // positive end 473 // positive end 461 #ifdef G4TWISTDEBUG 474 #ifdef G4TWISTDEBUG 462 G4cout << "add guess at +z/2 .. " << G4end 475 G4cout << "add guess at +z/2 .. " << G4endl ; 463 #endif 476 #endif 464 477 465 phi = fPhiTwist/2 ; 478 phi = fPhiTwist/2 ; 466 u = 0 ; 479 u = 0 ; 467 480 468 xbuftmp.phi = phi ; 481 xbuftmp.phi = phi ; 469 xbuftmp.u = u ; 482 xbuftmp.u = u ; 470 xbuftmp.areacode = sOutside ; 483 xbuftmp.areacode = sOutside ; 471 xbuftmp.distance = kInfinity ; 484 xbuftmp.distance = kInfinity ; 472 xbuftmp.isvalid = false ; 485 xbuftmp.isvalid = false ; 473 486 474 xbuf.push_back(xbuftmp) ; // store it to 487 xbuf.push_back(xbuftmp) ; // store it to xbuf 475 488 >> 489 476 #ifdef G4TWISTDEBUG 490 #ifdef G4TWISTDEBUG 477 G4cout << "add guess at -z/2 .. " << G4end 491 G4cout << "add guess at -z/2 .. " << G4endl ; 478 #endif 492 #endif 479 493 480 phi = -fPhiTwist/2 ; 494 phi = -fPhiTwist/2 ; 481 u = 0 ; 495 u = 0 ; 482 496 483 xbuftmp.phi = phi ; 497 xbuftmp.phi = phi ; 484 xbuftmp.u = u ; 498 xbuftmp.u = u ; 485 xbuftmp.areacode = sOutside ; 499 xbuftmp.areacode = sOutside ; 486 xbuftmp.distance = kInfinity ; 500 xbuftmp.distance = kInfinity ; 487 xbuftmp.isvalid = false ; 501 xbuftmp.isvalid = false ; 488 502 489 xbuf.push_back(xbuftmp) ; // store it to 503 xbuf.push_back(xbuftmp) ; // store it to xbuf 490 504 491 for ( std::size_t k = nxxtmp ; k<xbuf.size << 505 for ( size_t k = nxxtmp ; k<xbuf.size() ; k++ ) { 492 { << 506 493 #ifdef G4TWISTDEBUG 507 #ifdef G4TWISTDEBUG 494 G4cout << "Solution " << k << " : " 508 G4cout << "Solution " << k << " : " 495 << "reconstructed phiR = " << xbu 509 << "reconstructed phiR = " << xbuf[k].phi 496 << ", uR = " << xbuf[k].u << G4en 510 << ", uR = " << xbuf[k].u << G4endl ; 497 #endif 511 #endif 498 512 499 phi = xbuf[k].phi ; // get the stored v 513 phi = xbuf[k].phi ; // get the stored values for phi and u 500 u = xbuf[k].u ; 514 u = xbuf[k].u ; 501 515 502 IsConverged = false ; // no convergenc 516 IsConverged = false ; // no convergence at the beginning 503 517 504 for ( G4int i = 1 ; i<maxint ; ++i ) << 518 for ( G4int i = 1 ; i<maxint ; i++ ) { 505 { << 519 506 xxonsurface = SurfacePoint(phi,u) ; 520 xxonsurface = SurfacePoint(phi,u) ; 507 surfacenormal = NormAng(phi,u) ; 521 surfacenormal = NormAng(phi,u) ; 508 tmpdist = DistanceToPlaneWithV(p, v, x 522 tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 509 deltaX = ( tmpxx - xxonsurface ).mag() 523 deltaX = ( tmpxx - xxonsurface ).mag() ; 510 theta = std::fabs(std::acos(v*surfacen 524 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ; 511 if ( theta < 0.001 ) << 525 if ( theta < 0.001 ) { 512 { << 513 factor = 50 ; 526 factor = 50 ; 514 } 527 } 515 else << 528 else { 516 { << 517 factor = 1 ; 529 factor = 1 ; 518 } 530 } 519 531 520 #ifdef G4TWISTDEBUG 532 #ifdef G4TWISTDEBUG 521 G4cout << "Step i = " << i << ", dista << 533 G4cout << "Step i = " << i << ", distance = " << tmpdist << ", " << deltaX << G4endl ; 522 << tmpdist << ", " << deltaX << << 523 G4cout << "X = " << tmpxx << G4endl ; 534 G4cout << "X = " << tmpxx << G4endl ; 524 #endif 535 #endif 525 536 526 GetPhiUAtX(tmpxx, phi, u) ; // new poi << 537 GetPhiUAtX(tmpxx, phi, u) ; // the new point xx is accepted and phi/u replaced 527 << 538 528 #ifdef G4TWISTDEBUG 539 #ifdef G4TWISTDEBUG 529 G4cout << "approximated phi = " << phi 540 G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 530 #endif 541 #endif 531 542 532 if ( deltaX <= factor*ctol ) { IsConve 543 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; } >> 544 533 } // end iterative loop (i) 545 } // end iterative loop (i) >> 546 534 547 >> 548 // new code 21.09.05 O.Link 535 if ( std::fabs(tmpdist)<ctol ) tmpdist = 549 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ; 536 550 >> 551 537 #ifdef G4TWISTDEBUG 552 #ifdef G4TWISTDEBUG 538 G4cout << "refined solution " << phi << 553 G4cout << "refined solution " << phi << " , " << u << G4endl ; 539 G4cout << "distance = " << tmpdist << G4 554 G4cout << "distance = " << tmpdist << G4endl ; 540 G4cout << "local X = " << tmpxx << G4end 555 G4cout << "local X = " << tmpxx << G4endl ; 541 #endif 556 #endif 542 557 543 tmpisvalid = false ; // init 558 tmpisvalid = false ; // init 544 559 545 if ( IsConverged ) << 560 if ( IsConverged ) { 546 { << 561 547 if (validate == kValidateWithTol) << 562 if (validate == kValidateWithTol) { 548 { << 549 tmpareacode = GetAreaCode(tmpxx); 563 tmpareacode = GetAreaCode(tmpxx); 550 if (!IsOutside(tmpareacode)) << 564 if (!IsOutside(tmpareacode)) { 551 { << 552 if (tmpdist >= 0) tmpisvalid = tru 565 if (tmpdist >= 0) tmpisvalid = true; 553 } 566 } 554 } << 567 } else if (validate == kValidateWithoutTol) { 555 else if (validate == kValidateWithoutT << 556 { << 557 tmpareacode = GetAreaCode(tmpxx, fal 568 tmpareacode = GetAreaCode(tmpxx, false); 558 if (IsInside(tmpareacode)) << 569 if (IsInside(tmpareacode)) { 559 { << 560 if (tmpdist >= 0) tmpisvalid = tru 570 if (tmpdist >= 0) tmpisvalid = true; 561 } 571 } 562 } << 572 } else { // kDontValidate 563 else // kDontValidate << 564 { << 565 G4Exception("G4TwistedBoxSide::Dista 573 G4Exception("G4TwistedBoxSide::DistanceToSurface()", 566 "GeomSolids0001", FatalE 574 "GeomSolids0001", FatalException, 567 "Feature NOT implemented 575 "Feature NOT implemented !"); 568 } 576 } 569 577 570 } 578 } 571 else << 579 else { 572 { << 573 tmpdist = kInfinity; // no converg 580 tmpdist = kInfinity; // no convergence after 10 steps 574 tmpisvalid = false ; // solution i 581 tmpisvalid = false ; // solution is not vaild 575 } 582 } 576 << 583 >> 584 577 // store the found values 585 // store the found values 578 xbuf[k].xx = tmpxx ; 586 xbuf[k].xx = tmpxx ; 579 xbuf[k].distance = tmpdist ; 587 xbuf[k].distance = tmpdist ; 580 xbuf[k].areacode = tmpareacode ; 588 xbuf[k].areacode = tmpareacode ; 581 xbuf[k].isvalid = tmpisvalid ; 589 xbuf[k].isvalid = tmpisvalid ; 582 590 >> 591 583 } // end loop over physical solutions 592 } // end loop over physical solutions >> 593 >> 594 584 } // end less than 2 solutions 595 } // end less than 2 solutions 585 596 >> 597 586 // sort again 598 // sort again 587 std::sort(xbuf.begin() , xbuf.end(), Distanc 599 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting 588 600 589 // erase identical intersection (within kCar 601 // erase identical intersection (within kCarTolerance) 590 xbuf.erase(std::unique(xbuf.begin(),xbuf.end << 602 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) , xbuf.end() ) ; 591 603 592 #ifdef G4TWISTDEBUG 604 #ifdef G4TWISTDEBUG 593 G4cout << G4endl << "list xbuf after sorting 605 G4cout << G4endl << "list xbuf after sorting : " << G4endl ; 594 G4cout << G4endl << G4endl ; 606 G4cout << G4endl << G4endl ; 595 #endif 607 #endif 596 608 597 nxx = (G4int)xbuf.size() ; // determine nu << 609 nxx = xbuf.size() ; // determine number of solutions again. 598 610 599 for ( G4int i = 0 ; i<(G4int)xbuf.size() ; + << 611 for ( size_t i = 0 ; i<xbuf.size() ; i++ ) { 600 { << 612 601 distance[i] = xbuf[i].distance; 613 distance[i] = xbuf[i].distance; 602 gxx[i] = ComputeGlobalPoint(xbuf[i].x 614 gxx[i] = ComputeGlobalPoint(xbuf[i].xx); 603 areacode[i] = xbuf[i].areacode ; 615 areacode[i] = xbuf[i].areacode ; 604 isvalid[i] = xbuf[i].isvalid ; 616 isvalid[i] = xbuf[i].isvalid ; 605 617 606 fCurStatWithV.SetCurrentStatus(i, gxx[i], 618 fCurStatWithV.SetCurrentStatus(i, gxx[i], distance[i], areacode[i], 607 isvalid[i 619 isvalid[i], nxx, validate, &gp, &gv); >> 620 608 #ifdef G4TWISTDEBUG 621 #ifdef G4TWISTDEBUG 609 G4cout << "element Nr. " << i 622 G4cout << "element Nr. " << i 610 << ", local Intersection = " << xbu 623 << ", local Intersection = " << xbuf[i].xx 611 << ", distance = " << xbuf[i].dista 624 << ", distance = " << xbuf[i].distance 612 << ", u = " << xbuf[i].u 625 << ", u = " << xbuf[i].u 613 << ", phi = " << xbuf[i].phi 626 << ", phi = " << xbuf[i].phi 614 << ", isvalid = " << xbuf[i].isvali 627 << ", isvalid = " << xbuf[i].isvalid 615 << G4endl ; 628 << G4endl ; 616 #endif 629 #endif >> 630 617 } // end for( i ) loop 631 } // end for( i ) loop 618 632 >> 633 619 #ifdef G4TWISTDEBUG 634 #ifdef G4TWISTDEBUG 620 G4cout << "G4TwistTrapParallelSide finished 635 G4cout << "G4TwistTrapParallelSide finished " << G4endl ; 621 G4cout << nxx << " possible physical solutio 636 G4cout << nxx << " possible physical solutions found" << G4endl ; 622 for ( G4int k= 0 ; k< nxx ; ++k ) << 637 for ( G4int k= 0 ; k< nxx ; k++ ) { 623 { << 624 G4cout << "global intersection Point found 638 G4cout << "global intersection Point found: " << gxx[k] << G4endl ; 625 G4cout << "distance = " << distance[k] << 639 G4cout << "distance = " << distance[k] << G4endl ; 626 G4cout << "isvalid = " << isvalid[k] << G4 640 G4cout << "isvalid = " << isvalid[k] << G4endl ; 627 } 641 } 628 #endif 642 #endif 629 643 630 return nxx ; 644 return nxx ; >> 645 631 } 646 } 632 647 >> 648 >> 649 633 //============================================ 650 //===================================================================== 634 //* DistanceToSurface ------------------------ 651 //* DistanceToSurface ------------------------------------------------- 635 652 636 G4int G4TwistTrapParallelSide::DistanceToSurfa << 653 G4int G4TwistTrapParallelSide::DistanceToSurface(const G4ThreeVector &gp, 637 << 654 G4ThreeVector gxx[], 638 << 655 G4double distance[], 639 << 656 G4int areacode[]) 640 { 657 { 641 const G4double ctol = 0.5 * kCarTolerance; << 658 // to do >> 659 >> 660 static const G4double ctol = 0.5 * kCarTolerance; 642 661 643 fCurStat.ResetfDone(kDontValidate, &gp); 662 fCurStat.ResetfDone(kDontValidate, &gp); 644 663 645 if (fCurStat.IsDone()) << 664 if (fCurStat.IsDone()) { 646 { << 665 G4int i; 647 for (G4int i=0; i<fCurStat.GetNXX(); ++i << 666 for (i=0; i<fCurStat.GetNXX(); i++) { 648 { << 649 gxx[i] = fCurStat.GetXX(i); 667 gxx[i] = fCurStat.GetXX(i); 650 distance[i] = fCurStat.GetDistance(i) 668 distance[i] = fCurStat.GetDistance(i); 651 areacode[i] = fCurStat.GetAreacode(i) 669 areacode[i] = fCurStat.GetAreacode(i); 652 } 670 } 653 return fCurStat.GetNXX(); 671 return fCurStat.GetNXX(); 654 } << 672 } else { 655 else // initialize << 673 // initialize 656 { << 674 G4int i; 657 for (G4int i=0; i<G4VSURFACENXX; ++i) << 675 for (i=0; i<G4VSURFACENXX; i++) { 658 { << 659 distance[i] = kInfinity; 676 distance[i] = kInfinity; 660 areacode[i] = sOutside; 677 areacode[i] = sOutside; 661 gxx[i].set(kInfinity, kInfinity, kInf 678 gxx[i].set(kInfinity, kInfinity, kInfinity); 662 } 679 } 663 } 680 } 664 681 665 G4ThreeVector p = ComputeLocalPoint(gp); 682 G4ThreeVector p = ComputeLocalPoint(gp); 666 G4ThreeVector xx; // intersection point 683 G4ThreeVector xx; // intersection point 667 G4ThreeVector xxonsurface ; // interpolated 684 G4ThreeVector xxonsurface ; // interpolated intersection point 668 685 669 // the surfacenormal at that surface point 686 // the surfacenormal at that surface point 670 G4double phiR = 0 ; // 687 G4double phiR = 0 ; // 671 G4double uR = 0 ; 688 G4double uR = 0 ; 672 689 673 G4ThreeVector surfacenormal ; 690 G4ThreeVector surfacenormal ; 674 G4double deltaX ; 691 G4double deltaX ; 675 692 676 G4int maxint = 20 ; 693 G4int maxint = 20 ; 677 694 678 for ( G4int i = 1 ; i<maxint ; ++i ) << 695 for ( G4int i = 1 ; i<maxint ; i++ ) { 679 { << 696 680 xxonsurface = SurfacePoint(phiR,uR) ; 697 xxonsurface = SurfacePoint(phiR,uR) ; 681 surfacenormal = NormAng(phiR,uR) ; 698 surfacenormal = NormAng(phiR,uR) ; 682 distance[0] = DistanceToPlane(p, xxonsurf 699 distance[0] = DistanceToPlane(p, xxonsurface, surfacenormal, xx); // new XX 683 deltaX = ( xx - xxonsurface ).mag() ; 700 deltaX = ( xx - xxonsurface ).mag() ; 684 701 685 #ifdef G4TWISTDEBUG 702 #ifdef G4TWISTDEBUG 686 G4cout << "i = " << i << ", distance = " << 703 G4cout << "i = " << i << ", distance = " << distance[0] << ", " << deltaX << G4endl ; 687 << distance[0] << ", " << deltaX < << 688 G4cout << "X = " << xx << G4endl ; 704 G4cout << "X = " << xx << G4endl ; 689 #endif 705 #endif 690 706 691 // the new point xx is accepted and phi/p 707 // the new point xx is accepted and phi/psi replaced 692 GetPhiUAtX(xx, phiR, uR) ; 708 GetPhiUAtX(xx, phiR, uR) ; 693 709 694 if ( deltaX <= ctol ) { break ; } 710 if ( deltaX <= ctol ) { break ; } >> 711 695 } 712 } 696 713 697 // check validity of solution ( valid phi,p 714 // check validity of solution ( valid phi,psi ) 698 715 699 G4double halfphi = 0.5*fPhiTwist ; 716 G4double halfphi = 0.5*fPhiTwist ; 700 G4double uMax = GetBoundaryMax(phiR) ; 717 G4double uMax = GetBoundaryMax(phiR) ; 701 G4double uMin = GetBoundaryMin(phiR) ; 718 G4double uMin = GetBoundaryMin(phiR) ; 702 719 703 if ( phiR > halfphi ) phiR = halfphi ; 720 if ( phiR > halfphi ) phiR = halfphi ; 704 if ( phiR < -halfphi ) phiR = -halfphi ; 721 if ( phiR < -halfphi ) phiR = -halfphi ; 705 if ( uR > uMax ) uR = uMax ; 722 if ( uR > uMax ) uR = uMax ; 706 if ( uR < uMin ) uR = uMin ; 723 if ( uR < uMin ) uR = uMin ; 707 724 708 xxonsurface = SurfacePoint(phiR,uR) ; 725 xxonsurface = SurfacePoint(phiR,uR) ; 709 distance[0] = ( p - xx ).mag() ; 726 distance[0] = ( p - xx ).mag() ; 710 if ( distance[0] <= ctol ) { distance[0] = 727 if ( distance[0] <= ctol ) { distance[0] = 0 ; } 711 728 712 // end of validity 729 // end of validity 713 730 714 #ifdef G4TWISTDEBUG 731 #ifdef G4TWISTDEBUG 715 G4cout << "refined solution " << phiR << " 732 G4cout << "refined solution " << phiR << " , " << uR << " , " << G4endl ; 716 G4cout << "distance = " << distance[0] << G 733 G4cout << "distance = " << distance[0] << G4endl ; 717 G4cout << "X = " << xx << G4endl ; 734 G4cout << "X = " << xx << G4endl ; 718 #endif 735 #endif 719 736 720 G4bool isvalid = true; 737 G4bool isvalid = true; 721 gxx[0] = ComputeGlobalPoint(xx); 738 gxx[0] = ComputeGlobalPoint(xx); 722 739 723 #ifdef G4TWISTDEBUG 740 #ifdef G4TWISTDEBUG 724 G4cout << "intersection Point found: " << g 741 G4cout << "intersection Point found: " << gxx[0] << G4endl ; 725 G4cout << "distance = " << distance[0] << G 742 G4cout << "distance = " << distance[0] << G4endl ; 726 #endif 743 #endif 727 744 728 fCurStat.SetCurrentStatus(0, gxx[0], distan 745 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], 729 isvalid, 1, kDontV 746 isvalid, 1, kDontValidate, &gp); 730 return 1; 747 return 1; 731 } 748 } 732 749 >> 750 733 //============================================ 751 //===================================================================== 734 //* GetAreaCode ------------------------------ 752 //* GetAreaCode ------------------------------------------------------- 735 753 736 G4int G4TwistTrapParallelSide::GetAreaCode(con << 754 G4int G4TwistTrapParallelSide::GetAreaCode(const G4ThreeVector &xx, 737 G4b << 755 G4bool withTol) 738 { 756 { 739 // We must use the function in local coordi 757 // We must use the function in local coordinate system. 740 // See the description of DistanceToSurface 758 // See the description of DistanceToSurface(p,v). 741 759 742 const G4double ctol = 0.5 * kCarTolerance; << 760 static const G4double ctol = 0.5 * kCarTolerance; 743 761 744 G4double phi ; 762 G4double phi ; 745 G4double yprime ; 763 G4double yprime ; 746 GetPhiUAtX(xx, phi,yprime ) ; 764 GetPhiUAtX(xx, phi,yprime ) ; 747 765 748 G4double fXAxisMax = GetBoundaryMax(phi) ; 766 G4double fXAxisMax = GetBoundaryMax(phi) ; 749 G4double fXAxisMin = GetBoundaryMin(phi) ; 767 G4double fXAxisMin = GetBoundaryMin(phi) ; 750 768 751 #ifdef G4TWISTDEBUG 769 #ifdef G4TWISTDEBUG 752 G4cout << "GetAreaCode: phi = " << phi << G 770 G4cout << "GetAreaCode: phi = " << phi << G4endl ; 753 G4cout << "GetAreaCode: yprime = " << yprim 771 G4cout << "GetAreaCode: yprime = " << yprime << G4endl ; 754 G4cout << "Intervall is " << fXAxisMin << " 772 G4cout << "Intervall is " << fXAxisMin << " to " << fXAxisMax << G4endl ; 755 #endif 773 #endif 756 774 757 G4int areacode = sInside; 775 G4int areacode = sInside; 758 776 759 if (fAxis[0] == kXAxis && fAxis[1] == kZAxi << 777 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) { 760 { << 778 761 G4int zaxis = 1; 779 G4int zaxis = 1; 762 780 763 if (withTol) << 781 if (withTol) { 764 { << 782 765 G4bool isoutside = false; << 783 G4bool isoutside = false; 766 784 767 // test boundary of xaxis 785 // test boundary of xaxis 768 786 769 if (yprime < fXAxisMin + ctol) << 787 if (yprime < fXAxisMin + ctol) { 770 { << 771 areacode |= (sAxis0 & (sAxisX | sA 788 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary; 772 if (yprime <= fXAxisMin - ctol) is 789 if (yprime <= fXAxisMin - ctol) isoutside = true; 773 790 774 } << 791 } else if (yprime > fXAxisMax - ctol) { 775 else if (yprime > fXAxisMax - ctol) << 776 { << 777 areacode |= (sAxis0 & (sAxisX | sA 792 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary; 778 if (yprime >= fXAxisMax + ctol) i 793 if (yprime >= fXAxisMax + ctol) isoutside = true; 779 } 794 } 780 795 781 // test boundary of z-axis 796 // test boundary of z-axis 782 797 783 if (xx.z() < fAxisMin[zaxis] + ctol) << 798 if (xx.z() < fAxisMin[zaxis] + ctol) { 784 { << 785 areacode |= (sAxis1 & (sAxisZ | sA 799 areacode |= (sAxis1 & (sAxisZ | sAxisMin)); 786 800 787 if ((areacode & sBoundary) != 0) << 801 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner. 788 else areaco 802 else areacode |= sBoundary; 789 if (xx.z() <= fAxisMin[zaxis] - ct 803 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true; 790 804 791 } << 805 } else if (xx.z() > fAxisMax[zaxis] - ctol) { 792 else if (xx.z() > fAxisMax[zaxis] - c << 793 { << 794 areacode |= (sAxis1 & (sAxisZ | sA 806 areacode |= (sAxis1 & (sAxisZ | sAxisMax)); 795 807 796 if ((areacode & sBoundary) != 0) << 808 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner. 797 else areaco 809 else areacode |= sBoundary; 798 if (xx.z() >= fAxisMax[zaxis] + ct 810 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true; 799 } 811 } 800 812 801 // if isoutside = true, clear inside 813 // if isoutside = true, clear inside bit. 802 // if not on boundary, add axis infor 814 // if not on boundary, add axis information. 803 815 804 if (isoutside) << 816 if (isoutside) { 805 { << 806 G4int tmpareacode = areacode & (~s 817 G4int tmpareacode = areacode & (~sInside); 807 areacode = tmpareacode; 818 areacode = tmpareacode; 808 } << 819 } else if ((areacode & sBoundary) != sBoundary) { 809 else if ((areacode & sBoundary) != sB << 810 { << 811 areacode |= (sAxis0 & sAxisX) | (s 820 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ); 812 } 821 } 813 822 814 } << 823 } else { 815 else << 824 816 { << 817 // boundary of y-axis 825 // boundary of y-axis 818 826 819 if (yprime < fXAxisMin ) << 827 if (yprime < fXAxisMin ) { 820 { << 821 areacode |= (sAxis0 & (sAxisX | sA 828 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary; 822 } << 829 } else if (yprime > fXAxisMax) { 823 else if (yprime > fXAxisMax) << 824 { << 825 areacode |= (sAxis0 & (sAxisX | sA 830 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary; 826 } 831 } 827 832 828 // boundary of z-axis 833 // boundary of z-axis 829 834 830 if (xx.z() < fAxisMin[zaxis]) << 835 if (xx.z() < fAxisMin[zaxis]) { 831 { << 832 areacode |= (sAxis1 & (sAxisZ | sA 836 areacode |= (sAxis1 & (sAxisZ | sAxisMin)); 833 if ((areacode & sBoundary) != 0) << 837 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner. 834 else areaco 838 else areacode |= sBoundary; 835 839 836 } << 840 } else if (xx.z() > fAxisMax[zaxis]) { 837 else if (xx.z() > fAxisMax[zaxis]) << 838 { << 839 areacode |= (sAxis1 & (sAxisZ | sA 841 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ; 840 if ((areacode & sBoundary) != 0) << 842 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner. 841 else areaco 843 else areacode |= sBoundary; 842 } 844 } 843 845 844 if ((areacode & sBoundary) != sBounda << 846 if ((areacode & sBoundary) != sBoundary) { 845 { << 846 areacode |= (sAxis0 & sAxisX) | (s 847 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ); 847 } 848 } 848 } 849 } 849 return areacode; 850 return areacode; 850 } << 851 } else { 851 else << 852 { << 853 G4Exception("G4TwistTrapParallelSide::Ge 852 G4Exception("G4TwistTrapParallelSide::GetAreaCode()", 854 "GeomSolids0001", FatalExcep 853 "GeomSolids0001", FatalException, 855 "Feature NOT implemented !") 854 "Feature NOT implemented !"); 856 } 855 } 857 return areacode; 856 return areacode; 858 } 857 } 859 858 860 //============================================ 859 //===================================================================== 861 //* SetCorners() ----------------------------- 860 //* SetCorners() ------------------------------------------------------ 862 861 863 void G4TwistTrapParallelSide::SetCorners() 862 void G4TwistTrapParallelSide::SetCorners() 864 { 863 { 865 864 866 // Set Corner points in local coodinate. 865 // Set Corner points in local coodinate. 867 866 868 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis << 867 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) { 869 { << 868 870 G4double x, y, z; 869 G4double x, y, z; 871 870 872 // corner of Axis0min and Axis1min 871 // corner of Axis0min and Axis1min 873 872 874 x = -fdeltaX/2. + (-fDx2 + fDy1*fTAlph)*st << 873 x = -fdeltaX/2. + (-fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ; 875 + fDy1*std::sin(fPhiTwist/2.) ; << 874 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) + (fDx2 - fDy1*fTAlph)*std::sin(fPhiTwist/2.) ; 876 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/ << 877 + (fDx2 - fDy1*fTAlph)*std::sin(fPhiTwis << 878 z = -fDz ; 875 z = -fDz ; 879 876 880 SetCorner(sC0Min1Min, x, y, z); 877 SetCorner(sC0Min1Min, x, y, z); 881 878 882 // corner of Axis0max and Axis1min 879 // corner of Axis0max and Axis1min 883 880 884 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std << 881 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ; 885 + fDy1*std::sin(fPhiTwist/2.) ; << 882 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ; 886 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/ << 887 - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwis << 888 z = -fDz; 883 z = -fDz; 889 884 890 SetCorner(sC0Max1Min, x, y, z); 885 SetCorner(sC0Max1Min, x, y, z); 891 886 892 // corner of Axis0max and Axis1max 887 // corner of Axis0max and Axis1max 893 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std: << 888 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ; 894 - fDy2*std::sin(fPhiTwist/2.) ; << 889 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ; 895 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2 << 896 + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwis << 897 z = fDz ; 890 z = fDz ; 898 891 899 SetCorner(sC0Max1Max, x, y, z); 892 SetCorner(sC0Max1Max, x, y, z); 900 893 901 // corner of Axis0min and Axis1max 894 // corner of Axis0min and Axis1max 902 x = fdeltaX/2. + (-fDx4 + fDy2*fTAlph)*std << 895 x = fdeltaX/2. + (-fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ; 903 - fDy2*std::sin(fPhiTwist/2.) ; << 896 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (-fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ; 904 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2 << 905 + (-fDx4 + fDy2*fTAlph)*std::sin(fPhiTwi << 906 z = fDz ; 897 z = fDz ; 907 898 908 SetCorner(sC0Min1Max, x, y, z); 899 SetCorner(sC0Min1Max, x, y, z); 909 } << 900 910 else << 901 } else { 911 { << 902 912 G4Exception("G4TwistTrapParallelSide::SetC 903 G4Exception("G4TwistTrapParallelSide::SetCorners()", 913 "GeomSolids0001", FatalExcepti 904 "GeomSolids0001", FatalException, 914 "Method NOT implemented !"); 905 "Method NOT implemented !"); 915 } 906 } 916 } 907 } 917 908 918 //============================================ 909 //===================================================================== 919 //* SetBoundaries() -------------------------- 910 //* SetBoundaries() --------------------------------------------------- 920 911 921 void G4TwistTrapParallelSide::SetBoundaries() 912 void G4TwistTrapParallelSide::SetBoundaries() 922 { 913 { 923 // Set direction-unit vector of boundary-li 914 // Set direction-unit vector of boundary-lines in local coodinate. 924 // 915 // 925 916 926 G4ThreeVector direction; 917 G4ThreeVector direction; 927 918 928 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis << 919 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) { 929 { << 920 930 // sAxis0 & sAxisMin 921 // sAxis0 & sAxisMin 931 direction = GetCorner(sC0Min1Max) - GetCor 922 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min); 932 direction = direction.unit(); 923 direction = direction.unit(); 933 SetBoundary(sAxis0 & (sAxisX | sAxisMin), 924 SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction, 934 GetCorner(sC0Min1Min), sAxisZ) 925 GetCorner(sC0Min1Min), sAxisZ) ; 935 926 936 // sAxis0 & sAxisMax 927 // sAxis0 & sAxisMax 937 direction = GetCorner(sC0Max1Max) - GetCor 928 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min); 938 direction = direction.unit(); 929 direction = direction.unit(); 939 SetBoundary(sAxis0 & (sAxisX | sAxisMax), 930 SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction, 940 GetCorner(sC0Max1Min), sAxisZ) 931 GetCorner(sC0Max1Min), sAxisZ); 941 932 942 // sAxis1 & sAxisMin 933 // sAxis1 & sAxisMin 943 direction = GetCorner(sC0Max1Min) - GetCor 934 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min); 944 direction = direction.unit(); 935 direction = direction.unit(); 945 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), 936 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction, 946 GetCorner(sC0Min1Min), sAxisX) 937 GetCorner(sC0Min1Min), sAxisX); 947 938 948 // sAxis1 & sAxisMax 939 // sAxis1 & sAxisMax 949 direction = GetCorner(sC0Max1Max) - GetCor 940 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max); 950 direction = direction.unit(); 941 direction = direction.unit(); 951 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), 942 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction, 952 GetCorner(sC0Min1Max), sAxisX) 943 GetCorner(sC0Min1Max), sAxisX); >> 944 >> 945 } else { >> 946 >> 947 G4Exception("G4TwistTrapParallelSide::SetCorners()", >> 948 "GeomSolids0001", FatalException, >> 949 "Feature NOT implemented !"); 953 } 950 } 954 else << 951 955 { << 956 G4Exception("G4TwistTrapParallelSide::SetC << 957 "GeomSolids0001", FatalExcepti << 958 "Feature NOT implemented !"); << 959 } << 960 } 952 } 961 953 962 //============================================ 954 //===================================================================== 963 //* GetPhiUAtX() ----------------------------- 955 //* GetPhiUAtX() ------------------------------------------------------ 964 956 965 void 957 void 966 G4TwistTrapParallelSide::GetPhiUAtX( const G4T << 958 G4TwistTrapParallelSide::GetPhiUAtX( G4ThreeVector p, G4double &phi, G4double &u) 967 G4double& << 968 { 959 { 969 // find closest point XX on surface for a gi 960 // find closest point XX on surface for a given point p 970 // X0 is a point on the surface, d is the d << 961 // X0 is a point on the surface, d is the direction ( both for a fixed z = pz) 971 // ( both for a fixed z = pz) << 972 962 973 // phi is given by the z coordinate of p 963 // phi is given by the z coordinate of p 974 964 975 phi = p.z()/(2*fDz)*fPhiTwist ; 965 phi = p.z()/(2*fDz)*fPhiTwist ; 976 966 977 u = ((-(fdeltaX*phi) + fPhiTwist*p.x())* std << 967 u = ((-(fdeltaX*phi) + fPhiTwist*p.x())* std::cos(phi) + (-(fdeltaY*phi) + fPhiTwist*p.y())*std::sin(phi))/fPhiTwist ; 978 + (-(fdeltaY*phi) + fPhiTwist*p.y())*std:: << 968 979 } 969 } 980 970 981 //============================================ 971 //===================================================================== 982 //* ProjectPoint() --------------------------- 972 //* ProjectPoint() ---------------------------------------------------- 983 973 984 G4ThreeVector G4TwistTrapParallelSide::Project << 974 G4ThreeVector G4TwistTrapParallelSide::ProjectPoint(const G4ThreeVector &p, 985 975 G4bool isglobal) 986 { 976 { 987 // Get Rho at p.z() on Hyperbolic Surface. 977 // Get Rho at p.z() on Hyperbolic Surface. 988 G4ThreeVector tmpp; 978 G4ThreeVector tmpp; 989 if (isglobal) << 979 if (isglobal) { 990 { << 991 tmpp = fRot.inverse()*p - fTrans; 980 tmpp = fRot.inverse()*p - fTrans; 992 } << 981 } else { 993 else << 994 { << 995 tmpp = p; 982 tmpp = p; 996 } 983 } 997 984 998 G4double phi ; 985 G4double phi ; 999 G4double u ; 986 G4double u ; 1000 987 1001 GetPhiUAtX( tmpp, phi, u ) ; // calculate << 988 GetPhiUAtX( tmpp, phi, u ) ; // calculate (phi, u) for a point p close the surface 1002 989 1003 G4ThreeVector xx = SurfacePoint(phi,u) ; // << 990 G4ThreeVector xx = SurfacePoint(phi,u) ; // transform back to cartesian coordinates 1004 991 1005 if (isglobal) << 992 if (isglobal) { 1006 { << 1007 return (fRot * xx + fTrans); 993 return (fRot * xx + fTrans); 1008 } << 994 } else { 1009 else << 1010 { << 1011 return xx; 995 return xx; 1012 } 996 } 1013 } 997 } 1014 998 1015 //=========================================== 999 //===================================================================== 1016 //* GetFacets() ----------------------------- 1000 //* GetFacets() ------------------------------------------------------- 1017 1001 1018 void G4TwistTrapParallelSide::GetFacets( G4in 1002 void G4TwistTrapParallelSide::GetFacets( G4int k, G4int n, G4double xyz[][3], 1019 G4in 1003 G4int faces[][4], G4int iside ) 1020 { 1004 { >> 1005 1021 G4double phi ; 1006 G4double phi ; 1022 G4double z, u ; // the two parameters f 1007 G4double z, u ; // the two parameters for the surface equation 1023 G4ThreeVector p ; // a point on the surfac 1008 G4ThreeVector p ; // a point on the surface, given by (z,u) 1024 1009 1025 G4int nnode ; 1010 G4int nnode ; 1026 G4int nface ; 1011 G4int nface ; 1027 1012 1028 G4double umin, umax ; 1013 G4double umin, umax ; 1029 1014 1030 // calculate the (n-1)*(k-1) vertices 1015 // calculate the (n-1)*(k-1) vertices 1031 1016 1032 for ( G4int i = 0 ; i<n ; ++i ) << 1017 G4int i,j ; 1033 { << 1018 >> 1019 for ( i = 0 ; i<n ; i++ ) { >> 1020 1034 z = -fDz+i*(2.*fDz)/(n-1) ; 1021 z = -fDz+i*(2.*fDz)/(n-1) ; 1035 phi = z*fPhiTwist/(2*fDz) ; 1022 phi = z*fPhiTwist/(2*fDz) ; 1036 umin = GetBoundaryMin(phi) ; 1023 umin = GetBoundaryMin(phi) ; 1037 umax = GetBoundaryMax(phi) ; 1024 umax = GetBoundaryMax(phi) ; 1038 1025 1039 for ( G4int j = 0 ; j<k ; ++j ) << 1026 for ( j = 0 ; j<k ; j++ ) { 1040 { << 1027 1041 nnode = GetNode(i,j,k,n,iside) ; 1028 nnode = GetNode(i,j,k,n,iside) ; 1042 u = umax - j*(umax-umin)/(k-1) ; 1029 u = umax - j*(umax-umin)/(k-1) ; 1043 p = SurfacePoint(phi,u,true) ; // surf << 1030 p = SurfacePoint(phi,u,true) ; // surface point in global coordinate system 1044 1031 1045 xyz[nnode][0] = p.x() ; 1032 xyz[nnode][0] = p.x() ; 1046 xyz[nnode][1] = p.y() ; 1033 xyz[nnode][1] = p.y() ; 1047 xyz[nnode][2] = p.z() ; 1034 xyz[nnode][2] = p.z() ; 1048 1035 1049 if ( i<n-1 && j<k-1 ) // conterclock << 1036 if ( i<n-1 && j<k-1 ) { // conterclock wise filling 1050 { << 1037 1051 nface = GetFace(i,j,k,n,iside) ; 1038 nface = GetFace(i,j,k,n,iside) ; 1052 faces[nface][0] = GetEdgeVisibility(i << 1039 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * (GetNode(i ,j ,k,n,iside)+1) ; // fortran numbering 1053 * (GetNode(i ,j ,k, << 1040 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * (GetNode(i ,j+1,k,n,iside)+1) ; 1054 faces[nface][1] = GetEdgeVisibility(i << 1041 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * (GetNode(i+1,j+1,k,n,iside)+1) ; 1055 * (GetNode(i ,j+1,k, << 1042 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * (GetNode(i+1,j ,k,n,iside)+1) ; 1056 faces[nface][2] = GetEdgeVisibility(i << 1043 1057 * (GetNode(i+1,j+1,k, << 1058 faces[nface][3] = GetEdgeVisibility(i << 1059 * (GetNode(i+1,j ,k, << 1060 } 1044 } 1061 } 1045 } 1062 } 1046 } 1063 } 1047 } 1064 1048