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