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