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