| 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 // G4VTwistedFaceted implementation
27 //
28 // Author: 04-Nov-2004 - O.Link (Oliver.Link@cern.ch)
29 // --------------------------------------------------------------------
30
31 #include "G4VTwistedFaceted.hh"
32
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4VoxelLimits.hh"
36 #include "G4AffineTransform.hh"
37 #include "G4BoundingEnvelope.hh"
38 #include "G4SolidExtentList.hh"
39 #include "G4ClippablePolygon.hh"
40 #include "G4VPVParameterisation.hh"
41 #include "G4GeometryTolerance.hh"
42 #include "meshdefs.hh"
43
44 #include "G4VGraphicsScene.hh"
45 #include "G4Polyhedron.hh"
46 #include "G4VisExtent.hh"
47
48 #include "Randomize.hh"
49
50 #include "G4AutoLock.hh"
51
52 namespace
53 {
54 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
55 }
56
57
58 //=====================================================================
59 //* constructors ------------------------------------------------------
60
61 G4VTwistedFaceted::
62 G4VTwistedFaceted( const G4String& pname, // Name of instance
63 G4double PhiTwist, // twist angle
64 G4double pDz, // half z length
65 G4double pTheta, // direction between end planes
66 G4double pPhi, // defined by polar and azim. angles
67 G4double pDy1, // half y length at -pDz
68 G4double pDx1, // half x length at -pDz,-pDy
69 G4double pDx2, // half x length at -pDz,+pDy
70 G4double pDy2, // half y length at +pDz
71 G4double pDx3, // half x length at +pDz,-pDy
72 G4double pDx4, // half x length at +pDz,+pDy
73 G4double pAlph ) // tilt angle
74 : G4VSolid(pname),
75 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fSide0(nullptr),
76 fSide90(nullptr), fSide180(nullptr), fSide270(nullptr)
77 {
78
79 G4double pDytmp ;
80 G4double fDxUp ;
81 G4double fDxDown ;
82
83 fDx1 = pDx1 ;
84 fDx2 = pDx2 ;
85 fDx3 = pDx3 ;
86 fDx4 = pDx4 ;
87 fDy1 = pDy1 ;
88 fDy2 = pDy2 ;
89 fDz = pDz ;
90
91 G4double kAngTolerance
92 = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
93
94 // maximum values
95 //
96 fDxDown = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ;
97 fDxUp = ( fDx3 > fDx4 ? fDx3 : fDx4 ) ;
98 fDx = ( fDxUp > fDxDown ? fDxUp : fDxDown ) ;
99 fDy = ( fDy1 > fDy2 ? fDy1 : fDy2 ) ;
100
101 // planarity check
102 //
103 if ( fDx1 != fDx2 && fDx3 != fDx4 )
104 {
105 pDytmp = fDy1 * ( fDx3 - fDx4 ) / ( fDx1 - fDx2 ) ;
106 if ( std::fabs(pDytmp - fDy2) > kCarTolerance )
107 {
108 std::ostringstream message;
109 message << "Not planar surface in untwisted Trapezoid: "
110 << GetName() << G4endl
111 << "fDy2 is " << fDy2 << " but should be "
112 << pDytmp << ".";
113 G4Exception("G4VTwistedFaceted::G4VTwistedFaceted()", "GeomSolids0002",
114 FatalErrorInArgument, message);
115 }
116 }
117
118 #ifdef G4TWISTDEBUG
119 if ( fDx1 == fDx2 && fDx3 == fDx4 )
120 {
121 G4cout << "Trapezoid is a box" << G4endl ;
122 }
123
124 #endif
125
126 if ( ( fDx1 == fDx2 && fDx3 != fDx4 ) || ( fDx1 != fDx2 && fDx3 == fDx4 ) )
127 {
128 std::ostringstream message;
129 message << "Not planar surface in untwisted Trapezoid: "
130 << GetName() << G4endl
131 << "One endcap is rectangular, the other is a trapezoid." << G4endl
132 << "For planarity reasons they have to be rectangles or trapezoids "
133 << "on both sides.";
134 G4Exception("G4VTwistedFaceted::G4VTwistedFaceted()", "GeomSolids0002",
135 FatalErrorInArgument, message);
136 }
137
138 // twist angle
139 //
140 fPhiTwist = PhiTwist ;
141
142 // tilt angle
143 //
144 fAlph = pAlph ;
145 fTAlph = std::tan(fAlph) ;
146
147 fTheta = pTheta ;
148 fPhi = pPhi ;
149
150 // dx in surface equation
151 //
152 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
153
154 // dy in surface equation
155 //
156 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
157
158 if ( ( fDx1 <= 2*kCarTolerance)
159 || ( fDx2 <= 2*kCarTolerance)
160 || ( fDx3 <= 2*kCarTolerance)
161 || ( fDx4 <= 2*kCarTolerance)
162 || ( fDy1 <= 2*kCarTolerance)
163 || ( fDy2 <= 2*kCarTolerance)
164 || ( fDz <= 2*kCarTolerance)
165 || ( std::fabs(fPhiTwist) <= 2*kAngTolerance )
166 || ( std::fabs(fPhiTwist) >= pi/2 )
167 || ( std::fabs(fAlph) >= pi/2 )
168 || fTheta >= pi/2 || fTheta < 0
169 )
170 {
171 std::ostringstream message;
172 message << "Invalid dimensions. Too small, or twist angle too big: "
173 << GetName() << G4endl
174 << "fDx 1-4 = " << fDx1/cm << ", " << fDx2/cm << ", "
175 << fDx3/cm << ", " << fDx4/cm << " cm" << G4endl
176 << "fDy 1-2 = " << fDy1/cm << ", " << fDy2/cm << ", "
177 << " cm" << G4endl
178 << "fDz = " << fDz/cm << " cm" << G4endl
179 << " twistangle " << fPhiTwist/deg << " deg" << G4endl
180 << " phi,theta = " << fPhi/deg << ", " << fTheta/deg << " deg";
181 G4Exception("G4TwistedTrap::G4VTwistedFaceted()",
182 "GeomSolids0002", FatalErrorInArgument, message);
183 }
184 CreateSurfaces();
185 }
186
187
188 //=====================================================================
189 //* Fake default constructor ------------------------------------------
190
191 G4VTwistedFaceted::G4VTwistedFaceted( __void__& a )
192 : G4VSolid(a),
193 fLowerEndcap(nullptr), fUpperEndcap(nullptr),
194 fSide0(nullptr), fSide90(nullptr), fSide180(nullptr), fSide270(nullptr)
195 {
196 }
197
198
199 //=====================================================================
200 //* destructor --------------------------------------------------------
201
202 G4VTwistedFaceted::~G4VTwistedFaceted()
203 {
204 delete fLowerEndcap ;
205 delete fUpperEndcap ;
206
207 delete fSide0 ;
208 delete fSide90 ;
209 delete fSide180 ;
210 delete fSide270 ;
211 delete fpPolyhedron; fpPolyhedron = nullptr;
212 }
213
214
215 //=====================================================================
216 //* Copy constructor --------------------------------------------------
217
218 G4VTwistedFaceted::G4VTwistedFaceted(const G4VTwistedFaceted& rhs)
219 : G4VSolid(rhs),
220 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
221 fTheta(rhs.fTheta), fPhi(rhs.fPhi),
222 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fDy2(rhs.fDy2),
223 fDx3(rhs.fDx3), fDx4(rhs.fDx4), fDz(rhs.fDz), fDx(rhs.fDx), fDy(rhs.fDy),
224 fAlph(rhs.fAlph), fTAlph(rhs.fTAlph), fdeltaX(rhs.fdeltaX),
225 fdeltaY(rhs.fdeltaY), fPhiTwist(rhs.fPhiTwist), fLowerEndcap(nullptr),
226 fUpperEndcap(nullptr), fSide0(nullptr), fSide90(nullptr), fSide180(nullptr), fSide270(nullptr)
227 {
228 CreateSurfaces();
229 }
230
231
232 //=====================================================================
233 //* Assignment operator -----------------------------------------------
234
235 G4VTwistedFaceted& G4VTwistedFaceted::operator = (const G4VTwistedFaceted& rhs)
236 {
237 // Check assignment to self
238 //
239 if (this == &rhs) { return *this; }
240
241 // Copy base class data
242 //
243 G4VSolid::operator=(rhs);
244
245 // Copy data
246 //
247 fTheta = rhs.fTheta; fPhi = rhs.fPhi;
248 fDy1= rhs.fDy1; fDx1= rhs.fDx1; fDx2= rhs.fDx2; fDy2= rhs.fDy2;
249 fDx3= rhs.fDx3; fDx4= rhs.fDx4; fDz= rhs.fDz; fDx= rhs.fDx; fDy= rhs.fDy;
250 fAlph= rhs.fAlph; fTAlph= rhs.fTAlph; fdeltaX= rhs.fdeltaX;
251 fdeltaY= rhs.fdeltaY; fPhiTwist= rhs.fPhiTwist; fLowerEndcap= nullptr;
252 fUpperEndcap= nullptr; fSide0= nullptr; fSide90= nullptr; fSide180= nullptr; fSide270= nullptr;
253 fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
254 fRebuildPolyhedron = false;
255 delete fpPolyhedron; fpPolyhedron = nullptr;
256
257 CreateSurfaces();
258
259 return *this;
260 }
261
262
263 //=====================================================================
264 //* ComputeDimensions -------------------------------------------------
265
266 void G4VTwistedFaceted::ComputeDimensions(G4VPVParameterisation* ,
267 const G4int ,
268 const G4VPhysicalVolume* )
269 {
270 G4Exception("G4VTwistedFaceted::ComputeDimensions()",
271 "GeomSolids0001", FatalException,
272 "G4VTwistedFaceted does not support Parameterisation.");
273 }
274
275
276 //=====================================================================
277 //* Extent ------------------------------------------------------------
278
279 void G4VTwistedFaceted::BoundingLimits(G4ThreeVector& pMin,
280 G4ThreeVector& pMax) const
281 {
282 G4double cosPhi = std::cos(fPhi);
283 G4double sinPhi = std::sin(fPhi);
284 G4double tanTheta = std::tan(fTheta);
285 G4double tanAlpha = fTAlph;
286
287 G4double xmid1 = fDy1*tanAlpha;
288 G4double x1 = std::abs(xmid1 + fDx1);
289 G4double x2 = std::abs(xmid1 - fDx1);
290 G4double x3 = std::abs(xmid1 + fDx2);
291 G4double x4 = std::abs(xmid1 - fDx2);
292 G4double xmax1 = std::max(std::max(std::max(x1, x2), x3), x4);
293 G4double rmax1 = std::sqrt(xmax1*xmax1 + fDy1*fDy1);
294
295 G4double xmid2 = fDy2*tanAlpha;
296 G4double x5 = std::abs(xmid2 + fDx3);
297 G4double x6 = std::abs(xmid2 - fDx3);
298 G4double x7 = std::abs(xmid2 + fDx4);
299 G4double x8 = std::abs(xmid2 - fDx4);
300 G4double xmax2 = std::max(std::max(std::max(x5, x6), x7), x8);
301 G4double rmax2 = std::sqrt(xmax2*xmax2 + fDy2*fDy2);
302
303 G4double x0 = fDz*tanTheta*cosPhi;
304 G4double y0 = fDz*tanTheta*sinPhi;
305 G4double xmin = std::min(-x0 - rmax1, x0 - rmax2);
306 G4double ymin = std::min(-y0 - rmax1, y0 - rmax2);
307 G4double xmax = std::max(-x0 + rmax1, x0 + rmax2);
308 G4double ymax = std::max(-y0 + rmax1, y0 + rmax2);
309 pMin.set(xmin, ymin,-fDz);
310 pMax.set(xmax, ymax, fDz);
311 }
312
313
314 //=====================================================================
315 //* CalculateExtent ---------------------------------------------------
316
317 G4bool
318 G4VTwistedFaceted::CalculateExtent( const EAxis pAxis,
319 const G4VoxelLimits& pVoxelLimit,
320 const G4AffineTransform& pTransform,
321 G4double& pMin,
322 G4double& pMax ) const
323 {
324 G4ThreeVector bmin, bmax;
325
326 // Get bounding box
327 BoundingLimits(bmin,bmax);
328
329 // Find extent
330 G4BoundingEnvelope bbox(bmin,bmax);
331 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
332 }
333
334
335 //=====================================================================
336 //* Inside ------------------------------------------------------------
337
338 EInside G4VTwistedFaceted::Inside(const G4ThreeVector& p) const
339 {
340
341 EInside tmpin = kOutside ;
342
343 G4double phi = p.z()/(2*fDz) * fPhiTwist ; // rotate the point to z=0
344 G4double cphi = std::cos(-phi) ;
345 G4double sphi = std::sin(-phi) ;
346
347 G4double px = p.x() + fdeltaX * ( -phi/fPhiTwist) ; // shift
348 G4double py = p.y() + fdeltaY * ( -phi/fPhiTwist) ;
349 G4double pz = p.z() ;
350
351 G4double posx = px * cphi - py * sphi ; // rotation
352 G4double posy = px * sphi + py * cphi ;
353 G4double posz = pz ;
354
355 G4double xMin = Xcoef(posy,phi,fTAlph) - 2*Xcoef(posy,phi,0.) ;
356 G4double xMax = Xcoef(posy,phi,fTAlph) ;
357
358 G4double yMax = GetValueB(phi)/2. ; // b(phi)/2 is limit
359 G4double yMin = -yMax ;
360
361 #ifdef G4TWISTDEBUG
362
363 G4cout << "inside called: p = " << p << G4endl ;
364 G4cout << "fDx1 = " << fDx1 << G4endl ;
365 G4cout << "fDx2 = " << fDx2 << G4endl ;
366 G4cout << "fDx3 = " << fDx3 << G4endl ;
367 G4cout << "fDx4 = " << fDx4 << G4endl ;
368
369 G4cout << "fDy1 = " << fDy1 << G4endl ;
370 G4cout << "fDy2 = " << fDy2 << G4endl ;
371
372 G4cout << "fDz = " << fDz << G4endl ;
373
374 G4cout << "Tilt angle alpha = " << fAlph << G4endl ;
375 G4cout << "phi,theta = " << fPhi << " , " << fTheta << G4endl ;
376
377 G4cout << "Twist angle = " << fPhiTwist << G4endl ;
378
379 G4cout << "posx = " << posx << G4endl ;
380 G4cout << "posy = " << posy << G4endl ;
381 G4cout << "xMin = " << xMin << G4endl ;
382 G4cout << "xMax = " << xMax << G4endl ;
383 G4cout << "yMin = " << yMin << G4endl ;
384 G4cout << "yMax = " << yMax << G4endl ;
385
386 #endif
387
388
389 if ( posx <= xMax - kCarTolerance*0.5
390 && posx >= xMin + kCarTolerance*0.5 )
391 {
392 if ( posy <= yMax - kCarTolerance*0.5
393 && posy >= yMin + kCarTolerance*0.5 )
394 {
395 if (std::fabs(posz) <= fDz - kCarTolerance*0.5 ) tmpin = kInside ;
396 else if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) tmpin = kSurface ;
397 }
398 else if ( posy <= yMax + kCarTolerance*0.5
399 && posy >= yMin - kCarTolerance*0.5 )
400 {
401 if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) tmpin = kSurface ;
402 }
403 }
404 else if ( posx <= xMax + kCarTolerance*0.5
405 && posx >= xMin - kCarTolerance*0.5 )
406 {
407 if ( posy <= yMax + kCarTolerance*0.5
408 && posy >= yMin - kCarTolerance*0.5 )
409 {
410 if (std::fabs(posz) <= fDz + kCarTolerance*0.5) tmpin = kSurface ;
411 }
412 }
413
414 #ifdef G4TWISTDEBUG
415 G4cout << "inside = " << tmpin << G4endl ;
416 #endif
417
418 return tmpin;
419
420 }
421
422
423 //=====================================================================
424 //* SurfaceNormal -----------------------------------------------------
425
426 G4ThreeVector G4VTwistedFaceted::SurfaceNormal(const G4ThreeVector& p) const
427 {
428 //
429 // return the normal unit vector to the Hyperbolical Surface at a point
430 // p on (or nearly on) the surface
431 //
432 // Which of the three or four surfaces are we closest to?
433 //
434
435
436 G4double distance = kInfinity;
437
438 G4VTwistSurface* surfaces[6];
439
440 surfaces[0] = fSide0 ;
441 surfaces[1] = fSide90 ;
442 surfaces[2] = fSide180 ;
443 surfaces[3] = fSide270 ;
444 surfaces[4] = fLowerEndcap;
445 surfaces[5] = fUpperEndcap;
446
447 G4ThreeVector xx;
448 G4ThreeVector bestxx;
449 G4int i;
450 G4int besti = -1;
451 for (i=0; i< 6; i++)
452 {
453 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
454 if (tmpdistance < distance)
455 {
456 distance = tmpdistance;
457 bestxx = xx;
458 besti = i;
459 }
460 }
461
462 return surfaces[besti]->GetNormal(bestxx, true);
463 }
464
465
466 //=====================================================================
467 //* DistanceToIn (p, v) -----------------------------------------------
468
469 G4double G4VTwistedFaceted::DistanceToIn (const G4ThreeVector& p,
470 const G4ThreeVector& v ) const
471 {
472
473 // DistanceToIn (p, v):
474 // Calculate distance to surface of shape from `outside'
475 // along with the v, allowing for tolerance.
476 // The function returns kInfinity if no intersection or
477 // just grazing within tolerance.
478
479 //
480 // Calculate DistanceToIn(p,v)
481 //
482
483 EInside currentside = Inside(p);
484
485 if (currentside == kInside)
486 {
487 }
488 else if (currentside == kSurface)
489 {
490 // particle is just on a boundary.
491 // if the particle is entering to the volume, return 0
492 //
493 G4ThreeVector normal = SurfaceNormal(p);
494 if (normal*v < 0)
495 {
496 return 0;
497 }
498 }
499
500 // now, we can take smallest positive distance.
501
502 // Initialize
503 //
504 G4double distance = kInfinity;
505
506 // Find intersections and choose nearest one
507 //
508 G4VTwistSurface *surfaces[6];
509
510 surfaces[0] = fSide0;
511 surfaces[1] = fSide90 ;
512 surfaces[2] = fSide180 ;
513 surfaces[3] = fSide270 ;
514 surfaces[4] = fLowerEndcap;
515 surfaces[5] = fUpperEndcap;
516
517 G4ThreeVector xx;
518 G4ThreeVector bestxx;
519 for (const auto & surface : surfaces)
520 {
521 #ifdef G4TWISTDEBUG
522 G4cout << G4endl << "surface " << &surface - &*surfaces << ": " << G4endl << G4endl ;
523 #endif
524 G4double tmpdistance = surface->DistanceToIn(p, v, xx);
525 #ifdef G4TWISTDEBUG
526 G4cout << "Solid DistanceToIn : distance = " << tmpdistance << G4endl ;
527 G4cout << "intersection point = " << xx << G4endl ;
528 #endif
529 if (tmpdistance < distance)
530 {
531 distance = tmpdistance;
532 bestxx = xx;
533 }
534 }
535
536 #ifdef G4TWISTDEBUG
537 G4cout << "best distance = " << distance << G4endl ;
538 #endif
539
540 // timer.Stop();
541 return distance;
542 }
543
544
545 //=====================================================================
546 //* DistanceToIn (p) --------------------------------------------------
547
548 G4double G4VTwistedFaceted::DistanceToIn (const G4ThreeVector& p) const
549 {
550 // DistanceToIn(p):
551 // Calculate distance to surface of shape from `outside',
552 // allowing for tolerance
553 //
554
555 //
556 // Calculate DistanceToIn(p)
557 //
558
559 EInside currentside = Inside(p);
560
561 switch (currentside)
562 {
563 case (kInside) :
564 {
565 }
566
567 case (kSurface) :
568 {
569 return 0;
570 }
571
572 case (kOutside) :
573 {
574 // Initialize
575 //
576 G4double distance = kInfinity;
577
578 // Find intersections and choose nearest one
579 //
580 G4VTwistSurface* surfaces[6];
581
582 surfaces[0] = fSide0;
583 surfaces[1] = fSide90 ;
584 surfaces[2] = fSide180 ;
585 surfaces[3] = fSide270 ;
586 surfaces[4] = fLowerEndcap;
587 surfaces[5] = fUpperEndcap;
588
589 G4ThreeVector xx;
590 G4ThreeVector bestxx;
591 for (const auto & surface : surfaces)
592 {
593 G4double tmpdistance = surface->DistanceTo(p, xx);
594 if (tmpdistance < distance)
595 {
596 distance = tmpdistance;
597 bestxx = xx;
598 }
599 }
600 return distance;
601 }
602
603 default:
604 {
605 G4Exception("G4VTwistedFaceted::DistanceToIn(p)", "GeomSolids0003",
606 FatalException, "Unknown point location!");
607 }
608 } // switch end
609
610 return 0.;
611 }
612
613
614 //=====================================================================
615 //* DistanceToOut (p, v) ----------------------------------------------
616
617 G4double
618 G4VTwistedFaceted::DistanceToOut( const G4ThreeVector& p,
619 const G4ThreeVector& v,
620 const G4bool calcNorm,
621 G4bool* validNorm,
622 G4ThreeVector* norm ) const
623 {
624 // DistanceToOut (p, v):
625 // Calculate distance to surface of shape from `inside'
626 // along with the v, allowing for tolerance.
627 // The function returns kInfinity if no intersection or
628 // just grazing within tolerance.
629
630 //
631 // Calculate DistanceToOut(p,v)
632 //
633
634 EInside currentside = Inside(p);
635
636 if (currentside == kOutside)
637 {
638 }
639 else if (currentside == kSurface)
640 {
641 // particle is just on a boundary.
642 // if the particle is exiting from the volume, return 0
643 //
644 G4ThreeVector normal = SurfaceNormal(p);
645 if (normal*v > 0)
646 {
647 if (calcNorm)
648 {
649 *norm = normal;
650 *validNorm = true;
651 }
652 // timer.Stop();
653 return 0;
654 }
655 }
656
657 // now, we can take smallest positive distance.
658
659 // Initialize
660 G4double distance = kInfinity;
661
662 // find intersections and choose nearest one.
663 G4VTwistSurface *surfaces[6];
664
665 surfaces[0] = fSide0;
666 surfaces[1] = fSide90 ;
667 surfaces[2] = fSide180 ;
668 surfaces[3] = fSide270 ;
669 surfaces[4] = fLowerEndcap;
670 surfaces[5] = fUpperEndcap;
671
672 G4int besti = -1;
673 G4ThreeVector xx;
674 G4ThreeVector bestxx;
675 for (auto i=0; i<6 ; ++i)
676 {
677 G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
678 if (tmpdistance < distance)
679 {
680 distance = tmpdistance;
681 bestxx = xx;
682 besti = i;
683 }
684 }
685
686 if (calcNorm)
687 {
688 if (besti != -1)
689 {
690 *norm = (surfaces[besti]->GetNormal(p, true));
691 *validNorm = surfaces[besti]->IsValidNorm();
692 }
693 }
694
695 return distance;
696 }
697
698
699 //=====================================================================
700 //* DistanceToOut (p) -------------------------------------------------
701
702 G4double G4VTwistedFaceted::DistanceToOut( const G4ThreeVector& p ) const
703 {
704 // DistanceToOut(p):
705 // Calculate distance to surface of shape from `inside',
706 // allowing for tolerance
707
708 //
709 // Calculate DistanceToOut(p)
710 //
711
712 EInside currentside = Inside(p);
713 G4double retval = kInfinity;
714
715 switch (currentside)
716 {
717 case (kOutside) :
718 {
719 #ifdef G4SPECSDEBUG
720 G4int oldprc = G4cout.precision(16) ;
721 G4cout << G4endl ;
722 DumpInfo();
723 G4cout << "Position:" << G4endl << G4endl ;
724 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
725 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
726 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
727 G4cout.precision(oldprc) ;
728 G4Exception("G4VTwistedFaceted::DistanceToOut(p)", "GeomSolids1002",
729 JustWarning, "Point p is outside !?" );
730 #endif
731 break;
732 }
733 case (kSurface) :
734 {
735 retval = 0;
736 break;
737 }
738
739 case (kInside) :
740 {
741 // Initialize
742 //
743 G4double distance = kInfinity;
744
745 // find intersections and choose nearest one
746 //
747 G4VTwistSurface* surfaces[6];
748
749 surfaces[0] = fSide0;
750 surfaces[1] = fSide90 ;
751 surfaces[2] = fSide180 ;
752 surfaces[3] = fSide270 ;
753 surfaces[4] = fLowerEndcap;
754 surfaces[5] = fUpperEndcap;
755
756 G4ThreeVector xx;
757 G4ThreeVector bestxx;
758 for (const auto & surface : surfaces)
759 {
760 G4double tmpdistance = surface->DistanceTo(p, xx);
761 if (tmpdistance < distance)
762 {
763 distance = tmpdistance;
764 bestxx = xx;
765 }
766 }
767 retval = distance;
768 break;
769 }
770
771 default :
772 {
773 G4Exception("G4VTwistedFaceted::DistanceToOut(p)", "GeomSolids0003",
774 FatalException, "Unknown point location!");
775 break;
776 }
777 } // switch end
778
779 return retval;
780 }
781
782
783 //=====================================================================
784 //* StreamInfo --------------------------------------------------------
785
786 std::ostream& G4VTwistedFaceted::StreamInfo(std::ostream& os) const
787 {
788 //
789 // Stream object contents to an output stream
790 //
791 G4long oldprc = os.precision(16);
792 os << "-----------------------------------------------------------\n"
793 << " *** Dump for solid - " << GetName() << " ***\n"
794 << " ===================================================\n"
795 << " Solid type: G4VTwistedFaceted\n"
796 << " Parameters: \n"
797 << " polar angle theta = " << fTheta/degree << " deg" << G4endl
798 << " azimuthal angle phi = " << fPhi/degree << " deg" << G4endl
799 << " tilt angle alpha = " << fAlph/degree << " deg" << G4endl
800 << " TWIST angle = " << fPhiTwist/degree << " deg" << G4endl
801 << " Half length along y (lower endcap) = " << fDy1/cm << " cm"
802 << G4endl
803 << " Half length along x (lower endcap, bottom) = " << fDx1/cm << " cm"
804 << G4endl
805 << " Half length along x (lower endcap, top) = " << fDx2/cm << " cm"
806 << G4endl
807 << " Half length along y (upper endcap) = " << fDy2/cm << " cm"
808 << G4endl
809 << " Half length along x (upper endcap, bottom) = " << fDx3/cm << " cm"
810 << G4endl
811 << " Half length along x (upper endcap, top) = " << fDx4/cm << " cm"
812 << G4endl
813 << "-----------------------------------------------------------\n";
814 os.precision(oldprc);
815
816 return os;
817 }
818
819
820 //=====================================================================
821 //* DiscribeYourselfTo ------------------------------------------------
822
823 void G4VTwistedFaceted::DescribeYourselfTo (G4VGraphicsScene& scene) const
824 {
825 scene.AddSolid (*this);
826 }
827
828
829 //=====================================================================
830 //* GetExtent ---------------------------------------------------------
831
832 G4VisExtent G4VTwistedFaceted::GetExtent() const
833 {
834 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
835
836 return { -maxRad, maxRad ,
837 -maxRad, maxRad ,
838 -fDz, fDz };
839 }
840
841
842 //=====================================================================
843 //* CreateSurfaces ----------------------------------------------------
844
845 void G4VTwistedFaceted::CreateSurfaces()
846 {
847
848 // create 6 surfaces of TwistedTub.
849
850 if ( fDx1 == fDx2 && fDx3 == fDx4 ) // special case : Box
851 {
852 fSide0 = new G4TwistBoxSide("0deg", fPhiTwist, fDz, fTheta, fPhi,
853 fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 0.*deg);
854 fSide180 = new G4TwistBoxSide("180deg", fPhiTwist, fDz, fTheta, fPhi+pi,
855 fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 180.*deg);
856 }
857 else // default general case
858 {
859 fSide0 = new G4TwistTrapAlphaSide("0deg" ,fPhiTwist, fDz, fTheta,
860 fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
861 fSide180 = new G4TwistTrapAlphaSide("180deg", fPhiTwist, fDz, fTheta,
862 fPhi+pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*deg);
863 }
864
865 // create parallel sides
866 //
867 fSide90 = new G4TwistTrapParallelSide("90deg", fPhiTwist, fDz, fTheta,
868 fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
869 fSide270 = new G4TwistTrapParallelSide("270deg", fPhiTwist, fDz, fTheta,
870 fPhi+pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*deg);
871
872 // create endcaps
873 //
874 fUpperEndcap = new G4TwistTrapFlatSide("UpperCap",fPhiTwist, fDx3, fDx4, fDy2,
875 fDz, fAlph, fPhi, fTheta, 1 );
876 fLowerEndcap = new G4TwistTrapFlatSide("LowerCap",fPhiTwist, fDx1, fDx2, fDy1,
877 fDz, fAlph, fPhi, fTheta, -1 );
878
879 // Set neighbour surfaces
880
881 fSide0->SetNeighbours( fSide270 , fLowerEndcap , fSide90 , fUpperEndcap );
882 fSide90->SetNeighbours( fSide0 , fLowerEndcap , fSide180 , fUpperEndcap );
883 fSide180->SetNeighbours(fSide90 , fLowerEndcap , fSide270 , fUpperEndcap );
884 fSide270->SetNeighbours(fSide180 , fLowerEndcap , fSide0 , fUpperEndcap );
885 fUpperEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
886 fLowerEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
887
888 }
889
890 //=====================================================================
891 //* GetCubicVolume ----------------------------------------------------
892
893 G4double G4VTwistedFaceted::GetCubicVolume()
894 {
895 if(fCubicVolume == 0.)
896 {
897 fCubicVolume = ((fDx1 + fDx2 + fDx3 + fDx4)*(fDy1 + fDy2) +
898 (fDx4 + fDx3 - fDx2 - fDx1)*(fDy2 - fDy1)/3)*fDz;
899 }
900 return fCubicVolume;
901 }
902
903 //=====================================================================
904 //* GetLateralFaceArea ------------------------------------------------
905
906 G4double
907 G4VTwistedFaceted::GetLateralFaceArea(const G4TwoVector& p1,
908 const G4TwoVector& p2,
909 const G4TwoVector& p3,
910 const G4TwoVector& p4) const
911 {
912 constexpr G4int NSTEP = 100;
913 constexpr G4double dt = 1./NSTEP;
914
915 G4double h = 2.*fDz;
916 G4double hh = h*h;
917 G4double hTanTheta = h*std::tan(fTheta);
918 G4double x1 = p1.x();
919 G4double y1 = p1.y();
920 G4double x21 = p2.x() - p1.x();
921 G4double y21 = p2.y() - p1.y();
922 G4double x31 = p3.x() - p1.x();
923 G4double y31 = p3.y() - p1.y();
924 G4double x42 = p4.x() - p2.x();
925 G4double y42 = p4.y() - p2.y();
926 G4double x43 = p4.x() - p3.x();
927 G4double y43 = p4.y() - p3.y();
928
929 // check if face is plane (just in case)
930 G4double lmax = std::max(std::max(std::abs(x21),std::abs(y21)),
931 std::max(std::abs(x43),std::abs(y43)));
932 G4double eps = lmax*kCarTolerance;
933 if (std::abs(fPhiTwist) < kCarTolerance &&
934 std::abs(x21*y43 - y21*x43) < eps)
935 {
936 G4double x0 = hTanTheta*std::cos(fPhi);
937 G4double y0 = hTanTheta*std::sin(fPhi);
938 G4ThreeVector A(p4.x() - p1.x() + x0, p4.y() - p1.y() + y0, h);
939 G4ThreeVector B(p3.x() - p2.x() + x0, p3.y() - p2.y() + y0, h);
940 return (A.cross(B)).mag()*0.5;
941 }
942
943 // twisted face
944 G4double area = 0.;
945 for (G4int i = 0; i < NSTEP; ++i)
946 {
947 G4double t = (i + 0.5)*dt;
948 G4double I = x21 + (x42 - x31)*t;
949 G4double J = y21 + (y42 - y31)*t;
950 G4double II = I*I;
951 G4double JJ = J*J;
952 G4double IIJJ = hh*(I*I + J*J);
953
954 G4double ang = fPhi + fPhiTwist*(0.5 - t);
955 G4double A = fPhiTwist*(II + JJ) + x21*y43 - x43*y21;
956 G4double B = fPhiTwist*(I*(x1 + x31*t) + J*(y1 + y31*t)) +
957 hTanTheta*(I*std::sin(ang) - J*std::cos(ang)) +
958 (I*y31 - J*x31);
959
960 G4double invAA = 1./(A*A);
961 G4double sqrtAA = 2.*std::abs(A);
962 G4double invSqrtAA = 1./sqrtAA;
963
964 G4double aa = A*A;
965 G4double bb = 2.*A*B;
966 G4double cc = IIJJ + B*B;
967
968 G4double R1 = std::sqrt(aa + bb + cc);
969 G4double R0 = std::sqrt(cc);
970 G4double log1 = std::log(std::abs(sqrtAA*R1 + 2.*aa + bb));
971 G4double log0 = std::log(std::abs(sqrtAA*R0 + bb));
972
973 area += 0.5*R1 + 0.25*bb*invAA*(R1 - R0) + IIJJ*invSqrtAA*(log1 - log0);
974 }
975 return area*dt;
976 }
977
978 //=====================================================================
979 //* GetSurfaceArea ----------------------------------------------------
980
981 G4double G4VTwistedFaceted::GetSurfaceArea()
982 {
983 if (fSurfaceArea == 0.)
984 {
985 G4TwoVector vv[8];
986 vv[0] = G4TwoVector(-fDx1 - fDy1*fTAlph,-fDy1);
987 vv[1] = G4TwoVector( fDx1 - fDy1*fTAlph,-fDy1);
988 vv[2] = G4TwoVector(-fDx2 + fDy1*fTAlph, fDy1);
989 vv[3] = G4TwoVector( fDx2 + fDy1*fTAlph, fDy1);
990 vv[4] = G4TwoVector(-fDx3 - fDy2*fTAlph,-fDy2);
991 vv[5] = G4TwoVector( fDx3 - fDy2*fTAlph,-fDy2);
992 vv[6] = G4TwoVector(-fDx4 + fDy2*fTAlph, fDy2);
993 vv[7] = G4TwoVector( fDx4 + fDy2*fTAlph, fDy2);
994 fSurfaceArea = 2.*(fDy1*(fDx1 + fDx2) + fDy2*(fDx3 + fDx4)) +
995 GetLateralFaceArea(vv[0], vv[1], vv[4], vv[5]) +
996 GetLateralFaceArea(vv[1], vv[3], vv[5], vv[7]) +
997 GetLateralFaceArea(vv[3], vv[2], vv[7], vv[6]) +
998 GetLateralFaceArea(vv[2], vv[0], vv[6], vv[4]);
999 }
1000 return fSurfaceArea;
1001 }
1002
1003 //=====================================================================
1004 //* GetEntityType -----------------------------------------------------
1005
1006 G4GeometryType G4VTwistedFaceted::GetEntityType() const
1007 {
1008 return {"G4VTwistedFaceted"};
1009 }
1010
1011
1012 //=====================================================================
1013 //* GetPolyhedron -----------------------------------------------------
1014
1015 G4Polyhedron* G4VTwistedFaceted::GetPolyhedron() const
1016 {
1017 if (fpPolyhedron == nullptr ||
1018 fRebuildPolyhedron ||
1019 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1020 fpPolyhedron->GetNumberOfRotationSteps())
1021 {
1022 G4AutoLock l(&polyhedronMutex);
1023 delete fpPolyhedron;
1024 fpPolyhedron = CreatePolyhedron();
1025 fRebuildPolyhedron = false;
1026 l.unlock();
1027 }
1028
1029 return fpPolyhedron;
1030 }
1031
1032
1033 //=====================================================================
1034 //* GetPointInSolid ---------------------------------------------------
1035
1036 G4ThreeVector G4VTwistedFaceted::GetPointInSolid(G4double z) const
1037 {
1038
1039
1040 // this routine is only used for a test
1041 // can be deleted ...
1042
1043 if ( z == fDz ) z -= 0.1*fDz ;
1044 if ( z == -fDz ) z += 0.1*fDz ;
1045
1046 G4double phi = z/(2*fDz)*fPhiTwist ;
1047
1048 return { fdeltaX * phi/fPhiTwist, fdeltaY * phi/fPhiTwist, z };
1049 }
1050
1051
1052 //=====================================================================
1053 //* GetPointOnSurface -------------------------------------------------
1054
1055 G4ThreeVector G4VTwistedFaceted::GetPointOnSurface() const
1056 {
1057
1058 G4double phi = G4RandFlat::shoot(-fPhiTwist/2.,fPhiTwist/2.);
1059 G4double u , umin, umax ; // variable for twisted surfaces
1060 G4double y ; // variable for flat surface (top and bottom)
1061
1062 // Compute the areas. Attention: Only correct for trapezoids
1063 // where the twisting is done along the z-axis. In the general case
1064 // the computed surface area is more difficult. However this simplification
1065 // does not affect the tracking through the solid.
1066
1067 G4double a1 = fSide0->GetSurfaceArea();
1068 G4double a2 = fSide90->GetSurfaceArea();
1069 G4double a3 = fSide180->GetSurfaceArea() ;
1070 G4double a4 = fSide270->GetSurfaceArea() ;
1071 G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1072 G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1073
1074 #ifdef G4TWISTDEBUG
1075 G4cout << "Surface 0 deg = " << a1 << G4endl ;
1076 G4cout << "Surface 90 deg = " << a2 << G4endl ;
1077 G4cout << "Surface 180 deg = " << a3 << G4endl ;
1078 G4cout << "Surface 270 deg = " << a4 << G4endl ;
1079 G4cout << "Surface Lower = " << a5 << G4endl ;
1080 G4cout << "Surface Upper = " << a6 << G4endl ;
1081 #endif
1082
1083 G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1084
1085 if(chose < a1)
1086 {
1087 umin = fSide0->GetBoundaryMin(phi) ;
1088 umax = fSide0->GetBoundaryMax(phi) ;
1089 u = G4RandFlat::shoot(umin,umax) ;
1090
1091 return fSide0->SurfacePoint(phi, u, true) ; // point on 0deg surface
1092 }
1093
1094 else if( (chose >= a1) && (chose < a1 + a2 ) )
1095 {
1096 umin = fSide90->GetBoundaryMin(phi) ;
1097 umax = fSide90->GetBoundaryMax(phi) ;
1098
1099 u = G4RandFlat::shoot(umin,umax) ;
1100
1101 return fSide90->SurfacePoint(phi, u, true); // point on 90deg surface
1102 }
1103 else if( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1104 {
1105 umin = fSide180->GetBoundaryMin(phi) ;
1106 umax = fSide180->GetBoundaryMax(phi) ;
1107 u = G4RandFlat::shoot(umin,umax) ;
1108
1109 return fSide180->SurfacePoint(phi, u, true); // point on 180 deg surface
1110 }
1111 else if( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1112 {
1113 umin = fSide270->GetBoundaryMin(phi) ;
1114 umax = fSide270->GetBoundaryMax(phi) ;
1115 u = G4RandFlat::shoot(umin,umax) ;
1116 return fSide270->SurfacePoint(phi, u, true); // point on 270 deg surface
1117 }
1118 else if( (chose >= a1 + a2 + a3 + a4 ) && (chose < a1 + a2 + a3 + a4 + a5 ) )
1119 {
1120 y = G4RandFlat::shoot(-fDy1,fDy1) ;
1121 umin = fLowerEndcap->GetBoundaryMin(y) ;
1122 umax = fLowerEndcap->GetBoundaryMax(y) ;
1123 u = G4RandFlat::shoot(umin,umax) ;
1124
1125 return fLowerEndcap->SurfacePoint(u,y,true); // point on lower endcap
1126 }
1127 else
1128 {
1129 y = G4RandFlat::shoot(-fDy2,fDy2) ;
1130 umin = fUpperEndcap->GetBoundaryMin(y) ;
1131 umax = fUpperEndcap->GetBoundaryMax(y) ;
1132 u = G4RandFlat::shoot(umin,umax) ;
1133
1134 return fUpperEndcap->SurfacePoint(u,y,true) ; // point on upper endcap
1135
1136 }
1137 }
1138
1139
1140 //=====================================================================
1141 //* CreatePolyhedron --------------------------------------------------
1142
1143 G4Polyhedron* G4VTwistedFaceted::CreatePolyhedron () const
1144 {
1145 // number of meshes
1146 const G4int k =
1147 G4int(G4Polyhedron::GetNumberOfRotationSteps() *
1148 std::abs(fPhiTwist) / twopi) + 2;
1149 const G4int n = k;
1150
1151 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
1152 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
1153
1154 auto ph = new G4Polyhedron;
1155 typedef G4double G4double3[3];
1156 typedef G4int G4int4[4];
1157 auto xyz = new G4double3[nnodes]; // number of nodes
1158 auto faces = new G4int4[nfaces] ; // number of faces
1159
1160 fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
1161 fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
1162 fSide270->GetFacets(k,n,xyz,faces,2) ;
1163 fSide0->GetFacets(k,n,xyz,faces,3) ;
1164 fSide90->GetFacets(k,n,xyz,faces,4) ;
1165 fSide180->GetFacets(k,n,xyz,faces,5) ;
1166
1167 ph->createPolyhedron(nnodes,nfaces,xyz,faces);
1168
1169 return ph;
1170 }
1171