| Geant4 Cross Reference |
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 // Implementation for G4UOrb wrapper class
27 //
28 // 30.10.13 G.Cosmo, CERN/PH
29 // --------------------------------------------------------------------
30
31 #include "G4Orb.hh"
32 #include "G4UOrb.hh"
33
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
35
36 #include "G4TwoVector.hh"
37 #include "G4AffineTransform.hh"
38 #include "G4GeometryTolerance.hh"
39 #include "G4BoundingEnvelope.hh"
40
41 #include "G4VPVParameterisation.hh"
42 #include "G4PhysicalConstants.hh"
43
44 using namespace CLHEP;
45
46 ////////////////////////////////////////////////////////////////////////
47 //
48 // constructor - check positive radius
49 //
50
51 G4UOrb::G4UOrb( const G4String& pName, G4double pRmax )
52 : Base_t(pName, pRmax)
53 {
54 }
55
56 ///////////////////////////////////////////////////////////////////////
57 //
58 // Fake default constructor - sets only member data and allocates memory
59 // for usage restricted to object persistency.
60 //
61 G4UOrb::G4UOrb( __void__& a )
62 : Base_t(a)
63 {
64 }
65
66 /////////////////////////////////////////////////////////////////////
67 //
68 // Destructor
69
70 G4UOrb::~G4UOrb() = default;
71
72 //////////////////////////////////////////////////////////////////////////
73 //
74 // Copy constructor
75
76 G4UOrb::G4UOrb(const G4UOrb& rhs)
77 : Base_t(rhs)
78 {
79 }
80
81 //////////////////////////////////////////////////////////////////////////
82 //
83 // Assignment operator
84
85 G4UOrb& G4UOrb::operator = (const G4UOrb& rhs)
86 {
87 // Check assignment to self
88 //
89 if (this == &rhs) { return *this; }
90
91 // Copy base class data
92 //
93 Base_t::operator=(rhs);
94
95 return *this;
96 }
97
98 //////////////////////////////////////////////////////////////////////////
99 //
100 // Accessors & modifiers
101
102 G4double G4UOrb::GetRadius() const
103 {
104 return Base_t::GetRadius();
105 }
106
107 void G4UOrb::SetRadius(G4double newRmax)
108 {
109 Base_t::SetRadius(newRmax);
110 fRebuildPolyhedron = true;
111 }
112
113 G4double G4UOrb::GetRadialTolerance() const
114 {
115 return Base_t::GetRadialTolerance();
116 }
117
118 //////////////////////////////////////////////////////////////////////////
119 //
120 // Dispatch to parameterisation for replication mechanism dimension
121 // computation & modification.
122
123 void G4UOrb::ComputeDimensions( G4VPVParameterisation* p,
124 const G4int n,
125 const G4VPhysicalVolume* pRep )
126 {
127 p->ComputeDimensions(*(G4Orb*)this,n,pRep);
128 }
129
130 //////////////////////////////////////////////////////////////////////////
131 //
132 // Make a clone of the object
133
134 G4VSolid* G4UOrb::Clone() const
135 {
136 return new G4UOrb(*this);
137 }
138
139 //////////////////////////////////////////////////////////////////////////
140 //
141 // Get bounding box
142
143 void G4UOrb::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
144 {
145 G4double radius = GetRadius();
146 pMin.set(-radius,-radius,-radius);
147 pMax.set( radius, radius, radius);
148
149 // Check correctness of the bounding box
150 //
151 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
152 {
153 std::ostringstream message;
154 message << "Bad bounding box (min >= max) for solid: "
155 << GetName() << " !"
156 << "\npMin = " << pMin
157 << "\npMax = " << pMax;
158 G4Exception("G4UOrb::BoundingLimits()", "GeomMgt0001",
159 JustWarning, message);
160 StreamInfo(G4cout);
161 }
162 }
163
164 //////////////////////////////////////////////////////////////////////////
165 //
166 // Calculate extent under transform and specified limit
167
168 G4bool
169 G4UOrb::CalculateExtent(const EAxis pAxis,
170 const G4VoxelLimits& pVoxelLimit,
171 const G4AffineTransform& pTransform,
172 G4double& pMin, G4double& pMax) const
173 {
174 G4ThreeVector bmin, bmax;
175 G4bool exist;
176
177 // Get bounding box
178 BoundingLimits(bmin,bmax);
179
180 // Check bounding box
181 G4BoundingEnvelope bbox(bmin,bmax);
182 #ifdef G4BBOX_EXTENT
183 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
184 #endif
185 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
186 {
187 return exist = pMin < pMax;
188 }
189
190 // Find bounding envelope and calculate extent
191 //
192 static const G4int NTHETA = 8; // number of steps along Theta
193 static const G4int NPHI = 16; // number of steps along Phi
194 static const G4double sinHalfTheta = std::sin(halfpi/NTHETA);
195 static const G4double cosHalfTheta = std::cos(halfpi/NTHETA);
196 static const G4double sinHalfPhi = std::sin(pi/NPHI);
197 static const G4double cosHalfPhi = std::cos(pi/NPHI);
198 static const G4double sinStepTheta = 2.*sinHalfTheta*cosHalfTheta;
199 static const G4double cosStepTheta = 1. - 2.*sinHalfTheta*sinHalfTheta;
200 static const G4double sinStepPhi = 2.*sinHalfPhi*cosHalfPhi;
201 static const G4double cosStepPhi = 1. - 2.*sinHalfPhi*sinHalfPhi;
202
203 G4double radius = GetRadius();
204 G4double rtheta = radius/cosHalfTheta;
205 G4double rphi = rtheta/cosHalfPhi;
206
207 // set reference circle
208 G4TwoVector xy[NPHI];
209 G4double sinCurPhi = sinHalfPhi;
210 G4double cosCurPhi = cosHalfPhi;
211 for (auto & k : xy)
212 {
213 k.set(cosCurPhi,sinCurPhi);
214 G4double sinTmpPhi = sinCurPhi;
215 sinCurPhi = sinCurPhi*cosStepPhi + cosCurPhi*sinStepPhi;
216 cosCurPhi = cosCurPhi*cosStepPhi - sinTmpPhi*sinStepPhi;
217 }
218
219 // set bounding circles
220 G4ThreeVectorList circles[NTHETA];
221 for (auto & circle : circles) circle.resize(NPHI);
222
223 G4double sinCurTheta = sinHalfTheta;
224 G4double cosCurTheta = cosHalfTheta;
225 for (auto & circle : circles)
226 {
227 G4double z = rtheta*cosCurTheta;
228 G4double rho = rphi*sinCurTheta;
229 for (G4int k=0; k<NPHI; ++k)
230 {
231 circle[k].set(rho*xy[k].x(),rho*xy[k].y(),z);
232 }
233 G4double sinTmpTheta = sinCurTheta;
234 sinCurTheta = sinCurTheta*cosStepTheta + cosCurTheta*sinStepTheta;
235 cosCurTheta = cosCurTheta*cosStepTheta - sinTmpTheta*sinStepTheta;
236 }
237
238 // set envelope and calculate extent
239 std::vector<const G4ThreeVectorList *> polygons;
240 polygons.resize(NTHETA);
241 for (G4int i=0; i<NTHETA; ++i) polygons[i] = &circles[i];
242
243 G4BoundingEnvelope benv(bmin,bmax,polygons);
244 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
245 return exist;
246 }
247
248 //////////////////////////////////////////////////////////////////////////
249 //
250 // Create polyhedron for visualization
251
252 G4Polyhedron* G4UOrb::CreatePolyhedron() const
253 {
254 return new G4PolyhedronSphere(0., GetRadius(), 0., twopi, 0., pi);
255 }
256
257 #endif // G4GEOM_USE_USOLIDS
258