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