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