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