| 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 //
>> 27 // $Id: G4Trap.cc 101121 2016-11-07 09:18:01Z gcosmo $
>> 28 //
>> 29 // class G4Trap
>> 30 //
26 // Implementation for G4Trap class 31 // Implementation for G4Trap class
27 // 32 //
28 // 21.03.95 P.Kent: Modified for `tolerant' ge << 33 // History:
>> 34 //
>> 35 // 23.09.16 E.Tcherniaev: added Extent(pmin,pmax),
>> 36 // use G4BoundingEnvelope for CalculateExtent(),
>> 37 // removed CreateRotatedVertices()
>> 38 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal
>> 39 // 26.04.05 V.Grichine: new SurfaceNormal is default
>> 40 // 19.04.05 V.Grichine: bug fixed in G4Trap("name",G4ThreeVector[8] vp)
>> 41 // 12.12.04 V.Grichine: SurfaceNormal with edges/vertices
>> 42 // 15.11.04 V.Grichine: bug fixed in G4Trap("name",G4ThreeVector[8] vp)
>> 43 // 13.12.99 V.Grichine: bug fixed in DistanceToIn(p,v)
>> 44 // 19.11.99 V.Grichine: kUndef was added to Eside enum
>> 45 // 04.06.99 S.Giani: Fixed CalculateExtent in rotated case.
>> 46 // 08.12.97 J.Allison: Added "nominal" constructor and method SetAllParameters.
>> 47 // 01.11.96 V.Grichine: Costructor for Right Angular Wedge from STEP, G4Trd/Para
29 // 09.09.96 V.Grichine: Final modifications be 48 // 09.09.96 V.Grichine: Final modifications before to commit
30 // 08.12.97 J.Allison: Added "nominal" constru << 49 // 21.03.95 P.Kent: Modified for `tolerant' geometry
31 // 28.04.05 V.Grichine: new SurfaceNormal acco << 50 //
32 // 18.04.17 E.Tcherniaev: complete revision, s << 51 ////////////////////////////////////////////////////////////////////////////////////
33 // ------------------------------------------- <<
34 52
35 #include "G4Trap.hh" 53 #include "G4Trap.hh"
36 54
37 #if !defined(G4GEOM_USE_UTRAP) 55 #if !defined(G4GEOM_USE_UTRAP)
38 56
39 #include "globals.hh" 57 #include "globals.hh"
40 #include "G4GeomTools.hh" <<
41 58
42 #include "G4VoxelLimits.hh" 59 #include "G4VoxelLimits.hh"
43 #include "G4AffineTransform.hh" 60 #include "G4AffineTransform.hh"
44 #include "G4BoundingEnvelope.hh" 61 #include "G4BoundingEnvelope.hh"
45 62
46 #include "G4VPVParameterisation.hh" 63 #include "G4VPVParameterisation.hh"
47 64
48 #include "G4QuickRand.hh" << 65 #include "Randomize.hh"
49 66
50 #include "G4VGraphicsScene.hh" 67 #include "G4VGraphicsScene.hh"
51 #include "G4Polyhedron.hh" 68 #include "G4Polyhedron.hh"
52 69
53 using namespace CLHEP; 70 using namespace CLHEP;
54 71
>> 72 ////////////////////////////////////////////////////////////////////////
>> 73 //
>> 74 // Accuracy of coplanarity
>> 75
>> 76 const G4double kCoplanar_Tolerance = 1E-4 ;
>> 77
55 ////////////////////////////////////////////// 78 //////////////////////////////////////////////////////////////////////////
56 // 79 //
57 // Constructor - check and set half-widths as << 80 // Private enum: Not for external use
>> 81
>> 82 enum Eside {kUndef,ks0,ks1,ks2,ks3,kPZ,kMZ};
>> 83
>> 84 //////////////////////////////////////////////////////////////////////////
>> 85 //
>> 86 // Constructor - check and set half-widths as well as angles:
58 // final check of coplanarity 87 // final check of coplanarity
59 88
60 G4Trap::G4Trap( const G4String& pName, 89 G4Trap::G4Trap( const G4String& pName,
61 G4double pDz, 90 G4double pDz,
62 G4double pTheta, G4doubl 91 G4double pTheta, G4double pPhi,
63 G4double pDy1, G4double 92 G4double pDy1, G4double pDx1, G4double pDx2,
64 G4double pAlp1, 93 G4double pAlp1,
65 G4double pDy2, G4double 94 G4double pDy2, G4double pDx3, G4double pDx4,
66 G4double pAlp2 ) << 95 G4double pAlp2)
67 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 96 : G4CSGSolid(pName)
68 { 97 {
69 fDz = pDz; << 98 if ( pDz <= 0 || pDy1 <= 0 || pDx1 <= 0 ||
70 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi << 99 pDx2 <= 0 || pDy2 <= 0 || pDx3 <= 0 || pDx4 <= 0 )
71 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi << 100 {
>> 101 std::ostringstream message;
>> 102 message << "Invalid length parameters for Solid: " << GetName() << G4endl
>> 103 << " X - "
>> 104 << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
>> 105 << " Y - " << pDy1 << ", " << pDy2 << G4endl
>> 106 << " Z - " << pDz;
>> 107 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 108 FatalException, message);
>> 109 }
72 110
73 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp << 111 fDz=pDz;
74 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp << 112 fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
>> 113 fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
>> 114
>> 115 fDy1=pDy1;
>> 116 fDx1=pDx1;
>> 117 fDx2=pDx2;
>> 118 fTalpha1=std::tan(pAlp1);
>> 119
>> 120 fDy2=pDy2;
>> 121 fDx3=pDx3;
>> 122 fDx4=pDx4;
>> 123 fTalpha2=std::tan(pAlp2);
75 124
76 CheckParameters(); <<
77 MakePlanes(); 125 MakePlanes();
78 } 126 }
79 127
80 ////////////////////////////////////////////// << 128 ////////////////////////////////////////////////////////////////////////////
81 // 129 //
82 // Constructor - Design of trapezoid based on << 130 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
83 // which are its vertices. Checking of planari << 131 // which are its vertices. Checking of planarity with preparation of
84 // fPlanes[] and than calculation of other mem 132 // fPlanes[] and than calculation of other members
85 133
86 G4Trap::G4Trap( const G4String& pName, 134 G4Trap::G4Trap( const G4String& pName,
87 const G4ThreeVector pt[8] ) 135 const G4ThreeVector pt[8] )
88 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 136 : G4CSGSolid(pName)
89 { 137 {
>> 138 G4bool good;
>> 139
90 // Start with check of centering - the cente 140 // Start with check of centering - the center of gravity trap line
91 // should cross the origin of frame 141 // should cross the origin of frame
92 // << 142
93 if ( pt[0].z() >= 0 << 143 if (!( pt[0].z() < 0
94 || pt[0].z() != pt[1].z() << 144 && pt[0].z() == pt[1].z() && pt[0].z() == pt[2].z()
95 || pt[0].z() != pt[2].z() << 145 && pt[0].z() == pt[3].z()
96 || pt[0].z() != pt[3].z() << 146 && pt[4].z() > 0
97 << 147 && pt[4].z() == pt[5].z() && pt[4].z() == pt[6].z()
98 || pt[4].z() <= 0 << 148 && pt[4].z() == pt[7].z()
99 || pt[4].z() != pt[5].z() << 149 && std::fabs( pt[0].z() + pt[4].z() ) < kCarTolerance
100 || pt[4].z() != pt[6].z() << 150 && pt[0].y() == pt[1].y() && pt[2].y() == pt[3].y()
101 || pt[4].z() != pt[7].z() << 151 && pt[4].y() == pt[5].y() && pt[6].y() == pt[7].y()
102 << 152 && std::fabs( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) < kCarTolerance
103 || std::fabs( pt[0].z() + pt[4].z() ) << 153 && std::fabs( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() +
104 << 154 pt[2].x() + pt[3].x() + pt[6].x() + pt[7].x() ) < kCarTolerance ) )
105 || pt[0].y() != pt[1].y() <<
106 || pt[2].y() != pt[3].y() <<
107 || pt[4].y() != pt[5].y() <<
108 || pt[6].y() != pt[7].y() <<
109 <<
110 || std::fabs(pt[0].y()+pt[2].y()+pt[4] <<
111 || std::fabs(pt[0].x()+pt[1].x()+pt[4] <<
112 pt[2].x()+pt[3].x()+pt[6] <<
113 { 155 {
114 std::ostringstream message; 156 std::ostringstream message;
115 message << "Invalid vertice coordinates fo 157 message << "Invalid vertice coordinates for Solid: " << GetName();
116 G4Exception("G4Trap::G4Trap()", "GeomSolid 158 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
117 FatalException, message); 159 FatalException, message);
118 } 160 }
>> 161
>> 162 // Bottom side with normal approx. -Y
>> 163
>> 164 good = MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
119 165
120 // Set parameters << 166 if (!good)
121 // << 167 {
122 fDz = (pt[7]).z(); << 168 DumpInfo();
>> 169 G4Exception("G4Trap::G4Trap()", "GeomSolids0002", FatalException,
>> 170 "Face at ~-Y not planar.");
>> 171 }
123 172
>> 173 // Top side with normal approx. +Y
>> 174
>> 175 good = MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
>> 176
>> 177 if (!good)
>> 178 {
>> 179 std::ostringstream message;
>> 180 message << "Face at ~+Y not planar for Solid: " << GetName();
>> 181 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 182 FatalException, message);
>> 183 }
>> 184
>> 185 // Front side with normal approx. -X
>> 186
>> 187 good = MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
>> 188
>> 189 if (!good)
>> 190 {
>> 191 std::ostringstream message;
>> 192 message << "Face at ~-X not planar for Solid: " << GetName();
>> 193 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 194 FatalException, message);
>> 195 }
>> 196
>> 197 // Back side iwth normal approx. +X
>> 198
>> 199 good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
>> 200 if (!good)
>> 201 {
>> 202 std::ostringstream message;
>> 203 message << "Face at ~+X not planar for Solid: " << GetName();
>> 204 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 205 FatalException, message);
>> 206 }
>> 207 fDz = (pt[7]).z() ;
>> 208
124 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5; 209 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5;
125 fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5; 210 fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5;
126 fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5; 211 fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5;
127 fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]). 212 fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy1;
128 213
129 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5; 214 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5;
130 fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5; 215 fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5;
131 fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5; 216 fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5;
132 fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]). 217 fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy2;
133 218
134 fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx 219 fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx3)/fDz;
135 fTthetaSphi = ((pt[4]).y()+fDy2)/fDz; 220 fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;
136 <<
137 CheckParameters(); <<
138 MakePlanes(pt); <<
139 } 221 }
140 222
141 ////////////////////////////////////////////// << 223 //////////////////////////////////////////////////////////////////////////////
142 // 224 //
143 // Constructor for Right Angular Wedge from ST 225 // Constructor for Right Angular Wedge from STEP
144 226
145 G4Trap::G4Trap( const G4String& pName, 227 G4Trap::G4Trap( const G4String& pName,
146 G4double pZ, 228 G4double pZ,
147 G4double pY, 229 G4double pY,
148 G4double pX, G4double pL 230 G4double pX, G4double pLTX )
149 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 231 : G4CSGSolid(pName)
150 { 232 {
151 fDz = 0.5*pZ; fTthetaCphi = 0; fTthetaSphi << 233 G4bool good;
152 fDy1 = 0.5*pY; fDx1 = 0.5*pX; fDx2 = 0.5*pLT <<
153 fDy2 = fDy1; fDx3 = fDx1; fDx4 = fDx2; <<
154 234
155 CheckParameters(); << 235 if ( pZ<=0 || pY<=0 || pX<=0 || pLTX<=0 || pLTX>pX )
156 MakePlanes(); << 236 {
>> 237 std::ostringstream message;
>> 238 message << "Invalid length parameters for Solid: " << GetName();
>> 239 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 240 FatalException, message);
>> 241 }
>> 242
>> 243 fDz = 0.5*pZ ;
>> 244 fTthetaCphi = 0 ;
>> 245 fTthetaSphi = 0 ;
>> 246
>> 247 fDy1 = 0.5*pY;
>> 248 fDx1 = 0.5*pX ;
>> 249 fDx2 = 0.5*pLTX;
>> 250 fTalpha1 = 0.5*(pLTX - pX)/pY;
>> 251
>> 252 fDy2 = fDy1 ;
>> 253 fDx3 = fDx1;
>> 254 fDx4 = fDx2 ;
>> 255 fTalpha2 = fTalpha1 ;
>> 256
>> 257 G4ThreeVector pt[8] ;
>> 258
>> 259 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
>> 260 -fDz*fTthetaSphi-fDy1,-fDz);
>> 261 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
>> 262 -fDz*fTthetaSphi-fDy1,-fDz);
>> 263 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
>> 264 -fDz*fTthetaSphi+fDy1,-fDz);
>> 265 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
>> 266 -fDz*fTthetaSphi+fDy1,-fDz);
>> 267 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
>> 268 +fDz*fTthetaSphi-fDy2,+fDz);
>> 269 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
>> 270 +fDz*fTthetaSphi-fDy2,+fDz);
>> 271 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
>> 272 +fDz*fTthetaSphi+fDy2,+fDz);
>> 273 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
>> 274 +fDz*fTthetaSphi+fDy2,+fDz);
>> 275
>> 276 // Bottom side with normal approx. -Y
>> 277 //
>> 278 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
>> 279 if (!good)
>> 280 {
>> 281 std::ostringstream message;
>> 282 message << "Face at ~-Y not planar for Solid: " << GetName();
>> 283 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 284 FatalException, message);
>> 285 }
>> 286
>> 287 // Top side with normal approx. +Y
>> 288 //
>> 289 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
>> 290 if (!good)
>> 291 {
>> 292 std::ostringstream message;
>> 293 message << "Face at ~+Y not planar for Solid: " << GetName();
>> 294 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 295 FatalException, message);
>> 296 }
>> 297
>> 298 // Front side with normal approx. -X
>> 299 //
>> 300 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
>> 301 if (!good)
>> 302 {
>> 303 std::ostringstream message;
>> 304 message << "Face at ~-X not planar for Solid: " << GetName();
>> 305 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 306 FatalException, message);
>> 307 }
>> 308
>> 309 // Back side iwth normal approx. +X
>> 310 //
>> 311 good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
>> 312 if (!good)
>> 313 {
>> 314 std::ostringstream message;
>> 315 message << "Face at ~+X not planar for Solid: " << GetName();
>> 316 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 317 FatalException, message);
>> 318 }
157 } 319 }
158 320
159 ////////////////////////////////////////////// << 321 ///////////////////////////////////////////////////////////////////////////////
160 // 322 //
161 // Constructor for G4Trd 323 // Constructor for G4Trd
162 324
163 G4Trap::G4Trap( const G4String& pName, 325 G4Trap::G4Trap( const G4String& pName,
164 G4double pDx1, G4double 326 G4double pDx1, G4double pDx2,
165 G4double pDy1, G4double 327 G4double pDy1, G4double pDy2,
166 G4double pDz ) 328 G4double pDz )
167 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 329 : G4CSGSolid(pName)
168 { 330 {
169 fDz = pDz; fTthetaCphi = 0; fTthetaSphi = << 331 G4bool good;
170 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx1; fTalp <<
171 fDy2 = pDy2; fDx3 = pDx2; fDx4 = pDx2; fTalp <<
172 332
173 CheckParameters(); << 333 if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 )
174 MakePlanes(); << 334 {
>> 335 std::ostringstream message;
>> 336 message << "Invalid length parameters for Solid: " << GetName();
>> 337 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 338 FatalException, message);
>> 339 }
>> 340
>> 341 fDz = pDz;
>> 342 fTthetaCphi = 0 ;
>> 343 fTthetaSphi = 0 ;
>> 344
>> 345 fDy1 = pDy1 ;
>> 346 fDx1 = pDx1 ;
>> 347 fDx2 = pDx1 ;
>> 348 fTalpha1 = 0 ;
>> 349
>> 350 fDy2 = pDy2 ;
>> 351 fDx3 = pDx2 ;
>> 352 fDx4 = pDx2 ;
>> 353 fTalpha2 = 0 ;
>> 354
>> 355 G4ThreeVector pt[8] ;
>> 356
>> 357 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
>> 358 -fDz*fTthetaSphi-fDy1,-fDz);
>> 359 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
>> 360 -fDz*fTthetaSphi-fDy1,-fDz);
>> 361 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
>> 362 -fDz*fTthetaSphi+fDy1,-fDz);
>> 363 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
>> 364 -fDz*fTthetaSphi+fDy1,-fDz);
>> 365 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
>> 366 +fDz*fTthetaSphi-fDy2,+fDz);
>> 367 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
>> 368 +fDz*fTthetaSphi-fDy2,+fDz);
>> 369 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
>> 370 +fDz*fTthetaSphi+fDy2,+fDz);
>> 371 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
>> 372 +fDz*fTthetaSphi+fDy2,+fDz);
>> 373
>> 374 // Bottom side with normal approx. -Y
>> 375 //
>> 376 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
>> 377 if (!good)
>> 378 {
>> 379 std::ostringstream message;
>> 380 message << "Face at ~-Y not planar for Solid: " << GetName();
>> 381 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 382 FatalException, message);
>> 383 }
>> 384
>> 385 // Top side with normal approx. +Y
>> 386 //
>> 387 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
>> 388 if (!good)
>> 389 {
>> 390 std::ostringstream message;
>> 391 message << "Face at ~+Y not planar for Solid: " << GetName();
>> 392 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 393 FatalException, message);
>> 394 }
>> 395
>> 396 // Front side with normal approx. -X
>> 397 //
>> 398 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
>> 399 if (!good)
>> 400 {
>> 401 std::ostringstream message;
>> 402 message << "Face at ~-X not planar for Solid: " << GetName();
>> 403 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 404 FatalException, message);
>> 405 }
>> 406
>> 407 // Back side iwth normal approx. +X
>> 408 //
>> 409 good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
>> 410 if (!good)
>> 411 {
>> 412 std::ostringstream message;
>> 413 message << "Face at ~+X not planar for Solid: " << GetName();
>> 414 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 415 FatalException, message);
>> 416 }
175 } 417 }
176 418
177 ////////////////////////////////////////////// << 419 ////////////////////////////////////////////////////////////////////////////
178 // 420 //
179 // Constructor for G4Para 421 // Constructor for G4Para
180 422
181 G4Trap::G4Trap( const G4String& pName, 423 G4Trap::G4Trap( const G4String& pName,
182 G4double pDx, G4double p 424 G4double pDx, G4double pDy,
183 G4double pDz, 425 G4double pDz,
184 G4double pAlpha, 426 G4double pAlpha,
185 G4double pTheta, G4doubl 427 G4double pTheta, G4double pPhi )
186 : G4CSGSolid(pName), halfCarTolerance(0.5*kC << 428 : G4CSGSolid(pName)
187 { 429 {
188 fDz = pDz; << 430 G4bool good;
189 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi <<
190 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi <<
191 431
192 fDy1 = pDy; fDx1 = pDx; fDx2 = pDx; fTalpha1 << 432 if ( pDz<=0 || pDy<=0 || pDx<=0 )
193 fDy2 = pDy; fDx3 = pDx; fDx4 = pDx; fTalpha2 << 433 {
>> 434 std::ostringstream message;
>> 435 message << "Invalid length parameters for Solid: " << GetName();
>> 436 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 437 FatalException, message);
>> 438 }
194 439
195 CheckParameters(); << 440 fDz = pDz ;
196 MakePlanes(); << 441 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi) ;
>> 442 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi) ;
>> 443
>> 444 fDy1 = pDy ;
>> 445 fDx1 = pDx ;
>> 446 fDx2 = pDx ;
>> 447 fTalpha1 = std::tan(pAlpha) ;
>> 448
>> 449 fDy2 = pDy ;
>> 450 fDx3 = pDx ;
>> 451 fDx4 = pDx ;
>> 452 fTalpha2 = fTalpha1 ;
>> 453
>> 454 G4ThreeVector pt[8] ;
>> 455
>> 456 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
>> 457 -fDz*fTthetaSphi-fDy1,-fDz);
>> 458 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
>> 459 -fDz*fTthetaSphi-fDy1,-fDz);
>> 460 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
>> 461 -fDz*fTthetaSphi+fDy1,-fDz);
>> 462 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
>> 463 -fDz*fTthetaSphi+fDy1,-fDz);
>> 464 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
>> 465 +fDz*fTthetaSphi-fDy2,+fDz);
>> 466 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
>> 467 +fDz*fTthetaSphi-fDy2,+fDz);
>> 468 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
>> 469 +fDz*fTthetaSphi+fDy2,+fDz);
>> 470 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
>> 471 +fDz*fTthetaSphi+fDy2,+fDz);
>> 472
>> 473 // Bottom side with normal approx. -Y
>> 474 //
>> 475 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
>> 476 if (!good)
>> 477 {
>> 478 std::ostringstream message;
>> 479 message << "Face at ~-Y not planar for Solid: " << GetName();
>> 480 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 481 FatalException, message);
>> 482 }
>> 483
>> 484 // Top side with normal approx. +Y
>> 485 //
>> 486 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
>> 487 if (!good)
>> 488 {
>> 489 std::ostringstream message;
>> 490 message << "Face at ~+Y not planar for Solid: " << GetName();
>> 491 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 492 FatalException, message);
>> 493 }
>> 494
>> 495 // Front side with normal approx. -X
>> 496 //
>> 497 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
>> 498 if (!good)
>> 499 {
>> 500 std::ostringstream message;
>> 501 message << "Face at ~-X not planar for Solid: " << GetName();
>> 502 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 503 FatalException, message);
>> 504 }
>> 505
>> 506 // Back side iwth normal approx. +X
>> 507 //
>> 508 good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
>> 509 if (!good)
>> 510 {
>> 511 std::ostringstream message;
>> 512 message << "Face at ~+X not planar for Solid: " << GetName();
>> 513 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
>> 514 FatalException, message);
>> 515 }
197 } 516 }
198 517
199 ////////////////////////////////////////////// << 518 ///////////////////////////////////////////////////////////////////////////
200 // 519 //
201 // Nominal constructor for G4Trap whose parame 520 // Nominal constructor for G4Trap whose parameters are to be set by
202 // a G4VParamaterisation later. Check and set 521 // a G4VParamaterisation later. Check and set half-widths as well as
203 // angles: final check of coplanarity 522 // angles: final check of coplanarity
204 523
205 G4Trap::G4Trap( const G4String& pName ) 524 G4Trap::G4Trap( const G4String& pName )
206 : G4CSGSolid (pName), halfCarTolerance(0.5*k << 525 : G4CSGSolid (pName), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
207 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), <<
208 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.) 526 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
209 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.) 527 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
210 { 528 {
211 MakePlanes(); 529 MakePlanes();
212 } 530 }
213 531
214 ////////////////////////////////////////////// << 532 ///////////////////////////////////////////////////////////////////////
215 // 533 //
216 // Fake default constructor - sets only member 534 // Fake default constructor - sets only member data and allocates memory
217 // for usage restri 535 // for usage restricted to object persistency.
218 // 536 //
219 G4Trap::G4Trap( __void__& a ) 537 G4Trap::G4Trap( __void__& a )
220 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo << 538 : G4CSGSolid(a), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
221 fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.), <<
222 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.) 539 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
223 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.) 540 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
224 { 541 {
225 MakePlanes(); 542 MakePlanes();
226 } 543 }
227 544
228 ////////////////////////////////////////////// << 545 ////////////////////////////////////////////////////////////////////////
229 // 546 //
230 // Destructor 547 // Destructor
231 548
232 G4Trap::~G4Trap() = default; << 549 G4Trap::~G4Trap()
>> 550 {
>> 551 }
233 552
234 ////////////////////////////////////////////// 553 //////////////////////////////////////////////////////////////////////////
235 // 554 //
236 // Copy constructor 555 // Copy constructor
237 556
238 G4Trap::G4Trap(const G4Trap& rhs) 557 G4Trap::G4Trap(const G4Trap& rhs)
239 : G4CSGSolid(rhs), halfCarTolerance(rhs.half << 558 : G4CSGSolid(rhs), fDz(rhs.fDz),
240 fDz(rhs.fDz), fTthetaCphi(rhs.fTthetaCphi) << 559 fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
241 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.f 560 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
242 fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.f 561 fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
243 { 562 {
244 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs << 563 for (size_t i=0; i<4; ++i)
245 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs. << 564 {
246 fTrapType = rhs.fTrapType; << 565 fPlanes[i].a = rhs.fPlanes[i].a;
>> 566 fPlanes[i].b = rhs.fPlanes[i].b;
>> 567 fPlanes[i].c = rhs.fPlanes[i].c;
>> 568 fPlanes[i].d = rhs.fPlanes[i].d;
>> 569 }
247 } 570 }
248 571
249 ////////////////////////////////////////////// 572 //////////////////////////////////////////////////////////////////////////
250 // 573 //
251 // Assignment operator 574 // Assignment operator
252 575
253 G4Trap& G4Trap::operator = (const G4Trap& rhs) << 576 G4Trap& G4Trap::operator = (const G4Trap& rhs)
254 { 577 {
255 // Check assignment to self 578 // Check assignment to self
256 // 579 //
257 if (this == &rhs) { return *this; } 580 if (this == &rhs) { return *this; }
258 581
259 // Copy base class data 582 // Copy base class data
260 // 583 //
261 G4CSGSolid::operator=(rhs); 584 G4CSGSolid::operator=(rhs);
262 585
263 // Copy data 586 // Copy data
264 // 587 //
265 halfCarTolerance = rhs.halfCarTolerance; << 588 fDz = rhs.fDz;
266 fDz = rhs.fDz; fTthetaCphi = rhs.fTthetaCphi << 589 fTthetaCphi = rhs.fTthetaCphi; fTthetaSphi = rhs.fTthetaSphi;
267 fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs 590 fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1;
268 fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs 591 fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2;
269 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs << 592 for (size_t i=0; i<4; ++i)
270 for (G4int i=0; i<6; ++i) { fAreas[i] = rhs. << 593 {
271 fTrapType = rhs.fTrapType; << 594 fPlanes[i].a = rhs.fPlanes[i].a;
>> 595 fPlanes[i].b = rhs.fPlanes[i].b;
>> 596 fPlanes[i].c = rhs.fPlanes[i].c;
>> 597 fPlanes[i].d = rhs.fPlanes[i].d;
>> 598 }
>> 599
272 return *this; 600 return *this;
273 } 601 }
274 602
275 ////////////////////////////////////////////// << 603 ///////////////////////////////////////////////////////////////////////
276 // 604 //
277 // Set all parameters, as for constructor - ch 605 // Set all parameters, as for constructor - check and set half-widths
278 // as well as angles: final check of coplanari 606 // as well as angles: final check of coplanarity
279 607
280 void G4Trap::SetAllParameters ( G4double pDz, 608 void G4Trap::SetAllParameters ( G4double pDz,
281 G4double pThet 609 G4double pTheta,
282 G4double pPhi, 610 G4double pPhi,
283 G4double pDy1, 611 G4double pDy1,
284 G4double pDx1, 612 G4double pDx1,
285 G4double pDx2, 613 G4double pDx2,
286 G4double pAlp1 614 G4double pAlp1,
287 G4double pDy2, 615 G4double pDy2,
288 G4double pDx3, 616 G4double pDx3,
289 G4double pDx4, 617 G4double pDx4,
290 G4double pAlp2 618 G4double pAlp2 )
291 { 619 {
292 // Reset data of the base class << 620 if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 || pDx3<=0 || pDx4<=0 )
293 fCubicVolume = 0; << 621 {
294 fSurfaceArea = 0; << 622 std::ostringstream message;
>> 623 message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
>> 624 << " X - "
>> 625 << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
>> 626 << " Y - " << pDy1 << ", " << pDy2 << G4endl
>> 627 << " Z - " << pDz;
>> 628 G4Exception("G4Trap::SetAllParameters()", "GeomSolids0002",
>> 629 FatalException, message);
>> 630 }
>> 631 fCubicVolume= 0.;
>> 632 fSurfaceArea= 0.;
295 fRebuildPolyhedron = true; 633 fRebuildPolyhedron = true;
>> 634 fDz=pDz;
>> 635 fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
>> 636 fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
>> 637
>> 638 fDy1=pDy1;
>> 639 fDx1=pDx1;
>> 640 fDx2=pDx2;
>> 641 fTalpha1=std::tan(pAlp1);
>> 642
>> 643 fDy2=pDy2;
>> 644 fDx3=pDx3;
>> 645 fDx4=pDx4;
>> 646 fTalpha2=std::tan(pAlp2);
296 647
297 // Set parameters <<
298 fDz = pDz; <<
299 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi <<
300 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi <<
301 <<
302 fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalp <<
303 fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalp <<
304 <<
305 CheckParameters(); <<
306 MakePlanes(); 648 MakePlanes();
307 } 649 }
308 650
309 ////////////////////////////////////////////// 651 //////////////////////////////////////////////////////////////////////////
310 // 652 //
311 // Check length parameters << 653 // Checking of coplanarity
312 654
313 void G4Trap::CheckParameters() << 655 G4bool G4Trap::MakePlanes()
314 { 656 {
315 if (fDz<=0 || << 657 G4bool good = true;
316 fDy1<=0 || fDx1<=0 || fDx2<=0 || << 658
317 fDy2<=0 || fDx3<=0 || fDx4<=0) << 659 G4ThreeVector pt[8] ;
>> 660
>> 661 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
>> 662 -fDz*fTthetaSphi-fDy1,-fDz);
>> 663 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
>> 664 -fDz*fTthetaSphi-fDy1,-fDz);
>> 665 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
>> 666 -fDz*fTthetaSphi+fDy1,-fDz);
>> 667 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
>> 668 -fDz*fTthetaSphi+fDy1,-fDz);
>> 669 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
>> 670 +fDz*fTthetaSphi-fDy2,+fDz);
>> 671 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
>> 672 +fDz*fTthetaSphi-fDy2,+fDz);
>> 673 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
>> 674 +fDz*fTthetaSphi+fDy2,+fDz);
>> 675 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
>> 676 +fDz*fTthetaSphi+fDy2,+fDz);
>> 677
>> 678 // Bottom side with normal approx. -Y
>> 679 //
>> 680 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]) ;
>> 681 if (!good)
318 { 682 {
319 std::ostringstream message; 683 std::ostringstream message;
320 message << "Invalid Length Parameters for << 684 message << "Face at ~-Y not planar for Solid: " << GetName();
321 << "\n X - " <<fDx1<<", "<<fDx2<< << 685 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
322 << "\n Y - " <<fDy1<<", "<<fDy2 <<
323 << "\n Z - " <<fDz; <<
324 G4Exception("G4Trap::CheckParameters()", " <<
325 FatalException, message); 686 FatalException, message);
326 } 687 }
327 } <<
328 <<
329 ////////////////////////////////////////////// <<
330 // <<
331 // Compute vertices and set side planes <<
332 <<
333 void G4Trap::MakePlanes() <<
334 { <<
335 G4double DzTthetaCphi = fDz*fTthetaCphi; <<
336 G4double DzTthetaSphi = fDz*fTthetaSphi; <<
337 G4double Dy1Talpha1 = fDy1*fTalpha1; <<
338 G4double Dy2Talpha2 = fDy2*fTalpha2; <<
339 <<
340 G4ThreeVector pt[8] = <<
341 { <<
342 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx <<
343 G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx <<
344 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx <<
345 G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx <<
346 G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx <<
347 G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx <<
348 G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx <<
349 G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx <<
350 }; <<
351 688
352 MakePlanes(pt); << 689 // Top side with normal approx. +Y
353 } << 690 //
354 << 691 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
355 ////////////////////////////////////////////// << 692 if (!good)
356 // <<
357 // Set side planes, check planarity <<
358 <<
359 void G4Trap::MakePlanes(const G4ThreeVector pt <<
360 { <<
361 constexpr G4int iface[4][4] = { {0,4,5,1}, { <<
362 const static G4String side[4] = { "~-Y", "~+ <<
363 <<
364 for (G4int i=0; i<4; ++i) <<
365 { 693 {
366 if (MakePlane(pt[iface[i][0]], << 694 std::ostringstream message;
367 pt[iface[i][1]], << 695 message << "Face at ~+Y not planar for Solid: " << GetName();
368 pt[iface[i][2]], << 696 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
369 pt[iface[i][3]], << 697 FatalException, message);
370 fPlanes[i])) continue; << 698 }
371 699
372 // Non planar side face << 700 // Front side with normal approx. -X
373 G4ThreeVector normal(fPlanes[i].a,fPlanes[ << 701 //
374 G4double dmax = 0; << 702 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
375 for (G4int k=0; k<4; ++k) << 703 if (!good)
376 { << 704 {
377 G4double dist = normal.dot(pt[iface[i][k <<
378 if (std::abs(dist) > std::abs(dmax)) dma <<
379 } <<
380 std::ostringstream message; 705 std::ostringstream message;
381 message << "Side face " << side[i] << " is << 706 message << "Face at ~-X not planar for Solid: " << GetName();
382 << GetName() << "\nDiscrepancy: " << 707 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
383 StreamInfo(message); << 708 FatalException, message);
>> 709 }
>> 710
>> 711 // Back side iwth normal approx. +X
>> 712 //
>> 713 good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
>> 714 if ( !good )
>> 715 {
>> 716 std::ostringstream message;
>> 717 message << "Face at ~+X not planar for Solid: " << GetName();
384 G4Exception("G4Trap::MakePlanes()", "GeomS 718 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
385 FatalException, message); 719 FatalException, message);
386 } 720 }
387 721
388 // Re-compute parameters << 722 return good;
389 SetCachedValues(); <<
390 } 723 }
391 724
392 ////////////////////////////////////////////// << 725 //////////////////////////////////////////////////////////////////////////////
393 // 726 //
394 // Calculate the coef's of the plane p1->p2->p 727 // Calculate the coef's of the plane p1->p2->p3->p4->p1
395 // where the ThreeVectors 1-4 are in anti-cloc << 728 // where the ThreeVectors 1-4 are in anti-clockwise order when viewed from
396 // from infront of the plane (i.e. from normal << 729 // infront of the plane (i.e. from normal direction).
397 // 730 //
398 // Return true if the points are coplanar, fal << 731 // Return true if the ThreeVectors are coplanar + set coef;s
>> 732 // false if ThreeVectors are not coplanar
399 733
400 G4bool G4Trap::MakePlane( const G4ThreeVector& 734 G4bool G4Trap::MakePlane( const G4ThreeVector& p1,
401 const G4ThreeVector& 735 const G4ThreeVector& p2,
402 const G4ThreeVector& 736 const G4ThreeVector& p3,
403 const G4ThreeVector& 737 const G4ThreeVector& p4,
404 TrapSidePlane& 738 TrapSidePlane& plane )
405 { 739 {
406 G4ThreeVector normal = ((p4 - p2).cross(p3 - << 740 G4double a, b, c, sd;
407 if (std::abs(normal.x()) < DBL_EPSILON) norm << 741 G4ThreeVector v12, v13, v14, Vcross;
408 if (std::abs(normal.y()) < DBL_EPSILON) norm <<
409 if (std::abs(normal.z()) < DBL_EPSILON) norm <<
410 normal = normal.unit(); <<
411 <<
412 G4ThreeVector centre = (p1 + p2 + p3 + p4)*0 <<
413 plane.a = normal.x(); <<
414 plane.b = normal.y(); <<
415 plane.c = normal.z(); <<
416 plane.d = -normal.dot(centre); <<
417 <<
418 // compute distances and check planarity <<
419 G4double d1 = std::abs(normal.dot(p1) + plan <<
420 G4double d2 = std::abs(normal.dot(p2) + plan <<
421 G4double d3 = std::abs(normal.dot(p3) + plan <<
422 G4double d4 = std::abs(normal.dot(p4) + plan <<
423 G4double dmax = std::max(std::max(std::max(d <<
424 742
425 return dmax <= 1000 * kCarTolerance; << 743 G4bool good;
426 } <<
427 <<
428 ////////////////////////////////////////////// <<
429 // <<
430 // Recompute parameters using planes <<
431 744
432 void G4Trap::SetCachedValues() << 745 v12 = p2 - p1;
433 { << 746 v13 = p3 - p1;
434 // Set indeces << 747 v14 = p4 - p1;
435 constexpr G4int iface[6][4] = << 748 Vcross = v12.cross(v13);
436 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2, <<
437 <<
438 // Get vertices <<
439 G4ThreeVector pt[8]; <<
440 GetVertices(pt); <<
441 749
442 // Set face areas << 750 if (std::fabs(Vcross.dot(v14)/(Vcross.mag()*v14.mag())) > kCoplanar_Tolerance)
443 for (G4int i=0; i<6; ++i) <<
444 { 751 {
445 fAreas[i] = G4GeomTools::QuadAreaNormal(pt << 752 good = false;
446 pt <<
447 pt <<
448 pt <<
449 } <<
450 for (G4int i=1; i<6; ++i) { fAreas[i] += fAr <<
451 <<
452 // Define type of trapezoid <<
453 fTrapType = 0; <<
454 if (fPlanes[0].b == -1 && fPlanes[1].b == 1 <<
455 std::abs(fPlanes[0].a) < DBL_EPSILON && <<
456 std::abs(fPlanes[0].c) < DBL_EPSILON && <<
457 std::abs(fPlanes[1].a) < DBL_EPSILON && <<
458 std::abs(fPlanes[1].c) < DBL_EPSILON) <<
459 { <<
460 fTrapType = 1; // YZ section is a rectangl <<
461 if (std::abs(fPlanes[2].a + fPlanes[3].a) <<
462 std::abs(fPlanes[2].c - fPlanes[3].c) <<
463 fPlanes[2].b == 0 && <<
464 fPlanes[3].b == 0) <<
465 { <<
466 fTrapType = 2; // ... and XZ section is <<
467 fPlanes[2].a = -fPlanes[3].a; <<
468 fPlanes[2].c = fPlanes[3].c; <<
469 } <<
470 if (std::abs(fPlanes[2].a + fPlanes[3].a) <<
471 std::abs(fPlanes[2].b - fPlanes[3].b) <<
472 fPlanes[2].c == 0 && <<
473 fPlanes[3].c == 0) <<
474 { <<
475 fTrapType = 3; // ... and XY section is <<
476 fPlanes[2].a = -fPlanes[3].a; <<
477 fPlanes[2].b = fPlanes[3].b; <<
478 } <<
479 } 753 }
480 } << 754 else
481 <<
482 ////////////////////////////////////////////// <<
483 // <<
484 // Get volume <<
485 <<
486 G4double G4Trap::GetCubicVolume() <<
487 { <<
488 if (fCubicVolume == 0) <<
489 { <<
490 G4ThreeVector pt[8]; <<
491 GetVertices(pt); <<
492 <<
493 G4double dz = pt[4].z() - pt[0].z(); <<
494 G4double dy1 = pt[2].y() - pt[0].y(); <<
495 G4double dx1 = pt[1].x() - pt[0].x(); <<
496 G4double dx2 = pt[3].x() - pt[2].x(); <<
497 G4double dy2 = pt[6].y() - pt[4].y(); <<
498 G4double dx3 = pt[5].x() - pt[4].x(); <<
499 G4double dx4 = pt[7].x() - pt[6].x(); <<
500 <<
501 fCubicVolume = ((dx1 + dx2 + dx3 + dx4)*(d <<
502 (dx4 + dx3 - dx2 - dx1)*(d <<
503 } <<
504 return fCubicVolume; <<
505 } <<
506 <<
507 ////////////////////////////////////////////// <<
508 // <<
509 // Get surface area <<
510 <<
511 G4double G4Trap::GetSurfaceArea() <<
512 { <<
513 if (fSurfaceArea == 0) <<
514 { 755 {
515 G4ThreeVector pt[8]; << 756 // a,b,c correspond to the x/y/z components of the
516 G4int iface [6][4] = << 757 // normal vector to the plane
517 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2, << 758
518 << 759 // a = (p2.y()-p1.y())*(p1.z()+p2.z())+(p3.y()-p2.y())*(p2.z()+p3.z());
519 GetVertices(pt); << 760 // a += (p4.y()-p3.y())*(p3.z()+p4.z())+(p1.y()-p4.y())*(p4.z()+p1.z()); // ?
520 for (const auto & i : iface) << 761 // b = (p2.z()-p1.z())*(p1.x()+p2.x())+(p3.z()-p2.z())*(p2.x()+p3.x());
>> 762 // b += (p4.z()-p3.z())*(p3.x()+p4.x())+(p1.z()-p4.z())*(p4.x()+p1.x()); // ?
>> 763 // c = (p2.x()-p1.x())*(p1.y()+p2.y())+(p3.x()-p2.x())*(p2.y()+p3.y());
>> 764 // c += (p4.x()-p3.x())*(p3.y()+p4.y())+(p1.x()-p4.x())*(p4.y()+p1.y()); // ?
>> 765
>> 766 // Let create diagonals 4-2 and 3-1 than (4-2)x(3-1) provides
>> 767 // vector perpendicular to the plane directed to outside !!!
>> 768 // and a,b,c, = f(1,2,3,4) external relative to trap normal
>> 769
>> 770 a = +(p4.y() - p2.y())*(p3.z() - p1.z())
>> 771 - (p3.y() - p1.y())*(p4.z() - p2.z());
>> 772
>> 773 b = -(p4.x() - p2.x())*(p3.z() - p1.z())
>> 774 + (p3.x() - p1.x())*(p4.z() - p2.z());
>> 775
>> 776 c = +(p4.x() - p2.x())*(p3.y() - p1.y())
>> 777 - (p3.x() - p1.x())*(p4.y() - p2.y());
>> 778
>> 779 sd = std::sqrt( a*a + b*b + c*c ); // so now vector plane.(a,b,c) is unit
>> 780
>> 781 if( sd > 0 )
>> 782 {
>> 783 plane.a = a/sd;
>> 784 plane.b = b/sd;
>> 785 plane.c = c/sd;
>> 786 }
>> 787 else
521 { 788 {
522 fSurfaceArea += G4GeomTools::QuadAreaNor << 789 std::ostringstream message;
523 << 790 message << "Invalid parameters: norm.mod() <= 0, for Solid: "
524 << 791 << GetName();
525 << 792 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
>> 793 FatalException, message) ;
526 } 794 }
>> 795 // Calculate D: p1 in in plane so D=-n.p1.Vect()
>> 796
>> 797 plane.d = -( plane.a*p1.x() + plane.b*p1.y() + plane.c*p1.z() );
>> 798
>> 799 good = true;
527 } 800 }
528 return fSurfaceArea; << 801 return good;
529 } 802 }
530 803
531 ////////////////////////////////////////////// << 804 //////////////////////////////////////////////////////////////////////////////
532 // 805 //
533 // Dispatch to parameterisation for replicatio 806 // Dispatch to parameterisation for replication mechanism dimension
534 // computation & modification. 807 // computation & modification.
535 808
536 void G4Trap::ComputeDimensions( G4VPVPar 809 void G4Trap::ComputeDimensions( G4VPVParameterisation* p,
537 const G4int n, 810 const G4int n,
538 const G4VPhysi 811 const G4VPhysicalVolume* pRep )
539 { 812 {
540 p->ComputeDimensions(*this,n,pRep); 813 p->ComputeDimensions(*this,n,pRep);
541 } 814 }
542 815
543 ////////////////////////////////////////////// << 816 ////////////////////////////////////////////////////////////////////////
544 // 817 //
545 // Get bounding box 818 // Get bounding box
546 819
547 void G4Trap::BoundingLimits(G4ThreeVector& pMi << 820 void G4Trap::Extent(G4ThreeVector& pMin, G4ThreeVector& pMax) const
548 { 821 {
549 G4ThreeVector pt[8]; << 822 G4double dz = GetZHalfLength();
550 GetVertices(pt); << 823 G4double dx1 = GetXHalfLength1();
551 << 824 G4double dx2 = GetXHalfLength2();
552 G4double xmin = kInfinity, xmax = -kInfinity << 825 G4double dx3 = GetXHalfLength3();
553 G4double ymin = kInfinity, ymax = -kInfinity << 826 G4double dx4 = GetXHalfLength4();
554 for (const auto & i : pt) << 827 G4double dy1 = GetYHalfLength1();
555 { << 828 G4double dy2 = GetYHalfLength2();
556 G4double x = i.x(); << 829
557 if (x < xmin) xmin = x; << 830 G4double x0 = dz*fTthetaCphi;
558 if (x > xmax) xmax = x; << 831 G4double x1 = dy1*GetTanAlpha1();
559 G4double y = i.y(); << 832 G4double x2 = dy2*GetTanAlpha2();
560 if (y < ymin) ymin = y; << 833 G4double xmin =
561 if (y > ymax) ymax = y; << 834 std::min(
562 } << 835 std::min(
>> 836 std::min(-x0-x1-dx1,-x0+x1-dx2),x0-x2-dx3),x0+x2-dx4);
>> 837 G4double xmax =
>> 838 std::max(
>> 839 std::max(
>> 840 std::max(-x0-x1+dx1,-x0+x1+dx2),x0-x2+dx3),x0+x2+dx4);
>> 841
>> 842 G4double y0 = dz*fTthetaSphi;
>> 843 G4double ymin = std::min(-y0-dy1,y0-dy2);
>> 844 G4double ymax = std::max(-y0+dy1,y0+dy2);
563 845
564 G4double dz = GetZHalfLength(); <<
565 pMin.set(xmin,ymin,-dz); 846 pMin.set(xmin,ymin,-dz);
566 pMax.set(xmax,ymax, dz); 847 pMax.set(xmax,ymax, dz);
567 848
568 // Check correctness of the bounding box 849 // Check correctness of the bounding box
569 // 850 //
570 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 851 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
571 { 852 {
572 std::ostringstream message; 853 std::ostringstream message;
573 message << "Bad bounding box (min >= max) 854 message << "Bad bounding box (min >= max) for solid: "
574 << GetName() << " !" 855 << GetName() << " !"
575 << "\npMin = " << pMin 856 << "\npMin = " << pMin
576 << "\npMax = " << pMax; 857 << "\npMax = " << pMax;
577 G4Exception("G4Trap::BoundingLimits()", "G << 858 G4Exception("G4Trap::Extent()", "GeomMgt0001", JustWarning, message);
578 JustWarning, message); <<
579 DumpInfo(); 859 DumpInfo();
580 } 860 }
581 } 861 }
582 862
583 ////////////////////////////////////////////// << 863 ////////////////////////////////////////////////////////////////////////
584 // 864 //
585 // Calculate extent under transform and specif 865 // Calculate extent under transform and specified limit
586 866
587 G4bool G4Trap::CalculateExtent( const EAxis pA 867 G4bool G4Trap::CalculateExtent( const EAxis pAxis,
588 const G4VoxelL 868 const G4VoxelLimits& pVoxelLimit,
589 const G4Affine 869 const G4AffineTransform& pTransform,
590 G4double 870 G4double& pMin, G4double& pMax) const
591 { 871 {
592 G4ThreeVector bmin, bmax; 872 G4ThreeVector bmin, bmax;
593 G4bool exist; 873 G4bool exist;
594 874
595 // Check bounding box (bbox) 875 // Check bounding box (bbox)
596 // 876 //
597 BoundingLimits(bmin,bmax); << 877 Extent(bmin,bmax);
598 G4BoundingEnvelope bbox(bmin,bmax); 878 G4BoundingEnvelope bbox(bmin,bmax);
599 #ifdef G4BBOX_EXTENT 879 #ifdef G4BBOX_EXTENT
600 return bbox.CalculateExtent(pAxis,pVoxelLimi << 880 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
601 #endif 881 #endif
602 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox 882 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
603 { 883 {
604 return exist = pMin < pMax; << 884 return exist = (pMin < pMax) ? true : false;
605 } 885 }
606 886
607 // Set bounding envelope (benv) and calculat 887 // Set bounding envelope (benv) and calculate extent
608 // 888 //
609 G4ThreeVector pt[8]; << 889 G4double dz = GetZHalfLength();
610 GetVertices(pt); << 890 G4double dx1 = GetXHalfLength1();
>> 891 G4double dx2 = GetXHalfLength2();
>> 892 G4double dx3 = GetXHalfLength3();
>> 893 G4double dx4 = GetXHalfLength4();
>> 894 G4double dy1 = GetYHalfLength1();
>> 895 G4double dy2 = GetYHalfLength2();
>> 896
>> 897 G4double x0 = dz*fTthetaCphi;
>> 898 G4double x1 = dy1*GetTanAlpha1();
>> 899 G4double x2 = dy2*GetTanAlpha2();
>> 900 G4double y0 = dz*fTthetaSphi;
611 901
612 G4ThreeVectorList baseA(4), baseB(4); 902 G4ThreeVectorList baseA(4), baseB(4);
613 baseA[0] = pt[0]; << 903 baseA[0].set(-x0-x1-dx1,-y0-dy1,-dz);
614 baseA[1] = pt[1]; << 904 baseA[1].set(-x0-x1+dx1,-y0-dy1,-dz);
615 baseA[2] = pt[3]; << 905 baseA[2].set(-x0+x1+dx2,-y0+dy1,-dz);
616 baseA[3] = pt[2]; << 906 baseA[3].set(-x0+x1-dx2,-y0+dy1,-dz);
617 << 907
618 baseB[0] = pt[4]; << 908 baseB[0].set( x0-x2-dx3, y0-dy2, dz);
619 baseB[1] = pt[5]; << 909 baseB[1].set( x0-x2+dx3, y0-dy2, dz);
620 baseB[2] = pt[7]; << 910 baseB[2].set( x0+x2+dx4, y0+dy2, dz);
621 baseB[3] = pt[6]; << 911 baseB[3].set( x0+x2-dx4, y0+dy2, dz);
622 912
623 std::vector<const G4ThreeVectorList *> polyg 913 std::vector<const G4ThreeVectorList *> polygons(2);
624 polygons[0] = &baseA; 914 polygons[0] = &baseA;
625 polygons[1] = &baseB; 915 polygons[1] = &baseB;
626 916
627 G4BoundingEnvelope benv(bmin,bmax,polygons); 917 G4BoundingEnvelope benv(bmin,bmax,polygons);
628 exist = benv.CalculateExtent(pAxis,pVoxelLim 918 exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
629 return exist; 919 return exist;
630 } 920 }
631 921
632 ////////////////////////////////////////////// << 922 ////////////////////////////////////////////////////////////////////////
633 // 923 //
634 // Return whether point is inside/outside/on_s << 924 // Return whether point inside/outside/on surface, using tolerance
635 925
636 EInside G4Trap::Inside( const G4ThreeVector& p 926 EInside G4Trap::Inside( const G4ThreeVector& p ) const
637 { 927 {
638 switch (fTrapType) << 928 EInside in;
>> 929 G4double Dist;
>> 930 G4int i;
>> 931 if ( std::fabs(p.z()) <= fDz-kCarTolerance*0.5)
>> 932 {
>> 933 in = kInside;
>> 934
>> 935 for ( i = 0;i < 4;i++ )
>> 936 {
>> 937 Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
>> 938 +fPlanes[i].c*p.z() + fPlanes[i].d;
>> 939
>> 940 if (Dist > kCarTolerance*0.5) return in = kOutside;
>> 941 else if (Dist > -kCarTolerance*0.5) in = kSurface;
>> 942
>> 943 }
>> 944 }
>> 945 else if (std::fabs(p.z()) <= fDz+kCarTolerance*0.5)
639 { 946 {
640 case 0: // General case << 947 in = kSurface;
>> 948
>> 949 for ( i = 0; i < 4; i++ )
641 { 950 {
642 G4double dz = std::abs(p.z())-fDz; << 951 Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
643 G4double dy1 = fPlanes[0].b*p.y()+fPlane << 952 +fPlanes[i].c*p.z() + fPlanes[i].d;
644 G4double dy2 = fPlanes[1].b*p.y()+fPlane <<
645 G4double dy = std::max(dz,std::max(dy1,d <<
646 <<
647 G4double dx1 = fPlanes[2].a*p.x()+fPlane <<
648 + fPlanes[2].c*p.z()+fPlane <<
649 G4double dx2 = fPlanes[3].a*p.x()+fPlane <<
650 + fPlanes[3].c*p.z()+fPlane <<
651 G4double dist = std::max(dy,std::max(dx1 <<
652 <<
653 return (dist > halfCarTolerance) ? kOuts <<
654 ((dist > -halfCarTolerance) ? kSurface <<
655 } <<
656 case 1: // YZ section is a rectangle <<
657 { <<
658 G4double dz = std::abs(p.z())-fDz; <<
659 G4double dy = std::max(dz,std::abs(p.y() <<
660 G4double dx1 = fPlanes[2].a*p.x()+fPlane <<
661 + fPlanes[2].c*p.z()+fPlane <<
662 G4double dx2 = fPlanes[3].a*p.x()+fPlane <<
663 + fPlanes[3].c*p.z()+fPlane <<
664 G4double dist = std::max(dy,std::max(dx1 <<
665 <<
666 return (dist > halfCarTolerance) ? kOuts <<
667 ((dist > -halfCarTolerance) ? kSurface <<
668 } <<
669 case 2: // YZ section is a rectangle and <<
670 { // XZ section is an isosceles trap <<
671 G4double dz = std::abs(p.z())-fDz; <<
672 G4double dy = std::max(dz,std::abs(p.y() <<
673 G4double dx = fPlanes[3].a*std::abs(p.x( <<
674 + fPlanes[3].c*p.z()+fPlanes <<
675 G4double dist = std::max(dy,dx); <<
676 <<
677 return (dist > halfCarTolerance) ? kOuts <<
678 ((dist > -halfCarTolerance) ? kSurface <<
679 } <<
680 case 3: // YZ section is a rectangle and <<
681 { // XY section is an isosceles trap <<
682 G4double dz = std::abs(p.z())-fDz; <<
683 G4double dy = std::max(dz,std::abs(p.y() <<
684 G4double dx = fPlanes[3].a*std::abs(p.x( <<
685 + fPlanes[3].b*p.y()+fPlanes <<
686 G4double dist = std::max(dy,dx); <<
687 953
688 return (dist > halfCarTolerance) ? kOuts << 954 if (Dist > kCarTolerance*0.5) return in = kOutside;
689 ((dist > -halfCarTolerance) ? kSurface <<
690 } 955 }
691 } 956 }
692 return kOutside; << 957 else in = kOutside;
>> 958
>> 959 return in;
693 } 960 }
694 961
695 ////////////////////////////////////////////// << 962 /////////////////////////////////////////////////////////////////////////////
696 // 963 //
697 // Determine side, and return corresponding no << 964 // Calculate side nearest to p, and return normal
>> 965 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
698 966
699 G4ThreeVector G4Trap::SurfaceNormal( const G4T 967 G4ThreeVector G4Trap::SurfaceNormal( const G4ThreeVector& p ) const
700 { 968 {
701 G4double nx = 0, ny = 0, nz = 0; << 969 G4int i, noSurfaces = 0;
702 G4double dz = std::abs(p.z()) - fDz; << 970 G4double dist, distz, distx, disty, distmx, distmy, safe = kInfinity;
703 nz = std::copysign(G4double(std::abs(dz) <= << 971 G4double delta = 0.5*kCarTolerance;
>> 972 G4ThreeVector norm, sumnorm(0.,0.,0.);
704 973
705 switch (fTrapType) << 974 for (i = 0; i < 4; i++)
706 { 975 {
707 case 0: // General case << 976 dist = std::fabs(fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
708 { << 977 + fPlanes[i].c*p.z() + fPlanes[i].d);
709 for (G4int i=0; i<2; ++i) << 978 if ( dist < safe )
710 { <<
711 G4double dy = fPlanes[i].b*p.y() + fPl <<
712 if (std::abs(dy) > halfCarTolerance) c <<
713 ny = fPlanes[i].b; <<
714 nz += fPlanes[i].c; <<
715 break; <<
716 } <<
717 for (G4int i=2; i<4; ++i) <<
718 { <<
719 G4double dx = fPlanes[i].a*p.x() + <<
720 fPlanes[i].b*p.y() + fPl <<
721 if (std::abs(dx) > halfCarTolerance) c <<
722 nx = fPlanes[i].a; <<
723 ny += fPlanes[i].b; <<
724 nz += fPlanes[i].c; <<
725 break; <<
726 } <<
727 break; <<
728 } <<
729 case 1: // YZ section - rectangle <<
730 { <<
731 G4double dy = std::abs(p.y()) + fPlanes[ <<
732 ny = std::copysign(G4double(std::abs(dy) <<
733 for (G4int i=2; i<4; ++i) <<
734 { <<
735 G4double dx = fPlanes[i].a*p.x() + <<
736 fPlanes[i].b*p.y() + fPl <<
737 if (std::abs(dx) > halfCarTolerance) c <<
738 nx = fPlanes[i].a; <<
739 ny += fPlanes[i].b; <<
740 nz += fPlanes[i].c; <<
741 break; <<
742 } <<
743 break; <<
744 } <<
745 case 2: // YZ section - rectangle, XZ sect <<
746 { 979 {
747 G4double dy = std::abs(p.y()) + fPlanes[ << 980 safe = dist;
748 ny = std::copysign(G4double(std::abs(dy) <<
749 G4double dx = fPlanes[3].a*std::abs(p.x( <<
750 fPlanes[3].c*p.z() + fPlan <<
751 G4double k = std::abs(dx) <= halfCarTole <<
752 nx = std::copysign(k, p.x())*fPlanes[3] <<
753 nz += k*fPlanes[3].c; <<
754 break; <<
755 } <<
756 case 3: // YZ section - rectangle, XY sect <<
757 { <<
758 G4double dy = std::abs(p.y()) + fPlanes[ <<
759 ny = std::copysign(G4double(std::abs(dy) <<
760 G4double dx = fPlanes[3].a*std::abs(p.x( <<
761 fPlanes[3].b*p.y() + fPlan <<
762 G4double k = std::abs(dx) <= halfCarTole <<
763 nx = std::copysign(k, p.x())*fPlanes[3] <<
764 ny += k*fPlanes[3].b; <<
765 break; <<
766 } 981 }
767 } 982 }
>> 983 distz = std::fabs( std::fabs( p.z() ) - fDz );
768 984
769 // Return normal << 985 distmy = std::fabs( fPlanes[0].a*p.x() + fPlanes[0].b*p.y()
770 // << 986 + fPlanes[0].c*p.z() + fPlanes[0].d );
771 G4double mag2 = nx*nx + ny*ny + nz*nz; << 987
772 if (mag2 == 1) return { nx,ny,nz }; << 988 disty = std::fabs( fPlanes[1].a*p.x() + fPlanes[1].b*p.y()
773 else if (mag2 != 0) return G4ThreeVector(nx, << 989 + fPlanes[1].c*p.z() + fPlanes[1].d );
774 else << 990
>> 991 distmx = std::fabs( fPlanes[2].a*p.x() + fPlanes[2].b*p.y()
>> 992 + fPlanes[2].c*p.z() + fPlanes[2].d );
>> 993
>> 994 distx = std::fabs( fPlanes[3].a*p.x() + fPlanes[3].b*p.y()
>> 995 + fPlanes[3].c*p.z() + fPlanes[3].d );
>> 996
>> 997 G4ThreeVector nX = G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
>> 998 G4ThreeVector nmX = G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
>> 999 G4ThreeVector nY = G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
>> 1000 G4ThreeVector nmY = G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
>> 1001 G4ThreeVector nZ = G4ThreeVector(0.,0.,1.0);
>> 1002
>> 1003 if (distx <= delta)
>> 1004 {
>> 1005 noSurfaces ++;
>> 1006 sumnorm += nX;
>> 1007 }
>> 1008 if (distmx <= delta)
>> 1009 {
>> 1010 noSurfaces ++;
>> 1011 sumnorm += nmX;
>> 1012 }
>> 1013 if (disty <= delta)
>> 1014 {
>> 1015 noSurfaces ++;
>> 1016 sumnorm += nY;
>> 1017 }
>> 1018 if (distmy <= delta)
>> 1019 {
>> 1020 noSurfaces ++;
>> 1021 sumnorm += nmY;
>> 1022 }
>> 1023 if (distz <= delta)
>> 1024 {
>> 1025 noSurfaces ++;
>> 1026 if ( p.z() >= 0.) sumnorm += nZ;
>> 1027 else sumnorm -= nZ;
>> 1028 }
>> 1029 if ( noSurfaces == 0 )
775 { 1030 {
776 // Point is not on the surface <<
777 // <<
778 #ifdef G4CSGDEBUG 1031 #ifdef G4CSGDEBUG
779 std::ostringstream message; <<
780 G4long oldprc = message.precision(16); <<
781 message << "Point p is not on surface (!?) <<
782 << GetName() << G4endl; <<
783 message << "Position:\n"; <<
784 message << " p.x() = " << p.x()/mm << " <<
785 message << " p.y() = " << p.y()/mm << " <<
786 message << " p.z() = " << p.z()/mm << " <<
787 G4cout.precision(oldprc) ; <<
788 G4Exception("G4Trap::SurfaceNormal(p)", "G 1032 G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002",
789 JustWarning, message ); << 1033 JustWarning, "Point p is not on surface !?" );
790 DumpInfo(); << 1034 #endif
791 #endif << 1035 norm = ApproxSurfaceNormal(p);
792 return ApproxSurfaceNormal(p); <<
793 } 1036 }
>> 1037 else if ( noSurfaces == 1 ) norm = sumnorm;
>> 1038 else norm = sumnorm.unit();
>> 1039 return norm;
794 } 1040 }
795 1041
796 ////////////////////////////////////////////// << 1042 ////////////////////////////////////////////////////////////////////////////////////
797 // 1043 //
798 // Algorithm for SurfaceNormal() following the 1044 // Algorithm for SurfaceNormal() following the original specification
799 // for points not on the surface 1045 // for points not on the surface
800 1046
801 G4ThreeVector G4Trap::ApproxSurfaceNormal( con 1047 G4ThreeVector G4Trap::ApproxSurfaceNormal( const G4ThreeVector& p ) const
802 { 1048 {
803 G4double dist = -DBL_MAX; << 1049 G4double safe=kInfinity,Dist,safez;
804 G4int iside = 0; << 1050 G4int i,imin=0;
805 for (G4int i=0; i<4; ++i) << 1051 for (i=0;i<4;i++)
806 { << 1052 {
807 G4double d = fPlanes[i].a*p.x() + << 1053 Dist=std::fabs(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
808 fPlanes[i].b*p.y() + << 1054 +fPlanes[i].c*p.z()+fPlanes[i].d);
809 fPlanes[i].c*p.z() + fPlanes[ << 1055 if (Dist<safe)
810 if (d > dist) { dist = d; iside = i; } << 1056 {
>> 1057 safe=Dist;
>> 1058 imin=i;
>> 1059 }
>> 1060 }
>> 1061 safez=std::fabs(std::fabs(p.z())-fDz);
>> 1062 if (safe<safez)
>> 1063 {
>> 1064 return G4ThreeVector(fPlanes[imin].a,fPlanes[imin].b,fPlanes[imin].c);
811 } 1065 }
812 <<
813 G4double distz = std::abs(p.z()) - fDz; <<
814 if (dist > distz) <<
815 return { fPlanes[iside].a, fPlanes[iside]. <<
816 else 1066 else
817 return { 0, 0, (G4double)((p.z() < 0) ? -1 << 1067 {
>> 1068 if (p.z()>0)
>> 1069 {
>> 1070 return G4ThreeVector(0,0,1);
>> 1071 }
>> 1072 else
>> 1073 {
>> 1074 return G4ThreeVector(0,0,-1);
>> 1075 }
>> 1076 }
818 } 1077 }
819 1078
820 ////////////////////////////////////////////// << 1079 ////////////////////////////////////////////////////////////////////////////
821 // 1080 //
822 // Calculate distance to shape from outside << 1081 // Calculate distance to shape from outside - return kInfinity if no intersection
823 // - return kInfinity if no intersection << 1082 //
>> 1083 // ALGORITHM:
>> 1084 // For each component, calculate pair of minimum and maximum intersection
>> 1085 // values for which the particle is in the extent of the shape
>> 1086 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
824 1087
825 G4double G4Trap::DistanceToIn(const G4ThreeVec << 1088 G4double G4Trap::DistanceToIn( const G4ThreeVector& p,
826 const G4ThreeVec << 1089 const G4ThreeVector& v ) const
827 { 1090 {
828 // Z intersections <<
829 // <<
830 if ((std::abs(p.z()) - fDz) >= -halfCarToler <<
831 return kInfinity; <<
832 G4double invz = (-v.z() == 0) ? DBL_MAX : -1 <<
833 G4double dz = (invz < 0) ? fDz : -fDz; <<
834 G4double tzmin = (p.z() + dz)*invz; <<
835 G4double tzmax = (p.z() - dz)*invz; <<
836 1091
837 // Y intersections << 1092 G4double snxt; // snxt = default return value
>> 1093 G4double max,smax,smin;
>> 1094 G4double pdist,Comp,vdist;
>> 1095 G4int i;
838 // 1096 //
839 G4double tymin = 0, tymax = DBL_MAX; << 1097 // Z Intersection range
840 G4int i = 0; << 1098 //
841 for ( ; i<2; ++i) << 1099 if ( v.z() > 0 )
842 { 1100 {
843 G4double cosa = fPlanes[i].b*v.y() + fPlan << 1101 max = fDz - p.z() ;
844 G4double dist = fPlanes[i].b*p.y() + fPlan << 1102 if (max > 0.5*kCarTolerance)
845 if (dist >= -halfCarTolerance) <<
846 { 1103 {
847 if (cosa >= 0) return kInfinity; << 1104 smax = max/v.z();
848 G4double tmp = -dist/cosa; << 1105 smin = (-fDz-p.z())/v.z();
849 if (tymin < tmp) tymin = tmp; <<
850 } 1106 }
851 else if (cosa > 0) << 1107 else
852 { 1108 {
853 G4double tmp = -dist/cosa; << 1109 return snxt=kInfinity;
854 if (tymax > tmp) tymax = tmp; <<
855 } 1110 }
856 } 1111 }
857 << 1112 else if (v.z() < 0 )
858 // Z intersections << 1113 {
859 // << 1114 max = - fDz - p.z() ;
860 G4double txmin = 0, txmax = DBL_MAX; << 1115 if (max < -0.5*kCarTolerance )
861 for ( ; i<4; ++i) << 1116 {
>> 1117 smax=max/v.z();
>> 1118 smin=(fDz-p.z())/v.z();
>> 1119 }
>> 1120 else
>> 1121 {
>> 1122 return snxt=kInfinity;
>> 1123 }
>> 1124 }
>> 1125 else
862 { 1126 {
863 G4double cosa = fPlanes[i].a*v.x()+fPlanes << 1127 if (std::fabs(p.z())<fDz - 0.5*kCarTolerance) // Inside was <=fDz
864 G4double dist = fPlanes[i].a*p.x()+fPlanes <<
865 fPlanes[i].d; <<
866 if (dist >= -halfCarTolerance) <<
867 { 1128 {
868 if (cosa >= 0) return kInfinity; << 1129 smin=0;
869 G4double tmp = -dist/cosa; << 1130 smax=kInfinity;
870 if (txmin < tmp) txmin = tmp; <<
871 } 1131 }
872 else if (cosa > 0) << 1132 else
873 { 1133 {
874 G4double tmp = -dist/cosa; << 1134 return snxt=kInfinity;
875 if (txmax > tmp) txmax = tmp; <<
876 } 1135 }
877 } 1136 }
878 1137
879 // Find distance << 1138 for (i=0;i<4;i++)
>> 1139 {
>> 1140 pdist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
>> 1141 +fPlanes[i].c*p.z()+fPlanes[i].d;
>> 1142 Comp=fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
>> 1143 if ( pdist >= -0.5*kCarTolerance ) // was >0
>> 1144 {
>> 1145 //
>> 1146 // Outside the plane -> this is an extent entry distance
>> 1147 //
>> 1148 if (Comp >= 0) // was >0
>> 1149 {
>> 1150 return snxt=kInfinity ;
>> 1151 }
>> 1152 else
>> 1153 {
>> 1154 vdist=-pdist/Comp;
>> 1155 if (vdist>smin)
>> 1156 {
>> 1157 if (vdist<smax)
>> 1158 {
>> 1159 smin = vdist;
>> 1160 }
>> 1161 else
>> 1162 {
>> 1163 return snxt=kInfinity;
>> 1164 }
>> 1165 }
>> 1166 }
>> 1167 }
>> 1168 else
>> 1169 {
>> 1170 //
>> 1171 // Inside the plane -> couble be an extent exit distance (smax)
>> 1172 //
>> 1173 if (Comp>0) // Will leave extent
>> 1174 {
>> 1175 vdist=-pdist/Comp;
>> 1176 if (vdist<smax)
>> 1177 {
>> 1178 if (vdist>smin)
>> 1179 {
>> 1180 smax=vdist;
>> 1181 }
>> 1182 else
>> 1183 {
>> 1184 return snxt=kInfinity;
>> 1185 }
>> 1186 }
>> 1187 }
>> 1188 }
>> 1189 }
880 // 1190 //
881 G4double tmin = std::max(std::max(txmin,tymi << 1191 // Checks in non z plane intersections ensure smin<smax
882 G4double tmax = std::min(std::min(txmax,tyma << 1192 //
883 << 1193 if (smin >=0 )
884 if (tmax <= tmin + halfCarTolerance) return << 1194 {
885 return (tmin < halfCarTolerance ) ? 0. : tmi << 1195 snxt = smin ;
>> 1196 }
>> 1197 else
>> 1198 {
>> 1199 snxt = 0 ;
>> 1200 }
>> 1201 return snxt;
886 } 1202 }
887 1203
888 ////////////////////////////////////////////// << 1204 ///////////////////////////////////////////////////////////////////////////
889 // 1205 //
890 // Calculate exact shortest distance to any bo 1206 // Calculate exact shortest distance to any boundary from outside
891 // This is the best fast estimation of the sho 1207 // This is the best fast estimation of the shortest distance to trap
892 // - return 0 if point is inside << 1208 // - Returns 0 is ThreeVector inside
893 1209
894 G4double G4Trap::DistanceToIn( const G4ThreeVe 1210 G4double G4Trap::DistanceToIn( const G4ThreeVector& p ) const
895 { 1211 {
896 switch (fTrapType) << 1212 G4double safe=0.0,Dist;
>> 1213 G4int i;
>> 1214 safe=std::fabs(p.z())-fDz;
>> 1215 for (i=0;i<4;i++)
897 { 1216 {
898 case 0: // General case << 1217 Dist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
899 { << 1218 +fPlanes[i].c*p.z()+fPlanes[i].d;
900 G4double dz = std::abs(p.z())-fDz; << 1219 if (Dist > safe) safe=Dist;
901 G4double dy1 = fPlanes[0].b*p.y()+fPlane <<
902 G4double dy2 = fPlanes[1].b*p.y()+fPlane <<
903 G4double dy = std::max(dz,std::max(dy1,d <<
904 <<
905 G4double dx1 = fPlanes[2].a*p.x()+fPlane <<
906 + fPlanes[2].c*p.z()+fPlane <<
907 G4double dx2 = fPlanes[3].a*p.x()+fPlane <<
908 + fPlanes[3].c*p.z()+fPlane <<
909 G4double dist = std::max(dy,std::max(dx1 <<
910 return (dist > 0) ? dist : 0.; <<
911 } <<
912 case 1: // YZ section is a rectangle <<
913 { <<
914 G4double dz = std::abs(p.z())-fDz; <<
915 G4double dy = std::max(dz,std::abs(p.y() <<
916 G4double dx1 = fPlanes[2].a*p.x()+fPlane <<
917 + fPlanes[2].c*p.z()+fPlane <<
918 G4double dx2 = fPlanes[3].a*p.x()+fPlane <<
919 + fPlanes[3].c*p.z()+fPlane <<
920 G4double dist = std::max(dy,std::max(dx1 <<
921 return (dist > 0) ? dist : 0.; <<
922 } <<
923 case 2: // YZ section is a rectangle and <<
924 { // XZ section is an isosceles trap <<
925 G4double dz = std::abs(p.z())-fDz; <<
926 G4double dy = std::max(dz,std::abs(p.y() <<
927 G4double dx = fPlanes[3].a*std::abs(p.x( <<
928 + fPlanes[3].c*p.z()+fPlanes <<
929 G4double dist = std::max(dy,dx); <<
930 return (dist > 0) ? dist : 0.; <<
931 } <<
932 case 3: // YZ section is a rectangle and <<
933 { // XY section is an isosceles trap <<
934 G4double dz = std::abs(p.z())-fDz; <<
935 G4double dy = std::max(dz,std::abs(p.y() <<
936 G4double dx = fPlanes[3].a*std::abs(p.x( <<
937 + fPlanes[3].b*p.y()+fPlanes <<
938 G4double dist = std::max(dy,dx); <<
939 return (dist > 0) ? dist : 0.; <<
940 } <<
941 } 1220 }
942 return 0.; << 1221 if (safe<0) safe=0;
>> 1222 return safe;
943 } 1223 }
944 1224
945 ////////////////////////////////////////////// << 1225 /////////////////////////////////////////////////////////////////////////////////
946 // 1226 //
947 // Calculate distance to surface of shape from << 1227 // Calculate distance to surface of shape from inside
948 // find normal at exit point, if required << 1228 // Calculate distance to x/y/z planes - smallest is exiting distance
949 // - when leaving the surface, return 0 <<
950 1229
951 G4double G4Trap::DistanceToOut(const G4ThreeVe 1230 G4double G4Trap::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
952 const G4bool ca 1231 const G4bool calcNorm,
953 G4bool* v << 1232 G4bool *validNorm, G4ThreeVector *n) const
954 { 1233 {
955 // Z intersections << 1234 Eside side = kUndef;
>> 1235 G4double snxt; // snxt = return value
>> 1236 G4double pdist,Comp,vdist,max;
>> 1237 //
>> 1238 // Z Intersections
956 // 1239 //
957 if ((std::abs(p.z()) - fDz) >= -halfCarToler << 1240 if (v.z()>0)
958 { 1241 {
959 if (calcNorm) << 1242 max=fDz-p.z();
>> 1243 if (max>kCarTolerance/2)
>> 1244 {
>> 1245 snxt=max/v.z();
>> 1246 side=kPZ;
>> 1247 }
>> 1248 else
960 { 1249 {
961 *validNorm = true; << 1250 if (calcNorm)
962 n->set(0, 0, (p.z() < 0) ? -1 : 1); << 1251 {
>> 1252 *validNorm=true;
>> 1253 *n=G4ThreeVector(0,0,1);
>> 1254 }
>> 1255 return snxt=0;
963 } 1256 }
964 return 0; <<
965 } 1257 }
966 G4double vz = v.z(); << 1258 else if (v.z()<0)
967 G4double tmax = (vz == 0) ? DBL_MAX : (std:: << 1259 {
968 G4int iside = (vz < 0) ? -4 : -2; // little << 1260 max=-fDz-p.z();
>> 1261 if (max<-kCarTolerance/2)
>> 1262 {
>> 1263 snxt=max/v.z();
>> 1264 side=kMZ;
>> 1265 }
>> 1266 else
>> 1267 {
>> 1268 if (calcNorm)
>> 1269 {
>> 1270 *validNorm=true;
>> 1271 *n=G4ThreeVector(0,0,-1);
>> 1272 }
>> 1273 return snxt=0;
>> 1274 }
>> 1275 }
>> 1276 else
>> 1277 {
>> 1278 snxt=kInfinity;
>> 1279 }
969 1280
970 // Y intersections <<
971 // 1281 //
972 G4int i = 0; << 1282 // Intersections with planes[0] (expanded because of setting enum)
973 for ( ; i<2; ++i) << 1283 //
>> 1284 pdist=fPlanes[0].a*p.x()+fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
>> 1285 Comp=fPlanes[0].a*v.x()+fPlanes[0].b*v.y()+fPlanes[0].c*v.z();
>> 1286 if (pdist>0)
974 { 1287 {
975 G4double cosa = fPlanes[i].b*v.y() + fPlan << 1288 // Outside the plane
976 if (cosa > 0) << 1289 if (Comp>0)
977 { 1290 {
978 G4double dist = fPlanes[i].b*p.y() + fPl << 1291 // Leaving immediately
979 if (dist >= -halfCarTolerance) << 1292 if (calcNorm)
980 { 1293 {
981 if (calcNorm) << 1294 *validNorm=true;
982 { << 1295 *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
983 *validNorm = true; << 1296 }
984 n->set(0, fPlanes[i].b, fPlanes[i].c << 1297 return snxt=0;
985 } << 1298 }
986 return 0; << 1299 }
>> 1300 else if (pdist<-kCarTolerance/2)
>> 1301 {
>> 1302 // Inside the plane
>> 1303 if (Comp>0)
>> 1304 {
>> 1305 // Will leave extent
>> 1306 vdist=-pdist/Comp;
>> 1307 if (vdist<snxt)
>> 1308 {
>> 1309 snxt=vdist;
>> 1310 side=ks0;
987 } 1311 }
988 G4double tmp = -dist/cosa; << 1312 }
989 if (tmax > tmp) { tmax = tmp; iside = i; << 1313 }
>> 1314 else
>> 1315 {
>> 1316 // On surface
>> 1317 if (Comp>0)
>> 1318 {
>> 1319 if (calcNorm)
>> 1320 {
>> 1321 *validNorm=true;
>> 1322 *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
>> 1323 }
>> 1324 return snxt=0;
990 } 1325 }
991 } 1326 }
992 1327
993 // X intersections <<
994 // 1328 //
995 for ( ; i<4; ++i) << 1329 // Intersections with planes[1] (expanded because of setting enum)
>> 1330 //
>> 1331 pdist=fPlanes[1].a*p.x()+fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
>> 1332 Comp=fPlanes[1].a*v.x()+fPlanes[1].b*v.y()+fPlanes[1].c*v.z();
>> 1333 if (pdist>0)
996 { 1334 {
997 G4double cosa = fPlanes[i].a*v.x()+fPlanes << 1335 // Outside the plane
998 if (cosa > 0) << 1336 if (Comp>0)
999 { 1337 {
1000 G4double dist = fPlanes[i].a*p.x() + << 1338 // Leaving immediately
1001 fPlanes[i].b*p.y() + fP << 1339 if (calcNorm)
1002 if (dist >= -halfCarTolerance) <<
1003 { 1340 {
1004 if (calcNorm) << 1341 *validNorm=true;
1005 { << 1342 *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
1006 *validNorm = true; <<
1007 n->set(fPlanes[i].a, fPlanes[i].b, <<
1008 } <<
1009 return 0; <<
1010 } 1343 }
1011 G4double tmp = -dist/cosa; << 1344 return snxt=0;
1012 if (tmax > tmp) { tmax = tmp; iside = i << 1345 }
>> 1346 }
>> 1347 else if (pdist<-kCarTolerance/2)
>> 1348 {
>> 1349 // Inside the plane
>> 1350 if (Comp>0)
>> 1351 {
>> 1352 // Will leave extent
>> 1353 vdist=-pdist/Comp;
>> 1354 if (vdist<snxt)
>> 1355 {
>> 1356 snxt=vdist;
>> 1357 side=ks1;
>> 1358 }
>> 1359 }
>> 1360 }
>> 1361 else
>> 1362 {
>> 1363 // On surface
>> 1364 if (Comp>0)
>> 1365 {
>> 1366 if (calcNorm)
>> 1367 {
>> 1368 *validNorm=true;
>> 1369 *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
>> 1370 }
>> 1371 return snxt=0;
1013 } 1372 }
1014 } 1373 }
1015 1374
1016 // Set normal, if required, and return dist <<
1017 // 1375 //
>> 1376 // Intersections with planes[2] (expanded because of setting enum)
>> 1377 //
>> 1378 pdist=fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z()+fPlanes[2].d;
>> 1379 Comp=fPlanes[2].a*v.x()+fPlanes[2].b*v.y()+fPlanes[2].c*v.z();
>> 1380 if (pdist>0)
>> 1381 {
>> 1382 // Outside the plane
>> 1383 if (Comp>0)
>> 1384 {
>> 1385 // Leaving immediately
>> 1386 if (calcNorm)
>> 1387 {
>> 1388 *validNorm=true;
>> 1389 *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
>> 1390 }
>> 1391 return snxt=0;
>> 1392 }
>> 1393 }
>> 1394 else if (pdist<-kCarTolerance/2)
>> 1395 {
>> 1396 // Inside the plane
>> 1397 if (Comp>0)
>> 1398 {
>> 1399 // Will leave extent
>> 1400 vdist=-pdist/Comp;
>> 1401 if (vdist<snxt)
>> 1402 {
>> 1403 snxt=vdist;
>> 1404 side=ks2;
>> 1405 }
>> 1406 }
>> 1407 }
>> 1408 else
>> 1409 {
>> 1410 // On surface
>> 1411 if (Comp>0)
>> 1412 {
>> 1413 if (calcNorm)
>> 1414 {
>> 1415 *validNorm=true;
>> 1416 *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
>> 1417 }
>> 1418 return snxt=0;
>> 1419 }
>> 1420 }
>> 1421
>> 1422 //
>> 1423 // Intersections with planes[3] (expanded because of setting enum)
>> 1424 //
>> 1425 pdist=fPlanes[3].a*p.x()+fPlanes[3].b*p.y()+fPlanes[3].c*p.z()+fPlanes[3].d;
>> 1426 Comp=fPlanes[3].a*v.x()+fPlanes[3].b*v.y()+fPlanes[3].c*v.z();
>> 1427 if (pdist>0)
>> 1428 {
>> 1429 // Outside the plane
>> 1430 if (Comp>0)
>> 1431 {
>> 1432 // Leaving immediately
>> 1433 if (calcNorm)
>> 1434 {
>> 1435 *validNorm=true;
>> 1436 *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
>> 1437 }
>> 1438 return snxt=0;
>> 1439 }
>> 1440 }
>> 1441 else if (pdist<-kCarTolerance/2)
>> 1442 {
>> 1443 // Inside the plane
>> 1444 if (Comp>0)
>> 1445 {
>> 1446 // Will leave extent
>> 1447 vdist=-pdist/Comp;
>> 1448 if (vdist<snxt)
>> 1449 {
>> 1450 snxt=vdist;
>> 1451 side=ks3;
>> 1452 }
>> 1453 }
>> 1454 }
>> 1455 else
>> 1456 {
>> 1457 // On surface
>> 1458 if (Comp>0)
>> 1459 {
>> 1460 if (calcNorm)
>> 1461 {
>> 1462 *validNorm=true;
>> 1463 *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
>> 1464 }
>> 1465 return snxt=0;
>> 1466 }
>> 1467 }
>> 1468
>> 1469 // set normal
1018 if (calcNorm) 1470 if (calcNorm)
1019 { 1471 {
1020 *validNorm = true; << 1472 *validNorm=true;
1021 if (iside < 0) << 1473 switch(side)
1022 n->set(0, 0, iside + 3); // (-4+3)=-1, << 1474 {
1023 else << 1475 case ks0:
1024 n->set(fPlanes[iside].a, fPlanes[iside] << 1476 *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
>> 1477 break;
>> 1478 case ks1:
>> 1479 *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
>> 1480 break;
>> 1481 case ks2:
>> 1482 *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
>> 1483 break;
>> 1484 case ks3:
>> 1485 *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
>> 1486 break;
>> 1487 case kMZ:
>> 1488 *n=G4ThreeVector(0,0,-1);
>> 1489 break;
>> 1490 case kPZ:
>> 1491 *n=G4ThreeVector(0,0,1);
>> 1492 break;
>> 1493 default:
>> 1494 G4cout << G4endl;
>> 1495 DumpInfo();
>> 1496 std::ostringstream message;
>> 1497 G4int oldprc = message.precision(16);
>> 1498 message << "Undefined side for valid surface normal to solid."
>> 1499 << G4endl
>> 1500 << "Position:" << G4endl << G4endl
>> 1501 << "p.x() = " << p.x()/mm << " mm" << G4endl
>> 1502 << "p.y() = " << p.y()/mm << " mm" << G4endl
>> 1503 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
>> 1504 << "Direction:" << G4endl << G4endl
>> 1505 << "v.x() = " << v.x() << G4endl
>> 1506 << "v.y() = " << v.y() << G4endl
>> 1507 << "v.z() = " << v.z() << G4endl << G4endl
>> 1508 << "Proposed distance :" << G4endl << G4endl
>> 1509 << "snxt = " << snxt/mm << " mm" << G4endl;
>> 1510 message.precision(oldprc);
>> 1511 G4Exception("G4Trap::DistanceToOut(p,v,..)","GeomSolids1002",
>> 1512 JustWarning, message);
>> 1513 break;
>> 1514 }
1025 } 1515 }
1026 return tmax; << 1516 return snxt;
1027 } 1517 }
1028 1518
1029 ///////////////////////////////////////////// << 1519 //////////////////////////////////////////////////////////////////////////////
1030 // 1520 //
1031 // Calculate exact shortest distance to any b 1521 // Calculate exact shortest distance to any boundary from inside
1032 // - Returns 0 is ThreeVector outside 1522 // - Returns 0 is ThreeVector outside
1033 1523
1034 G4double G4Trap::DistanceToOut( const G4Three 1524 G4double G4Trap::DistanceToOut( const G4ThreeVector& p ) const
1035 { 1525 {
>> 1526 G4double safe=0.0,Dist;
>> 1527 G4int i;
>> 1528
1036 #ifdef G4CSGDEBUG 1529 #ifdef G4CSGDEBUG
1037 if( Inside(p) == kOutside ) 1530 if( Inside(p) == kOutside )
1038 { 1531 {
1039 std::ostringstream message; << 1532 G4int oldprc = G4cout.precision(16) ;
1040 G4long oldprc = message.precision(16); << 1533 G4cout << G4endl ;
1041 message << "Point p is outside (!?) of so << 1534 DumpInfo();
1042 message << "Position:\n"; << 1535 G4cout << "Position:" << G4endl << G4endl ;
1043 message << " p.x() = " << p.x()/mm << " << 1536 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1044 message << " p.y() = " << p.y()/mm << " << 1537 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1045 message << " p.z() = " << p.z()/mm << " << 1538 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1046 G4cout.precision(oldprc); << 1539 G4cout.precision(oldprc) ;
1047 G4Exception("G4Trap::DistanceToOut(p)", " << 1540 G4Exception("G4Trap::DistanceToOut(p)",
1048 JustWarning, message ); << 1541 "GeomSolids1002", JustWarning, "Point p is outside !?" );
1049 DumpInfo(); <<
1050 } 1542 }
1051 #endif 1543 #endif
1052 switch (fTrapType) << 1544
>> 1545 safe=fDz-std::fabs(p.z());
>> 1546 if (safe<0) safe=0;
>> 1547 else
1053 { 1548 {
1054 case 0: // General case << 1549 for (i=0;i<4;i++)
1055 { 1550 {
1056 G4double dz = std::abs(p.z())-fDz; << 1551 Dist=-(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
1057 G4double dy1 = fPlanes[0].b*p.y()+fPlan << 1552 +fPlanes[i].c*p.z()+fPlanes[i].d);
1058 G4double dy2 = fPlanes[1].b*p.y()+fPlan << 1553 if (Dist<safe) safe=Dist;
1059 G4double dy = std::max(dz,std::max(dy1, <<
1060 <<
1061 G4double dx1 = fPlanes[2].a*p.x()+fPlan <<
1062 + fPlanes[2].c*p.z()+fPlan <<
1063 G4double dx2 = fPlanes[3].a*p.x()+fPlan <<
1064 + fPlanes[3].c*p.z()+fPlan <<
1065 G4double dist = std::max(dy,std::max(dx <<
1066 return (dist < 0) ? -dist : 0.; <<
1067 } <<
1068 case 1: // YZ section is a rectangle <<
1069 { <<
1070 G4double dz = std::abs(p.z())-fDz; <<
1071 G4double dy = std::max(dz,std::abs(p.y( <<
1072 G4double dx1 = fPlanes[2].a*p.x()+fPlan <<
1073 + fPlanes[2].c*p.z()+fPlan <<
1074 G4double dx2 = fPlanes[3].a*p.x()+fPlan <<
1075 + fPlanes[3].c*p.z()+fPlan <<
1076 G4double dist = std::max(dy,std::max(dx <<
1077 return (dist < 0) ? -dist : 0.; <<
1078 } <<
1079 case 2: // YZ section is a rectangle and <<
1080 { // XZ section is an isosceles tra <<
1081 G4double dz = std::abs(p.z())-fDz; <<
1082 G4double dy = std::max(dz,std::abs(p.y( <<
1083 G4double dx = fPlanes[3].a*std::abs(p.x <<
1084 + fPlanes[3].c*p.z()+fPlane <<
1085 G4double dist = std::max(dy,dx); <<
1086 return (dist < 0) ? -dist : 0.; <<
1087 } <<
1088 case 3: // YZ section is a rectangle and <<
1089 { // XY section is an isosceles tra <<
1090 G4double dz = std::abs(p.z())-fDz; <<
1091 G4double dy = std::max(dz,std::abs(p.y( <<
1092 G4double dx = fPlanes[3].a*std::abs(p.x <<
1093 + fPlanes[3].b*p.y()+fPlane <<
1094 G4double dist = std::max(dy,dx); <<
1095 return (dist < 0) ? -dist : 0.; <<
1096 } 1554 }
>> 1555 if (safe<0) safe=0;
1097 } 1556 }
1098 return 0.; << 1557 return safe;
1099 } 1558 }
1100 1559
1101 ///////////////////////////////////////////// 1560 //////////////////////////////////////////////////////////////////////////
1102 // 1561 //
1103 // GetEntityType 1562 // GetEntityType
1104 1563
1105 G4GeometryType G4Trap::GetEntityType() const 1564 G4GeometryType G4Trap::GetEntityType() const
1106 { 1565 {
1107 return {"G4Trap"}; << 1566 return G4String("G4Trap");
1108 } <<
1109 <<
1110 ///////////////////////////////////////////// <<
1111 // <<
1112 // IsFaceted <<
1113 <<
1114 G4bool G4Trap::IsFaceted() const <<
1115 { <<
1116 return true; <<
1117 } 1567 }
1118 1568
1119 ///////////////////////////////////////////// 1569 //////////////////////////////////////////////////////////////////////////
1120 // 1570 //
1121 // Make a clone of the object 1571 // Make a clone of the object
1122 // 1572 //
1123 G4VSolid* G4Trap::Clone() const 1573 G4VSolid* G4Trap::Clone() const
1124 { 1574 {
1125 return new G4Trap(*this); 1575 return new G4Trap(*this);
1126 } 1576 }
1127 1577
1128 ///////////////////////////////////////////// 1578 //////////////////////////////////////////////////////////////////////////
1129 // 1579 //
1130 // Stream object contents to an output stream 1580 // Stream object contents to an output stream
1131 1581
1132 std::ostream& G4Trap::StreamInfo( std::ostrea 1582 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
1133 { 1583 {
1134 G4double phi = GetPhi(); << 1584 G4int oldprc = os.precision(16);
1135 G4double theta = GetTheta(); <<
1136 G4double alpha1 = GetAlpha1(); <<
1137 G4double alpha2 = GetAlpha2(); <<
1138 <<
1139 G4long oldprc = os.precision(16); <<
1140 os << "------------------------------------ 1585 os << "-----------------------------------------------------------\n"
1141 << " *** Dump for solid: " << GetName << 1586 << " *** Dump for solid - " << GetName() << " ***\n"
1142 << " ================================ 1587 << " ===================================================\n"
1143 << " Solid type: G4Trap\n" 1588 << " Solid type: G4Trap\n"
1144 << " Parameters:\n" << 1589 << " Parameters: \n"
1145 << " half length Z: " << fDz/mm << " << 1590 << " half length Z: " << fDz/mm << " mm \n"
1146 << " half length Y, face -Dz: " << fD << 1591 << " half length Y of face -fDz: " << fDy1/mm << " mm \n"
1147 << " half length X, face -Dz, side -D << 1592 << " half length X of side -fDy1, face -fDz: " << fDx1/mm << " mm \n"
1148 << " half length X, face -Dz, side +D << 1593 << " half length X of side +fDy1, face -fDz: " << fDx2/mm << " mm \n"
1149 << " half length Y, face +Dz: " << fD << 1594 << " half length Y of face +fDz: " << fDy2/mm << " mm \n"
1150 << " half length X, face +Dz, side -D << 1595 << " half length X of side -fDy2, face +fDz: " << fDx3/mm << " mm \n"
1151 << " half length X, face +Dz, side +D << 1596 << " half length X of side +fDy2, face +fDz: " << fDx4/mm << " mm \n"
1152 << " theta: " << theta/degree << " de << 1597 << " std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree << " degrees \n"
1153 << " phi: " << phi/degree << " degr << 1598 << " std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree << " degrees \n"
1154 << " alpha, face -Dz: " << alpha1/deg << 1599 << " std::tan(alpha), -fDz: " << fTalpha1/degree << " degrees \n"
1155 << " alpha, face +Dz: " << alpha2/deg << 1600 << " std::tan(alpha), +fDz: " << fTalpha2/degree << " degrees \n"
>> 1601 << " trap side plane equations:\n"
>> 1602 << " " << fPlanes[0].a << " X + " << fPlanes[0].b << " Y + "
>> 1603 << fPlanes[0].c << " Z + " << fPlanes[0].d << " = 0\n"
>> 1604 << " " << fPlanes[1].a << " X + " << fPlanes[1].b << " Y + "
>> 1605 << fPlanes[1].c << " Z + " << fPlanes[1].d << " = 0\n"
>> 1606 << " " << fPlanes[2].a << " X + " << fPlanes[2].b << " Y + "
>> 1607 << fPlanes[2].c << " Z + " << fPlanes[2].d << " = 0\n"
>> 1608 << " " << fPlanes[3].a << " X + " << fPlanes[3].b << " Y + "
>> 1609 << fPlanes[3].c << " Z + " << fPlanes[3].d << " = 0\n"
1156 << "------------------------------------ 1610 << "-----------------------------------------------------------\n";
1157 os.precision(oldprc); 1611 os.precision(oldprc);
1158 1612
1159 return os; 1613 return os;
1160 } 1614 }
1161 1615
1162 ///////////////////////////////////////////// << 1616 /////////////////////////////////////////////////////////////////////////
>> 1617 //
>> 1618 // GetPointOnPlane
1163 // 1619 //
1164 // Compute vertices from planes << 1620 // Auxiliary method for Get Point on Surface
1165 1621
1166 void G4Trap::GetVertices(G4ThreeVector pt[8]) << 1622 G4ThreeVector G4Trap::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1,
1167 { << 1623 G4ThreeVector p2, G4ThreeVector p3,
1168 for (G4int i=0; i<8; ++i) << 1624 G4double& area) const
1169 { << 1625 {
1170 G4int iy = (i==0 || i==1 || i==4 || i==5) << 1626 G4double lambda1, lambda2, chose, aOne, aTwo;
1171 G4int ix = (i==0 || i==2 || i==4 || i==6) << 1627 G4ThreeVector t, u, v, w, Area, normal;
1172 G4double z = (i < 4) ? -fDz : fDz; << 1628
1173 G4double y = -(fPlanes[iy].c*z + fPlanes[ << 1629 t = p1 - p0;
1174 G4double x = -(fPlanes[ix].b*y + fPlanes[ << 1630 u = p2 - p1;
1175 + fPlanes[ix].d)/fPlanes[i << 1631 v = p3 - p2;
1176 pt[i].set(x,y,z); << 1632 w = p0 - p3;
1177 } << 1633
>> 1634 Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
>> 1635 w.z()*v.x() - w.x()*v.z(),
>> 1636 w.x()*v.y() - w.y()*v.x());
>> 1637
>> 1638 aOne = 0.5*Area.mag();
>> 1639
>> 1640 Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
>> 1641 t.z()*u.x() - t.x()*u.z(),
>> 1642 t.x()*u.y() - t.y()*u.x());
>> 1643
>> 1644 aTwo = 0.5*Area.mag();
>> 1645
>> 1646 area = aOne + aTwo;
>> 1647
>> 1648 chose = G4RandFlat::shoot(0.,aOne+aTwo);
>> 1649
>> 1650 if( (chose>=0.) && (chose < aOne) )
>> 1651 {
>> 1652 lambda1 = G4RandFlat::shoot(0.,1.);
>> 1653 lambda2 = G4RandFlat::shoot(0.,lambda1);
>> 1654 return (p2+lambda1*v+lambda2*w);
>> 1655 }
>> 1656
>> 1657 // else
>> 1658
>> 1659 lambda1 = G4RandFlat::shoot(0.,1.);
>> 1660 lambda2 = G4RandFlat::shoot(0.,lambda1);
>> 1661
>> 1662 return (p0+lambda1*t+lambda2*u);
1178 } 1663 }
1179 1664
1180 ///////////////////////////////////////////// << 1665 ///////////////////////////////////////////////////////////////
1181 // 1666 //
1182 // Generate random point on the surface << 1667 // GetPointOnSurface
1183 1668
1184 G4ThreeVector G4Trap::GetPointOnSurface() con 1669 G4ThreeVector G4Trap::GetPointOnSurface() const
1185 { 1670 {
1186 // Set indeces << 1671 G4double aOne, aTwo, aThree, aFour, aFive, aSix, chose;
1187 constexpr G4int iface [6][4] = << 1672 G4ThreeVector One, Two, Three, Four, Five, Six, test;
1188 { {0,1,3,2}, {0,4,5,1}, {2,3,7,6}, {0,2,6 <<
1189 <<
1190 // Set vertices <<
1191 G4ThreeVector pt[8]; 1673 G4ThreeVector pt[8];
1192 GetVertices(pt); << 1674
1193 << 1675 pt[0] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
1194 // Select face << 1676 -fDz*fTthetaSphi-fDy1,-fDz);
1195 // << 1677 pt[1] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
1196 G4double select = fAreas[5]*G4QuickRand(); << 1678 -fDz*fTthetaSphi-fDy1,-fDz);
1197 G4int k = 5; << 1679 pt[2] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
1198 k -= (G4int)(select <= fAreas[4]); << 1680 -fDz*fTthetaSphi+fDy1,-fDz);
1199 k -= (G4int)(select <= fAreas[3]); << 1681 pt[3] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
1200 k -= (G4int)(select <= fAreas[2]); << 1682 -fDz*fTthetaSphi+fDy1,-fDz);
1201 k -= (G4int)(select <= fAreas[1]); << 1683 pt[4] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
1202 k -= (G4int)(select <= fAreas[0]); << 1684 +fDz*fTthetaSphi-fDy2,+fDz);
1203 << 1685 pt[5] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
1204 // Select sub-triangle << 1686 +fDz*fTthetaSphi-fDy2,+fDz);
1205 // << 1687 pt[6] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
1206 G4int i0 = iface[k][0]; << 1688 +fDz*fTthetaSphi+fDy2,+fDz);
1207 G4int i1 = iface[k][1]; << 1689 pt[7] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
1208 G4int i2 = iface[k][2]; << 1690 +fDz*fTthetaSphi+fDy2,+fDz);
1209 G4int i3 = iface[k][3]; << 1691
1210 G4double s2 = G4GeomTools::TriangleAreaNorm << 1692 // make sure we provide the points in a clockwise fashion
1211 if (select > fAreas[k] - s2) i0 = i2; << 1693
1212 << 1694 One = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
1213 // Generate point << 1695 Two = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
1214 // << 1696 Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
1215 G4double u = G4QuickRand(); << 1697 Four = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour);
1216 G4double v = G4QuickRand(); << 1698 Five = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
1217 if (u + v > 1.) { u = 1. - u; v = 1. - v; } << 1699 Six = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
1218 return (1.-u-v)*pt[i0] + u*pt[i1] + v*pt[i3 << 1700
>> 1701 chose = G4RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
>> 1702 if( (chose>=0.) && (chose<aOne) )
>> 1703 { return One; }
>> 1704 else if( (chose>=aOne) && (chose<aOne+aTwo) )
>> 1705 { return Two; }
>> 1706 else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) )
>> 1707 { return Three; }
>> 1708 else if( (chose>=aOne+aTwo+aThree) && (chose<aOne+aTwo+aThree+aFour) )
>> 1709 { return Four; }
>> 1710 else if( (chose>=aOne+aTwo+aThree+aFour)
>> 1711 && (chose<aOne+aTwo+aThree+aFour+aFive) )
>> 1712 { return Five; }
>> 1713 return Six;
1219 } 1714 }
1220 1715
1221 ///////////////////////////////////////////// 1716 //////////////////////////////////////////////////////////////////////////
1222 // 1717 //
1223 // Methods for visualisation 1718 // Methods for visualisation
1224 1719
1225 void G4Trap::DescribeYourselfTo ( G4VGraphics 1720 void G4Trap::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
1226 { 1721 {
1227 scene.AddSolid (*this); 1722 scene.AddSolid (*this);
1228 } 1723 }
1229 1724
1230 G4Polyhedron* G4Trap::CreatePolyhedron () con 1725 G4Polyhedron* G4Trap::CreatePolyhedron () const
1231 { 1726 {
1232 G4double phi = std::atan2(fTthetaSphi, fTth 1727 G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
1233 G4double alpha1 = std::atan(fTalpha1); 1728 G4double alpha1 = std::atan(fTalpha1);
1234 G4double alpha2 = std::atan(fTalpha2); 1729 G4double alpha2 = std::atan(fTalpha2);
1235 G4double theta = std::atan(std::sqrt(fTthet << 1730 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi+fTthetaSphi*fTthetaSphi));
1236 +fTthet <<
1237 1731
1238 return new G4PolyhedronTrap(fDz, theta, phi 1732 return new G4PolyhedronTrap(fDz, theta, phi,
1239 fDy1, fDx1, fDx 1733 fDy1, fDx1, fDx2, alpha1,
1240 fDy2, fDx3, fDx 1734 fDy2, fDx3, fDx4, alpha2);
1241 } 1735 }
1242 1736
1243 #endif 1737 #endif
1244 1738