ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TwistedTubs.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4TwistedTubs.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 // G4TwistedTubs implementation
27 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp), created.
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
30 // from original version in Jupiter-2.5.02 application.
31 // --------------------------------------------------------------------
32 
33 #include "G4TwistedTubs.hh"
34 
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4GeometryTolerance.hh"
38 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh"
40 #include "G4BoundingEnvelope.hh"
41 #include "G4ClippablePolygon.hh"
42 #include "G4VPVParameterisation.hh"
43 #include "meshdefs.hh"
44 
45 #include "G4VGraphicsScene.hh"
46 #include "G4Polyhedron.hh"
47 #include "G4VisExtent.hh"
48 
49 #include "Randomize.hh"
50 
51 #include "G4AutoLock.hh"
52 
53 namespace
54 {
55  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
56 }
57 
58 //=====================================================================
59 //* constructors ------------------------------------------------------
60 
61 G4TwistedTubs::G4TwistedTubs(const G4String& pname,
62  G4double twistedangle,
63  G4double endinnerrad,
64  G4double endouterrad,
65  G4double halfzlen,
66  G4double dphi)
67  : G4VSolid(pname), fDPhi(dphi),
68  fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
69  fFormerTwisted(0), fInnerHype(0), fOuterHype(0)
70 {
71  if (endinnerrad < DBL_MIN)
72  {
73  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
74  FatalErrorInArgument, "Invalid end-inner-radius!");
75  }
76 
77  G4double sinhalftwist = std::sin(0.5 * twistedangle);
78 
79  G4double endinnerradX = endinnerrad * sinhalftwist;
80  G4double innerrad = std::sqrt( endinnerrad * endinnerrad
81  - endinnerradX * endinnerradX );
82 
83  G4double endouterradX = endouterrad * sinhalftwist;
84  G4double outerrad = std::sqrt( endouterrad * endouterrad
85  - endouterradX * endouterradX );
86 
87  // temporary treatment!!
88  SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
89  CreateSurfaces();
90 }
91 
92 G4TwistedTubs::G4TwistedTubs(const G4String& pname,
93  G4double twistedangle,
94  G4double endinnerrad,
95  G4double endouterrad,
96  G4double halfzlen,
97  G4int nseg,
98  G4double totphi)
99  : G4VSolid(pname),
100  fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
101  fFormerTwisted(0), fInnerHype(0), fOuterHype(0)
102 {
103 
104  if (!nseg)
105  {
106  std::ostringstream message;
107  message << "Invalid number of segments." << G4endl
108  << " nseg = " << nseg;
109  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
110  FatalErrorInArgument, message);
111  }
112  if (totphi == DBL_MIN || endinnerrad < DBL_MIN)
113  {
114  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
115  FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
116  }
117 
118  G4double sinhalftwist = std::sin(0.5 * twistedangle);
119 
120  G4double endinnerradX = endinnerrad * sinhalftwist;
121  G4double innerrad = std::sqrt( endinnerrad * endinnerrad
122  - endinnerradX * endinnerradX );
123 
124  G4double endouterradX = endouterrad * sinhalftwist;
125  G4double outerrad = std::sqrt( endouterrad * endouterrad
126  - endouterradX * endouterradX );
127 
128  // temporary treatment!!
129  fDPhi = totphi / nseg;
130  SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
131  CreateSurfaces();
132 }
133 
134 G4TwistedTubs::G4TwistedTubs(const G4String& pname,
135  G4double twistedangle,
136  G4double innerrad,
137  G4double outerrad,
138  G4double negativeEndz,
139  G4double positiveEndz,
140  G4double dphi)
141  : G4VSolid(pname), fDPhi(dphi),
142  fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
143  fFormerTwisted(0), fInnerHype(0), fOuterHype(0)
144 {
145  if (innerrad < DBL_MIN)
146  {
147  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
148  FatalErrorInArgument, "Invalid end-inner-radius!");
149  }
150 
151  SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
152  CreateSurfaces();
153 }
154 
155 G4TwistedTubs::G4TwistedTubs(const G4String& pname,
156  G4double twistedangle,
157  G4double innerrad,
158  G4double outerrad,
159  G4double negativeEndz,
160  G4double positiveEndz,
161  G4int nseg,
162  G4double totphi)
163  : G4VSolid(pname),
164  fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
165  fFormerTwisted(0), fInnerHype(0), fOuterHype(0)
166 {
167  if (!nseg)
168  {
169  std::ostringstream message;
170  message << "Invalid number of segments." << G4endl
171  << " nseg = " << nseg;
172  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
173  FatalErrorInArgument, message);
174  }
175  if (totphi == DBL_MIN || innerrad < DBL_MIN)
176  {
177  G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
178  FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
179  }
180 
181  fDPhi = totphi / nseg;
182  SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
183  CreateSurfaces();
184 }
185 
186 //=====================================================================
187 //* Fake default constructor ------------------------------------------
188 
189 G4TwistedTubs::G4TwistedTubs( __void__& a )
190  : G4VSolid(a), fPhiTwist(0.), fInnerRadius(0.), fOuterRadius(0.), fDPhi(0.),
191  fZHalfLength(0.), fInnerStereo(0.), fOuterStereo(0.), fTanInnerStereo(0.),
192  fTanOuterStereo(0.), fKappa(0.), fInnerRadius2(0.), fOuterRadius2(0.),
193  fTanInnerStereo2(0.), fTanOuterStereo2(0.), fLowerEndcap(0), fUpperEndcap(0),
194  fLatterTwisted(0), fFormerTwisted(0), fInnerHype(0), fOuterHype(0)
195 {
196 }
197 
198 //=====================================================================
199 //* destructor --------------------------------------------------------
200 
201 G4TwistedTubs::~G4TwistedTubs()
202 {
203  if (fLowerEndcap) { delete fLowerEndcap; }
204  if (fUpperEndcap) { delete fUpperEndcap; }
205  if (fLatterTwisted) { delete fLatterTwisted; }
206  if (fFormerTwisted) { delete fFormerTwisted; }
207  if (fInnerHype) { delete fInnerHype; }
208  if (fOuterHype) { delete fOuterHype; }
209  if (fpPolyhedron) { delete fpPolyhedron; fpPolyhedron = nullptr; }
210 }
211 
212 //=====================================================================
213 //* Copy constructor --------------------------------------------------
214 
215 G4TwistedTubs::G4TwistedTubs(const G4TwistedTubs& rhs)
216  : G4VSolid(rhs), fPhiTwist(rhs.fPhiTwist),
217  fInnerRadius(rhs.fInnerRadius), fOuterRadius(rhs.fOuterRadius),
218  fDPhi(rhs.fDPhi), fZHalfLength(rhs.fZHalfLength),
219  fInnerStereo(rhs.fInnerStereo), fOuterStereo(rhs.fOuterStereo),
220  fTanInnerStereo(rhs.fTanInnerStereo), fTanOuterStereo(rhs.fTanOuterStereo),
221  fKappa(rhs.fKappa), fInnerRadius2(rhs.fInnerRadius2),
222  fOuterRadius2(rhs.fOuterRadius2), fTanInnerStereo2(rhs.fTanInnerStereo2),
223  fTanOuterStereo2(rhs.fTanOuterStereo2),
224  fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0), fFormerTwisted(0),
225  fInnerHype(0), fOuterHype(0),
226  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
227  fLastInside(rhs.fLastInside), fLastNormal(rhs.fLastNormal),
228  fLastDistanceToIn(rhs.fLastDistanceToIn),
229  fLastDistanceToOut(rhs.fLastDistanceToOut),
230  fLastDistanceToInWithV(rhs.fLastDistanceToInWithV),
231  fLastDistanceToOutWithV(rhs.fLastDistanceToOutWithV)
232 {
233  for (auto i=0; i<2; ++i)
234  {
235  fEndZ[i] = rhs.fEndZ[i];
236  fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
237  fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
238  fEndPhi[i] = rhs.fEndPhi[i];
239  fEndZ2[i] = rhs.fEndZ2[i];
240  }
241  CreateSurfaces();
242 }
243 
244 
245 //=====================================================================
246 //* Assignment operator -----------------------------------------------
247 
248 G4TwistedTubs& G4TwistedTubs::operator = (const G4TwistedTubs& rhs)
249 {
250  // Check assignment to self
251  //
252  if (this == &rhs) { return *this; }
253 
254  // Copy base class data
255  //
256  G4VSolid::operator=(rhs);
257 
258  // Copy data
259  //
260  fPhiTwist= rhs.fPhiTwist;
261  fInnerRadius= rhs.fInnerRadius; fOuterRadius= rhs.fOuterRadius;
262  fDPhi= rhs.fDPhi; fZHalfLength= rhs.fZHalfLength;
263  fInnerStereo= rhs.fInnerStereo; fOuterStereo= rhs.fOuterStereo;
264  fTanInnerStereo= rhs.fTanInnerStereo; fTanOuterStereo= rhs.fTanOuterStereo;
265  fKappa= rhs.fKappa; fInnerRadius2= rhs.fInnerRadius2;
266  fOuterRadius2= rhs.fOuterRadius2; fTanInnerStereo2= rhs.fTanInnerStereo2;
267  fTanOuterStereo2= rhs.fTanOuterStereo2;
268  fLowerEndcap= fUpperEndcap= fLatterTwisted= fFormerTwisted= 0;
269  fInnerHype= fOuterHype= 0;
270  fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
271  fLastInside= rhs.fLastInside; fLastNormal= rhs.fLastNormal;
272  fLastDistanceToIn= rhs.fLastDistanceToIn;
273  fLastDistanceToOut= rhs.fLastDistanceToOut;
274  fLastDistanceToInWithV= rhs.fLastDistanceToInWithV;
275  fLastDistanceToOutWithV= rhs.fLastDistanceToOutWithV;
276 
277  for (auto i=0; i<2; ++i)
278  {
279  fEndZ[i] = rhs.fEndZ[i];
280  fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
281  fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
282  fEndPhi[i] = rhs.fEndPhi[i];
283  fEndZ2[i] = rhs.fEndZ2[i];
284  }
285 
286  CreateSurfaces();
287  fRebuildPolyhedron = false;
288  delete fpPolyhedron; fpPolyhedron = nullptr;
289 
290  return *this;
291 }
292 
293 //=====================================================================
294 //* ComputeDimensions -------------------------------------------------
295 
296 void G4TwistedTubs::ComputeDimensions(G4VPVParameterisation* /* p */ ,
297  const G4int /* n */ ,
298  const G4VPhysicalVolume* /* pRep */ )
299 {
300  G4Exception("G4TwistedTubs::ComputeDimensions()",
301  "GeomSolids0001", FatalException,
302  "G4TwistedTubs does not support Parameterisation.");
303 }
304 
305 //=====================================================================
306 //* BoundingLimits ----------------------------------------------------
307 
308 void G4TwistedTubs::BoundingLimits(G4ThreeVector& pMin,
309  G4ThreeVector& pMax) const
310 {
311  G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
312  fEndOuterRadius[0] : fEndOuterRadius[1]);
313  pMin.set(-maxEndOuterRad,-maxEndOuterRad,-fZHalfLength);
314  pMax.set( maxEndOuterRad, maxEndOuterRad, fZHalfLength);
315 
316  // Check correctness of the bounding box
317  //
318  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
319  {
320  std::ostringstream message;
321  message << "Bad bounding box (min >= max) for solid: "
322  << GetName() << " !"
323  << "\npMin = " << pMin
324  << "\npMax = " << pMax;
325  G4Exception("G4TwistedTubs::BoundingLimits()", "GeomMgt0001",
326  JustWarning, message);
327  DumpInfo();
328  }
329 }
330 
331 //=====================================================================
332 //* CalculateExtent ---------------------------------------------------
333 
334 G4bool
335 G4TwistedTubs::CalculateExtent(const EAxis pAxis,
336  const G4VoxelLimits& pVoxelLimit,
337  const G4AffineTransform& pTransform,
338  G4double& pMin, G4double& pMax) const
339 {
340  G4ThreeVector bmin, bmax;
341 
342  // Get bounding box
343  BoundingLimits(bmin,bmax);
344 
345  // Find extent
346  G4BoundingEnvelope bbox(bmin,bmax);
347  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
348 }
349 
350 
351 //=====================================================================
352 //* Inside ------------------------------------------------------------
353 
354 EInside G4TwistedTubs::Inside(const G4ThreeVector& p) const
355 {
356 
357  const G4double halftol
358  = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
359  // static G4int timerid = -1;
360  // G4Timer timer(timerid, "G4TwistedTubs", "Inside");
361  // timer.Start();
362 
363  G4ThreeVector *tmpp;
364  EInside *tmpinside;
365  if (fLastInside.p == p)
366  {
367  return fLastInside.inside;
368  }
369  else
370  {
371  tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
372  tmpinside = const_cast<EInside*>(&(fLastInside.inside));
373  tmpp->set(p.x(), p.y(), p.z());
374  }
375 
376  EInside outerhypearea = ((G4TwistTubsHypeSide *)fOuterHype)->Inside(p);
377  G4double innerhyperho = ((G4TwistTubsHypeSide *)fInnerHype)->GetRhoAtPZ(p);
378  G4double distanceToOut = p.getRho() - innerhyperho; // +ve: inside
379 
380  if ((outerhypearea == kOutside) || (distanceToOut < -halftol))
381  {
382  *tmpinside = kOutside;
383  }
384  else if (outerhypearea == kSurface)
385  {
386  *tmpinside = kSurface;
387  }
388  else
389  {
390  if (distanceToOut <= halftol)
391  {
392  *tmpinside = kSurface;
393  }
394  else
395  {
396  *tmpinside = kInside;
397  }
398  }
399 
400  return fLastInside.inside;
401 }
402 
403 //=====================================================================
404 //* SurfaceNormal -----------------------------------------------------
405 
406 G4ThreeVector G4TwistedTubs::SurfaceNormal(const G4ThreeVector& p) const
407 {
408  //
409  // return the normal unit vector to the Hyperbolical Surface at a point
410  // p on (or nearly on) the surface
411  //
412  // Which of the three or four surfaces are we closest to?
413  //
414 
415  if (fLastNormal.p == p)
416  {
417  return fLastNormal.vec;
418  }
419  G4ThreeVector *tmpp =
420  const_cast<G4ThreeVector*>(&(fLastNormal.p));
421  G4ThreeVector *tmpnormal =
422  const_cast<G4ThreeVector*>(&(fLastNormal.vec));
423  G4VTwistSurface **tmpsurface =
424  const_cast<G4VTwistSurface**>(fLastNormal.surface);
425  tmpp->set(p.x(), p.y(), p.z());
426 
427  G4double distance = kInfinity;
428 
429  G4VTwistSurface *surfaces[6];
430  surfaces[0] = fLatterTwisted;
431  surfaces[1] = fFormerTwisted;
432  surfaces[2] = fInnerHype;
433  surfaces[3] = fOuterHype;
434  surfaces[4] = fLowerEndcap;
435  surfaces[5] = fUpperEndcap;
436 
437  G4ThreeVector xx;
438  G4ThreeVector bestxx;
439  G4int besti = -1;
440  for (auto i=0; i<6; ++i)
441  {
442  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
443  if (tmpdistance < distance)
444  {
445  distance = tmpdistance;
446  bestxx = xx;
447  besti = i;
448  }
449  }
450 
451  tmpsurface[0] = surfaces[besti];
452  *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
453 
454  return fLastNormal.vec;
455 }
456 
457 //=====================================================================
458 //* DistanceToIn (p, v) -----------------------------------------------
459 
460 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p,
461  const G4ThreeVector& v ) const
462 {
463 
464  // DistanceToIn (p, v):
465  // Calculate distance to surface of shape from `outside'
466  // along with the v, allowing for tolerance.
467  // The function returns kInfinity if no intersection or
468  // just grazing within tolerance.
469 
470  //
471  // checking last value
472  //
473 
474  G4ThreeVector* tmpp;
475  G4ThreeVector* tmpv;
476  G4double* tmpdist;
477  if ((fLastDistanceToInWithV.p == p) && (fLastDistanceToInWithV.vec == v))
478  {
479  return fLastDistanceToIn.value;
480  }
481  else
482  {
483  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
484  tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
485  tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
486  tmpp->set(p.x(), p.y(), p.z());
487  tmpv->set(v.x(), v.y(), v.z());
488  }
489 
490  //
491  // Calculate DistanceToIn(p,v)
492  //
493 
494  EInside currentside = Inside(p);
495 
496  if (currentside == kInside)
497  {
498  }
499  else
500  {
501  if (currentside == kSurface)
502  {
503  // particle is just on a boundary.
504  // If the particle is entering to the volume, return 0.
505  //
506  G4ThreeVector normal = SurfaceNormal(p);
507  if (normal*v < 0)
508  {
509  *tmpdist = 0.;
510  return fLastDistanceToInWithV.value;
511  }
512  }
513  }
514 
515  // now, we can take smallest positive distance.
516 
517  // Initialize
518  //
519  G4double distance = kInfinity;
520 
521  // find intersections and choose nearest one.
522  //
523  G4VTwistSurface* surfaces[6];
524  surfaces[0] = fLowerEndcap;
525  surfaces[1] = fUpperEndcap;
526  surfaces[2] = fLatterTwisted;
527  surfaces[3] = fFormerTwisted;
528  surfaces[4] = fInnerHype;
529  surfaces[5] = fOuterHype;
530 
531  G4ThreeVector xx;
532  G4ThreeVector bestxx;
533  for (auto i=0; i<6; ++i)
534  {
535  G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
536  if (tmpdistance < distance)
537  {
538  distance = tmpdistance;
539  bestxx = xx;
540  }
541  }
542  *tmpdist = distance;
543 
544  return fLastDistanceToInWithV.value;
545 }
546 
547 //=====================================================================
548 //* DistanceToIn (p) --------------------------------------------------
549 
550 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p) const
551 {
552  // DistanceToIn(p):
553  // Calculate distance to surface of shape from `outside',
554  // allowing for tolerance
555 
556  //
557  // checking last value
558  //
559 
560  G4ThreeVector* tmpp;
561  G4double* tmpdist;
562  if (fLastDistanceToIn.p == p)
563  {
564  return fLastDistanceToIn.value;
565  }
566  else
567  {
568  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
569  tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
570  tmpp->set(p.x(), p.y(), p.z());
571  }
572 
573  //
574  // Calculate DistanceToIn(p)
575  //
576 
577  EInside currentside = Inside(p);
578 
579  switch (currentside)
580  {
581  case (kInside) :
582  {}
583  case (kSurface) :
584  {
585  *tmpdist = 0.;
586  return fLastDistanceToIn.value;
587  }
588  case (kOutside) :
589  {
590  // Initialize
591  G4double distance = kInfinity;
592 
593  // find intersections and choose nearest one.
594  G4VTwistSurface *surfaces[6];
595  surfaces[0] = fLowerEndcap;
596  surfaces[1] = fUpperEndcap;
597  surfaces[2] = fLatterTwisted;
598  surfaces[3] = fFormerTwisted;
599  surfaces[4] = fInnerHype;
600  surfaces[5] = fOuterHype;
601 
602  G4ThreeVector xx;
603  G4ThreeVector bestxx;
604  for (auto i=0; i<6; ++i)
605  {
606  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
607  if (tmpdistance < distance)
608  {
609  distance = tmpdistance;
610  bestxx = xx;
611  }
612  }
613  *tmpdist = distance;
614  return fLastDistanceToIn.value;
615  }
616  default :
617  {
618  G4Exception("G4TwistedTubs::DistanceToIn(p)", "GeomSolids0003",
619  FatalException, "Unknown point location!");
620  }
621  } // switch end
622 
623  return kInfinity;
624 }
625 
626 //=====================================================================
627 //* DistanceToOut (p, v) ----------------------------------------------
628 
629 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p,
630  const G4ThreeVector& v,
631  const G4bool calcNorm,
632  G4bool *validNorm,
633  G4ThreeVector *norm ) const
634 {
635  // DistanceToOut (p, v):
636  // Calculate distance to surface of shape from `inside'
637  // along with the v, allowing for tolerance.
638  // The function returns kInfinity if no intersection or
639  // just grazing within tolerance.
640 
641  //
642  // checking last value
643  //
644 
645  G4ThreeVector* tmpp;
646  G4ThreeVector* tmpv;
647  G4double* tmpdist;
648  if ((fLastDistanceToOutWithV.p == p) && (fLastDistanceToOutWithV.vec == v) )
649  {
650  return fLastDistanceToOutWithV.value;
651  }
652  else
653  {
654  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
655  tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
656  tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
657  tmpp->set(p.x(), p.y(), p.z());
658  tmpv->set(v.x(), v.y(), v.z());
659  }
660 
661  //
662  // Calculate DistanceToOut(p,v)
663  //
664 
665  EInside currentside = Inside(p);
666 
667  if (currentside == kOutside)
668  {
669  }
670  else
671  {
672  if (currentside == kSurface)
673  {
674  // particle is just on a boundary.
675  // If the particle is exiting from the volume, return 0.
676  //
677  G4ThreeVector normal = SurfaceNormal(p);
678  G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
679  if (normal*v > 0)
680  {
681  if (calcNorm)
682  {
683  *norm = (blockedsurface->GetNormal(p, true));
684  *validNorm = blockedsurface->IsValidNorm();
685  }
686  *tmpdist = 0.;
687  return fLastDistanceToOutWithV.value;
688  }
689  }
690  }
691 
692  // now, we can take smallest positive distance.
693 
694  // Initialize
695  //
696  G4double distance = kInfinity;
697 
698  // find intersections and choose nearest one.
699  //
700  G4VTwistSurface* surfaces[6];
701  surfaces[0] = fLatterTwisted;
702  surfaces[1] = fFormerTwisted;
703  surfaces[2] = fInnerHype;
704  surfaces[3] = fOuterHype;
705  surfaces[4] = fLowerEndcap;
706  surfaces[5] = fUpperEndcap;
707 
708  G4int besti = -1;
709  G4ThreeVector xx;
710  G4ThreeVector bestxx;
711  for (auto i=0; i<6; ++i)
712  {
713  G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
714  if (tmpdistance < distance)
715  {
716  distance = tmpdistance;
717  bestxx = xx;
718  besti = i;
719  }
720  }
721 
722  if (calcNorm)
723  {
724  if (besti != -1)
725  {
726  *norm = (surfaces[besti]->GetNormal(p, true));
727  *validNorm = surfaces[besti]->IsValidNorm();
728  }
729  }
730 
731  *tmpdist = distance;
732 
733  return fLastDistanceToOutWithV.value;
734 }
735 
736 
737 //=====================================================================
738 //* DistanceToOut (p) ----------------------------------------------
739 
740 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p ) const
741 {
742  // DistanceToOut(p):
743  // Calculate distance to surface of shape from `inside',
744  // allowing for tolerance
745 
746  //
747  // checking last value
748  //
749 
750  G4ThreeVector* tmpp;
751  G4double* tmpdist;
752  if (fLastDistanceToOut.p == p)
753  {
754  return fLastDistanceToOut.value;
755  }
756  else
757  {
758  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
759  tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
760  tmpp->set(p.x(), p.y(), p.z());
761  }
762 
763  //
764  // Calculate DistanceToOut(p)
765  //
766 
767  EInside currentside = Inside(p);
768 
769  switch (currentside)
770  {
771  case (kOutside) :
772  {
773  }
774  case (kSurface) :
775  {
776  *tmpdist = 0.;
777  return fLastDistanceToOut.value;
778  }
779  case (kInside) :
780  {
781  // Initialize
782  G4double distance = kInfinity;
783 
784  // find intersections and choose nearest one.
785  G4VTwistSurface* surfaces[6];
786  surfaces[0] = fLatterTwisted;
787  surfaces[1] = fFormerTwisted;
788  surfaces[2] = fInnerHype;
789  surfaces[3] = fOuterHype;
790  surfaces[4] = fLowerEndcap;
791  surfaces[5] = fUpperEndcap;
792 
793  G4ThreeVector xx;
794  G4ThreeVector bestxx;
795  for (auto i=0; i<6; ++i)
796  {
797  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
798  if (tmpdistance < distance)
799  {
800  distance = tmpdistance;
801  bestxx = xx;
802  }
803  }
804  *tmpdist = distance;
805 
806  return fLastDistanceToOut.value;
807  }
808  default :
809  {
810  G4Exception("G4TwistedTubs::DistanceToOut(p)", "GeomSolids0003",
811  FatalException, "Unknown point location!");
812  }
813  } // switch end
814 
815  return 0.;
816 }
817 
818 //=====================================================================
819 //* StreamInfo --------------------------------------------------------
820 
821 std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
822 {
823  //
824  // Stream object contents to an output stream
825  //
826  G4int oldprc = os.precision(16);
827  os << "-----------------------------------------------------------\n"
828  << " *** Dump for solid - " << GetName() << " ***\n"
829  << " ===================================================\n"
830  << " Solid type: G4TwistedTubs\n"
831  << " Parameters: \n"
832  << " -ve end Z : " << fEndZ[0]/mm << " mm \n"
833  << " +ve end Z : " << fEndZ[1]/mm << " mm \n"
834  << " inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
835  << " inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
836  << " outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
837  << " outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
838  << " inner radius (z=0) : " << fInnerRadius/mm << " mm \n"
839  << " outer radius (z=0) : " << fOuterRadius/mm << " mm \n"
840  << " twisted angle : " << fPhiTwist/degree << " degrees \n"
841  << " inner stereo angle : " << fInnerStereo/degree << " degrees \n"
842  << " outer stereo angle : " << fOuterStereo/degree << " degrees \n"
843  << " phi-width of a piece : " << fDPhi/degree << " degrees \n"
844  << "-----------------------------------------------------------\n";
845  os.precision(oldprc);
846 
847  return os;
848 }
849 
850 
851 //=====================================================================
852 //* DiscribeYourselfTo ------------------------------------------------
853 
854 void G4TwistedTubs::DescribeYourselfTo (G4VGraphicsScene& scene) const
855 {
856  scene.AddSolid (*this);
857 }
858 
859 //=====================================================================
860 //* GetExtent ---------------------------------------------------------
861 
862 G4VisExtent G4TwistedTubs::GetExtent() const
863 {
864  // Define the sides of the box into which the G4Tubs instance would fit.
865 
866  G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
867  fEndOuterRadius[0] : fEndOuterRadius[1]);
868  return G4VisExtent( -maxEndOuterRad, maxEndOuterRad,
869  -maxEndOuterRad, maxEndOuterRad,
870  -fZHalfLength, fZHalfLength );
871 }
872 
873 //=====================================================================
874 //* CreatePolyhedron --------------------------------------------------
875 
876 G4Polyhedron* G4TwistedTubs::CreatePolyhedron () const
877 {
878  // number of meshes
879  //
880  G4double absPhiTwist = std::abs(fPhiTwist);
881  G4double dA = std::max(fDPhi,absPhiTwist);
882  const G4int k =
883  G4int(G4Polyhedron::GetNumberOfRotationSteps() * dA / twopi) + 2;
884  const G4int n =
885  G4int(G4Polyhedron::GetNumberOfRotationSteps() * absPhiTwist / twopi) + 2;
886 
887  const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
888  const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
889 
890  G4Polyhedron* ph = new G4Polyhedron;
891  typedef G4double G4double3[3];
892  typedef G4int G4int4[4];
893  G4double3* xyz = new G4double3[nnodes]; // number of nodes
894  G4int4* faces = new G4int4[nfaces] ; // number of faces
895  fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
896  fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
897  fInnerHype->GetFacets(k,n,xyz,faces,2) ;
898  fFormerTwisted->GetFacets(k,n,xyz,faces,3) ;
899  fOuterHype->GetFacets(k,n,xyz,faces,4) ;
900  fLatterTwisted->GetFacets(k,n,xyz,faces,5) ;
901 
902  ph->createPolyhedron(nnodes,nfaces,xyz,faces);
903 
904  delete[] xyz;
905  delete[] faces;
906 
907  return ph;
908 }
909 
910 //=====================================================================
911 //* GetPolyhedron -----------------------------------------------------
912 
913 G4Polyhedron* G4TwistedTubs::GetPolyhedron () const
914 {
915  if (fpPolyhedron == nullptr ||
916  fRebuildPolyhedron ||
917  fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
918  fpPolyhedron->GetNumberOfRotationSteps())
919  {
920  G4AutoLock l(&polyhedronMutex);
921  delete fpPolyhedron;
922  fpPolyhedron = CreatePolyhedron();
923  fRebuildPolyhedron = false;
924  l.unlock();
925  }
926  return fpPolyhedron;
927 }
928 
929 //=====================================================================
930 //* CreateSurfaces ----------------------------------------------------
931 
932 void G4TwistedTubs::CreateSurfaces()
933 {
934  // create 6 surfaces of TwistedTub
935 
936  G4ThreeVector x0(0, 0, fEndZ[0]);
937  G4ThreeVector n (0, 0, -1);
938 
939  fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap",
940  fEndInnerRadius, fEndOuterRadius,
941  fDPhi, fEndPhi, fEndZ, -1) ;
942 
943  fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap",
944  fEndInnerRadius, fEndOuterRadius,
945  fDPhi, fEndPhi, fEndZ, 1) ;
946 
947  G4RotationMatrix rotHalfDPhi;
948  rotHalfDPhi.rotateZ(0.5*fDPhi);
949 
950  fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
951  fEndInnerRadius, fEndOuterRadius,
952  fDPhi, fEndPhi, fEndZ,
953  fInnerRadius, fOuterRadius, fKappa,
954  1 ) ;
955  fFormerTwisted = new G4TwistTubsSide("FormerTwisted",
956  fEndInnerRadius, fEndOuterRadius,
957  fDPhi, fEndPhi, fEndZ,
958  fInnerRadius, fOuterRadius, fKappa,
959  -1 ) ;
960 
961  fInnerHype = new G4TwistTubsHypeSide("InnerHype",
962  fEndInnerRadius, fEndOuterRadius,
963  fDPhi, fEndPhi, fEndZ,
964  fInnerRadius, fOuterRadius,fKappa,
965  fTanInnerStereo, fTanOuterStereo, -1) ;
966  fOuterHype = new G4TwistTubsHypeSide("OuterHype",
967  fEndInnerRadius, fEndOuterRadius,
968  fDPhi, fEndPhi, fEndZ,
969  fInnerRadius, fOuterRadius,fKappa,
970  fTanInnerStereo, fTanOuterStereo, 1) ;
971 
972 
973  // set neighbour surfaces
974  //
975  fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
976  fOuterHype, fFormerTwisted);
977  fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
978  fOuterHype, fFormerTwisted);
979  fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
980  fOuterHype, fUpperEndcap);
981  fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
982  fOuterHype, fUpperEndcap);
983  fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
984  fFormerTwisted, fUpperEndcap);
985  fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
986  fFormerTwisted, fUpperEndcap);
987 }
988 
989 
990 //=====================================================================
991 //* GetEntityType -----------------------------------------------------
992 
993 G4GeometryType G4TwistedTubs::GetEntityType() const
994 {
995  return G4String("G4TwistedTubs");
996 }
997 
998 //=====================================================================
999 //* Clone -------------------------------------------------------------
1000 
1001 G4VSolid* G4TwistedTubs::Clone() const
1002 {
1003  return new G4TwistedTubs(*this);
1004 }
1005 
1006 //=====================================================================
1007 //* GetCubicVolume ----------------------------------------------------
1008 
1009 G4double G4TwistedTubs::GetCubicVolume()
1010 {
1011  if(fCubicVolume != 0.) {;}
1012  else { fCubicVolume = fDPhi*fZHalfLength*(fOuterRadius*fOuterRadius
1013  -fInnerRadius*fInnerRadius); }
1014  return fCubicVolume;
1015 }
1016 
1017 //=====================================================================
1018 //* GetSurfaceArea ----------------------------------------------------
1019 
1020 G4double G4TwistedTubs::GetSurfaceArea()
1021 {
1022  if(fSurfaceArea != 0.) {;}
1023  else { fSurfaceArea = G4VSolid::GetSurfaceArea(); }
1024  return fSurfaceArea;
1025 }
1026 
1027 //=====================================================================
1028 //* GetPointOnSurface -------------------------------------------------
1029 
1030 G4ThreeVector G4TwistedTubs::GetPointOnSurface() const
1031 {
1032 
1033  G4double z = G4RandFlat::shoot(fEndZ[0],fEndZ[1]);
1034  G4double phi , phimin, phimax ;
1035  G4double x , xmin, xmax ;
1036  G4double r , rmin, rmax ;
1037 
1038  G4double a1 = fOuterHype->GetSurfaceArea() ;
1039  G4double a2 = fInnerHype->GetSurfaceArea() ;
1040  G4double a3 = fLatterTwisted->GetSurfaceArea() ;
1041  G4double a4 = fFormerTwisted->GetSurfaceArea() ;
1042  G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1043  G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1044 
1045  G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1046 
1047  if(chose < a1)
1048  {
1049 
1050  phimin = fOuterHype->GetBoundaryMin(z) ;
1051  phimax = fOuterHype->GetBoundaryMax(z) ;
1052  phi = G4RandFlat::shoot(phimin,phimax) ;
1053 
1054  return fOuterHype->SurfacePoint(phi,z,true) ;
1055 
1056  }
1057  else if ( (chose >= a1) && (chose < a1 + a2 ) )
1058  {
1059 
1060  phimin = fInnerHype->GetBoundaryMin(z) ;
1061  phimax = fInnerHype->GetBoundaryMax(z) ;
1062  phi = G4RandFlat::shoot(phimin,phimax) ;
1063 
1064  return fInnerHype->SurfacePoint(phi,z,true) ;
1065 
1066  }
1067  else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1068  {
1069 
1070  xmin = fLatterTwisted->GetBoundaryMin(z) ;
1071  xmax = fLatterTwisted->GetBoundaryMax(z) ;
1072  x = G4RandFlat::shoot(xmin,xmax) ;
1073 
1074  return fLatterTwisted->SurfacePoint(x,z,true) ;
1075 
1076  }
1077  else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1078  {
1079 
1080  xmin = fFormerTwisted->GetBoundaryMin(z) ;
1081  xmax = fFormerTwisted->GetBoundaryMax(z) ;
1082  x = G4RandFlat::shoot(xmin,xmax) ;
1083 
1084  return fFormerTwisted->SurfacePoint(x,z,true) ;
1085  }
1086  else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
1087  {
1088  rmin = GetEndInnerRadius(0) ;
1089  rmax = GetEndOuterRadius(0) ;
1090  r = std::sqrt(G4RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
1091 
1092  phimin = fLowerEndcap->GetBoundaryMin(r) ;
1093  phimax = fLowerEndcap->GetBoundaryMax(r) ;
1094  phi = G4RandFlat::shoot(phimin,phimax) ;
1095 
1096  return fLowerEndcap->SurfacePoint(phi,r,true) ;
1097  }
1098  else
1099  {
1100  rmin = GetEndInnerRadius(1) ;
1101  rmax = GetEndOuterRadius(1) ;
1102  r = rmin + (rmax-rmin)*std::sqrt(G4RandFlat::shoot());
1103 
1104  phimin = fUpperEndcap->GetBoundaryMin(r) ;
1105  phimax = fUpperEndcap->GetBoundaryMax(r) ;
1106  phi = G4RandFlat::shoot(phimin,phimax) ;
1107 
1108  return fUpperEndcap->SurfacePoint(phi,r,true) ;
1109  }
1110 }