Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TwistTubsHypeSide.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 //
27 // $Id: G4TwistTubsHypeSide.cc 72937 2013-08-14 13:20:38Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 //
34 // G4TwistTubsHypeSide.cc
35 //
36 // Author:
37 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
38 //
39 // History:
40 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
41 // from original version in Jupiter-2.5.02 application.
42 // --------------------------------------------------------------------
43 
44 #include "G4TwistTubsHypeSide.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4GeometryTolerance.hh"
47 
48 //=====================================================================
49 //* constructors ------------------------------------------------------
50 
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), fDPhi(twopi)
68 {
69  if ( (axis0 == kZAxis) && (axis1 == kPhi) )
70  {
71  G4Exception("G4TwistTubsHypeSide::G4TwistTubsHypeSide()",
72  "GeomSolids0002", FatalErrorInArgument,
73  "Should swap axis0 and axis1!");
74  }
75 
76  fInside.gp.set(kInfinity, kInfinity, kInfinity);
77  fInside.inside = kOutside;
78  fIsValidNorm = false;
79 
80  SetCorners();
81  SetBoundaries();
82 
83 }
84 
86  G4double EndInnerRadius[2],
87  G4double EndOuterRadius[2],
88  G4double DPhi,
89  G4double EndPhi[2],
90  G4double EndZ[2],
91  G4double InnerRadius,
92  G4double OuterRadius,
93  G4double Kappa,
94  G4double TanInnerStereo,
95  G4double TanOuterStereo,
96  G4int handedness)
97  : G4VTwistSurface(name)
98 {
99 
100  fHandedness = handedness; // +z = +ve, -z = -ve
101  fAxis[0] = kPhi;
102  fAxis[1] = kZAxis;
103  fAxisMin[0] = kInfinity; // we cannot fix boundary min of Phi,
104  fAxisMax[0] = kInfinity; // because it depends on z.
105  fAxisMin[1] = EndZ[0];
106  fAxisMax[1] = EndZ[1];
107  fKappa = Kappa;
108  fDPhi = DPhi ;
109 
110  if (handedness < 0) { // inner hyperbolic surface
111  fTanStereo = TanInnerStereo;
112  fR0 = InnerRadius;
113  } else { // outer hyperbolic surface
114  fTanStereo = TanOuterStereo;
115  fR0 = OuterRadius;
116  }
117  fTan2Stereo = fTanStereo * fTanStereo;
118  fR02 = fR0 * fR0;
119 
120  fTrans.set(0, 0, 0);
121  fIsValidNorm = false;
122 
123  fInside.gp.set(kInfinity, kInfinity, kInfinity);
124  fInside.inside = kOutside;
125 
126  SetCorners(EndInnerRadius, EndOuterRadius, DPhi, EndPhi, EndZ) ;
127 
128  SetBoundaries();
129 }
130 
131 //=====================================================================
132 //* Fake default constructor ------------------------------------------
133 
135  : G4VTwistSurface(a), fKappa(0.), fTanStereo(0.), fTan2Stereo(0.),
136  fR0(0.), fR02(0.), fDPhi(0.)
137 {
138 }
139 
140 //=====================================================================
141 //* destructor --------------------------------------------------------
142 
144 {
145 }
146 
147 //=====================================================================
148 //* GetNormal ---------------------------------------------------------
149 
151  G4bool isGlobal)
152 {
153  // GetNormal returns a normal vector at a surface (or very close
154  // to surface) point at tmpxx.
155  // If isGlobal=true, it returns the normal in global coordinate.
156  //
157 
158  G4ThreeVector xx;
159  if (isGlobal) {
160  xx = ComputeLocalPoint(tmpxx);
161  if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
163  }
164  } else {
165  xx = tmpxx;
166  if (xx == fCurrentNormal.p) {
167  return fCurrentNormal.normal;
168  }
169  }
170 
171  fCurrentNormal.p = xx;
172 
173  G4ThreeVector normal( xx.x(), xx.y(), -xx.z() * fTan2Stereo);
174  normal *= fHandedness;
175  normal = normal.unit();
176 
177  if (isGlobal) {
179  } else {
180  fCurrentNormal.normal = normal;
181  }
182  return fCurrentNormal.normal;
183 }
184 
185 //=====================================================================
186 //* Inside() ----------------------------------------------------------
187 
189 {
190  // Inside returns
191  const G4double halftol
193 
194  if (fInside.gp == gp) {
195  return fInside.inside;
196  }
197  fInside.gp = gp;
198 
200 
201 
202  if (p.mag() < DBL_MIN) {
203  fInside.inside = kOutside;
204  return fInside.inside;
205  }
206 
207  G4double rhohype = GetRhoAtPZ(p);
208  G4double distanceToOut = fHandedness * (rhohype - p.getRho());
209  // +ve : inside
210 
211  if (distanceToOut < -halftol) {
212 
213  fInside.inside = kOutside;
214 
215  } else {
216 
217  G4int areacode = GetAreaCode(p);
218  if (IsOutside(areacode)) {
219  fInside.inside = kOutside;
220  } else if (IsBoundary(areacode)) {
221  fInside.inside = kSurface;
222  } else if (IsInside(areacode)) {
223  if (distanceToOut <= halftol) {
224  fInside.inside = kSurface;
225  } else {
226  fInside.inside = kInside;
227  }
228  } else {
229  G4cout << "WARNING - G4TwistTubsHypeSide::Inside()" << G4endl
230  << " Invalid option !" << G4endl
231  << " name, areacode, distanceToOut = "
232  << GetName() << ", " << std::hex << areacode << std::dec << ", "
233  << distanceToOut << G4endl;
234  }
235  }
236 
237  return fInside.inside;
238 }
239 
240 //=====================================================================
241 //* DistanceToSurface -------------------------------------------------
242 
244  const G4ThreeVector &gv,
245  G4ThreeVector gxx[],
246  G4double distance[],
247  G4int areacode[],
248  G4bool isvalid[],
249  EValidate validate)
250 {
251  //
252  // Decide if and where a line intersects with a hyperbolic
253  // surface (of infinite extent)
254  //
255  // Arguments:
256  // p - (in) Point on trajectory
257  // v - (in) Vector along trajectory
258  // r2 - (in) Square of radius at z = 0
259  // tan2phi - (in) std::tan(stereo)**2
260  // s - (out) Up to two points of intersection, where the
261  // intersection point is p + s*v, and if there are
262  // two intersections, s[0] < s[1]. May be negative.
263  // Returns:
264  // The number of intersections. If 0, the trajectory misses.
265  //
266  //
267  // Equation of a line:
268  //
269  // x = x0 + s*tx y = y0 + s*ty z = z0 + s*tz
270  //
271  // Equation of a hyperbolic surface:
272  //
273  // x**2 + y**2 = r**2 + (z*tanPhi)**2
274  //
275  // Solution is quadratic:
276  //
277  // a*s**2 + b*s + c = 0
278  //
279  // where:
280  //
281  // a = tx**2 + ty**2 - (tz*tanPhi)**2
282  //
283  // b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 )
284  //
285  // c = x0**2 + y0**2 - r**2 - (z0*tanPhi)**2
286  //
287 
288  fCurStatWithV.ResetfDone(validate, &gp, &gv);
289 
290  if (fCurStatWithV.IsDone()) {
291  G4int i;
292  for (i=0; i<fCurStatWithV.GetNXX(); i++) {
293  gxx[i] = fCurStatWithV.GetXX(i);
294  distance[i] = fCurStatWithV.GetDistance(i);
295  areacode[i] = fCurStatWithV.GetAreacode(i);
296  isvalid[i] = fCurStatWithV.IsValid(i);
297  }
298  return fCurStatWithV.GetNXX();
299  } else {
300  // initialize
301  G4int i;
302  for (i=0; i<2; i++) {
303  distance[i] = kInfinity;
304  areacode[i] = sOutside;
305  isvalid[i] = false;
306  gxx[i].set(kInfinity, kInfinity, kInfinity);
307  }
308  }
309 
312  G4ThreeVector xx[2];
313 
314  //
315  // special case! p is on origin.
316  //
317 
318  if (p.mag() == 0) {
319  // p is origin.
320  // unique solution of 2-dimension question in r-z plane
321  // Equations:
322  // r^2 = fR02 + z^2*fTan2Stere0
323  // r = beta*z
324  // where
325  // beta = vrho / vz
326  // Solution (z value of intersection point):
327  // xxz = +- std::sqrt (fR02 / (beta^2 - fTan2Stereo))
328  //
329 
330  G4double vz = v.z();
331  G4double absvz = std::fabs(vz);
332  G4double vrho = v.getRho();
333  G4double vslope = vrho/vz;
334  G4double vslope2 = vslope * vslope;
335  if (vrho == 0 || (vrho/absvz) <= (absvz*std::fabs(fTanStereo)/absvz)) {
336  // vz/vrho is bigger than slope of asymptonic line
337  distance[0] = kInfinity;
338  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
339  isvalid[0], 0, validate, &gp, &gv);
340  return 0;
341  }
342 
343  if (vz) {
344  G4double xxz = std::sqrt(fR02 / (vslope2 - fTan2Stereo))
345  * (vz / std::fabs(vz)) ;
346  G4double t = xxz / vz;
347  xx[0].set(t*v.x(), t*v.y(), xxz);
348  } else {
349  // p.z = 0 && v.z =0
350  xx[0].set(v.x()*fR0, v.y()*fR0, 0); // v is a unit vector.
351  }
352  distance[0] = xx[0].mag();
353  gxx[0] = ComputeGlobalPoint(xx[0]);
354 
355  if (validate == kValidateWithTol) {
356  areacode[0] = GetAreaCode(xx[0]);
357  if (!IsOutside(areacode[0])) {
358  if (distance[0] >= 0) isvalid[0] = true;
359  }
360  } else if (validate == kValidateWithoutTol) {
361  areacode[0] = GetAreaCode(xx[0], false);
362  if (IsInside(areacode[0])) {
363  if (distance[0] >= 0) isvalid[0] = true;
364  }
365  } else { // kDontValidate
366  areacode[0] = sInside;
367  if (distance[0] >= 0) isvalid[0] = true;
368  }
369 
370  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
371  isvalid[0], 1, validate, &gp, &gv);
372  return 1;
373  }
374 
375  //
376  // special case end.
377  //
378 
379  G4double a = v.x()*v.x() + v.y()*v.y() - v.z()*v.z()*fTan2Stereo;
380  G4double b = 2.0 * ( p.x() * v.x() + p.y() * v.y() - p.z() * v.z() * fTan2Stereo );
381  G4double c = p.x()*p.x() + p.y()*p.y() - fR02 - p.z()*p.z()*fTan2Stereo;
382  G4double D = b*b - 4*a*c; //discriminant
383  G4int vout = 0;
384 
385  if (std::fabs(a) < DBL_MIN) {
386  if (std::fabs(b) > DBL_MIN) { // single solution
387 
388  distance[0] = -c/b;
389  xx[0] = p + distance[0]*v;
390  gxx[0] = ComputeGlobalPoint(xx[0]);
391 
392  if (validate == kValidateWithTol) {
393  areacode[0] = GetAreaCode(xx[0]);
394  if (!IsOutside(areacode[0])) {
395  if (distance[0] >= 0) isvalid[0] = true;
396  }
397  } else if (validate == kValidateWithoutTol) {
398  areacode[0] = GetAreaCode(xx[0], false);
399  if (IsInside(areacode[0])) {
400  if (distance[0] >= 0) isvalid[0] = true;
401  }
402  } else { // kDontValidate
403  areacode[0] = sInside;
404  if (distance[0] >= 0) isvalid[0] = true;
405  }
406 
407  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
408  isvalid[0], 1, validate, &gp, &gv);
409  vout = 1;
410 
411  } else {
412  // if a=b=0 and c != 0, p is origin and v is parallel to asymptotic line.
413  // if a=b=c=0, p is on surface and v is paralell to stereo wire.
414  // return distance = infinity.
415 
416  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
417  isvalid[0], 0, validate, &gp, &gv);
418 
419  vout = 0;
420  }
421 
422  } else if (D > DBL_MIN) { // double solutions
423 
424  D = std::sqrt(D);
425  G4double factor = 0.5/a;
426  G4double tmpdist[2] = {kInfinity, kInfinity};
427  G4ThreeVector tmpxx[2] ;
428  G4int tmpareacode[2] = {sOutside, sOutside};
429  G4bool tmpisvalid[2] = {false, false};
430  G4int i;
431 
432  for (i=0; i<2; i++) {
433  tmpdist[i] = factor*(-b - D);
434  D = -D;
435  tmpxx[i] = p + tmpdist[i]*v;
436 
437  if (validate == kValidateWithTol) {
438  tmpareacode[i] = GetAreaCode(tmpxx[i]);
439  if (!IsOutside(tmpareacode[i])) {
440  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
441  continue;
442  }
443  } else if (validate == kValidateWithoutTol) {
444  tmpareacode[i] = GetAreaCode(tmpxx[i], false);
445  if (IsInside(tmpareacode[i])) {
446  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
447  continue;
448  }
449  } else { // kDontValidate
450  tmpareacode[i] = sInside;
451  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
452  continue;
453  }
454  }
455 
456  if (tmpdist[0] <= tmpdist[1]) {
457  distance[0] = tmpdist[0];
458  distance[1] = tmpdist[1];
459  xx[0] = tmpxx[0];
460  xx[1] = tmpxx[1];
461  gxx[0] = ComputeGlobalPoint(tmpxx[0]);
462  gxx[1] = ComputeGlobalPoint(tmpxx[1]);
463  areacode[0] = tmpareacode[0];
464  areacode[1] = tmpareacode[1];
465  isvalid[0] = tmpisvalid[0];
466  isvalid[1] = tmpisvalid[1];
467  } else {
468  distance[0] = tmpdist[1];
469  distance[1] = tmpdist[0];
470  xx[0] = tmpxx[1];
471  xx[1] = tmpxx[0];
472  gxx[0] = ComputeGlobalPoint(tmpxx[1]);
473  gxx[1] = ComputeGlobalPoint(tmpxx[0]);
474  areacode[0] = tmpareacode[1];
475  areacode[1] = tmpareacode[0];
476  isvalid[0] = tmpisvalid[1];
477  isvalid[1] = tmpisvalid[0];
478  }
479 
480  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
481  isvalid[0], 2, validate, &gp, &gv);
482  fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
483  isvalid[1], 2, validate, &gp, &gv);
484  vout = 2;
485 
486  } else {
487  // if D<0, no solution
488  // if D=0, just grazing the surfaces, return kInfinity
489 
490  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
491  isvalid[0], 0, validate, &gp, &gv);
492  vout = 0;
493  }
494  return vout;
495 }
496 
497 
498 //=====================================================================
499 //* DistanceToSurface -------------------------------------------------
500 
502  G4ThreeVector gxx[],
503  G4double distance[],
504  G4int areacode[])
505 {
506  // Find the approximate distance of a point of a hyperbolic surface.
507  // The distance must be an underestimate.
508  // It will also be nice (although not necessary) that the estimate is
509  // always finite no matter how close the point is.
510  //
511  // We arranged G4Hype::ApproxDistOutside and G4Hype::ApproxDistInside
512  // for this function. See these discriptions.
513 
514  const G4double halftol
516 
518 
519  if (fCurStat.IsDone()) {
520  for (G4int i=0; i<fCurStat.GetNXX(); i++) {
521  gxx[i] = fCurStat.GetXX(i);
522  distance[i] = fCurStat.GetDistance(i);
523  areacode[i] = fCurStat.GetAreacode(i);
524  }
525  return fCurStat.GetNXX();
526  } else {
527  // initialize
528  for (G4int i=0; i<2; i++) {
529  distance[i] = kInfinity;
530  areacode[i] = sOutside;
531  gxx[i].set(kInfinity, kInfinity, kInfinity);
532  }
533  }
534 
535 
537  G4ThreeVector xx;
538 
539  //
540  // special case!
541  // If p is on surface, return distance = 0 immediatery .
542  //
543  G4ThreeVector lastgxx[2];
544  for (G4int i=0; i<2; i++) {
545  lastgxx[i] = fCurStatWithV.GetXX(i);
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 
649 G4int G4TwistTubsHypeSide::GetAreaCode(const G4ThreeVector &xx,
650  G4bool withTol)
651 {
652  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  std::ostringstream message;
751  message << "Feature NOT implemented !" << G4endl
752  << " fAxis[0] = " << fAxis[0] << G4endl
753  << " fAxis[1] = " << fAxis[1];
754  G4Exception("G4TwistTubsHypeSide::GetAreaCode()",
755  "GeomSolids0001", FatalException, message);
756  }
757  return areacode;
758 }
759 
760 //=====================================================================
761 //* GetAreaCodeInPhi --------------------------------------------------
762 
763 G4int 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 
810 void 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  std::ostringstream message;
862  message << "Feature NOT implemented !" << G4endl
863  << " fAxis[0] = " << fAxis[0] << G4endl
864  << " fAxis[1] = " << fAxis[1];
865  G4Exception("G4TwistTubsHypeSide::SetCorners()",
866  "GeomSolids0001", FatalException, message);
867  }
868 }
869 
870 
871 //=====================================================================
872 //* SetCorners() ------------------------------------------------------
873 
874 void G4TwistTubsHypeSide::SetCorners()
875 {
876  G4Exception("G4TwistTubsHypeSide::SetCorners()",
877  "GeomSolids0001", FatalException,
878  "Method NOT implemented !");
879 }
880 
881 //=====================================================================
882 //* SetBoundaries() ---------------------------------------------------
883 
884 void 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,
898 
899  // sAxis0 & sAxisMax
900  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
901  direction = direction.unit();
902  SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
904 
905  // sAxis1 & sAxisMin
906  direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
907  direction = direction.unit();
908  SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
910 
911  // sAxis1 & sAxisMax
912  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
913  direction = direction.unit();
914  SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
916  } else {
917  std::ostringstream message;
918  message << "Feature NOT implemented !" << G4endl
919  << " fAxis[0] = " << fAxis[0] << G4endl
920  << " fAxis[1] = " << fAxis[1];
921  G4Exception("G4TwistTubsHypeSide::SetBoundaries()",
922  "GeomSolids0001", FatalException, message);
923  }
924 }
925 
926 //=====================================================================
927 //* GetFacets() -------------------------------------------------------
928 
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)*(k-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<k ; j++ )
950  {
951  nnode = GetNode(i,j,k,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)/(k-1) ;
958  } else { // outer hyperbolic surface
959  x = xmax - j*(xmax-xmin)/(k-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<k-1 ) { // clock wise filling
969 
970  nface = GetFace(i,j,k,n,iside) ;
971 
972  faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i ,j ,k,n,iside)+1) ;
973  faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
974  faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
975  faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
976 
977  }
978  }
979  }
980 }
static const G4int sAxisZ
void set(double x, double y, double z)
Definition: geomdefs.hh:54
G4int GetAreacode(G4int i) const
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
static const G4int sC0Min1Max
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
static const G4int sAxisPhi
virtual G4double GetBoundaryMin(G4double phi)
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
double x() const
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
G4SurfCurNormal fCurrentNormal
G4double z
Definition: TRTMaterials.hh:39
const char * p
Definition: xmltok.h:285
G4double fAxisMax[2]
G4ThreeVector GetCorner(G4int areacode) const
static const G4int sOutside
virtual G4ThreeVector GetBoundaryAtPZ(G4int areacode, const G4ThreeVector &p) const
const XML_Char * name
G4bool IsBoundary(G4int areacode, G4bool testbitmode=false) const
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
double getRho() const
virtual G4double GetBoundaryMax(G4double phi)
CurrentStatus fCurStat
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
int G4int
Definition: G4Types.hh:78
double z() const
G4ThreeVector fTrans
G4double fAxisMin[2]
static const G4int sC0Min1Min
G4GLOB_DLL std::ostream G4cout
G4bool IsOutside(G4int areacode) const
virtual EInside Inside(const G4ThreeVector &gp)
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
G4double GetDistance(G4int i) const
static const G4int sAxis1
static const G4int sC0Max1Max
static const G4int sBoundary
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
G4double DistanceToLine(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &d, G4ThreeVector &xx)
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
G4bool IsValid(G4int i) const
static const G4int sAxis0
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
static const G4int sAxisMin
const G4int n
static const G4int sInside
virtual G4double GetRhoAtPZ(const G4ThreeVector &p, G4bool isglobal=false) const
G4TwistTubsHypeSide(const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, const G4int handedness, const G4double kappa, const G4double tanstereo, const G4double r0, const EAxis axis0=kPhi, const EAxis axis1=kZAxis, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const G4int sCorner
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
EInside
Definition: geomdefs.hh:58
static const G4int sAxisMax
EAxis
Definition: geomdefs.hh:54
Hep3Vector unit() const
double y() const
G4ThreeVector GetXX(G4int i) const
virtual G4String GetName() const
#define DBL_MIN
Definition: templates.hh:75
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
virtual G4int AmIOnLeftSide(const G4ThreeVector &me, const G4ThreeVector &vec, G4bool withTol=true)
const XML_Char int len
#define G4endl
Definition: G4ios.hh:61
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
double G4double
Definition: G4Types.hh:76
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
double mag() const
static const G4int sC0Max1Min
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)
CurrentStatus fCurStatWithV
static G4GeometryTolerance * GetInstance()
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)