Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 //
27 // $Id: G4TwistedTubs.cc 105776 2017-08-17 08:09:09Z 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;
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  G4ThreeVector& pMax) const
333 {
334  G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
336  pMin.set(-maxEndOuterRad,-maxEndOuterRad,-fZHalfLength);
337  pMax.set( maxEndOuterRad, maxEndOuterRad, fZHalfLength);
338 
339  // Check correctness of the bounding box
340  //
341  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
342  {
343  std::ostringstream message;
344  message << "Bad bounding box (min >= max) for solid: "
345  << GetName() << " !"
346  << "\npMin = " << pMin
347  << "\npMax = " << pMax;
348  G4Exception("G4TwistedTubs::BoundingLimits()", "GeomMgt0001",
349  JustWarning, message);
350  DumpInfo();
351  }
352 }
353 
355 //
356 // Calculate extent under transform and specified limit
357 
358 G4bool
360  const G4VoxelLimits& pVoxelLimit,
361  const G4AffineTransform& pTransform,
362  G4double& pMin, G4double& pMax) const
363 {
364  G4ThreeVector bmin, bmax;
365 
366  // Get bounding box
367  BoundingLimits(bmin,bmax);
368 
369  // Find extent
370  G4BoundingEnvelope bbox(bmin,bmax);
371  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
372 }
373 
374 
375 //=====================================================================
376 //* Inside ------------------------------------------------------------
377 
379 {
380 
381  const G4double halftol
383  // static G4int timerid = -1;
384  // G4Timer timer(timerid, "G4TwistedTubs", "Inside");
385  // timer.Start();
386 
387  G4ThreeVector *tmpp;
388  EInside *tmpinside;
389  if (fLastInside.p == p)
390  {
391  return fLastInside.inside;
392  }
393  else
394  {
395  tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
396  tmpinside = const_cast<EInside*>(&(fLastInside.inside));
397  tmpp->set(p.x(), p.y(), p.z());
398  }
399 
400  EInside outerhypearea = ((G4TwistTubsHypeSide *)fOuterHype)->Inside(p);
401  G4double innerhyperho = ((G4TwistTubsHypeSide *)fInnerHype)->GetRhoAtPZ(p);
402  G4double distanceToOut = p.getRho() - innerhyperho; // +ve: inside
403 
404  if ((outerhypearea == kOutside) || (distanceToOut < -halftol))
405  {
406  *tmpinside = kOutside;
407  }
408  else if (outerhypearea == kSurface)
409  {
410  *tmpinside = kSurface;
411  }
412  else
413  {
414  if (distanceToOut <= halftol)
415  {
416  *tmpinside = kSurface;
417  }
418  else
419  {
420  *tmpinside = kInside;
421  }
422  }
423 
424  return fLastInside.inside;
425 }
426 
427 //=====================================================================
428 //* SurfaceNormal -----------------------------------------------------
429 
431 {
432  //
433  // return the normal unit vector to the Hyperbolical Surface at a point
434  // p on (or nearly on) the surface
435  //
436  // Which of the three or four surfaces are we closest to?
437  //
438 
439  if (fLastNormal.p == p)
440  {
441  return fLastNormal.vec;
442  }
443  G4ThreeVector *tmpp =
444  const_cast<G4ThreeVector*>(&(fLastNormal.p));
445  G4ThreeVector *tmpnormal =
446  const_cast<G4ThreeVector*>(&(fLastNormal.vec));
447  G4VTwistSurface **tmpsurface =
448  const_cast<G4VTwistSurface**>(fLastNormal.surface);
449  tmpp->set(p.x(), p.y(), p.z());
450 
451  G4double distance = kInfinity;
452 
453  G4VTwistSurface *surfaces[6];
454  surfaces[0] = fLatterTwisted;
455  surfaces[1] = fFormerTwisted;
456  surfaces[2] = fInnerHype;
457  surfaces[3] = fOuterHype;
458  surfaces[4] = fLowerEndcap;
459  surfaces[5] = fUpperEndcap;
460 
462  G4ThreeVector bestxx;
463  G4int i;
464  G4int besti = -1;
465  for (i=0; i< 6; i++)
466  {
467  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
468  if (tmpdistance < distance)
469  {
470  distance = tmpdistance;
471  bestxx = xx;
472  besti = i;
473  }
474  }
475 
476  tmpsurface[0] = surfaces[besti];
477  *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
478 
479  return fLastNormal.vec;
480 }
481 
482 //=====================================================================
483 //* DistanceToIn (p, v) -----------------------------------------------
484 
486  const G4ThreeVector& v ) const
487 {
488 
489  // DistanceToIn (p, v):
490  // Calculate distance to surface of shape from `outside'
491  // along with the v, allowing for tolerance.
492  // The function returns kInfinity if no intersection or
493  // just grazing within tolerance.
494 
495  //
496  // checking last value
497  //
498 
499  G4ThreeVector *tmpp;
500  G4ThreeVector *tmpv;
501  G4double *tmpdist;
502  if ((fLastDistanceToInWithV.p == p) && (fLastDistanceToInWithV.vec == v))
503  {
504  return fLastDistanceToIn.value;
505  }
506  else
507  {
508  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
509  tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
510  tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
511  tmpp->set(p.x(), p.y(), p.z());
512  tmpv->set(v.x(), v.y(), v.z());
513  }
514 
515  //
516  // Calculate DistanceToIn(p,v)
517  //
518 
519  EInside currentside = Inside(p);
520 
521  if (currentside == kInside)
522  {
523  }
524  else
525  {
526  if (currentside == kSurface)
527  {
528  // particle is just on a boundary.
529  // If the particle is entering to the volume, return 0.
530  //
532  if (normal*v < 0)
533  {
534  *tmpdist = 0;
536  }
537  }
538  }
539 
540  // now, we can take smallest positive distance.
541 
542  // Initialize
543  //
544  G4double distance = kInfinity;
545 
546  // find intersections and choose nearest one.
547  //
548  G4VTwistSurface *surfaces[6];
549  surfaces[0] = fLowerEndcap;
550  surfaces[1] = fUpperEndcap;
551  surfaces[2] = fLatterTwisted;
552  surfaces[3] = fFormerTwisted;
553  surfaces[4] = fInnerHype;
554  surfaces[5] = fOuterHype;
555 
557  G4ThreeVector bestxx;
558  G4int i;
559  for (i=0; i< 6; i++)
560  {
561  G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
562  if (tmpdistance < distance)
563  {
564  distance = tmpdistance;
565  bestxx = xx;
566  }
567  }
568  *tmpdist = distance;
569 
571 }
572 
573 //=====================================================================
574 //* DistanceToIn (p) --------------------------------------------------
575 
577 {
578  // DistanceToIn(p):
579  // Calculate distance to surface of shape from `outside',
580  // allowing for tolerance
581 
582  //
583  // checking last value
584  //
585 
586  G4ThreeVector *tmpp;
587  G4double *tmpdist;
588  if (fLastDistanceToIn.p == p)
589  {
590  return fLastDistanceToIn.value;
591  }
592  else
593  {
594  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
595  tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
596  tmpp->set(p.x(), p.y(), p.z());
597  }
598 
599  //
600  // Calculate DistanceToIn(p)
601  //
602 
603  EInside currentside = Inside(p);
604 
605  switch (currentside)
606  {
607  case (kInside) :
608  {}
609  case (kSurface) :
610  {
611  *tmpdist = 0.;
612  return fLastDistanceToIn.value;
613  }
614  case (kOutside) :
615  {
616  // Initialize
617  G4double distance = kInfinity;
618 
619  // find intersections and choose nearest one.
620  G4VTwistSurface *surfaces[6];
621  surfaces[0] = fLowerEndcap;
622  surfaces[1] = fUpperEndcap;
623  surfaces[2] = fLatterTwisted;
624  surfaces[3] = fFormerTwisted;
625  surfaces[4] = fInnerHype;
626  surfaces[5] = fOuterHype;
627 
628  G4int i;
630  G4ThreeVector bestxx;
631  for (i=0; i< 6; i++)
632  {
633  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
634  if (tmpdistance < distance)
635  {
636  distance = tmpdistance;
637  bestxx = xx;
638  }
639  }
640  *tmpdist = distance;
641  return fLastDistanceToIn.value;
642  }
643  default :
644  {
645  G4Exception("G4TwistedTubs::DistanceToIn(p)", "GeomSolids0003",
646  FatalException, "Unknown point location!");
647  }
648  } // switch end
649 
650  return kInfinity;
651 }
652 
653 //=====================================================================
654 //* DistanceToOut (p, v) ----------------------------------------------
655 
657  const G4ThreeVector& v,
658  const G4bool calcNorm,
659  G4bool *validNorm,
660  G4ThreeVector *norm ) const
661 {
662  // DistanceToOut (p, v):
663  // Calculate distance to surface of shape from `inside'
664  // along with the v, allowing for tolerance.
665  // The function returns kInfinity if no intersection or
666  // just grazing within tolerance.
667 
668  //
669  // checking last value
670  //
671 
672  G4ThreeVector *tmpp;
673  G4ThreeVector *tmpv;
674  G4double *tmpdist;
676  {
678  }
679  else
680  {
681  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
682  tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
683  tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
684  tmpp->set(p.x(), p.y(), p.z());
685  tmpv->set(v.x(), v.y(), v.z());
686  }
687 
688  //
689  // Calculate DistanceToOut(p,v)
690  //
691 
692  EInside currentside = Inside(p);
693 
694  if (currentside == kOutside)
695  {
696  }
697  else
698  {
699  if (currentside == kSurface)
700  {
701  // particle is just on a boundary.
702  // If the particle is exiting from the volume, return 0.
703  //
705  G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
706  if (normal*v > 0)
707  {
708  if (calcNorm)
709  {
710  *norm = (blockedsurface->GetNormal(p, true));
711  *validNorm = blockedsurface->IsValidNorm();
712  }
713  *tmpdist = 0.;
715  }
716  }
717  }
718 
719  // now, we can take smallest positive distance.
720 
721  // Initialize
722  //
723  G4double distance = kInfinity;
724 
725  // find intersections and choose nearest one.
726  //
727  G4VTwistSurface *surfaces[6];
728  surfaces[0] = fLatterTwisted;
729  surfaces[1] = fFormerTwisted;
730  surfaces[2] = fInnerHype;
731  surfaces[3] = fOuterHype;
732  surfaces[4] = fLowerEndcap;
733  surfaces[5] = fUpperEndcap;
734 
735  G4int i;
736  G4int besti = -1;
738  G4ThreeVector bestxx;
739  for (i=0; i< 6; i++)
740  {
741  G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
742  if (tmpdistance < distance)
743  {
744  distance = tmpdistance;
745  bestxx = xx;
746  besti = i;
747  }
748  }
749 
750  if (calcNorm)
751  {
752  if (besti != -1)
753  {
754  *norm = (surfaces[besti]->GetNormal(p, true));
755  *validNorm = surfaces[besti]->IsValidNorm();
756  }
757  }
758 
759  *tmpdist = distance;
760 
762 }
763 
764 
765 //=====================================================================
766 //* DistanceToOut (p) ----------------------------------------------
767 
769 {
770  // DistanceToOut(p):
771  // Calculate distance to surface of shape from `inside',
772  // allowing for tolerance
773 
774  //
775  // checking last value
776  //
777 
778  G4ThreeVector *tmpp;
779  G4double *tmpdist;
780  if (fLastDistanceToOut.p == p)
781  {
782  return fLastDistanceToOut.value;
783  }
784  else
785  {
786  tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
787  tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
788  tmpp->set(p.x(), p.y(), p.z());
789  }
790 
791  //
792  // Calculate DistanceToOut(p)
793  //
794 
795  EInside currentside = Inside(p);
796 
797  switch (currentside)
798  {
799  case (kOutside) :
800  {
801  }
802  case (kSurface) :
803  {
804  *tmpdist = 0.;
805  return fLastDistanceToOut.value;
806  }
807  case (kInside) :
808  {
809  // Initialize
810  G4double distance = kInfinity;
811 
812  // find intersections and choose nearest one.
813  G4VTwistSurface *surfaces[6];
814  surfaces[0] = fLatterTwisted;
815  surfaces[1] = fFormerTwisted;
816  surfaces[2] = fInnerHype;
817  surfaces[3] = fOuterHype;
818  surfaces[4] = fLowerEndcap;
819  surfaces[5] = fUpperEndcap;
820 
821  G4int i;
823  G4ThreeVector bestxx;
824  for (i=0; i< 6; i++)
825  {
826  G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
827  if (tmpdistance < distance)
828  {
829  distance = tmpdistance;
830  bestxx = xx;
831  }
832  }
833  *tmpdist = distance;
834 
835  return fLastDistanceToOut.value;
836  }
837  default :
838  {
839  G4Exception("G4TwistedTubs::DistanceToOut(p)", "GeomSolids0003",
840  FatalException, "Unknown point location!");
841  }
842  } // switch end
843 
844  return 0;
845 }
846 
847 //=====================================================================
848 //* StreamInfo --------------------------------------------------------
849 
850 std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
851 {
852  //
853  // Stream object contents to an output stream
854  //
855  G4int oldprc = os.precision(16);
856  os << "-----------------------------------------------------------\n"
857  << " *** Dump for solid - " << GetName() << " ***\n"
858  << " ===================================================\n"
859  << " Solid type: G4TwistedTubs\n"
860  << " Parameters: \n"
861  << " -ve end Z : " << fEndZ[0]/mm << " mm \n"
862  << " +ve end Z : " << fEndZ[1]/mm << " mm \n"
863  << " inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
864  << " inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
865  << " outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
866  << " outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
867  << " inner radius (z=0) : " << fInnerRadius/mm << " mm \n"
868  << " outer radius (z=0) : " << fOuterRadius/mm << " mm \n"
869  << " twisted angle : " << fPhiTwist/degree << " degrees \n"
870  << " inner stereo angle : " << fInnerStereo/degree << " degrees \n"
871  << " outer stereo angle : " << fOuterStereo/degree << " degrees \n"
872  << " phi-width of a piece : " << fDPhi/degree << " degrees \n"
873  << "-----------------------------------------------------------\n";
874  os.precision(oldprc);
875 
876  return os;
877 }
878 
879 
880 //=====================================================================
881 //* DiscribeYourselfTo ------------------------------------------------
882 
884 {
885  scene.AddSolid (*this);
886 }
887 
888 //=====================================================================
889 //* GetExtent ---------------------------------------------------------
890 
892 {
893  // Define the sides of the box into which the G4Tubs instance would fit.
894 
895  G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
897  return G4VisExtent( -maxEndOuterRad, maxEndOuterRad,
898  -maxEndOuterRad, maxEndOuterRad,
900 }
901 
902 //=====================================================================
903 //* CreatePolyhedron --------------------------------------------------
904 
906 {
907  // number of meshes
908  //
909  G4double absPhiTwist = std::abs(fPhiTwist);
910  G4double dA = std::max(fDPhi,absPhiTwist);
911  const G4int k =
913  const G4int n =
914  G4int(G4Polyhedron::GetNumberOfRotationSteps() * absPhiTwist / twopi) + 2;
915 
916  const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
917  const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
918 
919  G4Polyhedron *ph=new G4Polyhedron;
920  typedef G4double G4double3[3];
921  typedef G4int G4int4[4];
922  G4double3* xyz = new G4double3[nnodes]; // number of nodes
923  G4int4* faces = new G4int4[nfaces] ; // number of faces
924  fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
925  fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
926  fInnerHype->GetFacets(k,n,xyz,faces,2) ;
927  fFormerTwisted->GetFacets(k,n,xyz,faces,3) ;
928  fOuterHype->GetFacets(k,n,xyz,faces,4) ;
929  fLatterTwisted->GetFacets(k,n,xyz,faces,5) ;
930 
931  ph->createPolyhedron(nnodes,nfaces,xyz,faces);
932 
933  delete[] xyz;
934  delete[] faces;
935 
936  return ph;
937 }
938 
939 //=====================================================================
940 //* GetPolyhedron -----------------------------------------------------
941 
943 {
944  if (!fpPolyhedron ||
948  {
949  G4AutoLock l(&polyhedronMutex);
950  delete fpPolyhedron;
952  fRebuildPolyhedron = false;
953  l.unlock();
954  }
955  return fpPolyhedron;
956 }
957 
958 //=====================================================================
959 //* CreateSurfaces ----------------------------------------------------
960 
962 {
963  // create 6 surfaces of TwistedTub
964 
965  G4ThreeVector x0(0, 0, fEndZ[0]);
966  G4ThreeVector n (0, 0, -1);
967 
968  fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap",
970  fDPhi, fEndPhi, fEndZ, -1) ;
971 
972  fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap",
974  fDPhi, fEndPhi, fEndZ, 1) ;
975 
976  G4RotationMatrix rotHalfDPhi;
977  rotHalfDPhi.rotateZ(0.5*fDPhi);
978 
979  fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
981  fDPhi, fEndPhi, fEndZ,
983  1 ) ;
984  fFormerTwisted = new G4TwistTubsSide("FormerTwisted",
986  fDPhi, fEndPhi, fEndZ,
988  -1 ) ;
989 
990  fInnerHype = new G4TwistTubsHypeSide("InnerHype",
992  fDPhi, fEndPhi, fEndZ,
995  fOuterHype = new G4TwistTubsHypeSide("OuterHype",
997  fDPhi, fEndPhi, fEndZ,
1000 
1001 
1002  // set neighbour surfaces
1003  //
1016 }
1017 
1018 
1019 //=====================================================================
1020 //* GetEntityType -----------------------------------------------------
1021 
1023 {
1024  return G4String("G4TwistedTubs");
1025 }
1026 
1027 //=====================================================================
1028 //* Clone -------------------------------------------------------------
1029 
1031 {
1032  return new G4TwistedTubs(*this);
1033 }
1034 
1035 //=====================================================================
1036 //* GetCubicVolume ----------------------------------------------------
1037 
1039 {
1040  if(fCubicVolume != 0.) {;}
1043  return fCubicVolume;
1044 }
1045 
1046 //=====================================================================
1047 //* GetSurfaceArea ----------------------------------------------------
1048 
1050 {
1051  if(fSurfaceArea != 0.) {;}
1053  return fSurfaceArea;
1054 }
1055 
1056 //=====================================================================
1057 //* GetPointOnSurface -------------------------------------------------
1058 
1060 {
1061 
1063  G4double phi , phimin, phimax ;
1064  G4double x , xmin, xmax ;
1065  G4double r , rmin, rmax ;
1066 
1073 
1074  G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1075 
1076  if(chose < a1)
1077  {
1078 
1079  phimin = fOuterHype->GetBoundaryMin(z) ;
1080  phimax = fOuterHype->GetBoundaryMax(z) ;
1081  phi = G4RandFlat::shoot(phimin,phimax) ;
1082 
1083  return fOuterHype->SurfacePoint(phi,z,true) ;
1084 
1085  }
1086  else if ( (chose >= a1) && (chose < a1 + a2 ) )
1087  {
1088 
1089  phimin = fInnerHype->GetBoundaryMin(z) ;
1090  phimax = fInnerHype->GetBoundaryMax(z) ;
1091  phi = G4RandFlat::shoot(phimin,phimax) ;
1092 
1093  return fInnerHype->SurfacePoint(phi,z,true) ;
1094 
1095  }
1096  else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1097  {
1098 
1099  xmin = fLatterTwisted->GetBoundaryMin(z) ;
1100  xmax = fLatterTwisted->GetBoundaryMax(z) ;
1101  x = G4RandFlat::shoot(xmin,xmax) ;
1102 
1103  return fLatterTwisted->SurfacePoint(x,z,true) ;
1104 
1105  }
1106  else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1107  {
1108 
1109  xmin = fFormerTwisted->GetBoundaryMin(z) ;
1110  xmax = fFormerTwisted->GetBoundaryMax(z) ;
1111  x = G4RandFlat::shoot(xmin,xmax) ;
1112 
1113  return fFormerTwisted->SurfacePoint(x,z,true) ;
1114  }
1115  else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
1116  {
1117  rmin = GetEndInnerRadius(0) ;
1118  rmax = GetEndOuterRadius(0) ;
1119  r = std::sqrt(G4RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
1120 
1121  phimin = fLowerEndcap->GetBoundaryMin(r) ;
1122  phimax = fLowerEndcap->GetBoundaryMax(r) ;
1123  phi = G4RandFlat::shoot(phimin,phimax) ;
1124 
1125  return fLowerEndcap->SurfacePoint(phi,r,true) ;
1126  }
1127  else
1128  {
1129  rmin = GetEndInnerRadius(1) ;
1130  rmax = GetEndOuterRadius(1) ;
1131  r = rmin + (rmax-rmin)*std::sqrt(G4RandFlat::shoot());
1132 
1133  phimin = fUpperEndcap->GetBoundaryMin(r) ;
1134  phimax = fUpperEndcap->GetBoundaryMax(r) ;
1135  phi = G4RandFlat::shoot(phimin,phimax) ;
1136 
1137  return fUpperEndcap->SurfacePoint(phi,r,true) ;
1138  }
1139 }
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
Float_t x
Definition: compare.C:6
G4double fTanOuterStereo
void set(double x, double y, double z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double fEndPhi[2]
void ComputeDimensions(G4VPVParameterisation *, const G4int, const G4VPhysicalVolume *)
G4Polyhedron * fpPolyhedron
G4double GetEndOuterRadius() const
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
LastValue fLastDistanceToOut
G4double fPhiTwist
Double_t xx
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4Polyhedron * GetPolyhedron() const
G4double fOuterRadius
virtual void AddSolid(const G4Box &)=0
G4double GetSurfaceArea()
static const G4double kInfinity
Definition: geomdefs.hh:42
G4VTwistSurface * fInnerHype
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4double fInnerStereo
G4double fTanOuterStereo2
static constexpr double mm
Definition: G4SIunits.hh:115
#define G4endl
Definition: G4ios.hh:61
const char * p
Definition: xmltok.h:285
G4bool fRebuildPolyhedron
Double_t z
G4TwistedTubs(const G4String &pname, G4double twistedangle, G4double endinnerrad, G4double endouterrad, G4double halfzlen, G4double dphi)
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
virtual G4double GetSurfaceArea()
Definition: G4VSolid.cc:250
G4VTwistSurface * fLatterTwisted
double z() const
G4double fTanInnerStereo2
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4double GetRadialTolerance() const
G4double fTanInnerStereo
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcnorm=G4bool(false), G4bool *validnorm=0, G4ThreeVector *n=0) const
G4bool IsValidNorm() const
G4VTwistSurface * fLowerEndcap
G4GeometryType GetEntityType() const
virtual G4double GetSurfaceArea()=0
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal)=0
G4double fCubicVolume
G4VisExtent GetExtent() const
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:110
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4VTwistSurface * fFormerTwisted
LastState fLastInside
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
G4double fOuterStereo
virtual G4double DistanceTo(const G4ThreeVector &gp, G4ThreeVector &gxx)
G4VSolid * Clone() const
static constexpr double twopi
Definition: G4SIunits.hh:76
LastVector fLastNormal
G4double GetCubicVolume()
virtual G4double DistanceToIn(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
EInside Inside(const G4ThreeVector &p) const
virtual G4double GetBoundaryMax(G4double)=0
virtual G4double GetBoundaryMin(G4double)=0
G4VTwistSurface * fUpperEndcap
G4TwistedTubs & operator=(const G4TwistedTubs &rhs)
LastValue fLastDistanceToIn
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double fEndOuterRadius[2]
static G4int GetNumberOfRotationSteps()
G4Polyhedron * CreatePolyhedron() const
G4double fInnerRadius2
G4double GetEndInnerRadius() const
G4VTwistSurface ** surface
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
void SetNeighbours(G4VTwistSurface *axis0min, G4VTwistSurface *axis1min, G4VTwistSurface *axis0max, G4VTwistSurface *axis1max)
std::ostream & StreamInfo(std::ostream &os) const
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
int G4int
Definition: G4Types.hh:78
EInside
Definition: geomdefs.hh:58
LastValueWithDoubleVector fLastDistanceToInWithV
G4VTwistSurface * fOuterHype
EAxis
Definition: geomdefs.hh:54
double getRho() const
G4ThreeVector GetPointOnSurface() const
G4String GetName() const
void DumpInfo() const
Float_t norm
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
#define DBL_MIN
Definition: templates.hh:75
virtual ~G4TwistedTubs()
void SetFields(G4double phitwist, G4double innerrad, G4double outerrad, G4double negativeEndz, G4double positiveEndz)
HepRotation & rotateZ(double delta)
Definition: Rotation.cc:92
G4double fEndInnerRadius[2]
static constexpr double degree
Definition: G4SIunits.hh:144
double x() const
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)=0
static G4GeometryTolerance * GetInstance()
Char_t n[5]
T sqr(const T &x)
Definition: templates.hh:145
G4double fInnerRadius
G4double fZHalfLength
double y() const
G4double fSurfaceArea
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
virtual G4double DistanceToOut(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
G4double fEndZ[2]
G4double fEndZ2[2]
void CreateSurfaces()
LastValueWithDoubleVector fLastDistanceToOutWithV
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)=0
G4double fOuterRadius2
std::mutex G4Mutex
Definition: G4Threading.hh:84