| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 #include "ChromosomeFactory.hh"
28
29 #include "CylindricalChromosome.hh"
30 #include "EllipticalChromosome.hh"
31 #include "RodChromosome.hh"
32 #include "SphericalChromosome.hh"
33 #include "BoxChromosome.hh"
34 #include "VirtualChromosome.hh"
35
36 #include "G4PhysicalConstants.hh"
37 #include "G4ThreeVector.hh"
38 #include "G4UnitsTable.hh"
39 #include "Randomize.hh"
40
41 //....oooOO0OOooo........oooOO0OOooo........oo
42
43 VirtualChromosome* ChromosomeFactory::MakeChro
44
45 {
46 VirtualChromosome* chromosome = nullptr;
47 G4String chromosome_type = commands[0];
48 if (chromosome_type == CylindricalChromosome
49 // interpret the command for cylinder
50 // expect cyl rad height x y z unit rx ry
51 // rotations are in degrees
52 if (commands.size() == 7) {
53 G4double unit = G4UnitDefinition::GetVal
54 G4double radius = std::stod(commands[1])
55 G4double hgt = std::stod(commands[2]) *
56 G4ThreeVector center(std::stod(commands[
57 std::stod(commands[
58 chromosome = new CylindricalChromosome(n
59 }
60 else if (commands.size() == 10) // euler
61 {
62 G4double unit = G4UnitDefinition::GetVal
63 G4double radius = std::stod(commands[1])
64 G4double hgt = std::stod(commands[2]) *
65 G4ThreeVector center(std::stod(commands[
66 std::stod(commands[
67
68 // Rotations are first X, then Y, Then Z
69 G4RotationMatrix rot;
70 rot.rotateX(std::stod(commands[7]) * pi
71 rot.rotateY(std::stod(commands[8]) * pi
72 rot.rotateZ(std::stod(commands[9]) * pi
73
74 chromosome = new CylindricalChromosome(n
75 }
76 else {
77 G4cout << "The arguments for a cylinder
78 G4cout << "1) name cyl rad height x y
79 G4cout << "2) name cyl rad height x y
80 G4cout << "Note that rotations are in de
81 InvalidReading(chromosome_type);
82 }
83 }
84 else if (chromosome_type == RodChromosome::f
85 // interpret the command for rod
86 // expect cyl rad height x y z unit rx ry
87 // rotations are in degrees
88 if (commands.size() == 7) {
89 G4double unit = G4UnitDefinition::GetVal
90 G4double radius = std::stod(commands[1])
91 G4double hgt = std::stod(commands[2]) *
92 G4ThreeVector center(std::stod(commands[
93 std::stod(commands[
94 chromosome = new RodChromosome(name, cen
95 }
96 else if (commands.size() == 10) // euler
97 {
98 G4double unit = G4UnitDefinition::GetVal
99 G4double radius = std::stod(commands[1])
100 G4double hgt = std::stod(commands[2]) *
101 G4ThreeVector center(std::stod(commands[
102 std::stod(commands[
103
104 // Rotations are first X, then Y, Then Z
105 G4RotationMatrix rot;
106 rot.rotateX(std::stod(commands[7]) * pi
107 rot.rotateY(std::stod(commands[8]) * pi
108 rot.rotateZ(std::stod(commands[9]) * pi
109
110 chromosome = new RodChromosome(name, cen
111 }
112 else {
113 G4cout << "The arguments for a cylinder
114 G4cout << "1) name rod rad height x y
115 G4cout << "2) name rod rad height x y
116 G4cout << "Note that rotations are in de
117 InvalidReading(chromosome_type);
118 }
119 }
120 else if (chromosome_type == SphericalChromos
121 // interpret the command for sphere
122 // expect sphere rad x y z unit rx ry rz
123 // rotations are in degrees
124 if (commands.size() == 6) {
125 G4double unit = G4UnitDefinition::GetVal
126 G4double radius = std::stod(commands[1])
127 G4ThreeVector center(std::stod(commands[
128 std::stod(commands[
129 chromosome = new SphericalChromosome(nam
130 }
131 else if (commands.size() == 9) // euler a
132 {
133 G4double unit = G4UnitDefinition::GetVal
134 G4double radius = std::stod(commands[1])
135 G4ThreeVector center(std::stod(commands[
136 std::stod(commands[
137
138 // Rotations are first X, then Y, Then Z
139 G4RotationMatrix rot;
140 rot.rotateX(std::stod(commands[6]) * pi
141 rot.rotateY(std::stod(commands[7]) * pi
142 rot.rotateZ(std::stod(commands[8]) * pi
143
144 chromosome = new SphericalChromosome(nam
145 }
146 else {
147 G4cout << "The arguments for a sphere ar
148 G4cout << "1) name sphere rad x y z u
149 G4cout << "2) name sphere rad x y z u
150 G4cout << "Note that rotations are in de
151 InvalidReading(chromosome_type);
152 }
153 }
154 else if (chromosome_type == EllipticalChromo
155 // interpret the command for Ellipse
156 // expect ellipse sx sy sz x y z unit rx r
157 // rotations are in degrees
158 // sx sy sz are semi-major axes
159 if (commands.size() == 8) {
160 G4double unit = G4UnitDefinition::GetVal
161 G4double sx = std::stod(commands[1]) * u
162 G4double sy = std::stod(commands[2]) * u
163 G4double sz = std::stod(commands[3]) * u
164 G4ThreeVector center(std::stod(commands[
165 std::stod(commands[
166 chromosome = new EllipticalChromosome(na
167 }
168 else if (commands.size() == 11) // euler
169 {
170 G4double unit = G4UnitDefinition::GetVal
171 G4double sx = std::stod(commands[1]) * u
172 G4double sy = std::stod(commands[2]) * u
173 G4double sz = std::stod(commands[3]) * u
174 G4ThreeVector center(std::stod(commands[
175 std::stod(commands[
176
177 // Rotations are first X, then Y, Then Z
178 G4RotationMatrix rot;
179 rot.rotateX(std::stod(commands[8]) * pi
180 rot.rotateY(std::stod(commands[9]) * pi
181 rot.rotateZ(std::stod(commands[10]) * pi
182
183 chromosome = new EllipticalChromosome(na
184 }
185 else {
186 G4cout << "The arguments for a ellipse a
187 G4cout << "1) name ellipse sx sy sz x
188 G4cout << "2) name ellipse sx sy sz x
189 G4cout << "Note that rotations are in de
190 G4cout << "Note that dimensions (sx, sy,
191 InvalidReading(chromosome_type);
192 }
193 }
194 // added by sara
195
196 else if(chromosome_type == BoxChromosome::fS
197 {
198 // interpret the command for Box
199 // expect ellipse xdim ydim zdim x y z uni
200 // rotations are in degrees
201 // xdim ydim zdim are edges of box
202 if(commands.size() == 8)
203 {
204 G4double unit = G4UnitDefinition::GetVal
205 G4double xdim = std::stod(commands[1])
206 G4double ydim = std::stod(commands[2])
207 G4double zdim = std::stod(commands[3])
208 G4ThreeVector center(std::stod(commands[
209 std::stod(commands[
210 std::stod(commands[
211 chromosome = new BoxChromosome(name, cen
212 }
213 else if(commands.size() == 11) // euler a
214 {
215 G4double unit = G4UnitDefinition::GetVal
216 G4double xdim = std::stod(commands[1])
217 G4double ydim = std::stod(commands[2])
218 G4double zdim = std::stod(commands[3])
219 G4ThreeVector center(std::stod(commands[
220 std::stod(commands[
221 std::stod(commands[
222
223 // Rotations are first X, then Y, Then Z
224 G4RotationMatrix rot;
225 rot.rotateX(std::stod(commands[8]) * pi
226 rot.rotateY(std::stod(commands[9]) * pi
227 rot.rotateZ(std::stod(commands[10]) * pi
228
229 chromosome = new BoxChromosome(name, cen
230 }
231 else
232 {
233 G4cout << "The arguments for a box are:"
234 G4cout << "1) name box xdim ydim zdim
235 G4cout << "2) name box xdim ydim zdim
236 G4cout << "Note that rotations are in de
237 G4cout << "Note that dimensions (xdim, y
238 << G4endl;
239 InvalidReading(chromosome_type);
240 }
241 }
242 else
243 {
244 chromosome = nullptr;
245 InvalidReading(chromosome_type);
246 }
247 return chromosome;
248 }
249
250 //....oooOO0OOooo........oooOO0OOooo........oo
251 void ChromosomeFactory::InvalidReading(const G
252 {
253 G4ExceptionDescription errmsg;
254 errmsg << "Chromosome type: " << chromosome_
255 G4Exception("ChromosomeFactory::MakeChromoso
256 }
257 //....oooOO0OOooo........oooOO0OOooo........oo
258
259 void ChromosomeFactory::Test()
260 {
261 G4cout << "---------------------------------
262 << "--------------------" << G4endl;
263 G4cout << "Chromosome Test" << G4endl;
264 // populate vector for tests
265 std::vector<std::vector<G4String>*> tests;
266 G4double r;
267 for (G4int ii = 0; ii != 50; ii++) {
268 // No rotation constructor test
269 auto vec = new std::vector<G4String>();
270 r = G4UniformRand();
271 if (r < 0.25) {
272 vec->push_back("rod");
273 vec->push_back(std::to_string(10 * G4Uni
274 vec->push_back(std::to_string(5 * G4Unif
275 }
276 else if (r < 0.5) {
277 vec->push_back("sphere");
278 vec->push_back(std::to_string(10 * G4Uni
279 }
280 else if (r < 0.75) {
281 vec->push_back("cyl");
282 vec->push_back(std::to_string(10 * G4Uni
283 vec->push_back(std::to_string(5 * G4Unif
284 }
285 else {
286 vec->push_back("ellipse");
287 vec->push_back(std::to_string(10 * G4Uni
288 vec->push_back(std::to_string(5 * G4Unif
289 vec->push_back(std::to_string(20 * G4Uni
290 }
291
292 vec->push_back(std::to_string(10 * (G4Unif
293 vec->push_back(std::to_string(10 * (G4Unif
294 vec->push_back(std::to_string(10 * (G4Unif
295 vec->push_back("um");
296 tests.push_back(vec);
297 }
298 for (G4int ii = 0; ii != 50; ++ii) {
299 // Test with rotation
300 auto vec = new std::vector<G4String>();
301 r = G4UniformRand();
302 if (r < 0.25) {
303 vec->push_back("rod");
304 vec->push_back(std::to_string(10 * G4Uni
305 vec->push_back(std::to_string(5 * G4Unif
306 }
307 else if (r < 0.5) {
308 vec->push_back("sphere");
309 vec->push_back(std::to_string(10 * G4Uni
310 }
311 else if (r < 0.75) {
312 vec->push_back("cyl");
313 vec->push_back(std::to_string(10 * G4Uni
314 vec->push_back(std::to_string(5 * G4Unif
315 }
316 else {
317 vec->push_back("ellipse");
318 vec->push_back(std::to_string(10 * G4Uni
319 vec->push_back(std::to_string(5 * G4Unif
320 vec->push_back(std::to_string(20 * G4Uni
321 }
322
323 vec->push_back(std::to_string(10 * (G4Unif
324 vec->push_back(std::to_string(10 * (G4Unif
325 vec->push_back(std::to_string(10 * (G4Unif
326
327 vec->push_back("um");
328
329 vec->push_back(std::to_string(720 * (G4Uni
330 vec->push_back(std::to_string(720 * (G4Uni
331 vec->push_back(std::to_string(720 * (G4Uni
332
333 tests.push_back(vec);
334 }
335
336 VirtualChromosome* chromo;
337 for (auto& test : tests) {
338 chromo = this->MakeChromosome("test", (*te
339 G4int passes = 0;
340 G4int n = 1000;
341 for (G4int jj = 0; jj != n; jj++) {
342 G4ThreeVector point = chromo->RandomPoin
343 if (chromo->PointInChromosome(point)) {
344 passes++;
345 }
346 else {
347 G4cout << point << G4endl;
348 }
349 }
350 if (passes == n) {
351 G4cout << "Chromosome Test Passed for "
352 << " test points" << G4endl;
353 }
354 else {
355 G4cout << "Chromosome Test Failed for "
356 << " test points (only " << passe
357 chromo->Print();
358 }
359 delete chromo;
360 }
361
362 for (auto& test : tests) {
363 delete test;
364 }
365 G4cout << "---------------------------------
366 << "--------------------" << G4endl;
367 }
368
369 //....oooOO0OOooo........oooOO0OOooo........oo