| 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 // G4VTwistedFaceted implementation << 26 // $Id$
>> 27 //
>> 28 //
>> 29 // --------------------------------------------------------------------
>> 30 // GEANT 4 class source file
>> 31 //
>> 32 //
>> 33 // G4VTwistedFaceted.cc
>> 34 //
>> 35 // Author:
>> 36 //
>> 37 // 04-Nov-2004 - O.Link (Oliver.Link@cern.ch)
27 // 38 //
28 // Author: 04-Nov-2004 - O.Link (Oliver.Link@c <<
29 // ------------------------------------------- 39 // --------------------------------------------------------------------
30 40
31 #include "G4VTwistedFaceted.hh" 41 #include "G4VTwistedFaceted.hh"
32 42
33 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh"
35 #include "G4VoxelLimits.hh" 45 #include "G4VoxelLimits.hh"
36 #include "G4AffineTransform.hh" 46 #include "G4AffineTransform.hh"
37 #include "G4BoundingEnvelope.hh" <<
38 #include "G4SolidExtentList.hh" 47 #include "G4SolidExtentList.hh"
39 #include "G4ClippablePolygon.hh" 48 #include "G4ClippablePolygon.hh"
40 #include "G4VPVParameterisation.hh" 49 #include "G4VPVParameterisation.hh"
41 #include "G4GeometryTolerance.hh" 50 #include "G4GeometryTolerance.hh"
42 #include "meshdefs.hh" 51 #include "meshdefs.hh"
43 52
44 #include "G4VGraphicsScene.hh" 53 #include "G4VGraphicsScene.hh"
45 #include "G4Polyhedron.hh" 54 #include "G4Polyhedron.hh"
46 #include "G4VisExtent.hh" 55 #include "G4VisExtent.hh"
>> 56 #include "G4NURBS.hh"
>> 57 #include "G4NURBStube.hh"
>> 58 #include "G4NURBScylinder.hh"
>> 59 #include "G4NURBStubesector.hh"
47 60
48 #include "Randomize.hh" 61 #include "Randomize.hh"
49 62
50 #include "G4AutoLock.hh" <<
51 <<
52 namespace <<
53 { <<
54 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE <<
55 } <<
56 <<
57 <<
58 //============================================ 63 //=====================================================================
59 //* constructors ----------------------------- 64 //* constructors ------------------------------------------------------
60 65
61 G4VTwistedFaceted:: 66 G4VTwistedFaceted::
62 G4VTwistedFaceted( const G4String& pname, << 67 G4VTwistedFaceted( const G4String &pname, // Name of instance
63 G4double PhiTwist, 68 G4double PhiTwist, // twist angle
64 G4double pDz, 69 G4double pDz, // half z length
65 G4double pTheta, // 70 G4double pTheta, // direction between end planes
66 G4double pPhi, // 71 G4double pPhi, // defined by polar and azim. angles
67 G4double pDy1, // 72 G4double pDy1, // half y length at -pDz
68 G4double pDx1, // 73 G4double pDx1, // half x length at -pDz,-pDy
69 G4double pDx2, // 74 G4double pDx2, // half x length at -pDz,+pDy
70 G4double pDy2, // 75 G4double pDy2, // half y length at +pDz
71 G4double pDx3, // 76 G4double pDx3, // half x length at +pDz,-pDy
72 G4double pDx4, // 77 G4double pDx4, // half x length at +pDz,+pDy
73 G4double pAlph ) // << 78 G4double pAlph // tilt angle
74 : G4VSolid(pname), << 79 )
75 fLowerEndcap(nullptr), fUpperEndcap(nullpt << 80 : G4VSolid(pname),
76 fSide90(nullptr), fSide180(nullptr), fSide << 81 fLowerEndcap(0), fUpperEndcap(0), fSide0(0),
>> 82 fSide90(0), fSide180(0), fSide270(0),
>> 83 fSurfaceArea(0.), fpPolyhedron(0)
77 { 84 {
78 85
79 G4double pDytmp ; 86 G4double pDytmp ;
80 G4double fDxUp ; 87 G4double fDxUp ;
81 G4double fDxDown ; 88 G4double fDxDown ;
82 89
83 fDx1 = pDx1 ; 90 fDx1 = pDx1 ;
84 fDx2 = pDx2 ; 91 fDx2 = pDx2 ;
85 fDx3 = pDx3 ; 92 fDx3 = pDx3 ;
86 fDx4 = pDx4 ; 93 fDx4 = pDx4 ;
87 fDy1 = pDy1 ; 94 fDy1 = pDy1 ;
88 fDy2 = pDy2 ; 95 fDy2 = pDy2 ;
89 fDz = pDz ; 96 fDz = pDz ;
90 97
91 G4double kAngTolerance 98 G4double kAngTolerance
92 = G4GeometryTolerance::GetInstance()->GetA 99 = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
93 100
94 // maximum values 101 // maximum values
95 // 102 //
96 fDxDown = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ; 103 fDxDown = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ;
97 fDxUp = ( fDx3 > fDx4 ? fDx3 : fDx4 ) ; 104 fDxUp = ( fDx3 > fDx4 ? fDx3 : fDx4 ) ;
98 fDx = ( fDxUp > fDxDown ? fDxUp : fDxDow 105 fDx = ( fDxUp > fDxDown ? fDxUp : fDxDown ) ;
99 fDy = ( fDy1 > fDy2 ? fDy1 : fDy2 ) ; 106 fDy = ( fDy1 > fDy2 ? fDy1 : fDy2 ) ;
100 107
101 // planarity check 108 // planarity check
102 // 109 //
103 if ( fDx1 != fDx2 && fDx3 != fDx4 ) 110 if ( fDx1 != fDx2 && fDx3 != fDx4 )
104 { 111 {
105 pDytmp = fDy1 * ( fDx3 - fDx4 ) / ( fDx1 - 112 pDytmp = fDy1 * ( fDx3 - fDx4 ) / ( fDx1 - fDx2 ) ;
106 if ( std::fabs(pDytmp - fDy2) > kCarTolera 113 if ( std::fabs(pDytmp - fDy2) > kCarTolerance )
107 { 114 {
108 std::ostringstream message; 115 std::ostringstream message;
109 message << "Not planar surface in untwis 116 message << "Not planar surface in untwisted Trapezoid: "
110 << GetName() << G4endl 117 << GetName() << G4endl
111 << "fDy2 is " << fDy2 << " but s 118 << "fDy2 is " << fDy2 << " but should be "
112 << pDytmp << "."; 119 << pDytmp << ".";
113 G4Exception("G4VTwistedFaceted::G4VTwist 120 G4Exception("G4VTwistedFaceted::G4VTwistedFaceted()", "GeomSolids0002",
114 FatalErrorInArgument, messag 121 FatalErrorInArgument, message);
115 } 122 }
116 } 123 }
117 124
118 #ifdef G4TWISTDEBUG 125 #ifdef G4TWISTDEBUG
119 if ( fDx1 == fDx2 && fDx3 == fDx4 ) 126 if ( fDx1 == fDx2 && fDx3 == fDx4 )
120 { 127 {
121 G4cout << "Trapezoid is a box" << G4endl 128 G4cout << "Trapezoid is a box" << G4endl ;
122 } 129 }
123 130
124 #endif 131 #endif
125 132
126 if ( ( fDx1 == fDx2 && fDx3 != fDx4 ) || ( 133 if ( ( fDx1 == fDx2 && fDx3 != fDx4 ) || ( fDx1 != fDx2 && fDx3 == fDx4 ) )
127 { 134 {
128 std::ostringstream message; 135 std::ostringstream message;
129 message << "Not planar surface in untwiste 136 message << "Not planar surface in untwisted Trapezoid: "
130 << GetName() << G4endl 137 << GetName() << G4endl
131 << "One endcap is rectangular, the 138 << "One endcap is rectangular, the other is a trapezoid." << G4endl
132 << "For planarity reasons they hav 139 << "For planarity reasons they have to be rectangles or trapezoids "
133 << "on both sides."; 140 << "on both sides.";
134 G4Exception("G4VTwistedFaceted::G4VTwisted 141 G4Exception("G4VTwistedFaceted::G4VTwistedFaceted()", "GeomSolids0002",
135 FatalErrorInArgument, message) 142 FatalErrorInArgument, message);
136 } 143 }
137 144
138 // twist angle 145 // twist angle
139 // 146 //
140 fPhiTwist = PhiTwist ; 147 fPhiTwist = PhiTwist ;
141 148
142 // tilt angle 149 // tilt angle
143 // 150 //
144 fAlph = pAlph ; 151 fAlph = pAlph ;
145 fTAlph = std::tan(fAlph) ; 152 fTAlph = std::tan(fAlph) ;
146 153
147 fTheta = pTheta ; 154 fTheta = pTheta ;
148 fPhi = pPhi ; 155 fPhi = pPhi ;
149 156
150 // dx in surface equation 157 // dx in surface equation
151 // 158 //
152 fdeltaX = 2 * fDz * std::tan(fTheta) * std:: 159 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
153 160
154 // dy in surface equation 161 // dy in surface equation
155 // 162 //
156 fdeltaY = 2 * fDz * std::tan(fTheta) * std:: 163 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
157 164
158 if ( ( fDx1 <= 2*kCarTolerance) << 165 if ( ! ( ( fDx1 > 2*kCarTolerance)
159 || ( fDx2 <= 2*kCarTolerance) << 166 && ( fDx2 > 2*kCarTolerance)
160 || ( fDx3 <= 2*kCarTolerance) << 167 && ( fDx3 > 2*kCarTolerance)
161 || ( fDx4 <= 2*kCarTolerance) << 168 && ( fDx4 > 2*kCarTolerance)
162 || ( fDy1 <= 2*kCarTolerance) << 169 && ( fDy1 > 2*kCarTolerance)
163 || ( fDy2 <= 2*kCarTolerance) << 170 && ( fDy2 > 2*kCarTolerance)
164 || ( fDz <= 2*kCarTolerance) << 171 && ( fDz > 2*kCarTolerance)
165 || ( std::fabs(fPhiTwist) <= 2*kAngTo << 172 && ( std::fabs(fPhiTwist) > 2*kAngTolerance )
166 || ( std::fabs(fPhiTwist) >= pi/2 ) << 173 && ( std::fabs(fPhiTwist) < pi/2 )
167 || ( std::fabs(fAlph) >= pi/2 ) << 174 && ( std::fabs(fAlph) < pi/2 )
168 || fTheta >= pi/2 || fTheta < 0 << 175 && ( fTheta < pi/2 && fTheta >= 0 ) )
169 ) 176 )
170 { 177 {
171 std::ostringstream message; 178 std::ostringstream message;
172 message << "Invalid dimensions. Too small, 179 message << "Invalid dimensions. Too small, or twist angle too big: "
173 << GetName() << G4endl 180 << GetName() << G4endl
174 << "fDx 1-4 = " << fDx1/cm << ", " 181 << "fDx 1-4 = " << fDx1/cm << ", " << fDx2/cm << ", "
175 << fDx3/cm << ", " << fDx4/cm << " 182 << fDx3/cm << ", " << fDx4/cm << " cm" << G4endl
176 << "fDy 1-2 = " << fDy1/cm << ", " 183 << "fDy 1-2 = " << fDy1/cm << ", " << fDy2/cm << ", "
177 << " cm" << G4endl 184 << " cm" << G4endl
178 << "fDz = " << fDz/cm << " cm" << 185 << "fDz = " << fDz/cm << " cm" << G4endl
179 << " twistangle " << fPhiTwist/deg 186 << " twistangle " << fPhiTwist/deg << " deg" << G4endl
180 << " phi,theta = " << fPhi/deg << 187 << " phi,theta = " << fPhi/deg << ", " << fTheta/deg << " deg";
181 G4Exception("G4TwistedTrap::G4VTwistedFace 188 G4Exception("G4TwistedTrap::G4VTwistedFaceted()",
182 "GeomSolids0002", FatalErrorIn 189 "GeomSolids0002", FatalErrorInArgument, message);
183 } 190 }
184 CreateSurfaces(); 191 CreateSurfaces();
>> 192 fCubicVolume = 2 * fDz * ( ( fDx1 + fDx2 ) * fDy1 + ( fDx3 + fDx4 ) * fDy2 );
185 } 193 }
186 194
187 195
188 //============================================ 196 //=====================================================================
189 //* Fake default constructor ----------------- 197 //* Fake default constructor ------------------------------------------
190 198
191 G4VTwistedFaceted::G4VTwistedFaceted( __void__ 199 G4VTwistedFaceted::G4VTwistedFaceted( __void__& a )
192 : G4VSolid(a), << 200 : G4VSolid(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
193 fLowerEndcap(nullptr), fUpperEndcap(nullpt << 201 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fDx(0.), fDy(0.), fAlph(0.),
194 fSide0(nullptr), fSide90(nullptr), fSide18 << 202 fTAlph(0.), fdeltaX(0.), fdeltaY(0.), fPhiTwist(0.),
>> 203 fLowerEndcap(0), fUpperEndcap(0), fSide0(0), fSide90(0), fSide180(0),
>> 204 fSide270(0), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
195 { 205 {
196 } 206 }
197 207
198 <<
199 //============================================ 208 //=====================================================================
200 //* destructor ------------------------------- 209 //* destructor --------------------------------------------------------
201 210
202 G4VTwistedFaceted::~G4VTwistedFaceted() 211 G4VTwistedFaceted::~G4VTwistedFaceted()
203 { 212 {
204 delete fLowerEndcap ; << 213 if (fLowerEndcap) delete fLowerEndcap ;
205 delete fUpperEndcap ; << 214 if (fUpperEndcap) delete fUpperEndcap ;
206 215
207 delete fSide0 ; << 216 if (fSide0) delete fSide0 ;
208 delete fSide90 ; << 217 if (fSide90) delete fSide90 ;
209 delete fSide180 ; << 218 if (fSide180) delete fSide180 ;
210 delete fSide270 ; << 219 if (fSide270) delete fSide270 ;
211 delete fpPolyhedron; fpPolyhedron = nullptr; << 220 if (fpPolyhedron) delete fpPolyhedron;
212 } 221 }
213 222
214 <<
215 //============================================ 223 //=====================================================================
216 //* Copy constructor ------------------------- 224 //* Copy constructor --------------------------------------------------
217 225
218 G4VTwistedFaceted::G4VTwistedFaceted(const G4V 226 G4VTwistedFaceted::G4VTwistedFaceted(const G4VTwistedFaceted& rhs)
219 : G4VSolid(rhs), << 227 : G4VSolid(rhs), fTheta(rhs.fTheta), fPhi(rhs.fPhi),
220 fCubicVolume(rhs.fCubicVolume), fSurfaceAr <<
221 fTheta(rhs.fTheta), fPhi(rhs.fPhi), <<
222 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.f 228 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fDy2(rhs.fDy2),
223 fDx3(rhs.fDx3), fDx4(rhs.fDx4), fDz(rhs.fD 229 fDx3(rhs.fDx3), fDx4(rhs.fDx4), fDz(rhs.fDz), fDx(rhs.fDx), fDy(rhs.fDy),
224 fAlph(rhs.fAlph), fTAlph(rhs.fTAlph), fdel 230 fAlph(rhs.fAlph), fTAlph(rhs.fTAlph), fdeltaX(rhs.fdeltaX),
225 fdeltaY(rhs.fdeltaY), fPhiTwist(rhs.fPhiTw << 231 fdeltaY(rhs.fdeltaY), fPhiTwist(rhs.fPhiTwist), fLowerEndcap(0),
226 fUpperEndcap(nullptr), fSide0(nullptr), fS << 232 fUpperEndcap(0), fSide0(0), fSide90(0), fSide180(0), fSide270(0),
>> 233 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
>> 234 fpPolyhedron(0),
>> 235 fLastInside(rhs.fLastInside), fLastNormal(rhs.fLastNormal),
>> 236 fLastDistanceToIn(rhs.fLastDistanceToIn),
>> 237 fLastDistanceToOut(rhs.fLastDistanceToOut),
>> 238 fLastDistanceToInWithV(rhs.fLastDistanceToInWithV),
>> 239 fLastDistanceToOutWithV(rhs.fLastDistanceToOutWithV)
227 { 240 {
228 CreateSurfaces(); 241 CreateSurfaces();
229 } 242 }
230 243
231 244
232 //============================================ 245 //=====================================================================
233 //* Assignment operator ---------------------- 246 //* Assignment operator -----------------------------------------------
234 247
235 G4VTwistedFaceted& G4VTwistedFaceted::operator 248 G4VTwistedFaceted& G4VTwistedFaceted::operator = (const G4VTwistedFaceted& rhs)
236 { 249 {
237 // Check assignment to self 250 // Check assignment to self
238 // 251 //
239 if (this == &rhs) { return *this; } 252 if (this == &rhs) { return *this; }
240 253
241 // Copy base class data 254 // Copy base class data
242 // 255 //
243 G4VSolid::operator=(rhs); 256 G4VSolid::operator=(rhs);
244 257
245 // Copy data 258 // Copy data
246 // 259 //
247 fTheta = rhs.fTheta; fPhi = rhs.fPhi; 260 fTheta = rhs.fTheta; fPhi = rhs.fPhi;
248 fDy1= rhs.fDy1; fDx1= rhs.fDx1; fDx2= rhs.f 261 fDy1= rhs.fDy1; fDx1= rhs.fDx1; fDx2= rhs.fDx2; fDy2= rhs.fDy2;
249 fDx3= rhs.fDx3; fDx4= rhs.fDx4; fDz= rhs.fD 262 fDx3= rhs.fDx3; fDx4= rhs.fDx4; fDz= rhs.fDz; fDx= rhs.fDx; fDy= rhs.fDy;
250 fAlph= rhs.fAlph; fTAlph= rhs.fTAlph; fdelt 263 fAlph= rhs.fAlph; fTAlph= rhs.fTAlph; fdeltaX= rhs.fdeltaX;
251 fdeltaY= rhs.fdeltaY; fPhiTwist= rhs.fPhiTw << 264 fdeltaY= rhs.fdeltaY; fPhiTwist= rhs.fPhiTwist; fLowerEndcap= 0;
252 fUpperEndcap= nullptr; fSide0= nullptr; fSi << 265 fUpperEndcap= 0; fSide0= 0; fSide90= 0; fSide180= 0; fSide270= 0;
253 fCubicVolume= rhs.fCubicVolume; fSurfaceAre 266 fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
254 fRebuildPolyhedron = false; << 267 fpPolyhedron= 0;
255 delete fpPolyhedron; fpPolyhedron = nullptr << 268 fLastInside= rhs.fLastInside; fLastNormal= rhs.fLastNormal;
>> 269 fLastDistanceToIn= rhs.fLastDistanceToIn;
>> 270 fLastDistanceToOut= rhs.fLastDistanceToOut;
>> 271 fLastDistanceToInWithV= rhs.fLastDistanceToInWithV;
>> 272 fLastDistanceToOutWithV= rhs.fLastDistanceToOutWithV;
256 273
257 CreateSurfaces(); 274 CreateSurfaces();
258 275
259 return *this; 276 return *this;
260 } 277 }
261 278
262 <<
263 //============================================ 279 //=====================================================================
264 //* ComputeDimensions ------------------------ 280 //* ComputeDimensions -------------------------------------------------
265 281
266 void G4VTwistedFaceted::ComputeDimensions(G4VP 282 void G4VTwistedFaceted::ComputeDimensions(G4VPVParameterisation* ,
267 cons 283 const G4int ,
268 cons 284 const G4VPhysicalVolume* )
269 { 285 {
270 G4Exception("G4VTwistedFaceted::ComputeDimen 286 G4Exception("G4VTwistedFaceted::ComputeDimensions()",
271 "GeomSolids0001", FatalException 287 "GeomSolids0001", FatalException,
272 "G4VTwistedFaceted does not supp 288 "G4VTwistedFaceted does not support Parameterisation.");
273 } 289 }
274 290
275 <<
276 //============================================ 291 //=====================================================================
277 //* Extent ----------------------------------- << 292 //* CalculateExtent ---------------------------------------------------
278 293
279 void G4VTwistedFaceted::BoundingLimits(G4Three << 294 G4bool
280 G4Three << 295 G4VTwistedFaceted::CalculateExtent( const EAxis pAxis,
>> 296 const G4VoxelLimits &pVoxelLimit,
>> 297 const G4AffineTransform &pTransform,
>> 298 G4double &pMin,
>> 299 G4double &pMax ) const
281 { 300 {
282 G4double cosPhi = std::cos(fPhi); << 301 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
283 G4double sinPhi = std::sin(fPhi); <<
284 G4double tanTheta = std::tan(fTheta); <<
285 G4double tanAlpha = fTAlph; <<
286 <<
287 G4double xmid1 = fDy1*tanAlpha; <<
288 G4double x1 = std::abs(xmid1 + fDx1); <<
289 G4double x2 = std::abs(xmid1 - fDx1); <<
290 G4double x3 = std::abs(xmid1 + fDx2); <<
291 G4double x4 = std::abs(xmid1 - fDx2); <<
292 G4double xmax1 = std::max(std::max(std::max( <<
293 G4double rmax1 = std::sqrt(xmax1*xmax1 + fDy <<
294 <<
295 G4double xmid2 = fDy2*tanAlpha; <<
296 G4double x5 = std::abs(xmid2 + fDx3); <<
297 G4double x6 = std::abs(xmid2 - fDx3); <<
298 G4double x7 = std::abs(xmid2 + fDx4); <<
299 G4double x8 = std::abs(xmid2 - fDx4); <<
300 G4double xmax2 = std::max(std::max(std::max( <<
301 G4double rmax2 = std::sqrt(xmax2*xmax2 + fDy <<
302 <<
303 G4double x0 = fDz*tanTheta*cosPhi; <<
304 G4double y0 = fDz*tanTheta*sinPhi; <<
305 G4double xmin = std::min(-x0 - rmax1, x0 - r <<
306 G4double ymin = std::min(-y0 - rmax1, y0 - r <<
307 G4double xmax = std::max(-x0 + rmax1, x0 + r <<
308 G4double ymax = std::max(-y0 + rmax1, y0 + r <<
309 pMin.set(xmin, ymin,-fDz); <<
310 pMax.set(xmax, ymax, fDz); <<
311 } <<
312 302
>> 303 if (!pTransform.IsRotated())
>> 304 {
>> 305 // Special case handling for unrotated boxes
>> 306 // Compute x/y/z mins and maxs respecting limits, with early returns
>> 307 // if outside limits. Then switch() on pAxis
>> 308
>> 309 G4double xoffset,xMin,xMax;
>> 310 G4double yoffset,yMin,yMax;
>> 311 G4double zoffset,zMin,zMax;
>> 312
>> 313 xoffset = pTransform.NetTranslation().x() ;
>> 314 xMin = xoffset - maxRad ;
>> 315 xMax = xoffset + maxRad ;
313 316
314 //============================================ << 317 if (pVoxelLimit.IsXLimited())
315 //* CalculateExtent -------------------------- << 318 {
>> 319 if ( xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance ||
>> 320 xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance ) return false;
>> 321 else
>> 322 {
>> 323 if (xMin < pVoxelLimit.GetMinXExtent())
>> 324 {
>> 325 xMin = pVoxelLimit.GetMinXExtent() ;
>> 326 }
>> 327 if (xMax > pVoxelLimit.GetMaxXExtent())
>> 328 {
>> 329 xMax = pVoxelLimit.GetMaxXExtent() ;
>> 330 }
>> 331 }
>> 332 }
>> 333 yoffset = pTransform.NetTranslation().y() ;
>> 334 yMin = yoffset - maxRad ;
>> 335 yMax = yoffset + maxRad ;
>> 336
>> 337 if (pVoxelLimit.IsYLimited())
>> 338 {
>> 339 if ( yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance ||
>> 340 yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance ) return false;
>> 341 else
>> 342 {
>> 343 if (yMin < pVoxelLimit.GetMinYExtent())
>> 344 {
>> 345 yMin = pVoxelLimit.GetMinYExtent() ;
>> 346 }
>> 347 if (yMax > pVoxelLimit.GetMaxYExtent())
>> 348 {
>> 349 yMax = pVoxelLimit.GetMaxYExtent() ;
>> 350 }
>> 351 }
>> 352 }
>> 353 zoffset = pTransform.NetTranslation().z() ;
>> 354 zMin = zoffset - fDz ;
>> 355 zMax = zoffset + fDz ;
>> 356
>> 357 if (pVoxelLimit.IsZLimited())
>> 358 {
>> 359 if ( zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance ||
>> 360 zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance ) return false;
>> 361 else
>> 362 {
>> 363 if (zMin < pVoxelLimit.GetMinZExtent())
>> 364 {
>> 365 zMin = pVoxelLimit.GetMinZExtent() ;
>> 366 }
>> 367 if (zMax > pVoxelLimit.GetMaxZExtent())
>> 368 {
>> 369 zMax = pVoxelLimit.GetMaxZExtent() ;
>> 370 }
>> 371 }
>> 372 }
>> 373 switch (pAxis)
>> 374 {
>> 375 case kXAxis:
>> 376 pMin = xMin ;
>> 377 pMax = xMax ;
>> 378 break ;
>> 379 case kYAxis:
>> 380 pMin=yMin;
>> 381 pMax=yMax;
>> 382 break;
>> 383 case kZAxis:
>> 384 pMin=zMin;
>> 385 pMax=zMax;
>> 386 break;
>> 387 default:
>> 388 break;
>> 389 }
>> 390 pMin -= kCarTolerance ;
>> 391 pMax += kCarTolerance ;
>> 392
>> 393 return true;
>> 394 }
>> 395 else // General rotated case - create and clip mesh to boundaries
>> 396 {
>> 397 G4bool existsAfterClip = false ;
>> 398 G4ThreeVectorList* vertices ;
316 399
317 G4bool << 400 pMin = +kInfinity ;
318 G4VTwistedFaceted::CalculateExtent( const EAxi << 401 pMax = -kInfinity ;
319 const G4Vo << 402
320 const G4Af << 403 // Calculate rotated vertex coordinates
321 G4do << 404
322 G4do << 405 vertices = CreateRotatedVertices(pTransform) ;
323 { << 406 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
324 G4ThreeVector bmin, bmax; << 407 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax) ;
325 << 408 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
326 // Get bounding box << 409
327 BoundingLimits(bmin,bmax); << 410 if (pVoxelLimit.IsLimited(pAxis) == false)
328 << 411 {
329 // Find extent << 412 if ( pMin != kInfinity || pMax != -kInfinity )
330 G4BoundingEnvelope bbox(bmin,bmax); << 413 {
331 return bbox.CalculateExtent(pAxis,pVoxelLimi << 414 existsAfterClip = true ;
>> 415
>> 416 // Add 2*tolerance to avoid precision troubles
>> 417
>> 418 pMin -= kCarTolerance;
>> 419 pMax += kCarTolerance;
>> 420 }
>> 421 }
>> 422 else
>> 423 {
>> 424 G4ThreeVector clipCentre(
>> 425 ( pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
>> 426 ( pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
>> 427 ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 );
>> 428
>> 429 if ( pMin != kInfinity || pMax != -kInfinity )
>> 430 {
>> 431 existsAfterClip = true ;
>> 432
>> 433 // Check to see if endpoints are in the solid
>> 434
>> 435 clipCentre(pAxis) = pVoxelLimit.GetMinExtent(pAxis);
>> 436
>> 437 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
>> 438 != kOutside)
>> 439 {
>> 440 pMin = pVoxelLimit.GetMinExtent(pAxis);
>> 441 }
>> 442 else
>> 443 {
>> 444 pMin -= kCarTolerance;
>> 445 }
>> 446 clipCentre(pAxis) = pVoxelLimit.GetMaxExtent(pAxis);
>> 447
>> 448 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
>> 449 != kOutside)
>> 450 {
>> 451 pMax = pVoxelLimit.GetMaxExtent(pAxis);
>> 452 }
>> 453 else
>> 454 {
>> 455 pMax += kCarTolerance;
>> 456 }
>> 457 }
>> 458
>> 459 // Check for case where completely enveloping clipping volume
>> 460 // If point inside then we are confident that the solid completely
>> 461 // envelopes the clipping volume. Hence set min/max extents according
>> 462 // to clipping volume extents along the specified axis.
>> 463
>> 464 else if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
>> 465 != kOutside)
>> 466 {
>> 467 existsAfterClip = true ;
>> 468 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
>> 469 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
>> 470 }
>> 471 }
>> 472 delete vertices;
>> 473 return existsAfterClip;
>> 474 }
>> 475
>> 476
332 } 477 }
333 478
>> 479 G4ThreeVectorList* G4VTwistedFaceted::
>> 480 CreateRotatedVertices(const G4AffineTransform& pTransform) const
>> 481 {
>> 482
>> 483 G4ThreeVectorList* vertices = new G4ThreeVectorList();
>> 484
>> 485 if (vertices)
>> 486 {
>> 487 vertices->reserve(8);
>> 488
>> 489 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
>> 490
>> 491 G4ThreeVector vertex0(-maxRad,-maxRad,-fDz) ;
>> 492 G4ThreeVector vertex1(maxRad,-maxRad,-fDz) ;
>> 493 G4ThreeVector vertex2(maxRad,maxRad,-fDz) ;
>> 494 G4ThreeVector vertex3(-maxRad,maxRad,-fDz) ;
>> 495 G4ThreeVector vertex4(-maxRad,-maxRad,fDz) ;
>> 496 G4ThreeVector vertex5(maxRad,-maxRad,fDz) ;
>> 497 G4ThreeVector vertex6(maxRad,maxRad,fDz) ;
>> 498 G4ThreeVector vertex7(-maxRad,maxRad,fDz) ;
>> 499
>> 500 vertices->push_back(pTransform.TransformPoint(vertex0));
>> 501 vertices->push_back(pTransform.TransformPoint(vertex1));
>> 502 vertices->push_back(pTransform.TransformPoint(vertex2));
>> 503 vertices->push_back(pTransform.TransformPoint(vertex3));
>> 504 vertices->push_back(pTransform.TransformPoint(vertex4));
>> 505 vertices->push_back(pTransform.TransformPoint(vertex5));
>> 506 vertices->push_back(pTransform.TransformPoint(vertex6));
>> 507 vertices->push_back(pTransform.TransformPoint(vertex7));
>> 508 }
>> 509 else
>> 510 {
>> 511 DumpInfo();
>> 512 G4Exception("G4VTwistedFaceted::CreateRotatedVertices()",
>> 513 "GeomSolids0003", FatalException,
>> 514 "Error in allocation of vertices. Out of memory !");
>> 515 }
>> 516 return vertices;
>> 517 }
334 518
335 //============================================ 519 //=====================================================================
336 //* Inside ----------------------------------- 520 //* Inside ------------------------------------------------------------
337 521
338 EInside G4VTwistedFaceted::Inside(const G4Thre 522 EInside G4VTwistedFaceted::Inside(const G4ThreeVector& p) const
339 { 523 {
340 524
341 EInside tmpin = kOutside ; << 525 G4ThreeVector *tmpp;
>> 526 EInside *tmpin;
>> 527 if (fLastInside.p == p) {
>> 528 return fLastInside.inside;
>> 529 } else {
>> 530 tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
>> 531 tmpin = const_cast<EInside*>(&(fLastInside.inside));
>> 532 tmpp->set(p.x(), p.y(), p.z());
>> 533 }
>> 534
>> 535 *tmpin = kOutside ;
342 536
343 G4double phi = p.z()/(2*fDz) * fPhiTwist ; 537 G4double phi = p.z()/(2*fDz) * fPhiTwist ; // rotate the point to z=0
344 G4double cphi = std::cos(-phi) ; 538 G4double cphi = std::cos(-phi) ;
345 G4double sphi = std::sin(-phi) ; 539 G4double sphi = std::sin(-phi) ;
346 540
347 G4double px = p.x() + fdeltaX * ( -phi/fPh 541 G4double px = p.x() + fdeltaX * ( -phi/fPhiTwist) ; // shift
348 G4double py = p.y() + fdeltaY * ( -phi/fPh 542 G4double py = p.y() + fdeltaY * ( -phi/fPhiTwist) ;
349 G4double pz = p.z() ; 543 G4double pz = p.z() ;
350 544
351 G4double posx = px * cphi - py * sphi ; 545 G4double posx = px * cphi - py * sphi ; // rotation
352 G4double posy = px * sphi + py * cphi ; 546 G4double posy = px * sphi + py * cphi ;
353 G4double posz = pz ; 547 G4double posz = pz ;
354 548
355 G4double xMin = Xcoef(posy,phi,fTAlph) - 2* 549 G4double xMin = Xcoef(posy,phi,fTAlph) - 2*Xcoef(posy,phi,0.) ;
356 G4double xMax = Xcoef(posy,phi,fTAlph) ; 550 G4double xMax = Xcoef(posy,phi,fTAlph) ;
357 551
358 G4double yMax = GetValueB(phi)/2. ; // b(p 552 G4double yMax = GetValueB(phi)/2. ; // b(phi)/2 is limit
359 G4double yMin = -yMax ; 553 G4double yMin = -yMax ;
360 554
361 #ifdef G4TWISTDEBUG 555 #ifdef G4TWISTDEBUG
362 556
363 G4cout << "inside called: p = " << p << G4e 557 G4cout << "inside called: p = " << p << G4endl ;
364 G4cout << "fDx1 = " << fDx1 << G4endl ; 558 G4cout << "fDx1 = " << fDx1 << G4endl ;
365 G4cout << "fDx2 = " << fDx2 << G4endl ; 559 G4cout << "fDx2 = " << fDx2 << G4endl ;
366 G4cout << "fDx3 = " << fDx3 << G4endl ; 560 G4cout << "fDx3 = " << fDx3 << G4endl ;
367 G4cout << "fDx4 = " << fDx4 << G4endl ; 561 G4cout << "fDx4 = " << fDx4 << G4endl ;
368 562
369 G4cout << "fDy1 = " << fDy1 << G4endl ; 563 G4cout << "fDy1 = " << fDy1 << G4endl ;
370 G4cout << "fDy2 = " << fDy2 << G4endl ; 564 G4cout << "fDy2 = " << fDy2 << G4endl ;
371 565
372 G4cout << "fDz = " << fDz << G4endl ; 566 G4cout << "fDz = " << fDz << G4endl ;
373 567
374 G4cout << "Tilt angle alpha = " << fAlph << 568 G4cout << "Tilt angle alpha = " << fAlph << G4endl ;
375 G4cout << "phi,theta = " << fPhi << " , " < 569 G4cout << "phi,theta = " << fPhi << " , " << fTheta << G4endl ;
376 570
377 G4cout << "Twist angle = " << fPhiTwist << 571 G4cout << "Twist angle = " << fPhiTwist << G4endl ;
378 572
379 G4cout << "posx = " << posx << G4endl ; 573 G4cout << "posx = " << posx << G4endl ;
380 G4cout << "posy = " << posy << G4endl ; 574 G4cout << "posy = " << posy << G4endl ;
381 G4cout << "xMin = " << xMin << G4endl ; 575 G4cout << "xMin = " << xMin << G4endl ;
382 G4cout << "xMax = " << xMax << G4endl ; 576 G4cout << "xMax = " << xMax << G4endl ;
383 G4cout << "yMin = " << yMin << G4endl ; 577 G4cout << "yMin = " << yMin << G4endl ;
384 G4cout << "yMax = " << yMax << G4endl ; 578 G4cout << "yMax = " << yMax << G4endl ;
385 579
386 #endif 580 #endif
387 581
388 582
389 if ( posx <= xMax - kCarTolerance*0.5 583 if ( posx <= xMax - kCarTolerance*0.5
390 && posx >= xMin + kCarTolerance*0.5 ) 584 && posx >= xMin + kCarTolerance*0.5 )
391 { 585 {
392 if ( posy <= yMax - kCarTolerance*0.5 586 if ( posy <= yMax - kCarTolerance*0.5
393 && posy >= yMin + kCarTolerance*0.5 ) 587 && posy >= yMin + kCarTolerance*0.5 )
394 { 588 {
395 if (std::fabs(posz) <= fDz - kCarTo << 589 if (std::fabs(posz) <= fDz - kCarTolerance*0.5 ) *tmpin = kInside ;
396 else if (std::fabs(posz) <= fDz + kCarTo << 590 else if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) *tmpin = kSurface ;
397 } 591 }
398 else if ( posy <= yMax + kCarTolerance*0.5 592 else if ( posy <= yMax + kCarTolerance*0.5
399 && posy >= yMin - kCarTolerance*0.5 593 && posy >= yMin - kCarTolerance*0.5 )
400 { 594 {
401 if (std::fabs(posz) <= fDz + kCarToleran << 595 if (std::fabs(posz) <= fDz + kCarTolerance*0.5 ) *tmpin = kSurface ;
402 } 596 }
403 } 597 }
404 else if ( posx <= xMax + kCarTolerance*0.5 598 else if ( posx <= xMax + kCarTolerance*0.5
405 && posx >= xMin - kCarTolerance*0.5 ) 599 && posx >= xMin - kCarTolerance*0.5 )
406 { 600 {
407 if ( posy <= yMax + kCarTolerance*0.5 601 if ( posy <= yMax + kCarTolerance*0.5
408 && posy >= yMin - kCarTolerance*0.5 ) 602 && posy >= yMin - kCarTolerance*0.5 )
409 { 603 {
410 if (std::fabs(posz) <= fDz + kCarToleran << 604 if (std::fabs(posz) <= fDz + kCarTolerance*0.5) *tmpin = kSurface ;
411 } 605 }
412 } 606 }
413 607
414 #ifdef G4TWISTDEBUG 608 #ifdef G4TWISTDEBUG
415 G4cout << "inside = " << tmpin << G4endl ; << 609 G4cout << "inside = " << fLastInside.inside << G4endl ;
416 #endif 610 #endif
417 611
418 return tmpin; << 612 return fLastInside.inside;
419 613
420 } 614 }
421 615
422 <<
423 //============================================ 616 //=====================================================================
424 //* SurfaceNormal ---------------------------- 617 //* SurfaceNormal -----------------------------------------------------
425 618
426 G4ThreeVector G4VTwistedFaceted::SurfaceNormal 619 G4ThreeVector G4VTwistedFaceted::SurfaceNormal(const G4ThreeVector& p) const
427 { 620 {
428 // 621 //
429 // return the normal unit vector to the Hyp 622 // return the normal unit vector to the Hyperbolical Surface at a point
430 // p on (or nearly on) the surface 623 // p on (or nearly on) the surface
431 // 624 //
432 // Which of the three or four surfaces are 625 // Which of the three or four surfaces are we closest to?
433 // 626 //
434 627
>> 628 if (fLastNormal.p == p)
>> 629 {
>> 630 return fLastNormal.vec;
>> 631 }
>> 632
>> 633 G4ThreeVector *tmpp = const_cast<G4ThreeVector*>(&(fLastNormal.p));
>> 634 G4ThreeVector *tmpnormal = const_cast<G4ThreeVector*>(&(fLastNormal.vec));
>> 635 G4VTwistSurface **tmpsurface = const_cast<G4VTwistSurface**>(fLastNormal.surface);
>> 636 tmpp->set(p.x(), p.y(), p.z());
435 637
436 G4double distance = kInfinity; << 638 G4double distance = kInfinity;
437 639
438 G4VTwistSurface* surfaces[6]; << 640 G4VTwistSurface *surfaces[6];
439 641
440 surfaces[0] = fSide0 ; 642 surfaces[0] = fSide0 ;
441 surfaces[1] = fSide90 ; 643 surfaces[1] = fSide90 ;
442 surfaces[2] = fSide180 ; 644 surfaces[2] = fSide180 ;
443 surfaces[3] = fSide270 ; 645 surfaces[3] = fSide270 ;
444 surfaces[4] = fLowerEndcap; 646 surfaces[4] = fLowerEndcap;
445 surfaces[5] = fUpperEndcap; 647 surfaces[5] = fUpperEndcap;
446 648
447 G4ThreeVector xx; 649 G4ThreeVector xx;
448 G4ThreeVector bestxx; 650 G4ThreeVector bestxx;
449 G4int i; 651 G4int i;
450 G4int besti = -1; 652 G4int besti = -1;
451 for (i=0; i< 6; i++) 653 for (i=0; i< 6; i++)
452 { 654 {
453 G4double tmpdistance = surfaces[i]->Dist 655 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
454 if (tmpdistance < distance) 656 if (tmpdistance < distance)
455 { 657 {
456 distance = tmpdistance; 658 distance = tmpdistance;
457 bestxx = xx; 659 bestxx = xx;
458 besti = i; 660 besti = i;
459 } 661 }
460 } 662 }
461 663
462 return surfaces[besti]->GetNormal(bestxx, t << 664 tmpsurface[0] = surfaces[besti];
>> 665 *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
>> 666
>> 667 return fLastNormal.vec;
463 } 668 }
464 669
465 <<
466 //============================================ 670 //=====================================================================
467 //* DistanceToIn (p, v) ---------------------- 671 //* DistanceToIn (p, v) -----------------------------------------------
468 672
469 G4double G4VTwistedFaceted::DistanceToIn (cons 673 G4double G4VTwistedFaceted::DistanceToIn (const G4ThreeVector& p,
470 cons 674 const G4ThreeVector& v ) const
471 { 675 {
472 676
473 // DistanceToIn (p, v): 677 // DistanceToIn (p, v):
474 // Calculate distance to surface of shape f 678 // Calculate distance to surface of shape from `outside'
475 // along with the v, allowing for tolerance 679 // along with the v, allowing for tolerance.
476 // The function returns kInfinity if no int 680 // The function returns kInfinity if no intersection or
477 // just grazing within tolerance. 681 // just grazing within tolerance.
478 682
479 // 683 //
>> 684 // checking last value
>> 685 //
>> 686
>> 687 G4ThreeVector *tmpp;
>> 688 G4ThreeVector *tmpv;
>> 689 G4double *tmpdist;
>> 690 if (fLastDistanceToInWithV.p == p && fLastDistanceToInWithV.vec == v)
>> 691 {
>> 692 return fLastDistanceToIn.value;
>> 693 }
>> 694 else
>> 695 {
>> 696 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
>> 697 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
>> 698 tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
>> 699 tmpp->set(p.x(), p.y(), p.z());
>> 700 tmpv->set(v.x(), v.y(), v.z());
>> 701 }
>> 702
>> 703 //
480 // Calculate DistanceToIn(p,v) 704 // Calculate DistanceToIn(p,v)
481 // 705 //
482 706
483 EInside currentside = Inside(p); 707 EInside currentside = Inside(p);
484 708
485 if (currentside == kInside) 709 if (currentside == kInside)
486 { 710 {
487 } 711 }
488 else if (currentside == kSurface) 712 else if (currentside == kSurface)
489 { 713 {
490 // particle is just on a boundary. 714 // particle is just on a boundary.
491 // if the particle is entering to the vol 715 // if the particle is entering to the volume, return 0
492 // 716 //
493 G4ThreeVector normal = SurfaceNormal(p); 717 G4ThreeVector normal = SurfaceNormal(p);
494 if (normal*v < 0) 718 if (normal*v < 0)
495 { 719 {
496 return 0; << 720 *tmpdist = 0;
>> 721 return fLastDistanceToInWithV.value;
497 } 722 }
498 } 723 }
499 724
500 // now, we can take smallest positive dista 725 // now, we can take smallest positive distance.
501 726
502 // Initialize 727 // Initialize
503 // 728 //
504 G4double distance = kInfinity; << 729 G4double distance = kInfinity;
505 730
506 // Find intersections and choose nearest on 731 // Find intersections and choose nearest one
507 // 732 //
508 G4VTwistSurface *surfaces[6]; 733 G4VTwistSurface *surfaces[6];
509 734
510 surfaces[0] = fSide0; 735 surfaces[0] = fSide0;
511 surfaces[1] = fSide90 ; 736 surfaces[1] = fSide90 ;
512 surfaces[2] = fSide180 ; 737 surfaces[2] = fSide180 ;
513 surfaces[3] = fSide270 ; 738 surfaces[3] = fSide270 ;
514 surfaces[4] = fLowerEndcap; 739 surfaces[4] = fLowerEndcap;
515 surfaces[5] = fUpperEndcap; 740 surfaces[5] = fUpperEndcap;
516 741
517 G4ThreeVector xx; 742 G4ThreeVector xx;
518 G4ThreeVector bestxx; 743 G4ThreeVector bestxx;
519 for (const auto & surface : surfaces) << 744 G4int i;
>> 745 for (i=0; i < 6 ; i++)
520 { 746 {
521 #ifdef G4TWISTDEBUG 747 #ifdef G4TWISTDEBUG
522 G4cout << G4endl << "surface " << &surfa << 748 G4cout << G4endl << "surface " << i << ": " << G4endl << G4endl ;
523 #endif 749 #endif
524 G4double tmpdistance = surface->Distance << 750 G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
525 #ifdef G4TWISTDEBUG 751 #ifdef G4TWISTDEBUG
526 G4cout << "Solid DistanceToIn : distance 752 G4cout << "Solid DistanceToIn : distance = " << tmpdistance << G4endl ;
527 G4cout << "intersection point = " << xx 753 G4cout << "intersection point = " << xx << G4endl ;
528 #endif 754 #endif
529 if (tmpdistance < distance) 755 if (tmpdistance < distance)
530 { 756 {
531 distance = tmpdistance; 757 distance = tmpdistance;
532 bestxx = xx; 758 bestxx = xx;
533 } 759 }
534 } 760 }
535 761
536 #ifdef G4TWISTDEBUG 762 #ifdef G4TWISTDEBUG
537 G4cout << "best distance = " << distance << 763 G4cout << "best distance = " << distance << G4endl ;
538 #endif 764 #endif
539 765
>> 766 *tmpdist = distance;
540 // timer.Stop(); 767 // timer.Stop();
541 return distance; << 768 return fLastDistanceToInWithV.value;
542 } 769 }
543 770
544 771
545 //============================================ <<
546 //* DistanceToIn (p) ------------------------- <<
547 <<
548 G4double G4VTwistedFaceted::DistanceToIn (cons 772 G4double G4VTwistedFaceted::DistanceToIn (const G4ThreeVector& p) const
549 { 773 {
550 // DistanceToIn(p): 774 // DistanceToIn(p):
551 // Calculate distance to surface of shape f 775 // Calculate distance to surface of shape from `outside',
552 // allowing for tolerance 776 // allowing for tolerance
553 // 777 //
554 778
555 // 779 //
>> 780 // checking last value
>> 781 //
>> 782
>> 783 G4ThreeVector *tmpp;
>> 784 G4double *tmpdist;
>> 785 if (fLastDistanceToIn.p == p)
>> 786 {
>> 787 return fLastDistanceToIn.value;
>> 788 }
>> 789 else
>> 790 {
>> 791 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
>> 792 tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
>> 793 tmpp->set(p.x(), p.y(), p.z());
>> 794 }
>> 795
>> 796 //
556 // Calculate DistanceToIn(p) 797 // Calculate DistanceToIn(p)
557 // 798 //
558 799
559 EInside currentside = Inside(p); 800 EInside currentside = Inside(p);
560 801
561 switch (currentside) 802 switch (currentside)
562 { 803 {
563 case (kInside) : 804 case (kInside) :
564 { 805 {
565 } 806 }
566 807
567 case (kSurface) : 808 case (kSurface) :
568 { 809 {
569 return 0; << 810 *tmpdist = 0.;
>> 811 return fLastDistanceToIn.value;
570 } 812 }
571 813
572 case (kOutside) : 814 case (kOutside) :
573 { 815 {
574 // Initialize 816 // Initialize
575 // 817 //
576 G4double distance = kInfinity; << 818 G4double distance = kInfinity;
577 819
578 // Find intersections and choose near 820 // Find intersections and choose nearest one
579 // 821 //
580 G4VTwistSurface* surfaces[6]; << 822 G4VTwistSurface *surfaces[6];
581 823
582 surfaces[0] = fSide0; 824 surfaces[0] = fSide0;
583 surfaces[1] = fSide90 ; 825 surfaces[1] = fSide90 ;
584 surfaces[2] = fSide180 ; 826 surfaces[2] = fSide180 ;
585 surfaces[3] = fSide270 ; 827 surfaces[3] = fSide270 ;
586 surfaces[4] = fLowerEndcap; 828 surfaces[4] = fLowerEndcap;
587 surfaces[5] = fUpperEndcap; 829 surfaces[5] = fUpperEndcap;
588 830
>> 831 G4int i;
589 G4ThreeVector xx; 832 G4ThreeVector xx;
590 G4ThreeVector bestxx; 833 G4ThreeVector bestxx;
591 for (const auto & surface : surfaces) << 834 for (i=0; i< 6; i++)
592 { 835 {
593 G4double tmpdistance = surface->Di << 836 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
594 if (tmpdistance < distance) 837 if (tmpdistance < distance)
595 { 838 {
596 distance = tmpdistance; 839 distance = tmpdistance;
597 bestxx = xx; 840 bestxx = xx;
598 } 841 }
599 } 842 }
600 return distance; << 843 *tmpdist = distance;
>> 844 return fLastDistanceToIn.value;
601 } 845 }
602 846
603 default: << 847 default :
604 { 848 {
605 G4Exception("G4VTwistedFaceted::Dista 849 G4Exception("G4VTwistedFaceted::DistanceToIn(p)", "GeomSolids0003",
606 FatalException, "Unknown 850 FatalException, "Unknown point location!");
607 } 851 }
608 } // switch end 852 } // switch end
609 853
610 return 0.; << 854 return 0;
611 } 855 }
612 856
613 857
614 //============================================ 858 //=====================================================================
615 //* DistanceToOut (p, v) --------------------- 859 //* DistanceToOut (p, v) ----------------------------------------------
616 860
617 G4double 861 G4double
618 G4VTwistedFaceted::DistanceToOut( const G4Thre 862 G4VTwistedFaceted::DistanceToOut( const G4ThreeVector& p,
619 const G4Thre 863 const G4ThreeVector& v,
620 const G4bool 864 const G4bool calcNorm,
621 G4bool << 865 G4bool *validNorm,
622 G4Thre << 866 G4ThreeVector *norm ) const
623 { 867 {
624 // DistanceToOut (p, v): 868 // DistanceToOut (p, v):
625 // Calculate distance to surface of shape f 869 // Calculate distance to surface of shape from `inside'
626 // along with the v, allowing for tolerance 870 // along with the v, allowing for tolerance.
627 // The function returns kInfinity if no int 871 // The function returns kInfinity if no intersection or
628 // just grazing within tolerance. 872 // just grazing within tolerance.
629 873
630 // 874 //
>> 875 // checking last value
>> 876 //
>> 877
>> 878 G4ThreeVector *tmpp;
>> 879 G4ThreeVector *tmpv;
>> 880 G4double *tmpdist;
>> 881 if (fLastDistanceToOutWithV.p == p && fLastDistanceToOutWithV.vec == v )
>> 882 {
>> 883 return fLastDistanceToOutWithV.value;
>> 884 }
>> 885 else
>> 886 {
>> 887 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
>> 888 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
>> 889 tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
>> 890 tmpp->set(p.x(), p.y(), p.z());
>> 891 tmpv->set(v.x(), v.y(), v.z());
>> 892 }
>> 893
>> 894 //
631 // Calculate DistanceToOut(p,v) 895 // Calculate DistanceToOut(p,v)
632 // 896 //
633 897
634 EInside currentside = Inside(p); 898 EInside currentside = Inside(p);
635 899
636 if (currentside == kOutside) 900 if (currentside == kOutside)
637 { 901 {
638 } 902 }
639 else if (currentside == kSurface) 903 else if (currentside == kSurface)
640 { 904 {
641 // particle is just on a boundary. 905 // particle is just on a boundary.
642 // if the particle is exiting from the v 906 // if the particle is exiting from the volume, return 0
643 // 907 //
644 G4ThreeVector normal = SurfaceNormal(p); 908 G4ThreeVector normal = SurfaceNormal(p);
>> 909 G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
645 if (normal*v > 0) 910 if (normal*v > 0)
646 { 911 {
647 if (calcNorm) 912 if (calcNorm)
648 { 913 {
649 *norm = normal; << 914 *norm = (blockedsurface->GetNormal(p, true));
650 *validNorm = true; << 915 *validNorm = blockedsurface->IsValidNorm();
651 } 916 }
>> 917 *tmpdist = 0.;
652 // timer.Stop(); 918 // timer.Stop();
653 return 0; << 919 return fLastDistanceToOutWithV.value;
654 } 920 }
655 } 921 }
656 922
657 // now, we can take smallest positive dista 923 // now, we can take smallest positive distance.
658 924
659 // Initialize 925 // Initialize
660 G4double distance = kInfinity; << 926 G4double distance = kInfinity;
661 927
662 // find intersections and choose nearest on 928 // find intersections and choose nearest one.
663 G4VTwistSurface *surfaces[6]; 929 G4VTwistSurface *surfaces[6];
664 930
665 surfaces[0] = fSide0; 931 surfaces[0] = fSide0;
666 surfaces[1] = fSide90 ; 932 surfaces[1] = fSide90 ;
667 surfaces[2] = fSide180 ; 933 surfaces[2] = fSide180 ;
668 surfaces[3] = fSide270 ; 934 surfaces[3] = fSide270 ;
669 surfaces[4] = fLowerEndcap; 935 surfaces[4] = fLowerEndcap;
670 surfaces[5] = fUpperEndcap; 936 surfaces[5] = fUpperEndcap;
671 937
>> 938 G4int i;
672 G4int besti = -1; 939 G4int besti = -1;
673 G4ThreeVector xx; 940 G4ThreeVector xx;
674 G4ThreeVector bestxx; 941 G4ThreeVector bestxx;
675 for (auto i=0; i<6 ; ++i) << 942 for (i=0; i< 6 ; i++) {
676 { <<
677 G4double tmpdistance = surfaces[i]->Dist 943 G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
678 if (tmpdistance < distance) 944 if (tmpdistance < distance)
679 { 945 {
680 distance = tmpdistance; 946 distance = tmpdistance;
681 bestxx = xx; 947 bestxx = xx;
682 besti = i; 948 besti = i;
683 } 949 }
684 } 950 }
685 951
686 if (calcNorm) 952 if (calcNorm)
687 { 953 {
688 if (besti != -1) 954 if (besti != -1)
689 { 955 {
690 *norm = (surfaces[besti]->GetNormal(p 956 *norm = (surfaces[besti]->GetNormal(p, true));
691 *validNorm = surfaces[besti]->IsValid 957 *validNorm = surfaces[besti]->IsValidNorm();
692 } 958 }
693 } 959 }
694 960
695 return distance; << 961 *tmpdist = distance;
>> 962 // timer.Stop();
>> 963 return fLastDistanceToOutWithV.value;
696 } 964 }
697 965
698 966
699 //============================================ <<
700 //* DistanceToOut (p) ------------------------ <<
701 <<
702 G4double G4VTwistedFaceted::DistanceToOut( con 967 G4double G4VTwistedFaceted::DistanceToOut( const G4ThreeVector& p ) const
703 { 968 {
704 // DistanceToOut(p): 969 // DistanceToOut(p):
705 // Calculate distance to surface of shape f 970 // Calculate distance to surface of shape from `inside',
706 // allowing for tolerance 971 // allowing for tolerance
>> 972
>> 973 //
>> 974 // checking last value
>> 975 //
>> 976
>> 977 G4ThreeVector *tmpp;
>> 978 G4double *tmpdist;
>> 979
>> 980 if (fLastDistanceToOut.p == p)
>> 981 {
>> 982 return fLastDistanceToOut.value;
>> 983 }
>> 984 else
>> 985 {
>> 986 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
>> 987 tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
>> 988 tmpp->set(p.x(), p.y(), p.z());
>> 989 }
707 990
708 // 991 //
709 // Calculate DistanceToOut(p) 992 // Calculate DistanceToOut(p)
710 // 993 //
711 994
712 EInside currentside = Inside(p); 995 EInside currentside = Inside(p);
713 G4double retval = kInfinity; 996 G4double retval = kInfinity;
714 997
715 switch (currentside) 998 switch (currentside)
716 { 999 {
717 case (kOutside) : 1000 case (kOutside) :
718 { 1001 {
719 #ifdef G4SPECSDEBUG 1002 #ifdef G4SPECSDEBUG
720 G4int oldprc = G4cout.precision(16) ; 1003 G4int oldprc = G4cout.precision(16) ;
721 G4cout << G4endl ; 1004 G4cout << G4endl ;
722 DumpInfo(); 1005 DumpInfo();
723 G4cout << "Position:" << G4endl << G4 1006 G4cout << "Position:" << G4endl << G4endl ;
724 G4cout << "p.x() = " << p.x()/mm << 1007 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
725 G4cout << "p.y() = " << p.y()/mm << 1008 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
726 G4cout << "p.z() = " << p.z()/mm << 1009 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
727 G4cout.precision(oldprc) ; 1010 G4cout.precision(oldprc) ;
728 G4Exception("G4VTwistedFaceted::Distan 1011 G4Exception("G4VTwistedFaceted::DistanceToOut(p)", "GeomSolids1002",
729 JustWarning, "Point p is o 1012 JustWarning, "Point p is outside !?" );
730 #endif 1013 #endif
731 break; 1014 break;
732 } 1015 }
733 case (kSurface) : 1016 case (kSurface) :
734 { 1017 {
735 retval = 0; << 1018 *tmpdist = 0.;
>> 1019 retval = fLastDistanceToOut.value;
736 break; 1020 break;
737 } 1021 }
738 1022
739 case (kInside) : 1023 case (kInside) :
740 { 1024 {
741 // Initialize 1025 // Initialize
742 // 1026 //
743 G4double distance = kInfinity; << 1027 G4double distance = kInfinity;
744 1028
745 // find intersections and choose neare 1029 // find intersections and choose nearest one
746 // 1030 //
747 G4VTwistSurface* surfaces[6]; << 1031 G4VTwistSurface *surfaces[6];
748 1032
749 surfaces[0] = fSide0; 1033 surfaces[0] = fSide0;
750 surfaces[1] = fSide90 ; 1034 surfaces[1] = fSide90 ;
751 surfaces[2] = fSide180 ; 1035 surfaces[2] = fSide180 ;
752 surfaces[3] = fSide270 ; 1036 surfaces[3] = fSide270 ;
753 surfaces[4] = fLowerEndcap; 1037 surfaces[4] = fLowerEndcap;
754 surfaces[5] = fUpperEndcap; 1038 surfaces[5] = fUpperEndcap;
755 1039
>> 1040 G4int i;
756 G4ThreeVector xx; 1041 G4ThreeVector xx;
757 G4ThreeVector bestxx; 1042 G4ThreeVector bestxx;
758 for (const auto & surface : surfaces) << 1043 for (i=0; i< 6; i++)
759 { 1044 {
760 G4double tmpdistance = surface->Dist << 1045 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
761 if (tmpdistance < distance) 1046 if (tmpdistance < distance)
762 { 1047 {
763 distance = tmpdistance; 1048 distance = tmpdistance;
764 bestxx = xx; 1049 bestxx = xx;
765 } 1050 }
766 } 1051 }
767 retval = distance; << 1052 *tmpdist = distance;
>> 1053
>> 1054 retval = fLastDistanceToOut.value;
768 break; 1055 break;
769 } 1056 }
770 1057
771 default : 1058 default :
772 { 1059 {
773 G4Exception("G4VTwistedFaceted::Distan 1060 G4Exception("G4VTwistedFaceted::DistanceToOut(p)", "GeomSolids0003",
774 FatalException, "Unknown p 1061 FatalException, "Unknown point location!");
775 break; 1062 break;
776 } 1063 }
777 } // switch end 1064 } // switch end
778 1065
779 return retval; 1066 return retval;
780 } 1067 }
781 1068
782 1069
783 //============================================ 1070 //=====================================================================
784 //* StreamInfo ------------------------------- 1071 //* StreamInfo --------------------------------------------------------
785 1072
786 std::ostream& G4VTwistedFaceted::StreamInfo(st 1073 std::ostream& G4VTwistedFaceted::StreamInfo(std::ostream& os) const
787 { 1074 {
788 // 1075 //
789 // Stream object contents to an output strea 1076 // Stream object contents to an output stream
790 // 1077 //
791 G4long oldprc = os.precision(16); << 1078 G4int oldprc = os.precision(16);
792 os << "------------------------------------- 1079 os << "-----------------------------------------------------------\n"
793 << " *** Dump for solid - " << GetName 1080 << " *** Dump for solid - " << GetName() << " ***\n"
794 << " ================================= 1081 << " ===================================================\n"
795 << " Solid type: G4VTwistedFaceted\n" 1082 << " Solid type: G4VTwistedFaceted\n"
796 << " Parameters: \n" 1083 << " Parameters: \n"
797 << " polar angle theta = " << fTheta/ 1084 << " polar angle theta = " << fTheta/degree << " deg" << G4endl
798 << " azimuthal angle phi = " << fPhi/de 1085 << " azimuthal angle phi = " << fPhi/degree << " deg" << G4endl
799 << " tilt angle alpha = " << fAlph/de 1086 << " tilt angle alpha = " << fAlph/degree << " deg" << G4endl
800 << " TWIST angle = " << fPhiTwis 1087 << " TWIST angle = " << fPhiTwist/degree << " deg" << G4endl
801 << " Half length along y (lower endcap) 1088 << " Half length along y (lower endcap) = " << fDy1/cm << " cm"
802 << G4endl 1089 << G4endl
803 << " Half length along x (lower endcap, 1090 << " Half length along x (lower endcap, bottom) = " << fDx1/cm << " cm"
804 << G4endl 1091 << G4endl
805 << " Half length along x (lower endcap, 1092 << " Half length along x (lower endcap, top) = " << fDx2/cm << " cm"
806 << G4endl 1093 << G4endl
807 << " Half length along y (upper endcap) 1094 << " Half length along y (upper endcap) = " << fDy2/cm << " cm"
808 << G4endl 1095 << G4endl
809 << " Half length along x (upper endcap, 1096 << " Half length along x (upper endcap, bottom) = " << fDx3/cm << " cm"
810 << G4endl 1097 << G4endl
811 << " Half length along x (upper endcap, 1098 << " Half length along x (upper endcap, top) = " << fDx4/cm << " cm"
812 << G4endl 1099 << G4endl
813 << "------------------------------------- 1100 << "-----------------------------------------------------------\n";
814 os.precision(oldprc); 1101 os.precision(oldprc);
815 1102
816 return os; 1103 return os;
817 } 1104 }
818 1105
819 1106
820 //============================================ 1107 //=====================================================================
821 //* DiscribeYourselfTo ----------------------- 1108 //* DiscribeYourselfTo ------------------------------------------------
822 1109
823 void G4VTwistedFaceted::DescribeYourselfTo (G4 1110 void G4VTwistedFaceted::DescribeYourselfTo (G4VGraphicsScene& scene) const
824 { 1111 {
825 scene.AddSolid (*this); 1112 scene.AddSolid (*this);
826 } 1113 }
827 1114
828 <<
829 //============================================ 1115 //=====================================================================
830 //* GetExtent -------------------------------- 1116 //* GetExtent ---------------------------------------------------------
831 1117
832 G4VisExtent G4VTwistedFaceted::GetExtent() con 1118 G4VisExtent G4VTwistedFaceted::GetExtent() const
833 { 1119 {
834 G4double maxRad = std::sqrt( fDx*fDx + fDy*f 1120 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
835 1121
836 return { -maxRad, maxRad , << 1122 return G4VisExtent(-maxRad, maxRad ,
837 -maxRad, maxRad , << 1123 -maxRad, maxRad ,
838 -fDz, fDz }; << 1124 -fDz, fDz );
>> 1125 }
>> 1126
>> 1127
>> 1128 //=====================================================================
>> 1129 //* CreateNUBS --------------------------------------------------------
>> 1130
>> 1131 G4NURBS* G4VTwistedFaceted::CreateNURBS () const
>> 1132 {
>> 1133 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
>> 1134
>> 1135 return new G4NURBStube(maxRad, maxRad, fDz);
>> 1136 // Tube for now!!!
839 } 1137 }
840 1138
841 1139
842 //============================================ 1140 //=====================================================================
843 //* CreateSurfaces --------------------------- 1141 //* CreateSurfaces ----------------------------------------------------
844 1142
845 void G4VTwistedFaceted::CreateSurfaces() 1143 void G4VTwistedFaceted::CreateSurfaces()
846 { 1144 {
847 1145
848 // create 6 surfaces of TwistedTub. 1146 // create 6 surfaces of TwistedTub.
849 1147
850 if ( fDx1 == fDx2 && fDx3 == fDx4 ) // sp 1148 if ( fDx1 == fDx2 && fDx3 == fDx4 ) // special case : Box
851 { 1149 {
852 fSide0 = new G4TwistBoxSide("0deg", fPhi << 1150 fSide0 = new G4TwistBoxSide("0deg", fPhiTwist, fDz, fTheta, fPhi,
853 fDy1, fDx1, fDx1, fD << 1151 fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 0.*deg);
854 fSide180 = new G4TwistBoxSide("180deg", fP 1152 fSide180 = new G4TwistBoxSide("180deg", fPhiTwist, fDz, fTheta, fPhi+pi,
855 fDy1, fDx1, fDx1, fD << 1153 fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 180.*deg);
856 } 1154 }
857 else // default general case 1155 else // default general case
858 { 1156 {
859 fSide0 = new G4TwistTrapAlphaSide("0deg" 1157 fSide0 = new G4TwistTrapAlphaSide("0deg" ,fPhiTwist, fDz, fTheta,
860 fPhi, fDy1, fDx1, fDx2, 1158 fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
861 fSide180 = new G4TwistTrapAlphaSide("180de 1159 fSide180 = new G4TwistTrapAlphaSide("180deg", fPhiTwist, fDz, fTheta,
862 fPhi+pi, fDy1, fDx2, fDx1, fD 1160 fPhi+pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*deg);
863 } 1161 }
864 1162
865 // create parallel sides 1163 // create parallel sides
866 // 1164 //
867 fSide90 = new G4TwistTrapParallelSide("90deg 1165 fSide90 = new G4TwistTrapParallelSide("90deg", fPhiTwist, fDz, fTheta,
868 fPhi, fDy1, fDx1, fDx2, 1166 fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
869 fSide270 = new G4TwistTrapParallelSide("270d 1167 fSide270 = new G4TwistTrapParallelSide("270deg", fPhiTwist, fDz, fTheta,
870 fPhi+pi, fDy1, fDx2, fDx1, fD 1168 fPhi+pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*deg);
871 1169
872 // create endcaps 1170 // create endcaps
873 // 1171 //
874 fUpperEndcap = new G4TwistTrapFlatSide("Uppe 1172 fUpperEndcap = new G4TwistTrapFlatSide("UpperCap",fPhiTwist, fDx3, fDx4, fDy2,
875 fDz, fAlph 1173 fDz, fAlph, fPhi, fTheta, 1 );
876 fLowerEndcap = new G4TwistTrapFlatSide("Lowe 1174 fLowerEndcap = new G4TwistTrapFlatSide("LowerCap",fPhiTwist, fDx1, fDx2, fDy1,
877 fDz, fAlph 1175 fDz, fAlph, fPhi, fTheta, -1 );
878 1176
879 // Set neighbour surfaces 1177 // Set neighbour surfaces
880 1178
881 fSide0->SetNeighbours( fSide270 , fLowerEnd 1179 fSide0->SetNeighbours( fSide270 , fLowerEndcap , fSide90 , fUpperEndcap );
882 fSide90->SetNeighbours( fSide0 , fLowerEnd 1180 fSide90->SetNeighbours( fSide0 , fLowerEndcap , fSide180 , fUpperEndcap );
883 fSide180->SetNeighbours(fSide90 , fLowerEnd 1181 fSide180->SetNeighbours(fSide90 , fLowerEndcap , fSide270 , fUpperEndcap );
884 fSide270->SetNeighbours(fSide180 , fLowerEnd 1182 fSide270->SetNeighbours(fSide180 , fLowerEndcap , fSide0 , fUpperEndcap );
885 fUpperEndcap->SetNeighbours( fSide180, fSide 1183 fUpperEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
886 fLowerEndcap->SetNeighbours( fSide180, fSide 1184 fLowerEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
887 1185
888 } 1186 }
889 1187
890 //============================================ <<
891 //* GetCubicVolume --------------------------- <<
892 <<
893 G4double G4VTwistedFaceted::GetCubicVolume() <<
894 { <<
895 if(fCubicVolume == 0.) <<
896 { <<
897 fCubicVolume = ((fDx1 + fDx2 + fDx3 + fDx4 <<
898 (fDx4 + fDx3 - fDx2 - fDx1 <<
899 } <<
900 return fCubicVolume; <<
901 } <<
902 <<
903 //============================================ <<
904 //* GetLateralFaceArea ----------------------- <<
905 <<
906 G4double <<
907 G4VTwistedFaceted::GetLateralFaceArea(const G4 <<
908 const G4 <<
909 const G4 <<
910 const G4 <<
911 { <<
912 constexpr G4int NSTEP = 100; <<
913 constexpr G4double dt = 1./NSTEP; <<
914 <<
915 G4double h = 2.*fDz; <<
916 G4double hh = h*h; <<
917 G4double hTanTheta = h*std::tan(fTheta); <<
918 G4double x1 = p1.x(); <<
919 G4double y1 = p1.y(); <<
920 G4double x21 = p2.x() - p1.x(); <<
921 G4double y21 = p2.y() - p1.y(); <<
922 G4double x31 = p3.x() - p1.x(); <<
923 G4double y31 = p3.y() - p1.y(); <<
924 G4double x42 = p4.x() - p2.x(); <<
925 G4double y42 = p4.y() - p2.y(); <<
926 G4double x43 = p4.x() - p3.x(); <<
927 G4double y43 = p4.y() - p3.y(); <<
928 <<
929 // check if face is plane (just in case) <<
930 G4double lmax = std::max(std::max(std::abs(x <<
931 std::max(std::abs(x <<
932 G4double eps = lmax*kCarTolerance; <<
933 if (std::abs(fPhiTwist) < kCarTolerance && <<
934 std::abs(x21*y43 - y21*x43) < eps) <<
935 { <<
936 G4double x0 = hTanTheta*std::cos(fPhi); <<
937 G4double y0 = hTanTheta*std::sin(fPhi); <<
938 G4ThreeVector A(p4.x() - p1.x() + x0, p4.y <<
939 G4ThreeVector B(p3.x() - p2.x() + x0, p3.y <<
940 return (A.cross(B)).mag()*0.5; <<
941 } <<
942 <<
943 // twisted face <<
944 G4double area = 0.; <<
945 for (G4int i = 0; i < NSTEP; ++i) <<
946 { <<
947 G4double t = (i + 0.5)*dt; <<
948 G4double I = x21 + (x42 - x31)*t; <<
949 G4double J = y21 + (y42 - y31)*t; <<
950 G4double II = I*I; <<
951 G4double JJ = J*J; <<
952 G4double IIJJ = hh*(I*I + J*J); <<
953 <<
954 G4double ang = fPhi + fPhiTwist*(0.5 - t); <<
955 G4double A = fPhiTwist*(II + JJ) + x21*y43 <<
956 G4double B = fPhiTwist*(I*(x1 + x31*t) + J <<
957 hTanTheta*(I*std::sin(ang) - J*std::cos( <<
958 (I*y31 - J*x31); <<
959 <<
960 G4double invAA = 1./(A*A); <<
961 G4double sqrtAA = 2.*std::abs(A); <<
962 G4double invSqrtAA = 1./sqrtAA; <<
963 <<
964 G4double aa = A*A; <<
965 G4double bb = 2.*A*B; <<
966 G4double cc = IIJJ + B*B; <<
967 <<
968 G4double R1 = std::sqrt(aa + bb + cc); <<
969 G4double R0 = std::sqrt(cc); <<
970 G4double log1 = std::log(std::abs(sqrtAA*R <<
971 G4double log0 = std::log(std::abs(sqrtAA*R <<
972 <<
973 area += 0.5*R1 + 0.25*bb*invAA*(R1 - R0) + <<
974 } <<
975 return area*dt; <<
976 } <<
977 <<
978 //============================================ <<
979 //* GetSurfaceArea --------------------------- <<
980 <<
981 G4double G4VTwistedFaceted::GetSurfaceArea() <<
982 { <<
983 if (fSurfaceArea == 0.) <<
984 { <<
985 G4TwoVector vv[8]; <<
986 vv[0] = G4TwoVector(-fDx1 - fDy1*fTAlph,-f <<
987 vv[1] = G4TwoVector( fDx1 - fDy1*fTAlph,-f <<
988 vv[2] = G4TwoVector(-fDx2 + fDy1*fTAlph, f <<
989 vv[3] = G4TwoVector( fDx2 + fDy1*fTAlph, f <<
990 vv[4] = G4TwoVector(-fDx3 - fDy2*fTAlph,-f <<
991 vv[5] = G4TwoVector( fDx3 - fDy2*fTAlph,-f <<
992 vv[6] = G4TwoVector(-fDx4 + fDy2*fTAlph, f <<
993 vv[7] = G4TwoVector( fDx4 + fDy2*fTAlph, f <<
994 fSurfaceArea = 2.*(fDy1*(fDx1 + fDx2) + fD <<
995 GetLateralFaceArea(vv[0], vv[1], vv[4], <<
996 GetLateralFaceArea(vv[1], vv[3], vv[5], <<
997 GetLateralFaceArea(vv[3], vv[2], vv[7], <<
998 GetLateralFaceArea(vv[2], vv[0], vv[6], <<
999 } <<
1000 return fSurfaceArea; <<
1001 } <<
1002 1188
1003 //=========================================== 1189 //=====================================================================
1004 //* GetEntityType --------------------------- 1190 //* GetEntityType -----------------------------------------------------
1005 1191
1006 G4GeometryType G4VTwistedFaceted::GetEntityTy 1192 G4GeometryType G4VTwistedFaceted::GetEntityType() const
1007 { 1193 {
1008 return {"G4VTwistedFaceted"}; << 1194 return G4String("G4VTwistedFaceted");
1009 } 1195 }
1010 1196
1011 1197
1012 //=========================================== 1198 //=====================================================================
1013 //* GetPolyhedron --------------------------- 1199 //* GetPolyhedron -----------------------------------------------------
1014 1200
1015 G4Polyhedron* G4VTwistedFaceted::GetPolyhedro 1201 G4Polyhedron* G4VTwistedFaceted::GetPolyhedron() const
1016 { 1202 {
1017 if (fpPolyhedron == nullptr || << 1203 if (!fpPolyhedron ||
1018 fRebuildPolyhedron || <<
1019 fpPolyhedron->GetNumberOfRotationStepsA 1204 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1020 fpPolyhedron->GetNumberOfRotationSteps( 1205 fpPolyhedron->GetNumberOfRotationSteps())
1021 { 1206 {
1022 G4AutoLock l(&polyhedronMutex); <<
1023 delete fpPolyhedron; 1207 delete fpPolyhedron;
1024 fpPolyhedron = CreatePolyhedron(); 1208 fpPolyhedron = CreatePolyhedron();
1025 fRebuildPolyhedron = false; <<
1026 l.unlock(); <<
1027 } 1209 }
1028 1210
1029 return fpPolyhedron; 1211 return fpPolyhedron;
1030 } 1212 }
1031 1213
1032 1214
1033 //=========================================== 1215 //=====================================================================
1034 //* GetPointInSolid ------------------------- 1216 //* GetPointInSolid ---------------------------------------------------
1035 1217
1036 G4ThreeVector G4VTwistedFaceted::GetPointInSo 1218 G4ThreeVector G4VTwistedFaceted::GetPointInSolid(G4double z) const
1037 { 1219 {
1038 1220
1039 1221
1040 // this routine is only used for a test 1222 // this routine is only used for a test
1041 // can be deleted ... 1223 // can be deleted ...
1042 1224
1043 if ( z == fDz ) z -= 0.1*fDz ; 1225 if ( z == fDz ) z -= 0.1*fDz ;
1044 if ( z == -fDz ) z += 0.1*fDz ; 1226 if ( z == -fDz ) z += 0.1*fDz ;
1045 1227
1046 G4double phi = z/(2*fDz)*fPhiTwist ; 1228 G4double phi = z/(2*fDz)*fPhiTwist ;
1047 1229
1048 return { fdeltaX * phi/fPhiTwist, fdeltaY * << 1230 return G4ThreeVector(fdeltaX * phi/fPhiTwist, fdeltaY * phi/fPhiTwist, z ) ;
1049 } 1231 }
1050 1232
1051 1233
1052 //=========================================== 1234 //=====================================================================
1053 //* GetPointOnSurface ----------------------- 1235 //* GetPointOnSurface -------------------------------------------------
1054 1236
1055 G4ThreeVector G4VTwistedFaceted::GetPointOnSu 1237 G4ThreeVector G4VTwistedFaceted::GetPointOnSurface() const
1056 { 1238 {
1057 1239
1058 G4double phi = G4RandFlat::shoot(-fPhiTwist << 1240 G4double phi = G4RandFlat::shoot(-fPhiTwist/2.,fPhiTwist/2.);
1059 G4double u , umin, umax ; // variable for 1241 G4double u , umin, umax ; // variable for twisted surfaces
1060 G4double y ; // variable for 1242 G4double y ; // variable for flat surface (top and bottom)
1061 1243
1062 // Compute the areas. Attention: Only corre 1244 // Compute the areas. Attention: Only correct for trapezoids
1063 // where the twisting is done along the z-a 1245 // where the twisting is done along the z-axis. In the general case
1064 // the computed surface area is more diffic 1246 // the computed surface area is more difficult. However this simplification
1065 // does not affect the tracking through the 1247 // does not affect the tracking through the solid.
1066 1248
1067 G4double a1 = fSide0->GetSurfaceArea(); << 1249 G4double a1 = fSide0->GetSurfaceArea();
1068 G4double a2 = fSide90->GetSurfaceArea(); << 1250 G4double a2 = fSide90->GetSurfaceArea();
1069 G4double a3 = fSide180->GetSurfaceArea() ; << 1251 G4double a3 = fSide180->GetSurfaceArea() ;
1070 G4double a4 = fSide270->GetSurfaceArea() ; << 1252 G4double a4 = fSide270->GetSurfaceArea() ;
1071 G4double a5 = fLowerEndcap->GetSurfaceArea( << 1253 G4double a5 = fLowerEndcap->GetSurfaceArea() ;
1072 G4double a6 = fUpperEndcap->GetSurfaceArea( << 1254 G4double a6 = fUpperEndcap->GetSurfaceArea() ;
1073 1255
1074 #ifdef G4TWISTDEBUG 1256 #ifdef G4TWISTDEBUG
1075 G4cout << "Surface 0 deg = " << a1 << G4e 1257 G4cout << "Surface 0 deg = " << a1 << G4endl ;
1076 G4cout << "Surface 90 deg = " << a2 << G4e 1258 G4cout << "Surface 90 deg = " << a2 << G4endl ;
1077 G4cout << "Surface 180 deg = " << a3 << G4e 1259 G4cout << "Surface 180 deg = " << a3 << G4endl ;
1078 G4cout << "Surface 270 deg = " << a4 << G4e 1260 G4cout << "Surface 270 deg = " << a4 << G4endl ;
1079 G4cout << "Surface Lower = " << a5 << G4e 1261 G4cout << "Surface Lower = " << a5 << G4endl ;
1080 G4cout << "Surface Upper = " << a6 << G4e 1262 G4cout << "Surface Upper = " << a6 << G4endl ;
1081 #endif 1263 #endif
1082 1264
1083 G4double chose = G4RandFlat::shoot(0.,a1 + 1265 G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1084 1266
1085 if(chose < a1) 1267 if(chose < a1)
1086 { 1268 {
>> 1269
1087 umin = fSide0->GetBoundaryMin(phi) ; 1270 umin = fSide0->GetBoundaryMin(phi) ;
1088 umax = fSide0->GetBoundaryMax(phi) ; 1271 umax = fSide0->GetBoundaryMax(phi) ;
1089 u = G4RandFlat::shoot(umin,umax) ; 1272 u = G4RandFlat::shoot(umin,umax) ;
1090 1273
1091 return fSide0->SurfacePoint(phi, u, true 1274 return fSide0->SurfacePoint(phi, u, true) ; // point on 0deg surface
1092 } 1275 }
1093 1276
1094 else if( (chose >= a1) && (chose < a1 + a2 1277 else if( (chose >= a1) && (chose < a1 + a2 ) )
1095 { 1278 {
>> 1279
1096 umin = fSide90->GetBoundaryMin(phi) ; 1280 umin = fSide90->GetBoundaryMin(phi) ;
1097 umax = fSide90->GetBoundaryMax(phi) ; 1281 umax = fSide90->GetBoundaryMax(phi) ;
1098 1282
1099 u = G4RandFlat::shoot(umin,umax) ; 1283 u = G4RandFlat::shoot(umin,umax) ;
1100 1284
1101 return fSide90->SurfacePoint(phi, u, true 1285 return fSide90->SurfacePoint(phi, u, true); // point on 90deg surface
1102 } 1286 }
>> 1287
1103 else if( (chose >= a1 + a2 ) && (chose < a1 1288 else if( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1104 { 1289 {
>> 1290
1105 umin = fSide180->GetBoundaryMin(phi) ; 1291 umin = fSide180->GetBoundaryMin(phi) ;
1106 umax = fSide180->GetBoundaryMax(phi) ; 1292 umax = fSide180->GetBoundaryMax(phi) ;
1107 u = G4RandFlat::shoot(umin,umax) ; 1293 u = G4RandFlat::shoot(umin,umax) ;
1108 1294
1109 return fSide180->SurfacePoint(phi, u, tru << 1295 return fSide180->SurfacePoint(phi, u, true); // point on 180 deg surface
1110 } 1296 }
>> 1297
1111 else if( (chose >= a1 + a2 + a3 ) && (chos 1298 else if( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1112 { 1299 {
>> 1300
1113 umin = fSide270->GetBoundaryMin(phi) ; 1301 umin = fSide270->GetBoundaryMin(phi) ;
1114 umax = fSide270->GetBoundaryMax(phi) ; 1302 umax = fSide270->GetBoundaryMax(phi) ;
1115 u = G4RandFlat::shoot(umin,umax) ; 1303 u = G4RandFlat::shoot(umin,umax) ;
>> 1304
1116 return fSide270->SurfacePoint(phi, u, tru 1305 return fSide270->SurfacePoint(phi, u, true); // point on 270 deg surface
1117 } 1306 }
>> 1307
1118 else if( (chose >= a1 + a2 + a3 + a4 ) && 1308 else if( (chose >= a1 + a2 + a3 + a4 ) && (chose < a1 + a2 + a3 + a4 + a5 ) )
1119 { 1309 {
>> 1310
1120 y = G4RandFlat::shoot(-fDy1,fDy1) ; 1311 y = G4RandFlat::shoot(-fDy1,fDy1) ;
1121 umin = fLowerEndcap->GetBoundaryMin(y) ; 1312 umin = fLowerEndcap->GetBoundaryMin(y) ;
1122 umax = fLowerEndcap->GetBoundaryMax(y) ; 1313 umax = fLowerEndcap->GetBoundaryMax(y) ;
1123 u = G4RandFlat::shoot(umin,umax) ; 1314 u = G4RandFlat::shoot(umin,umax) ;
1124 1315
1125 return fLowerEndcap->SurfacePoint(u,y,tru 1316 return fLowerEndcap->SurfacePoint(u,y,true); // point on lower endcap
1126 } 1317 }
1127 else << 1318 else {
1128 { << 1319
1129 y = G4RandFlat::shoot(-fDy2,fDy2) ; 1320 y = G4RandFlat::shoot(-fDy2,fDy2) ;
1130 umin = fUpperEndcap->GetBoundaryMin(y) ; 1321 umin = fUpperEndcap->GetBoundaryMin(y) ;
1131 umax = fUpperEndcap->GetBoundaryMax(y) ; 1322 umax = fUpperEndcap->GetBoundaryMax(y) ;
1132 u = G4RandFlat::shoot(umin,umax) ; 1323 u = G4RandFlat::shoot(umin,umax) ;
1133 1324
1134 return fUpperEndcap->SurfacePoint(u,y,tru 1325 return fUpperEndcap->SurfacePoint(u,y,true) ; // point on upper endcap
1135 1326
1136 } 1327 }
1137 } 1328 }
1138 1329
1139 1330
1140 //=========================================== 1331 //=====================================================================
1141 //* CreatePolyhedron ------------------------ 1332 //* CreatePolyhedron --------------------------------------------------
1142 1333
1143 G4Polyhedron* G4VTwistedFaceted::CreatePolyhe 1334 G4Polyhedron* G4VTwistedFaceted::CreatePolyhedron () const
1144 { 1335 {
1145 // number of meshes 1336 // number of meshes
1146 const G4int k = 1337 const G4int k =
1147 G4int(G4Polyhedron::GetNumberOfRotationStep << 1338 G4int(G4Polyhedron::GetNumberOfRotationSteps() * fPhiTwist / twopi) + 2;
1148 std::abs(fPhiTwist) / twopi) + 2; <<
1149 const G4int n = k; 1339 const G4int n = k;
1150 1340
1151 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k 1341 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
1152 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1 1342 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
1153 1343
1154 auto ph = new G4Polyhedron; << 1344 G4Polyhedron *ph=new G4Polyhedron;
1155 typedef G4double G4double3[3]; 1345 typedef G4double G4double3[3];
1156 typedef G4int G4int4[4]; 1346 typedef G4int G4int4[4];
1157 auto xyz = new G4double3[nnodes]; // numbe << 1347 G4double3* xyz = new G4double3[nnodes]; // number of nodes
1158 auto faces = new G4int4[nfaces] ; // numbe << 1348 G4int4* faces = new G4int4[nfaces] ; // number of faces
1159 1349
1160 fLowerEndcap->GetFacets(k,k,xyz,faces,0) ; 1350 fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
1161 fUpperEndcap->GetFacets(k,k,xyz,faces,1) ; 1351 fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
1162 fSide270->GetFacets(k,n,xyz,faces,2) ; 1352 fSide270->GetFacets(k,n,xyz,faces,2) ;
1163 fSide0->GetFacets(k,n,xyz,faces,3) ; 1353 fSide0->GetFacets(k,n,xyz,faces,3) ;
1164 fSide90->GetFacets(k,n,xyz,faces,4) ; 1354 fSide90->GetFacets(k,n,xyz,faces,4) ;
1165 fSide180->GetFacets(k,n,xyz,faces,5) ; 1355 fSide180->GetFacets(k,n,xyz,faces,5) ;
1166 1356
1167 ph->createPolyhedron(nnodes,nfaces,xyz,face 1357 ph->createPolyhedron(nnodes,nfaces,xyz,faces);
1168 1358
1169 return ph; 1359 return ph;
1170 } 1360 }
1171 1361