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