Context Navigation

source: trunk/source/geometry/solids/specific/src/G4VTwistedFaceted.cc@ 1127

Last change on this file since 1127 was 1058, checked in by garnier, 17 years ago
file release beta
File size: 39.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: