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