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