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