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