Geant4  10.02
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
const G4double x[NPOINTSGL]
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