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