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