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