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