| 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 // G4TwistTubsHypeSide 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 "G4TwistTubsHypeSide.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4GeometryTolerance.hh"
36
37 //=====================================================================
38 //* constructors ------------------------------------------------------
39
40 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const G4String& name,
41 const G4RotationMatrix& rot,
42 const G4ThreeVector& tlate,
43 const G4int handedness,
44 const G4double kappa,
45 const G4double tanstereo,
46 const G4double r0,
47 const EAxis axis0,
48 const EAxis axis1,
49 G4double axis0min,
50 G4double axis1min,
51 G4double axis0max,
52 G4double axis1max )
53 : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
54 axis0min, axis1min, axis0max, axis1max),
55 fKappa(kappa), fTanStereo(tanstereo),
56 fTan2Stereo(tanstereo*tanstereo), fR0(r0), fR02(r0*r0), fDPhi(twopi)
57 {
58 if ( (axis0 == kZAxis) && (axis1 == kPhi) )
59 {
60 G4Exception("G4TwistTubsHypeSide::G4TwistTubsHypeSide()",
61 "GeomSolids0002", FatalErrorInArgument,
62 "Should swap axis0 and axis1!");
63 }
64 fInside.gp.set(kInfinity, kInfinity, kInfinity);
65 fInside.inside = kOutside;
66 fIsValidNorm = false;
67 SetCorners();
68 SetBoundaries();
69 }
70
71 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const G4String& name,
72 G4double EndInnerRadius[2],
73 G4double EndOuterRadius[2],
74 G4double DPhi,
75 G4double EndPhi[2],
76 G4double EndZ[2],
77 G4double InnerRadius,
78 G4double OuterRadius,
79 G4double Kappa,
80 G4double TanInnerStereo,
81 G4double TanOuterStereo,
82 G4int handedness)
83 : G4VTwistSurface(name)
84 {
85
86 fHandedness = handedness; // +z = +ve, -z = -ve
87 fAxis[0] = kPhi;
88 fAxis[1] = kZAxis;
89 fAxisMin[0] = kInfinity; // we cannot fix boundary min of Phi,
90 fAxisMax[0] = kInfinity; // because it depends on z.
91 fAxisMin[1] = EndZ[0];
92 fAxisMax[1] = EndZ[1];
93 fKappa = Kappa;
94 fDPhi = DPhi ;
95
96 if (handedness < 0) // inner hyperbolic surface
97 {
98 fTanStereo = TanInnerStereo;
99 fR0 = InnerRadius;
100 }
101 else // outer hyperbolic surface
102 {
103 fTanStereo = TanOuterStereo;
104 fR0 = OuterRadius;
105 }
106 fTan2Stereo = fTanStereo * fTanStereo;
107 fR02 = fR0 * fR0;
108
109 fTrans.set(0, 0, 0);
110 fIsValidNorm = false;
111
112 fInside.gp.set(kInfinity, kInfinity, kInfinity);
113 fInside.inside = kOutside;
114
115 SetCorners(EndInnerRadius, EndOuterRadius, DPhi, EndPhi, EndZ) ;
116
117 SetBoundaries();
118 }
119
120 //=====================================================================
121 //* Fake default constructor ------------------------------------------
122
123 G4TwistTubsHypeSide::G4TwistTubsHypeSide( __void__& a )
124 : G4VTwistSurface(a)
125 {
126 }
127
128 //=====================================================================
129 //* destructor --------------------------------------------------------
130
131 G4TwistTubsHypeSide::~G4TwistTubsHypeSide() = default;
132
133 //=====================================================================
134 //* GetNormal ---------------------------------------------------------
135
136 G4ThreeVector G4TwistTubsHypeSide::GetNormal(const G4ThreeVector& tmpxx,
137 G4bool isGlobal)
138 {
139 // GetNormal returns a normal vector at a surface (or very close
140 // to surface) point at tmpxx.
141 // If isGlobal=true, it returns the normal in global coordinate.
142
143 G4ThreeVector xx;
144 if (isGlobal)
145 {
146 xx = ComputeLocalPoint(tmpxx);
147 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
148 {
149 return ComputeGlobalDirection(fCurrentNormal.normal);
150 }
151 }
152 else
153 {
154 xx = tmpxx;
155 if (xx == fCurrentNormal.p)
156 {
157 return fCurrentNormal.normal;
158 }
159 }
160
161 fCurrentNormal.p = xx;
162
163 G4ThreeVector normal( xx.x(), xx.y(), -xx.z() * fTan2Stereo);
164 normal *= fHandedness;
165 normal = normal.unit();
166
167 if (isGlobal)
168 {
169 fCurrentNormal.normal = ComputeLocalDirection(normal);
170 }
171 else
172 {
173 fCurrentNormal.normal = normal;
174 }
175 return fCurrentNormal.normal;
176 }
177
178 //=====================================================================
179 //* Inside() ----------------------------------------------------------
180
181 EInside G4TwistTubsHypeSide::Inside(const G4ThreeVector& gp)
182 {
183 // Inside returns
184 const G4double halftol
185 = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
186
187 if (fInside.gp == gp)
188 {
189 return fInside.inside;
190 }
191 fInside.gp = gp;
192
193 G4ThreeVector p = ComputeLocalPoint(gp);
194
195
196 if (p.mag() < DBL_MIN)
197 {
198 fInside.inside = kOutside;
199 return fInside.inside;
200 }
201
202 G4double rhohype = GetRhoAtPZ(p);
203 G4double distanceToOut = fHandedness * (rhohype - p.getRho());
204 // +ve : inside
205
206 if (distanceToOut < -halftol)
207 {
208 fInside.inside = kOutside;
209 }
210 else
211 {
212 G4int areacode = GetAreaCode(p);
213 if (IsOutside(areacode))
214 {
215 fInside.inside = kOutside;
216 }
217 else if (IsBoundary(areacode))
218 {
219 fInside.inside = kSurface;
220 }
221 else if (IsInside(areacode))
222 {
223 if (distanceToOut <= halftol)
224 {
225 fInside.inside = kSurface;
226 }
227 else
228 {
229 fInside.inside = kInside;
230 }
231 }
232 else
233 {
234 G4cout << "WARNING - G4TwistTubsHypeSide::Inside()" << G4endl
235 << " Invalid option !" << G4endl
236 << " name, areacode, distanceToOut = "
237 << GetName() << ", " << std::hex << areacode
238 << std::dec << ", " << distanceToOut << G4endl;
239 }
240 }
241
242 return fInside.inside;
243 }
244
245 //=====================================================================
246 //* DistanceToSurface -------------------------------------------------
247
248 G4int G4TwistTubsHypeSide::DistanceToSurface(const G4ThreeVector& gp,
249 const G4ThreeVector& gv,
250 G4ThreeVector gxx[],
251 G4double distance[],
252 G4int areacode[],
253 G4bool isvalid[],
254 EValidate validate)
255 {
256 // Decide if and where a line intersects with a hyperbolic
257 // surface (of infinite extent)
258 //
259 // Arguments:
260 // p - (in) Point on trajectory
261 // v - (in) Vector along trajectory
262 // r2 - (in) Square of radius at z = 0
263 // tan2phi - (in) std::tan(stereo)**2
264 // s - (out) Up to two points of intersection, where the
265 // intersection point is p + s*v, and if there are
266 // two intersections, s[0] < s[1]. May be negative.
267 // Returns:
268 // The number of intersections. If 0, the trajectory misses.
269 //
270 //
271 // Equation of a line:
272 //
273 // x = x0 + s*tx y = y0 + s*ty z = z0 + s*tz
274 //
275 // Equation of a hyperbolic surface:
276 //
277 // x**2 + y**2 = r**2 + (z*tanPhi)**2
278 //
279 // Solution is quadratic:
280 //
281 // a*s**2 + b*s + c = 0
282 //
283 // where:
284 //
285 // a = tx**2 + ty**2 - (tz*tanPhi)**2
286 //
287 // b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 )
288 //
289 // c = x0**2 + y0**2 - r**2 - (z0*tanPhi)**2
290 //
291
292 fCurStatWithV.ResetfDone(validate, &gp, &gv);
293
294 if (fCurStatWithV.IsDone())
295 {
296 for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
297 {
298 gxx[i] = fCurStatWithV.GetXX(i);
299 distance[i] = fCurStatWithV.GetDistance(i);
300 areacode[i] = fCurStatWithV.GetAreacode(i);
301 isvalid[i] = fCurStatWithV.IsValid(i);
302 }
303 return fCurStatWithV.GetNXX();
304 }
305 else // initialize
306 {
307 for (auto i=0; i<2; ++i)
308 {
309 distance[i] = kInfinity;
310 areacode[i] = sOutside;
311 isvalid[i] = false;
312 gxx[i].set(kInfinity, kInfinity, kInfinity);
313 }
314 }
315
316 G4ThreeVector p = ComputeLocalPoint(gp);
317 G4ThreeVector v = ComputeLocalDirection(gv);
318 G4ThreeVector xx[2];
319
320 //
321 // special case! p is on origin.
322 //
323
324 if (p.mag() == 0)
325 {
326 // p is origin.
327 // unique solution of 2-dimension question in r-z plane
328 // Equations:
329 // r^2 = fR02 + z^2*fTan2Stere0
330 // r = beta*z
331 // where
332 // beta = vrho / vz
333 // Solution (z value of intersection point):
334 // xxz = +- std::sqrt (fR02 / (beta^2 - fTan2Stereo))
335 //
336
337 G4double vz = v.z();
338 G4double absvz = std::fabs(vz);
339 G4double vrho = v.getRho();
340 G4double vslope = vrho/vz;
341 G4double vslope2 = vslope * vslope;
342 if (vrho == 0 || (vrho/absvz) <= (absvz*std::fabs(fTanStereo)/absvz))
343 {
344 // vz/vrho is bigger than slope of asymptonic line
345 distance[0] = kInfinity;
346 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
347 isvalid[0], 0, validate, &gp, &gv);
348 return 0;
349 }
350
351 if (vz != 0.0)
352 {
353 G4double xxz = std::sqrt(fR02 / (vslope2 - fTan2Stereo))
354 * (vz / std::fabs(vz)) ;
355 G4double t = xxz / vz;
356 xx[0].set(t*v.x(), t*v.y(), xxz);
357 }
358 else
359 {
360 // p.z = 0 && v.z =0
361 xx[0].set(v.x()*fR0, v.y()*fR0, 0); // v is a unit vector.
362 }
363 distance[0] = xx[0].mag();
364 gxx[0] = ComputeGlobalPoint(xx[0]);
365
366 if (validate == kValidateWithTol)
367 {
368 areacode[0] = GetAreaCode(xx[0]);
369 if (!IsOutside(areacode[0]))
370 {
371 if (distance[0] >= 0) isvalid[0] = true;
372 }
373 }
374 else if (validate == kValidateWithoutTol)
375 {
376 areacode[0] = GetAreaCode(xx[0], false);
377 if (IsInside(areacode[0]))
378 {
379 if (distance[0] >= 0) isvalid[0] = true;
380 }
381 }
382 else // kDontValidate
383 {
384 areacode[0] = sInside;
385 if (distance[0] >= 0) isvalid[0] = true;
386 }
387
388 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
389 isvalid[0], 1, validate, &gp, &gv);
390 return 1;
391 }
392
393 //
394 // special case end.
395 //
396
397 G4double a = v.x()*v.x() + v.y()*v.y() - v.z()*v.z()*fTan2Stereo;
398 G4double b = 2.0
399 * ( p.x() * v.x() + p.y() * v.y() - p.z() * v.z() * fTan2Stereo );
400 G4double c = p.x()*p.x() + p.y()*p.y() - fR02 - p.z()*p.z()*fTan2Stereo;
401 G4double D = b*b - 4*a*c; //discriminant
402 G4int vout = 0;
403
404 if (std::fabs(a) < DBL_MIN)
405 {
406 if (std::fabs(b) > DBL_MIN) // single solution
407 {
408 distance[0] = -c/b;
409 xx[0] = p + distance[0]*v;
410 gxx[0] = ComputeGlobalPoint(xx[0]);
411
412 if (validate == kValidateWithTol)
413 {
414 areacode[0] = GetAreaCode(xx[0]);
415 if (!IsOutside(areacode[0]))
416 {
417 if (distance[0] >= 0) isvalid[0] = true;
418 }
419 }
420 else if (validate == kValidateWithoutTol)
421 {
422 areacode[0] = GetAreaCode(xx[0], false);
423 if (IsInside(areacode[0]))
424 {
425 if (distance[0] >= 0) isvalid[0] = true;
426 }
427 }
428 else // kDontValidate
429 {
430 areacode[0] = sInside;
431 if (distance[0] >= 0) isvalid[0] = true;
432 }
433
434 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
435 isvalid[0], 1, validate, &gp, &gv);
436 vout = 1;
437 }
438 else
439 {
440 // if a=b=0 and c != 0, p is origin and v is parallel to asymptotic line
441 // if a=b=c=0, p is on surface and v is paralell to stereo wire.
442 // return distance = infinity
443
444 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
445 isvalid[0], 0, validate, &gp, &gv);
446 vout = 0;
447 }
448 }
449 else if (D > DBL_MIN) // double solutions
450 {
451 D = std::sqrt(D);
452 G4double factor = 0.5/a;
453 G4double tmpdist[2] = {kInfinity, kInfinity};
454 G4ThreeVector tmpxx[2] ;
455 G4int tmpareacode[2] = {sOutside, sOutside};
456 G4bool tmpisvalid[2] = {false, false};
457
458 for (auto i=0; i<2; ++i)
459 {
460 tmpdist[i] = factor*(-b - D);
461 D = -D;
462 tmpxx[i] = p + tmpdist[i]*v;
463
464 if (validate == kValidateWithTol)
465 {
466 tmpareacode[i] = GetAreaCode(tmpxx[i]);
467 if (!IsOutside(tmpareacode[i]))
468 {
469 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
470 continue;
471 }
472 }
473 else if (validate == kValidateWithoutTol)
474 {
475 tmpareacode[i] = GetAreaCode(tmpxx[i], false);
476 if (IsInside(tmpareacode[i]))
477 {
478 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
479 continue;
480 }
481 }
482 else // kDontValidate
483 {
484 tmpareacode[i] = sInside;
485 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
486 continue;
487 }
488 }
489
490 if (tmpdist[0] <= tmpdist[1])
491 {
492 distance[0] = tmpdist[0];
493 distance[1] = tmpdist[1];
494 xx[0] = tmpxx[0];
495 xx[1] = tmpxx[1];
496 gxx[0] = ComputeGlobalPoint(tmpxx[0]);
497 gxx[1] = ComputeGlobalPoint(tmpxx[1]);
498 areacode[0] = tmpareacode[0];
499 areacode[1] = tmpareacode[1];
500 isvalid[0] = tmpisvalid[0];
501 isvalid[1] = tmpisvalid[1];
502 }
503 else
504 {
505 distance[0] = tmpdist[1];
506 distance[1] = tmpdist[0];
507 xx[0] = tmpxx[1];
508 xx[1] = tmpxx[0];
509 gxx[0] = ComputeGlobalPoint(tmpxx[1]);
510 gxx[1] = ComputeGlobalPoint(tmpxx[0]);
511 areacode[0] = tmpareacode[1];
512 areacode[1] = tmpareacode[0];
513 isvalid[0] = tmpisvalid[1];
514 isvalid[1] = tmpisvalid[0];
515 }
516
517 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
518 isvalid[0], 2, validate, &gp, &gv);
519 fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
520 isvalid[1], 2, validate, &gp, &gv);
521 vout = 2;
522 }
523 else
524 {
525 // if D<0, no solution
526 // if D=0, just grazing the surfaces, return kInfinity
527
528 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
529 isvalid[0], 0, validate, &gp, &gv);
530 vout = 0;
531 }
532 return vout;
533 }
534
535 //=====================================================================
536 //* DistanceToSurface -------------------------------------------------
537
538 G4int G4TwistTubsHypeSide::DistanceToSurface(const G4ThreeVector& gp,
539 G4ThreeVector gxx[],
540 G4double distance[],
541 G4int areacode[])
542 {
543 // Find the approximate distance of a point of a hyperbolic surface.
544 // The distance must be an underestimate.
545 // It will also be nice (although not necessary) that the estimate is
546 // always finite no matter how close the point is.
547 //
548 // We arranged G4Hype::ApproxDistOutside and G4Hype::ApproxDistInside
549 // for this function. See these discriptions.
550
551 const G4double halftol
552 = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
553
554 fCurStat.ResetfDone(kDontValidate, &gp);
555
556 if (fCurStat.IsDone())
557 {
558 for (G4int i=0; i<fCurStat.GetNXX(); ++i)
559 {
560 gxx[i] = fCurStat.GetXX(i);
561 distance[i] = fCurStat.GetDistance(i);
562 areacode[i] = fCurStat.GetAreacode(i);
563 }
564 return fCurStat.GetNXX();
565 }
566 else // initialize
567 {
568 for (auto i=0; i<2; ++i)
569 {
570 distance[i] = kInfinity;
571 areacode[i] = sOutside;
572 gxx[i].set(kInfinity, kInfinity, kInfinity);
573 }
574 }
575
576
577 G4ThreeVector p = ComputeLocalPoint(gp);
578 G4ThreeVector xx;
579
580 //
581 // special case!
582 // If p is on surface, return distance = 0 immediatery .
583 //
584 G4ThreeVector lastgxx[2];
585 for (auto i=0; i<2; ++i)
586 {
587 lastgxx[i] = fCurStatWithV.GetXX(i);
588 }
589
590 if ((gp - lastgxx[0]).mag() < halftol || (gp - lastgxx[1]).mag() < halftol)
591 {
592 // last winner, or last poststep point is on the surface.
593 xx = p;
594 gxx[0] = gp;
595 distance[0] = 0;
596
597 G4bool isvalid = true;
598 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
599 isvalid, 1, kDontValidate, &gp);
600
601 return 1;
602 }
603 //
604 // special case end
605 //
606
607 G4double prho = p.getRho();
608 G4double pz = std::fabs(p.z()); // use symmetry
609 G4double r1 = std::sqrt(fR02 + pz * pz * fTan2Stereo);
610
611 G4ThreeVector pabsz(p.x(), p.y(), pz);
612
613 if (prho > r1 + halftol) // p is outside of Hyperbolic surface
614 {
615 // First point xx1
616 G4double t = r1 / prho;
617 G4ThreeVector xx1(t * pabsz.x(), t * pabsz.y() , pz);
618
619 // Second point xx2
620 G4double z2 = (prho * fTanStereo + pz) / (1 + fTan2Stereo);
621 G4double r2 = std::sqrt(fR02 + z2 * z2 * fTan2Stereo);
622 t = r2 / prho;
623 G4ThreeVector xx2(t * pabsz.x(), t * pabsz.y() , z2);
624
625 G4double len = (xx2 - xx1).mag();
626 if (len < DBL_MIN)
627 {
628 // xx2 = xx1?? I guess we
629 // must have really bracketed the normal
630 distance[0] = (pabsz - xx1).mag();
631 xx = xx1;
632 }
633 else
634 {
635 distance[0] = DistanceToLine(pabsz, xx1, (xx2 - xx1) , xx);
636 }
637
638 }
639 else if (prho < r1 - halftol) // p is inside of Hyperbolic surface.
640 {
641 // First point xx1
642 G4double t;
643 G4ThreeVector xx1;
644 if (prho < DBL_MIN)
645 {
646 xx1.set(r1, 0. , pz);
647 }
648 else
649 {
650 t = r1 / prho;
651 xx1.set(t * pabsz.x(), t * pabsz.y() , pz);
652 }
653
654 // dr, dz is tangential vector of Hyparbolic surface at xx1
655 // dr = r, dz = z*tan2stereo
656 G4double dr = pz * fTan2Stereo;
657 G4double dz = r1;
658 G4double tanbeta = dr / dz;
659 G4double pztanbeta = pz * tanbeta;
660
661 // Second point xx2
662 // xx2 is intersection between x-axis and tangential vector
663 G4double r2 = r1 - pztanbeta;
664 G4ThreeVector xx2;
665 if (prho < DBL_MIN)
666 {
667 xx2.set(r2, 0. , 0.);
668 }
669 else
670 {
671 t = r2 / prho;
672 xx2.set(t * pabsz.x(), t * pabsz.y() , 0.);
673 }
674
675 G4ThreeVector d = xx2 - xx1;
676 distance[0] = DistanceToLine(pabsz, xx1, d, xx);
677
678 }
679 else // p is on Hyperbolic surface.
680 {
681 distance[0] = 0;
682 xx.set(p.x(), p.y(), pz);
683 }
684
685 if (p.z() < 0)
686 {
687 G4ThreeVector tmpxx(xx.x(), xx.y(), -xx.z());
688 xx = tmpxx;
689 }
690
691 gxx[0] = ComputeGlobalPoint(xx);
692 areacode[0] = sInside;
693 G4bool isvalid = true;
694 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
695 isvalid, 1, kDontValidate, &gp);
696 return 1;
697 }
698
699 //=====================================================================
700 //* GetAreaCode -------------------------------------------------------
701
702 G4int G4TwistTubsHypeSide::GetAreaCode(const G4ThreeVector& xx,
703 G4bool withTol)
704 {
705 const G4double ctol = 0.5 * kCarTolerance;
706 G4int areacode = sInside;
707
708 if ((fAxis[0] == kPhi && fAxis[1] == kZAxis))
709 {
710 G4int zaxis = 1;
711
712 if (withTol)
713 {
714 G4bool isoutside = false;
715 G4int phiareacode = GetAreaCodeInPhi(xx);
716 G4bool isoutsideinphi = IsOutside(phiareacode);
717
718 // test boundary of phiaxis
719
720 if ((phiareacode & sAxisMin) == sAxisMin)
721 {
722 areacode |= (sAxis0 & (sAxisPhi | sAxisMin)) | sBoundary;
723 if (isoutsideinphi) isoutside = true;
724 }
725 else if ((phiareacode & sAxisMax) == sAxisMax)
726 {
727 areacode |= (sAxis0 & (sAxisPhi | sAxisMax)) | sBoundary;
728 if (isoutsideinphi) isoutside = true;
729 }
730
731 // test boundary of zaxis
732
733 if (xx.z() < fAxisMin[zaxis] + ctol)
734 {
735 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
736 if ((areacode & sBoundary) != 0) areacode |= sCorner; // xx is on corner.
737 else areacode |= sBoundary;
738
739 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
740
741 }
742 else if (xx.z() > fAxisMax[zaxis] - ctol)
743 {
744 areacode |= (sAxis1 & (sAxisZ | sAxisMax));
745 if ((areacode & sBoundary) != 0) areacode |= sCorner; // xx is on corner.
746 else areacode |= sBoundary;
747
748 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
749 }
750
751 // if isoutside = true, clear sInside bit.
752 // if not on boundary, add boundary information.
753
754 if (isoutside)
755 {
756 G4int tmpareacode = areacode & (~sInside);
757 areacode = tmpareacode;
758 }
759 else if ((areacode & sBoundary) != sBoundary)
760 {
761 areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
762 }
763
764 return areacode;
765 }
766 else
767 {
768 G4int phiareacode = GetAreaCodeInPhi(xx, false);
769
770 // test boundary of z-axis
771
772 if (xx.z() < fAxisMin[zaxis])
773 {
774 areacode |= (sAxis1 & (sAxisZ | sAxisMin)) | sBoundary;
775
776 }
777 else if (xx.z() > fAxisMax[zaxis])
778 {
779 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) | sBoundary;
780 }
781
782 // boundary of phi-axis
783
784 if (phiareacode == sAxisMin)
785 {
786 areacode |= (sAxis0 & (sAxisPhi | sAxisMin));
787 if ((areacode & sBoundary) != 0) areacode |= sCorner; // xx is on corner.
788 else areacode |= sBoundary;
789
790 }
791 else if (phiareacode == sAxisMax)
792 {
793 areacode |= (sAxis0 & (sAxisPhi | sAxisMax));
794 if ((areacode & sBoundary) != 0) areacode |= sCorner; // xx is on corner.
795 else areacode |= sBoundary;
796 }
797
798 // if not on boundary, add boundary information.
799
800 if ((areacode & sBoundary) != sBoundary)
801 {
802 areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
803 }
804 return areacode;
805 }
806 }
807 else
808 {
809 std::ostringstream message;
810 message << "Feature NOT implemented !" << G4endl
811 << " fAxis[0] = " << fAxis[0] << G4endl
812 << " fAxis[1] = " << fAxis[1];
813 G4Exception("G4TwistTubsHypeSide::GetAreaCode()",
814 "GeomSolids0001", FatalException, message);
815 }
816 return areacode;
817 }
818
819 //=====================================================================
820 //* GetAreaCodeInPhi --------------------------------------------------
821
822 G4int G4TwistTubsHypeSide::GetAreaCodeInPhi(const G4ThreeVector& xx,
823 G4bool withTol)
824 {
825
826 G4ThreeVector lowerlimit; // lower phi-boundary limit at z = xx.z()
827 G4ThreeVector upperlimit; // upper phi-boundary limit at z = xx.z()
828 lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxisMin, xx);
829 upperlimit = GetBoundaryAtPZ(sAxis0 & sAxisMax, xx);
830
831 G4int areacode = sInside;
832 G4bool isoutside = false;
833
834 if (withTol)
835 {
836 if (AmIOnLeftSide(xx, lowerlimit) >= 0) // xx is on lowerlimit
837 {
838 areacode |= (sAxisMin | sBoundary);
839 if (AmIOnLeftSide(xx, lowerlimit) > 0) isoutside = true;
840
841 }
842 else if (AmIOnLeftSide(xx, upperlimit) <= 0) // xx is on upperlimit
843 {
844 areacode |= (sAxisMax | sBoundary);
845 if (AmIOnLeftSide(xx, upperlimit) < 0) isoutside = true;
846 }
847
848 // if isoutside = true, clear inside bit.
849
850 if (isoutside)
851 {
852 G4int tmpareacode = areacode & (~sInside);
853 areacode = tmpareacode;
854 }
855 }
856 else
857 {
858 if (AmIOnLeftSide(xx, lowerlimit, false) >= 0)
859 {
860 areacode |= (sAxisMin | sBoundary);
861 }
862 else if (AmIOnLeftSide(xx, upperlimit, false) <= 0)
863 {
864 areacode |= (sAxisMax | sBoundary);
865 }
866 }
867
868 return areacode;
869 }
870
871 //=====================================================================
872 //* SetCorners(EndInnerRadius, EndOuterRadius,DPhi,EndPhi,EndZ) -------
873
874 void G4TwistTubsHypeSide::SetCorners( G4double EndInnerRadius[2],
875 G4double EndOuterRadius[2],
876 G4double DPhi,
877 G4double endPhi[2],
878 G4double endZ[2] )
879 {
880 // Set Corner points in local coodinate.
881
882 if (fAxis[0] == kPhi && fAxis[1] == kZAxis) {
883
884 G4double endRad[2];
885 G4double halfdphi = 0.5*DPhi;
886
887 for (auto i=0; i<2; ++i) // i=0,1 : -ve z, +ve z
888 {
889 endRad[i] = (fHandedness == 1 ? EndOuterRadius[i] : EndInnerRadius[i]);
890 }
891
892 G4int zmin = 0 ; // at -ve z
893 G4int zmax = 1 ; // at +ve z
894
895 G4double x, y, z;
896
897 // corner of Axis0min and Axis1min
898 x = endRad[zmin]*std::cos(endPhi[zmin] - halfdphi);
899 y = endRad[zmin]*std::sin(endPhi[zmin] - halfdphi);
900 z = endZ[zmin];
901 SetCorner(sC0Min1Min, x, y, z);
902
903 // corner of Axis0max and Axis1min
904 x = endRad[zmin]*std::cos(endPhi[zmin] + halfdphi);
905 y = endRad[zmin]*std::sin(endPhi[zmin] + halfdphi);
906 z = endZ[zmin];
907 SetCorner(sC0Max1Min, x, y, z);
908
909 // corner of Axis0max and Axis1max
910 x = endRad[zmax]*std::cos(endPhi[zmax] + halfdphi);
911 y = endRad[zmax]*std::sin(endPhi[zmax] + halfdphi);
912 z = endZ[zmax];
913 SetCorner(sC0Max1Max, x, y, z);
914
915 // corner of Axis0min and Axis1max
916 x = endRad[zmax]*std::cos(endPhi[zmax] - halfdphi);
917 y = endRad[zmax]*std::sin(endPhi[zmax] - halfdphi);
918 z = endZ[zmax];
919 SetCorner(sC0Min1Max, x, y, z);
920
921 }
922 else
923 {
924 std::ostringstream message;
925 message << "Feature NOT implemented !" << G4endl
926 << " fAxis[0] = " << fAxis[0] << G4endl
927 << " fAxis[1] = " << fAxis[1];
928 G4Exception("G4TwistTubsHypeSide::SetCorners()",
929 "GeomSolids0001", FatalException, message);
930 }
931 }
932
933 //=====================================================================
934 //* SetCorners() ------------------------------------------------------
935
936 void G4TwistTubsHypeSide::SetCorners()
937 {
938 G4Exception("G4TwistTubsHypeSide::SetCorners()",
939 "GeomSolids0001", FatalException,
940 "Method NOT implemented !");
941 }
942
943 //=====================================================================
944 //* SetBoundaries() ---------------------------------------------------
945
946 void G4TwistTubsHypeSide::SetBoundaries()
947 {
948 // Set direction-unit vector of phi-boundary-lines in local coodinate.
949 // sAxis0 must be kPhi.
950 // This fanction set lower phi-boundary and upper phi-boundary.
951
952 if (fAxis[0] == kPhi && fAxis[1] == kZAxis)
953 {
954 G4ThreeVector direction;
955 // sAxis0 & sAxisMin
956 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
957 direction = direction.unit();
958 SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction,
959 GetCorner(sC0Min1Min), sAxisZ);
960
961 // sAxis0 & sAxisMax
962 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
963 direction = direction.unit();
964 SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
965 GetCorner(sC0Max1Min), sAxisZ);
966
967 // sAxis1 & sAxisMin
968 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
969 direction = direction.unit();
970 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
971 GetCorner(sC0Min1Min), sAxisPhi);
972
973 // sAxis1 & sAxisMax
974 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
975 direction = direction.unit();
976 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
977 GetCorner(sC0Min1Max), sAxisPhi);
978 }
979 else
980 {
981 std::ostringstream message;
982 message << "Feature NOT implemented !" << G4endl
983 << " fAxis[0] = " << fAxis[0] << G4endl
984 << " fAxis[1] = " << fAxis[1];
985 G4Exception("G4TwistTubsHypeSide::SetBoundaries()",
986 "GeomSolids0001", FatalException, message);
987 }
988 }
989
990 //=====================================================================
991 //* GetFacets() -------------------------------------------------------
992
993 void G4TwistTubsHypeSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
994 G4int faces[][4], G4int iside )
995 {
996 G4double z ; // the two parameters for the surface equation
997 G4double x,xmin,xmax ;
998
999 G4ThreeVector p ; // a point on the surface, given by (z,u)
1000
1001 G4int nnode ;
1002 G4int nface ;
1003
1004 // calculate the (n-1)*(k-1) vertices
1005
1006 for ( G4int i = 0 ; i<n ; ++i )
1007 {
1008 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
1009
1010 for ( G4int j = 0 ; j<k ; ++j )
1011 {
1012 nnode = GetNode(i,j,k,n,iside) ;
1013
1014 xmin = GetBoundaryMin(z) ;
1015 xmax = GetBoundaryMax(z) ;
1016
1017 if (fHandedness < 0) // inner hyperbolic surface
1018 {
1019 x = xmin + j*(xmax-xmin)/(k-1) ;
1020 }
1021 else // outer hyperbolic surface
1022 {
1023 x = xmax - j*(xmax-xmin)/(k-1) ;
1024 }
1025
1026 p = SurfacePoint(x,z,true) ; // surface point in global coord.system
1027
1028 xyz[nnode][0] = p.x() ;
1029 xyz[nnode][1] = p.y() ;
1030 xyz[nnode][2] = p.z() ;
1031
1032 if ( i<n-1 && j<k-1 ) // clock wise filling
1033 {
1034 nface = GetFace(i,j,k,n,iside) ;
1035 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1)
1036 * ( GetNode(i ,j ,k,n,iside)+1) ;
1037 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1)
1038 * ( GetNode(i+1,j ,k,n,iside)+1) ;
1039 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1)
1040 * ( GetNode(i+1,j+1,k,n,iside)+1) ;
1041 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1)
1042 * ( GetNode(i ,j+1,k,n,iside)+1) ;
1043 }
1044 }
1045 }
1046 }
1047