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