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