| 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: G4Para.cc,v 1.14 2003/06/16 16:53:37 gunter Exp $
>> 25 // GEANT4 tag $Name: geant4-05-02-patch-01 $
>> 26 //
>> 27 // class G4Para
>> 28 //
26 // Implementation for G4Para class 29 // Implementation for G4Para class
27 // 30 //
>> 31 // History:
28 // 21.03.95 P.Kent: Modified for `tolerant' ge 32 // 21.03.95 P.Kent: Modified for `tolerant' geom
29 // 31.10.96 V.Grichine: Modifications accordin 33 // 31.10.96 V.Grichine: Modifications according G4Box/Tubs before to commit
30 // 28.04.05 V.Grichine: new SurfaceNormal acco << 34 // 18.11.99 V.Grichine: kUndef was added to ESide
31 // 29.05.17 E.Tcherniaev: complete revision, s << 35 // 14.02.02 V.Grichine: bug fixed in Inside according to proposal of D.Wright
32 ////////////////////////////////////////////// << 36 // ********************************************************************
33 37
34 #include "G4Para.hh" 38 #include "G4Para.hh"
35 39
36 #if !defined(G4GEOM_USE_UPARA) <<
37 <<
38 #include "G4VoxelLimits.hh" 40 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh" 41 #include "G4AffineTransform.hh"
40 #include "G4BoundingEnvelope.hh" <<
41 #include "Randomize.hh" <<
42 42
43 #include "G4VPVParameterisation.hh" 43 #include "G4VPVParameterisation.hh"
44 44
45 #include "G4VGraphicsScene.hh" 45 #include "G4VGraphicsScene.hh"
>> 46 #include "G4Polyhedron.hh"
>> 47 #include "G4NURBS.hh"
>> 48 #include "G4NURBSbox.hh"
46 49
47 using namespace CLHEP; << 50 // Private enum: Not for external use
>> 51
>> 52 enum ESide {kUndef,kPX,kMX,kPY,kMY,kPZ,kMZ};
48 53
49 ////////////////////////////////////////////// << 54 // used internally for normal routine
50 // <<
51 // Constructor - set & check half widths <<
52 55
53 G4Para::G4Para(const G4String& pName, << 56 enum ENSide {kNZ,kNX,kNY};
54 G4double pDx, G4double pD <<
55 G4double pAlpha, G4double <<
56 : G4CSGSolid(pName), halfCarTolerance(0.5*kC <<
57 { <<
58 SetAllParameters(pDx, pDy, pDz, pAlpha, pThe <<
59 fRebuildPolyhedron = false; // default valu <<
60 } <<
61 57
62 ////////////////////////////////////////////// << 58 /////////////////////////////////////////////////////////////////////
63 // 59 //
64 // Constructor - design of trapezoid based on << 60 // Constructor - check and set half-widths
65 61
66 G4Para::G4Para( const G4String& pName, << 62 void G4Para::SetAllParameters( G4double pDx, G4double pDy, G4double pDz,
67 const G4ThreeVector pt[8] ) << 63 G4double pAlpha, G4double pTheta, G4double pPhi )
68 : G4CSGSolid(pName), halfCarTolerance(0.5*kC <<
69 { 64 {
70 // Find dimensions and trigonometric values << 65 if (pDx>0&&pDy>0&&pDz>0)
71 // << 66 {
72 fDx = (pt[3].x() - pt[2].x())*0.5; << 67 fDx=pDx;
73 fDy = (pt[2].y() - pt[1].y())*0.5; << 68 fDy=pDy;
74 fDz = pt[7].z(); << 69 fDz=pDz;
75 CheckParameters(); // check dimensions << 70 fTalpha=tan(pAlpha);
76 << 71 fTthetaCphi=tan(pTheta)*cos(pPhi);
77 fTalpha = (pt[2].x() + pt[3].x() - pt[1].x() << 72 fTthetaSphi=tan(pTheta)*sin(pPhi);
78 fTthetaCphi = (pt[4].x() + fDy*fTalpha + fDx << 73 }
79 fTthetaSphi = (pt[4].y() + fDy)/fDz; << 74 else
80 MakePlanes(); << 75 {
81 << 76 G4cout << "ERROR - G4Para()::SetAllParameters(): " << GetName() << G4endl
82 // Recompute vertices << 77 << " Invalid dimensions ! - "
83 // << 78 << pDx << ", " << pDy << ", " << pDz << G4endl;
84 G4ThreeVector v[8]; << 79 G4cerr << "ERROR - G4Para()::SetAllParameters(): " << GetName() << G4endl
85 G4double DyTalpha = fDy*fTalpha; << 80 << " Invalid dimensions ! - "
86 G4double DzTthetaSphi = fDz*fTthetaSphi; << 81 << pDx << ", " << pDy << ", " << pDz << G4endl;
87 G4double DzTthetaCphi = fDz*fTthetaCphi; << 82 G4Exception("G4Para::SetAllParameters() - Invalid Length Parameters");
88 v[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTthe <<
89 v[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTthe <<
90 v[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTthe <<
91 v[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTthe <<
92 v[4].set( DzTthetaCphi-DyTalpha-fDx, DzTthe <<
93 v[5].set( DzTthetaCphi-DyTalpha+fDx, DzTthe <<
94 v[6].set( DzTthetaCphi+DyTalpha-fDx, DzTthe <<
95 v[7].set( DzTthetaCphi+DyTalpha+fDx, DzTthe <<
96 <<
97 // Compare with original vertices <<
98 // <<
99 for (G4int i=0; i<8; ++i) <<
100 { <<
101 G4double delx = std::abs(pt[i].x() - v[i]. <<
102 G4double dely = std::abs(pt[i].y() - v[i]. <<
103 G4double delz = std::abs(pt[i].z() - v[i]. <<
104 G4double discrepancy = std::max(std::max(d <<
105 if (discrepancy > 0.1*kCarTolerance) <<
106 { <<
107 std::ostringstream message; <<
108 G4long oldprc = message.precision(16); <<
109 message << "Invalid vertice coordinates <<
110 << "\nVertix #" << i << ", discr <<
111 << "\n original : " << pt[i] <<
112 << "\n recomputed : " << v[i]; <<
113 G4cout.precision(oldprc); <<
114 G4Exception("G4Para::G4Para()", "GeomSol <<
115 FatalException, message); <<
116 <<
117 } <<
118 } 83 }
119 } 84 }
120 85
121 ////////////////////////////////////////////// << 86 ///////////////////////////////////////////////////////////////////////////
122 // <<
123 // Fake default constructor - sets only member <<
124 // for usage restri <<
125 <<
126 G4Para::G4Para( __void__& a ) <<
127 : G4CSGSolid(a), halfCarTolerance(0.5*kCarTo <<
128 { <<
129 SetAllParameters(1., 1., 1., 0., 0., 0.); <<
130 fRebuildPolyhedron = false; // default value <<
131 } <<
132 <<
133 ////////////////////////////////////////////// <<
134 // <<
135 // Destructor <<
136 <<
137 G4Para::~G4Para() = default; <<
138 <<
139 ////////////////////////////////////////////// <<
140 // <<
141 // Copy constructor <<
142 <<
143 G4Para::G4Para(const G4Para& rhs) <<
144 : G4CSGSolid(rhs), halfCarTolerance(rhs.half <<
145 fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), <<
146 fTthetaCphi(rhs.fTthetaCphi),fTthetaSphi(r <<
147 { <<
148 for (G4int i=0; i<4; ++i) { fPlanes[i] = rhs <<
149 } <<
150 <<
151 ////////////////////////////////////////////// <<
152 // <<
153 // Assignment operator <<
154 <<
155 G4Para& G4Para::operator = (const G4Para& rhs) <<
156 { <<
157 // Check assignment to self <<
158 // <<
159 if (this == &rhs) { return *this; } <<
160 <<
161 // Copy base class data <<
162 // <<
163 G4CSGSolid::operator=(rhs); <<
164 <<
165 // Copy data <<
166 // <<
167 halfCarTolerance = rhs.halfCarTolerance; <<
168 fDx = rhs.fDx; <<
169 fDy = rhs.fDy; <<
170 fDz = rhs.fDz; <<
171 fTalpha = rhs.fTalpha; <<
172 fTthetaCphi = rhs.fTthetaCphi; <<
173 fTthetaSphi = rhs.fTthetaSphi; <<
174 for (G4int i=0; i<4; ++i) { fPlanes[i] = rh <<
175 <<
176 return *this; <<
177 } <<
178 <<
179 ////////////////////////////////////////////// <<
180 // <<
181 // Set all parameters, as for constructor - se <<
182 <<
183 void G4Para::SetAllParameters(G4double pDx, G4 <<
184 G4double pAlpha, <<
185 { <<
186 // Reset data of the base class <<
187 fCubicVolume = 0; <<
188 fSurfaceArea = 0; <<
189 fRebuildPolyhedron = true; <<
190 <<
191 // Set parameters <<
192 fDx = pDx; <<
193 fDy = pDy; <<
194 fDz = pDz; <<
195 fTalpha = std::tan(pAlpha); <<
196 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi <<
197 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi <<
198 <<
199 CheckParameters(); <<
200 MakePlanes(); <<
201 } <<
202 <<
203 ////////////////////////////////////////////// <<
204 // 87 //
205 // Check dimensions <<
206 88
207 void G4Para::CheckParameters() << 89 G4Para::G4Para(const G4String& pName,
>> 90 G4double pDx, G4double pDy, G4double pDz,
>> 91 G4double pAlpha, G4double pTheta, G4double pPhi)
>> 92 : G4CSGSolid(pName)
208 { 93 {
209 if (fDx < 2*kCarTolerance || << 94 if (pDx>0&&pDy>0&&pDz>0)
210 fDy < 2*kCarTolerance || << 95 {
211 fDz < 2*kCarTolerance) << 96 SetAllParameters( pDx, pDy, pDz, pAlpha, pTheta, pPhi);
212 { << 97 }
213 std::ostringstream message; << 98 else
214 message << "Invalid (too small or negative << 99 {
215 << GetName() << 100 G4cout << "ERROR - G4Para()::G4Para(): " << GetName() << G4endl
216 << "\n X - " << fDx << 101 << " Invalid dimensions ! - "
217 << "\n Y - " << fDy << 102 << pDx << ", " << pDy << ", " << pDz << G4endl;
218 << "\n Z - " << fDz; << 103 G4cerr << "ERROR - G4Para()::G4Para(): " << GetName() << G4endl
219 G4Exception("G4Para::CheckParameters()", " << 104 << " Invalid dimensions ! - "
220 FatalException, message); << 105 << pDx << ", " << pDy << ", " << pDz << G4endl;
>> 106 G4Exception("G4Para::G4Para() - Invalid Length Parameters");
221 } 107 }
222 } 108 }
223 109
224 ////////////////////////////////////////////// << 110 ////////////////////////////////////////////////////////////////////////
225 // 111 //
226 // Set side planes << 112 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters,
>> 113 // which are its vertices. Checking of planarity with preparation of
>> 114 // fPlanes[] and than calculation of other members
227 115
228 void G4Para::MakePlanes() << 116 G4Para::G4Para( const G4String& pName,
>> 117 const G4ThreeVector pt[8] )
>> 118 : G4CSGSolid(pName)
229 { 119 {
230 G4ThreeVector vx(1, 0, 0); << 120 if ( pt[0].z()<0 && pt[0].z()==pt[1].z() && pt[0].z()==pt[2].z() &&
231 G4ThreeVector vy(fTalpha, 1, 0); << 121 pt[0].z()==pt[3].z() && pt[4].z()>0 && pt[4].z()==pt[5].z() &&
232 G4ThreeVector vz(fTthetaCphi, fTthetaSphi, 1 << 122 pt[4].z()==pt[6].z() && pt[4].z()==pt[7].z() &&
233 << 123 (pt[0].z()+pt[4].z())==0 &&
234 // Set -Y & +Y planes << 124 pt[0].y()==pt[1].y() && pt[2].y()==pt[3].y() &&
235 // << 125 pt[4].y()==pt[5].y() && pt[6].y()==pt[7].y() &&
236 G4ThreeVector ynorm = (vx.cross(vz)).unit(); << 126 (pt[0].y()+pt[2].y()+pt[4].y()+pt[6].y())==0 )
237 << 127 {
238 fPlanes[0].a = 0.; << 128 fDz = (pt[7]).z() ;
239 fPlanes[0].b = ynorm.y(); << 129
240 fPlanes[0].c = ynorm.z(); << 130 fDy = ((pt[2]).y()-(pt[1]).y())*0.5 ;
241 fPlanes[0].d = fPlanes[0].b*fDy; // point (0 << 131 fDx = ((pt[1]).x()-(pt[0]).x())*0.5 ;
242 << 132 fDx = ((pt[3]).x()-(pt[2]).x())*0.5 ;
243 fPlanes[1].a = 0.; << 133 fTalpha = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy ;
244 fPlanes[1].b = -fPlanes[0].b; << 134
245 fPlanes[1].c = -fPlanes[0].c; << 135 // fDy = ((pt[6]).y()-(pt[5]).y())*0.5 ;
246 fPlanes[1].d = fPlanes[0].d; << 136 // fDx = ((pt[5]).x()-(pt[4]).x())*0.5 ;
247 << 137 // fDx = ((pt[7]).x()-(pt[6]).x())*0.5 ;
248 // Set -X & +X planes << 138 // fTalpha = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy ;
249 // <<
250 G4ThreeVector xnorm = (vz.cross(vy)).unit(); <<
251 <<
252 fPlanes[2].a = xnorm.x(); <<
253 fPlanes[2].b = xnorm.y(); <<
254 fPlanes[2].c = xnorm.z(); <<
255 fPlanes[2].d = fPlanes[2].a*fDx; // point (f <<
256 139
257 fPlanes[3].a = -fPlanes[2].a; << 140 fTthetaCphi = ((pt[4]).x()+fDy*fTalpha+fDx)/fDz ;
258 fPlanes[3].b = -fPlanes[2].b; << 141 fTthetaSphi = ((pt[4]).y()+fDy)/fDz ;
259 fPlanes[3].c = -fPlanes[2].c; <<
260 fPlanes[3].d = fPlanes[2].d; <<
261 } <<
262 <<
263 ////////////////////////////////////////////// <<
264 // <<
265 // Get volume <<
266 <<
267 G4double G4Para::GetCubicVolume() <<
268 { <<
269 // It is like G4Box, since para transformati <<
270 if (fCubicVolume == 0) <<
271 { <<
272 fCubicVolume = 8*fDx*fDy*fDz; <<
273 } 142 }
274 return fCubicVolume; << 143 else
>> 144 {
>> 145 G4cout << "ERROR - G4Para()::G4Para(): " << GetName() << G4endl
>> 146 << " Invalid dimensions !" << G4endl;
>> 147 G4cerr << "ERROR - G4Para()::G4Para(): " << GetName() << G4endl
>> 148 << " Invalid dimensions !" << G4endl;
>> 149 G4Exception("G4Para::G4Para() - Invalid vertice coordinates");
>> 150 }
275 } 151 }
276 152
277 ////////////////////////////////////////////// 153 //////////////////////////////////////////////////////////////////////////
278 // 154 //
279 // Get surface area <<
280 155
281 G4double G4Para::GetSurfaceArea() << 156 G4Para::~G4Para()
282 { 157 {
283 if(fSurfaceArea == 0) <<
284 { <<
285 G4ThreeVector vx(fDx, 0, 0); <<
286 G4ThreeVector vy(fDy*fTalpha, fDy, 0); <<
287 G4ThreeVector vz(fDz*fTthetaCphi, fDz*fTth <<
288 <<
289 G4double sxy = fDx*fDy; // (vx.cross(vy)). <<
290 G4double sxz = (vx.cross(vz)).mag(); <<
291 G4double syz = (vy.cross(vz)).mag(); <<
292 <<
293 fSurfaceArea = 8*(sxy+sxz+syz); <<
294 } <<
295 return fSurfaceArea; <<
296 } 158 }
297 159
298 ////////////////////////////////////////////// 160 //////////////////////////////////////////////////////////////////////////
299 // 161 //
300 // Dispatch to parameterisation for replicatio 162 // Dispatch to parameterisation for replication mechanism dimension
301 // computation & modification << 163 // computation & modification.
302 164
303 void G4Para::ComputeDimensions( G4VPVPara 165 void G4Para::ComputeDimensions( G4VPVParameterisation* p,
304 const G4int n, 166 const G4int n,
305 const G4VPhysi 167 const G4VPhysicalVolume* pRep )
306 { 168 {
307 p->ComputeDimensions(*this,n,pRep); 169 p->ComputeDimensions(*this,n,pRep);
308 } 170 }
309 171
310 ////////////////////////////////////////////// <<
311 // <<
312 // Get bounding box <<
313 172
314 void G4Para::BoundingLimits(G4ThreeVector& pMi << 173 //////////////////////////////////////////////////////////////
315 { <<
316 G4double dz = GetZHalfLength(); <<
317 G4double dx = GetXHalfLength(); <<
318 G4double dy = GetYHalfLength(); <<
319 <<
320 G4double x0 = dz*fTthetaCphi; <<
321 G4double x1 = dy*GetTanAlpha(); <<
322 G4double xmin = <<
323 std::min( <<
324 std::min( <<
325 std::min(-x0-x1-dx,-x0+x1-dx),x0-x1-dx),x0 <<
326 G4double xmax = <<
327 std::max( <<
328 std::max( <<
329 std::max(-x0-x1+dx,-x0+x1+dx),x0-x1+dx),x0 <<
330 <<
331 G4double y0 = dz*fTthetaSphi; <<
332 G4double ymin = std::min(-y0-dy,y0-dy); <<
333 G4double ymax = std::max(-y0+dy,y0+dy); <<
334 <<
335 pMin.set(xmin,ymin,-dz); <<
336 pMax.set(xmax,ymax, dz); <<
337 <<
338 // Check correctness of the bounding box <<
339 // <<
340 if (pMin.x() >= pMax.x() || pMin.y() >= pMax <<
341 { <<
342 std::ostringstream message; <<
343 message << "Bad bounding box (min >= max) <<
344 << GetName() << " !" <<
345 << "\npMin = " << pMin <<
346 << "\npMax = " << pMax; <<
347 G4Exception("G4Para::BoundingLimits()", "G <<
348 JustWarning, message); <<
349 DumpInfo(); <<
350 } <<
351 } <<
352 <<
353 ////////////////////////////////////////////// <<
354 // 174 //
355 // Calculate extent under transform and specif 175 // Calculate extent under transform and specified limit
356 176
357 G4bool G4Para::CalculateExtent( const EAxis pA 177 G4bool G4Para::CalculateExtent( const EAxis pAxis,
358 const G4VoxelL 178 const G4VoxelLimits& pVoxelLimit,
359 const G4Affine 179 const G4AffineTransform& pTransform,
360 G4double& 180 G4double& pMin, G4double& pMax ) const
361 { 181 {
362 G4ThreeVector bmin, bmax; << 182 G4bool flag;
363 G4bool exist; <<
364 183
365 // Check bounding box (bbox) << 184 if (!pTransform.IsRotated())
366 // << 185 {
367 BoundingLimits(bmin,bmax); << 186 // Special case handling for unrotated trapezoids
368 G4BoundingEnvelope bbox(bmin,bmax); << 187 // Compute z/x/y/ mins and maxs respecting limits, with early returns
369 #ifdef G4BBOX_EXTENT << 188 // if outside limits. Then switch() on pAxis
370 return bbox.CalculateExtent(pAxis,pVoxelLimi << 189
371 #endif << 190 G4int i ;
372 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 191 G4double xoffset,xMin,xMax;
373 { << 192 G4double yoffset,yMin,yMax;
374 return exist = pMin < pMax; << 193 G4double zoffset,zMin,zMax;
375 } << 194 G4double temp[8] ; // some points for intersection with zMin/zMax
>> 195
>> 196 xoffset=pTransform.NetTranslation().x();
>> 197 yoffset=pTransform.NetTranslation().y();
>> 198 zoffset=pTransform.NetTranslation().z();
>> 199
>> 200 G4ThreeVector pt[8]; // vertices after translation
>> 201 pt[0]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha-fDx,
>> 202 yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
>> 203 pt[1]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha+fDx,
>> 204 yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
>> 205 pt[2]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha-fDx,
>> 206 yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
>> 207 pt[3]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha+fDx,
>> 208 yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
>> 209 pt[4]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha-fDx,
>> 210 yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
>> 211 pt[5]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha+fDx,
>> 212 yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
>> 213 pt[6]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha-fDx,
>> 214 yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
>> 215 pt[7]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha+fDx,
>> 216 yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
>> 217 zMin=zoffset-fDz;
>> 218 zMax=zoffset+fDz;
>> 219 if ( pVoxelLimit.IsZLimited() )
>> 220 {
>> 221 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
>> 222 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
>> 223 {
>> 224 return false;
>> 225 }
>> 226 else
>> 227 {
>> 228 if (zMin<pVoxelLimit.GetMinZExtent())
>> 229 {
>> 230 zMin=pVoxelLimit.GetMinZExtent();
>> 231 }
>> 232 if (zMax>pVoxelLimit.GetMaxZExtent())
>> 233 {
>> 234 zMax=pVoxelLimit.GetMaxZExtent();
>> 235 }
>> 236 }
>> 237 }
376 238
377 // Set bounding envelope (benv) and calculat << 239 temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())
378 // << 240 *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
379 G4double dz = GetZHalfLength(); << 241 temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())
380 G4double dx = GetXHalfLength(); << 242 *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
381 G4double dy = GetYHalfLength(); << 243 temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())
>> 244 *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
>> 245 temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())
>> 246 *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
>> 247 yMax = yoffset - fabs(fDz*fTthetaSphi) - fDy - fDy ;
>> 248 yMin = -yMax ;
>> 249 for(i=0;i<4;i++)
>> 250 {
>> 251 if(temp[i] > yMax) yMax = temp[i] ;
>> 252 if(temp[i] < yMin) yMin = temp[i] ;
>> 253 }
>> 254
>> 255 if (pVoxelLimit.IsYLimited())
>> 256 {
>> 257 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
>> 258 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
>> 259 {
>> 260 return false;
>> 261 }
>> 262 else
>> 263 {
>> 264 if (yMin<pVoxelLimit.GetMinYExtent())
>> 265 {
>> 266 yMin=pVoxelLimit.GetMinYExtent();
>> 267 }
>> 268 if (yMax>pVoxelLimit.GetMaxYExtent())
>> 269 {
>> 270 yMax=pVoxelLimit.GetMaxYExtent();
>> 271 }
>> 272 }
>> 273 }
382 274
383 G4double x0 = dz*fTthetaCphi; << 275 temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())
384 G4double x1 = dy*GetTanAlpha(); << 276 *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
385 G4double y0 = dz*fTthetaSphi; << 277 temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())
>> 278 *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
>> 279 temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())
>> 280 *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
>> 281 temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())
>> 282 *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
>> 283 temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())
>> 284 *(zMin-pt[3].z())/(pt[7].z()-pt[3].z()) ;
>> 285 temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())
>> 286 *(zMax-pt[3].z())/(pt[7].z()-pt[3].z()) ;
>> 287 temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())
>> 288 *(zMin-pt[1].z())/(pt[5].z()-pt[1].z()) ;
>> 289 temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())
>> 290 *(zMax-pt[1].z())/(pt[5].z()-pt[1].z()) ;
>> 291
>> 292 xMax = xoffset - fabs(fDz*fTthetaCphi) - fDx - fDx -fDx - fDx;
>> 293 xMin = -xMax ;
>> 294 for(i=0;i<8;i++)
>> 295 {
>> 296 if(temp[i] > xMax) xMax = temp[i] ;
>> 297 if(temp[i] < xMin) xMin = temp[i] ;
>> 298 }
>> 299 // xMax/Min = f(yMax/Min) ?
>> 300 if (pVoxelLimit.IsXLimited())
>> 301 {
>> 302 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
>> 303 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
>> 304 {
>> 305 return false;
>> 306 }
>> 307 else
>> 308 {
>> 309 if (xMin<pVoxelLimit.GetMinXExtent())
>> 310 {
>> 311 xMin=pVoxelLimit.GetMinXExtent();
>> 312 }
>> 313 if (xMax>pVoxelLimit.GetMaxXExtent())
>> 314 {
>> 315 xMax=pVoxelLimit.GetMaxXExtent();
>> 316 }
>> 317 }
>> 318 }
>> 319
>> 320 switch (pAxis)
>> 321 {
>> 322 case kXAxis:
>> 323 pMin=xMin;
>> 324 pMax=xMax;
>> 325 break;
>> 326 case kYAxis:
>> 327 pMin=yMin;
>> 328 pMax=yMax;
>> 329 break;
>> 330 case kZAxis:
>> 331 pMin=zMin;
>> 332 pMax=zMax;
>> 333 break;
>> 334 default:
>> 335 break;
>> 336 }
386 337
387 G4ThreeVectorList baseA(4), baseB(4); << 338 pMin-=kCarTolerance;
388 baseA[0].set(-x0-x1-dx,-y0-dy,-dz); << 339 pMax+=kCarTolerance;
389 baseA[1].set(-x0-x1+dx,-y0-dy,-dz); << 340 flag = true;
390 baseA[2].set(-x0+x1+dx,-y0+dy,-dz); << 341 }
391 baseA[3].set(-x0+x1-dx,-y0+dy,-dz); << 342 else
>> 343 {
>> 344 // General rotated case - create and clip mesh to boundaries
392 345
393 baseB[0].set(+x0-x1-dx, y0-dy, dz); << 346 G4bool existsAfterClip=false;
394 baseB[1].set(+x0-x1+dx, y0-dy, dz); << 347 G4ThreeVectorList *vertices;
395 baseB[2].set(+x0+x1+dx, y0+dy, dz); <<
396 baseB[3].set(+x0+x1-dx, y0+dy, dz); <<
397 348
398 std::vector<const G4ThreeVectorList *> polyg << 349 pMin=+kInfinity;
399 polygons[0] = &baseA; << 350 pMax=-kInfinity;
400 polygons[1] = &baseB; <<
401 351
402 G4BoundingEnvelope benv(bmin,bmax,polygons); << 352 // Calculate rotated vertex coordinates
403 exist = benv.CalculateExtent(pAxis,pVoxelLim << 353
404 return exist; << 354 vertices=CreateRotatedVertices(pTransform);
>> 355 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
>> 356 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
>> 357 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
>> 358
>> 359 if (pMin!=kInfinity||pMax!=-kInfinity)
>> 360 {
>> 361 existsAfterClip=true;
>> 362
>> 363 // Add 2*tolerance to avoid precision troubles
>> 364 //
>> 365 pMin-=kCarTolerance;
>> 366 pMax+=kCarTolerance;
>> 367 }
>> 368 else
>> 369 {
>> 370 // Check for case where completely enveloping clipping volume
>> 371 // If point inside then we are confident that the solid completely
>> 372 // envelopes the clipping volume. Hence set min/max extents according
>> 373 // to clipping volume extents along the specified axis.
>> 374
>> 375 G4ThreeVector clipCentre(
>> 376 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
>> 377 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
>> 378 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
>> 379
>> 380 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
>> 381 {
>> 382 existsAfterClip=true;
>> 383 pMin=pVoxelLimit.GetMinExtent(pAxis);
>> 384 pMax=pVoxelLimit.GetMaxExtent(pAxis);
>> 385 }
>> 386 }
>> 387 delete vertices ; // 'new' in the function called
>> 388 flag = existsAfterClip ;
>> 389 }
>> 390 return flag;
405 } 391 }
406 392
407 ////////////////////////////////////////////// << 393 /////////////////////////////////////////////////////////////////////////////
408 // <<
409 // Determine where is point p, inside/on_surfa <<
410 // 394 //
>> 395 // Check in p is inside/on surface/outside solid
411 396
412 EInside G4Para::Inside( const G4ThreeVector& p 397 EInside G4Para::Inside( const G4ThreeVector& p ) const
413 { 398 {
414 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 399 G4double xt, yt, yt1;
415 G4double dx = std::abs(xx) + fPlanes[2].d; << 400 EInside in = kOutside;
416 401
417 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 402 yt1 = p.y() - fTthetaSphi*p.z();
418 G4double dy = std::abs(yy) + fPlanes[0].d; << 403 yt = fabs(yt1) ;
419 G4double dxy = std::max(dx,dy); <<
420 404
421 G4double dz = std::abs(p.z())-fDz; << 405 // xt = fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt );
422 G4double dist = std::max(dxy,dz); <<
423 406
424 if (dist > halfCarTolerance) return kOutside << 407 xt = fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt1 );
425 return (dist > -halfCarTolerance) ? kSurface << 408
>> 409 if ( fabs( p.z() ) <= fDz - kCarTolerance*0.5)
>> 410 {
>> 411 if (yt <= fDy - kCarTolerance*0.5)
>> 412 {
>> 413 if ( xt <= fDx - kCarTolerance*0.5 ) in = kInside;
>> 414 else if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
>> 415 }
>> 416 else if ( yt <= fDy + kCarTolerance*0.5)
>> 417 {
>> 418 if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
>> 419 }
>> 420 }
>> 421 else if ( fabs(p.z()) <= fDz + kCarTolerance*0.5 )
>> 422 {
>> 423 if ( yt <= fDy + kCarTolerance*0.5)
>> 424 {
>> 425 if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
>> 426 }
>> 427 }
>> 428 return in;
426 } 429 }
427 430
428 ////////////////////////////////////////////// << 431 ///////////////////////////////////////////////////////////////////////////
429 // 432 //
430 // Determine side where point is, and return c << 433 // Calculate side nearest to p, and return normal
>> 434 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
431 435
432 G4ThreeVector G4Para::SurfaceNormal( const G4T 436 G4ThreeVector G4Para::SurfaceNormal( const G4ThreeVector& p ) const
433 { 437 {
434 G4int nsurf = 0; // number of surfaces where << 438 ENSide side;
>> 439 G4ThreeVector norm;
>> 440 G4double distx,disty,distz;
>> 441 G4double newpx,newpy,xshift;
>> 442 G4double calpha,salpha; // Sin/Cos(alpha) - needed to recalc G4Parameter
>> 443 G4double tntheta,cosntheta; // tan and cos of normal's theta component
>> 444 G4double ycomp;
435 445
436 // Check Z faces << 446 newpx=p.x()-fTthetaCphi*p.z();
437 // << 447 newpy=p.y()-fTthetaSphi*p.z();
438 G4double nz = 0; <<
439 G4double dz = std::abs(p.z()) - fDz; <<
440 if (std::abs(dz) <= halfCarTolerance) <<
441 { <<
442 nz = (p.z() < 0) ? -1 : 1; <<
443 ++nsurf; <<
444 } <<
445 448
446 // Check Y faces << 449 calpha=1/sqrt(1+fTalpha*fTalpha);
447 // << 450 if (fTalpha)
448 G4double ny = 0; <<
449 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. <<
450 if (std::abs(fPlanes[0].d + yy) <= halfCarTo <<
451 { 451 {
452 ny = fPlanes[0].b; << 452 salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
453 nz += fPlanes[0].c; <<
454 ++nsurf; <<
455 } 453 }
456 else if (std::abs(fPlanes[1].d - yy) <= half << 454 else
457 { 455 {
458 ny = fPlanes[1].b; << 456 salpha=0;
459 nz += fPlanes[1].c; <<
460 ++nsurf; <<
461 } 457 }
462 458
463 // Check X faces << 459 xshift=newpx*calpha+newpy*salpha;
464 // << 460
465 G4double nx = 0; << 461 distx=fabs(fabs(xshift)-fDx*calpha);
466 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 462 disty=fabs(fabs(newpy)-fDy);
467 if (std::abs(fPlanes[2].d + xx) <= halfCarTo << 463 distz=fabs(fabs(p.z())-fDz);
468 { << 464
469 nx = fPlanes[2].a; << 465 if (distx<disty)
470 ny += fPlanes[2].b; <<
471 nz += fPlanes[2].c; <<
472 ++nsurf; <<
473 } <<
474 else if (std::abs(fPlanes[3].d - xx) <= half <<
475 { 466 {
476 nx = fPlanes[3].a; << 467 if (distx<distz) side=kNX;
477 ny += fPlanes[3].b; << 468 else side=kNZ;
478 nz += fPlanes[3].c; <<
479 ++nsurf; <<
480 } 469 }
481 <<
482 // Return normal <<
483 // <<
484 if (nsurf == 1) return {nx,ny,nz}; <<
485 else if (nsurf != 0) return G4ThreeVector(nx <<
486 else 470 else
487 { 471 {
488 // Point is not on the surface << 472 if (disty<distz) side=kNY;
489 // << 473 else side=kNZ;
490 #ifdef G4CSGDEBUG <<
491 std::ostringstream message; <<
492 G4int oldprc = message.precision(16); <<
493 message << "Point p is not on surface (!?) <<
494 << GetName() << G4endl; <<
495 message << "Position:\n"; <<
496 message << " p.x() = " << p.x()/mm << " <<
497 message << " p.y() = " << p.y()/mm << " <<
498 message << " p.z() = " << p.z()/mm << " <<
499 G4cout.precision(oldprc) ; <<
500 G4Exception("G4Para::SurfaceNormal(p)", "G <<
501 JustWarning, message ); <<
502 DumpInfo(); <<
503 #endif <<
504 return ApproxSurfaceNormal(p); <<
505 } 474 }
506 } <<
507 475
508 ////////////////////////////////////////////// << 476 switch (side)
509 // <<
510 // Algorithm for SurfaceNormal() following the <<
511 // for points not on the surface <<
512 <<
513 G4ThreeVector G4Para::ApproxSurfaceNormal( con <<
514 { <<
515 G4double dist = -DBL_MAX; <<
516 G4int iside = 0; <<
517 for (G4int i=0; i<4; ++i) <<
518 { 477 {
519 G4double d = fPlanes[i].a*p.x() + << 478 case kNX:
520 fPlanes[i].b*p.y() + << 479 tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
521 fPlanes[i].c*p.z() + fPlanes[ << 480 if (xshift<0)
522 if (d > dist) { dist = d; iside = i; } << 481 {
>> 482 cosntheta=-1/sqrt(1+tntheta*tntheta);
>> 483 }
>> 484 else
>> 485 {
>> 486 cosntheta=1/sqrt(1+tntheta*tntheta);
>> 487 }
>> 488 norm=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
>> 489 break;
>> 490 case kNY:
>> 491 if (newpy<0)
>> 492 {
>> 493 ycomp=-1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 494 }
>> 495 else
>> 496 {
>> 497 ycomp=1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 498 }
>> 499 norm=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
>> 500 break;
>> 501 case kNZ: // Closest to Z
>> 502 if (p.z()>=0)
>> 503 {
>> 504 norm=G4ThreeVector(0,0,1);
>> 505 }
>> 506 else
>> 507 {
>> 508 norm=G4ThreeVector(0,0,-1);
>> 509 }
>> 510 break;
523 } 511 }
524 << 512 return norm;
525 G4double distz = std::abs(p.z()) - fDz; <<
526 if (dist > distz) <<
527 return { fPlanes[iside].a, fPlanes[iside]. <<
528 else <<
529 return { 0, 0, (G4double)((p.z() < 0) ? -1 <<
530 } 513 }
531 514
532 ////////////////////////////////////////////// << 515 //////////////////////////////////////////////////////////////////////////////
533 // 516 //
534 // Calculate distance to shape from outside 517 // Calculate distance to shape from outside
535 // - return kInfinity if no intersection << 518 // - return kInfinity if no intersection
536 << 519 //
537 G4double G4Para::DistanceToIn(const G4ThreeVec << 520 // ALGORITHM:
538 const G4ThreeVec << 521 // For each component, calculate pair of minimum and maximum intersection
539 { << 522 // values for which the particle is in the extent of the shape
540 // Z intersections << 523 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
>> 524 // - Z plane intersectin uses tolerance
>> 525 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
>> 526 // (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
>> 527 // cases)
>> 528 // - Note: XZ and YZ planes each divide space into four regions,
>> 529 // characterised by ss1 ss2
>> 530
>> 531 G4double G4Para::DistanceToIn( const G4ThreeVector& p,
>> 532 const G4ThreeVector& v ) const
>> 533 {
>> 534 G4double snxt; // snxt = default return value
>> 535 G4double smin,smax;
>> 536 G4double tmin,tmax;
>> 537 G4double yt,vy,xt,vx;
>> 538 G4double max;
541 // 539 //
542 if ((std::abs(p.z()) - fDz) >= -halfCarToler << 540 // Z Intersection range
543 return kInfinity; <<
544 G4double invz = (-v.z() == 0) ? DBL_MAX : -1 <<
545 G4double dz = (invz < 0) ? fDz : -fDz; <<
546 G4double tzmin = (p.z() + dz)*invz; <<
547 G4double tzmax = (p.z() - dz)*invz; <<
548 <<
549 // Y intersections <<
550 // 541 //
551 G4double tmin0 = tzmin, tmax0 = tzmax; << 542 if (v.z()>0)
552 G4double cos0 = fPlanes[0].b*v.y() + fPlanes <<
553 G4double disy = fPlanes[0].b*p.y() + fPlanes <<
554 G4double dis0 = fPlanes[0].d + disy; <<
555 if (dis0 >= -halfCarTolerance) <<
556 { 543 {
557 if (cos0 >= 0) return kInfinity; << 544 max=fDz-p.z();
558 G4double tmp = -dis0/cos0; << 545 if (max>kCarTolerance*0.5)
559 if (tmin0 < tmp) tmin0 = tmp; << 546 {
>> 547 smax=max/v.z();
>> 548 smin=(-fDz-p.z())/v.z();
>> 549 }
>> 550 else
>> 551 {
>> 552 return snxt=kInfinity;
>> 553 }
>> 554 }
>> 555 else if (v.z()<0)
>> 556 {
>> 557 max=-fDz-p.z();
>> 558 if (max<-kCarTolerance*0.5)
>> 559 {
>> 560 smax=max/v.z();
>> 561 smin=(fDz-p.z())/v.z();
>> 562 }
>> 563 else
>> 564 {
>> 565 return snxt=kInfinity;
>> 566 }
560 } 567 }
561 else if (cos0 > 0) << 568 else
562 { 569 {
563 G4double tmp = -dis0/cos0; << 570 if (fabs(p.z())<=fDz) // Inside
564 if (tmax0 > tmp) tmax0 = tmp; << 571 {
>> 572 smin=0;
>> 573 smax=kInfinity;
>> 574 }
>> 575 else
>> 576 {
>> 577 return snxt=kInfinity;
>> 578 }
565 } 579 }
>> 580
>> 581 //
>> 582 // Y G4Parallel planes intersection
>> 583 //
566 584
567 G4double tmin1 = tmin0, tmax1 = tmax0; << 585 yt=p.y()-fTthetaSphi*p.z();
568 G4double cos1 = -cos0; << 586 vy=v.y()-fTthetaSphi*v.z();
569 G4double dis1 = fPlanes[1].d - disy; << 587
570 if (dis1 >= -halfCarTolerance) << 588 if (vy>0)
571 { 589 {
572 if (cos1 >= 0) return kInfinity; << 590 max=fDy-yt;
573 G4double tmp = -dis1/cos1; << 591 if (max>kCarTolerance*0.5)
574 if (tmin1 < tmp) tmin1 = tmp; << 592 {
>> 593 tmax=max/vy;
>> 594 tmin=(-fDy-yt)/vy;
>> 595 }
>> 596 else
>> 597 {
>> 598 return snxt=kInfinity;
>> 599 }
575 } 600 }
576 else if (cos1 > 0) << 601 else if (vy<0)
577 { 602 {
578 G4double tmp = -dis1/cos1; << 603 max=-fDy-yt;
579 if (tmax1 > tmp) tmax1 = tmp; << 604 if (max<-kCarTolerance*0.5)
>> 605 {
>> 606 tmax=max/vy;
>> 607 tmin=(fDy-yt)/vy;
>> 608 }
>> 609 else
>> 610 {
>> 611 return snxt=kInfinity;
>> 612 }
>> 613 }
>> 614 else
>> 615 {
>> 616 if (fabs(yt)<=fDy)
>> 617 {
>> 618 tmin=0;
>> 619 tmax=kInfinity;
>> 620 }
>> 621 else
>> 622 {
>> 623 return snxt=kInfinity;
>> 624 }
580 } 625 }
581 626
582 // X intersections << 627 // Re-Calc valid intersection range
583 // 628 //
584 G4double tmin2 = tmin1, tmax2 = tmax1; << 629 if (tmin>smin) smin=tmin;
585 G4double cos2 = fPlanes[2].a*v.x() + fPlanes << 630 if (tmax<smax) smax=tmax;
586 G4double disx = fPlanes[2].a*p.x() + fPlanes << 631 if (smax<=smin)
587 G4double dis2 = fPlanes[2].d + disx; <<
588 if (dis2 >= -halfCarTolerance) <<
589 { 632 {
590 if (cos2 >= 0) return kInfinity; << 633 return snxt=kInfinity;
591 G4double tmp = -dis2/cos2; <<
592 if (tmin2 < tmp) tmin2 = tmp; <<
593 } 634 }
594 else if (cos2 > 0) << 635 else
595 { 636 {
596 G4double tmp = -dis2/cos2; << 637 //
597 if (tmax2 > tmp) tmax2 = tmp; << 638 // X G4Parallel planes intersection
>> 639 //
>> 640 xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
>> 641 vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
>> 642 if (vx>0)
>> 643 {
>> 644 max=fDx-xt;
>> 645 if (max>kCarTolerance*0.5)
>> 646 {
>> 647 tmax=max/vx;
>> 648 tmin=(-fDx-xt)/vx;
>> 649 }
>> 650 else
>> 651 {
>> 652 return snxt=kInfinity;
>> 653 }
>> 654 }
>> 655 else if (vx<0)
>> 656 {
>> 657 max=-fDx-xt;
>> 658 if (max<-kCarTolerance*0.5)
>> 659 {
>> 660 tmax=max/vx;
>> 661 tmin=(fDx-xt)/vx;
>> 662 }
>> 663 else
>> 664 {
>> 665 return snxt=kInfinity;
>> 666 }
>> 667 }
>> 668 else
>> 669 {
>> 670 if (fabs(xt)<=fDx)
>> 671 {
>> 672 tmin=0;
>> 673 tmax=kInfinity;
>> 674 }
>> 675 else
>> 676 {
>> 677 return snxt=kInfinity;
>> 678 }
>> 679 }
>> 680 if (tmin>smin) smin=tmin;
>> 681 if (tmax<smax) smax=tmax;
598 } 682 }
599 683
600 G4double tmin3 = tmin2, tmax3 = tmax2; << 684 if (smax>0&&smin<smax)
601 G4double cos3 = -cos2; <<
602 G4double dis3 = fPlanes[3].d - disx; <<
603 if (dis3 >= -halfCarTolerance) <<
604 { 685 {
605 if (cos3 >= 0) return kInfinity; << 686 if (smin>0)
606 G4double tmp = -dis3/cos3; << 687 {
607 if (tmin3 < tmp) tmin3 = tmp; << 688 snxt=smin;
>> 689 }
>> 690 else
>> 691 {
>> 692 snxt=0;
>> 693 }
608 } 694 }
609 else if (cos3 > 0) << 695 else
610 { 696 {
611 G4double tmp = -dis3/cos3; << 697 snxt=kInfinity;
612 if (tmax3 > tmp) tmax3 = tmp; <<
613 } 698 }
614 << 699 return snxt;
615 // Find distance <<
616 // <<
617 G4double tmin = tmin3, tmax = tmax3; <<
618 if (tmax <= tmin + halfCarTolerance) return <<
619 return (tmin < halfCarTolerance ) ? 0. : tmi <<
620 } 700 }
621 701
622 ////////////////////////////////////////////// << 702 ////////////////////////////////////////////////////////////////////////////
623 // 703 //
624 // Calculate exact shortest distance to any bo 704 // Calculate exact shortest distance to any boundary from outside
625 // - returns 0 is point inside << 705 // - Returns 0 is point inside
626 706
627 G4double G4Para::DistanceToIn( const G4ThreeVe 707 G4double G4Para::DistanceToIn( const G4ThreeVector& p ) const
628 { 708 {
629 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. << 709 G4double safe;
630 G4double dx = std::abs(xx) + fPlanes[2].d; << 710 G4double distz1,distz2,disty1,disty2,distx1,distx2;
>> 711 G4double trany,cosy,tranx,cosx;
>> 712
>> 713 // Z planes
>> 714 //
>> 715 distz1=p.z()-fDz;
>> 716 distz2=-fDz-p.z();
>> 717 if (distz1>distz2)
>> 718 {
>> 719 safe=distz1;
>> 720 }
>> 721 else
>> 722 {
>> 723 safe=distz2;
>> 724 }
631 725
632 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. << 726 trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
633 G4double dy = std::abs(yy) + fPlanes[0].d; <<
634 G4double dxy = std::max(dx,dy); <<
635 727
636 G4double dz = std::abs(p.z())-fDz; << 728 // Transformed x into `box' system
637 G4double dist = std::max(dxy,dz); << 729 //
>> 730 cosy=1.0/sqrt(1.0+fTthetaSphi*fTthetaSphi);
>> 731 disty1=(trany-fDy)*cosy;
>> 732 disty2=(-fDy-trany)*cosy;
>> 733
>> 734 if (disty1>safe) safe=disty1;
>> 735 if (disty2>safe) safe=disty2;
638 736
639 return (dist > 0) ? dist : 0.; << 737 tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
>> 738 cosx=1.0/sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
>> 739 distx1=(tranx-fDx)*cosx;
>> 740 distx2=(-fDx-tranx)*cosx;
>> 741
>> 742 if (distx1>safe) safe=distx1;
>> 743 if (distx2>safe) safe=distx2;
>> 744
>> 745 if (safe<0) safe=0;
>> 746 return safe;
640 } 747 }
641 748
642 ////////////////////////////////////////////// 749 //////////////////////////////////////////////////////////////////////////
643 // 750 //
644 // Calculate distance to surface of shape from << 751 // Calculate distance to surface of shape from inside
645 // find normal at exit point << 752 // Calculate distance to x/y/z planes - smallest is exiting distance
646 // - when leaving the surface, return 0 <<
647 753
648 G4double G4Para::DistanceToOut(const G4ThreeVe 754 G4double G4Para::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
649 const G4bool ca 755 const G4bool calcNorm,
650 G4bool* v << 756 G4bool *validNorm, G4ThreeVector *n) const
651 { 757 {
652 // Z intersections << 758 ESide side = kUndef;
>> 759 G4double snxt; // snxt = return value
>> 760 G4double max,tmax;
>> 761 G4double yt,vy,xt,vx;
>> 762
>> 763 G4double ycomp,calpha,salpha,tntheta,cosntheta;
>> 764
>> 765 //
>> 766 // Z Intersections
653 // 767 //
654 if ((std::abs(p.z()) - fDz) >= -halfCarToler << 768
>> 769 if (v.z()>0)
655 { 770 {
656 if (calcNorm) << 771 max=fDz-p.z();
>> 772 if (max>kCarTolerance*0.5)
657 { 773 {
658 *validNorm = true; << 774 snxt=max/v.z();
659 n->set(0, 0, (p.z() < 0) ? -1 : 1); << 775 side=kPZ;
>> 776 }
>> 777 else
>> 778 {
>> 779 if (calcNorm)
>> 780 {
>> 781 *validNorm=true;
>> 782 *n=G4ThreeVector(0,0,1);
>> 783 }
>> 784 return snxt=0;
660 } 785 }
661 return 0.; <<
662 } 786 }
663 G4double vz = v.z(); << 787 else if (v.z()<0)
664 G4double tmax = (vz == 0) ? DBL_MAX : (std:: <<
665 G4int iside = (vz < 0) ? -4 : -2; // little <<
666 <<
667 // Y intersections <<
668 // <<
669 G4double cos0 = fPlanes[0].b*v.y() + fPlanes <<
670 if (cos0 > 0) <<
671 { 788 {
672 G4double dis0 = fPlanes[0].b*p.y() + fPlan << 789 max=-fDz-p.z();
673 if (dis0 >= -halfCarTolerance) << 790 if (max<-kCarTolerance*0.5)
>> 791 {
>> 792 snxt=max/v.z();
>> 793 side=kMZ;
>> 794 }
>> 795 else
674 { 796 {
675 if (calcNorm) 797 if (calcNorm)
676 { 798 {
677 *validNorm = true; << 799 *validNorm=true;
678 n->set(0, fPlanes[0].b, fPlanes[0].c); << 800 *n=G4ThreeVector(0,0,-1);
679 } 801 }
680 return 0.; << 802 return snxt=0;
681 } 803 }
682 G4double tmp = -dis0/cos0; <<
683 if (tmax > tmp) { tmax = tmp; iside = 0; } <<
684 } 804 }
>> 805 else
>> 806 {
>> 807 snxt=kInfinity;
>> 808 }
>> 809
>> 810 //
>> 811 // Y plane intersection
>> 812 //
>> 813
>> 814 yt=p.y()-fTthetaSphi*p.z();
>> 815 vy=v.y()-fTthetaSphi*v.z();
685 816
686 G4double cos1 = -cos0; << 817 if (vy>0)
687 if (cos1 > 0) <<
688 { 818 {
689 G4double dis1 = fPlanes[1].b*p.y() + fPlan << 819 max=fDy-yt;
690 if (dis1 >= -halfCarTolerance) << 820 if (max>kCarTolerance*0.5)
>> 821 {
>> 822 tmax=max/vy;
>> 823 if (tmax<snxt)
>> 824 {
>> 825 snxt=tmax;
>> 826 side=kPY;
>> 827 }
>> 828 }
>> 829 else
691 { 830 {
692 if (calcNorm) 831 if (calcNorm)
>> 832 {
>> 833 *validNorm=true; // Leaving via plus Y
>> 834 ycomp=1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 835 *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
>> 836 }
>> 837 return snxt=0;
>> 838 }
>> 839 }
>> 840 else if (vy<0)
>> 841 {
>> 842 max=-fDy-yt;
>> 843 if (max<-kCarTolerance*0.5)
>> 844 {
>> 845 tmax=max/vy;
>> 846 if (tmax<snxt)
693 { 847 {
694 *validNorm = true; << 848 snxt=tmax;
695 n->set(0, fPlanes[1].b, fPlanes[1].c); << 849 side=kMY;
696 } 850 }
697 return 0.; <<
698 } 851 }
699 G4double tmp = -dis1/cos1; << 852 else
700 if (tmax > tmp) { tmax = tmp; iside = 1; } << 853 {
>> 854 if (calcNorm)
>> 855 {
>> 856 *validNorm=true; // Leaving via minus Y
>> 857 ycomp=-1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 858 *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
>> 859 }
>> 860 return snxt=0;
>> 861 }
701 } 862 }
702 863
703 // X intersections <<
704 // 864 //
705 G4double cos2 = fPlanes[2].a*v.x() + fPlanes << 865 // X plane intersection
706 if (cos2 > 0) << 866 //
>> 867
>> 868 xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
>> 869 vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
>> 870 if (vx>0)
707 { 871 {
708 G4double dis2 = fPlanes[2].a*p.x()+fPlanes << 872 max=fDx-xt;
709 if (dis2 >= -halfCarTolerance) << 873 if (max>kCarTolerance*0.5)
>> 874 {
>> 875 tmax=max/vx;
>> 876 if (tmax<snxt)
>> 877 {
>> 878 snxt=tmax;
>> 879 side=kPX;
>> 880 }
>> 881 }
>> 882 else
710 { 883 {
711 if (calcNorm) 884 if (calcNorm)
712 { 885 {
713 *validNorm = true; << 886 *validNorm=true; // Leaving via plus X
714 n->set(fPlanes[2].a, fPlanes[2].b, fP << 887 calpha=1/sqrt(1+fTalpha*fTalpha);
>> 888 if (fTalpha)
>> 889 {
>> 890 salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
>> 891 }
>> 892 else
>> 893 {
>> 894 salpha=0;
>> 895 }
>> 896 tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
>> 897 cosntheta=1/sqrt(1+tntheta*tntheta);
>> 898 *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
715 } 899 }
716 return 0.; << 900 return snxt=0;
717 } 901 }
718 G4double tmp = -dis2/cos2; <<
719 if (tmax > tmp) { tmax = tmp; iside = 2; } <<
720 } 902 }
721 << 903 else if (vx<0)
722 G4double cos3 = -cos2; <<
723 if (cos3 > 0) <<
724 { 904 {
725 G4double dis3 = fPlanes[3].a*p.x()+fPlanes << 905 max=-fDx-xt;
726 if (dis3 >= -halfCarTolerance) << 906 if (max<-kCarTolerance*0.5)
>> 907 {
>> 908 tmax=max/vx;
>> 909 if (tmax<snxt)
>> 910 {
>> 911 snxt=tmax;
>> 912 side=kMX;
>> 913 }
>> 914 }
>> 915 else
727 { 916 {
728 if (calcNorm) 917 if (calcNorm)
729 { 918 {
730 *validNorm = true; << 919 *validNorm=true; // Leaving via minus X
731 n->set(fPlanes[3].a, fPlanes[3].b, fP << 920 calpha=1/sqrt(1+fTalpha*fTalpha);
>> 921 if (fTalpha)
>> 922 {
>> 923 salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
>> 924 }
>> 925 else
>> 926 {
>> 927 salpha=0;
>> 928 }
>> 929 tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
>> 930 cosntheta=-1/sqrt(1+tntheta*tntheta);
>> 931 *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
>> 932 return snxt=0;
732 } 933 }
733 return 0.; <<
734 } 934 }
735 G4double tmp = -dis3/cos3; <<
736 if (tmax > tmp) { tmax = tmp; iside = 3; } <<
737 } 935 }
738 936
739 // Set normal, if required, and return dista << 937 if (calcNorm)
740 // <<
741 if (calcNorm) <<
742 { 938 {
743 *validNorm = true; << 939 *validNorm=true;
744 if (iside < 0) << 940 switch (side)
745 n->set(0, 0, iside + 3); // (-4+3)=-1, ( << 941 {
746 else << 942 case kMZ:
747 n->set(fPlanes[iside].a, fPlanes[iside]. << 943 *n=G4ThreeVector(0,0,-1);
>> 944 break;
>> 945 case kPZ:
>> 946 *n=G4ThreeVector(0,0,1);
>> 947 break;
>> 948 case kMY:
>> 949 ycomp=-1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 950 *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
>> 951 break;
>> 952 case kPY:
>> 953 ycomp=1/sqrt(1+fTthetaSphi*fTthetaSphi);
>> 954 *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
>> 955 break;
>> 956 case kMX:
>> 957 calpha=1/sqrt(1+fTalpha*fTalpha);
>> 958 if (fTalpha)
>> 959 {
>> 960 salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
>> 961 }
>> 962 else
>> 963 {
>> 964 salpha=0;
>> 965 }
>> 966 tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
>> 967 cosntheta=-1/sqrt(1+tntheta*tntheta);
>> 968 *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
>> 969 break;
>> 970 case kPX:
>> 971 calpha=1/sqrt(1+fTalpha*fTalpha);
>> 972 if (fTalpha)
>> 973 {
>> 974 salpha=-calpha/fTalpha; // NOTE: actually use MINUS sin(alpha)
>> 975 }
>> 976 else
>> 977 {
>> 978 salpha=0;
>> 979 }
>> 980 tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
>> 981 cosntheta=1/sqrt(1+tntheta*tntheta);
>> 982 *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
>> 983 break;
>> 984 default:
>> 985 DumpInfo();
>> 986 G4Exception("G4Para::DistanceToOut() - Invalid enum");
>> 987 break;
>> 988 }
748 } 989 }
749 return tmax; << 990 return snxt;
750 } 991 }
751 992
752 ////////////////////////////////////////////// << 993 /////////////////////////////////////////////////////////////////////////////
753 // 994 //
754 // Calculate exact shortest distance to any bo 995 // Calculate exact shortest distance to any boundary from inside
755 // - returns 0 is point outside << 996 // - Returns 0 is point outside
756 997
757 G4double G4Para::DistanceToOut( const G4ThreeV 998 G4double G4Para::DistanceToOut( const G4ThreeVector& p ) const
758 { 999 {
>> 1000 G4double safe;
>> 1001 G4double distz1,distz2,disty1,disty2,distx1,distx2;
>> 1002 G4double trany,cosy,tranx,cosx;
>> 1003
759 #ifdef G4CSGDEBUG 1004 #ifdef G4CSGDEBUG
760 if( Inside(p) == kOutside ) 1005 if( Inside(p) == kOutside )
761 { 1006 {
762 std::ostringstream message; << 1007 G4cout.precision(16) ;
763 G4int oldprc = message.precision(16); << 1008 G4cout << G4endl ;
764 message << "Point p is outside (!?) of sol << 1009 DumpInfo();
765 message << "Position:\n"; << 1010 G4cout << "Position:" << G4endl << G4endl ;
766 message << " p.x() = " << p.x()/mm << " << 1011 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
767 message << " p.y() = " << p.y()/mm << " << 1012 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
768 message << " p.z() = " << p.z()/mm << " << 1013 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
769 G4cout.precision(oldprc) ; << 1014 G4cout << "G4Para::DistanceToOut(p) - point p is outside ?!" << G4endl ;
770 G4Exception("G4Para::DistanceToOut(p)", "G << 1015 G4cerr << "G4Para::DistanceToOut(p) - point p is outside ?!" << G4endl ;
771 JustWarning, message ); << 1016 }
772 DumpInfo(); <<
773 } <<
774 #endif 1017 #endif
775 G4double xx = fPlanes[2].a*p.x()+fPlanes[2]. <<
776 G4double dx = std::abs(xx) + fPlanes[2].d; <<
777 <<
778 G4double yy = fPlanes[0].b*p.y()+fPlanes[0]. <<
779 G4double dy = std::abs(yy) + fPlanes[0].d; <<
780 G4double dxy = std::max(dx,dy); <<
781 1018
782 G4double dz = std::abs(p.z())-fDz; << 1019 // Z planes
783 G4double dist = std::max(dxy,dz); << 1020 //
>> 1021 distz1=fDz-p.z();
>> 1022 distz2=fDz+p.z();
>> 1023 if (distz1<distz2)
>> 1024 {
>> 1025 safe=distz1;
>> 1026 }
>> 1027 else
>> 1028 {
>> 1029 safe=distz2;
>> 1030 }
784 1031
785 return (dist < 0) ? -dist : 0.; << 1032 trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
786 } <<
787 1033
788 ////////////////////////////////////////////// << 1034 // Transformed x into `box' system
789 // << 1035 //
790 // GetEntityType << 1036 cosy=1.0/sqrt(1.0+fTthetaSphi*fTthetaSphi);
>> 1037 disty1=(fDy-trany)*cosy;
>> 1038 disty2=(fDy+trany)*cosy;
>> 1039
>> 1040 if (disty1<safe) safe=disty1;
>> 1041 if (disty2<safe) safe=disty2;
791 1042
792 G4GeometryType G4Para::GetEntityType() const << 1043 tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
793 { << 1044 cosx=1.0/sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
794 return {"G4Para"}; << 1045 distx1=(fDx-tranx)*cosx;
>> 1046 distx2=(fDx+tranx)*cosx;
>> 1047
>> 1048 if (distx1<safe) safe=distx1;
>> 1049 if (distx2<safe) safe=distx2;
>> 1050
>> 1051 if (safe<0) safe=0;
>> 1052 return safe;
795 } 1053 }
796 1054
797 ////////////////////////////////////////////// << 1055 ////////////////////////////////////////////////////////////////////////////////
798 // 1056 //
799 // IsFaceted << 1057 // Create a List containing the transformed vertices
800 << 1058 // Ordering [0-3] -fDz cross section
801 G4bool G4Para::IsFaceted() const << 1059 // [4-7] +fDz cross section such that [0] is below [4],
802 { << 1060 // [1] below [5] etc.
803 return true; << 1061 // Note:
>> 1062 // Caller has deletion resposibility
>> 1063
>> 1064 G4ThreeVectorList*
>> 1065 G4Para::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
>> 1066 {
>> 1067 G4ThreeVectorList *vertices;
>> 1068 vertices=new G4ThreeVectorList();
>> 1069 vertices->reserve(8);
>> 1070 if (vertices)
>> 1071 {
>> 1072 G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy*fTalpha-fDx,
>> 1073 -fDz*fTthetaSphi-fDy, -fDz);
>> 1074 G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy*fTalpha+fDx,
>> 1075 -fDz*fTthetaSphi-fDy, -fDz);
>> 1076 G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy*fTalpha-fDx,
>> 1077 -fDz*fTthetaSphi+fDy, -fDz);
>> 1078 G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy*fTalpha+fDx,
>> 1079 -fDz*fTthetaSphi+fDy, -fDz);
>> 1080 G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy*fTalpha-fDx,
>> 1081 +fDz*fTthetaSphi-fDy, +fDz);
>> 1082 G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy*fTalpha+fDx,
>> 1083 +fDz*fTthetaSphi-fDy, +fDz);
>> 1084 G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy*fTalpha-fDx,
>> 1085 +fDz*fTthetaSphi+fDy, +fDz);
>> 1086 G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy*fTalpha+fDx,
>> 1087 +fDz*fTthetaSphi+fDy, +fDz);
>> 1088
>> 1089 vertices->push_back(pTransform.TransformPoint(vertex0));
>> 1090 vertices->push_back(pTransform.TransformPoint(vertex1));
>> 1091 vertices->push_back(pTransform.TransformPoint(vertex2));
>> 1092 vertices->push_back(pTransform.TransformPoint(vertex3));
>> 1093 vertices->push_back(pTransform.TransformPoint(vertex4));
>> 1094 vertices->push_back(pTransform.TransformPoint(vertex5));
>> 1095 vertices->push_back(pTransform.TransformPoint(vertex6));
>> 1096 vertices->push_back(pTransform.TransformPoint(vertex7));
>> 1097 }
>> 1098 else
>> 1099 {
>> 1100 DumpInfo();
>> 1101 G4Exception("G4Para::CreateRotatedVertices() - Out of memory !");
>> 1102 }
>> 1103 return vertices;
804 } 1104 }
805 1105
806 ////////////////////////////////////////////// 1106 //////////////////////////////////////////////////////////////////////////
807 // 1107 //
808 // Make a clone of the object << 1108 // GetEntityType
809 // << 1109
810 G4VSolid* G4Para::Clone() const << 1110 G4GeometryType G4Para::GetEntityType() const
811 { 1111 {
812 return new G4Para(*this); << 1112 return G4String("G4Para");
813 } 1113 }
814 1114
815 ////////////////////////////////////////////// 1115 //////////////////////////////////////////////////////////////////////////
816 // 1116 //
817 // Stream object contents to an output stream 1117 // Stream object contents to an output stream
818 1118
819 std::ostream& G4Para::StreamInfo( std::ostream 1119 std::ostream& G4Para::StreamInfo( std::ostream& os ) const
820 { 1120 {
821 G4double alpha = std::atan(fTalpha); <<
822 G4double theta = std::atan(std::sqrt(fTtheta <<
823 fTtheta <<
824 G4double phi = std::atan2(fTthetaSphi,fTth <<
825 <<
826 G4long oldprc = os.precision(16); <<
827 os << "------------------------------------- 1121 os << "-----------------------------------------------------------\n"
828 << " *** Dump for solid - " << GetName 1122 << " *** Dump for solid - " << GetName() << " ***\n"
829 << " ================================= 1123 << " ===================================================\n"
830 << " Solid type: G4Para\n" 1124 << " Solid type: G4Para\n"
831 << " Parameters:\n" << 1125 << " Parameters: \n"
832 << " half length X: " << fDx/mm << " m << 1126 << " half length X: " << fDx/mm << " mm \n"
833 << " half length Y: " << fDy/mm << " m << 1127 << " half length Y: " << fDy/mm << " mm \n"
834 << " half length Z: " << fDz/mm << " m << 1128 << " half length Z: " << fDz/mm << " mm \n"
835 << " alpha: " << alpha/degree << "degr << 1129 << " tan(alpha) : " << fTalpha/degree << " degrees \n"
836 << " theta: " << theta/degree << "degr << 1130 << " tan(theta)*cos(phi): " << fTthetaCphi/degree << " degrees \n"
837 << " phi: " << phi/degree << "degrees\ << 1131 << " tan(theta)*sin(phi): " << fTthetaSphi/degree << " degrees \n"
838 << "------------------------------------- 1132 << "-----------------------------------------------------------\n";
839 os.precision(oldprc); <<
840 1133
841 return os; 1134 return os;
842 } 1135 }
843 1136
844 ////////////////////////////////////////////// << 1137 ////////////////////////////////////////////////////////////////////////////
845 // <<
846 // Return a point randomly and uniformly selec <<
847 <<
848 G4ThreeVector G4Para::GetPointOnSurface() cons <<
849 { <<
850 G4double DyTalpha = fDy*fTalpha; <<
851 G4double DzTthetaSphi = fDz*fTthetaSphi; <<
852 G4double DzTthetaCphi = fDz*fTthetaCphi; <<
853 <<
854 // Set vertices <<
855 // <<
856 G4ThreeVector pt[8]; <<
857 pt[0].set(-DzTthetaCphi-DyTalpha-fDx, -DzTth <<
858 pt[1].set(-DzTthetaCphi-DyTalpha+fDx, -DzTth <<
859 pt[2].set(-DzTthetaCphi+DyTalpha-fDx, -DzTth <<
860 pt[3].set(-DzTthetaCphi+DyTalpha+fDx, -DzTth <<
861 pt[4].set( DzTthetaCphi-DyTalpha-fDx, DzTth <<
862 pt[5].set( DzTthetaCphi-DyTalpha+fDx, DzTth <<
863 pt[6].set( DzTthetaCphi+DyTalpha-fDx, DzTth <<
864 pt[7].set( DzTthetaCphi+DyTalpha+fDx, DzTth <<
865 <<
866 // Set areas (-Z, -Y, +Y, -X, +X, +Z) <<
867 // <<
868 G4ThreeVector vx(fDx, 0, 0); <<
869 G4ThreeVector vy(DyTalpha, fDy, 0); <<
870 G4ThreeVector vz(DzTthetaCphi, DzTthetaSphi, <<
871 <<
872 G4double sxy = fDx*fDy; // (vx.cross(vy)).ma <<
873 G4double sxz = (vx.cross(vz)).mag(); <<
874 G4double syz = (vy.cross(vz)).mag(); <<
875 <<
876 G4double sface[6] = { sxy, syz, syz, sxz, sx <<
877 for (G4int i=1; i<6; ++i) { sface[i] += sfac <<
878 <<
879 // Select face <<
880 // <<
881 G4double select = sface[5]*G4UniformRand(); <<
882 G4int k = 5; <<
883 if (select <= sface[4]) k = 4; <<
884 if (select <= sface[3]) k = 3; <<
885 if (select <= sface[2]) k = 2; <<
886 if (select <= sface[1]) k = 1; <<
887 if (select <= sface[0]) k = 0; <<
888 <<
889 // Generate point <<
890 // <<
891 G4int ip[6][3] = {{0,1,2}, {0,4,1}, {2,3,6}, <<
892 G4double u = G4UniformRand(); <<
893 G4double v = G4UniformRand(); <<
894 return (1.-u-v)*pt[ip[k][0]] + u*pt[ip[k][1] <<
895 } <<
896 <<
897 ////////////////////////////////////////////// <<
898 // 1138 //
899 // Methods for visualisation 1139 // Methods for visualisation
900 1140
901 void G4Para::DescribeYourselfTo ( G4VGraphicsS 1141 void G4Para::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
902 { 1142 {
903 scene.AddSolid (*this); << 1143 scene.AddThis (*this);
904 } 1144 }
905 1145
906 G4Polyhedron* G4Para::CreatePolyhedron () cons 1146 G4Polyhedron* G4Para::CreatePolyhedron () const
907 { 1147 {
908 G4double phi = std::atan2(fTthetaSphi, fTthe << 1148 G4double phi = atan2(fTthetaSphi, fTthetaCphi);
909 G4double alpha = std::atan(fTalpha); << 1149 G4double alpha = atan(fTalpha);
910 G4double theta = std::atan(std::sqrt(fTtheta << 1150 G4double theta = atan(sqrt(fTthetaCphi*fTthetaCphi
911 fTtheta << 1151 +fTthetaSphi*fTthetaSphi));
912 1152
913 return new G4PolyhedronPara(fDx, fDy, fDz, a 1153 return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi);
914 } 1154 }
915 #endif << 1155
>> 1156 G4NURBS* G4Para::CreateNURBS () const
>> 1157 {
>> 1158 // return new G4NURBSbox (fDx, fDy, fDz);
>> 1159 return 0 ;
>> 1160 }
>> 1161
>> 1162 //
>> 1163 //
>> 1164 /////////////////////////// End of G4Para.cc ///////////////////////////
916 1165