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