source: trunk/source/geometry/solids/specific/src/G4TwistTubsHypeSide.cc@ 1202

Last change on this file since 1202 was 1058, checked in by garnier, 17 years ago

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