Geant4  10.01.p03
G4TwistTrapFlatSide.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: G4TwistTrapFlatSide.cc 66356 2012-12-18 09:02:32Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 //
34 // G4TwistTrapFlatSide.cc
35 //
36 // Author:
37 // 30-Aug-2002 - O.Link (Oliver.Link@cern.ch)
38 //
39 // --------------------------------------------------------------------
40 
41 #include "G4TwistTrapFlatSide.hh"
42 
43 //=====================================================================
44 //* constructors ------------------------------------------------------
45 
47  G4double PhiTwist,
48  G4double pDx1,
49  G4double pDx2,
50  G4double pDy,
51  G4double pDz,
52  G4double pAlpha,
53  G4double pPhi,
54  G4double pTheta,
55  G4int handedness)
56 
57  : G4VTwistSurface(name)
58 {
59  fHandedness = handedness; // +z = +ve, -z = -ve
60 
61  fDx1 = pDx1 ;
62  fDx2 = pDx2 ;
63  fDy = pDy ;
64  fDz = pDz ;
65  fAlpha = pAlpha ;
66  fTAlph = std::tan(fAlpha) ;
67  fPhi = pPhi ;
68  fTheta = pTheta ;
69 
70  fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
71  // dx in surface equation
72  fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
73  // dy in surface equation
74 
75  fPhiTwist = PhiTwist ;
76 
77  fCurrentNormal.normal.set( 0, 0, (fHandedness < 0 ? -1 : 1));
78  // Unit vector, in local coordinate system
79  fRot.rotateZ( fHandedness > 0
80  ? 0.5 * fPhiTwist
81  : -0.5 * fPhiTwist );
82 
83  fTrans.set(
84  fHandedness > 0 ? 0.5*fdeltaX : -0.5*fdeltaX ,
85  fHandedness > 0 ? 0.5*fdeltaY : -0.5*fdeltaY ,
86  fHandedness > 0 ? fDz : -fDz ) ;
87 
88  fIsValidNorm = true;
89 
90 
91  fAxis[0] = kXAxis ;
92  fAxis[1] = kYAxis ;
93  fAxisMin[0] = kInfinity ; // x-Axis cannot be fixed, because it
94  fAxisMax[0] = kInfinity ; // depends on y
95  fAxisMin[1] = -fDy ; // y - axis
96  fAxisMax[1] = fDy ;
97 
98  SetCorners();
99  SetBoundaries();
100 }
101 
102 
103 //=====================================================================
104 //* Fake default constructor ------------------------------------------
105 
107  : G4VTwistSurface(a), fDx1(0.), fDx2(0.), fDy(0.), fDz(0.), fPhiTwist(0.),
108  fAlpha(0.), fTAlph(0.), fPhi(0.), fTheta(0.), fdeltaX(0.), fdeltaY(0.)
109 {
110 }
111 
112 
113 //=====================================================================
114 //* destructor --------------------------------------------------------
115 
117 {
118 }
119 
120 //=====================================================================
121 //* GetNormal ---------------------------------------------------------
122 
124  G4bool isGlobal)
125 {
126  if (isGlobal) {
128  } else {
129  return fCurrentNormal.normal;
130  }
131 }
132 
133 //=====================================================================
134 //* DistanceToSurface(p, v) -------------------------------------------
135 
137  const G4ThreeVector &gv,
138  G4ThreeVector gxx[],
139  G4double distance[],
140  G4int areacode[],
141  G4bool isvalid[],
142  EValidate validate)
143 {
144  fCurStatWithV.ResetfDone(validate, &gp, &gv);
145 
146  if (fCurStatWithV.IsDone()) {
147  G4int i;
148  for (i=0; i<fCurStatWithV.GetNXX(); i++) {
149  gxx[i] = fCurStatWithV.GetXX(i);
150  distance[i] = fCurStatWithV.GetDistance(i);
151  areacode[i] = fCurStatWithV.GetAreacode(i);
152  isvalid[i] = fCurStatWithV.IsValid(i);
153  }
154  return fCurStatWithV.GetNXX();
155  } else {
156  // initialize
157  G4int i;
158  for (i=0; i<2; i++) {
159  distance[i] = kInfinity;
160  areacode[i] = sOutside;
161  isvalid[i] = false;
162  gxx[i].set(kInfinity, kInfinity, kInfinity);
163  }
164  }
165 
168 
169  //
170  // special case!
171  // if p is on surface, distance = 0.
172  //
173 
174  if (std::fabs(p.z()) == 0.) { // if p is on the plane
175  distance[0] = 0;
176  G4ThreeVector xx = p;
177  gxx[0] = ComputeGlobalPoint(xx);
178 
179  if (validate == kValidateWithTol) {
180  areacode[0] = GetAreaCode(xx);
181  if (!IsOutside(areacode[0])) {
182  isvalid[0] = true;
183  }
184  } else if (validate == kValidateWithoutTol) {
185  areacode[0] = GetAreaCode(xx, false);
186  if (IsInside(areacode[0])) {
187  isvalid[0] = true;
188  }
189  } else { // kDontValidate
190  areacode[0] = sInside;
191  isvalid[0] = true;
192  }
193 
194  return 1;
195  }
196  //
197  // special case end
198  //
199 
200  if (v.z() == 0) {
201 
202  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
203  isvalid[0], 0, validate, &gp, &gv);
204  return 0;
205  }
206 
207  distance[0] = - (p.z() / v.z());
208 
209  G4ThreeVector xx = p + distance[0]*v;
210  gxx[0] = ComputeGlobalPoint(xx);
211 
212  if (validate == kValidateWithTol) {
213  areacode[0] = GetAreaCode(xx);
214  if (!IsOutside(areacode[0])) {
215  if (distance[0] >= 0) isvalid[0] = true;
216  }
217  } else if (validate == kValidateWithoutTol) {
218  areacode[0] = GetAreaCode(xx, false);
219  if (IsInside(areacode[0])) {
220  if (distance[0] >= 0) isvalid[0] = true;
221  }
222  } else { // kDontValidate
223  areacode[0] = sInside;
224  if (distance[0] >= 0) isvalid[0] = true;
225  }
226 
227 
228  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
229  isvalid[0], 1, validate, &gp, &gv);
230 
231 #ifdef G4TWISTDEBUG
232  G4cerr << "ERROR - G4TwistTrapFlatSide::DistanceToSurface(p,v)" << G4endl;
233  G4cerr << " Name : " << GetName() << G4endl;
234  G4cerr << " xx : " << xx << G4endl;
235  G4cerr << " gxx[0] : " << gxx[0] << G4endl;
236  G4cerr << " dist[0] : " << distance[0] << G4endl;
237  G4cerr << " areacode[0] : " << areacode[0] << G4endl;
238  G4cerr << " isvalid[0] : " << isvalid[0] << G4endl;
239 #endif
240  return 1;
241 }
242 
243 //=====================================================================
244 //* DistanceToSurface(p) ----------------------------------------------
245 
247  G4ThreeVector gxx[],
248  G4double distance[],
249  G4int areacode[])
250 {
251  // Calculate distance to plane in local coordinate,
252  // then return distance and global intersection points.
253  //
254 
256 
257  if (fCurStat.IsDone()) {
258  G4int i;
259  for (i=0; i<fCurStat.GetNXX(); i++) {
260  gxx[i] = fCurStat.GetXX(i);
261  distance[i] = fCurStat.GetDistance(i);
262  areacode[i] = fCurStat.GetAreacode(i);
263  }
264  return fCurStat.GetNXX();
265  } else {
266  // initialize
267  G4int i;
268  for (i=0; i<2; i++) {
269  distance[i] = kInfinity;
270  areacode[i] = sOutside;
271  gxx[i].set(kInfinity, kInfinity, kInfinity);
272  }
273  }
274 
276  G4ThreeVector xx;
277 
278  // The plane is placed on origin with making its normal
279  // parallel to z-axis.
280  if (std::fabs(p.z()) <= 0.5 * kCarTolerance)
281  { // if p is on the plane, return 1
282  distance[0] = 0;
283  xx = p;
284  } else {
285  distance[0] = std::fabs(p.z());
286  xx.set(p.x(), p.y(), 0);
287  }
288 
289  gxx[0] = ComputeGlobalPoint(xx);
290  areacode[0] = sInside;
291  G4bool isvalid = true;
292  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
293  isvalid, 1, kDontValidate, &gp);
294  return 1;
295 
296 }
297 
299  G4bool withTol)
300 {
301 
302  static const G4double ctol = 0.5 * kCarTolerance;
303  G4int areacode = sInside;
304 
305  if (fAxis[0] == kXAxis && fAxis[1] == kYAxis) {
306 
307  G4int yaxis = 1;
308 
309  G4double wmax = xAxisMax(xx.y(), fTAlph ) ;
310  G4double wmin = -xAxisMax(xx.y(), -fTAlph ) ;
311 
312  if (withTol) {
313 
314  G4bool isoutside = false;
315 
316  // test boundary of x-axis
317 
318  if (xx.x() < wmin + ctol) {
319  areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
320  if (xx.x() <= wmin - ctol) isoutside = true;
321 
322  } else if (xx.x() > wmax - ctol) {
323  areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
324  if (xx.x() >= wmax + ctol) isoutside = true;
325  }
326 
327  // test boundary of y-axis
328 
329  if (xx.y() < fAxisMin[yaxis] + ctol) {
330  areacode |= (sAxis1 & (sAxisY | sAxisMin));
331 
332  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
333  else areacode |= sBoundary;
334  if (xx.y() <= fAxisMin[yaxis] - ctol) isoutside = true;
335 
336  } else if (xx.y() > fAxisMax[yaxis] - ctol) {
337  areacode |= (sAxis1 & (sAxisY | sAxisMax));
338 
339  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
340  else areacode |= sBoundary;
341  if (xx.y() >= fAxisMax[yaxis] + ctol) isoutside = true;
342  }
343 
344  // if isoutside = true, clear inside bit.
345  // if not on boundary, add axis information.
346 
347  if (isoutside) {
348  G4int tmpareacode = areacode & (~sInside);
349  areacode = tmpareacode;
350  } else if ((areacode & sBoundary) != sBoundary) {
351  areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisY);
352  }
353 
354  } else {
355 
356  // boundary of x-axis
357 
358  if (xx.x() < wmin ) {
359  areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
360  } else if (xx.x() > wmax) {
361  areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
362  }
363 
364  // boundary of y-axis
365 
366  if (xx.y() < fAxisMin[yaxis]) {
367  areacode |= (sAxis1 & (sAxisY | sAxisMin));
368  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
369  else areacode |= sBoundary;
370 
371  } else if (xx.y() > fAxisMax[yaxis]) {
372  areacode |= (sAxis1 & (sAxisY | sAxisMax)) ;
373  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
374  else areacode |= sBoundary;
375  }
376 
377  if ((areacode & sBoundary) != sBoundary) {
378  areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisY);
379  }
380  }
381  return areacode;
382  } else {
383  G4Exception("G4TwistTrapFlatSide::GetAreaCode()",
384  "GeomSolids0001", FatalException,
385  "Feature NOT implemented !");
386  }
387 
388  return areacode;
389 }
390 
391 
392 //=====================================================================
393 //* SetCorners --------------------------------------------------------
394 
396 {
397  // Set Corner points in local coodinate.
398 
399  if (fAxis[0] == kXAxis && fAxis[1] == kYAxis) {
400 
401  G4double x, y, z;
402 
403  // corner of Axis0min and Axis1min
404  x = -fDx1 + fDy * fTAlph ;
405  y = -fDy ;
406  z = 0 ;
407  SetCorner(sC0Min1Min, x, y, z);
408 
409  // corner of Axis0max and Axis1min
410  x = fDx1 + fDy * fTAlph ;
411  y = -fDy ;
412  z = 0 ;
413  SetCorner(sC0Max1Min, x, y, z);
414 
415  // corner of Axis0max and Axis1max
416  x = fDx2 - fDy * fTAlph ;
417  y = fDy ;
418  z = 0 ;
419  SetCorner(sC0Max1Max, x, y, z);
420 
421  // corner of Axis0min and Axis1max
422  x = -fDx2 - fDy * fTAlph ;
423  y = fDy ;
424  z = 0 ;
425  SetCorner(sC0Min1Max, x, y, z);
426 
427  } else {
428  std::ostringstream message;
429  message << "Feature NOT implemented !" << G4endl
430  << " fAxis[0] = " << fAxis[0] << G4endl
431  << " fAxis[1] = " << fAxis[1];
432  G4Exception("G4TwistTrapFlatSide::SetCorners()",
433  "GeomSolids0001", FatalException, message);
434  }
435 }
436 
437 //=====================================================================
438 //* SetBoundaries() ---------------------------------------------------
439 
441 {
442  // Set direction-unit vector of phi-boundary-lines in local coodinate.
443  // Don't call the function twice.
444 
445  G4ThreeVector direction ;
446 
447  if (fAxis[0] == kXAxis && fAxis[1] == kYAxis) {
448 
449  // sAxis0 & sAxisMin
450  direction = - ( GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min) ) ;
451  direction = direction.unit();
452  SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction,
454 
455  // sAxis0 & sAxisMax
456  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min) ; // inverse
457  direction = direction.unit();
458  SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction,
460 
461  // sAxis1 & sAxisMin
462  direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
463  direction = direction.unit();
464  SetBoundary(sAxis1 & (sAxisY | sAxisMin), direction,
466 
467  // sAxis1 & sAxisMax
468  direction = - ( GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max) ) ;
469  direction = direction.unit();
470  SetBoundary(sAxis1 & (sAxisY | sAxisMax), direction,
472 
473  } else {
474  std::ostringstream message;
475  message << "Feature NOT implemented !" << G4endl
476  << " fAxis[0] = " << fAxis[0] << G4endl
477  << " fAxis[1] = " << fAxis[1];
478  G4Exception("G4TwistTrapFlatSide::SetCorners()",
479  "GeomSolids0001", FatalException, message);
480  }
481 }
482 
483 //=====================================================================
484 //* GetFacets() -------------------------------------------------------
485 
487  G4int faces[][4], G4int iside )
488 {
489 
490  G4double x,y ; // the two parameters for the surface equation
491  G4ThreeVector p ; // a point on the surface, given by (z,u)
492 
493  G4int nnode ;
494  G4int nface ;
495 
496  G4double xmin,xmax ;
497 
498  // calculate the (n-1)*(k-1) vertices
499 
500  G4int i,j ;
501 
502  for ( i = 0 ; i<n ; i++ ) {
503 
504  y = -fDy + i*(2*fDy)/(n-1) ;
505 
506  for ( j = 0 ; j<k ; j++ ) {
507 
508  xmin = GetBoundaryMin(y) ;
509  xmax = GetBoundaryMax(y) ;
510  x = xmin + j*(xmax-xmin)/(k-1) ;
511 
512  nnode = GetNode(i,j,k,n,iside) ;
513  p = SurfacePoint(x,y,true) ; // surface point in global coordinate system
514 
515  xyz[nnode][0] = p.x() ;
516  xyz[nnode][1] = p.y() ;
517  xyz[nnode][2] = p.z() ;
518 
519  if ( i<n-1 && j<k-1 ) {
520 
521  nface = GetFace(i,j,k,n,iside) ;
522 
523  if (fHandedness < 0) { // lower side
524  faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i ,j ,k,n,iside)+1) ;
525  faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
526  faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
527  faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
528  } else { // upper side
529  faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * ( GetNode(i ,j ,k,n,iside)+1) ;
530  faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
531  faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
532  faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
533  }
534 
535  }
536  }
537  }
538 }
virtual G4ThreeVector SurfacePoint(G4double x, G4double y, G4bool isGlobal=false)
G4int GetAreacode(G4int i) const
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
static const G4int sC0Min1Max
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
G4SurfCurNormal fCurrentNormal
G4double z
Definition: TRTMaterials.hh:39
G4String name
Definition: TRTMaterials.hh:40
G4double fAxisMax[2]
G4ThreeVector GetCorner(G4int areacode) const
static const G4int sOutside
G4double a
Definition: TRTMaterials.hh:39
CurrentStatus fCurStat
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
int G4int
Definition: G4Types.hh:78
G4RotationMatrix fRot
static const G4int sAxisX
G4ThreeVector fTrans
G4double fAxisMin[2]
virtual G4ThreeVector GetNormal(const G4ThreeVector &, G4bool isGlobal=false)
static const G4int sC0Min1Min
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
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
static const G4int sAxisY
const G4int n
static const G4int sInside
virtual G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)
G4TwistTrapFlatSide(const G4String &name, G4double PhiTwist, G4double pDx1, G4double pDx2, G4double pDy, G4double pDz, G4double pAlpha, G4double pPhi, G4double pTheta, G4int handedness)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const G4int sCorner
virtual G4double GetBoundaryMin(G4double u)
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sAxisMax
virtual G4double GetBoundaryMax(G4double u)
G4ThreeVector GetXX(G4int i) const
virtual G4String GetName() const
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
virtual void SetBoundaries()
#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
CurrentStatus fCurStatWithV
G4GLOB_DLL std::ostream G4cerr
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
double xAxisMax(G4double u, G4double fTanAlpha) const