Context Navigation

source: trunk/source/geometry/solids/specific/src/G4TwistTubsFlatSide.cc @ 1202

Last change on this file since 1202 was 1058, checked in by garnier, 15 years ago
file release beta
File size: 19.8 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	//
27	// $Id: G4TwistTubsFlatSide.cc,v 1.7 2007/05/23 09:31:02 gcosmo Exp $
28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	//
30	//
31	// --------------------------------------------------------------------
32	// GEANT 4 class source file
33	//
34	//
35	// G4TwistTubsFlatSide.cc
36	//
37	// Author:
38	// 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
39	//
40	// History:
41	// 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
42	// from original version in Jupiter-2.5.02 application.
43	// --------------------------------------------------------------------
44
45	#include "G4TwistTubsFlatSide.hh"
46	#include "G4GeometryTolerance.hh"
47
48	//=====================================================================
49	//* constructors ------------------------------------------------------
50
51	G4TwistTubsFlatSide::G4TwistTubsFlatSide(const G4String &name,
52	const G4RotationMatrix &rot,
53	const G4ThreeVector &tlate,
54	const G4ThreeVector &n,
55	const EAxis axis0 ,
56	const EAxis axis1 ,
57	G4double axis0min,
58	G4double axis1min,
59	G4double axis0max,
60	G4double axis1max )
61	: G4VTwistSurface(name, rot, tlate, 0, axis0, axis1,
62	axis0min, axis1min, axis0max, axis1max)
63	{
64	if (axis0 == kPhi && axis1 == kRho) {
65	G4Exception("G4TwistTubsFlatSide::G4TwistTubsFlatSide()", "InvalidSetup",
66	FatalException, "Should swap axis0 and axis1!");
67	}
68
69	G4ThreeVector normal = rot.inverse()*n;
70	fCurrentNormal.normal = normal.unit(); // in local coordinate system
71	fIsValidNorm = true;
72
73	SetCorners();
74	SetBoundaries();
75
76	fSurfaceArea = 1 ; // not yet implemented. This is NOT a problem for tracking
77
78	}
79
80
81
82	G4TwistTubsFlatSide::G4TwistTubsFlatSide( const G4String &name,
83	G4double EndInnerRadius[2],
84	G4double EndOuterRadius[2],
85	G4double DPhi,
86	G4double EndPhi[2],
87	G4double EndZ[2],
88	G4int handedness )
89	: G4VTwistSurface(name)
90	{
91	fHandedness = handedness; // +z = +ve, -z = -ve
92	fAxis[0] = kRho; // in local coordinate system
93	fAxis[1] = kPhi;
94	G4int i = (handedness < 0 ? 0 : 1);
95	fAxisMin[0] = EndInnerRadius[i]; // Inner-hype radius at z=0
96	fAxisMax[0] = EndOuterRadius[i]; // Outer-hype radius at z=0
97	fAxisMin[1] = -0.5*DPhi;
98	fAxisMax[1] = -fAxisMin[1];
99	fCurrentNormal.normal.set(0, 0, (fHandedness < 0 ? -1 : 1));
100	// Unit vector, in local coordinate system
101	fRot.rotateZ(EndPhi[i]);
102	fTrans.set(0, 0, EndZ[i]);
103	fIsValidNorm = true;
104
105	SetCorners();
106	SetBoundaries();
107
108	fSurfaceArea = 0.5DPhi (EndOuterRadius[i]*EndOuterRadius[i]
109	- EndInnerRadius[i]*EndInnerRadius[i] ) ;
110
111	}
112
113
114	//=====================================================================
115	//* Fake default constructor ------------------------------------------
116
117	G4TwistTubsFlatSide::G4TwistTubsFlatSide( __void__& a )
118	: G4VTwistSurface(a)
119	{
120	}
121
122
123	//=====================================================================
124	//* destructor --------------------------------------------------------
125
126	G4TwistTubsFlatSide::~G4TwistTubsFlatSide()
127	{
128	}
129
130	//=====================================================================
131	//* GetNormal ---------------------------------------------------------
132
133	G4ThreeVector G4TwistTubsFlatSide::GetNormal(const G4ThreeVector & /* xx */ ,
134	G4bool isGlobal)
135	{
136	if (isGlobal) {
137	return ComputeGlobalDirection(fCurrentNormal.normal);
138	} else {
139	return fCurrentNormal.normal;
140	}
141	}
142
143	//=====================================================================
144	//* DistanceToSurface(p, v) -------------------------------------------
145
146	G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector &gp,
147	const G4ThreeVector &gv,
148	G4ThreeVector gxx[],
149	G4double distance[],
150	G4int areacode[],
151	G4bool isvalid[],
152	EValidate validate)
153	{
154	fCurStatWithV.ResetfDone(validate, &gp, &gv);
155
156	if (fCurStatWithV.IsDone()) {
157	G4int i;
158	for (i=0; i<fCurStatWithV.GetNXX(); i++) {
159	gxx[i] = fCurStatWithV.GetXX(i);
160	distance[i] = fCurStatWithV.GetDistance(i);
161	areacode[i] = fCurStatWithV.GetAreacode(i);
162	isvalid[i] = fCurStatWithV.IsValid(i);
163	}
164	return fCurStatWithV.GetNXX();
165	} else {
166	// initialize
167	G4int i;
168	for (i=0; i<2; i++) {
169	distance[i] = kInfinity;
170	areacode[i] = sOutside;
171	isvalid[i] = false;
172	gxx[i].set(kInfinity, kInfinity, kInfinity);
173	}
174	}
175
176	G4ThreeVector p = ComputeLocalPoint(gp);
177	G4ThreeVector v = ComputeLocalDirection(gv);
178
179	//
180	// special case!
181	// if p is on surface, distance = 0.
182	//
183
184	if (std::fabs(p.z()) == 0.) { // if p is on the plane
185	distance[0] = 0;
186	G4ThreeVector xx = p;
187	gxx[0] = ComputeGlobalPoint(xx);
188
189	if (validate == kValidateWithTol) {
190	areacode[0] = GetAreaCode(xx);
191	if (!IsOutside(areacode[0])) {
192	isvalid[0] = true;
193	}
194	} else if (validate == kValidateWithoutTol) {
195	areacode[0] = GetAreaCode(xx, false);
196	if (IsInside(areacode[0])) {
197	isvalid[0] = true;
198	}
199	} else { // kDontValidate
200	areacode[0] = sInside;
201	isvalid[0] = true;
202	}
203
204	return 1;
205	}
206	//
207	// special case end
208	//
209
210	if (v.z() == 0) {
211
212	fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
213	isvalid[0], 0, validate, &gp, &gv);
214	return 0;
215	}
216
217	distance[0] = - (p.z() / v.z());
218
219	G4ThreeVector xx = p + distance[0]*v;
220	gxx[0] = ComputeGlobalPoint(xx);
221
222	if (validate == kValidateWithTol) {
223	areacode[0] = GetAreaCode(xx);
224	if (!IsOutside(areacode[0])) {
225	if (distance[0] >= 0) isvalid[0] = true;
226	}
227	} else if (validate == kValidateWithoutTol) {
228	areacode[0] = GetAreaCode(xx, false);
229	if (IsInside(areacode[0])) {
230	if (distance[0] >= 0) isvalid[0] = true;
231	}
232	} else { // kDontValidate
233	areacode[0] = sInside;
234	if (distance[0] >= 0) isvalid[0] = true;
235	}
236
237	fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
238	isvalid[0], 1, validate, &gp, &gv);
239
240	#ifdef G4TWISTDEBUG
241	G4cerr << "ERROR - G4TwistTubsFlatSide::DistanceToSurface(p,v)" << G4endl;
242	G4cerr << " Name : " << GetName() << G4endl;
243	G4cerr << " xx : " << xx << G4endl;
244	G4cerr << " gxx[0] : " << gxx[0] << G4endl;
245	G4cerr << " dist[0] : " << distance[0] << G4endl;
246	G4cerr << " areacode[0] : " << areacode[0] << G4endl;
247	G4cerr << " isvalid[0] : " << isvalid[0] << G4endl;
248	#endif
249	return 1;
250	}
251
252	//=====================================================================
253	//* DistanceToSurface(p) ----------------------------------------------
254
255	G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector &gp,
256	G4ThreeVector gxx[],
257	G4double distance[],
258	G4int areacode[])
259	{
260	// Calculate distance to plane in local coordinate,
261	// then return distance and global intersection points.
262	//
263
264	fCurStat.ResetfDone(kDontValidate, &gp);
265
266	if (fCurStat.IsDone()) {
267	G4int i;
268	for (i=0; i<fCurStat.GetNXX(); i++) {
269	gxx[i] = fCurStat.GetXX(i);
270	distance[i] = fCurStat.GetDistance(i);
271	areacode[i] = fCurStat.GetAreacode(i);
272	}
273	return fCurStat.GetNXX();
274	} else {
275	// initialize
276	G4int i;
277	for (i=0; i<2; i++) {
278	distance[i] = kInfinity;
279	areacode[i] = sOutside;
280	gxx[i].set(kInfinity, kInfinity, kInfinity);
281	}
282	}
283
284	G4ThreeVector p = ComputeLocalPoint(gp);
285	G4ThreeVector xx;
286
287	// The plane is placed on origin with making its normal
288	// parallel to z-axis.
289	if (std::fabs(p.z()) <= 0.5 * kCarTolerance) { // if p is on the plane, return 1
290	distance[0] = 0;
291	xx = p;
292	} else {
293	distance[0] = std::fabs(p.z());
294	xx.set(p.x(), p.y(), 0);
295	}
296
297	gxx[0] = ComputeGlobalPoint(xx);
298	areacode[0] = sInside;
299	G4bool isvalid = true;
300	fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
301	isvalid, 1, kDontValidate, &gp);
302	return 1;
303
304	}
305
306	//=====================================================================
307	//* GetAreaCode -------------------------------------------------------
308
309	G4int G4TwistTubsFlatSide::GetAreaCode(const G4ThreeVector &xx,
310	G4bool withTol)
311	{
312
313	static const G4double rtol
314	= 0.5*G4GeometryTolerance::GetInstance()->GetRadialTolerance();
315
316	G4int areacode = sInside;
317
318	if (fAxis[0] == kRho && fAxis[1] == kPhi) {
319	G4int rhoaxis = 0;
320	// G4int phiaxis = 0;
321
322	G4ThreeVector dphimin; // direction of phi-minimum boundary
323	G4ThreeVector dphimax; // direction of phi-maximum boundary
324	dphimin = GetCorner(sC0Max1Min);
325	dphimax = GetCorner(sC0Max1Max);
326
327	if (withTol) {
328
329	G4bool isoutside = false;
330
331	// test boundary of rho-axis
332
333	if (xx.getRho() <= fAxisMin[rhoaxis] + rtol) {
334
335	areacode \|= (sAxis0 & (sAxisRho \| sAxisMin)) \| sBoundary; // rho-min
336	if (xx.getRho() < fAxisMin[rhoaxis] - rtol) isoutside = true;
337
338	} else if (xx.getRho() >= fAxisMax[rhoaxis] - rtol) {
339
340	areacode \|= (sAxis0 & (sAxisRho \| sAxisMax)) \| sBoundary; // rho-max
341	if (xx.getRho() > fAxisMax[rhoaxis] + rtol) isoutside = true;
342
343	}
344
345	// test boundary of phi-axis
346
347	if (AmIOnLeftSide(xx, dphimin) >= 0) { // xx is on dphimin
348
349	areacode \|= (sAxis1 & (sAxisPhi \| sAxisMin));
350	if (areacode & sBoundary) areacode \|= sCorner; // xx is on the corner.
351	else areacode \|= sBoundary;
352
353	if (AmIOnLeftSide(xx, dphimin) > 0) isoutside = true;
354
355	} else if (AmIOnLeftSide(xx, dphimax) <= 0) { // xx is on dphimax
356
357	areacode \|= (sAxis1 & (sAxisPhi \| sAxisMax));
358	if (areacode & sBoundary) areacode \|= sCorner; // xx is on the corner.
359	else areacode \|= sBoundary;
360
361	if (AmIOnLeftSide(xx, dphimax) < 0) isoutside = true;
362
363	}
364
365	// if isoutside = true, clear inside bit.
366	// if not on boundary, add axis information.
367
368	if (isoutside) {
369	G4int tmpareacode = areacode & (~sInside);
370	areacode = tmpareacode;
371	} else if ((areacode & sBoundary) != sBoundary) {
372	areacode \|= (sAxis0 & sAxisRho) \| (sAxis1 & sAxisPhi);
373	}
374
375	} else {
376
377	// out of boundary of rho-axis
378
379	if (xx.getRho() < fAxisMin[rhoaxis]) {
380	areacode \|= (sAxis0 & (sAxisRho \| sAxisMin)) \| sBoundary;
381	} else if (xx.getRho() > fAxisMax[rhoaxis]) {
382	areacode \|= (sAxis0 & (sAxisRho \| sAxisMax)) \| sBoundary;
383	}
384
385	// out of boundary of phi-axis
386
387	if (AmIOnLeftSide(xx, dphimin, false) >= 0) { // xx is leftside or
388	areacode \|= (sAxis1 & (sAxisPhi \| sAxisMin)) ; // boundary of dphimin
389	if (areacode & sBoundary) areacode \|= sCorner; // xx is on the corner.
390	else areacode \|= sBoundary;
391
392	} else if (AmIOnLeftSide(xx, dphimax, false) <= 0) { // xx is rightside or
393	areacode \|= (sAxis1 & (sAxisPhi \| sAxisMax)) ; // boundary of dphimax
394	if (areacode & sBoundary) areacode \|= sCorner; // xx is on the corner.
395	else areacode \|= sBoundary;
396
397	}
398
399	if ((areacode & sBoundary) != sBoundary) {
400	areacode \|= (sAxis0 & sAxisRho) \| (sAxis1 & sAxisPhi);
401	}
402
403	}
404	return areacode;
405	} else {
406
407	G4cerr << "ERROR - G4TwistTubsFlatSide::GetAreaCode()" << G4endl
408	<< " fAxis[0] = " << fAxis[0] << G4endl
409	<< " fAxis[1] = " << fAxis[1] << G4endl;
410	G4Exception("G4TwistTubsFlatSide::GetAreaCode()", "NotImplemented",
411	FatalException, "Feature NOT implemented !");
412	}
413	return areacode;
414	}
415
416
417	//=====================================================================
418	//* SetCorners --------------------------------------------------------
419
420	void G4TwistTubsFlatSide::SetCorners()
421	{
422	// Set Corner points in local coodinate.
423
424	if (fAxis[0] == kRho && fAxis[1] == kPhi) {
425
426	G4int rhoaxis = 0; // kRho
427	G4int phiaxis = 1; // kPhi
428
429	G4double x, y, z;
430	// corner of Axis0min and Axis1min
431	x = fAxisMin[rhoaxis]*std::cos(fAxisMin[phiaxis]);
432	y = fAxisMin[rhoaxis]*std::sin(fAxisMin[phiaxis]);
433	z = 0;
434	SetCorner(sC0Min1Min, x, y, z);
435	// corner of Axis0max and Axis1min
436	x = fAxisMax[rhoaxis]*std::cos(fAxisMin[phiaxis]);
437	y = fAxisMax[rhoaxis]*std::sin(fAxisMin[phiaxis]);
438	z = 0;
439	SetCorner(sC0Max1Min, x, y, z);
440	// corner of Axis0max and Axis1max
441	x = fAxisMax[rhoaxis]*std::cos(fAxisMax[phiaxis]);
442	y = fAxisMax[rhoaxis]*std::sin(fAxisMax[phiaxis]);
443	z = 0;
444	SetCorner(sC0Max1Max, x, y, z);
445	// corner of Axis0min and Axis1max
446	x = fAxisMin[rhoaxis]*std::cos(fAxisMax[phiaxis]);
447	y = fAxisMin[rhoaxis]*std::sin(fAxisMax[phiaxis]);
448	z = 0;
449	SetCorner(sC0Min1Max, x, y, z);
450
451	} else {
452	G4cerr << "ERROR - G4TwistTubsFlatSide::SetCorners()" << G4endl
453	<< " fAxis[0] = " << fAxis[0] << G4endl
454	<< " fAxis[1] = " << fAxis[1] << G4endl;
455	G4Exception("G4TwistTubsFlatSide::SetCorners()", "NotImplemented",
456	FatalException, "Feature NOT implemented !");
457	}
458	}
459
460	//=====================================================================
461	//* SetBoundaries() ---------------------------------------------------
462
463	void G4TwistTubsFlatSide::SetBoundaries()
464	{
465	// Set direction-unit vector of phi-boundary-lines in local coodinate.
466	// Don't call the function twice.
467
468	if (fAxis[0] == kRho && fAxis[1] == kPhi) {
469
470	G4ThreeVector direction;
471	// sAxis0 & sAxisMin
472	direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
473	direction = direction.unit();
474	SetBoundary(sAxis0 & (sAxisPhi \| sAxisMin), direction,
475	GetCorner(sC0Min1Min), sAxisPhi);
476
477	// sAxis0 & sAxisMax
478	direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
479	direction = direction.unit();
480	SetBoundary(sAxis0 & (sAxisPhi \| sAxisMax), direction,
481	GetCorner(sC0Max1Min), sAxisPhi);
482
483	// sAxis1 & sAxisMin
484	direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
485	direction = direction.unit();
486	SetBoundary(sAxis1 & (sAxisRho \| sAxisMin), direction,
487	GetCorner(sC0Min1Min), sAxisRho);
488
489	// sAxis1 & sAxisMax
490	direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
491	direction = direction.unit();
492	SetBoundary(sAxis1 & (sAxisRho \| sAxisMax), direction,
493	GetCorner(sC0Min1Max), sAxisPhi);
494	} else {
495	G4cerr << "ERROR - G4TwistTubsFlatSide::SetBoundaries()" << G4endl
496	<< " fAxis[0] = " << fAxis[0] << G4endl
497	<< " fAxis[1] = " << fAxis[1] << G4endl;
498	G4Exception("G4TwistTubsFlatSide::SetBoundaries()", "NotImplemented",
499	FatalException, "Feature NOT implemented !");
500	}
501	}
502
503	//=====================================================================
504	//* GetFacets() -------------------------------------------------------
505
506	void G4TwistTubsFlatSide::GetFacets( G4int m, G4int n, G4double xyz[][3],
507	G4int faces[][4], G4int iside )
508	{
509
510	G4ThreeVector p ;
511
512	G4double rmin = fAxisMin[0] ;
513	G4double rmax = fAxisMax[0] ;
514	G4double phimin, phimax ;
515
516	G4double r,phi ;
517
518	G4int i,j ;
519
520	G4int nnode,nface ;
521
522	for ( i = 0 ; i<n ; i++ ) {
523
524	r = rmin + i*(rmax-rmin)/(n-1) ;
525
526	phimin = GetBoundaryMin(r) ;
527	phimax = GetBoundaryMax(r) ;
528
529	for ( j = 0 ; j<m ; j++ )
530	{
531	phi = phimin + j*(phimax-phimin)/(m-1) ;
532
533	nnode = GetNode(i,j,m,n,iside) ;
534	p = SurfacePoint(phi,r,true) ; // surface point in global coord.system
535
536	xyz[nnode][0] = p.x() ;
537	xyz[nnode][1] = p.y() ;
538	xyz[nnode][2] = p.z() ;
539
540	if ( i<n-1 && j<m-1 ) { // conterclock wise filling
541
542	nface = GetFace(i,j,m,n,iside) ;
543
544	if (fHandedness < 0) { // lower side
545	faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,-1) * ( GetNode(i ,j ,m,n,iside)+1) ;
546	faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,-1) * ( GetNode(i ,j+1,m,n,iside)+1) ;
547	faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,-1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
548	faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,-1) * ( GetNode(i+1,j ,m,n,iside)+1) ;
549	} else { // upper side
550	faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,1) * ( GetNode(i ,j ,m,n,iside)+1) ;
551	faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,1) * ( GetNode(i+1,j ,m,n,iside)+1) ;
552	faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
553	faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,1) * ( GetNode(i ,j+1,m,n,iside)+1) ;
554
555	}
556
557
558
559	}
560	}
561	}
562	}

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

Download in other formats: