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