1 package org.lcsim.detector.solids;
2
3 import hep.physics.vec.BasicHep3Vector;
4 import java.util.Random;
5 import junit.framework.TestCase;
6 import org.lcsim.detector.DetectorElement;
7 import org.lcsim.detector.IPhysicalVolume;
8 import org.lcsim.detector.IPhysicalVolumeNavigator;
9 import org.lcsim.detector.IRotation3D;
10 import org.lcsim.detector.ITransform3D;
11 import org.lcsim.detector.ITranslation3D;
12 import org.lcsim.detector.LogicalVolume;
13 import org.lcsim.detector.PhysicalVolume;
14 import org.lcsim.detector.PhysicalVolumeNavigatorStore;
15 import org.lcsim.detector.RotationGeant;
16 import org.lcsim.detector.Transform3D;
17 import org.lcsim.detector.Translation3D;
18 import org.lcsim.detector.converter.heprep.DetectorElementToHepRepConverter;
19 import org.lcsim.detector.material.IMaterial;
20 import org.lcsim.detector.material.MaterialElement;
21 import org.lcsim.util.test.TestUtil.TestOutputFile;
22
23 import static java.lang.Math.PI;
24 import static java.lang.Math.abs;
25 import static java.lang.Math.atan;
26 import static java.lang.Math.cos;
27
28 import static org.lcsim.units.clhep.SystemOfUnits.m;
29
30
31
32
33
34
35
36 public class RightIsoscelesTrapezoidTest extends TestCase
37 {
38
39 private boolean debug = false;
40 private boolean writeHeprep = false;
41 private static IMaterial dummymat = new MaterialElement("dummymat", 1, 1, 1.0);
42 IPhysicalVolumeNavigator nav;
43 IPhysicalVolume world;
44
45 public RightIsoscelesTrapezoidTest(String testName)
46 {
47 super(testName);
48 }
49
50 public void testRightIsoscelesTrapezoid() throws Exception
51 {
52 String name = "testSolid";
53
54 double b = 1.0;
55 double t = 1.0;
56 double height = 1.0;
57 double thickness = 1.0;
58
59 RightIsoscelesTrapezoid rit1 = new RightIsoscelesTrapezoid(name, b, t, height, thickness);
60
61
62
63 Box box = new Box(name, b, height, thickness);
64 if (debug)
65 System.out.println(rit1);
66
67 assertEquals(rit1.getCubicVolume(), box.getCubicVolume());
68
69
70
71
72
73
74
75
76 double r = 5.;
77
78 double sphereVolume = 4 * PI * r * r * r / 3.;
79 double polyVolume = rit1.getCubicVolume();
80
81 int nmax = 100000;
82 int nIn = 0;
83 Random ran = new Random();
84 int i = 0;
85 long startTime = System.nanoTime();
86 long endTime;
87
88 try
89 {
90 while (i < nmax)
91 {
92 double x = r * (2. * ran.nextDouble() - 1.);
93 double y = r * (2. * ran.nextDouble() - 1.);
94 double z = r * (2. * ran.nextDouble() - 1.);
95
96 if ((x * x + y * y + z * z) < r * r)
97 {
98 ++i;
99 BasicHep3Vector pos = new BasicHep3Vector(x, y, z);
100 Inside in = rit1.inside(pos);
101 Inside inBox = box.inside(pos);
102 assertTrue(in.equals(inBox));
103 if (Inside.INSIDE.compareTo(in) == 0)
104 ++nIn;
105 }
106 }
107 } finally
108 {
109 endTime = System.nanoTime();
110 }
111 long duration = endTime - startTime;
112 if (debug)
113 System.out.println("duration= " + duration);
114 double meas = (double) nIn / nmax;
115 double pred = polyVolume / sphereVolume;
116 if (debug)
117 System.out.println("nmax= " + nmax + " nIn= " + nIn + " ratio= " + meas);
118 if (debug)
119 System.out.println("volume ratio= " + pred);
120 double err = (meas - pred) / pred;
121 if (debug)
122 System.out.println("(meas-pred)/pred= " + err);
123
124 assertTrue(abs(err) < .1);
125
126
127
128 RightIsoscelesTrapezoid rit2 = new RightIsoscelesTrapezoid(name, b, 2. * t, height, thickness);
129 if (debug)
130 System.out.println(rit2);
131
132 nmax = 100000;
133 nIn = 0;
134 i = 0;
135 startTime = System.nanoTime();
136
137 try
138 {
139 while (i < nmax)
140 {
141 double x = r * (2. * ran.nextDouble() - 1.);
142 double y = r * (2. * ran.nextDouble() - 1.);
143 double z = r * (2. * ran.nextDouble() - 1.);
144
145 if ((x * x + y * y + z * z) < r * r)
146 {
147 ++i;
148 BasicHep3Vector pos = new BasicHep3Vector(x, y, z);
149 Inside in = rit1.inside(pos);
150 if (Inside.INSIDE.compareTo(in) == 0)
151 ++nIn;
152 }
153 }
154 } finally
155 {
156 endTime = System.nanoTime();
157 }
158 duration = endTime - startTime;
159 if (debug)
160 System.out.println("duration= " + duration);
161 meas = (double) nIn / nmax;
162 pred = polyVolume / sphereVolume;
163 if (debug)
164 System.out.println("nmax= " + nmax + " nIn= " + nIn + " ratio= " + meas);
165 if (debug)
166 System.out.println("volume ratio= " + pred);
167 err = (meas - pred) / pred;
168 if (debug)
169 System.out.println("(meas-pred)/pred= " + err);
170
171 assertTrue(abs(err) < .1);
172
173
174
175 double skinThickness = .1*b;
176 double phi = atan((t-b)/2*height);
177 double l = skinThickness/cos(phi);
178 double newTop = 2.*t-l;
179 double newBottom = b-l;
180 RightIsoscelesTrapezoid rit3 = new RightIsoscelesTrapezoid(name, newBottom, newTop, height, thickness);
181
182
183
184 if (writeHeprep)
185 {
186 world = createWorld();
187 nav = PhysicalVolumeNavigatorStore.getInstance().createDefault(world);
188
189 IRotation3D rotation = new RotationGeant(0, 0, 0);
190 ITranslation3D translation = new Translation3D(0, 0, 0);
191 ITransform3D transform = new Transform3D(translation, rotation);
192 LogicalVolume lvTest = new LogicalVolume("lvtest", rit2, dummymat);
193 new PhysicalVolume(
194 transform,
195 "pvtest",
196 lvTest,
197 world.getLogicalVolume(),
198 0);
199 new DetectorElement("detest", null, "/pvtest");
200
201
202 LogicalVolume lvTest2 = new LogicalVolume("lvtest2", rit3, dummymat);
203 new PhysicalVolume(
204 transform,
205 "pvtest2",
206 lvTest2,
207 world.getLogicalVolume(),
208 0);
209 new DetectorElement("detest2", null, "/pvtest2");
210
211
212
213
214 DetectorElementToHepRepConverter.writeHepRep(new TestOutputFile("RightIsoscelesTrapezoidTest.heprep").getAbsolutePath());
215
216 }
217
218 }
219
220 final IPhysicalVolume createWorld()
221 {
222 Box boxWorld = new Box(
223 "world_box",
224 10.0 * m,
225 10.0 * m,
226 10.0 * m);
227
228 LogicalVolume lvWorld =
229 new LogicalVolume(
230 "world",
231 boxWorld,
232 dummymat);
233
234 IPhysicalVolume pvTop =
235 new PhysicalVolume(
236 null,
237 "world",
238 lvWorld,
239 null,
240 0);
241
242 return pvTop;
243 }
244 }