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