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

Last change on this file since 1347 was 1347, checked in by garnier, 13 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.