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