| 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-patch-02 $
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 << 394
312 G4double totalphi = ephi - sphi; << 395 //
>> 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)
>> 411 {
>> 412 v4 = transform.TransformPoint(
>> 413 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, 0));
>> 414 }
>> 415
>> 416 // Inner hyperbolic surface
313 417
314 // Find bounding box << 418 G4double zInnerSplit = 0.;
315 if (dphi <= 0 || totalphi >= CLHEP::twopi) << 419 if (splitInner)
316 { 420 {
317 pMin.set(-rmax,-rmax, zmin); << 421 v2 = transform.TransformPoint(
318 pMax.set( rmax, rmax, zmax); << 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
>> 477 {
>> 478 AddPolyToExtent( v0, v1, w1, w0, voxelLimit, axis, extentList );
>> 479 }
>> 480
>> 481 // Inner hyperbolic surface
>> 482
>> 483 if (splitInner)
331 { 484 {
332 std::ostringstream message; << 485 w2 = transform.TransformPoint(
333 message << "Bad bounding box (min >= max) << 486 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
334 << GetName() << " !" << 487 + fZHalfLength));
335 << "\npMin = " << pMin << 488 w3 = transform.TransformPoint(
336 << "\npMax = " << pMax; << 489 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
337 G4Exception("G4TwistedTubs::BoundingLimits << 490 - fZHalfLength));
338 JustWarning, message); << 491
339 DumpInfo(); << 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
340 } 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;
>> 543
>> 544 phiPoly.AddVertexInOrder( v0 );
>> 545 phiPoly.AddVertexInOrder( v1 );
>> 546 phiPoly.AddVertexInOrder( w1 );
>> 547 phiPoly.AddVertexInOrder( w0 );
345 548
346 G4bool << 549 if (phiPoly.PartialClip( voxelLimit, axis ))
347 G4TwistedTubs::CalculateExtent(const EAxis pAx << 550 {
348 const G4VoxelLi << 551 phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
349 const G4AffineT << 552 extentList.AddSurface( phiPoly );
350 G4double& << 553 }
351 { <<
352 G4ThreeVector bmin, bmax; <<
353 <<
354 // Get bounding box <<
355 BoundingLimits(bmin,bmax); <<
356 <<
357 // Find extent <<
358 G4BoundingEnvelope bbox(bmin,bmax); <<
359 return bbox.CalculateExtent(pAxis,pVoxelLimi <<
360 } 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
>> 569
374 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 for (i=0; i< 6; i++)
499 { 744 {
500 G4double tmpdistance = surface->Distance << 745 G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
501 if (tmpdistance < distance) 746 if (tmpdistance < distance)
502 { 747 {
503 distance = tmpdistance; 748 distance = tmpdistance;
504 bestxx = xx; 749 bestxx = xx;
505 } 750 }
506 } 751 }
507 return distance; << 752 *tmpdist = distance;
>> 753
>> 754 return fLastDistanceToInWithV.value;
508 } 755 }
509 756
510 //============================================ 757 //=====================================================================
511 //* DistanceToIn (p) ------------------------- 758 //* DistanceToIn (p) --------------------------------------------------
512 759
513 G4double G4TwistedTubs::DistanceToIn (const G4 760 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p) const
514 { 761 {
515 // DistanceToIn(p): 762 // DistanceToIn(p):
516 // Calculate distance to surface of shape f 763 // Calculate distance to surface of shape from `outside',
517 // allowing for tolerance 764 // allowing for tolerance
>> 765
>> 766 //
>> 767 // checking last value
>> 768 //
>> 769
>> 770 G4ThreeVector *tmpp;
>> 771 G4double *tmpdist;
>> 772 if (fLastDistanceToIn.p == p)
>> 773 {
>> 774 return fLastDistanceToIn.value;
>> 775 }
>> 776 else
>> 777 {
>> 778 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
>> 779 tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
>> 780 tmpp->set(p.x(), p.y(), p.z());
>> 781 }
518 782
519 // 783 //
520 // Calculate DistanceToIn(p) 784 // Calculate DistanceToIn(p)
521 // 785 //
522 786
523 EInside currentside = Inside(p); 787 EInside currentside = Inside(p);
524 788
525 switch (currentside) 789 switch (currentside)
526 { 790 {
527 case (kInside) : 791 case (kInside) :
528 {} 792 {}
529 case (kSurface) : 793 case (kSurface) :
530 { 794 {
531 return 0; << 795 *tmpdist = 0.;
>> 796 return fLastDistanceToIn.value;
532 } 797 }
533 case (kOutside) : 798 case (kOutside) :
534 { 799 {
535 // Initialize 800 // Initialize
536 G4double distance = kInfinity; << 801 G4double distance = kInfinity;
537 802
538 // find intersections and choose near 803 // find intersections and choose nearest one.
539 G4VTwistSurface *surfaces[6]; 804 G4VTwistSurface *surfaces[6];
540 surfaces[0] = fLowerEndcap; 805 surfaces[0] = fLowerEndcap;
541 surfaces[1] = fUpperEndcap; 806 surfaces[1] = fUpperEndcap;
542 surfaces[2] = fLatterTwisted; 807 surfaces[2] = fLatterTwisted;
543 surfaces[3] = fFormerTwisted; 808 surfaces[3] = fFormerTwisted;
544 surfaces[4] = fInnerHype; 809 surfaces[4] = fInnerHype;
545 surfaces[5] = fOuterHype; 810 surfaces[5] = fOuterHype;
546 811
>> 812 G4int i;
547 G4ThreeVector xx; 813 G4ThreeVector xx;
548 G4ThreeVector bestxx; 814 G4ThreeVector bestxx;
549 for (const auto & surface : surfaces) << 815 for (i=0; i< 6; i++)
550 { 816 {
551 G4double tmpdistance = surface->Di << 817 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
552 if (tmpdistance < distance) 818 if (tmpdistance < distance)
553 { 819 {
554 distance = tmpdistance; 820 distance = tmpdistance;
555 bestxx = xx; 821 bestxx = xx;
556 } 822 }
557 } 823 }
558 return distance; << 824 *tmpdist = distance;
>> 825 return fLastDistanceToIn.value;
559 } 826 }
560 default : 827 default :
561 { 828 {
562 G4Exception("G4TwistedTubs::DistanceT << 829 G4Exception("G4TwistedTubs::DistanceToIn(p)", "InvalidCondition",
563 FatalException, "Unknown 830 FatalException, "Unknown point location!");
564 } 831 }
565 } // switch end 832 } // switch end
566 833
567 return kInfinity; 834 return kInfinity;
568 } 835 }
569 836
570 //============================================ 837 //=====================================================================
571 //* DistanceToOut (p, v) --------------------- 838 //* DistanceToOut (p, v) ----------------------------------------------
572 839
573 G4double G4TwistedTubs::DistanceToOut( const G 840 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p,
574 const G 841 const G4ThreeVector& v,
575 const G 842 const G4bool calcNorm,
576 G4bool 843 G4bool *validNorm,
577 G4Three 844 G4ThreeVector *norm ) const
578 { 845 {
579 // DistanceToOut (p, v): 846 // DistanceToOut (p, v):
580 // Calculate distance to surface of shape f 847 // Calculate distance to surface of shape from `inside'
581 // along with the v, allowing for tolerance 848 // along with the v, allowing for tolerance.
582 // The function returns kInfinity if no int 849 // The function returns kInfinity if no intersection or
583 // just grazing within tolerance. 850 // just grazing within tolerance.
584 851
585 // 852 //
>> 853 // checking last value
>> 854 //
>> 855
>> 856 G4ThreeVector *tmpp;
>> 857 G4ThreeVector *tmpv;
>> 858 G4double *tmpdist;
>> 859 if ((fLastDistanceToOutWithV.p == p) && (fLastDistanceToOutWithV.vec == v) )
>> 860 {
>> 861 return fLastDistanceToOutWithV.value;
>> 862 }
>> 863 else
>> 864 {
>> 865 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
>> 866 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
>> 867 tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
>> 868 tmpp->set(p.x(), p.y(), p.z());
>> 869 tmpv->set(v.x(), v.y(), v.z());
>> 870 }
>> 871
>> 872 //
586 // Calculate DistanceToOut(p,v) 873 // Calculate DistanceToOut(p,v)
587 // 874 //
588 875
589 EInside currentside = Inside(p); 876 EInside currentside = Inside(p);
>> 877
590 if (currentside == kOutside) 878 if (currentside == kOutside)
591 { 879 {
592 } 880 }
593 else 881 else
594 { 882 {
595 if (currentside == kSurface) 883 if (currentside == kSurface)
596 { 884 {
597 // particle is just on a boundary. 885 // particle is just on a boundary.
598 // If the particle is exiting from the 886 // If the particle is exiting from the volume, return 0.
599 // 887 //
600 G4ThreeVector normal = SurfaceNormal(p) 888 G4ThreeVector normal = SurfaceNormal(p);
>> 889 G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
601 if (normal*v > 0) 890 if (normal*v > 0)
602 { 891 {
603 if (calcNorm) 892 if (calcNorm)
604 { 893 {
605 *norm = normal; << 894 *norm = (blockedsurface->GetNormal(p, true));
606 *validNorm = true; << 895 *validNorm = blockedsurface->IsValidNorm();
607 } 896 }
608 return 0; << 897 *tmpdist = 0.;
>> 898 return fLastDistanceToOutWithV.value;
609 } 899 }
610 } 900 }
611 } 901 }
612 902
613 // now, we can take smallest positive dista 903 // now, we can take smallest positive distance.
614 904
615 // Initialize 905 // Initialize
616 // 906 //
617 G4double distance = kInfinity; << 907 G4double distance = kInfinity;
618 908
619 // find intersections and choose nearest on 909 // find intersections and choose nearest one.
620 // 910 //
621 G4VTwistSurface* surfaces[6]; << 911 G4VTwistSurface *surfaces[6];
622 surfaces[0] = fLatterTwisted; 912 surfaces[0] = fLatterTwisted;
623 surfaces[1] = fFormerTwisted; 913 surfaces[1] = fFormerTwisted;
624 surfaces[2] = fInnerHype; 914 surfaces[2] = fInnerHype;
625 surfaces[3] = fOuterHype; 915 surfaces[3] = fOuterHype;
626 surfaces[4] = fLowerEndcap; 916 surfaces[4] = fLowerEndcap;
627 surfaces[5] = fUpperEndcap; 917 surfaces[5] = fUpperEndcap;
628 918
>> 919 G4int i;
629 G4int besti = -1; 920 G4int besti = -1;
630 G4ThreeVector xx; 921 G4ThreeVector xx;
631 G4ThreeVector bestxx; 922 G4ThreeVector bestxx;
632 for (auto i=0; i<6; ++i) << 923 for (i=0; i< 6; i++)
633 { 924 {
634 G4double tmpdistance = surfaces[i]->Dist 925 G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
635 if (tmpdistance < distance) 926 if (tmpdistance < distance)
636 { 927 {
637 distance = tmpdistance; 928 distance = tmpdistance;
638 bestxx = xx; 929 bestxx = xx;
639 besti = i; 930 besti = i;
640 } 931 }
641 } 932 }
642 933
643 if (calcNorm) 934 if (calcNorm)
644 { 935 {
645 if (besti != -1) 936 if (besti != -1)
646 { 937 {
647 *norm = (surfaces[besti]->GetNormal(p 938 *norm = (surfaces[besti]->GetNormal(p, true));
648 *validNorm = surfaces[besti]->IsValid 939 *validNorm = surfaces[besti]->IsValidNorm();
649 } 940 }
650 } 941 }
651 942
652 return distance; << 943 *tmpdist = distance;
>> 944
>> 945 return fLastDistanceToOutWithV.value;
653 } 946 }
654 947
655 948
656 //============================================ 949 //=====================================================================
657 //* DistanceToOut (p) ------------------------ 950 //* DistanceToOut (p) ----------------------------------------------
658 951
659 G4double G4TwistedTubs::DistanceToOut( const G 952 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p ) const
660 { 953 {
661 // DistanceToOut(p): 954 // DistanceToOut(p):
662 // Calculate distance to surface of shape f 955 // Calculate distance to surface of shape from `inside',
663 // allowing for tolerance 956 // allowing for tolerance
664 957
665 // 958 //
>> 959 // checking last value
>> 960 //
>> 961
>> 962 G4ThreeVector *tmpp;
>> 963 G4double *tmpdist;
>> 964 if (fLastDistanceToOut.p == p)
>> 965 {
>> 966 return fLastDistanceToOut.value;
>> 967 }
>> 968 else
>> 969 {
>> 970 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
>> 971 tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
>> 972 tmpp->set(p.x(), p.y(), p.z());
>> 973 }
>> 974
>> 975 //
666 // Calculate DistanceToOut(p) 976 // Calculate DistanceToOut(p)
667 // 977 //
668 978
669 EInside currentside = Inside(p); 979 EInside currentside = Inside(p);
670 980
671 switch (currentside) 981 switch (currentside)
672 { 982 {
673 case (kOutside) : 983 case (kOutside) :
674 { 984 {
675 } 985 }
676 case (kSurface) : 986 case (kSurface) :
677 { 987 {
678 return 0; << 988 *tmpdist = 0.;
>> 989 return fLastDistanceToOut.value;
679 } 990 }
680 case (kInside) : 991 case (kInside) :
681 { 992 {
682 // Initialize 993 // Initialize
683 G4double distance = kInfinity; 994 G4double distance = kInfinity;
684 995
685 // find intersections and choose near 996 // find intersections and choose nearest one.
686 G4VTwistSurface* surfaces[6]; << 997 G4VTwistSurface *surfaces[6];
687 surfaces[0] = fLatterTwisted; 998 surfaces[0] = fLatterTwisted;
688 surfaces[1] = fFormerTwisted; 999 surfaces[1] = fFormerTwisted;
689 surfaces[2] = fInnerHype; 1000 surfaces[2] = fInnerHype;
690 surfaces[3] = fOuterHype; 1001 surfaces[3] = fOuterHype;
691 surfaces[4] = fLowerEndcap; 1002 surfaces[4] = fLowerEndcap;
692 surfaces[5] = fUpperEndcap; 1003 surfaces[5] = fUpperEndcap;
693 1004
>> 1005 G4int i;
694 G4ThreeVector xx; 1006 G4ThreeVector xx;
695 G4ThreeVector bestxx; 1007 G4ThreeVector bestxx;
696 for (const auto & surface : surfaces) << 1008 for (i=0; i< 6; i++)
697 { 1009 {
698 G4double tmpdistance = surface->Di << 1010 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
699 if (tmpdistance < distance) 1011 if (tmpdistance < distance)
700 { 1012 {
701 distance = tmpdistance; 1013 distance = tmpdistance;
702 bestxx = xx; 1014 bestxx = xx;
703 } 1015 }
704 } 1016 }
705 return distance; << 1017 *tmpdist = distance;
>> 1018
>> 1019 return fLastDistanceToOut.value;
706 } 1020 }
707 default : 1021 default :
708 { 1022 {
709 G4Exception("G4TwistedTubs::DistanceT << 1023 G4Exception("G4TwistedTubs::DistanceToOut(p)", "InvalidCondition",
710 FatalException, "Unknown 1024 FatalException, "Unknown point location!");
711 } 1025 }
712 } // switch end 1026 } // switch end
713 1027
714 return 0.; << 1028 return 0;
715 } 1029 }
716 1030
717 //============================================ 1031 //=====================================================================
718 //* StreamInfo ------------------------------- 1032 //* StreamInfo --------------------------------------------------------
719 1033
720 std::ostream& G4TwistedTubs::StreamInfo(std::o 1034 std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
721 { 1035 {
722 // 1036 //
723 // Stream object contents to an output strea 1037 // Stream object contents to an output stream
724 // 1038 //
725 G4long oldprc = os.precision(16); <<
726 os << "------------------------------------- 1039 os << "-----------------------------------------------------------\n"
727 << " *** Dump for solid - " << GetName 1040 << " *** Dump for solid - " << GetName() << " ***\n"
728 << " ================================= 1041 << " ===================================================\n"
729 << " Solid type: G4TwistedTubs\n" 1042 << " Solid type: G4TwistedTubs\n"
730 << " Parameters: \n" 1043 << " Parameters: \n"
731 << " -ve end Z : " << fEn 1044 << " -ve end Z : " << fEndZ[0]/mm << " mm \n"
732 << " +ve end Z : " << fEn 1045 << " +ve end Z : " << fEndZ[1]/mm << " mm \n"
733 << " inner end radius(-ve z): " << fEn 1046 << " inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
734 << " inner end radius(+ve z): " << fEn 1047 << " inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
735 << " outer end radius(-ve z): " << fEn 1048 << " outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
736 << " outer end radius(+ve z): " << fEn 1049 << " outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
737 << " inner radius (z=0) : " << fIn 1050 << " inner radius (z=0) : " << fInnerRadius/mm << " mm \n"
738 << " outer radius (z=0) : " << fOu 1051 << " outer radius (z=0) : " << fOuterRadius/mm << " mm \n"
739 << " twisted angle : " << fPh 1052 << " twisted angle : " << fPhiTwist/degree << " degrees \n"
740 << " inner stereo angle : " << fIn 1053 << " inner stereo angle : " << fInnerStereo/degree << " degrees \n"
741 << " outer stereo angle : " << fOu 1054 << " outer stereo angle : " << fOuterStereo/degree << " degrees \n"
742 << " phi-width of a piece : " << fDP 1055 << " phi-width of a piece : " << fDPhi/degree << " degrees \n"
743 << "------------------------------------- 1056 << "-----------------------------------------------------------\n";
744 os.precision(oldprc); <<
745 1057
746 return os; 1058 return os;
747 } 1059 }
748 1060
749 1061
750 //============================================ 1062 //=====================================================================
751 //* DiscribeYourselfTo ----------------------- 1063 //* DiscribeYourselfTo ------------------------------------------------
752 1064
753 void G4TwistedTubs::DescribeYourselfTo (G4VGra 1065 void G4TwistedTubs::DescribeYourselfTo (G4VGraphicsScene& scene) const
754 { 1066 {
755 scene.AddSolid (*this); 1067 scene.AddSolid (*this);
756 } 1068 }
757 1069
758 //============================================ 1070 //=====================================================================
759 //* GetExtent -------------------------------- 1071 //* GetExtent ---------------------------------------------------------
760 1072
761 G4VisExtent G4TwistedTubs::GetExtent() const << 1073 G4VisExtent G4TwistedTubs::GetExtent() const
762 { 1074 {
763 // Define the sides of the box into which th 1075 // Define the sides of the box into which the G4Tubs instance would fit.
764 // << 1076
765 G4ThreeVector pmin,pmax; << 1077 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
766 BoundingLimits(pmin,pmax); << 1078 return G4VisExtent( -maxEndOuterRad, maxEndOuterRad,
767 return { pmin.x(),pmax.x(), << 1079 -maxEndOuterRad, maxEndOuterRad,
768 pmin.y(),pmax.y(), << 1080 -fZHalfLength, fZHalfLength );
769 pmin.z(),pmax.z() }; <<
770 } 1081 }
771 1082
772 //============================================ 1083 //=====================================================================
773 //* CreatePolyhedron ------------------------- 1084 //* CreatePolyhedron --------------------------------------------------
774 1085
775 G4Polyhedron* G4TwistedTubs::CreatePolyhedron 1086 G4Polyhedron* G4TwistedTubs::CreatePolyhedron () const
776 { 1087 {
777 // number of meshes 1088 // number of meshes
778 // 1089 //
779 G4double absPhiTwist = std::abs(fPhiTwist); << 1090 G4double dA = std::max(fDPhi,fPhiTwist);
780 G4double dA = std::max(fDPhi,absPhiTwist); << 1091 const G4int m =
781 const G4int k = <<
782 G4int(G4Polyhedron::GetNumberOfRotationSte 1092 G4int(G4Polyhedron::GetNumberOfRotationSteps() * dA / twopi) + 2;
783 const G4int n = 1093 const G4int n =
784 G4int(G4Polyhedron::GetNumberOfRotationSte << 1094 G4int(G4Polyhedron::GetNumberOfRotationSteps() * fPhiTwist / twopi) + 2;
785 1095
786 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ; << 1096 const G4int nnodes = 4*(m-1)*(n-2) + 2*m*m ;
787 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1) << 1097 const G4int nfaces = 4*(m-1)*(n-1) + 2*(m-1)*(m-1) ;
788 1098
789 auto ph = new G4Polyhedron; << 1099 G4Polyhedron *ph=new G4Polyhedron;
790 typedef G4double G4double3[3]; 1100 typedef G4double G4double3[3];
791 typedef G4int G4int4[4]; 1101 typedef G4int G4int4[4];
792 auto xyz = new G4double3[nnodes]; // number << 1102 G4double3* xyz = new G4double3[nnodes]; // number of nodes
793 auto faces = new G4int4[nfaces] ; // number << 1103 G4int4* faces = new G4int4[nfaces] ; // number of faces
794 fLowerEndcap->GetFacets(k,k,xyz,faces,0) ; << 1104 fLowerEndcap->GetFacets(m,m,xyz,faces,0) ;
795 fUpperEndcap->GetFacets(k,k,xyz,faces,1) ; << 1105 fUpperEndcap->GetFacets(m,m,xyz,faces,1) ;
796 fInnerHype->GetFacets(k,n,xyz,faces,2) ; << 1106 fInnerHype->GetFacets(m,n,xyz,faces,2) ;
797 fFormerTwisted->GetFacets(k,n,xyz,faces,3) ; << 1107 fFormerTwisted->GetFacets(m,n,xyz,faces,3) ;
798 fOuterHype->GetFacets(k,n,xyz,faces,4) ; << 1108 fOuterHype->GetFacets(m,n,xyz,faces,4) ;
799 fLatterTwisted->GetFacets(k,n,xyz,faces,5) ; << 1109 fLatterTwisted->GetFacets(m,n,xyz,faces,5) ;
800 1110
801 ph->createPolyhedron(nnodes,nfaces,xyz,faces 1111 ph->createPolyhedron(nnodes,nfaces,xyz,faces);
802 1112
803 delete[] xyz; 1113 delete[] xyz;
804 delete[] faces; 1114 delete[] faces;
805 1115
806 return ph; 1116 return ph;
807 } 1117 }
808 1118
809 //============================================ 1119 //=====================================================================
>> 1120 //* CreateNUBS --------------------------------------------------------
>> 1121
>> 1122 G4NURBS* G4TwistedTubs::CreateNURBS () const
>> 1123 {
>> 1124 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
>> 1125 G4double maxEndInnerRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
>> 1126 return new G4NURBStube(maxEndInnerRad, maxEndOuterRad, fZHalfLength);
>> 1127 // Tube for now!!!
>> 1128 }
>> 1129
>> 1130 //=====================================================================
810 //* GetPolyhedron ---------------------------- 1131 //* GetPolyhedron -----------------------------------------------------
811 1132
812 G4Polyhedron* G4TwistedTubs::GetPolyhedron () 1133 G4Polyhedron* G4TwistedTubs::GetPolyhedron () const
813 { 1134 {
814 if (fpPolyhedron == nullptr || << 1135 if ((!fpPolyhedron) ||
815 fRebuildPolyhedron || << 1136 (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
816 fpPolyhedron->GetNumberOfRotationStepsAt << 1137 fpPolyhedron->GetNumberOfRotationSteps()))
817 fpPolyhedron->GetNumberOfRotationSteps() <<
818 { 1138 {
819 G4AutoLock l(&polyhedronMutex); <<
820 delete fpPolyhedron; 1139 delete fpPolyhedron;
821 fpPolyhedron = CreatePolyhedron(); 1140 fpPolyhedron = CreatePolyhedron();
822 fRebuildPolyhedron = false; <<
823 l.unlock(); <<
824 } 1141 }
825 return fpPolyhedron; 1142 return fpPolyhedron;
826 } 1143 }
827 1144
828 //============================================ 1145 //=====================================================================
829 //* CreateSurfaces --------------------------- 1146 //* CreateSurfaces ----------------------------------------------------
830 1147
831 void G4TwistedTubs::CreateSurfaces() 1148 void G4TwistedTubs::CreateSurfaces()
832 { 1149 {
833 // create 6 surfaces of TwistedTub 1150 // create 6 surfaces of TwistedTub
834 1151
835 G4ThreeVector x0(0, 0, fEndZ[0]); 1152 G4ThreeVector x0(0, 0, fEndZ[0]);
836 G4ThreeVector n (0, 0, -1); 1153 G4ThreeVector n (0, 0, -1);
837 1154
838 fLowerEndcap = new G4TwistTubsFlatSide("Low 1155 fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap",
839 fEndInnerR 1156 fEndInnerRadius, fEndOuterRadius,
840 fDPhi, fEn 1157 fDPhi, fEndPhi, fEndZ, -1) ;
841 1158
842 fUpperEndcap = new G4TwistTubsFlatSide("Upp 1159 fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap",
843 fEndInnerR 1160 fEndInnerRadius, fEndOuterRadius,
844 fDPhi, fEn 1161 fDPhi, fEndPhi, fEndZ, 1) ;
845 1162
846 G4RotationMatrix rotHalfDPhi; 1163 G4RotationMatrix rotHalfDPhi;
847 rotHalfDPhi.rotateZ(0.5*fDPhi); 1164 rotHalfDPhi.rotateZ(0.5*fDPhi);
848 1165
849 fLatterTwisted = new G4TwistTubsSide("Latte 1166 fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
850 fEndI 1167 fEndInnerRadius, fEndOuterRadius,
851 fDPhi 1168 fDPhi, fEndPhi, fEndZ,
852 fInne 1169 fInnerRadius, fOuterRadius, fKappa,
853 1 ) ; 1170 1 ) ;
854 fFormerTwisted = new G4TwistTubsSide("Forme 1171 fFormerTwisted = new G4TwistTubsSide("FormerTwisted",
855 fEndI 1172 fEndInnerRadius, fEndOuterRadius,
856 fDPhi 1173 fDPhi, fEndPhi, fEndZ,
857 fInne 1174 fInnerRadius, fOuterRadius, fKappa,
858 -1 ) 1175 -1 ) ;
859 1176
860 fInnerHype = new G4TwistTubsHypeSide("Inner 1177 fInnerHype = new G4TwistTubsHypeSide("InnerHype",
861 fEndIn 1178 fEndInnerRadius, fEndOuterRadius,
862 fDPhi, 1179 fDPhi, fEndPhi, fEndZ,
863 fInner 1180 fInnerRadius, fOuterRadius,fKappa,
864 fTanIn 1181 fTanInnerStereo, fTanOuterStereo, -1) ;
865 fOuterHype = new G4TwistTubsHypeSide("Outer 1182 fOuterHype = new G4TwistTubsHypeSide("OuterHype",
866 fEndIn 1183 fEndInnerRadius, fEndOuterRadius,
867 fDPhi, 1184 fDPhi, fEndPhi, fEndZ,
868 fInner 1185 fInnerRadius, fOuterRadius,fKappa,
869 fTanIn 1186 fTanInnerStereo, fTanOuterStereo, 1) ;
870 1187
871 1188
872 // set neighbour surfaces 1189 // set neighbour surfaces
873 // 1190 //
874 fLowerEndcap->SetNeighbours(fInnerHype, fLa 1191 fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
875 fOuterHype, fFo 1192 fOuterHype, fFormerTwisted);
876 fUpperEndcap->SetNeighbours(fInnerHype, fLa 1193 fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
877 fOuterHype, fFo 1194 fOuterHype, fFormerTwisted);
878 fLatterTwisted->SetNeighbours(fInnerHype, f 1195 fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
879 fOuterHype, f 1196 fOuterHype, fUpperEndcap);
880 fFormerTwisted->SetNeighbours(fInnerHype, f 1197 fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
881 fOuterHype, f 1198 fOuterHype, fUpperEndcap);
882 fInnerHype->SetNeighbours(fLatterTwisted, f 1199 fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
883 fFormerTwisted, f 1200 fFormerTwisted, fUpperEndcap);
884 fOuterHype->SetNeighbours(fLatterTwisted, f 1201 fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
885 fFormerTwisted, f 1202 fFormerTwisted, fUpperEndcap);
886 } 1203 }
887 1204
888 1205
889 //============================================ 1206 //=====================================================================
890 //* GetEntityType ---------------------------- 1207 //* GetEntityType -----------------------------------------------------
891 1208
892 G4GeometryType G4TwistedTubs::GetEntityType() 1209 G4GeometryType G4TwistedTubs::GetEntityType() const
893 { 1210 {
894 return {"G4TwistedTubs"}; << 1211 return G4String("G4TwistedTubs");
895 } 1212 }
896 1213
897 //============================================ 1214 //=====================================================================
898 //* Clone ------------------------------------ 1215 //* Clone -------------------------------------------------------------
899 1216
900 G4VSolid* G4TwistedTubs::Clone() const 1217 G4VSolid* G4TwistedTubs::Clone() const
901 { 1218 {
902 return new G4TwistedTubs(*this); 1219 return new G4TwistedTubs(*this);
903 } 1220 }
904 1221
905 //============================================ 1222 //=====================================================================
906 //* GetCubicVolume --------------------------- 1223 //* GetCubicVolume ----------------------------------------------------
907 1224
908 G4double G4TwistedTubs::GetCubicVolume() 1225 G4double G4TwistedTubs::GetCubicVolume()
909 { 1226 {
910 if (fCubicVolume == 0.) << 1227 if(fCubicVolume != 0.) {;}
911 { << 1228 else { fCubicVolume = fDPhi*fZHalfLength*(fOuterRadius*fOuterRadius
912 G4double DPhi = GetDPhi(); << 1229 -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; 1230 return fCubicVolume;
928 } 1231 }
929 1232
930 //============================================ 1233 //=====================================================================
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 --------------------------- 1234 //* GetSurfaceArea ----------------------------------------------------
980 1235
981 G4double G4TwistedTubs::GetSurfaceArea() 1236 G4double G4TwistedTubs::GetSurfaceArea()
982 { 1237 {
983 if (fSurfaceArea == 0.) << 1238 if(fSurfaceArea != 0.) {;}
984 { << 1239 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; 1240 return fSurfaceArea;
1019 } 1241 }
1020 1242
1021 //=========================================== 1243 //=====================================================================
1022 //* GetPointOnSurface ----------------------- 1244 //* GetPointOnSurface -------------------------------------------------
1023 1245
1024 G4ThreeVector G4TwistedTubs::GetPointOnSurfac 1246 G4ThreeVector G4TwistedTubs::GetPointOnSurface() const
1025 { 1247 {
1026 1248
1027 G4double z = G4RandFlat::shoot(fEndZ[0],fEn << 1249 G4double z = G4RandFlat::shoot(fEndZ[0],fEndZ[1]);
1028 G4double phi , phimin, phimax ; 1250 G4double phi , phimin, phimax ;
1029 G4double x , xmin, xmax ; 1251 G4double x , xmin, xmax ;
1030 G4double r , rmin, rmax ; 1252 G4double r , rmin, rmax ;
1031 1253
1032 G4double a1 = fOuterHype->GetSurfaceArea() 1254 G4double a1 = fOuterHype->GetSurfaceArea() ;
1033 G4double a2 = fInnerHype->GetSurfaceArea() 1255 G4double a2 = fInnerHype->GetSurfaceArea() ;
1034 G4double a3 = fLatterTwisted->GetSurfaceAre 1256 G4double a3 = fLatterTwisted->GetSurfaceArea() ;
1035 G4double a4 = fFormerTwisted->GetSurfaceAre 1257 G4double a4 = fFormerTwisted->GetSurfaceArea() ;
1036 G4double a5 = fLowerEndcap->GetSurfaceArea( 1258 G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1037 G4double a6 = fUpperEndcap->GetSurfaceArea( 1259 G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1038 1260
1039 G4double chose = G4RandFlat::shoot(0.,a1 + 1261 G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1040 1262
1041 if(chose < a1) 1263 if(chose < a1)
1042 { 1264 {
1043 1265
1044 phimin = fOuterHype->GetBoundaryMin(z) ; 1266 phimin = fOuterHype->GetBoundaryMin(z) ;
1045 phimax = fOuterHype->GetBoundaryMax(z) ; 1267 phimax = fOuterHype->GetBoundaryMax(z) ;
1046 phi = G4RandFlat::shoot(phimin,phimax) ; 1268 phi = G4RandFlat::shoot(phimin,phimax) ;
1047 1269
1048 return fOuterHype->SurfacePoint(phi,z,tru 1270 return fOuterHype->SurfacePoint(phi,z,true) ;
1049 1271
1050 } 1272 }
1051 else if ( (chose >= a1) && (chose < a1 + a2 1273 else if ( (chose >= a1) && (chose < a1 + a2 ) )
1052 { 1274 {
1053 1275
1054 phimin = fInnerHype->GetBoundaryMin(z) ; 1276 phimin = fInnerHype->GetBoundaryMin(z) ;
1055 phimax = fInnerHype->GetBoundaryMax(z) ; 1277 phimax = fInnerHype->GetBoundaryMax(z) ;
1056 phi = G4RandFlat::shoot(phimin,phimax) ; 1278 phi = G4RandFlat::shoot(phimin,phimax) ;
1057 1279
1058 return fInnerHype->SurfacePoint(phi,z,tru 1280 return fInnerHype->SurfacePoint(phi,z,true) ;
1059 1281
1060 } 1282 }
1061 else if ( (chose >= a1 + a2 ) && (chose < a 1283 else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1062 { 1284 {
1063 1285
1064 xmin = fLatterTwisted->GetBoundaryMin(z) 1286 xmin = fLatterTwisted->GetBoundaryMin(z) ;
1065 xmax = fLatterTwisted->GetBoundaryMax(z) 1287 xmax = fLatterTwisted->GetBoundaryMax(z) ;
1066 x = G4RandFlat::shoot(xmin,xmax) ; 1288 x = G4RandFlat::shoot(xmin,xmax) ;
1067 1289
1068 return fLatterTwisted->SurfacePoint(x,z,t 1290 return fLatterTwisted->SurfacePoint(x,z,true) ;
1069 1291
1070 } 1292 }
1071 else if ( (chose >= a1 + a2 + a3 ) && (cho 1293 else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1072 { 1294 {
1073 1295
1074 xmin = fFormerTwisted->GetBoundaryMin(z) 1296 xmin = fFormerTwisted->GetBoundaryMin(z) ;
1075 xmax = fFormerTwisted->GetBoundaryMax(z) 1297 xmax = fFormerTwisted->GetBoundaryMax(z) ;
1076 x = G4RandFlat::shoot(xmin,xmax) ; 1298 x = G4RandFlat::shoot(xmin,xmax) ;
1077 1299
1078 return fFormerTwisted->SurfacePoint(x,z,t 1300 return fFormerTwisted->SurfacePoint(x,z,true) ;
1079 } << 1301
>> 1302 }
1080 else if( (chose >= a1 + a2 + a3 + a4 )&&(c 1303 else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
1081 { 1304 {
>> 1305
1082 rmin = GetEndInnerRadius(0) ; 1306 rmin = GetEndInnerRadius(0) ;
1083 rmax = GetEndOuterRadius(0) ; 1307 rmax = GetEndOuterRadius(0) ;
1084 r = std::sqrt(G4RandFlat::shoot()*(sqr(rm << 1308 r = G4RandFlat::shoot(rmin,rmax) ;
1085 1309
1086 phimin = fLowerEndcap->GetBoundaryMin(r) 1310 phimin = fLowerEndcap->GetBoundaryMin(r) ;
1087 phimax = fLowerEndcap->GetBoundaryMax(r) 1311 phimax = fLowerEndcap->GetBoundaryMax(r) ;
1088 phi = G4RandFlat::shoot(phimin,phimax) 1312 phi = G4RandFlat::shoot(phimin,phimax) ;
1089 1313
1090 return fLowerEndcap->SurfacePoint(phi,r,t 1314 return fLowerEndcap->SurfacePoint(phi,r,true) ;
>> 1315
1091 } 1316 }
1092 else 1317 else
1093 { 1318 {
1094 rmin = GetEndInnerRadius(1) ; 1319 rmin = GetEndInnerRadius(1) ;
1095 rmax = GetEndOuterRadius(1) ; 1320 rmax = GetEndOuterRadius(1) ;
1096 r = rmin + (rmax-rmin)*std::sqrt(G4RandFl << 1321 r = G4RandFlat::shoot(rmin,rmax) ;
1097 1322
1098 phimin = fUpperEndcap->GetBoundaryMin(r) 1323 phimin = fUpperEndcap->GetBoundaryMin(r) ;
1099 phimax = fUpperEndcap->GetBoundaryMax(r) 1324 phimax = fUpperEndcap->GetBoundaryMax(r) ;
1100 phi = G4RandFlat::shoot(phimin,phimax) 1325 phi = G4RandFlat::shoot(phimin,phimax) ;
1101 1326
1102 return fUpperEndcap->SurfacePoint(phi,r,t 1327 return fUpperEndcap->SurfacePoint(phi,r,true) ;
1103 } 1328 }
1104 } 1329 }
1105 1330