| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4TwistedTubs implementation <<
27 // 23 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu << 24 // $Id: G4TwistedTubs.cc,v 1.12 2005/04/04 11:56:59 gcosmo Exp $
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), << 25 // GEANT4 tag $Name: geant4-07-01 $
30 // from original version in Jupi << 26 //
>> 27 //
>> 28 // --------------------------------------------------------------------
>> 29 // GEANT 4 class source file
>> 30 //
>> 31 //
>> 32 // G4TwistedSurface.cc
>> 33 //
>> 34 // Author:
>> 35 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
>> 36 //
>> 37 // History:
>> 38 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
>> 39 // from original version in Jupiter-2.5.02 application.
31 // ------------------------------------------- 40 // --------------------------------------------------------------------
32 41
33 #include "G4TwistedTubs.hh" 42 #include "G4TwistedTubs.hh"
34 43
35 #include "G4GeomTools.hh" <<
36 #include "G4PhysicalConstants.hh" <<
37 #include "G4SystemOfUnits.hh" <<
38 #include "G4GeometryTolerance.hh" <<
39 #include "G4VoxelLimits.hh" 44 #include "G4VoxelLimits.hh"
40 #include "G4AffineTransform.hh" 45 #include "G4AffineTransform.hh"
41 #include "G4BoundingEnvelope.hh" << 46 #include "G4SolidExtentList.hh"
42 #include "G4ClippablePolygon.hh" 47 #include "G4ClippablePolygon.hh"
43 #include "G4VPVParameterisation.hh" 48 #include "G4VPVParameterisation.hh"
44 #include "meshdefs.hh" 49 #include "meshdefs.hh"
45 50
46 #include "G4VGraphicsScene.hh" 51 #include "G4VGraphicsScene.hh"
47 #include "G4Polyhedron.hh" 52 #include "G4Polyhedron.hh"
48 #include "G4VisExtent.hh" 53 #include "G4VisExtent.hh"
49 << 54 #include "G4NURBS.hh"
50 #include "Randomize.hh" << 55 #include "G4NURBStube.hh"
51 << 56 #include "G4NURBScylinder.hh"
52 #include "G4AutoLock.hh" << 57 #include "G4NURBStubesector.hh"
53 <<
54 namespace <<
55 { <<
56 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE <<
57 } <<
58 <<
59 58
60 //============================================ 59 //=====================================================================
61 //* constructors ----------------------------- 60 //* constructors ------------------------------------------------------
62 61
63 G4TwistedTubs::G4TwistedTubs(const G4String& p << 62 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
64 G4double t 63 G4double twistedangle,
65 G4double e 64 G4double endinnerrad,
66 G4double e 65 G4double endouterrad,
67 G4double h 66 G4double halfzlen,
68 G4double d 67 G4double dphi)
69 : G4VSolid(pname), fDPhi(dphi), 68 : G4VSolid(pname), fDPhi(dphi),
70 fLowerEndcap(nullptr), fUpperEndcap(nullp << 69 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
71 fFormerTwisted(nullptr), fInnerHype(nullp << 70 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
72 { << 71 fCubicVolume(0.), fpPolyhedron(0)
73 if (endinnerrad < DBL_MIN) << 72 {
74 { << 73 if (endinnerrad < DBL_MIN) {
75 G4Exception("G4TwistedTubs::G4TwistedTub << 74 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
76 FatalErrorInArgument, "Inval << 75 FatalException, "Invalid end-inner-radius!");
77 } 76 }
78 77
79 G4double sinhalftwist = std::sin(0.5 * twis 78 G4double sinhalftwist = std::sin(0.5 * twistedangle);
80 79
81 G4double endinnerradX = endinnerrad * sinha 80 G4double endinnerradX = endinnerrad * sinhalftwist;
82 G4double innerrad = std::sqrt( endinner 81 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
83 - endinnerrad 82 - endinnerradX * endinnerradX );
84 83
85 G4double endouterradX = endouterrad * sinha 84 G4double endouterradX = endouterrad * sinhalftwist;
86 G4double outerrad = std::sqrt( endouter 85 G4double outerrad = std::sqrt( endouterrad * endouterrad
87 - endouterrad 86 - endouterradX * endouterradX );
88 87
89 // temporary treatment!! 88 // temporary treatment!!
90 SetFields(twistedangle, innerrad, outerrad, 89 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
91 CreateSurfaces(); 90 CreateSurfaces();
92 } 91 }
93 92
94 G4TwistedTubs::G4TwistedTubs(const G4String& p << 93 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
95 G4double t 94 G4double twistedangle,
96 G4double e 95 G4double endinnerrad,
97 G4double e 96 G4double endouterrad,
98 G4double h 97 G4double halfzlen,
99 G4int n 98 G4int nseg,
100 G4double t 99 G4double totphi)
101 : G4VSolid(pname), 100 : G4VSolid(pname),
102 fLowerEndcap(nullptr), fUpperEndcap(nullp << 101 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
103 fFormerTwisted(nullptr), fInnerHype(nullp << 102 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
>> 103 fCubicVolume(0.), fpPolyhedron(0)
104 { 104 {
105 105
106 if (nseg == 0) << 106 if (!nseg) G4cerr << "ERROR - G4TwistedTubs::G4TwistedTubs()" << G4endl
107 { << 107 << " Invalid nseg. nseg = " << nseg << G4endl;
108 std::ostringstream message; << 108 if (totphi == DBL_MIN || endinnerrad < DBL_MIN) {
109 message << "Invalid number of segments." << 109 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
110 << " nseg = " << nseg; << 110 FatalException, "Invalid total-phi or end-inner-radius!");
111 G4Exception("G4TwistedTubs::G4TwistedTub << 111 }
112 FatalErrorInArgument, messag <<
113 } <<
114 if (totphi == DBL_MIN || endinnerrad < DBL_ <<
115 { <<
116 G4Exception("G4TwistedTubs::G4TwistedTub <<
117 FatalErrorInArgument, "Invalid <<
118 } <<
119 112
120 G4double sinhalftwist = std::sin(0.5 * twis 113 G4double sinhalftwist = std::sin(0.5 * twistedangle);
121 114
122 G4double endinnerradX = endinnerrad * sinha 115 G4double endinnerradX = endinnerrad * sinhalftwist;
123 G4double innerrad = std::sqrt( endinner 116 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
124 - endinnerrad 117 - endinnerradX * endinnerradX );
125 118
126 G4double endouterradX = endouterrad * sinha 119 G4double endouterradX = endouterrad * sinhalftwist;
127 G4double outerrad = std::sqrt( endouter 120 G4double outerrad = std::sqrt( endouterrad * endouterrad
128 - endouterrad 121 - endouterradX * endouterradX );
129 122
130 // temporary treatment!! 123 // temporary treatment!!
131 fDPhi = totphi / nseg; 124 fDPhi = totphi / nseg;
132 SetFields(twistedangle, innerrad, outerrad, 125 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
133 CreateSurfaces(); 126 CreateSurfaces();
134 } 127 }
135 128
136 G4TwistedTubs::G4TwistedTubs(const G4String& p << 129 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
137 G4double t 130 G4double twistedangle,
138 G4double i 131 G4double innerrad,
139 G4double o 132 G4double outerrad,
140 G4double n 133 G4double negativeEndz,
141 G4double p 134 G4double positiveEndz,
142 G4double d 135 G4double dphi)
143 : G4VSolid(pname), fDPhi(dphi), 136 : G4VSolid(pname), fDPhi(dphi),
144 fLowerEndcap(nullptr), fUpperEndcap(nullp << 137 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
145 fFormerTwisted(nullptr), fInnerHype(nullp << 138 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
146 { << 139 fCubicVolume(0.), fpPolyhedron(0)
147 if (innerrad < DBL_MIN) << 140 {
148 { << 141 if (innerrad < DBL_MIN) {
149 G4Exception("G4TwistedTubs::G4TwistedTub << 142 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
150 FatalErrorInArgument, "Inval << 143 FatalException, "Invalid end-inner-radius!");
151 } 144 }
152 145
153 SetFields(twistedangle, innerrad, outerrad, 146 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
154 CreateSurfaces(); 147 CreateSurfaces();
155 } 148 }
156 149
157 G4TwistedTubs::G4TwistedTubs(const G4String& p << 150 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
158 G4double t 151 G4double twistedangle,
159 G4double i 152 G4double innerrad,
160 G4double o 153 G4double outerrad,
161 G4double n 154 G4double negativeEndz,
162 G4double p 155 G4double positiveEndz,
163 G4int n 156 G4int nseg,
164 G4double t 157 G4double totphi)
165 : G4VSolid(pname), 158 : G4VSolid(pname),
166 fLowerEndcap(nullptr), fUpperEndcap(nullp << 159 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
167 fFormerTwisted(nullptr), fInnerHype(nullp << 160 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
168 { << 161 fCubicVolume(0.), fpPolyhedron(0)
169 if (nseg == 0) << 162 {
170 { << 163 if (!nseg) G4cerr << "ERROR - G4TwistedTubs::G4TwistedTubs()" << G4endl
171 std::ostringstream message; << 164 << " Invalid nseg. nseg = " << nseg << G4endl;
172 message << "Invalid number of segments." << 165 if (totphi == DBL_MIN || innerrad < DBL_MIN) {
173 << " nseg = " << nseg; << 166 G4Exception("G4TwistedTubs::G4TwistedTubs()", "InvalidSetup",
174 G4Exception("G4TwistedTubs::G4TwistedTub << 167 FatalException, "Invalid total-phi or end-inner-radius!");
175 FatalErrorInArgument, messag <<
176 } <<
177 if (totphi == DBL_MIN || innerrad < DBL_MIN <<
178 { <<
179 G4Exception("G4TwistedTubs::G4TwistedTub <<
180 FatalErrorInArgument, "Invalid <<
181 } 168 }
182 169
183 fDPhi = totphi / nseg; 170 fDPhi = totphi / nseg;
184 SetFields(twistedangle, innerrad, outerrad, 171 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
185 CreateSurfaces(); 172 CreateSurfaces();
186 } 173 }
187 174
188 //============================================ 175 //=====================================================================
189 //* Fake default constructor ----------------- <<
190 <<
191 G4TwistedTubs::G4TwistedTubs( __void__& a ) <<
192 : G4VSolid(a), <<
193 fLowerEndcap(nullptr), fUpperEndcap(nullpt <<
194 fLatterTwisted(nullptr), fFormerTwisted(nu <<
195 fInnerHype(nullptr), fOuterHype(nullptr) <<
196 { <<
197 } <<
198 <<
199 //============================================ <<
200 //* destructor ------------------------------- 176 //* destructor --------------------------------------------------------
201 177
202 G4TwistedTubs::~G4TwistedTubs() 178 G4TwistedTubs::~G4TwistedTubs()
203 { 179 {
204 delete fLowerEndcap; << 180 if (fLowerEndcap) delete fLowerEndcap;
205 delete fUpperEndcap; << 181 if (fUpperEndcap) delete fUpperEndcap;
206 delete fLatterTwisted; << 182 if (fLatterTwisted) delete fLatterTwisted;
207 delete fFormerTwisted; << 183 if (fFormerTwisted) delete fFormerTwisted;
208 delete fInnerHype; << 184 if (fInnerHype) delete fInnerHype;
209 delete fOuterHype; << 185 if (fOuterHype) delete fOuterHype;
210 delete fpPolyhedron; fpPolyhedron = nullptr << 186 if (fpPolyhedron) delete fpPolyhedron;
211 } <<
212 <<
213 //============================================ <<
214 //* Copy constructor ------------------------- <<
215 <<
216 G4TwistedTubs::G4TwistedTubs(const G4TwistedTu <<
217 : G4VSolid(rhs), fPhiTwist(rhs.fPhiTwist), <<
218 fInnerRadius(rhs.fInnerRadius), fOuterRadi <<
219 fDPhi(rhs.fDPhi), fZHalfLength(rhs.fZHalfL <<
220 fInnerStereo(rhs.fInnerStereo), fOuterSter <<
221 fTanInnerStereo(rhs.fTanInnerStereo), fTan <<
222 fKappa(rhs.fKappa), fInnerRadius2(rhs.fInn <<
223 fOuterRadius2(rhs.fOuterRadius2), fTanInne <<
224 fTanOuterStereo2(rhs.fTanOuterStereo2), <<
225 fLowerEndcap(nullptr), fUpperEndcap(nullpt <<
226 fInnerHype(nullptr), fOuterHype(nullptr), <<
227 fCubicVolume(rhs.fCubicVolume), fSurfaceAr <<
228 { <<
229 for (auto i=0; i<2; ++i) <<
230 { <<
231 fEndZ[i] = rhs.fEndZ[i]; <<
232 fEndInnerRadius[i] = rhs.fEndInnerRadius[i <<
233 fEndOuterRadius[i] = rhs.fEndOuterRadius[i <<
234 fEndPhi[i] = rhs.fEndPhi[i]; <<
235 fEndZ2[i] = rhs.fEndZ2[i]; <<
236 } <<
237 CreateSurfaces(); <<
238 } <<
239 <<
240 <<
241 //============================================ <<
242 //* Assignment operator ---------------------- <<
243 <<
244 G4TwistedTubs& G4TwistedTubs::operator = (cons <<
245 { <<
246 // Check assignment to self <<
247 // <<
248 if (this == &rhs) { return *this; } <<
249 <<
250 // Copy base class data <<
251 // <<
252 G4VSolid::operator=(rhs); <<
253 <<
254 // Copy data <<
255 // <<
256 fPhiTwist= rhs.fPhiTwist; <<
257 fInnerRadius= rhs.fInnerRadius; fOuterRadiu <<
258 fDPhi= rhs.fDPhi; fZHalfLength= rhs.fZHalfL <<
259 fInnerStereo= rhs.fInnerStereo; fOuterStere <<
260 fTanInnerStereo= rhs.fTanInnerStereo; fTanO <<
261 fKappa= rhs.fKappa; fInnerRadius2= rhs.fInn <<
262 fOuterRadius2= rhs.fOuterRadius2; fTanInner <<
263 fTanOuterStereo2= rhs.fTanOuterStereo2; <<
264 fLowerEndcap= fUpperEndcap= fLatterTwisted= <<
265 fInnerHype= fOuterHype= nullptr; <<
266 fCubicVolume= rhs.fCubicVolume; fSurfaceAre <<
267 <<
268 for (auto i=0; i<2; ++i) <<
269 { <<
270 fEndZ[i] = rhs.fEndZ[i]; <<
271 fEndInnerRadius[i] = rhs.fEndInnerRadius[ <<
272 fEndOuterRadius[i] = rhs.fEndOuterRadius[ <<
273 fEndPhi[i] = rhs.fEndPhi[i]; <<
274 fEndZ2[i] = rhs.fEndZ2[i]; <<
275 } <<
276 <<
277 CreateSurfaces(); <<
278 fRebuildPolyhedron = false; <<
279 delete fpPolyhedron; fpPolyhedron = nullptr <<
280 <<
281 return *this; <<
282 } 187 }
283 188
284 //============================================ 189 //=====================================================================
285 //* ComputeDimensions ------------------------ 190 //* ComputeDimensions -------------------------------------------------
286 191
287 void G4TwistedTubs::ComputeDimensions(G4VPVPar 192 void G4TwistedTubs::ComputeDimensions(G4VPVParameterisation* /* p */ ,
288 const G4 193 const G4int /* n */ ,
289 const G4 194 const G4VPhysicalVolume* /* pRep */ )
290 { 195 {
291 G4Exception("G4TwistedTubs::ComputeDimension 196 G4Exception("G4TwistedTubs::ComputeDimensions()",
292 "GeomSolids0001", FatalException << 197 "NotSupported", FatalException,
293 "G4TwistedTubs does not support 198 "G4TwistedTubs does not support Parameterisation.");
294 } 199 }
295 200
>> 201
296 //============================================ 202 //=====================================================================
297 //* BoundingLimits --------------------------- << 203 //* CalculateExtent ---------------------------------------------------
298 204
299 void G4TwistedTubs::BoundingLimits(G4ThreeVect << 205 G4bool G4TwistedTubs::CalculateExtent( const EAxis axis,
300 G4ThreeVect << 206 const G4VoxelLimits &voxelLimit,
>> 207 const G4AffineTransform &transform,
>> 208 G4double &min,
>> 209 G4double &max ) const
301 { 210 {
302 // Find bounding tube << 211
303 G4double rmin = GetInnerRadius(); << 212 G4SolidExtentList extentList( axis, voxelLimit );
304 G4double rmax = GetEndOuterRadius(); << 213 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
305 << 214 fEndOuterRadius[0] : fEndOuterRadius[1]);
306 G4double zmin = std::min(GetEndZ(0), GetEndZ << 215 G4double maxEndInnerRad = (fEndInnerRadius[0] > fEndInnerRadius[1] ?
307 G4double zmax = std::max(GetEndZ(0), GetEndZ << 216 fEndInnerRadius[0] : fEndInnerRadius[1]);
308 << 217 G4double maxphi = (std::fabs(fEndPhi[0]) > std::fabs(fEndPhi[1]) ?
309 G4double dphi = 0.5*GetDPhi(); << 218 std::fabs(fEndPhi[0]) : std::fabs(fEndPhi[1]));
310 G4double sphi = std::min(GetEndPhi(0), GetEn << 219
311 G4double ephi = std::max(GetEndPhi(0), GetEn << 220 //
312 G4double totalphi = ephi - sphi; << 221 // Choose phi size of our segment(s) based on constants as
313 << 222 // defined in meshdefs.hh
314 // Find bounding box << 223 //
315 if (dphi <= 0 || totalphi >= CLHEP::twopi) << 224 // G4int numPhi = kMaxMeshSections;
316 { << 225 G4double sigPhi = 2*maxphi + fDPhi;
317 pMin.set(-rmax,-rmax, zmin); << 226 G4double rFudge = 1.0/std::cos(0.5*sigPhi);
318 pMax.set( rmax, rmax, zmax); << 227 G4double fudgeEndOuterRad = rFudge * maxEndOuterRad;
>> 228
>> 229 //
>> 230 // We work around in phi building polygons along the way.
>> 231 // As a reasonable compromise between accuracy and
>> 232 // complexity (=cpu time), the following facets are chosen:
>> 233 //
>> 234 // 1. If fOuterRadius/maxEndOuterRad > 0.95, approximate
>> 235 // the outer surface as a cylinder, and use one
>> 236 // rectangular polygon (0-1) to build its mesh.
>> 237 //
>> 238 // Otherwise, use two trapazoidal polygons that
>> 239 // meet at z = 0 (0-4-1)
>> 240 //
>> 241 // 2. If there is no inner surface, then use one
>> 242 // polygon for each entire endcap. (0) and (1)
>> 243 //
>> 244 // Otherwise, use a trapazoidal polygon for each
>> 245 // phi segment of each endcap. (0-2) and (1-3)
>> 246 //
>> 247 // 3. For the inner surface, if fInnerRadius/maxEndInnerRad > 0.95,
>> 248 // approximate the inner surface as a cylinder of
>> 249 // radius fInnerRadius and use one rectangular polygon
>> 250 // to build each phi segment of its mesh. (2-3)
>> 251 //
>> 252 // Otherwise, use one rectangular polygon centered
>> 253 // at z = 0 (5-6) and two connecting trapazoidal polygons
>> 254 // for each phi segment (2-5) and (3-6).
>> 255 //
>> 256
>> 257 G4bool splitOuter = (fOuterRadius/maxEndOuterRad < 0.95);
>> 258 G4bool splitInner = (fInnerRadius/maxEndInnerRad < 0.95);
>> 259
>> 260 //
>> 261 // Vertex assignments (v and w arrays)
>> 262 // [0] and [1] are mandatory
>> 263 // the rest are optional
>> 264 //
>> 265 // + -
>> 266 // [0]------[4]------[1] <--- outer radius
>> 267 // | |
>> 268 // | |
>> 269 // [2]---[5]---[6]---[3] <--- inner radius
>> 270 //
>> 271
>> 272 G4ClippablePolygon endPoly1, endPoly2;
>> 273
>> 274 G4double phimax = maxphi + 0.5*fDPhi;
>> 275 G4double phimin = - phimax;
>> 276
>> 277 G4ThreeVector v0, v1, v2, v3, v4, v5, v6; // -ve phi verticies for polygon
>> 278 G4ThreeVector w0, w1, w2, w3, w4, w5, w6; // +ve phi verticies for polygon
>> 279
>> 280 //
>> 281 // decide verticies of -ve phi boundary
>> 282 //
>> 283
>> 284 G4double cosPhi = std::cos(phimin);
>> 285 G4double sinPhi = std::sin(phimin);
>> 286
>> 287 // Outer hyperbolic surface
>> 288
>> 289 v0 = transform.TransformPoint(
>> 290 G4ThreeVector(fudgeEndOuterRad * cosPhi, fudgeEndOuterRad * sinPhi,
>> 291 + fZHalfLength));
>> 292 v1 = transform.TransformPoint(
>> 293 G4ThreeVector(fudgeEndOuterRad * cosPhi, fudgeEndOuterRad * sinPhi,
>> 294 - fZHalfLength));
>> 295 if (splitOuter) {
>> 296 v4 = transform.TransformPoint(
>> 297 G4ThreeVector(fudgeEndOuterRad * cosPhi, fudgeEndOuterRad * sinPhi, 0));
319 } 298 }
320 else << 299
321 { << 300 // Inner hyperbolic surface
322 G4TwoVector vmin,vmax; << 301
323 G4GeomTools::DiskExtent(rmin, rmax, sphi, << 302 G4double zInnerSplit = 0.;
324 pMin.set(vmin.x(), vmin.y(), zmin); << 303 if (splitInner) {
325 pMax.set(vmax.x(), vmax.y(), zmax); << 304
>> 305 v2 = transform.TransformPoint(
>> 306 G4ThreeVector(maxEndInnerRad * cosPhi, maxEndInnerRad * sinPhi,
>> 307 + fZHalfLength));
>> 308 v3 = transform.TransformPoint(
>> 309 G4ThreeVector(maxEndInnerRad * cosPhi, maxEndInnerRad * sinPhi,
>> 310 - fZHalfLength));
>> 311
>> 312 // Find intersection of tangential line of inner
>> 313 // surface at z = fZHalfLength and line r=fInnerRadius.
>> 314 G4double dr = fZHalfLength * fTanInnerStereo2;
>> 315 G4double dz = maxEndInnerRad;
>> 316 zInnerSplit = fZHalfLength + (fInnerRadius - maxEndInnerRad) * dz / dr;
>> 317
>> 318 // Build associated vertices
>> 319 v5 = transform.TransformPoint(
>> 320 G4ThreeVector(fInnerRadius * cosPhi, fInnerRadius * sinPhi,
>> 321 + zInnerSplit));
>> 322 v6 = transform.TransformPoint(
>> 323 G4ThreeVector(fInnerRadius * cosPhi, fInnerRadius * sinPhi,
>> 324 - zInnerSplit));
>> 325 } else {
>> 326
>> 327 v2 = transform.TransformPoint(
>> 328 G4ThreeVector(fInnerRadius * cosPhi, fInnerRadius * sinPhi,
>> 329 + fZHalfLength));
>> 330 v3 = transform.TransformPoint(
>> 331 G4ThreeVector(fInnerRadius * cosPhi, fInnerRadius * sinPhi,
>> 332 - fZHalfLength));
326 } 333 }
327 334
328 // Check correctness of the bounding box <<
329 // 335 //
330 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 336 // decide vertices of +ve phi boundary
331 { << 337 //
332 std::ostringstream message; << 338
333 message << "Bad bounding box (min >= max) << 339 cosPhi = std::cos(phimax);
334 << GetName() << " !" << 340 sinPhi = std::sin(phimax);
335 << "\npMin = " << pMin << 341
336 << "\npMax = " << pMax; << 342 // Outer hyperbolic surface
337 G4Exception("G4TwistedTubs::BoundingLimits << 343
338 JustWarning, message); << 344 w0 = transform.TransformPoint(
339 DumpInfo(); << 345 G4ThreeVector(fudgeEndOuterRad * cosPhi, fudgeEndOuterRad * sinPhi,
>> 346 + fZHalfLength));
>> 347 w1 = transform.TransformPoint(
>> 348 G4ThreeVector(fudgeEndOuterRad * cosPhi, fudgeEndOuterRad * sinPhi,
>> 349 - fZHalfLength));
>> 350 if (splitOuter) {
>> 351 G4double r = rFudge*fOuterRadius;
>> 352
>> 353 w4 = transform.TransformPoint(G4ThreeVector( r*cosPhi, r*sinPhi, 0 ));
>> 354
>> 355 AddPolyToExtent( v0, v4, w4, w0, voxelLimit, axis, extentList );
>> 356 AddPolyToExtent( v4, v1, w1, w4, voxelLimit, axis, extentList );
>> 357 } else {
>> 358 AddPolyToExtent( v0, v1, w1, w0, voxelLimit, axis, extentList );
>> 359 }
>> 360
>> 361 // Inner hyperbolic surface
>> 362
>> 363 if (splitInner) {
>> 364
>> 365 w2 = transform.TransformPoint(
>> 366 G4ThreeVector(maxEndInnerRad * cosPhi, maxEndInnerRad * sinPhi,
>> 367 + fZHalfLength));
>> 368 w3 = transform.TransformPoint(
>> 369 G4ThreeVector(maxEndInnerRad * cosPhi, maxEndInnerRad * sinPhi,
>> 370 - fZHalfLength));
>> 371
>> 372 w5 = transform.TransformPoint(
>> 373 G4ThreeVector(fInnerRadius * cosPhi, fInnerRadius * sinPhi,
>> 374 + zInnerSplit));
>> 375 w6 = transform.TransformPoint(
>> 376 G4ThreeVector(fInnerRadius * cosPhi, fInnerRadius * sinPhi,
>> 377 - zInnerSplit));
>> 378
>> 379 AddPolyToExtent( v3, v6, w6, w3, voxelLimit, axis, extentList );
>> 380 AddPolyToExtent( v6, v5, w5, w6, voxelLimit, axis, extentList );
>> 381 AddPolyToExtent( v5, v2, w2, w5, voxelLimit, axis, extentList );
>> 382
>> 383 } else {
>> 384 w2 = transform.TransformPoint(
>> 385 G4ThreeVector(fInnerRadius * cosPhi, fInnerRadius * sinPhi,
>> 386 + fZHalfLength));
>> 387 w3 = transform.TransformPoint(
>> 388 G4ThreeVector(fInnerRadius * cosPhi, fInnerRadius * sinPhi,
>> 389 - fZHalfLength));
>> 390
>> 391 AddPolyToExtent( v3, v2, w2, w3, voxelLimit, axis, extentList );
340 } 392 }
>> 393
>> 394 //
>> 395 // Endplate segments
>> 396 //
>> 397 AddPolyToExtent( v1, v3, w3, w1, voxelLimit, axis, extentList );
>> 398 AddPolyToExtent( v2, v0, w0, w2, voxelLimit, axis, extentList );
>> 399
>> 400 //
>> 401 // Return min/max value
>> 402 //
>> 403 return extentList.GetExtent( min, max );
341 } 404 }
342 405
>> 406
343 //============================================ 407 //=====================================================================
344 //* CalculateExtent -------------------------- << 408 //* AddPolyToExtent ---------------------------------------------------
>> 409
>> 410 void G4TwistedTubs::AddPolyToExtent( const G4ThreeVector &v0,
>> 411 const G4ThreeVector &v1,
>> 412 const G4ThreeVector &w1,
>> 413 const G4ThreeVector &w0,
>> 414 const G4VoxelLimits &voxelLimit,
>> 415 const EAxis axis,
>> 416 G4SolidExtentList &extentList )
>> 417 {
>> 418 // Utility function for CalculateExtent
>> 419 //
>> 420 G4ClippablePolygon phiPoly;
>> 421
>> 422 phiPoly.AddVertexInOrder( v0 );
>> 423 phiPoly.AddVertexInOrder( v1 );
>> 424 phiPoly.AddVertexInOrder( w1 );
>> 425 phiPoly.AddVertexInOrder( w0 );
345 426
346 G4bool << 427 if (phiPoly.PartialClip( voxelLimit, axis )) {
347 G4TwistedTubs::CalculateExtent(const EAxis pAx << 428 phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
348 const G4VoxelLi << 429 extentList.AddSurface( phiPoly );
349 const G4AffineT << 430 }
350 G4double& <<
351 { <<
352 G4ThreeVector bmin, bmax; <<
353 <<
354 // Get bounding box <<
355 BoundingLimits(bmin,bmax); <<
356 <<
357 // Find extent <<
358 G4BoundingEnvelope bbox(bmin,bmax); <<
359 return bbox.CalculateExtent(pAxis,pVoxelLimi <<
360 } 431 }
361 432
362 433
363 //============================================ 434 //=====================================================================
364 //* Inside ----------------------------------- 435 //* Inside ------------------------------------------------------------
365 436
366 EInside G4TwistedTubs::Inside(const G4ThreeVec 437 EInside G4TwistedTubs::Inside(const G4ThreeVector& p) const
367 { 438 {
368 439
369 const G4double halftol << 440 static const G4double halftol = 0.5 * kRadTolerance;
370 = 0.5 * G4GeometryTolerance::GetInstance( <<
371 // static G4int timerid = -1; 441 // static G4int timerid = -1;
372 // G4Timer timer(timerid, "G4TwistedTubs", 442 // G4Timer timer(timerid, "G4TwistedTubs", "Inside");
373 // timer.Start(); 443 // timer.Start();
>> 444
374 445
>> 446
>> 447 G4ThreeVector *tmpp;
>> 448 EInside *tmpinside;
>> 449 if (fLastInside.p == p) {
>> 450
>> 451 return fLastInside.inside;
>> 452
>> 453 } else {
>> 454 tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
>> 455 tmpinside = const_cast<EInside*>(&(fLastInside.inside));
>> 456 tmpp->set(p.x(), p.y(), p.z());
>> 457 }
375 458
376 EInside outerhypearea = ((G4TwistTubsHypeS << 459
377 G4double innerhyperho = ((G4TwistTubsHypeS << 460 EInside outerhypearea = ((G4HyperbolicSurface *)fOuterHype)->Inside(p);
>> 461 G4double innerhyperho = ((G4HyperbolicSurface *)fInnerHype)->GetRhoAtPZ(p);
378 G4double distanceToOut = p.getRho() - inner 462 G4double distanceToOut = p.getRho() - innerhyperho; // +ve: inside
379 EInside tmpinside; << 463
380 if ((outerhypearea == kOutside) || (distanc << 464 if (outerhypearea == kOutside || distanceToOut < -halftol) {
381 { << 465 *tmpinside = kOutside;
382 tmpinside = kOutside; << 466 } else if (outerhypearea == kSurface) {
383 } << 467 *tmpinside = kSurface;
384 else if (outerhypearea == kSurface) << 468 } else {
385 { << 469 if (distanceToOut <= halftol) {
386 tmpinside = kSurface; << 470 *tmpinside = kSurface;
387 } << 471 } else {
388 else << 472 *tmpinside = kInside;
389 { <<
390 if (distanceToOut <= halftol) <<
391 { <<
392 tmpinside = kSurface; <<
393 } <<
394 else <<
395 { <<
396 tmpinside = kInside; <<
397 } 473 }
398 } 474 }
>> 475
>> 476 return fLastInside.inside;
399 477
400 return tmpinside; <<
401 } 478 }
402 479
403 //============================================ 480 //=====================================================================
404 //* SurfaceNormal ---------------------------- 481 //* SurfaceNormal -----------------------------------------------------
405 482
406 G4ThreeVector G4TwistedTubs::SurfaceNormal(con 483 G4ThreeVector G4TwistedTubs::SurfaceNormal(const G4ThreeVector& p) const
407 { 484 {
408 // 485 //
409 // return the normal unit vector to the Hyp 486 // return the normal unit vector to the Hyperbolical Surface at a point
410 // p on (or nearly on) the surface 487 // p on (or nearly on) the surface
411 // 488 //
412 // Which of the three or four surfaces are 489 // Which of the three or four surfaces are we closest to?
413 // 490 //
414 491
>> 492 // static G4int timerid = -1;
>> 493 // G4Timer timer(timerid, "G4TwistedTubs", "SurfaceNormal");
>> 494 // timer.Start();
>> 495
>> 496
>> 497
>> 498
>> 499 if (fLastNormal.p == p) {
>> 500
>> 501 return fLastNormal.vec;
>> 502 }
>> 503 G4ThreeVector *tmpp = const_cast<G4ThreeVector*>(&(fLastNormal.p));
>> 504 G4ThreeVector *tmpnormal = const_cast<G4ThreeVector*>(&(fLastNormal.vec));
>> 505 G4VSurface **tmpsurface = const_cast<G4VSurface**>(fLastNormal.surface);
>> 506 tmpp->set(p.x(), p.y(), p.z());
415 507
416 G4double distance = kInfinity; 508 G4double distance = kInfinity;
417 509
418 G4VTwistSurface *surfaces[6]; << 510 G4VSurface *surfaces[6];
419 surfaces[0] = fLatterTwisted; 511 surfaces[0] = fLatterTwisted;
420 surfaces[1] = fFormerTwisted; 512 surfaces[1] = fFormerTwisted;
421 surfaces[2] = fInnerHype; 513 surfaces[2] = fInnerHype;
422 surfaces[3] = fOuterHype; 514 surfaces[3] = fOuterHype;
423 surfaces[4] = fLowerEndcap; 515 surfaces[4] = fLowerEndcap;
424 surfaces[5] = fUpperEndcap; 516 surfaces[5] = fUpperEndcap;
425 517
426 G4ThreeVector xx; 518 G4ThreeVector xx;
427 G4ThreeVector bestxx; 519 G4ThreeVector bestxx;
>> 520 G4int i;
428 G4int besti = -1; 521 G4int besti = -1;
429 for (auto i=0; i<6; ++i) << 522 for (i=0; i< 6; i++) {
430 { <<
431 G4double tmpdistance = surfaces[i]->Dist 523 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
432 if (tmpdistance < distance) << 524 if (tmpdistance < distance) {
433 { <<
434 distance = tmpdistance; 525 distance = tmpdistance;
435 bestxx = xx; 526 bestxx = xx;
436 besti = i; 527 besti = i;
437 } 528 }
438 } 529 }
439 530
440 return surfaces[besti]->GetNormal(bestxx, tr << 531 tmpsurface[0] = surfaces[besti];
>> 532 *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
>> 533
>> 534 return fLastNormal.vec;
441 } 535 }
442 536
443 //============================================ 537 //=====================================================================
444 //* DistanceToIn (p, v) ---------------------- 538 //* DistanceToIn (p, v) -----------------------------------------------
445 539
446 G4double G4TwistedTubs::DistanceToIn (const G4 540 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p,
447 const G4 541 const G4ThreeVector& v ) const
448 { 542 {
449 543
450 // DistanceToIn (p, v): 544 // DistanceToIn (p, v):
451 // Calculate distance to surface of shape f 545 // Calculate distance to surface of shape from `outside'
452 // along with the v, allowing for tolerance 546 // along with the v, allowing for tolerance.
453 // The function returns kInfinity if no int 547 // The function returns kInfinity if no intersection or
454 // just grazing within tolerance. 548 // just grazing within tolerance.
455 549
456 // 550 //
>> 551 // checking last value
>> 552 //
>> 553
>> 554 G4ThreeVector *tmpp;
>> 555 G4ThreeVector *tmpv;
>> 556 G4double *tmpdist;
>> 557 if (fLastDistanceToInWithV.p == p && fLastDistanceToInWithV.vec == v) {
>> 558
>> 559 return fLastDistanceToIn.value;
>> 560
>> 561 } else {
>> 562 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
>> 563 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
>> 564 tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
>> 565 tmpp->set(p.x(), p.y(), p.z());
>> 566 tmpv->set(v.x(), v.y(), v.z());
>> 567 }
>> 568
>> 569 //
457 // Calculate DistanceToIn(p,v) 570 // Calculate DistanceToIn(p,v)
458 // 571 //
459 572
460 EInside currentside = Inside(p); 573 EInside currentside = Inside(p);
461 574
462 if (currentside == kInside) << 575 if (currentside == kInside) {
463 { << 576
464 } << 577 } else if (currentside == kSurface) {
465 else << 578 // particle is just on a boundary.
466 { << 579 // if the particle is entering to the volume, return 0.
467 if (currentside == kSurface) << 580 G4ThreeVector normal = SurfaceNormal(p);
468 { << 581 if (normal*v < 0) {
469 // particle is just on a boundary. << 582
470 // If the particle is entering to the v << 583 *tmpdist = 0;
471 // << 584 return fLastDistanceToInWithV.value;
472 G4ThreeVector normal = SurfaceNormal(p) << 585 }
473 if (normal*v < 0) <<
474 { <<
475 return 0; <<
476 } <<
477 } <<
478 } 586 }
479 587
480 // now, we can take smallest positive dista 588 // now, we can take smallest positive distance.
481 589
482 // Initialize 590 // Initialize
483 // <<
484 G4double distance = kInfinity; 591 G4double distance = kInfinity;
485 592
486 // find intersections and choose nearest on 593 // find intersections and choose nearest one.
487 // << 594 G4VSurface *surfaces[6];
488 G4VTwistSurface* surfaces[6]; <<
489 surfaces[0] = fLowerEndcap; 595 surfaces[0] = fLowerEndcap;
490 surfaces[1] = fUpperEndcap; 596 surfaces[1] = fUpperEndcap;
491 surfaces[2] = fLatterTwisted; 597 surfaces[2] = fLatterTwisted;
492 surfaces[3] = fFormerTwisted; 598 surfaces[3] = fFormerTwisted;
493 surfaces[4] = fInnerHype; 599 surfaces[4] = fInnerHype;
494 surfaces[5] = fOuterHype; 600 surfaces[5] = fOuterHype;
495 601
496 G4ThreeVector xx; 602 G4ThreeVector xx;
497 G4ThreeVector bestxx; 603 G4ThreeVector bestxx;
498 for (const auto & surface : surfaces) << 604 G4int i;
499 { << 605 G4int besti = -1;
500 G4double tmpdistance = surface->Distance << 606 for (i=0; i< 6; i++) {
501 if (tmpdistance < distance) << 607 G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
502 { << 608 if (tmpdistance < distance) {
503 distance = tmpdistance; 609 distance = tmpdistance;
504 bestxx = xx; 610 bestxx = xx;
>> 611 besti = i;
505 } 612 }
506 } 613 }
507 return distance; << 614
>> 615 *tmpdist = distance;
>> 616 // timer.Stop();
>> 617 return fLastDistanceToInWithV.value;
508 } 618 }
509 619
510 //============================================ 620 //=====================================================================
511 //* DistanceToIn (p) ------------------------- 621 //* DistanceToIn (p) --------------------------------------------------
512 622
513 G4double G4TwistedTubs::DistanceToIn (const G4 623 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p) const
514 { 624 {
515 // DistanceToIn(p): 625 // DistanceToIn(p):
516 // Calculate distance to surface of shape f 626 // Calculate distance to surface of shape from `outside',
517 // allowing for tolerance 627 // allowing for tolerance
>> 628 //
>> 629
>> 630 //
>> 631 // checking last value
>> 632 //
>> 633
>> 634 G4ThreeVector *tmpp;
>> 635 G4double *tmpdist;
>> 636 if (fLastDistanceToIn.p == p) {
>> 637
>> 638 return fLastDistanceToIn.value;
>> 639 } else {
>> 640 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
>> 641 tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
>> 642 tmpp->set(p.x(), p.y(), p.z());
>> 643 }
518 644
519 // 645 //
520 // Calculate DistanceToIn(p) 646 // Calculate DistanceToIn(p)
521 // 647 //
522 648
523 EInside currentside = Inside(p); 649 EInside currentside = Inside(p);
524 650
525 switch (currentside) << 651 switch (currentside) {
526 { << 652
527 case (kInside) : << 653 case (kInside) : {
528 {} << 654
529 case (kSurface) : << 655 }
530 { << 656
531 return 0; << 657 case (kSurface) : {
>> 658 *tmpdist = 0.;
>> 659 return fLastDistanceToIn.value;
532 } 660 }
533 case (kOutside) : << 661
534 { << 662 case (kOutside) : {
535 // Initialize 663 // Initialize
536 G4double distance = kInfinity; << 664 G4double distance = kInfinity;
537 665
538 // find intersections and choose near 666 // find intersections and choose nearest one.
539 G4VTwistSurface *surfaces[6]; << 667 G4VSurface *surfaces[6];
540 surfaces[0] = fLowerEndcap; 668 surfaces[0] = fLowerEndcap;
541 surfaces[1] = fUpperEndcap; 669 surfaces[1] = fUpperEndcap;
542 surfaces[2] = fLatterTwisted; 670 surfaces[2] = fLatterTwisted;
543 surfaces[3] = fFormerTwisted; 671 surfaces[3] = fFormerTwisted;
544 surfaces[4] = fInnerHype; 672 surfaces[4] = fInnerHype;
545 surfaces[5] = fOuterHype; 673 surfaces[5] = fOuterHype;
546 674
>> 675 G4int i;
>> 676 G4int besti = -1;
547 G4ThreeVector xx; 677 G4ThreeVector xx;
548 G4ThreeVector bestxx; 678 G4ThreeVector bestxx;
549 for (const auto & surface : surfaces) << 679 for (i=0; i< 6; i++) {
550 { << 680 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
551 G4double tmpdistance = surface->Di << 681 if (tmpdistance < distance) {
552 if (tmpdistance < distance) <<
553 { <<
554 distance = tmpdistance; 682 distance = tmpdistance;
555 bestxx = xx; 683 bestxx = xx;
>> 684 besti = i;
556 } 685 }
557 } 686 }
558 return distance; << 687
>> 688
>> 689 *tmpdist = distance;
>> 690 return fLastDistanceToIn.value;
559 } 691 }
560 default : << 692
561 { << 693 default : {
562 G4Exception("G4TwistedTubs::DistanceT << 694 G4Exception("G4TwistedTubs::DistanceToIn(p)", "InvalidCondition",
563 FatalException, "Unknown 695 FatalException, "Unknown point location!");
564 } 696 }
565 } // switch end 697 } // switch end
566 <<
567 return kInfinity; 698 return kInfinity;
568 } 699 }
569 700
570 //============================================ 701 //=====================================================================
571 //* DistanceToOut (p, v) --------------------- 702 //* DistanceToOut (p, v) ----------------------------------------------
572 703
573 G4double G4TwistedTubs::DistanceToOut( const G 704 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p,
574 const G 705 const G4ThreeVector& v,
575 const G 706 const G4bool calcNorm,
576 G4bool << 707 G4bool *validNorm, G4ThreeVector *norm ) const
577 G4Three <<
578 { 708 {
579 // DistanceToOut (p, v): 709 // DistanceToOut (p, v):
580 // Calculate distance to surface of shape f 710 // Calculate distance to surface of shape from `inside'
581 // along with the v, allowing for tolerance 711 // along with the v, allowing for tolerance.
582 // The function returns kInfinity if no int 712 // The function returns kInfinity if no intersection or
583 // just grazing within tolerance. 713 // just grazing within tolerance.
>> 714 //
>> 715
>> 716
>> 717 //
>> 718 // checking last value
>> 719 //
>> 720
>> 721 G4ThreeVector *tmpp;
>> 722 G4ThreeVector *tmpv;
>> 723 G4double *tmpdist;
>> 724 if (fLastDistanceToOutWithV.p == p && fLastDistanceToOutWithV.vec == v ) {
>> 725
>> 726 return fLastDistanceToOutWithV.value;
>> 727 } else {
>> 728 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
>> 729 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
>> 730 tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
>> 731 tmpp->set(p.x(), p.y(), p.z());
>> 732 tmpv->set(v.x(), v.y(), v.z());
>> 733 }
584 734
585 // 735 //
586 // Calculate DistanceToOut(p,v) 736 // Calculate DistanceToOut(p,v)
587 // 737 //
588 738
589 EInside currentside = Inside(p); 739 EInside currentside = Inside(p);
590 if (currentside == kOutside) << 740
591 { << 741 if (currentside == kOutside) {
592 } << 742
593 else << 743 } else if (currentside == kSurface) {
594 { << 744 // particle is just on a boundary.
595 if (currentside == kSurface) << 745 // if the particle is exiting from the volume, return 0.
596 { << 746 G4ThreeVector normal = SurfaceNormal(p);
597 // particle is just on a boundary. << 747 G4VSurface *blockedsurface = fLastNormal.surface[0];
598 // If the particle is exiting from the << 748 if (normal*v > 0) {
599 // << 749
600 G4ThreeVector normal = SurfaceNormal(p) << 750 if (calcNorm) {
601 if (normal*v > 0) << 751 *norm = (blockedsurface->GetNormal(p, true));
602 { << 752 *validNorm = blockedsurface->IsValidNorm();
603 if (calcNorm) << 753 }
604 { << 754 *tmpdist = 0.;
605 *norm = normal; << 755 return fLastDistanceToOutWithV.value;
606 *validNorm = true; << 756 }
607 } <<
608 return 0; <<
609 } <<
610 } <<
611 } 757 }
612 758
613 // now, we can take smallest positive dista 759 // now, we can take smallest positive distance.
614 760
615 // Initialize 761 // Initialize
616 // << 762 G4double distance = kInfinity;
617 G4double distance = kInfinity; <<
618 763
619 // find intersections and choose nearest on 764 // find intersections and choose nearest one.
620 // << 765 G4VSurface *surfaces[6];
621 G4VTwistSurface* surfaces[6]; <<
622 surfaces[0] = fLatterTwisted; 766 surfaces[0] = fLatterTwisted;
623 surfaces[1] = fFormerTwisted; 767 surfaces[1] = fFormerTwisted;
624 surfaces[2] = fInnerHype; 768 surfaces[2] = fInnerHype;
625 surfaces[3] = fOuterHype; 769 surfaces[3] = fOuterHype;
626 surfaces[4] = fLowerEndcap; 770 surfaces[4] = fLowerEndcap;
627 surfaces[5] = fUpperEndcap; 771 surfaces[5] = fUpperEndcap;
628 772
>> 773 G4int i;
629 G4int besti = -1; 774 G4int besti = -1;
630 G4ThreeVector xx; 775 G4ThreeVector xx;
631 G4ThreeVector bestxx; 776 G4ThreeVector bestxx;
632 for (auto i=0; i<6; ++i) << 777 for (i=0; i< 6; i++) {
633 { <<
634 G4double tmpdistance = surfaces[i]->Dist 778 G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
635 if (tmpdistance < distance) << 779 if (tmpdistance < distance) {
636 { <<
637 distance = tmpdistance; 780 distance = tmpdistance;
638 bestxx = xx; 781 bestxx = xx;
639 besti = i; 782 besti = i;
640 } 783 }
641 } 784 }
642 785
643 if (calcNorm) << 786 if (calcNorm) {
644 { << 787 if (besti != -1) {
645 if (besti != -1) <<
646 { <<
647 *norm = (surfaces[besti]->GetNormal(p 788 *norm = (surfaces[besti]->GetNormal(p, true));
648 *validNorm = surfaces[besti]->IsValid 789 *validNorm = surfaces[besti]->IsValidNorm();
649 } 790 }
650 } 791 }
651 792
652 return distance; << 793 *tmpdist = distance;
>> 794 // timer.Stop();
>> 795 return fLastDistanceToOutWithV.value;
653 } 796 }
654 797
655 798
656 //============================================ 799 //=====================================================================
657 //* DistanceToOut (p) ------------------------ 800 //* DistanceToOut (p) ----------------------------------------------
658 801
659 G4double G4TwistedTubs::DistanceToOut( const G 802 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p ) const
660 { 803 {
661 // DistanceToOut(p): 804 // DistanceToOut(p):
662 // Calculate distance to surface of shape f 805 // Calculate distance to surface of shape from `inside',
663 // allowing for tolerance 806 // allowing for tolerance
>> 807 //
>> 808
664 809
665 // 810 //
>> 811 // checking last value
>> 812 //
>> 813
>> 814 G4ThreeVector *tmpp;
>> 815 G4double *tmpdist;
>> 816 if (fLastDistanceToOut.p == p) {
>> 817 return fLastDistanceToOut.value;
>> 818 } else {
>> 819 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
>> 820 tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
>> 821 tmpp->set(p.x(), p.y(), p.z());
>> 822 }
>> 823
>> 824 //
666 // Calculate DistanceToOut(p) 825 // Calculate DistanceToOut(p)
667 // 826 //
668 827
669 EInside currentside = Inside(p); 828 EInside currentside = Inside(p);
670 829
671 switch (currentside) << 830 switch (currentside) {
672 { << 831 case (kOutside) : {
673 case (kOutside) : << 832
674 { <<
675 } 833 }
676 case (kSurface) : << 834 case (kSurface) : {
677 { << 835 *tmpdist = 0.;
678 return 0; << 836 return fLastDistanceToOut.value;
679 } 837 }
680 case (kInside) : << 838
681 { << 839 case (kInside) : {
682 // Initialize 840 // Initialize
683 G4double distance = kInfinity; 841 G4double distance = kInfinity;
684 842
685 // find intersections and choose near 843 // find intersections and choose nearest one.
686 G4VTwistSurface* surfaces[6]; << 844 G4VSurface *surfaces[6];
687 surfaces[0] = fLatterTwisted; 845 surfaces[0] = fLatterTwisted;
688 surfaces[1] = fFormerTwisted; 846 surfaces[1] = fFormerTwisted;
689 surfaces[2] = fInnerHype; 847 surfaces[2] = fInnerHype;
690 surfaces[3] = fOuterHype; 848 surfaces[3] = fOuterHype;
691 surfaces[4] = fLowerEndcap; 849 surfaces[4] = fLowerEndcap;
692 surfaces[5] = fUpperEndcap; 850 surfaces[5] = fUpperEndcap;
693 851
>> 852 G4int i;
>> 853 G4int besti = -1;
694 G4ThreeVector xx; 854 G4ThreeVector xx;
695 G4ThreeVector bestxx; 855 G4ThreeVector bestxx;
696 for (const auto & surface : surfaces) << 856 for (i=0; i< 6; i++) {
697 { << 857 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
698 G4double tmpdistance = surface->Di << 858 if (tmpdistance < distance) {
699 if (tmpdistance < distance) <<
700 { <<
701 distance = tmpdistance; 859 distance = tmpdistance;
702 bestxx = xx; 860 bestxx = xx;
>> 861 besti = i;
703 } 862 }
704 } 863 }
705 return distance; << 864
>> 865
>> 866 *tmpdist = distance;
>> 867
>> 868 return fLastDistanceToOut.value;
706 } 869 }
707 default : << 870
708 { << 871 default : {
709 G4Exception("G4TwistedTubs::DistanceT << 872 G4Exception("G4TwistedTubs::DistanceToOut(p)", "InvalidCondition",
710 FatalException, "Unknown 873 FatalException, "Unknown point location!");
711 } 874 }
712 } // switch end 875 } // switch end
713 << 876 return 0;
714 return 0.; <<
715 } 877 }
716 878
717 //============================================ 879 //=====================================================================
718 //* StreamInfo ------------------------------- 880 //* StreamInfo --------------------------------------------------------
719 881
720 std::ostream& G4TwistedTubs::StreamInfo(std::o 882 std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
721 { 883 {
722 // 884 //
723 // Stream object contents to an output strea 885 // Stream object contents to an output stream
724 // 886 //
725 G4long oldprc = os.precision(16); <<
726 os << "------------------------------------- 887 os << "-----------------------------------------------------------\n"
727 << " *** Dump for solid - " << GetName 888 << " *** Dump for solid - " << GetName() << " ***\n"
728 << " ================================= 889 << " ===================================================\n"
729 << " Solid type: G4TwistedTubs\n" 890 << " Solid type: G4TwistedTubs\n"
730 << " Parameters: \n" 891 << " Parameters: \n"
731 << " -ve end Z : " << fEn 892 << " -ve end Z : " << fEndZ[0]/mm << " mm \n"
732 << " +ve end Z : " << fEn 893 << " +ve end Z : " << fEndZ[1]/mm << " mm \n"
733 << " inner end radius(-ve z): " << fEn 894 << " inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
734 << " inner end radius(+ve z): " << fEn 895 << " inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
735 << " outer end radius(-ve z): " << fEn 896 << " outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
736 << " outer end radius(+ve z): " << fEn 897 << " outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
737 << " inner radius (z=0) : " << fIn 898 << " inner radius (z=0) : " << fInnerRadius/mm << " mm \n"
738 << " outer radius (z=0) : " << fOu 899 << " outer radius (z=0) : " << fOuterRadius/mm << " mm \n"
739 << " twisted angle : " << fPh 900 << " twisted angle : " << fPhiTwist/degree << " degrees \n"
740 << " inner stereo angle : " << fIn 901 << " inner stereo angle : " << fInnerStereo/degree << " degrees \n"
741 << " outer stereo angle : " << fOu 902 << " outer stereo angle : " << fOuterStereo/degree << " degrees \n"
742 << " phi-width of a piece : " << fDP 903 << " phi-width of a piece : " << fDPhi/degree << " degrees \n"
743 << "------------------------------------- 904 << "-----------------------------------------------------------\n";
744 os.precision(oldprc); <<
745 905
746 return os; 906 return os;
747 } 907 }
748 908
749 909
750 //============================================ 910 //=====================================================================
751 //* DiscribeYourselfTo ----------------------- 911 //* DiscribeYourselfTo ------------------------------------------------
752 912
753 void G4TwistedTubs::DescribeYourselfTo (G4VGra 913 void G4TwistedTubs::DescribeYourselfTo (G4VGraphicsScene& scene) const
754 { 914 {
755 scene.AddSolid (*this); 915 scene.AddSolid (*this);
756 } 916 }
757 917
758 //============================================ 918 //=====================================================================
759 //* GetExtent -------------------------------- 919 //* GetExtent ---------------------------------------------------------
760 920
761 G4VisExtent G4TwistedTubs::GetExtent() const << 921 G4VisExtent G4TwistedTubs::GetExtent() const
762 { 922 {
763 // Define the sides of the box into which th 923 // Define the sides of the box into which the G4Tubs instance would fit.
764 // << 924 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
765 G4ThreeVector pmin,pmax; << 925 return G4VisExtent( -maxEndOuterRad, maxEndOuterRad,
766 BoundingLimits(pmin,pmax); << 926 -maxEndOuterRad, maxEndOuterRad,
767 return { pmin.x(),pmax.x(), << 927 -fZHalfLength, fZHalfLength );
768 pmin.y(),pmax.y(), <<
769 pmin.z(),pmax.z() }; <<
770 } 928 }
771 929
772 //============================================ 930 //=====================================================================
773 //* CreatePolyhedron ------------------------- 931 //* CreatePolyhedron --------------------------------------------------
774 932
775 G4Polyhedron* G4TwistedTubs::CreatePolyhedron 933 G4Polyhedron* G4TwistedTubs::CreatePolyhedron () const
776 { 934 {
777 // number of meshes << 935 // Tube for now!!!
778 // 936 //
779 G4double absPhiTwist = std::abs(fPhiTwist); << 937 return new G4PolyhedronTubs (fInnerRadius, fOuterRadius, fZHalfLength,
780 G4double dA = std::max(fDPhi,absPhiTwist); << 938 fEndPhi[0], fDPhi);
781 const G4int k = << 939 }
782 G4int(G4Polyhedron::GetNumberOfRotationSte <<
783 const G4int n = <<
784 G4int(G4Polyhedron::GetNumberOfRotationSte <<
785 <<
786 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ; <<
787 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1) <<
788 <<
789 auto ph = new G4Polyhedron; <<
790 typedef G4double G4double3[3]; <<
791 typedef G4int G4int4[4]; <<
792 auto xyz = new G4double3[nnodes]; // number <<
793 auto faces = new G4int4[nfaces] ; // number <<
794 fLowerEndcap->GetFacets(k,k,xyz,faces,0) ; <<
795 fUpperEndcap->GetFacets(k,k,xyz,faces,1) ; <<
796 fInnerHype->GetFacets(k,n,xyz,faces,2) ; <<
797 fFormerTwisted->GetFacets(k,n,xyz,faces,3) ; <<
798 fOuterHype->GetFacets(k,n,xyz,faces,4) ; <<
799 fLatterTwisted->GetFacets(k,n,xyz,faces,5) ; <<
800 <<
801 ph->createPolyhedron(nnodes,nfaces,xyz,faces <<
802 940
803 delete[] xyz; << 941 //=====================================================================
804 delete[] faces; << 942 //* CreateNUBS --------------------------------------------------------
805 943
806 return ph; << 944 G4NURBS* G4TwistedTubs::CreateNURBS () const
>> 945 {
>> 946 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
>> 947 G4double maxEndInnerRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
>> 948 return new G4NURBStube(maxEndInnerRad, maxEndOuterRad, fZHalfLength);
>> 949 // Tube for now!!!
807 } 950 }
808 951
809 //============================================ 952 //=====================================================================
810 //* GetPolyhedron ---------------------------- 953 //* GetPolyhedron -----------------------------------------------------
811 954
812 G4Polyhedron* G4TwistedTubs::GetPolyhedron () 955 G4Polyhedron* G4TwistedTubs::GetPolyhedron () const
813 { 956 {
814 if (fpPolyhedron == nullptr || << 957 if (!fpPolyhedron ||
815 fRebuildPolyhedron || <<
816 fpPolyhedron->GetNumberOfRotationStepsAt 958 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
817 fpPolyhedron->GetNumberOfRotationSteps() 959 fpPolyhedron->GetNumberOfRotationSteps())
818 { << 960 {
819 G4AutoLock l(&polyhedronMutex); << 961 delete fpPolyhedron;
820 delete fpPolyhedron; << 962 fpPolyhedron = CreatePolyhedron();
821 fpPolyhedron = CreatePolyhedron(); << 963 }
822 fRebuildPolyhedron = false; <<
823 l.unlock(); <<
824 } <<
825 return fpPolyhedron; 964 return fpPolyhedron;
826 } 965 }
827 966
828 //============================================ 967 //=====================================================================
829 //* CreateSurfaces --------------------------- 968 //* CreateSurfaces ----------------------------------------------------
830 969
831 void G4TwistedTubs::CreateSurfaces() 970 void G4TwistedTubs::CreateSurfaces()
832 { 971 {
833 // create 6 surfaces of TwistedTub << 972
834 << 973 // create 6 surfaces of TwistedTub.
835 G4ThreeVector x0(0, 0, fEndZ[0]); 974 G4ThreeVector x0(0, 0, fEndZ[0]);
836 G4ThreeVector n (0, 0, -1); 975 G4ThreeVector n (0, 0, -1);
837 976
838 fLowerEndcap = new G4TwistTubsFlatSide("Low << 977 fLowerEndcap = new G4FlatSurface("LowerEndcap",
839 fEndInnerR 978 fEndInnerRadius, fEndOuterRadius,
840 fDPhi, fEn 979 fDPhi, fEndPhi, fEndZ, -1) ;
841 980
842 fUpperEndcap = new G4TwistTubsFlatSide("Upp << 981 fUpperEndcap = new G4FlatSurface("UpperEndcap",
843 fEndInnerR 982 fEndInnerRadius, fEndOuterRadius,
844 fDPhi, fEn 983 fDPhi, fEndPhi, fEndZ, 1) ;
845 984
846 G4RotationMatrix rotHalfDPhi; 985 G4RotationMatrix rotHalfDPhi;
847 rotHalfDPhi.rotateZ(0.5*fDPhi); 986 rotHalfDPhi.rotateZ(0.5*fDPhi);
848 987
849 fLatterTwisted = new G4TwistTubsSide("Latte << 988 fLatterTwisted = new G4TwistedSurface("LatterTwisted",
850 fEndI 989 fEndInnerRadius, fEndOuterRadius,
851 fDPhi 990 fDPhi, fEndPhi, fEndZ,
852 fInne 991 fInnerRadius, fOuterRadius, fKappa,
853 1 ) ; 992 1 ) ;
854 fFormerTwisted = new G4TwistTubsSide("Forme << 993 fFormerTwisted = new G4TwistedSurface("FormerTwisted",
855 fEndI 994 fEndInnerRadius, fEndOuterRadius,
856 fDPhi 995 fDPhi, fEndPhi, fEndZ,
857 fInne 996 fInnerRadius, fOuterRadius, fKappa,
858 -1 ) 997 -1 ) ;
859 998
860 fInnerHype = new G4TwistTubsHypeSide("Inner << 999 fInnerHype = new G4HyperbolicSurface("InnerHype",
861 fEndIn 1000 fEndInnerRadius, fEndOuterRadius,
862 fDPhi, 1001 fDPhi, fEndPhi, fEndZ,
863 fInner 1002 fInnerRadius, fOuterRadius,fKappa,
864 fTanIn 1003 fTanInnerStereo, fTanOuterStereo, -1) ;
865 fOuterHype = new G4TwistTubsHypeSide("Outer << 1004 fOuterHype = new G4HyperbolicSurface("OuterHype",
866 fEndIn 1005 fEndInnerRadius, fEndOuterRadius,
867 fDPhi, 1006 fDPhi, fEndPhi, fEndZ,
868 fInner 1007 fInnerRadius, fOuterRadius,fKappa,
869 fTanIn 1008 fTanInnerStereo, fTanOuterStereo, 1) ;
870 1009
871 1010
872 // set neighbour surfaces << 1011 // set neighbour surfaces.
873 // <<
874 fLowerEndcap->SetNeighbours(fInnerHype, fLa 1012 fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
875 fOuterHype, fFo 1013 fOuterHype, fFormerTwisted);
876 fUpperEndcap->SetNeighbours(fInnerHype, fLa 1014 fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
877 fOuterHype, fFo 1015 fOuterHype, fFormerTwisted);
878 fLatterTwisted->SetNeighbours(fInnerHype, f 1016 fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
879 fOuterHype, f 1017 fOuterHype, fUpperEndcap);
880 fFormerTwisted->SetNeighbours(fInnerHype, f 1018 fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
881 fOuterHype, f 1019 fOuterHype, fUpperEndcap);
882 fInnerHype->SetNeighbours(fLatterTwisted, f 1020 fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
883 fFormerTwisted, f 1021 fFormerTwisted, fUpperEndcap);
884 fOuterHype->SetNeighbours(fLatterTwisted, f 1022 fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
885 fFormerTwisted, f 1023 fFormerTwisted, fUpperEndcap);
886 } 1024 }
887 1025
888 1026
889 //============================================ 1027 //=====================================================================
890 //* GetEntityType ---------------------------- 1028 //* GetEntityType -----------------------------------------------------
891 1029
892 G4GeometryType G4TwistedTubs::GetEntityType() 1030 G4GeometryType G4TwistedTubs::GetEntityType() const
893 { 1031 {
894 return {"G4TwistedTubs"}; << 1032 return G4String("G4TwistedTubs");
895 } <<
896 <<
897 //============================================ <<
898 //* Clone ------------------------------------ <<
899 <<
900 G4VSolid* G4TwistedTubs::Clone() const <<
901 { <<
902 return new G4TwistedTubs(*this); <<
903 } 1033 }
904 1034
905 //============================================ 1035 //=====================================================================
906 //* GetCubicVolume --------------------------- 1036 //* GetCubicVolume ----------------------------------------------------
907 1037
908 G4double G4TwistedTubs::GetCubicVolume() 1038 G4double G4TwistedTubs::GetCubicVolume()
909 { 1039 {
910 if (fCubicVolume == 0.) << 1040 if(fCubicVolume != 0.) ;
911 { << 1041 else fCubicVolume = G4VSolid::GetCubicVolume();
912 G4double DPhi = GetDPhi(); <<
913 G4double Z0 = GetEndZ(0); <<
914 G4double Z1 = GetEndZ(1); <<
915 G4double Ain = GetInnerRadius(); <<
916 G4double Aout = GetOuterRadius(); <<
917 G4double R0in = GetEndInnerRadius(0); <<
918 G4double R1in = GetEndInnerRadius(1); <<
919 G4double R0out = GetEndOuterRadius(0); <<
920 G4double R1out = GetEndOuterRadius(1); <<
921 <<
922 // V_hyperboloid = pi*h*(2*a*a + R*R)/3 <<
923 fCubicVolume = (2.*(Z1 - Z0)*(Aout + Ain)* <<
924 + Z1*(R1out + R1in)*(R1out <<
925 - Z0*(R0out + R0in)*(R0out <<
926 } <<
927 return fCubicVolume; 1042 return fCubicVolume;
928 } <<
929 <<
930 //============================================ <<
931 //* GetLateralArea --------------------------- <<
932 <<
933 G4double <<
934 G4TwistedTubs::GetLateralArea(G4double a, G4do <<
935 { <<
936 if (z == 0) return 0.; <<
937 G4double h = std::abs(z); <<
938 G4double area = h*a; <<
939 if (std::abs(a - r) > kCarTolerance) <<
940 { <<
941 G4double aa = a*a; <<
942 G4double hh = h*h; <<
943 G4double rr = r*r; <<
944 G4double cc = aa*hh/(rr - aa); <<
945 G4double k = std::sqrt(aa + cc)/cc; <<
946 G4double kh = k*h; <<
947 area = 0.5*a*(h*std::sqrt(1. + kh*kh) + st <<
948 } <<
949 return GetDPhi()*area; <<
950 } <<
951 <<
952 //============================================ <<
953 //* GetPhiCutArea ---------------------------- <<
954 <<
955 G4double <<
956 G4TwistedTubs::GetPhiCutArea(G4double a, G4dou <<
957 { <<
958 if (GetDPhi() >= CLHEP::twopi || r <= 0 || z <<
959 G4double h = std::abs(z); <<
960 G4double area = h*a; <<
961 if (GetPhiTwist() > kCarTolerance) <<
962 { <<
963 G4double sinw = std::sin(0.5*GetPhiTwist() <<
964 G4double p = sinw*r/h; <<
965 G4double q = sinw*r/a; <<
966 G4double pp = p*p; <<
967 G4double qq = q*q; <<
968 G4double pq = p*q; <<
969 G4double sqroot = std::sqrt(pp + qq + 1); <<
970 area = (pq*sqroot + <<
971 0.5*p*(pp + 3.)*std::atanh(q/sqroo <<
972 0.5*q*(qq + 3.)*std::atanh(p/sqroo <<
973 std::atan(sqroot/(pq)) - CLHEP::ha <<
974 } <<
975 return area; <<
976 } <<
977 <<
978 //============================================ <<
979 //* GetSurfaceArea --------------------------- <<
980 <<
981 G4double G4TwistedTubs::GetSurfaceArea() <<
982 { <<
983 if (fSurfaceArea == 0.) <<
984 { <<
985 G4double dphi = GetDPhi(); <<
986 G4double Ainn = GetInnerRadius(); <<
987 G4double Aout = GetOuterRadius(); <<
988 G4double Rinn0 = GetEndInnerRadius(0); <<
989 G4double Rout0 = GetEndOuterRadius(0); <<
990 G4double Rinn1 = GetEndInnerRadius(1); <<
991 G4double Rout1 = GetEndOuterRadius(1); <<
992 G4double z0 = GetEndZ(0); <<
993 G4double z1 = GetEndZ(1); <<
994 <<
995 G4double base0 = 0.5*dphi*(Rout0*Rout0 - R <<
996 G4double inner0 = GetLateralArea(Ainn, Rin <<
997 G4double outer0 = GetLateralArea(Aout, Rou <<
998 G4double cut0 = <<
999 GetPhiCutArea(Aout, Rout0, z0) - GetPhiC <<
1000 <<
1001 G4double base1 = base0; <<
1002 G4double inner1 = inner0; <<
1003 G4double outer1 = outer0; <<
1004 G4double cut1 = cut0; <<
1005 if (std::abs(z0) != std::abs(z1)) <<
1006 { <<
1007 base1 = 0.5*dphi*(Rout1*Rout1 - Rinn1*R <<
1008 inner1 = GetLateralArea(Ainn, Rinn1, z1 <<
1009 outer1 = GetLateralArea(Aout, Rout1, z1 <<
1010 cut1 = <<
1011 GetPhiCutArea(Aout, Rout1, z1) - GetPhi <<
1012 } <<
1013 fSurfaceArea = base0 + base1 + <<
1014 ((z0*z1 < 0) ? <<
1015 (inner0 + inner1 + outer0 + outer1 + 2. <<
1016 std::abs(inner0 - inner1 + outer0 - out <<
1017 } <<
1018 return fSurfaceArea; <<
1019 } <<
1020 <<
1021 //=========================================== <<
1022 //* GetPointOnSurface ----------------------- <<
1023 <<
1024 G4ThreeVector G4TwistedTubs::GetPointOnSurfac <<
1025 { <<
1026 <<
1027 G4double z = G4RandFlat::shoot(fEndZ[0],fEn <<
1028 G4double phi , phimin, phimax ; <<
1029 G4double x , xmin, xmax ; <<
1030 G4double r , rmin, rmax ; <<
1031 <<
1032 G4double a1 = fOuterHype->GetSurfaceArea() <<
1033 G4double a2 = fInnerHype->GetSurfaceArea() <<
1034 G4double a3 = fLatterTwisted->GetSurfaceAre <<
1035 G4double a4 = fFormerTwisted->GetSurfaceAre <<
1036 G4double a5 = fLowerEndcap->GetSurfaceArea( <<
1037 G4double a6 = fUpperEndcap->GetSurfaceArea( <<
1038 <<
1039 G4double chose = G4RandFlat::shoot(0.,a1 + <<
1040 <<
1041 if(chose < a1) <<
1042 { <<
1043 <<
1044 phimin = fOuterHype->GetBoundaryMin(z) ; <<
1045 phimax = fOuterHype->GetBoundaryMax(z) ; <<
1046 phi = G4RandFlat::shoot(phimin,phimax) ; <<
1047 <<
1048 return fOuterHype->SurfacePoint(phi,z,tru <<
1049 <<
1050 } <<
1051 else if ( (chose >= a1) && (chose < a1 + a2 <<
1052 { <<
1053 <<
1054 phimin = fInnerHype->GetBoundaryMin(z) ; <<
1055 phimax = fInnerHype->GetBoundaryMax(z) ; <<
1056 phi = G4RandFlat::shoot(phimin,phimax) ; <<
1057 <<
1058 return fInnerHype->SurfacePoint(phi,z,tru <<
1059 <<
1060 } <<
1061 else if ( (chose >= a1 + a2 ) && (chose < a <<
1062 { <<
1063 <<
1064 xmin = fLatterTwisted->GetBoundaryMin(z) <<
1065 xmax = fLatterTwisted->GetBoundaryMax(z) <<
1066 x = G4RandFlat::shoot(xmin,xmax) ; <<
1067 <<
1068 return fLatterTwisted->SurfacePoint(x,z,t <<
1069 <<
1070 } <<
1071 else if ( (chose >= a1 + a2 + a3 ) && (cho <<
1072 { <<
1073 <<
1074 xmin = fFormerTwisted->GetBoundaryMin(z) <<
1075 xmax = fFormerTwisted->GetBoundaryMax(z) <<
1076 x = G4RandFlat::shoot(xmin,xmax) ; <<
1077 <<
1078 return fFormerTwisted->SurfacePoint(x,z,t <<
1079 } <<
1080 else if( (chose >= a1 + a2 + a3 + a4 )&&(c <<
1081 { <<
1082 rmin = GetEndInnerRadius(0) ; <<
1083 rmax = GetEndOuterRadius(0) ; <<
1084 r = std::sqrt(G4RandFlat::shoot()*(sqr(rm <<
1085 <<
1086 phimin = fLowerEndcap->GetBoundaryMin(r) <<
1087 phimax = fLowerEndcap->GetBoundaryMax(r) <<
1088 phi = G4RandFlat::shoot(phimin,phimax) <<
1089 <<
1090 return fLowerEndcap->SurfacePoint(phi,r,t <<
1091 } <<
1092 else <<
1093 { <<
1094 rmin = GetEndInnerRadius(1) ; <<
1095 rmax = GetEndOuterRadius(1) ; <<
1096 r = rmin + (rmax-rmin)*std::sqrt(G4RandFl <<
1097 <<
1098 phimin = fUpperEndcap->GetBoundaryMin(r) <<
1099 phimax = fUpperEndcap->GetBoundaryMax(r) <<
1100 phi = G4RandFlat::shoot(phimin,phimax) <<
1101 <<
1102 return fUpperEndcap->SurfacePoint(phi,r,t <<
1103 } <<
1104 } 1043 }
1105 1044