| 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.7 2010-07-12 15:25:37 gcosmo Exp $
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), << 28 // GEANT4 tag $Name: not supported by cvs2svn $
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 { <<
55 G4Exception("G4TwistTubsSide::G4TwistTub 67 G4Exception("G4TwistTubsSide::G4TwistTubsSide()", "GeomSolids0002",
56 FatalErrorInArgument, "Shoul 68 FatalErrorInArgument, "Should swap axis0 and axis1!");
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), fKappa(0.)
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; 265 G4int vout = 0;
260 266
261 if (std::fabs(a) < DBL_MIN) << 267 if (std::fabs(a) < DBL_MIN) {
262 { << 268 if (std::fabs(b) > DBL_MIN) {
263 if (std::fabs(b) > DBL_MIN) << 269
264 { <<
265 // single solution 270 // single solution
266 271
267 distance[0] = - c / b; 272 distance[0] = - c / b;
268 xx[0] = p + distance[0]*v; 273 xx[0] = p + distance[0]*v;
269 gxx[0] = ComputeGlobalPoint(xx[0 274 gxx[0] = ComputeGlobalPoint(xx[0]);
270 275
271 if (validate == kValidateWithTol) << 276 if (validate == kValidateWithTol) {
272 { <<
273 areacode[0] = GetAreaCode(xx[0]); 277 areacode[0] = GetAreaCode(xx[0]);
274 if (!IsOutside(areacode[0])) << 278 if (!IsOutside(areacode[0])) {
275 { <<
276 if (distance[0] >= 0) isvalid[0 279 if (distance[0] >= 0) isvalid[0] = true;
277 } 280 }
278 } << 281 } else if (validate == kValidateWithoutTol) {
279 else if (validate == kValidateWithout <<
280 { <<
281 areacode[0] = GetAreaCode(xx[0], f 282 areacode[0] = GetAreaCode(xx[0], false);
282 if (IsInside(areacode[0])) << 283 if (IsInside(areacode[0])) {
283 { <<
284 if (distance[0] >= 0) isvalid[0 284 if (distance[0] >= 0) isvalid[0] = true;
285 } 285 }
286 } << 286 } else { // kDontValidate
287 else // kDontValidate <<
288 { <<
289 // we must omit x(rho,z) = rho(z=0 287 // we must omit x(rho,z) = rho(z=0) < 0
290 if (xx[0].x() > 0) << 288 if (xx[0].x() > 0) {
291 { <<
292 areacode[0] = sInside; 289 areacode[0] = sInside;
293 if (distance[0] >= 0) isvalid[0 290 if (distance[0] >= 0) isvalid[0] = true;
294 } << 291 } else {
295 else <<
296 { <<
297 distance[0] = kInfinity; 292 distance[0] = kInfinity;
298 fCurStatWithV.SetCurrentStatus( 293 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
299 294 areacode[0], isvalid[0],
300 295 0, validate, &gp, &gv);
301 return vout; 296 return vout;
302 } 297 }
303 } 298 }
304 299
305 fCurStatWithV.SetCurrentStatus(0, gxx 300 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
306 isvali 301 isvalid[0], 1, validate, &gp, &gv);
307 vout = 1; 302 vout = 1;
308 } << 303
309 else << 304 } else {
310 { <<
311 // if a=b=0 , v.y=0 and (v.x=0 && p.x 305 // 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 306 // if v.x=0 && p.x=0, no intersection unless p is on z-axis
313 // (in that case, v is paralell to 307 // (in that case, v is paralell to surface).
314 // if v.z=0 && p.z=0, no intersect 308 // if v.z=0 && p.z=0, no intersection unless p is on x-axis
315 // (in that case, v is paralell to 309 // (in that case, v is paralell to surface).
316 // return distance = infinity. 310 // return distance = infinity.
317 311
318 fCurStatWithV.SetCurrentStatus(0, gxx 312 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
319 isvali 313 isvalid[0], 0, validate, &gp, &gv);
320 } 314 }
321 } << 315
322 else if (D > DBL_MIN) << 316 } else if (D > DBL_MIN) {
323 { << 317
324 // double solutions 318 // double solutions
325 319
326 D = std::sqrt(D); 320 D = std::sqrt(D);
327 G4double factor = 0.5/a; 321 G4double factor = 0.5/a;
328 G4double tmpdist[2] = {kInfinity, k 322 G4double tmpdist[2] = {kInfinity, kInfinity};
329 G4ThreeVector tmpxx[2]; 323 G4ThreeVector tmpxx[2];
330 G4int tmpareacode[2] = {sOutside 324 G4int tmpareacode[2] = {sOutside, sOutside};
331 G4bool tmpisvalid[2] = {false, f 325 G4bool tmpisvalid[2] = {false, false};
>> 326 G4int i;
332 327
333 for (auto i=0; i<2; ++i) << 328 for (i=0; i<2; i++) {
334 { <<
335 G4double bminusD = - b - D; 329 G4double bminusD = - b - D;
336 330
337 // protection against round off error 331 // protection against round off error
338 //G4double protection = 1.0e-6; 332 //G4double protection = 1.0e-6;
339 G4double protection = 0; 333 G4double protection = 0;
340 if ( b * D < 0 && std::fabs(bminusD / << 334 if ( b * D < 0 && std::fabs(bminusD / D) < protection ) {
341 { <<
342 G4double acovbb = (a*c)/(b*b); 335 G4double acovbb = (a*c)/(b*b);
343 tmpdist[i] = - c/b * ( 1 - acovbb 336 tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
344 } << 337 } else {
345 else <<
346 { <<
347 tmpdist[i] = factor * bminusD; 338 tmpdist[i] = factor * bminusD;
348 } 339 }
349 340
350 D = -D; 341 D = -D;
351 tmpxx[i] = p + tmpdist[i]*v; 342 tmpxx[i] = p + tmpdist[i]*v;
352 343
353 if (validate == kValidateWithTol) << 344 if (validate == kValidateWithTol) {
354 { <<
355 tmpareacode[i] = GetAreaCode(tmpxx 345 tmpareacode[i] = GetAreaCode(tmpxx[i]);
356 if (!IsOutside(tmpareacode[i])) << 346 if (!IsOutside(tmpareacode[i])) {
357 { <<
358 if (tmpdist[i] >= 0) tmpisvalid 347 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
359 continue; 348 continue;
360 } 349 }
361 } << 350 } else if (validate == kValidateWithoutTol) {
362 else if (validate == kValidateWithout <<
363 { <<
364 tmpareacode[i] = GetAreaCode(tmpxx 351 tmpareacode[i] = GetAreaCode(tmpxx[i], false);
365 if (IsInside(tmpareacode[i])) << 352 if (IsInside(tmpareacode[i])) {
366 { <<
367 if (tmpdist[i] >= 0) tmpisvalid 353 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
368 continue; 354 continue;
369 } 355 }
370 } << 356 } else { // kDontValidate
371 else // kDontValidate <<
372 { <<
373 // we must choose x(rho,z) = rho(z 357 // we must choose x(rho,z) = rho(z=0) > 0
374 if (tmpxx[i].x() > 0) << 358 if (tmpxx[i].x() > 0) {
375 { <<
376 tmpareacode[i] = sInside; 359 tmpareacode[i] = sInside;
377 if (tmpdist[i] >= 0) tmpisvalid 360 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
378 continue; 361 continue;
379 } else { 362 } else {
380 tmpdist[i] = kInfinity; 363 tmpdist[i] = kInfinity;
381 continue; 364 continue;
382 } 365 }
383 } 366 }
384 } 367 }
385 368
386 if (tmpdist[0] <= tmpdist[1]) << 369 if (tmpdist[0] <= tmpdist[1]) {
387 { <<
388 distance[0] = tmpdist[0]; 370 distance[0] = tmpdist[0];
389 distance[1] = tmpdist[1]; 371 distance[1] = tmpdist[1];
390 xx[0] = tmpxx[0]; 372 xx[0] = tmpxx[0];
391 xx[1] = tmpxx[1]; 373 xx[1] = tmpxx[1];
392 gxx[0] = ComputeGlobalPoint(tmpx 374 gxx[0] = ComputeGlobalPoint(tmpxx[0]);
393 gxx[1] = ComputeGlobalPoint(tmpx 375 gxx[1] = ComputeGlobalPoint(tmpxx[1]);
394 areacode[0] = tmpareacode[0]; 376 areacode[0] = tmpareacode[0];
395 areacode[1] = tmpareacode[1]; 377 areacode[1] = tmpareacode[1];
396 isvalid[0] = tmpisvalid[0]; 378 isvalid[0] = tmpisvalid[0];
397 isvalid[1] = tmpisvalid[1]; 379 isvalid[1] = tmpisvalid[1];
398 } << 380 } else {
399 else <<
400 { <<
401 distance[0] = tmpdist[1]; 381 distance[0] = tmpdist[1];
402 distance[1] = tmpdist[0]; 382 distance[1] = tmpdist[0];
403 xx[0] = tmpxx[1]; 383 xx[0] = tmpxx[1];
404 xx[1] = tmpxx[0]; 384 xx[1] = tmpxx[0];
405 gxx[0] = ComputeGlobalPoint(tmpx 385 gxx[0] = ComputeGlobalPoint(tmpxx[1]);
406 gxx[1] = ComputeGlobalPoint(tmpx 386 gxx[1] = ComputeGlobalPoint(tmpxx[0]);
407 areacode[0] = tmpareacode[1]; 387 areacode[0] = tmpareacode[1];
408 areacode[1] = tmpareacode[0]; 388 areacode[1] = tmpareacode[0];
409 isvalid[0] = tmpisvalid[1]; 389 isvalid[0] = tmpisvalid[1];
410 isvalid[1] = tmpisvalid[0]; 390 isvalid[1] = tmpisvalid[0];
411 } 391 }
412 392
413 fCurStatWithV.SetCurrentStatus(0, gxx[0] 393 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
414 isvalid[0 394 isvalid[0], 2, validate, &gp, &gv);
415 fCurStatWithV.SetCurrentStatus(1, gxx[1] 395 fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
416 isvalid[1 396 isvalid[1], 2, validate, &gp, &gv);
417 397
418 // protection against roundoff error 398 // protection against roundoff error
419 399
420 for (G4int k=0; k<2; ++k) << 400 for (G4int k=0; k<2; k++) {
421 { <<
422 if (!isvalid[k]) continue; 401 if (!isvalid[k]) continue;
423 402
424 G4ThreeVector xxonsurface(xx[k].x(), 403 G4ThreeVector xxonsurface(xx[k].x(), fKappa * std::fabs(xx[k].x())
425 404 * xx[k].z() , xx[k].z());
426 G4double deltaY = (xx[k] - xxo 405 G4double deltaY = (xx[k] - xxonsurface).mag();
427 406
428 if ( deltaY > 0.5*kCarTolerance ) << 407 if ( deltaY > 0.5*kCarTolerance ) {
429 { << 408
430 G4int maxcount = 10; 409 G4int maxcount = 10;
431 G4int l; 410 G4int l;
432 G4double lastdeltaY = deltaY; << 411 G4double lastdeltaY = deltaY;
433 for (l=0; l<maxcount; ++l) << 412 for (l=0; l<maxcount; l++) {
434 { <<
435 G4ThreeVector surfacenormal = Get 413 G4ThreeVector surfacenormal = GetNormal(xxonsurface);
436 distance[k] = DistanceToPlaneWith 414 distance[k] = DistanceToPlaneWithV(p, v, xxonsurface,
437 415 surfacenormal, xx[k]);
438 deltaY = (xx[k] - xxonsurfac 416 deltaY = (xx[k] - xxonsurface).mag();
439 if (deltaY > lastdeltaY) { } // << 417 if (deltaY > lastdeltaY) {
>> 418
>> 419 }
440 gxx[k] = ComputeGlobalPoint( 420 gxx[k] = ComputeGlobalPoint(xx[k]);
441 421
442 if (deltaY <= 0.5*kCarTolerance) << 422 if (deltaY <= 0.5*kCarTolerance) {
443 xxonsurface.set(xx[k].x(), << 423
444 fKappa * std::fab << 424 break;
445 xx[k].z()); << 425 }
446 } << 426 xxonsurface.set(xx[k].x(),
447 if (l == maxcount) << 427 fKappa * std::fabs(xx[k].x()) * xx[k].z(),
448 { << 428 xx[k].z());
449 std::ostringstream message; << 429 }
450 message << "Exceeded maxloop coun << 430 if (l == maxcount) {
451 << " maxloop count << 431 std::ostringstream message;
452 G4Exception("G4TwistTubsFlatSide: << 432 message << "Exceeded maxloop count!" << G4endl
453 "GeomSolids0003", Fa << 433 << " maxloop count " << maxcount;
454 } << 434 G4Exception("G4TwistTubsFlatSide::DistanceToSurface(p,v)",
>> 435 "GeomSolids0003", FatalException, message);
>> 436 }
455 } 437 }
>> 438
456 } 439 }
457 vout = 2; 440 vout = 2;
458 } << 441 } else {
459 else <<
460 { <<
461 // if D<0, no solution 442 // if D<0, no solution
462 // if D=0, just grazing the surfaces, re 443 // if D=0, just grazing the surfaces, return kInfinity
463 444
464 fCurStatWithV.SetCurrentStatus(0, gxx[0] 445 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
465 isvalid[0 446 isvalid[0], 0, validate, &gp, &gv);
466 } 447 }
467 448
468 return vout; 449 return vout;
469 } 450 }
470 451
471 //============================================ 452 //=====================================================================
472 //* DistanceToSurface ------------------------ 453 //* DistanceToSurface -------------------------------------------------
473 454
474 G4int G4TwistTubsSide::DistanceToSurface(const << 455 G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector &gp,
475 << 456 G4ThreeVector gxx[],
476 << 457 G4double distance[],
477 << 458 G4int areacode[])
478 { 459 {
479 fCurStat.ResetfDone(kDontValidate, &gp); 460 fCurStat.ResetfDone(kDontValidate, &gp);
480 if (fCurStat.IsDone()) << 461 G4int i = 0;
481 { << 462 if (fCurStat.IsDone()) {
482 for (G4int i=0; i<fCurStat.GetNXX(); ++i << 463 for (i=0; i<fCurStat.GetNXX(); i++) {
483 { <<
484 gxx[i] = fCurStat.GetXX(i); 464 gxx[i] = fCurStat.GetXX(i);
485 distance[i] = fCurStat.GetDistance(i) 465 distance[i] = fCurStat.GetDistance(i);
486 areacode[i] = fCurStat.GetAreacode(i) 466 areacode[i] = fCurStat.GetAreacode(i);
487 } 467 }
488 return fCurStat.GetNXX(); 468 return fCurStat.GetNXX();
489 } << 469 } else {
490 else // initialize << 470 // initialize
491 { << 471 for (i=0; i<2; i++) {
492 for (auto i=0; i<2; ++i) <<
493 { <<
494 distance[i] = kInfinity; 472 distance[i] = kInfinity;
495 areacode[i] = sOutside; 473 areacode[i] = sOutside;
496 gxx[i].set(kInfinity, kInfinity, kInf 474 gxx[i].set(kInfinity, kInfinity, kInfinity);
497 } 475 }
498 } 476 }
499 477
500 const G4double halftol = 0.5 * kCarToleranc << 478 static const G4double halftol = 0.5 * kCarTolerance;
501 479
502 G4ThreeVector p = ComputeLocalPoint( 480 G4ThreeVector p = ComputeLocalPoint(gp);
503 G4ThreeVector xx; 481 G4ThreeVector xx;
504 G4int parity = (fKappa >= 0 ? 1 : 482 G4int parity = (fKappa >= 0 ? 1 : -1);
505 483
506 // 484 //
507 // special case! 485 // special case!
508 // If p is on surface, or 486 // If p is on surface, or
509 // p is on z-axis, 487 // p is on z-axis,
510 // return here immediatery. 488 // return here immediatery.
511 // 489 //
512 490
513 G4ThreeVector lastgxx[2]; 491 G4ThreeVector lastgxx[2];
514 for (auto i=0; i<2; ++i) << 492 for (i=0; i<2; i++) {
515 { <<
516 lastgxx[i] = fCurStatWithV.GetXX(i); 493 lastgxx[i] = fCurStatWithV.GetXX(i);
517 } 494 }
518 495
519 if ((gp - lastgxx[0]).mag() < halftol 496 if ((gp - lastgxx[0]).mag() < halftol
520 || (gp - lastgxx[1]).mag() < halftol) << 497 || (gp - lastgxx[1]).mag() < halftol) {
521 { <<
522 // last winner, or last poststep point i 498 // last winner, or last poststep point is on the surface.
523 xx = p; 499 xx = p;
524 distance[0] = 0; 500 distance[0] = 0;
525 gxx[0] = gp; 501 gxx[0] = gp;
526 502
527 G4bool isvalid = true; 503 G4bool isvalid = true;
528 fCurStat.SetCurrentStatus(0, gxx[0], dis 504 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
529 isvalid, 1, kDont 505 isvalid, 1, kDontValidate, &gp);
530 return 1; 506 return 1;
531 } 507 }
532 508
533 if (p.getRho() == 0) << 509 if (p.getRho() == 0) {
534 { <<
535 // p is on z-axis. Namely, p is on twist 510 // p is on z-axis. Namely, p is on twisted surface (invalid area).
536 // We must return here, however, returni 511 // We must return here, however, returning distance to x-minimum
537 // boundary is better than return 0-dist 512 // boundary is better than return 0-distance.
538 // 513 //
539 G4bool isvalid = true; 514 G4bool isvalid = true;
540 if (fAxis[0] == kXAxis && fAxis[1] == kZ << 515 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
541 { <<
542 distance[0] = DistanceToBoundary(sAxi 516 distance[0] = DistanceToBoundary(sAxis0 & sAxisMin, xx, p);
543 areacode[0] = sInside; 517 areacode[0] = sInside;
544 } << 518 } else {
545 else <<
546 { <<
547 distance[0] = 0; 519 distance[0] = 0;
548 xx.set(0., 0., 0.); 520 xx.set(0., 0., 0.);
549 } 521 }
550 gxx[0] = ComputeGlobalPoint(xx); 522 gxx[0] = ComputeGlobalPoint(xx);
551 fCurStat.SetCurrentStatus(0, gxx[0], dis 523 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
552 isvalid, 0, kD 524 isvalid, 0, kDontValidate, &gp);
553 return 1; 525 return 1;
554 } 526 }
555 527
556 // 528 //
557 // special case end 529 // special case end
558 // 530 //
559 531
560 // set corner points of quadrangle try area 532 // set corner points of quadrangle try area ...
561 533
562 G4ThreeVector A; // foot of normal from p 534 G4ThreeVector A; // foot of normal from p to boundary of sAxis0 & sAxisMin
563 G4ThreeVector C; // foot of normal from p 535 G4ThreeVector C; // foot of normal from p to boundary of sAxis0 & sAxisMax
564 G4ThreeVector B; // point on boundary 536 G4ThreeVector B; // point on boundary sAxis0 & sAxisMax at z = A.z()
565 G4ThreeVector D; // point on boundary 537 G4ThreeVector D; // point on boundary sAxis0 & sAxisMin at z = C.z()
566 538
567 // G4double distToA; // distance from 539 // G4double distToA; // distance from p to A
568 DistanceToBoundary(sAxis0 & sAxisMin, A, p) 540 DistanceToBoundary(sAxis0 & sAxisMin, A, p);
569 // G4double distToC; // distance from 541 // G4double distToC; // distance from p to C
570 DistanceToBoundary(sAxis0 & sAxisMax, C, p) 542 DistanceToBoundary(sAxis0 & sAxisMax, C, p);
571 543
572 // is p.z between a.z and c.z? 544 // is p.z between a.z and c.z?
573 // p.z must be bracketed a.z and c.z. 545 // p.z must be bracketed a.z and c.z.
574 if (A.z() > C.z()) << 546 if (A.z() > C.z()) {
575 { << 547 if (p.z() > A.z()) {
576 if (p.z() > A.z()) <<
577 { <<
578 A = GetBoundaryAtPZ(sAxis0 & sAxisMin 548 A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
579 } << 549 } else if (p.z() < C.z()) {
580 else if (p.z() < C.z()) <<
581 { <<
582 C = GetBoundaryAtPZ(sAxis0 & sAxisMax 550 C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
583 } 551 }
584 } << 552 } else {
585 else << 553 if (p.z() > C.z()) {
586 { <<
587 if (p.z() > C.z()) <<
588 { <<
589 C = GetBoundaryAtPZ(sAxis0 & sAxisMax 554 C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
590 } << 555 } else if (p.z() < A.z()) {
591 else if (p.z() < A.z()) <<
592 { <<
593 A = GetBoundaryAtPZ(sAxis0 & sAxisMin 556 A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
594 } 557 }
595 } 558 }
596 559
597 G4ThreeVector d[2]; // direction vecto 560 G4ThreeVector d[2]; // direction vectors of boundary
598 G4ThreeVector x0[2]; // foot of normal 561 G4ThreeVector x0[2]; // foot of normal from line to p
599 G4int btype[2]; // boundary type 562 G4int btype[2]; // boundary type
600 563
601 for (auto i=0; i<2; ++i) << 564 for (i=0; i<2; i++) {
602 { << 565 if (i == 0) {
603 if (i == 0) <<
604 { <<
605 GetBoundaryParameters((sAxis0 & sAxis 566 GetBoundaryParameters((sAxis0 & sAxisMax), d[i], x0[i], btype[i]);
606 B = x0[i] + ((A.z() - x0[i].z()) / d[ 567 B = x0[i] + ((A.z() - x0[i].z()) / d[i].z()) * d[i];
607 // x0 + t*d , d is direction unit vec 568 // x0 + t*d , d is direction unit vector.
608 } << 569 } else {
609 else <<
610 { <<
611 GetBoundaryParameters((sAxis0 & sAxis 570 GetBoundaryParameters((sAxis0 & sAxisMin), d[i], x0[i], btype[i]);
612 D = x0[i] + ((C.z() - x0[i].z()) / d[ 571 D = x0[i] + ((C.z() - x0[i].z()) / d[i].z()) * d[i];
613 } 572 }
614 } 573 }
615 574
616 // In order to set correct diagonal, swap A 575 // In order to set correct diagonal, swap A and D, C and B if needed.
617 G4ThreeVector pt(p.x(), p.y(), 0.); 576 G4ThreeVector pt(p.x(), p.y(), 0.);
618 G4double rc = std::fabs(p.x()); 577 G4double rc = std::fabs(p.x());
619 G4ThreeVector surfacevector(rc, rc * fKappa 578 G4ThreeVector surfacevector(rc, rc * fKappa * p.z(), 0.);
620 G4int pside = AmIOnLeftSide(pt, sur 579 G4int pside = AmIOnLeftSide(pt, surfacevector);
621 G4double test = (A.z() - C.z()) * par 580 G4double test = (A.z() - C.z()) * parity * pside;
622 581
623 if (test == 0) << 582 if (test == 0) {
624 { << 583 if (pside == 0) {
625 if (pside == 0) <<
626 { <<
627 // p is on surface. 584 // p is on surface.
628 xx = p; 585 xx = p;
629 distance[0] = 0; 586 distance[0] = 0;
630 gxx[0] = gp; 587 gxx[0] = gp;
631 588
632 G4bool isvalid = true; 589 G4bool isvalid = true;
633 fCurStat.SetCurrentStatus(0, gxx[0], 590 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
634 isvalid, 1, << 591 isvalid, 1, kDontValidate, &gp);
635 return 1; 592 return 1;
636 } << 593 } else {
637 else <<
638 { <<
639 // A.z = C.z(). return distance to li 594 // A.z = C.z(). return distance to line.
640 d[0] = C - A; 595 d[0] = C - A;
641 distance[0] = DistanceToLine(p, A, d[ 596 distance[0] = DistanceToLine(p, A, d[0], xx);
642 areacode[0] = sInside; 597 areacode[0] = sInside;
643 gxx[0] = ComputeGlobalPoint(xx); 598 gxx[0] = ComputeGlobalPoint(xx);
644 G4bool isvalid = true; 599 G4bool isvalid = true;
645 fCurStat.SetCurrentStatus(0, gxx[0], 600 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
646 isvalid, 1, << 601 isvalid, 1, kDontValidate, &gp);
647 return 1; 602 return 1;
648 } 603 }
649 } << 604
650 else if (test < 0) // wrong diagonal. vect << 605 } else if (test < 0) {
651 { // swap A and D, C and << 606
>> 607 // wrong diagonal. vector AC is crossing the surface!
>> 608 // swap A and D, C and B
652 G4ThreeVector tmp; 609 G4ThreeVector tmp;
653 tmp = A; 610 tmp = A;
654 A = D; 611 A = D;
655 D = tmp; 612 D = tmp;
656 tmp = C; 613 tmp = C;
657 C = B; 614 C = B;
658 B = tmp; 615 B = tmp;
659 616
660 } << 617 } else {
661 else // correct diagonal. nothing to do. << 618 // correct diagonal. nothing to do.
662 { <<
663 } 619 }
664 620
665 // Now, we chose correct diagonal. << 621 // Now, we chose correct diaglnal.
666 // First try. divide quadrangle into double 622 // First try. divide quadrangle into double triangle by diagonal and
667 // calculate distance to both surfaces. 623 // calculate distance to both surfaces.
668 624
669 G4ThreeVector xxacb; // foot of normal fr 625 G4ThreeVector xxacb; // foot of normal from plane ACB to p
670 G4ThreeVector nacb; // normal of plane A 626 G4ThreeVector nacb; // normal of plane ACD
671 G4ThreeVector xxcad; // foot of normal fr 627 G4ThreeVector xxcad; // foot of normal from plane CAD to p
672 G4ThreeVector ncad; // normal of plane C 628 G4ThreeVector ncad; // normal of plane CAD
673 G4ThreeVector AB(A.x(), A.y(), 0); 629 G4ThreeVector AB(A.x(), A.y(), 0);
674 G4ThreeVector DC(C.x(), C.y(), 0); 630 G4ThreeVector DC(C.x(), C.y(), 0);
675 631
676 G4double distToACB = G4VTwistSurface::Dista << 632 G4double distToACB = G4VTwistSurface::DistanceToPlane(p, A, C-A, AB, xxacb, nacb) * parity;
677 << 633 G4double distToCAD = G4VTwistSurface::DistanceToPlane(p, C, C-A, DC, xxcad, ncad) * parity;
678 G4double distToCAD = G4VTwistSurface::Dista << 634
679 <<
680 // if calculated distance = 0, return 635 // if calculated distance = 0, return
681 636
682 if (std::fabs(distToACB) <= halftol || std: << 637 if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol) {
683 { <<
684 xx = (std::fabs(distToACB) < std::fabs(d 638 xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
685 areacode[0] = sInside; 639 areacode[0] = sInside;
686 gxx[0] = ComputeGlobalPoint(xx); 640 gxx[0] = ComputeGlobalPoint(xx);
687 distance[0] = 0; 641 distance[0] = 0;
688 G4bool isvalid = true; 642 G4bool isvalid = true;
689 fCurStat.SetCurrentStatus(0, gxx[0], dis 643 fCurStat.SetCurrentStatus(0, gxx[0], distance[0] , areacode[0],
690 isvalid, 1, kD 644 isvalid, 1, kDontValidate, &gp);
691 return 1; 645 return 1;
692 } 646 }
693 647
694 if (distToACB * distToCAD > 0 && distToACB << 648 if (distToACB * distToCAD > 0 && distToACB < 0) {
695 { <<
696 // both distToACB and distToCAD are nega 649 // both distToACB and distToCAD are negative.
697 // divide quadrangle into double triangl 650 // divide quadrangle into double triangle by diagonal
698 G4ThreeVector normal; 651 G4ThreeVector normal;
699 distance[0] = DistanceToPlane(p, A, B, C 652 distance[0] = DistanceToPlane(p, A, B, C, D, parity, xx, normal);
700 } << 653 } else {
701 else << 654 if (distToACB * distToCAD > 0) {
702 { <<
703 if (distToACB * distToCAD > 0) <<
704 { <<
705 // both distToACB and distToCAD are p 655 // both distToACB and distToCAD are positive.
706 // Take smaller one. 656 // Take smaller one.
707 if (distToACB <= distToCAD) << 657 if (distToACB <= distToCAD) {
708 { <<
709 distance[0] = distToACB; 658 distance[0] = distToACB;
710 xx = xxacb; 659 xx = xxacb;
711 } << 660 } else {
712 else <<
713 { <<
714 distance[0] = distToCAD; 661 distance[0] = distToCAD;
715 xx = xxcad; 662 xx = xxcad;
716 } 663 }
717 } << 664 } else {
718 else <<
719 { <<
720 // distToACB * distToCAD is negative. 665 // distToACB * distToCAD is negative.
721 // take positive one 666 // take positive one
722 if (distToACB > 0) << 667 if (distToACB > 0) {
723 { <<
724 distance[0] = distToACB; 668 distance[0] = distToACB;
725 xx = xxacb; 669 xx = xxacb;
726 } << 670 } else {
727 else <<
728 { <<
729 distance[0] = distToCAD; 671 distance[0] = distToCAD;
730 xx = xxcad; 672 xx = xxcad;
731 } 673 }
732 } 674 }
>> 675
733 } 676 }
734 areacode[0] = sInside; 677 areacode[0] = sInside;
735 gxx[0] = ComputeGlobalPoint(xx); 678 gxx[0] = ComputeGlobalPoint(xx);
736 G4bool isvalid = true; 679 G4bool isvalid = true;
737 fCurStat.SetCurrentStatus(0, gxx[0], distan 680 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
738 isvalid, 1, kDont 681 isvalid, 1, kDontValidate, &gp);
739 return 1; 682 return 1;
740 } 683 }
741 684
742 //============================================ 685 //=====================================================================
743 //* DistanceToPlane -------------------------- 686 //* DistanceToPlane ---------------------------------------------------
744 687
745 G4double G4TwistTubsSide::DistanceToPlane(cons << 688 G4double G4TwistTubsSide::DistanceToPlane(const G4ThreeVector &p,
746 cons << 689 const G4ThreeVector &A,
747 cons << 690 const G4ThreeVector &B,
748 cons << 691 const G4ThreeVector &C,
749 cons << 692 const G4ThreeVector &D,
750 cons << 693 const G4int parity,
751 << 694 G4ThreeVector &xx,
752 << 695 G4ThreeVector &n)
753 { 696 {
754 const G4double halftol = 0.5 * kCarToleranc << 697 static const G4double halftol = 0.5 * kCarTolerance;
755 698
756 G4ThreeVector M = 0.5*(A + B); 699 G4ThreeVector M = 0.5*(A + B);
757 G4ThreeVector N = 0.5*(C + D); 700 G4ThreeVector N = 0.5*(C + D);
758 G4ThreeVector xxanm; // foot of normal fro 701 G4ThreeVector xxanm; // foot of normal from p to plane ANM
759 G4ThreeVector nanm; // normal of plane AN 702 G4ThreeVector nanm; // normal of plane ANM
760 G4ThreeVector xxcmn; // foot of normal fro 703 G4ThreeVector xxcmn; // foot of normal from p to plane CMN
761 G4ThreeVector ncmn; // normal of plane CM 704 G4ThreeVector ncmn; // normal of plane CMN
762 705
763 G4double distToanm = G4VTwistSurface::Dista << 706 G4double distToanm = G4VTwistSurface::DistanceToPlane(p, A, (N - A), (M - A), xxanm, nanm) * parity;
764 << 707 G4double distTocmn = G4VTwistSurface::DistanceToPlane(p, C, (M - C), (N - C), xxcmn, ncmn) * parity;
765 G4double distTocmn = G4VTwistSurface::Dista << 708
766 <<
767 #ifdef G4SPECSDEBUG <<
768 // if p is behind of both surfaces, abort. 709 // if p is behind of both surfaces, abort.
769 if (distToanm * distTocmn > 0 && distToanm << 710 if (distToanm * distTocmn > 0 && distToanm < 0) {
770 { <<
771 G4Exception("G4TwistTubsSide::DistanceToP 711 G4Exception("G4TwistTubsSide::DistanceToPlane()",
772 "GeomSolids0003", FatalExcept 712 "GeomSolids0003", FatalException,
773 "Point p is behind the surfac 713 "Point p is behind the surfaces.");
774 } 714 }
775 #endif << 715
776 // if p is on surface, return 0. 716 // if p is on surface, return 0.
777 if (std::fabs(distToanm) <= halftol) << 717 if (std::fabs(distToanm) <= halftol) {
778 { <<
779 xx = xxanm; 718 xx = xxanm;
780 n = nanm * parity; 719 n = nanm * parity;
781 return 0; 720 return 0;
782 } << 721 } else if (std::fabs(distTocmn) <= halftol) {
783 else if (std::fabs(distTocmn) <= halftol) <<
784 { <<
785 xx = xxcmn; 722 xx = xxcmn;
786 n = ncmn * parity; 723 n = ncmn * parity;
787 return 0; 724 return 0;
788 } 725 }
789 726
790 if (distToanm <= distTocmn) << 727 if (distToanm <= distTocmn) {
791 { << 728 if (distToanm > 0) {
792 if (distToanm > 0) <<
793 { <<
794 // both distanses are positive. take 729 // both distanses are positive. take smaller one.
795 xx = xxanm; 730 xx = xxanm;
796 n = nanm * parity; 731 n = nanm * parity;
797 return distToanm; 732 return distToanm;
798 } << 733 } else {
799 else <<
800 { <<
801 // take -ve distance and call the fun 734 // take -ve distance and call the function recursively.
802 return DistanceToPlane(p, A, M, N, D, 735 return DistanceToPlane(p, A, M, N, D, parity, xx, n);
803 } 736 }
804 } << 737 } else {
805 else << 738 if (distTocmn > 0) {
806 { <<
807 if (distTocmn > 0) <<
808 { <<
809 // both distanses are positive. take 739 // both distanses are positive. take smaller one.
810 xx = xxcmn; 740 xx = xxcmn;
811 n = ncmn * parity; 741 n = ncmn * parity;
812 return distTocmn; 742 return distTocmn;
813 } << 743 } else {
814 else <<
815 { <<
816 // take -ve distance and call the fun 744 // take -ve distance and call the function recursively.
817 return DistanceToPlane(p, C, N, M, B, 745 return DistanceToPlane(p, C, N, M, B, parity, xx, n);
818 } 746 }
819 } 747 }
820 } 748 }
821 749
822 //============================================ 750 //=====================================================================
823 //* GetAreaCode ------------------------------ 751 //* GetAreaCode -------------------------------------------------------
824 752
825 G4int G4TwistTubsSide::GetAreaCode(const G4Thr << 753 G4int G4TwistTubsSide::GetAreaCode(const G4ThreeVector &xx,
826 G4boo << 754 G4bool withTol)
827 { 755 {
828 // We must use the function in local coordi 756 // We must use the function in local coordinate system.
829 // See the description of DistanceToSurface 757 // See the description of DistanceToSurface(p,v).
830 758
831 const G4double ctol = 0.5 * kCarTolerance; << 759 static const G4double ctol = 0.5 * kCarTolerance;
832 G4int areacode = sInside; 760 G4int areacode = sInside;
833 761
834 if (fAxis[0] == kXAxis && fAxis[1] == kZAxi << 762 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
835 { <<
836 G4int xaxis = 0; 763 G4int xaxis = 0;
837 G4int zaxis = 1; 764 G4int zaxis = 1;
838 765
839 if (withTol) << 766 if (withTol) {
840 { << 767
841 G4bool isoutside = false; 768 G4bool isoutside = false;
842 769
843 // test boundary of xaxis 770 // test boundary of xaxis
844 771
845 if (xx.x() < fAxisMin[xaxis] + ctol) << 772 if (xx.x() < fAxisMin[xaxis] + ctol) {
846 { <<
847 areacode |= (sAxis0 & (sAxisX | sA 773 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
848 if (xx.x() <= fAxisMin[xaxis] - ct 774 if (xx.x() <= fAxisMin[xaxis] - ctol) isoutside = true;
849 775
850 } << 776 } else if (xx.x() > fAxisMax[xaxis] - ctol) {
851 else if (xx.x() > fAxisMax[xaxis] - c <<
852 { <<
853 areacode |= (sAxis0 & (sAxisX | sA 777 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
854 if (xx.x() >= fAxisMax[xaxis] + ct 778 if (xx.x() >= fAxisMax[xaxis] + ctol) isoutside = true;
855 } 779 }
856 780
857 // test boundary of z-axis 781 // test boundary of z-axis
858 782
859 if (xx.z() < fAxisMin[zaxis] + ctol) << 783 if (xx.z() < fAxisMin[zaxis] + ctol) {
860 { <<
861 areacode |= (sAxis1 & (sAxisZ | sA 784 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
862 785
863 if ((areacode & sBoundary) != 0) << 786 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
864 else areaco 787 else areacode |= sBoundary;
865 if (xx.z() <= fAxisMin[zaxis] - ct 788 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
866 789
867 } << 790 } else if (xx.z() > fAxisMax[zaxis] - ctol) {
868 else if (xx.z() > fAxisMax[zaxis] - c <<
869 { <<
870 areacode |= (sAxis1 & (sAxisZ | sA 791 areacode |= (sAxis1 & (sAxisZ | sAxisMax));
871 792
872 if ((areacode & sBoundary) != 0) << 793 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
873 else areaco 794 else areacode |= sBoundary;
874 if (xx.z() >= fAxisMax[zaxis] + ct 795 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
875 } 796 }
876 797
877 // if isoutside = true, clear inside 798 // if isoutside = true, clear inside bit.
878 // if not on boundary, add axis infor 799 // if not on boundary, add axis information.
879 800
880 if (isoutside) << 801 if (isoutside) {
881 { <<
882 G4int tmpareacode = areacode & (~s 802 G4int tmpareacode = areacode & (~sInside);
883 areacode = tmpareacode; 803 areacode = tmpareacode;
884 } << 804 } else if ((areacode & sBoundary) != sBoundary) {
885 else if ((areacode & sBoundary) != sB <<
886 { <<
887 areacode |= (sAxis0 & sAxisX) | (s 805 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
888 } << 806 }
889 } << 807
890 else << 808 } else {
891 { << 809
892 // boundary of x-axis 810 // boundary of x-axis
893 811
894 if (xx.x() < fAxisMin[xaxis] ) << 812 if (xx.x() < fAxisMin[xaxis] ) {
895 { <<
896 areacode |= (sAxis0 & (sAxisX | sA 813 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
897 } << 814 } else if (xx.x() > fAxisMax[xaxis]) {
898 else if (xx.x() > fAxisMax[xaxis]) <<
899 { <<
900 areacode |= (sAxis0 & (sAxisX | sA 815 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
901 } 816 }
902 817
903 // boundary of z-axis 818 // boundary of z-axis
904 819
905 if (xx.z() < fAxisMin[zaxis]) << 820 if (xx.z() < fAxisMin[zaxis]) {
906 { <<
907 areacode |= (sAxis1 & (sAxisZ | sA 821 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
908 if ((areacode & sBoundary) != 0) << 822 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
909 else areaco 823 else areacode |= sBoundary;
910 824
911 } << 825 } else if (xx.z() > fAxisMax[zaxis]) {
912 else if (xx.z() > fAxisMax[zaxis]) <<
913 { <<
914 areacode |= (sAxis1 & (sAxisZ | sA 826 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
915 if ((areacode & sBoundary) != 0) << 827 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
916 else areaco 828 else areacode |= sBoundary;
917 } 829 }
918 830
919 if ((areacode & sBoundary) != sBounda << 831 if ((areacode & sBoundary) != sBoundary) {
920 { <<
921 areacode |= (sAxis0 & sAxisX) | (s 832 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
922 } 833 }
923 } 834 }
924 return areacode; 835 return areacode;
925 } << 836 } else {
926 else <<
927 { <<
928 G4Exception("G4TwistTubsSide::GetAreaCod 837 G4Exception("G4TwistTubsSide::GetAreaCode()",
929 "GeomSolids0001", FatalExcep 838 "GeomSolids0001", FatalException,
930 "Feature NOT implemented !") 839 "Feature NOT implemented !");
931 } 840 }
932 return areacode; 841 return areacode;
933 } 842 }
934 843
935 //============================================ 844 //=====================================================================
936 //* SetCorners( arglist ) -------------------- 845 //* SetCorners( arglist ) -------------------------------------------------
937 846
938 void G4TwistTubsSide::SetCorners( G4double end << 847 void G4TwistTubsSide::SetCorners(
939 G4double end << 848 G4double endInnerRad[2],
940 G4double end << 849 G4double endOuterRad[2],
941 G4double end << 850 G4double endPhi[2],
>> 851 G4double endZ[2])
942 { 852 {
943 // Set Corner points in local coodinate. 853 // Set Corner points in local coodinate.
944 854
945 if (fAxis[0] == kXAxis && fAxis[1] == kZAxi << 855 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
946 { << 856
947 G4int zmin = 0 ; // at -ve z 857 G4int zmin = 0 ; // at -ve z
948 G4int zmax = 1 ; // at +ve z 858 G4int zmax = 1 ; // at +ve z
949 859
950 G4double x, y, z; 860 G4double x, y, z;
951 861
952 // corner of Axis0min and Axis1min 862 // corner of Axis0min and Axis1min
953 x = endInnerRad[zmin]*std::cos(endPhi[zm 863 x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
954 y = endInnerRad[zmin]*std::sin(endPhi[zm 864 y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
955 z = endZ[zmin]; 865 z = endZ[zmin];
956 SetCorner(sC0Min1Min, x, y, z); 866 SetCorner(sC0Min1Min, x, y, z);
957 867
958 // corner of Axis0max and Axis1min 868 // corner of Axis0max and Axis1min
959 x = endOuterRad[zmin]*std::cos(endPhi[zm 869 x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
960 y = endOuterRad[zmin]*std::sin(endPhi[zm 870 y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
961 z = endZ[zmin]; 871 z = endZ[zmin];
962 SetCorner(sC0Max1Min, x, y, z); 872 SetCorner(sC0Max1Min, x, y, z);
963 873
964 // corner of Axis0max and Axis1max 874 // corner of Axis0max and Axis1max
965 x = endOuterRad[zmax]*std::cos(endPhi[zm 875 x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
966 y = endOuterRad[zmax]*std::sin(endPhi[zm 876 y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
967 z = endZ[zmax]; 877 z = endZ[zmax];
968 SetCorner(sC0Max1Max, x, y, z); 878 SetCorner(sC0Max1Max, x, y, z);
969 879
970 // corner of Axis0min and Axis1max 880 // corner of Axis0min and Axis1max
971 x = endInnerRad[zmax]*std::cos(endPhi[zm 881 x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
972 y = endInnerRad[zmax]*std::sin(endPhi[zm 882 y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
973 z = endZ[zmax]; 883 z = endZ[zmax];
974 SetCorner(sC0Min1Max, x, y, z); 884 SetCorner(sC0Min1Max, x, y, z);
975 885
976 } << 886 } else {
977 else <<
978 { <<
979 std::ostringstream message; 887 std::ostringstream message;
980 message << "Feature NOT implemented !" < 888 message << "Feature NOT implemented !" << G4endl
981 << " fAxis[0] = " << fAxi 889 << " fAxis[0] = " << fAxis[0] << G4endl
982 << " fAxis[1] = " << fAxi 890 << " fAxis[1] = " << fAxis[1];
983 G4Exception("G4TwistTubsSide::SetCorners 891 G4Exception("G4TwistTubsSide::SetCorners()",
984 "GeomSolids0001", FatalExcep 892 "GeomSolids0001", FatalException, message);
985 } 893 }
986 } 894 }
987 895
988 //============================================ 896 //=====================================================================
989 //* SetCorners() ----------------------------- 897 //* SetCorners() ------------------------------------------------------
990 898
991 void G4TwistTubsSide::SetCorners() 899 void G4TwistTubsSide::SetCorners()
992 { 900 {
993 G4Exception("G4TwistTubsSide::SetCorners()" 901 G4Exception("G4TwistTubsSide::SetCorners()",
994 "GeomSolids0001", FatalExceptio 902 "GeomSolids0001", FatalException,
995 "Method NOT implemented !"); 903 "Method NOT implemented !");
996 } 904 }
997 905
998 //============================================ 906 //=====================================================================
999 //* SetBoundaries() -------------------------- 907 //* SetBoundaries() ---------------------------------------------------
1000 908
1001 void G4TwistTubsSide::SetBoundaries() 909 void G4TwistTubsSide::SetBoundaries()
1002 { 910 {
1003 // Set direction-unit vector of boundary-l 911 // Set direction-unit vector of boundary-lines in local coodinate.
1004 // 912 //
1005 G4ThreeVector direction; 913 G4ThreeVector direction;
1006 914
1007 if (fAxis[0] == kXAxis && fAxis[1] == kZAx << 915 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
1008 { << 916
1009 // sAxis0 & sAxisMin 917 // sAxis0 & sAxisMin
1010 direction = GetCorner(sC0Min1Max) - Get 918 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
1011 direction = direction.unit(); 919 direction = direction.unit();
1012 SetBoundary(sAxis0 & (sAxisX | sAxisMin 920 SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction,
1013 GetCorner(sC0Min1Min), sAxi 921 GetCorner(sC0Min1Min), sAxisZ) ;
1014 922
1015 // sAxis0 & sAxisMax 923 // sAxis0 & sAxisMax
1016 direction = GetCorner(sC0Max1Max) - Get 924 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
1017 direction = direction.unit(); 925 direction = direction.unit();
1018 SetBoundary(sAxis0 & (sAxisX | sAxisMax 926 SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction,
1019 GetCorner(sC0Max1Min), sAxi 927 GetCorner(sC0Max1Min), sAxisZ);
1020 928
1021 // sAxis1 & sAxisMin 929 // sAxis1 & sAxisMin
1022 direction = GetCorner(sC0Max1Min) - Get 930 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
1023 direction = direction.unit(); 931 direction = direction.unit();
1024 SetBoundary(sAxis1 & (sAxisZ | sAxisMin 932 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
1025 GetCorner(sC0Min1Min), sAxi 933 GetCorner(sC0Min1Min), sAxisX);
1026 934
1027 // sAxis1 & sAxisMax 935 // sAxis1 & sAxisMax
1028 direction = GetCorner(sC0Max1Max) - Get 936 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
1029 direction = direction.unit(); 937 direction = direction.unit();
1030 SetBoundary(sAxis1 & (sAxisZ | sAxisMax 938 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
1031 GetCorner(sC0Min1Max), sAxi 939 GetCorner(sC0Min1Max), sAxisX);
1032 940
1033 } << 941 } else {
1034 else <<
1035 { <<
1036 std::ostringstream message; 942 std::ostringstream message;
1037 message << "Feature NOT implemented !" 943 message << "Feature NOT implemented !" << G4endl
1038 << " fAxis[0] = " << fAx 944 << " fAxis[0] = " << fAxis[0] << G4endl
1039 << " fAxis[1] = " << fAx 945 << " fAxis[1] = " << fAxis[1];
1040 G4Exception("G4TwistTubsSide::SetCorner 946 G4Exception("G4TwistTubsSide::SetCorners()",
1041 "GeomSolids0001", FatalExce 947 "GeomSolids0001", FatalException, message);
1042 } 948 }
1043 } 949 }
1044 950
1045 //=========================================== 951 //=====================================================================
1046 //* GetFacets() ----------------------------- 952 //* GetFacets() -------------------------------------------------------
1047 953
1048 void G4TwistTubsSide::GetFacets( G4int k, G4i << 954 void G4TwistTubsSide::GetFacets( G4int m, G4int n, G4double xyz[][3],
1049 G4int faces[ 955 G4int faces[][4], G4int iside )
1050 { 956 {
>> 957
1051 G4double z ; // the two parameters for 958 G4double z ; // the two parameters for the surface equation
1052 G4double x,xmin,xmax ; 959 G4double x,xmin,xmax ;
1053 960
1054 G4ThreeVector p ; // a point on the surfac 961 G4ThreeVector p ; // a point on the surface, given by (z,u)
1055 962
1056 G4int nnode ; 963 G4int nnode ;
1057 G4int nface ; 964 G4int nface ;
1058 965
1059 // calculate the (n-1)*(k-1) vertices << 966 // calculate the (n-1)*(m-1) vertices
1060 967
1061 for ( G4int i = 0 ; i<n ; ++i ) << 968 G4int i,j ;
>> 969
>> 970 for ( i = 0 ; i<n ; i++ )
1062 { 971 {
>> 972
1063 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin 973 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
1064 974
1065 for ( G4int j = 0 ; j<k ; ++j ) << 975 for ( j = 0 ; j<m ; j++ ) {
1066 { << 976
1067 nnode = GetNode(i,j,k,n,iside) ; << 977 nnode = GetNode(i,j,m,n,iside) ;
1068 978
1069 xmin = GetBoundaryMin(z) ; 979 xmin = GetBoundaryMin(z) ;
1070 xmax = GetBoundaryMax(z) ; 980 xmax = GetBoundaryMax(z) ;
1071 981
1072 if (fHandedness < 0) << 982 if (fHandedness < 0) {
1073 { << 983 x = xmin + j*(xmax-xmin)/(m-1) ;
1074 x = xmin + j*(xmax-xmin)/(k-1) ; << 984 } else {
1075 } << 985 x = xmax - j*(xmax-xmin)/(m-1) ;
1076 else <<
1077 { <<
1078 x = xmax - j*(xmax-xmin)/(k-1) ; <<
1079 } 986 }
1080 987
1081 p = SurfacePoint(x,z,true) ; // surfac 988 p = SurfacePoint(x,z,true) ; // surface point in global coord.system
1082 989
1083 xyz[nnode][0] = p.x() ; 990 xyz[nnode][0] = p.x() ;
1084 xyz[nnode][1] = p.y() ; 991 xyz[nnode][1] = p.y() ;
1085 xyz[nnode][2] = p.z() ; 992 xyz[nnode][2] = p.z() ;
1086 993
1087 if ( i<n-1 && j<k-1 ) // clock wise f << 994 if ( i<n-1 && j<m-1 ) { // clock wise filling
1088 { << 995
1089 nface = GetFace(i,j,k,n,iside) ; << 996 nface = GetFace(i,j,m,n,iside) ;
1090 << 997
1091 faces[nface][0] = GetEdgeVisibility(i,j,k,n << 998 faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,1) * ( GetNode(i ,j ,m,n,iside)+1) ;
1092 * ( GetNode(i ,j ,k << 999 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 << 1000 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 << 1001 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 << 1002
1096 * ( GetNode(i+1,j+1,k <<
1097 faces[nface][3] = GetEdgeVisibility(i,j,k,n <<
1098 * ( GetNode(i ,j+1,k <<
1099 } 1003 }
1100 } 1004 }
1101 } 1005 }
1102 } 1006 }
1103 1007