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