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