Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4TwistTubsSide.cc
Go to the documentation of this file.
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 ([email protected]), created.
29// 13-Nov-2003 - O.Link ([email protected]), Integration in Geant4
30// from original version in Jupiter-2.5.02 application.
31// --------------------------------------------------------------------
32
33#include "G4TwistTubsSide.hh"
34
35//=====================================================================
36//* constructors ------------------------------------------------------
37
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
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;
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
99{
100}
101
102
103//=====================================================================
104//* destructor --------------------------------------------------------
105
107
108//=====================================================================
109//* GetNormal ---------------------------------------------------------
110
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 {
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 {
143 }
144 else
145 {
146 fCurrentNormal.normal = normal.unit();
147 }
148 return fCurrentNormal.normal;
149}
150
151//=====================================================================
152//* DistanceToSurface -------------------------------------------------
153
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
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
475 G4ThreeVector gxx[],
476 G4double distance[],
477 G4int areacode[])
478{
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
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
569 // G4double distToC; // distance from p to C
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 {
579 }
580 else if (p.z() < C.z())
581 {
583 }
584 }
585 else
586 {
587 if (p.z() > C.z())
588 {
590 }
591 else if (p.z() < A.z())
592 {
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
745G4double 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
825G4int 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
938void 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
991void G4TwistTubsSide::SetCorners()
992{
993 G4Exception("G4TwistTubsSide::SetCorners()",
994 "GeomSolids0001", FatalException,
995 "Method NOT implemented !");
996}
997
998//=====================================================================
999//* SetBoundaries() ---------------------------------------------------
1000
1001void 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,
1014
1015 // sAxis0 & sAxisMax
1016 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
1017 direction = direction.unit();
1018 SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction,
1020
1021 // sAxis1 & sAxisMin
1022 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
1023 direction = direction.unit();
1024 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
1026
1027 // sAxis1 & sAxisMax
1028 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
1029 direction = direction.unit();
1030 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
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
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}
G4double C(G4double temp)
G4double B(G4double temperature)
G4double D(G4double temp)
@ FatalException
@ FatalErrorInArgument
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
#define M(row, col)
double G4double
Definition G4Types.hh:83
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
const G4double A[17]
#define G4endl
Definition G4ios.hh:67
double z() const
Hep3Vector unit() const
double x() const
double y() const
Hep3Vector cross(const Hep3Vector &) const
double getRho() const
void set(double x, double y, double z)
HepRotation & rotateZ(double delta)
Definition Rotation.cc:87
G4double GetBoundaryMax(G4double phi) override
G4double GetBoundaryMin(G4double phi) override
G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false) override
~G4TwistTubsSide() override
void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside) override
G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false) override
G4TwistTubsSide(const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, G4int handedness, const G4double kappa, const EAxis axis0=kXAxis, const EAxis axis1=kZAxis, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol) override
G4double GetDistance(G4int i) const
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=nullptr)
G4ThreeVector GetXX(G4int i) const
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=nullptr)
virtual G4int AmIOnLeftSide(const G4ThreeVector &me, const G4ThreeVector &vec, G4bool withTol=true)
static const G4int sC0Min1Min
static const G4int sC0Min1Max
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sOutside
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
static const G4int sAxisMax
static const G4int sAxis0
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
G4RotationMatrix fRot
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
G4double DistanceToLine(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &d, G4ThreeVector &xx)
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
static const G4int sAxisMin
static const G4int sC0Max1Max
static const G4int sAxis1
virtual G4ThreeVector GetBoundaryAtPZ(G4int areacode, const G4ThreeVector &p) const
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
G4ThreeVector fTrans
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
G4ThreeVector GetCorner(G4int areacode) const
static const G4int sBoundary
static const G4int sAxisZ
G4bool IsOutside(G4int areacode) const
virtual G4double DistanceToBoundary(G4int areacode, G4ThreeVector &xx, const G4ThreeVector &p)
static const G4int sCorner
static const G4int sC0Max1Min
static const G4int sInside
CurrentStatus fCurStatWithV
static const G4int sAxisX
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
G4SurfCurNormal fCurrentNormal
virtual void GetBoundaryParameters(const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
CurrentStatus fCurStat
#define N
Definition crc32.c:57
EAxis
Definition geomdefs.hh:54
@ kXAxis
Definition geomdefs.hh:55
@ kZAxis
Definition geomdefs.hh:57
#define DBL_MIN
Definition templates.hh:54