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