| 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 // G4TwistTubsSide implementation 26 // G4TwistTubsSide implementation
27 // 27 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu 28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp), created.
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), 29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
30 // from original version in Jupi 30 // from original version in Jupiter-2.5.02 application.
31 // ------------------------------------------- 31 // --------------------------------------------------------------------
32 32
33 #include "G4TwistTubsSide.hh" 33 #include "G4TwistTubsSide.hh"
34 34
35 //============================================ 35 //=====================================================================
36 //* constructors ----------------------------- 36 //* constructors ------------------------------------------------------
37 37
38 G4TwistTubsSide::G4TwistTubsSide(const G4Strin 38 G4TwistTubsSide::G4TwistTubsSide(const G4String& name,
39 const G4Rotat 39 const G4RotationMatrix& rot,
40 const G4Three 40 const G4ThreeVector& tlate,
41 G4int 41 G4int handedness,
42 const G4doubl 42 const G4double kappa,
43 const EAxis 43 const EAxis axis0,
44 const EAxis 44 const EAxis axis1,
45 G4doubl 45 G4double axis0min,
46 G4doubl 46 G4double axis1min,
47 G4doubl 47 G4double axis0max,
48 G4doubl 48 G4double axis1max)
49 : G4VTwistSurface(name, rot, tlate, handedn 49 : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
50 axis0min, axis1min, axis0 50 axis0min, axis1min, axis0max, axis1max),
51 fKappa(kappa) 51 fKappa(kappa)
52 { 52 {
53 if (axis0 == kZAxis && axis1 == kXAxis) 53 if (axis0 == kZAxis && axis1 == kXAxis)
54 { 54 {
55 G4Exception("G4TwistTubsSide::G4TwistTub 55 G4Exception("G4TwistTubsSide::G4TwistTubsSide()", "GeomSolids0002",
56 FatalErrorInArgument, "Shoul 56 FatalErrorInArgument, "Should swap axis0 and axis1!");
57 } 57 }
58 fIsValidNorm = false; 58 fIsValidNorm = false;
59 SetCorners(); 59 SetCorners();
60 SetBoundaries(); 60 SetBoundaries();
61 } 61 }
62 62
63 G4TwistTubsSide::G4TwistTubsSide(const G4Strin 63 G4TwistTubsSide::G4TwistTubsSide(const G4String& name,
64 G4doubl 64 G4double EndInnerRadius[2],
65 G4doubl 65 G4double EndOuterRadius[2],
66 G4doubl 66 G4double DPhi,
67 G4doubl 67 G4double EndPhi[2],
68 G4doubl 68 G4double EndZ[2],
69 G4doubl 69 G4double InnerRadius,
70 G4doubl 70 G4double OuterRadius,
71 G4doubl 71 G4double Kappa,
72 G4int 72 G4int handedness)
73 : G4VTwistSurface(name) 73 : G4VTwistSurface(name)
74 { 74 {
75 fHandedness = handedness; // +z = +ve, -z 75 fHandedness = handedness; // +z = +ve, -z = -ve
76 fAxis[0] = kXAxis; // in local coordinat 76 fAxis[0] = kXAxis; // in local coordinate system
77 fAxis[1] = kZAxis; 77 fAxis[1] = kZAxis;
78 fAxisMin[0] = InnerRadius; // Inner-hype r 78 fAxisMin[0] = InnerRadius; // Inner-hype radius at z=0
79 fAxisMax[0] = OuterRadius; // Outer-hype r 79 fAxisMax[0] = OuterRadius; // Outer-hype radius at z=0
80 fAxisMin[1] = EndZ[0]; 80 fAxisMin[1] = EndZ[0];
81 fAxisMax[1] = EndZ[1]; 81 fAxisMax[1] = EndZ[1];
82 82
83 fKappa = Kappa; 83 fKappa = Kappa;
84 fRot.rotateZ( fHandedness > 0 84 fRot.rotateZ( fHandedness > 0
85 ? -0.5*DPhi 85 ? -0.5*DPhi
86 : 0.5*DPhi ); 86 : 0.5*DPhi );
87 fTrans.set(0, 0, 0); 87 fTrans.set(0, 0, 0);
88 fIsValidNorm = false; 88 fIsValidNorm = false;
89 89
90 SetCorners( EndInnerRadius, EndOuterRadius, 90 SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
91 SetBoundaries(); 91 SetBoundaries();
92 } 92 }
93 93
94 //============================================ 94 //=====================================================================
95 //* Fake default constructor ----------------- 95 //* Fake default constructor ------------------------------------------
96 96
97 G4TwistTubsSide::G4TwistTubsSide( __void__& a 97 G4TwistTubsSide::G4TwistTubsSide( __void__& a )
98 : G4VTwistSurface(a) << 98 : G4VTwistSurface(a), fKappa(0.)
99 { 99 {
100 } 100 }
101 101
102 102
103 //============================================ 103 //=====================================================================
104 //* destructor ------------------------------- 104 //* destructor --------------------------------------------------------
105 105
106 G4TwistTubsSide::~G4TwistTubsSide() = default; << 106 G4TwistTubsSide::~G4TwistTubsSide()
>> 107 {
>> 108 }
107 109
108 //============================================ 110 //=====================================================================
109 //* GetNormal -------------------------------- 111 //* GetNormal ---------------------------------------------------------
110 112
111 G4ThreeVector G4TwistTubsSide::GetNormal(const 113 G4ThreeVector G4TwistTubsSide::GetNormal(const G4ThreeVector& tmpxx,
112 114 G4bool isGlobal)
113 { 115 {
114 // GetNormal returns a normal vector at a s 116 // GetNormal returns a normal vector at a surface (or very close
115 // to surface) point at tmpxx. 117 // to surface) point at tmpxx.
116 // If isGlobal=true, it returns the normal 118 // If isGlobal=true, it returns the normal in global coordinate.
117 // 119 //
118 G4ThreeVector xx; 120 G4ThreeVector xx;
119 if (isGlobal) 121 if (isGlobal)
120 { 122 {
121 xx = ComputeLocalPoint(tmpxx); 123 xx = ComputeLocalPoint(tmpxx);
122 if ((xx - fCurrentNormal.p).mag() < 0.5 124 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
123 { 125 {
124 return ComputeGlobalDirection(fCurren 126 return ComputeGlobalDirection(fCurrentNormal.normal);
125 } 127 }
126 } 128 }
127 else 129 else
128 { 130 {
129 xx = tmpxx; 131 xx = tmpxx;
130 if (xx == fCurrentNormal.p) 132 if (xx == fCurrentNormal.p)
131 { 133 {
132 return fCurrentNormal.normal; 134 return fCurrentNormal.normal;
133 } 135 }
134 } 136 }
135 137
136 G4ThreeVector er(1, fKappa * xx.z(), 0); 138 G4ThreeVector er(1, fKappa * xx.z(), 0);
137 G4ThreeVector ez(0, fKappa * xx.x(), 1); 139 G4ThreeVector ez(0, fKappa * xx.x(), 1);
138 G4ThreeVector normal = fHandedness*(er.cros 140 G4ThreeVector normal = fHandedness*(er.cross(ez));
139 141
140 if (isGlobal) 142 if (isGlobal)
141 { 143 {
142 fCurrentNormal.normal = ComputeGlobalDir 144 fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
143 } 145 }
144 else 146 else
145 { 147 {
146 fCurrentNormal.normal = normal.unit(); 148 fCurrentNormal.normal = normal.unit();
147 } 149 }
148 return fCurrentNormal.normal; 150 return fCurrentNormal.normal;
149 } 151 }
150 152
151 //============================================ 153 //=====================================================================
152 //* DistanceToSurface ------------------------ 154 //* DistanceToSurface -------------------------------------------------
153 155
154 G4int G4TwistTubsSide::DistanceToSurface(const 156 G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector& gp,
155 const 157 const G4ThreeVector& gv,
156 158 G4ThreeVector gxx[],
157 159 G4double distance[],
158 160 G4int areacode[],
159 161 G4bool isvalid[],
160 162 EValidate validate)
161 { 163 {
162 // Coordinate system: 164 // Coordinate system:
163 // 165 //
164 // The coordinate system is so chosen th 166 // The coordinate system is so chosen that the intersection of
165 // the twisted surface with the z=0 plan 167 // the twisted surface with the z=0 plane coincides with the
166 // x-axis. 168 // x-axis.
167 // Rotation matrix from this coordinate 169 // Rotation matrix from this coordinate system (local system)
168 // to global system is saved in fRot fie 170 // to global system is saved in fRot field.
169 // So the (global) particle position and 171 // So the (global) particle position and (global) velocity vectors,
170 // p and v, should be rotated fRot.inver 172 // p and v, should be rotated fRot.inverse() in order to convert
171 // to local vectors. 173 // to local vectors.
172 // 174 //
173 // Equation of a twisted surface: 175 // Equation of a twisted surface:
174 // 176 //
175 // x(rho(z=0), z) = rho(z=0) 177 // x(rho(z=0), z) = rho(z=0)
176 // y(rho(z=0), z) = rho(z=0)*K*z 178 // y(rho(z=0), z) = rho(z=0)*K*z
177 // z(rho(z=0), z) = z 179 // z(rho(z=0), z) = z
178 // with 180 // with
179 // K = std::tan(fPhiTwist/2)/fZHalfLe 181 // K = std::tan(fPhiTwist/2)/fZHalfLen
180 // 182 //
181 // Equation of a line: 183 // Equation of a line:
182 // 184 //
183 // gxx = p + t*v 185 // gxx = p + t*v
184 // with 186 // with
185 // p = fRot.inverse()*gp 187 // p = fRot.inverse()*gp
186 // v = fRot.inverse()*gv 188 // v = fRot.inverse()*gv
187 // 189 //
188 // Solution for intersection: 190 // Solution for intersection:
189 // 191 //
190 // Required time for crossing is given b 192 // Required time for crossing is given by solving the
191 // following quadratic equation: 193 // following quadratic equation:
192 // 194 //
193 // a*t^2 + b*t + c = 0 195 // a*t^2 + b*t + c = 0
194 // 196 //
195 // where 197 // where
196 // 198 //
197 // a = K*v_x*v_z 199 // a = K*v_x*v_z
198 // b = K*(v_x*p_z + v_z*p_x) - v_y 200 // b = K*(v_x*p_z + v_z*p_x) - v_y
199 // c = K*p_x*p_z - p_y 201 // c = K*p_x*p_z - p_y
200 // 202 //
201 // Out of the possible two solutions you 203 // Out of the possible two solutions you must choose
202 // the one that gives a positive rho(z=0 204 // the one that gives a positive rho(z=0).
203 // 205 //
204 // 206 //
205 207
206 fCurStatWithV.ResetfDone(validate, &gp, &gv 208 fCurStatWithV.ResetfDone(validate, &gp, &gv);
207 209
208 if (fCurStatWithV.IsDone()) 210 if (fCurStatWithV.IsDone())
209 { 211 {
210 for (G4int i=0; i<fCurStatWithV.GetNXX() 212 for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
211 { 213 {
212 gxx[i] = fCurStatWithV.GetXX(i); 214 gxx[i] = fCurStatWithV.GetXX(i);
213 distance[i] = fCurStatWithV.GetDistan 215 distance[i] = fCurStatWithV.GetDistance(i);
214 areacode[i] = fCurStatWithV.GetAreaco 216 areacode[i] = fCurStatWithV.GetAreacode(i);
215 isvalid[i] = fCurStatWithV.IsValid(i 217 isvalid[i] = fCurStatWithV.IsValid(i);
216 } 218 }
217 return fCurStatWithV.GetNXX(); 219 return fCurStatWithV.GetNXX();
218 } 220 }
219 else // initialize 221 else // initialize
220 { 222 {
221 for (auto i=0; i<2; ++i) 223 for (auto i=0; i<2; ++i)
222 { 224 {
223 distance[i] = kInfinity; 225 distance[i] = kInfinity;
224 areacode[i] = sOutside; 226 areacode[i] = sOutside;
225 isvalid[i] = false; 227 isvalid[i] = false;
226 gxx[i].set(kInfinity, kInfinity, kInf 228 gxx[i].set(kInfinity, kInfinity, kInfinity);
227 } 229 }
228 } 230 }
229 231
230 G4ThreeVector p = ComputeLocalPoint(gp); 232 G4ThreeVector p = ComputeLocalPoint(gp);
231 G4ThreeVector v = ComputeLocalDirection(gv) 233 G4ThreeVector v = ComputeLocalDirection(gv);
232 G4ThreeVector xx[2]; 234 G4ThreeVector xx[2];
233 235
234 // 236 //
235 // special case! 237 // special case!
236 // p is origin or 238 // p is origin or
237 // 239 //
238 240
239 G4double absvz = std::fabs(v.z()); 241 G4double absvz = std::fabs(v.z());
240 242
241 if ((absvz<DBL_MIN) && (std::fabs(p.x() * v 243 if ((absvz<DBL_MIN) && (std::fabs(p.x() * v.y() - p.y() * v.x())<DBL_MIN))
242 { 244 {
243 // no intersection 245 // no intersection
244 246
245 isvalid[0] = false; 247 isvalid[0] = false;
246 fCurStat.SetCurrentStatus(0, gxx[0], dis 248 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
247 isvalid[0], 0, 249 isvalid[0], 0, validate, &gp, &gv);
248 return 0; 250 return 0;
249 } 251 }
250 252
251 // 253 //
252 // special case end 254 // special case end
253 // 255 //
254 256
255 G4double a = fKappa * v.x() * v.z(); 257 G4double a = fKappa * v.x() * v.z();
256 G4double b = fKappa * (v.x() * p.z() + v.z( 258 G4double b = fKappa * (v.x() * p.z() + v.z() * p.x()) - v.y();
257 G4double c = fKappa * p.x() * p.z() - p.y() 259 G4double c = fKappa * p.x() * p.z() - p.y();
258 G4double D = b * b - 4 * a * c; 260 G4double D = b * b - 4 * a * c; // discriminant
259 G4int vout = 0; 261 G4int vout = 0;
260 262
261 if (std::fabs(a) < DBL_MIN) 263 if (std::fabs(a) < DBL_MIN)
262 { 264 {
263 if (std::fabs(b) > DBL_MIN) 265 if (std::fabs(b) > DBL_MIN)
264 { 266 {
265 // single solution 267 // single solution
266 268
267 distance[0] = - c / b; 269 distance[0] = - c / b;
268 xx[0] = p + distance[0]*v; 270 xx[0] = p + distance[0]*v;
269 gxx[0] = ComputeGlobalPoint(xx[0 271 gxx[0] = ComputeGlobalPoint(xx[0]);
270 272
271 if (validate == kValidateWithTol) 273 if (validate == kValidateWithTol)
272 { 274 {
273 areacode[0] = GetAreaCode(xx[0]); 275 areacode[0] = GetAreaCode(xx[0]);
274 if (!IsOutside(areacode[0])) 276 if (!IsOutside(areacode[0]))
275 { 277 {
276 if (distance[0] >= 0) isvalid[0 278 if (distance[0] >= 0) isvalid[0] = true;
277 } 279 }
278 } 280 }
279 else if (validate == kValidateWithout 281 else if (validate == kValidateWithoutTol)
280 { 282 {
281 areacode[0] = GetAreaCode(xx[0], f 283 areacode[0] = GetAreaCode(xx[0], false);
282 if (IsInside(areacode[0])) 284 if (IsInside(areacode[0]))
283 { 285 {
284 if (distance[0] >= 0) isvalid[0 286 if (distance[0] >= 0) isvalid[0] = true;
285 } 287 }
286 } 288 }
287 else // kDontValidate 289 else // kDontValidate
288 { 290 {
289 // we must omit x(rho,z) = rho(z=0 291 // we must omit x(rho,z) = rho(z=0) < 0
290 if (xx[0].x() > 0) 292 if (xx[0].x() > 0)
291 { 293 {
292 areacode[0] = sInside; 294 areacode[0] = sInside;
293 if (distance[0] >= 0) isvalid[0 295 if (distance[0] >= 0) isvalid[0] = true;
294 } 296 }
295 else 297 else
296 { 298 {
297 distance[0] = kInfinity; 299 distance[0] = kInfinity;
298 fCurStatWithV.SetCurrentStatus( 300 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
299 301 areacode[0], isvalid[0],
300 302 0, validate, &gp, &gv);
301 return vout; 303 return vout;
302 } 304 }
303 } 305 }
304 306
305 fCurStatWithV.SetCurrentStatus(0, gxx 307 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
306 isvali 308 isvalid[0], 1, validate, &gp, &gv);
307 vout = 1; 309 vout = 1;
308 } 310 }
309 else 311 else
310 { 312 {
311 // if a=b=0 , v.y=0 and (v.x=0 && p.x 313 // if a=b=0 , v.y=0 and (v.x=0 && p.x=0) or (v.z=0 && p.z=0) .
312 // if v.x=0 && p.x=0, no intersect 314 // if v.x=0 && p.x=0, no intersection unless p is on z-axis
313 // (in that case, v is paralell to 315 // (in that case, v is paralell to surface).
314 // if v.z=0 && p.z=0, no intersect 316 // if v.z=0 && p.z=0, no intersection unless p is on x-axis
315 // (in that case, v is paralell to 317 // (in that case, v is paralell to surface).
316 // return distance = infinity. 318 // return distance = infinity.
317 319
318 fCurStatWithV.SetCurrentStatus(0, gxx 320 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
319 isvali 321 isvalid[0], 0, validate, &gp, &gv);
320 } 322 }
321 } 323 }
322 else if (D > DBL_MIN) 324 else if (D > DBL_MIN)
323 { 325 {
324 // double solutions 326 // double solutions
325 327
326 D = std::sqrt(D); 328 D = std::sqrt(D);
327 G4double factor = 0.5/a; 329 G4double factor = 0.5/a;
328 G4double tmpdist[2] = {kInfinity, k 330 G4double tmpdist[2] = {kInfinity, kInfinity};
329 G4ThreeVector tmpxx[2]; 331 G4ThreeVector tmpxx[2];
330 G4int tmpareacode[2] = {sOutside 332 G4int tmpareacode[2] = {sOutside, sOutside};
331 G4bool tmpisvalid[2] = {false, f 333 G4bool tmpisvalid[2] = {false, false};
332 334
333 for (auto i=0; i<2; ++i) 335 for (auto i=0; i<2; ++i)
334 { 336 {
335 G4double bminusD = - b - D; 337 G4double bminusD = - b - D;
336 338
337 // protection against round off error 339 // protection against round off error
338 //G4double protection = 1.0e-6; 340 //G4double protection = 1.0e-6;
339 G4double protection = 0; 341 G4double protection = 0;
340 if ( b * D < 0 && std::fabs(bminusD / 342 if ( b * D < 0 && std::fabs(bminusD / D) < protection )
341 { 343 {
342 G4double acovbb = (a*c)/(b*b); 344 G4double acovbb = (a*c)/(b*b);
343 tmpdist[i] = - c/b * ( 1 - acovbb 345 tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
344 } 346 }
345 else 347 else
346 { 348 {
347 tmpdist[i] = factor * bminusD; 349 tmpdist[i] = factor * bminusD;
348 } 350 }
349 351
350 D = -D; 352 D = -D;
351 tmpxx[i] = p + tmpdist[i]*v; 353 tmpxx[i] = p + tmpdist[i]*v;
352 354
353 if (validate == kValidateWithTol) 355 if (validate == kValidateWithTol)
354 { 356 {
355 tmpareacode[i] = GetAreaCode(tmpxx 357 tmpareacode[i] = GetAreaCode(tmpxx[i]);
356 if (!IsOutside(tmpareacode[i])) 358 if (!IsOutside(tmpareacode[i]))
357 { 359 {
358 if (tmpdist[i] >= 0) tmpisvalid 360 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
359 continue; 361 continue;
360 } 362 }
361 } 363 }
362 else if (validate == kValidateWithout 364 else if (validate == kValidateWithoutTol)
363 { 365 {
364 tmpareacode[i] = GetAreaCode(tmpxx 366 tmpareacode[i] = GetAreaCode(tmpxx[i], false);
365 if (IsInside(tmpareacode[i])) 367 if (IsInside(tmpareacode[i]))
366 { 368 {
367 if (tmpdist[i] >= 0) tmpisvalid 369 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
368 continue; 370 continue;
369 } 371 }
370 } 372 }
371 else // kDontValidate 373 else // kDontValidate
372 { 374 {
373 // we must choose x(rho,z) = rho(z 375 // we must choose x(rho,z) = rho(z=0) > 0
374 if (tmpxx[i].x() > 0) 376 if (tmpxx[i].x() > 0)
375 { 377 {
376 tmpareacode[i] = sInside; 378 tmpareacode[i] = sInside;
377 if (tmpdist[i] >= 0) tmpisvalid 379 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
378 continue; 380 continue;
379 } else { 381 } else {
380 tmpdist[i] = kInfinity; 382 tmpdist[i] = kInfinity;
381 continue; 383 continue;
382 } 384 }
383 } 385 }
384 } 386 }
385 387
386 if (tmpdist[0] <= tmpdist[1]) 388 if (tmpdist[0] <= tmpdist[1])
387 { 389 {
388 distance[0] = tmpdist[0]; 390 distance[0] = tmpdist[0];
389 distance[1] = tmpdist[1]; 391 distance[1] = tmpdist[1];
390 xx[0] = tmpxx[0]; 392 xx[0] = tmpxx[0];
391 xx[1] = tmpxx[1]; 393 xx[1] = tmpxx[1];
392 gxx[0] = ComputeGlobalPoint(tmpx 394 gxx[0] = ComputeGlobalPoint(tmpxx[0]);
393 gxx[1] = ComputeGlobalPoint(tmpx 395 gxx[1] = ComputeGlobalPoint(tmpxx[1]);
394 areacode[0] = tmpareacode[0]; 396 areacode[0] = tmpareacode[0];
395 areacode[1] = tmpareacode[1]; 397 areacode[1] = tmpareacode[1];
396 isvalid[0] = tmpisvalid[0]; 398 isvalid[0] = tmpisvalid[0];
397 isvalid[1] = tmpisvalid[1]; 399 isvalid[1] = tmpisvalid[1];
398 } 400 }
399 else 401 else
400 { 402 {
401 distance[0] = tmpdist[1]; 403 distance[0] = tmpdist[1];
402 distance[1] = tmpdist[0]; 404 distance[1] = tmpdist[0];
403 xx[0] = tmpxx[1]; 405 xx[0] = tmpxx[1];
404 xx[1] = tmpxx[0]; 406 xx[1] = tmpxx[0];
405 gxx[0] = ComputeGlobalPoint(tmpx 407 gxx[0] = ComputeGlobalPoint(tmpxx[1]);
406 gxx[1] = ComputeGlobalPoint(tmpx 408 gxx[1] = ComputeGlobalPoint(tmpxx[0]);
407 areacode[0] = tmpareacode[1]; 409 areacode[0] = tmpareacode[1];
408 areacode[1] = tmpareacode[0]; 410 areacode[1] = tmpareacode[0];
409 isvalid[0] = tmpisvalid[1]; 411 isvalid[0] = tmpisvalid[1];
410 isvalid[1] = tmpisvalid[0]; 412 isvalid[1] = tmpisvalid[0];
411 } 413 }
412 414
413 fCurStatWithV.SetCurrentStatus(0, gxx[0] 415 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
414 isvalid[0 416 isvalid[0], 2, validate, &gp, &gv);
415 fCurStatWithV.SetCurrentStatus(1, gxx[1] 417 fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
416 isvalid[1 418 isvalid[1], 2, validate, &gp, &gv);
417 419
418 // protection against roundoff error 420 // protection against roundoff error
419 421
420 for (G4int k=0; k<2; ++k) 422 for (G4int k=0; k<2; ++k)
421 { 423 {
422 if (!isvalid[k]) continue; 424 if (!isvalid[k]) continue;
423 425
424 G4ThreeVector xxonsurface(xx[k].x(), 426 G4ThreeVector xxonsurface(xx[k].x(), fKappa * std::fabs(xx[k].x())
425 427 * xx[k].z() , xx[k].z());
426 G4double deltaY = (xx[k] - xxo 428 G4double deltaY = (xx[k] - xxonsurface).mag();
427 429
428 if ( deltaY > 0.5*kCarTolerance ) 430 if ( deltaY > 0.5*kCarTolerance )
429 { 431 {
430 G4int maxcount = 10; 432 G4int maxcount = 10;
431 G4int l; 433 G4int l;
432 G4double lastdeltaY = deltaY; 434 G4double lastdeltaY = deltaY;
433 for (l=0; l<maxcount; ++l) 435 for (l=0; l<maxcount; ++l)
434 { 436 {
435 G4ThreeVector surfacenormal = Get 437 G4ThreeVector surfacenormal = GetNormal(xxonsurface);
436 distance[k] = DistanceToPlaneWith 438 distance[k] = DistanceToPlaneWithV(p, v, xxonsurface,
437 439 surfacenormal, xx[k]);
438 deltaY = (xx[k] - xxonsurfac 440 deltaY = (xx[k] - xxonsurface).mag();
439 if (deltaY > lastdeltaY) { } // 441 if (deltaY > lastdeltaY) { } // ???
440 gxx[k] = ComputeGlobalPoint( 442 gxx[k] = ComputeGlobalPoint(xx[k]);
441 443
442 if (deltaY <= 0.5*kCarTolerance) 444 if (deltaY <= 0.5*kCarTolerance) break;
443 xxonsurface.set(xx[k].x(), 445 xxonsurface.set(xx[k].x(),
444 fKappa * std::fab 446 fKappa * std::fabs(xx[k].x()) * xx[k].z(),
445 xx[k].z()); 447 xx[k].z());
446 } 448 }
447 if (l == maxcount) 449 if (l == maxcount)
448 { 450 {
449 std::ostringstream message; 451 std::ostringstream message;
450 message << "Exceeded maxloop coun 452 message << "Exceeded maxloop count!" << G4endl
451 << " maxloop count 453 << " maxloop count " << maxcount;
452 G4Exception("G4TwistTubsFlatSide: 454 G4Exception("G4TwistTubsFlatSide::DistanceToSurface(p,v)",
453 "GeomSolids0003", Fa 455 "GeomSolids0003", FatalException, message);
454 } 456 }
455 } 457 }
456 } 458 }
457 vout = 2; 459 vout = 2;
458 } 460 }
459 else 461 else
460 { 462 {
461 // if D<0, no solution 463 // if D<0, no solution
462 // if D=0, just grazing the surfaces, re 464 // if D=0, just grazing the surfaces, return kInfinity
463 465
464 fCurStatWithV.SetCurrentStatus(0, gxx[0] 466 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
465 isvalid[0 467 isvalid[0], 0, validate, &gp, &gv);
466 } 468 }
467 469
468 return vout; 470 return vout;
469 } 471 }
470 472
471 //============================================ 473 //=====================================================================
472 //* DistanceToSurface ------------------------ 474 //* DistanceToSurface -------------------------------------------------
473 475
474 G4int G4TwistTubsSide::DistanceToSurface(const 476 G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector& gp,
475 477 G4ThreeVector gxx[],
476 478 G4double distance[],
477 479 G4int areacode[])
478 { 480 {
479 fCurStat.ResetfDone(kDontValidate, &gp); 481 fCurStat.ResetfDone(kDontValidate, &gp);
480 if (fCurStat.IsDone()) 482 if (fCurStat.IsDone())
481 { 483 {
482 for (G4int i=0; i<fCurStat.GetNXX(); ++i 484 for (G4int i=0; i<fCurStat.GetNXX(); ++i)
483 { 485 {
484 gxx[i] = fCurStat.GetXX(i); 486 gxx[i] = fCurStat.GetXX(i);
485 distance[i] = fCurStat.GetDistance(i) 487 distance[i] = fCurStat.GetDistance(i);
486 areacode[i] = fCurStat.GetAreacode(i) 488 areacode[i] = fCurStat.GetAreacode(i);
487 } 489 }
488 return fCurStat.GetNXX(); 490 return fCurStat.GetNXX();
489 } 491 }
490 else // initialize 492 else // initialize
491 { 493 {
492 for (auto i=0; i<2; ++i) 494 for (auto i=0; i<2; ++i)
493 { 495 {
494 distance[i] = kInfinity; 496 distance[i] = kInfinity;
495 areacode[i] = sOutside; 497 areacode[i] = sOutside;
496 gxx[i].set(kInfinity, kInfinity, kInf 498 gxx[i].set(kInfinity, kInfinity, kInfinity);
497 } 499 }
498 } 500 }
499 501
500 const G4double halftol = 0.5 * kCarToleranc 502 const G4double halftol = 0.5 * kCarTolerance;
501 503
502 G4ThreeVector p = ComputeLocalPoint( 504 G4ThreeVector p = ComputeLocalPoint(gp);
503 G4ThreeVector xx; 505 G4ThreeVector xx;
504 G4int parity = (fKappa >= 0 ? 1 : 506 G4int parity = (fKappa >= 0 ? 1 : -1);
505 507
506 // 508 //
507 // special case! 509 // special case!
508 // If p is on surface, or 510 // If p is on surface, or
509 // p is on z-axis, 511 // p is on z-axis,
510 // return here immediatery. 512 // return here immediatery.
511 // 513 //
512 514
513 G4ThreeVector lastgxx[2]; 515 G4ThreeVector lastgxx[2];
514 for (auto i=0; i<2; ++i) 516 for (auto i=0; i<2; ++i)
515 { 517 {
516 lastgxx[i] = fCurStatWithV.GetXX(i); 518 lastgxx[i] = fCurStatWithV.GetXX(i);
517 } 519 }
518 520
519 if ((gp - lastgxx[0]).mag() < halftol 521 if ((gp - lastgxx[0]).mag() < halftol
520 || (gp - lastgxx[1]).mag() < halftol) 522 || (gp - lastgxx[1]).mag() < halftol)
521 { 523 {
522 // last winner, or last poststep point i 524 // last winner, or last poststep point is on the surface.
523 xx = p; 525 xx = p;
524 distance[0] = 0; 526 distance[0] = 0;
525 gxx[0] = gp; 527 gxx[0] = gp;
526 528
527 G4bool isvalid = true; 529 G4bool isvalid = true;
528 fCurStat.SetCurrentStatus(0, gxx[0], dis 530 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
529 isvalid, 1, kDont 531 isvalid, 1, kDontValidate, &gp);
530 return 1; 532 return 1;
531 } 533 }
532 534
533 if (p.getRho() == 0) 535 if (p.getRho() == 0)
534 { 536 {
535 // p is on z-axis. Namely, p is on twist 537 // p is on z-axis. Namely, p is on twisted surface (invalid area).
536 // We must return here, however, returni 538 // We must return here, however, returning distance to x-minimum
537 // boundary is better than return 0-dist 539 // boundary is better than return 0-distance.
538 // 540 //
539 G4bool isvalid = true; 541 G4bool isvalid = true;
540 if (fAxis[0] == kXAxis && fAxis[1] == kZ 542 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis)
541 { 543 {
542 distance[0] = DistanceToBoundary(sAxi 544 distance[0] = DistanceToBoundary(sAxis0 & sAxisMin, xx, p);
543 areacode[0] = sInside; 545 areacode[0] = sInside;
544 } 546 }
545 else 547 else
546 { 548 {
547 distance[0] = 0; 549 distance[0] = 0;
548 xx.set(0., 0., 0.); 550 xx.set(0., 0., 0.);
549 } 551 }
550 gxx[0] = ComputeGlobalPoint(xx); 552 gxx[0] = ComputeGlobalPoint(xx);
551 fCurStat.SetCurrentStatus(0, gxx[0], dis 553 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
552 isvalid, 0, kD 554 isvalid, 0, kDontValidate, &gp);
553 return 1; 555 return 1;
554 } 556 }
555 557
556 // 558 //
557 // special case end 559 // special case end
558 // 560 //
559 561
560 // set corner points of quadrangle try area 562 // set corner points of quadrangle try area ...
561 563
562 G4ThreeVector A; // foot of normal from p 564 G4ThreeVector A; // foot of normal from p to boundary of sAxis0 & sAxisMin
563 G4ThreeVector C; // foot of normal from p 565 G4ThreeVector C; // foot of normal from p to boundary of sAxis0 & sAxisMax
564 G4ThreeVector B; // point on boundary 566 G4ThreeVector B; // point on boundary sAxis0 & sAxisMax at z = A.z()
565 G4ThreeVector D; // point on boundary 567 G4ThreeVector D; // point on boundary sAxis0 & sAxisMin at z = C.z()
566 568
567 // G4double distToA; // distance from 569 // G4double distToA; // distance from p to A
568 DistanceToBoundary(sAxis0 & sAxisMin, A, p) 570 DistanceToBoundary(sAxis0 & sAxisMin, A, p);
569 // G4double distToC; // distance from 571 // G4double distToC; // distance from p to C
570 DistanceToBoundary(sAxis0 & sAxisMax, C, p) 572 DistanceToBoundary(sAxis0 & sAxisMax, C, p);
571 573
572 // is p.z between a.z and c.z? 574 // is p.z between a.z and c.z?
573 // p.z must be bracketed a.z and c.z. 575 // p.z must be bracketed a.z and c.z.
574 if (A.z() > C.z()) 576 if (A.z() > C.z())
575 { 577 {
576 if (p.z() > A.z()) 578 if (p.z() > A.z())
577 { 579 {
578 A = GetBoundaryAtPZ(sAxis0 & sAxisMin 580 A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
579 } 581 }
580 else if (p.z() < C.z()) 582 else if (p.z() < C.z())
581 { 583 {
582 C = GetBoundaryAtPZ(sAxis0 & sAxisMax 584 C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
583 } 585 }
584 } 586 }
585 else 587 else
586 { 588 {
587 if (p.z() > C.z()) 589 if (p.z() > C.z())
588 { 590 {
589 C = GetBoundaryAtPZ(sAxis0 & sAxisMax 591 C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
590 } 592 }
591 else if (p.z() < A.z()) 593 else if (p.z() < A.z())
592 { 594 {
593 A = GetBoundaryAtPZ(sAxis0 & sAxisMin 595 A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
594 } 596 }
595 } 597 }
596 598
597 G4ThreeVector d[2]; // direction vecto 599 G4ThreeVector d[2]; // direction vectors of boundary
598 G4ThreeVector x0[2]; // foot of normal 600 G4ThreeVector x0[2]; // foot of normal from line to p
599 G4int btype[2]; // boundary type 601 G4int btype[2]; // boundary type
600 602
601 for (auto i=0; i<2; ++i) 603 for (auto i=0; i<2; ++i)
602 { 604 {
603 if (i == 0) 605 if (i == 0)
604 { 606 {
605 GetBoundaryParameters((sAxis0 & sAxis 607 GetBoundaryParameters((sAxis0 & sAxisMax), d[i], x0[i], btype[i]);
606 B = x0[i] + ((A.z() - x0[i].z()) / d[ 608 B = x0[i] + ((A.z() - x0[i].z()) / d[i].z()) * d[i];
607 // x0 + t*d , d is direction unit vec 609 // x0 + t*d , d is direction unit vector.
608 } 610 }
609 else 611 else
610 { 612 {
611 GetBoundaryParameters((sAxis0 & sAxis 613 GetBoundaryParameters((sAxis0 & sAxisMin), d[i], x0[i], btype[i]);
612 D = x0[i] + ((C.z() - x0[i].z()) / d[ 614 D = x0[i] + ((C.z() - x0[i].z()) / d[i].z()) * d[i];
613 } 615 }
614 } 616 }
615 617
616 // In order to set correct diagonal, swap A 618 // In order to set correct diagonal, swap A and D, C and B if needed.
617 G4ThreeVector pt(p.x(), p.y(), 0.); 619 G4ThreeVector pt(p.x(), p.y(), 0.);
618 G4double rc = std::fabs(p.x()); 620 G4double rc = std::fabs(p.x());
619 G4ThreeVector surfacevector(rc, rc * fKappa 621 G4ThreeVector surfacevector(rc, rc * fKappa * p.z(), 0.);
620 G4int pside = AmIOnLeftSide(pt, sur 622 G4int pside = AmIOnLeftSide(pt, surfacevector);
621 G4double test = (A.z() - C.z()) * par 623 G4double test = (A.z() - C.z()) * parity * pside;
622 624
623 if (test == 0) 625 if (test == 0)
624 { 626 {
625 if (pside == 0) 627 if (pside == 0)
626 { 628 {
627 // p is on surface. 629 // p is on surface.
628 xx = p; 630 xx = p;
629 distance[0] = 0; 631 distance[0] = 0;
630 gxx[0] = gp; 632 gxx[0] = gp;
631 633
632 G4bool isvalid = true; 634 G4bool isvalid = true;
633 fCurStat.SetCurrentStatus(0, gxx[0], 635 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
634 isvalid, 1, 636 isvalid, 1, kDontValidate, &gp);
635 return 1; 637 return 1;
636 } 638 }
637 else 639 else
638 { 640 {
639 // A.z = C.z(). return distance to li 641 // A.z = C.z(). return distance to line.
640 d[0] = C - A; 642 d[0] = C - A;
641 distance[0] = DistanceToLine(p, A, d[ 643 distance[0] = DistanceToLine(p, A, d[0], xx);
642 areacode[0] = sInside; 644 areacode[0] = sInside;
643 gxx[0] = ComputeGlobalPoint(xx); 645 gxx[0] = ComputeGlobalPoint(xx);
644 G4bool isvalid = true; 646 G4bool isvalid = true;
645 fCurStat.SetCurrentStatus(0, gxx[0], 647 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
646 isvalid, 1, 648 isvalid, 1, kDontValidate, &gp);
647 return 1; 649 return 1;
648 } 650 }
649 } 651 }
650 else if (test < 0) // wrong diagonal. vect 652 else if (test < 0) // wrong diagonal. vector AC is crossing the surface!
651 { // swap A and D, C and 653 { // swap A and D, C and B
652 G4ThreeVector tmp; 654 G4ThreeVector tmp;
653 tmp = A; 655 tmp = A;
654 A = D; 656 A = D;
655 D = tmp; 657 D = tmp;
656 tmp = C; 658 tmp = C;
657 C = B; 659 C = B;
658 B = tmp; 660 B = tmp;
659 661
660 } 662 }
661 else // correct diagonal. nothing to do. 663 else // correct diagonal. nothing to do.
662 { 664 {
663 } 665 }
664 666
665 // Now, we chose correct diagonal. 667 // Now, we chose correct diagonal.
666 // First try. divide quadrangle into double 668 // First try. divide quadrangle into double triangle by diagonal and
667 // calculate distance to both surfaces. 669 // calculate distance to both surfaces.
668 670
669 G4ThreeVector xxacb; // foot of normal fr 671 G4ThreeVector xxacb; // foot of normal from plane ACB to p
670 G4ThreeVector nacb; // normal of plane A 672 G4ThreeVector nacb; // normal of plane ACD
671 G4ThreeVector xxcad; // foot of normal fr 673 G4ThreeVector xxcad; // foot of normal from plane CAD to p
672 G4ThreeVector ncad; // normal of plane C 674 G4ThreeVector ncad; // normal of plane CAD
673 G4ThreeVector AB(A.x(), A.y(), 0); 675 G4ThreeVector AB(A.x(), A.y(), 0);
674 G4ThreeVector DC(C.x(), C.y(), 0); 676 G4ThreeVector DC(C.x(), C.y(), 0);
675 677
676 G4double distToACB = G4VTwistSurface::Dista 678 G4double distToACB = G4VTwistSurface::DistanceToPlane(p, A, C-A, AB,
677 679 xxacb, nacb) * parity;
678 G4double distToCAD = G4VTwistSurface::Dista 680 G4double distToCAD = G4VTwistSurface::DistanceToPlane(p, C, C-A, DC,
679 681 xxcad, ncad) * parity;
680 // if calculated distance = 0, return 682 // if calculated distance = 0, return
681 683
682 if (std::fabs(distToACB) <= halftol || std: 684 if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol)
683 { 685 {
684 xx = (std::fabs(distToACB) < std::fabs(d 686 xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
685 areacode[0] = sInside; 687 areacode[0] = sInside;
686 gxx[0] = ComputeGlobalPoint(xx); 688 gxx[0] = ComputeGlobalPoint(xx);
687 distance[0] = 0; 689 distance[0] = 0;
688 G4bool isvalid = true; 690 G4bool isvalid = true;
689 fCurStat.SetCurrentStatus(0, gxx[0], dis 691 fCurStat.SetCurrentStatus(0, gxx[0], distance[0] , areacode[0],
690 isvalid, 1, kD 692 isvalid, 1, kDontValidate, &gp);
691 return 1; 693 return 1;
692 } 694 }
693 695
694 if (distToACB * distToCAD > 0 && distToACB 696 if (distToACB * distToCAD > 0 && distToACB < 0)
695 { 697 {
696 // both distToACB and distToCAD are nega 698 // both distToACB and distToCAD are negative.
697 // divide quadrangle into double triangl 699 // divide quadrangle into double triangle by diagonal
698 G4ThreeVector normal; 700 G4ThreeVector normal;
699 distance[0] = DistanceToPlane(p, A, B, C 701 distance[0] = DistanceToPlane(p, A, B, C, D, parity, xx, normal);
700 } 702 }
701 else 703 else
702 { 704 {
703 if (distToACB * distToCAD > 0) 705 if (distToACB * distToCAD > 0)
704 { 706 {
705 // both distToACB and distToCAD are p 707 // both distToACB and distToCAD are positive.
706 // Take smaller one. 708 // Take smaller one.
707 if (distToACB <= distToCAD) 709 if (distToACB <= distToCAD)
708 { 710 {
709 distance[0] = distToACB; 711 distance[0] = distToACB;
710 xx = xxacb; 712 xx = xxacb;
711 } 713 }
712 else 714 else
713 { 715 {
714 distance[0] = distToCAD; 716 distance[0] = distToCAD;
715 xx = xxcad; 717 xx = xxcad;
716 } 718 }
717 } 719 }
718 else 720 else
719 { 721 {
720 // distToACB * distToCAD is negative. 722 // distToACB * distToCAD is negative.
721 // take positive one 723 // take positive one
722 if (distToACB > 0) 724 if (distToACB > 0)
723 { 725 {
724 distance[0] = distToACB; 726 distance[0] = distToACB;
725 xx = xxacb; 727 xx = xxacb;
726 } 728 }
727 else 729 else
728 { 730 {
729 distance[0] = distToCAD; 731 distance[0] = distToCAD;
730 xx = xxcad; 732 xx = xxcad;
731 } 733 }
732 } 734 }
733 } 735 }
734 areacode[0] = sInside; 736 areacode[0] = sInside;
735 gxx[0] = ComputeGlobalPoint(xx); 737 gxx[0] = ComputeGlobalPoint(xx);
736 G4bool isvalid = true; 738 G4bool isvalid = true;
737 fCurStat.SetCurrentStatus(0, gxx[0], distan 739 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
738 isvalid, 1, kDont 740 isvalid, 1, kDontValidate, &gp);
739 return 1; 741 return 1;
740 } 742 }
741 743
742 //============================================ 744 //=====================================================================
743 //* DistanceToPlane -------------------------- 745 //* DistanceToPlane ---------------------------------------------------
744 746
745 G4double G4TwistTubsSide::DistanceToPlane(cons 747 G4double G4TwistTubsSide::DistanceToPlane(const G4ThreeVector& p,
746 cons 748 const G4ThreeVector& A,
747 cons 749 const G4ThreeVector& B,
748 cons 750 const G4ThreeVector& C,
749 cons 751 const G4ThreeVector& D,
750 cons 752 const G4int parity,
751 753 G4ThreeVector& xx,
752 754 G4ThreeVector& n)
753 { 755 {
754 const G4double halftol = 0.5 * kCarToleranc 756 const G4double halftol = 0.5 * kCarTolerance;
755 757
756 G4ThreeVector M = 0.5*(A + B); 758 G4ThreeVector M = 0.5*(A + B);
757 G4ThreeVector N = 0.5*(C + D); 759 G4ThreeVector N = 0.5*(C + D);
758 G4ThreeVector xxanm; // foot of normal fro 760 G4ThreeVector xxanm; // foot of normal from p to plane ANM
759 G4ThreeVector nanm; // normal of plane AN 761 G4ThreeVector nanm; // normal of plane ANM
760 G4ThreeVector xxcmn; // foot of normal fro 762 G4ThreeVector xxcmn; // foot of normal from p to plane CMN
761 G4ThreeVector ncmn; // normal of plane CM 763 G4ThreeVector ncmn; // normal of plane CMN
762 764
763 G4double distToanm = G4VTwistSurface::Dista 765 G4double distToanm = G4VTwistSurface::DistanceToPlane(p, A, (N - A), (M - A),
764 766 xxanm, nanm) * parity;
765 G4double distTocmn = G4VTwistSurface::Dista 767 G4double distTocmn = G4VTwistSurface::DistanceToPlane(p, C, (M - C), (N - C),
766 768 xxcmn, ncmn) * parity;
767 #ifdef G4SPECSDEBUG 769 #ifdef G4SPECSDEBUG
768 // if p is behind of both surfaces, abort. 770 // if p is behind of both surfaces, abort.
769 if (distToanm * distTocmn > 0 && distToanm 771 if (distToanm * distTocmn > 0 && distToanm < 0)
770 { 772 {
771 G4Exception("G4TwistTubsSide::DistanceToP 773 G4Exception("G4TwistTubsSide::DistanceToPlane()",
772 "GeomSolids0003", FatalExcept 774 "GeomSolids0003", FatalException,
773 "Point p is behind the surfac 775 "Point p is behind the surfaces.");
774 } 776 }
775 #endif 777 #endif
776 // if p is on surface, return 0. 778 // if p is on surface, return 0.
777 if (std::fabs(distToanm) <= halftol) 779 if (std::fabs(distToanm) <= halftol)
778 { 780 {
779 xx = xxanm; 781 xx = xxanm;
780 n = nanm * parity; 782 n = nanm * parity;
781 return 0; 783 return 0;
782 } 784 }
783 else if (std::fabs(distTocmn) <= halftol) 785 else if (std::fabs(distTocmn) <= halftol)
784 { 786 {
785 xx = xxcmn; 787 xx = xxcmn;
786 n = ncmn * parity; 788 n = ncmn * parity;
787 return 0; 789 return 0;
788 } 790 }
789 791
790 if (distToanm <= distTocmn) 792 if (distToanm <= distTocmn)
791 { 793 {
792 if (distToanm > 0) 794 if (distToanm > 0)
793 { 795 {
794 // both distanses are positive. take 796 // both distanses are positive. take smaller one.
795 xx = xxanm; 797 xx = xxanm;
796 n = nanm * parity; 798 n = nanm * parity;
797 return distToanm; 799 return distToanm;
798 } 800 }
799 else 801 else
800 { 802 {
801 // take -ve distance and call the fun 803 // take -ve distance and call the function recursively.
802 return DistanceToPlane(p, A, M, N, D, 804 return DistanceToPlane(p, A, M, N, D, parity, xx, n);
803 } 805 }
804 } 806 }
805 else 807 else
806 { 808 {
807 if (distTocmn > 0) 809 if (distTocmn > 0)
808 { 810 {
809 // both distanses are positive. take 811 // both distanses are positive. take smaller one.
810 xx = xxcmn; 812 xx = xxcmn;
811 n = ncmn * parity; 813 n = ncmn * parity;
812 return distTocmn; 814 return distTocmn;
813 } 815 }
814 else 816 else
815 { 817 {
816 // take -ve distance and call the fun 818 // take -ve distance and call the function recursively.
817 return DistanceToPlane(p, C, N, M, B, 819 return DistanceToPlane(p, C, N, M, B, parity, xx, n);
818 } 820 }
819 } 821 }
820 } 822 }
821 823
822 //============================================ 824 //=====================================================================
823 //* GetAreaCode ------------------------------ 825 //* GetAreaCode -------------------------------------------------------
824 826
825 G4int G4TwistTubsSide::GetAreaCode(const G4Thr 827 G4int G4TwistTubsSide::GetAreaCode(const G4ThreeVector& xx,
826 G4boo 828 G4bool withTol)
827 { 829 {
828 // We must use the function in local coordi 830 // We must use the function in local coordinate system.
829 // See the description of DistanceToSurface 831 // See the description of DistanceToSurface(p,v).
830 832
831 const G4double ctol = 0.5 * kCarTolerance; 833 const G4double ctol = 0.5 * kCarTolerance;
832 G4int areacode = sInside; 834 G4int areacode = sInside;
833 835
834 if (fAxis[0] == kXAxis && fAxis[1] == kZAxi 836 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis)
835 { 837 {
836 G4int xaxis = 0; 838 G4int xaxis = 0;
837 G4int zaxis = 1; 839 G4int zaxis = 1;
838 840
839 if (withTol) 841 if (withTol)
840 { 842 {
841 G4bool isoutside = false; 843 G4bool isoutside = false;
842 844
843 // test boundary of xaxis 845 // test boundary of xaxis
844 846
845 if (xx.x() < fAxisMin[xaxis] + ctol) 847 if (xx.x() < fAxisMin[xaxis] + ctol)
846 { 848 {
847 areacode |= (sAxis0 & (sAxisX | sA 849 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
848 if (xx.x() <= fAxisMin[xaxis] - ct 850 if (xx.x() <= fAxisMin[xaxis] - ctol) isoutside = true;
849 851
850 } 852 }
851 else if (xx.x() > fAxisMax[xaxis] - c 853 else if (xx.x() > fAxisMax[xaxis] - ctol)
852 { 854 {
853 areacode |= (sAxis0 & (sAxisX | sA 855 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
854 if (xx.x() >= fAxisMax[xaxis] + ct 856 if (xx.x() >= fAxisMax[xaxis] + ctol) isoutside = true;
855 } 857 }
856 858
857 // test boundary of z-axis 859 // test boundary of z-axis
858 860
859 if (xx.z() < fAxisMin[zaxis] + ctol) 861 if (xx.z() < fAxisMin[zaxis] + ctol)
860 { 862 {
861 areacode |= (sAxis1 & (sAxisZ | sA 863 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
862 864
863 if ((areacode & sBoundary) != 0) << 865 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner
864 else areaco 866 else areacode |= sBoundary;
865 if (xx.z() <= fAxisMin[zaxis] - ct 867 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
866 868
867 } 869 }
868 else if (xx.z() > fAxisMax[zaxis] - c 870 else if (xx.z() > fAxisMax[zaxis] - ctol)
869 { 871 {
870 areacode |= (sAxis1 & (sAxisZ | sA 872 areacode |= (sAxis1 & (sAxisZ | sAxisMax));
871 873
872 if ((areacode & sBoundary) != 0) << 874 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner
873 else areaco 875 else areacode |= sBoundary;
874 if (xx.z() >= fAxisMax[zaxis] + ct 876 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
875 } 877 }
876 878
877 // if isoutside = true, clear inside 879 // if isoutside = true, clear inside bit.
878 // if not on boundary, add axis infor 880 // if not on boundary, add axis information.
879 881
880 if (isoutside) 882 if (isoutside)
881 { 883 {
882 G4int tmpareacode = areacode & (~s 884 G4int tmpareacode = areacode & (~sInside);
883 areacode = tmpareacode; 885 areacode = tmpareacode;
884 } 886 }
885 else if ((areacode & sBoundary) != sB 887 else if ((areacode & sBoundary) != sBoundary)
886 { 888 {
887 areacode |= (sAxis0 & sAxisX) | (s 889 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
888 } 890 }
889 } 891 }
890 else 892 else
891 { 893 {
892 // boundary of x-axis 894 // boundary of x-axis
893 895
894 if (xx.x() < fAxisMin[xaxis] ) 896 if (xx.x() < fAxisMin[xaxis] )
895 { 897 {
896 areacode |= (sAxis0 & (sAxisX | sA 898 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
897 } 899 }
898 else if (xx.x() > fAxisMax[xaxis]) 900 else if (xx.x() > fAxisMax[xaxis])
899 { 901 {
900 areacode |= (sAxis0 & (sAxisX | sA 902 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
901 } 903 }
902 904
903 // boundary of z-axis 905 // boundary of z-axis
904 906
905 if (xx.z() < fAxisMin[zaxis]) 907 if (xx.z() < fAxisMin[zaxis])
906 { 908 {
907 areacode |= (sAxis1 & (sAxisZ | sA 909 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
908 if ((areacode & sBoundary) != 0) << 910 if (areacode & sBoundary) areacode |= sCorner; // xx is oncorner
909 else areaco 911 else areacode |= sBoundary;
910 912
911 } 913 }
912 else if (xx.z() > fAxisMax[zaxis]) 914 else if (xx.z() > fAxisMax[zaxis])
913 { 915 {
914 areacode |= (sAxis1 & (sAxisZ | sA 916 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
915 if ((areacode & sBoundary) != 0) << 917 if (areacode & sBoundary) areacode |= sCorner; // xx is on corner
916 else areaco 918 else areacode |= sBoundary;
917 } 919 }
918 920
919 if ((areacode & sBoundary) != sBounda 921 if ((areacode & sBoundary) != sBoundary)
920 { 922 {
921 areacode |= (sAxis0 & sAxisX) | (s 923 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
922 } 924 }
923 } 925 }
924 return areacode; 926 return areacode;
925 } 927 }
926 else 928 else
927 { 929 {
928 G4Exception("G4TwistTubsSide::GetAreaCod 930 G4Exception("G4TwistTubsSide::GetAreaCode()",
929 "GeomSolids0001", FatalExcep 931 "GeomSolids0001", FatalException,
930 "Feature NOT implemented !") 932 "Feature NOT implemented !");
931 } 933 }
932 return areacode; 934 return areacode;
933 } 935 }
934 936
935 //============================================ 937 //=====================================================================
936 //* SetCorners( arglist ) -------------------- 938 //* SetCorners( arglist ) -------------------------------------------------
937 939
938 void G4TwistTubsSide::SetCorners( G4double end 940 void G4TwistTubsSide::SetCorners( G4double endInnerRad[2],
939 G4double end 941 G4double endOuterRad[2],
940 G4double end 942 G4double endPhi[2],
941 G4double end 943 G4double endZ[2] )
942 { 944 {
943 // Set Corner points in local coodinate. 945 // Set Corner points in local coodinate.
944 946
945 if (fAxis[0] == kXAxis && fAxis[1] == kZAxi 947 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis)
946 { 948 {
947 G4int zmin = 0 ; // at -ve z 949 G4int zmin = 0 ; // at -ve z
948 G4int zmax = 1 ; // at +ve z 950 G4int zmax = 1 ; // at +ve z
949 951
950 G4double x, y, z; 952 G4double x, y, z;
951 953
952 // corner of Axis0min and Axis1min 954 // corner of Axis0min and Axis1min
953 x = endInnerRad[zmin]*std::cos(endPhi[zm 955 x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
954 y = endInnerRad[zmin]*std::sin(endPhi[zm 956 y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
955 z = endZ[zmin]; 957 z = endZ[zmin];
956 SetCorner(sC0Min1Min, x, y, z); 958 SetCorner(sC0Min1Min, x, y, z);
957 959
958 // corner of Axis0max and Axis1min 960 // corner of Axis0max and Axis1min
959 x = endOuterRad[zmin]*std::cos(endPhi[zm 961 x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
960 y = endOuterRad[zmin]*std::sin(endPhi[zm 962 y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
961 z = endZ[zmin]; 963 z = endZ[zmin];
962 SetCorner(sC0Max1Min, x, y, z); 964 SetCorner(sC0Max1Min, x, y, z);
963 965
964 // corner of Axis0max and Axis1max 966 // corner of Axis0max and Axis1max
965 x = endOuterRad[zmax]*std::cos(endPhi[zm 967 x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
966 y = endOuterRad[zmax]*std::sin(endPhi[zm 968 y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
967 z = endZ[zmax]; 969 z = endZ[zmax];
968 SetCorner(sC0Max1Max, x, y, z); 970 SetCorner(sC0Max1Max, x, y, z);
969 971
970 // corner of Axis0min and Axis1max 972 // corner of Axis0min and Axis1max
971 x = endInnerRad[zmax]*std::cos(endPhi[zm 973 x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
972 y = endInnerRad[zmax]*std::sin(endPhi[zm 974 y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
973 z = endZ[zmax]; 975 z = endZ[zmax];
974 SetCorner(sC0Min1Max, x, y, z); 976 SetCorner(sC0Min1Max, x, y, z);
975 977
976 } 978 }
977 else 979 else
978 { 980 {
979 std::ostringstream message; 981 std::ostringstream message;
980 message << "Feature NOT implemented !" < 982 message << "Feature NOT implemented !" << G4endl
981 << " fAxis[0] = " << fAxi 983 << " fAxis[0] = " << fAxis[0] << G4endl
982 << " fAxis[1] = " << fAxi 984 << " fAxis[1] = " << fAxis[1];
983 G4Exception("G4TwistTubsSide::SetCorners 985 G4Exception("G4TwistTubsSide::SetCorners()",
984 "GeomSolids0001", FatalExcep 986 "GeomSolids0001", FatalException, message);
985 } 987 }
986 } 988 }
987 989
988 //============================================ 990 //=====================================================================
989 //* SetCorners() ----------------------------- 991 //* SetCorners() ------------------------------------------------------
990 992
991 void G4TwistTubsSide::SetCorners() 993 void G4TwistTubsSide::SetCorners()
992 { 994 {
993 G4Exception("G4TwistTubsSide::SetCorners()" 995 G4Exception("G4TwistTubsSide::SetCorners()",
994 "GeomSolids0001", FatalExceptio 996 "GeomSolids0001", FatalException,
995 "Method NOT implemented !"); 997 "Method NOT implemented !");
996 } 998 }
997 999
998 //============================================ 1000 //=====================================================================
999 //* SetBoundaries() -------------------------- 1001 //* SetBoundaries() ---------------------------------------------------
1000 1002
1001 void G4TwistTubsSide::SetBoundaries() 1003 void G4TwistTubsSide::SetBoundaries()
1002 { 1004 {
1003 // Set direction-unit vector of boundary-l 1005 // Set direction-unit vector of boundary-lines in local coodinate.
1004 // 1006 //
1005 G4ThreeVector direction; 1007 G4ThreeVector direction;
1006 1008
1007 if (fAxis[0] == kXAxis && fAxis[1] == kZAx 1009 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis)
1008 { 1010 {
1009 // sAxis0 & sAxisMin 1011 // sAxis0 & sAxisMin
1010 direction = GetCorner(sC0Min1Max) - Get 1012 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
1011 direction = direction.unit(); 1013 direction = direction.unit();
1012 SetBoundary(sAxis0 & (sAxisX | sAxisMin 1014 SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction,
1013 GetCorner(sC0Min1Min), sAxi 1015 GetCorner(sC0Min1Min), sAxisZ) ;
1014 1016
1015 // sAxis0 & sAxisMax 1017 // sAxis0 & sAxisMax
1016 direction = GetCorner(sC0Max1Max) - Get 1018 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
1017 direction = direction.unit(); 1019 direction = direction.unit();
1018 SetBoundary(sAxis0 & (sAxisX | sAxisMax 1020 SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction,
1019 GetCorner(sC0Max1Min), sAxi 1021 GetCorner(sC0Max1Min), sAxisZ);
1020 1022
1021 // sAxis1 & sAxisMin 1023 // sAxis1 & sAxisMin
1022 direction = GetCorner(sC0Max1Min) - Get 1024 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
1023 direction = direction.unit(); 1025 direction = direction.unit();
1024 SetBoundary(sAxis1 & (sAxisZ | sAxisMin 1026 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
1025 GetCorner(sC0Min1Min), sAxi 1027 GetCorner(sC0Min1Min), sAxisX);
1026 1028
1027 // sAxis1 & sAxisMax 1029 // sAxis1 & sAxisMax
1028 direction = GetCorner(sC0Max1Max) - Get 1030 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
1029 direction = direction.unit(); 1031 direction = direction.unit();
1030 SetBoundary(sAxis1 & (sAxisZ | sAxisMax 1032 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
1031 GetCorner(sC0Min1Max), sAxi 1033 GetCorner(sC0Min1Max), sAxisX);
1032 1034
1033 } 1035 }
1034 else 1036 else
1035 { 1037 {
1036 std::ostringstream message; 1038 std::ostringstream message;
1037 message << "Feature NOT implemented !" 1039 message << "Feature NOT implemented !" << G4endl
1038 << " fAxis[0] = " << fAx 1040 << " fAxis[0] = " << fAxis[0] << G4endl
1039 << " fAxis[1] = " << fAx 1041 << " fAxis[1] = " << fAxis[1];
1040 G4Exception("G4TwistTubsSide::SetCorner 1042 G4Exception("G4TwistTubsSide::SetCorners()",
1041 "GeomSolids0001", FatalExce 1043 "GeomSolids0001", FatalException, message);
1042 } 1044 }
1043 } 1045 }
1044 1046
1045 //=========================================== 1047 //=====================================================================
1046 //* GetFacets() ----------------------------- 1048 //* GetFacets() -------------------------------------------------------
1047 1049
1048 void G4TwistTubsSide::GetFacets( G4int k, G4i 1050 void G4TwistTubsSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
1049 G4int faces[ 1051 G4int faces[][4], G4int iside )
1050 { 1052 {
1051 G4double z ; // the two parameters for 1053 G4double z ; // the two parameters for the surface equation
1052 G4double x,xmin,xmax ; 1054 G4double x,xmin,xmax ;
1053 1055
1054 G4ThreeVector p ; // a point on the surfac 1056 G4ThreeVector p ; // a point on the surface, given by (z,u)
1055 1057
1056 G4int nnode ; 1058 G4int nnode ;
1057 G4int nface ; 1059 G4int nface ;
1058 1060
1059 // calculate the (n-1)*(k-1) vertices 1061 // calculate the (n-1)*(k-1) vertices
1060 1062
1061 for ( G4int i = 0 ; i<n ; ++i ) 1063 for ( G4int i = 0 ; i<n ; ++i )
1062 { 1064 {
1063 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin 1065 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
1064 1066
1065 for ( G4int j = 0 ; j<k ; ++j ) 1067 for ( G4int j = 0 ; j<k ; ++j )
1066 { 1068 {
1067 nnode = GetNode(i,j,k,n,iside) ; 1069 nnode = GetNode(i,j,k,n,iside) ;
1068 1070
1069 xmin = GetBoundaryMin(z) ; 1071 xmin = GetBoundaryMin(z) ;
1070 xmax = GetBoundaryMax(z) ; 1072 xmax = GetBoundaryMax(z) ;
1071 1073
1072 if (fHandedness < 0) 1074 if (fHandedness < 0)
1073 { 1075 {
1074 x = xmin + j*(xmax-xmin)/(k-1) ; 1076 x = xmin + j*(xmax-xmin)/(k-1) ;
1075 } 1077 }
1076 else 1078 else
1077 { 1079 {
1078 x = xmax - j*(xmax-xmin)/(k-1) ; 1080 x = xmax - j*(xmax-xmin)/(k-1) ;
1079 } 1081 }
1080 1082
1081 p = SurfacePoint(x,z,true) ; // surfac 1083 p = SurfacePoint(x,z,true) ; // surface point in global coord.system
1082 1084
1083 xyz[nnode][0] = p.x() ; 1085 xyz[nnode][0] = p.x() ;
1084 xyz[nnode][1] = p.y() ; 1086 xyz[nnode][1] = p.y() ;
1085 xyz[nnode][2] = p.z() ; 1087 xyz[nnode][2] = p.z() ;
1086 1088
1087 if ( i<n-1 && j<k-1 ) // clock wise f 1089 if ( i<n-1 && j<k-1 ) // clock wise filling
1088 { 1090 {
1089 nface = GetFace(i,j,k,n,iside) ; 1091 nface = GetFace(i,j,k,n,iside) ;
1090 1092
1091 faces[nface][0] = GetEdgeVisibility(i,j,k,n 1093 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1)
1092 * ( GetNode(i ,j ,k 1094 * ( GetNode(i ,j ,k,n,iside)+1) ;
1093 faces[nface][1] = GetEdgeVisibility(i,j,k,n 1095 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1)
1094 * ( GetNode(i+1,j ,k 1096 * ( GetNode(i+1,j ,k,n,iside)+1) ;
1095 faces[nface][2] = GetEdgeVisibility(i,j,k,n 1097 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1)
1096 * ( GetNode(i+1,j+1,k 1098 * ( GetNode(i+1,j+1,k,n,iside)+1) ;
1097 faces[nface][3] = GetEdgeVisibility(i,j,k,n 1099 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1)
1098 * ( GetNode(i ,j+1,k 1100 * ( GetNode(i ,j+1,k,n,iside)+1) ;
1099 } 1101 }
1100 } 1102 }
1101 } 1103 }
1102 } 1104 }
1103 1105