source: trunk/source/geometry/solids/test/testSurfaceArea.cc@ 1354

Last change on this file since 1354 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 9.7 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: testSurfaceArea.cc,v 1.2 2007/02/12 11:29:23 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// --------------------------------------------------------------------
31// by Hans Dierckx to test surface Area
32
33// Validation of the Monte Carlo method to calculate arbitrary surface areas.
34// For the different CSG solids available, the surface is calculated
35// analytically and by the Monte Carlo method (after conversion into a
36// boolean intersection solid). Both results are compared afterwards.
37
38#include <assert.h>
39#include <cmath>
40#include <ctime>
41
42#include "globals.hh"
43#include "geomdefs.hh"
44
45#include "G4ThreeVector.hh"
46#include "G4TwoVector.hh"
47
48#include "G4Box.hh"
49#include "G4Tubs.hh"
50#include "G4Trap.hh"
51#include "G4Trd.hh"
52#include "G4Para.hh"
53#include "G4Sphere.hh"
54#include "G4Orb.hh"
55#include "G4Torus.hh"
56#include "G4Cons.hh"
57#include "G4BooleanSolid.hh"
58#include "G4IntersectionSolid.hh"
59#include "G4ExtrudedSolid.hh"
60
61#include "G4RotationMatrix.hh"
62#include "G4AffineTransform.hh"
63#include "G4VoxelLimits.hh"
64
65#include "Randomize.hh"
66
67G4bool testG4Surf()
68{
69 G4double surf, MCsurf;
70 G4int tic, toc;
71
72 G4Box w("Bigger Box",5,5,5);
73 G4cout << "\n\nCalculating Box:" << G4endl;
74 G4Box box("Test Box",1,0.5,3);
75 // fDx fDy fDz
76 surf = box.GetSurfaceArea();
77 G4cout <<"Surface = " << surf << G4endl;
78 G4IntersectionSolid box2("Box 2", &box, &w);
79 tic=clock();
80 MCsurf= box2.GetSurfaceArea();
81 toc = clock()-tic;
82 G4cout << "Elapsed time= "<< static_cast<G4double>(toc)/CLOCKS_PER_SEC << G4endl;
83 G4cout <<"MC result = " << MCsurf << G4endl;
84 G4cout << "deviation = " << (surf-MCsurf)/surf*100 << " %"<< G4endl;
85 G4cout <<"*******************" << G4endl;
86
87
88
89
90 G4cout << " Calculating Tube:" << G4endl;
91 // G4Tubs cyl("Test Cyl",1,4,3.5,pi/5,3*pi/2);
92 G4double R= std::sqrt(100./16./pi);
93 G4Tubs cyl("Test Cyl",0.5*R,1.5*R,1.5*R,0,twopi);
94 //pRmin pRmax pDz pSPhi pDPhi
95 surf = cyl.GetSurfaceArea();
96 G4cout << cyl.GetCubicVolume() << G4endl;
97 G4cout <<"Surface = " << surf << G4endl;
98 G4IntersectionSolid cyl2("Box 2", &cyl, &w);
99 tic=clock();
100 MCsurf=cyl2.GetSurfaceArea();
101 toc = clock()-tic;
102 G4cout << "Elapsed time= "<< static_cast<G4double>(toc)/CLOCKS_PER_SEC << G4endl;
103 G4cout <<"MC result = " << MCsurf << G4endl;
104 G4cout << "deviation = " << (surf-MCsurf)/surf*100 << " %"<< G4endl;
105 G4cout <<"*******************" << G4endl;
106
107 G4cout << " Calculating Cone:" << G4endl;
108 G4Cons cone("Test Cone",0.5,1,2,2.5,4,0,4*pi/3);
109 //pRmin1 pRmax1 pRmin2, pRmax2, pDz, pSphi pDphi
110 surf = cone.GetSurfaceArea();
111 G4cout <<"Surface = " << surf << G4endl;
112 G4IntersectionSolid cone2("Box 2", &cone, &w);
113 tic=clock();
114 MCsurf=cone2.GetSurfaceArea();
115 toc = clock()-tic;
116 G4cout << "Elapsed time= "<< static_cast<G4double>(toc)/CLOCKS_PER_SEC << G4endl;
117 G4cout <<"MC result = " << MCsurf << G4endl;
118 G4cout << "deviation = " << (surf-MCsurf)/surf*100 << " %"<< G4endl;
119 G4cout <<"*******************" << G4endl;
120
121
122 G4cout << " Calculating General G4Trap:" << G4endl;
123 G4Trap trap("Test Cyl",3, 20/180*pi,5/180*pi, 2, 1.5, 2, pi/18, 0.8 , 0.5, 0.7, pi/18);
124 //Dz,Theta, Phi, Dy1, Dx1, Dx2, alpha1, Dy2, Dx3,Dx4, alpha2
125 surf = trap.GetSurfaceArea();
126 G4cout <<"Surface = " << surf << G4endl;
127 G4IntersectionSolid trap2("Trap 2", &trap, &w);
128 tic=clock();
129 MCsurf=trap2.GetSurfaceArea();
130 toc = clock()-tic;
131 G4cout << "Elapsed time= "<< static_cast<G4double>(toc)/CLOCKS_PER_SEC << G4endl;
132 G4cout <<"MC result = " << MCsurf << G4endl;
133 G4cout << "deviation = " << (surf-MCsurf)/surf*100 << " %"<< G4endl;
134 G4cout <<"*******************" << G4endl;
135
136
137
138 G4cout << " Calculating Parallelepiped:" << G4endl;
139 G4Para para("Test Para", 1.5, 2, 3, pi/18, pi/9, 5/180*pi);
140 // Dx,Dy,Dz, alph, theta, phi
141 surf = para.GetSurfaceArea();
142 G4cout <<"Surface = " << surf << G4endl;
143 G4IntersectionSolid para2("Para 2", &para, &w);
144 tic=clock();
145 MCsurf=para2.GetSurfaceArea();
146 toc = clock()-tic;
147 G4cout << "Elapsed time= "<< static_cast<G4double>(toc)/CLOCKS_PER_SEC << G4endl;
148 G4cout <<"MC result = " << MCsurf << G4endl;
149 G4cout << "deviation = " << (surf-MCsurf)/surf*100 << " %"<< G4endl;
150 G4cout <<"*******************" << G4endl;
151
152
153
154 G4cout << " Calculating Trapezoid:" << G4endl;
155 // G4Trd trd("Test Trd", 3,1, 4,1.5, 3);
156 G4double a = std::sqrt(2.5);
157 G4Trd trd("Test Trd", a,2*a, 2*a,a,std::sqrt(3.)*0.5*a);
158 //Dx1,Dx2,Dy1,Dy2,Dz
159 surf = trd.GetSurfaceArea();
160 G4cout << trd.GetCubicVolume() << G4endl;
161 G4cout <<"Surface = " << surf << G4endl;
162 G4IntersectionSolid trd2("Trd 2", &trd, &w);
163 tic=clock();
164 MCsurf=trd2.GetSurfaceArea();
165 toc = clock()-tic;
166 G4cout << "Elapsed time= "<< static_cast<G4double>(toc)/CLOCKS_PER_SEC << G4endl;
167 G4cout <<"MC result = " << MCsurf << G4endl;
168 G4cout << "deviation = " << (surf-MCsurf)/surf*100 << " %"<< G4endl;
169 G4cout <<"*******************" << G4endl;
170
171
172
173 G4cout << " Calculating Sphere:" << G4endl;
174 G4Sphere sp("Test Sphere",2,3, 3*pi/4,pi, pi/6, pi);
175 //Rmin, Rmax, Sphi, Dphi, STheta, Dtheta
176 surf = sp.GetSurfaceArea();
177 G4cout <<"Surface = " << surf << G4endl;
178 G4IntersectionSolid sp2("Trd 2", &sp, &w);
179 tic=clock();
180 MCsurf=sp2.GetSurfaceArea();
181 toc = clock()-tic;
182 G4cout << "Elapsed time= "<< static_cast<G4double>(toc)/CLOCKS_PER_SEC << G4endl;
183 G4cout <<"MC result = " << MCsurf << G4endl;
184 G4cout << "deviation = " << (surf-MCsurf)/surf*100 << " %"<< G4endl;
185 G4cout <<"*******************" << G4endl;
186
187
188
189 G4cout << " Calculating Torus:" << G4endl;
190 G4Torus tor("Test Torus", 0.2,2,2.5, 0,0.2);
191 //Rmin Rmax Rtor Sphi Dphi
192 surf = tor.GetSurfaceArea();
193 G4cout <<"Surface = " << surf << G4endl;
194 G4IntersectionSolid tor2("Tor 2", &tor, &w);
195 tic=clock();
196 MCsurf=tor2.GetSurfaceArea();
197 toc = clock()-tic;
198 G4cout << "Elapsed time= "<< static_cast<G4double>(toc)/CLOCKS_PER_SEC << G4endl;
199 G4cout <<"MC result = " << MCsurf << G4endl;
200 G4cout << "deviation = " << (surf-MCsurf)/surf*100 << " %"<< G4endl;
201 G4cout <<"*******************" << G4endl;
202
203
204
205 G4cout << " Calculating Orb:" << G4endl;
206 G4Orb orb("Test Orb",2);
207 surf = orb.GetSurfaceArea();
208 G4cout <<"Surface = " << surf << G4endl;
209 G4IntersectionSolid orb2("Box 2", &orb, &w);
210 tic=clock();
211 MCsurf=orb2.GetSurfaceArea();
212 toc = clock()-tic;
213 G4cout << "Elapsed time= "<< static_cast<G4double>(toc)/CLOCKS_PER_SEC << G4endl;
214 G4cout <<"MC result = " << MCsurf << G4endl;
215 G4cout << "deviation = " << (surf-MCsurf)/surf*100 << " %"<< G4endl;
216 G4cout <<"*******************" << G4endl;
217
218
219
220 G4cout << " Calculating Extruded-Solid (against box):" << G4endl;
221 // Box above defined as Xtru ...
222 std::vector<G4TwoVector> polygon;
223 polygon.push_back(G4TwoVector(-1.,-0.5));
224 polygon.push_back(G4TwoVector(-1., 0.5));
225 polygon.push_back(G4TwoVector( 1., 0.5));
226 polygon.push_back(G4TwoVector( 1.,-0.5));
227
228 G4ExtrudedSolid xtru("xtru-box", polygon, 3,
229 G4TwoVector(), 1.0, G4TwoVector(), 1.0);
230
231 surf = xtru.GetSurfaceArea();
232 G4cout <<"Surface = " << surf << G4endl;
233 tic=clock();
234 MCsurf=box.GetSurfaceArea();
235 toc = clock()-tic;
236 G4cout << "Elapsed time= "<< static_cast<G4double>(toc)/CLOCKS_PER_SEC << G4endl;
237 G4cout <<"MC result = " << MCsurf << G4endl;
238 G4cout << "deviation = " << (surf-MCsurf)/surf*100 << " %"<< G4endl;
239 G4cout <<"*******************" << G4endl;
240
241 /////////////////////////////////////////////////////
242 return true;
243}
244
245int main()
246{
247
248// initialise random generator:
249CLHEP::RanluxEngine defaultEngine(1234567,4);
250CLHEP::HepRandom::setTheEngine(&defaultEngine);
251G4int seed = time(NULL);
252CLHEP::HepRandom::setTheSeed(seed);
253#ifdef NDEBUG
254 G4Exception("FAIL: *** Assertions must be compiled in! ***");
255#endif
256 assert(testG4Surf());
257 return 0;
258}
Note: See TracBrowser for help on using the repository browser.