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