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