Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 106567 2017-10-13 09:45:14Z 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 absPhiTwist = std::abs(fPhiTwist);
908  G4double dA = std::max(fDPhi,absPhiTwist);
909  const G4int k =
911  const G4int n =
912  G4int(G4Polyhedron::GetNumberOfRotationSteps() * absPhiTwist / twopi) + 2;
913 
914  const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
915  const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
916 
917  G4Polyhedron *ph=new G4Polyhedron;
918  typedef G4double G4double3[3];
919  typedef G4int G4int4[4];
920  G4double3* xyz = new G4double3[nnodes]; // number of nodes
921  G4int4* faces = new G4int4[nfaces] ; // number of faces
922  fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
923  fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
924  fInnerHype->GetFacets(k,n,xyz,faces,2) ;
925  fFormerTwisted->GetFacets(k,n,xyz,faces,3) ;
926  fOuterHype->GetFacets(k,n,xyz,faces,4) ;
927  fLatterTwisted->GetFacets(k,n,xyz,faces,5) ;
928 
929  ph->createPolyhedron(nnodes,nfaces,xyz,faces);
930 
931  delete[] xyz;
932  delete[] faces;
933 
934  return ph;
935 }
936 
937 //=====================================================================
938 //* GetPolyhedron -----------------------------------------------------
939 
941 {
942  if (!fpPolyhedron ||
943  fRebuildPolyhedron ||
945  fpPolyhedron->GetNumberOfRotationSteps())
946  {
947  G4AutoLock l(&polyhedronMutex);
948  delete fpPolyhedron;
949  fpPolyhedron = CreatePolyhedron();
950  fRebuildPolyhedron = false;
951  l.unlock();
952  }
953  return fpPolyhedron;
954 }
955 
956 //=====================================================================
957 //* CreateSurfaces ----------------------------------------------------
958 
959 void G4TwistedTubs::CreateSurfaces()
960 {
961  // create 6 surfaces of TwistedTub
962 
963  G4ThreeVector x0(0, 0, fEndZ[0]);
964  G4ThreeVector n (0, 0, -1);
965 
966  fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap",
967  fEndInnerRadius, fEndOuterRadius,
968  fDPhi, fEndPhi, fEndZ, -1) ;
969 
970  fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap",
971  fEndInnerRadius, fEndOuterRadius,
972  fDPhi, fEndPhi, fEndZ, 1) ;
973 
974  G4RotationMatrix rotHalfDPhi;
975  rotHalfDPhi.rotateZ(0.5*fDPhi);
976 
977  fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
978  fEndInnerRadius, fEndOuterRadius,
979  fDPhi, fEndPhi, fEndZ,
980  fInnerRadius, fOuterRadius, fKappa,
981  1 ) ;
982  fFormerTwisted = new G4TwistTubsSide("FormerTwisted",
983  fEndInnerRadius, fEndOuterRadius,
984  fDPhi, fEndPhi, fEndZ,
985  fInnerRadius, fOuterRadius, fKappa,
986  -1 ) ;
987 
988  fInnerHype = new G4TwistTubsHypeSide("InnerHype",
989  fEndInnerRadius, fEndOuterRadius,
990  fDPhi, fEndPhi, fEndZ,
991  fInnerRadius, fOuterRadius,fKappa,
992  fTanInnerStereo, fTanOuterStereo, -1) ;
993  fOuterHype = new G4TwistTubsHypeSide("OuterHype",
994  fEndInnerRadius, fEndOuterRadius,
995  fDPhi, fEndPhi, fEndZ,
996  fInnerRadius, fOuterRadius,fKappa,
997  fTanInnerStereo, fTanOuterStereo, 1) ;
998 
999 
1000  // set neighbour surfaces
1001  //
1002  fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
1003  fOuterHype, fFormerTwisted);
1004  fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
1005  fOuterHype, fFormerTwisted);
1006  fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
1007  fOuterHype, fUpperEndcap);
1008  fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
1009  fOuterHype, fUpperEndcap);
1010  fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
1011  fFormerTwisted, fUpperEndcap);
1012  fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
1013  fFormerTwisted, fUpperEndcap);
1014 }
1015 
1016 
1017 //=====================================================================
1018 //* GetEntityType -----------------------------------------------------
1019 
1021 {
1022  return G4String("G4TwistedTubs");
1023 }
1024 
1025 //=====================================================================
1026 //* Clone -------------------------------------------------------------
1027 
1029 {
1030  return new G4TwistedTubs(*this);
1031 }
1032 
1033 //=====================================================================
1034 //* GetCubicVolume ----------------------------------------------------
1035 
1037 {
1038  if(fCubicVolume != 0.) {;}
1039  else { fCubicVolume = fDPhi*fZHalfLength*(fOuterRadius*fOuterRadius
1040  -fInnerRadius*fInnerRadius); }
1041  return fCubicVolume;
1042 }
1043 
1044 //=====================================================================
1045 //* GetSurfaceArea ----------------------------------------------------
1046 
1048 {
1049  if(fSurfaceArea != 0.) {;}
1050  else { fSurfaceArea = G4VSolid::GetSurfaceArea(); }
1051  return fSurfaceArea;
1052 }
1053 
1054 //=====================================================================
1055 //* GetPointOnSurface -------------------------------------------------
1056 
1058 {
1059 
1060  G4double z = G4RandFlat::shoot(fEndZ[0],fEndZ[1]);
1061  G4double phi , phimin, phimax ;
1062  G4double x , xmin, xmax ;
1063  G4double r , rmin, rmax ;
1064 
1065  G4double a1 = fOuterHype->GetSurfaceArea() ;
1066  G4double a2 = fInnerHype->GetSurfaceArea() ;
1067  G4double a3 = fLatterTwisted->GetSurfaceArea() ;
1068  G4double a4 = fFormerTwisted->GetSurfaceArea() ;
1069  G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1070  G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1071 
1072  G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1073 
1074  if(chose < a1)
1075  {
1076 
1077  phimin = fOuterHype->GetBoundaryMin(z) ;
1078  phimax = fOuterHype->GetBoundaryMax(z) ;
1079  phi = G4RandFlat::shoot(phimin,phimax) ;
1080 
1081  return fOuterHype->SurfacePoint(phi,z,true) ;
1082 
1083  }
1084  else if ( (chose >= a1) && (chose < a1 + a2 ) )
1085  {
1086 
1087  phimin = fInnerHype->GetBoundaryMin(z) ;
1088  phimax = fInnerHype->GetBoundaryMax(z) ;
1089  phi = G4RandFlat::shoot(phimin,phimax) ;
1090 
1091  return fInnerHype->SurfacePoint(phi,z,true) ;
1092 
1093  }
1094  else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1095  {
1096 
1097  xmin = fLatterTwisted->GetBoundaryMin(z) ;
1098  xmax = fLatterTwisted->GetBoundaryMax(z) ;
1099  x = G4RandFlat::shoot(xmin,xmax) ;
1100 
1101  return fLatterTwisted->SurfacePoint(x,z,true) ;
1102 
1103  }
1104  else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1105  {
1106 
1107  xmin = fFormerTwisted->GetBoundaryMin(z) ;
1108  xmax = fFormerTwisted->GetBoundaryMax(z) ;
1109  x = G4RandFlat::shoot(xmin,xmax) ;
1110 
1111  return fFormerTwisted->SurfacePoint(x,z,true) ;
1112  }
1113  else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
1114  {
1115  rmin = GetEndInnerRadius(0) ;
1116  rmax = GetEndOuterRadius(0) ;
1117  r = std::sqrt(G4RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
1118 
1119  phimin = fLowerEndcap->GetBoundaryMin(r) ;
1120  phimax = fLowerEndcap->GetBoundaryMax(r) ;
1121  phi = G4RandFlat::shoot(phimin,phimax) ;
1122 
1123  return fLowerEndcap->SurfacePoint(phi,r,true) ;
1124  }
1125  else
1126  {
1127  rmin = GetEndInnerRadius(1) ;
1128  rmax = GetEndOuterRadius(1) ;
1129  r = rmin + (rmax-rmin)*std::sqrt(G4RandFlat::shoot());
1130 
1131  phimin = fUpperEndcap->GetBoundaryMin(r) ;
1132  phimax = fUpperEndcap->GetBoundaryMax(r) ;
1133  phi = G4RandFlat::shoot(phimin,phimax) ;
1134 
1135  return fUpperEndcap->SurfacePoint(phi,r,true) ;
1136  }
1137 }
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
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
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
tuple x
Definition: test.py:50
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
const G4int n
G4VSolid * Clone() const
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
virtual G4double GetBoundaryMax(G4double)=0
tuple v
Definition: test.py:18
string pname
Definition: eplot.py:33
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
tuple z
Definition: test.py:28
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)