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