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