| Geant4 Cross Reference |
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 #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........oooOO0OOooo........oooOO0OOooo......
42
43 VirtualChromosome* ChromosomeFactory::MakeChromosome(const G4String& name,
44 const std::vector<G4String>& commands)
45 {
46 VirtualChromosome* chromosome = nullptr;
47 G4String chromosome_type = commands[0];
48 if (chromosome_type == CylindricalChromosome::fShape) {
49 // interpret the command for cylinder
50 // expect cyl rad height x y z unit rx ry rz
51 // rotations are in degrees
52 if (commands.size() == 7) {
53 G4double unit = G4UnitDefinition::GetValueOf(commands[6]);
54 G4double radius = std::stod(commands[1]) * unit;
55 G4double hgt = std::stod(commands[2]) * unit;
56 G4ThreeVector center(std::stod(commands[3]) * unit, std::stod(commands[4]) * unit,
57 std::stod(commands[5]) * unit);
58 chromosome = new CylindricalChromosome(name, center, radius, hgt);
59 }
60 else if (commands.size() == 10) // euler angles given
61 {
62 G4double unit = G4UnitDefinition::GetValueOf(commands[6]);
63 G4double radius = std::stod(commands[1]) * unit;
64 G4double hgt = std::stod(commands[2]) * unit;
65 G4ThreeVector center(std::stod(commands[3]) * unit, std::stod(commands[4]) * unit,
66 std::stod(commands[5]) * unit);
67
68 // Rotations are first X, then Y, Then Z
69 G4RotationMatrix rot;
70 rot.rotateX(std::stod(commands[7]) * pi / 180);
71 rot.rotateY(std::stod(commands[8]) * pi / 180);
72 rot.rotateZ(std::stod(commands[9]) * pi / 180);
73
74 chromosome = new CylindricalChromosome(name, center, radius, hgt, rot);
75 }
76 else {
77 G4cout << "The arguments for a cylinder are:" << G4endl;
78 G4cout << "1) name cyl rad height x y z unit" << G4endl;
79 G4cout << "2) name cyl rad height x y z unit rx ry rz" << G4endl;
80 G4cout << "Note that rotations are in degrees" << G4endl;
81 InvalidReading(chromosome_type);
82 }
83 }
84 else if (chromosome_type == RodChromosome::fShape) {
85 // interpret the command for rod
86 // expect cyl rad height x y z unit rx ry rz
87 // rotations are in degrees
88 if (commands.size() == 7) {
89 G4double unit = G4UnitDefinition::GetValueOf(commands[6]);
90 G4double radius = std::stod(commands[1]) * unit;
91 G4double hgt = std::stod(commands[2]) * unit;
92 G4ThreeVector center(std::stod(commands[3]) * unit, std::stod(commands[4]) * unit,
93 std::stod(commands[5]) * unit);
94 chromosome = new RodChromosome(name, center, radius, hgt);
95 }
96 else if (commands.size() == 10) // euler angles given
97 {
98 G4double unit = G4UnitDefinition::GetValueOf(commands[6]);
99 G4double radius = std::stod(commands[1]) * unit;
100 G4double hgt = std::stod(commands[2]) * unit;
101 G4ThreeVector center(std::stod(commands[3]) * unit, std::stod(commands[4]) * unit,
102 std::stod(commands[5]) * unit);
103
104 // Rotations are first X, then Y, Then Z
105 G4RotationMatrix rot;
106 rot.rotateX(std::stod(commands[7]) * pi / 180);
107 rot.rotateY(std::stod(commands[8]) * pi / 180);
108 rot.rotateZ(std::stod(commands[9]) * pi / 180);
109
110 chromosome = new RodChromosome(name, center, radius, hgt, rot);
111 }
112 else {
113 G4cout << "The arguments for a cylinder are:" << G4endl;
114 G4cout << "1) name rod rad height x y z unit" << G4endl;
115 G4cout << "2) name rod rad height x y z unit rx ry rz" << G4endl;
116 G4cout << "Note that rotations are in degrees" << G4endl;
117 InvalidReading(chromosome_type);
118 }
119 }
120 else if (chromosome_type == SphericalChromosome::fShape) {
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::GetValueOf(commands[5]);
126 G4double radius = std::stod(commands[1]) * unit;
127 G4ThreeVector center(std::stod(commands[2]) * unit, std::stod(commands[3]) * unit,
128 std::stod(commands[4]) * unit);
129 chromosome = new SphericalChromosome(name, center, radius);
130 }
131 else if (commands.size() == 9) // euler angles given
132 {
133 G4double unit = G4UnitDefinition::GetValueOf(commands[5]);
134 G4double radius = std::stod(commands[1]) * unit;
135 G4ThreeVector center(std::stod(commands[2]) * unit, std::stod(commands[3]) * unit,
136 std::stod(commands[4]) * unit);
137
138 // Rotations are first X, then Y, Then Z
139 G4RotationMatrix rot;
140 rot.rotateX(std::stod(commands[6]) * pi / 180);
141 rot.rotateY(std::stod(commands[7]) * pi / 180);
142 rot.rotateZ(std::stod(commands[8]) * pi / 180);
143
144 chromosome = new SphericalChromosome(name, center, radius, rot);
145 }
146 else {
147 G4cout << "The arguments for a sphere are:" << G4endl;
148 G4cout << "1) name sphere rad x y z unit" << G4endl;
149 G4cout << "2) name sphere rad x y z unit rx ry rz" << G4endl;
150 G4cout << "Note that rotations are in degrees" << G4endl;
151 InvalidReading(chromosome_type);
152 }
153 }
154 else if (chromosome_type == EllipticalChromosome::fShape) {
155 // interpret the command for Ellipse
156 // expect ellipse sx sy sz x y z unit rx ry rz
157 // rotations are in degrees
158 // sx sy sz are semi-major axes
159 if (commands.size() == 8) {
160 G4double unit = G4UnitDefinition::GetValueOf(commands[7]);
161 G4double sx = std::stod(commands[1]) * unit;
162 G4double sy = std::stod(commands[2]) * unit;
163 G4double sz = std::stod(commands[3]) * unit;
164 G4ThreeVector center(std::stod(commands[4]) * unit, std::stod(commands[5]) * unit,
165 std::stod(commands[6]) * unit);
166 chromosome = new EllipticalChromosome(name, center, sx, sy, sz);
167 }
168 else if (commands.size() == 11) // euler angles given
169 {
170 G4double unit = G4UnitDefinition::GetValueOf(commands[7]);
171 G4double sx = std::stod(commands[1]) * unit;
172 G4double sy = std::stod(commands[2]) * unit;
173 G4double sz = std::stod(commands[3]) * unit;
174 G4ThreeVector center(std::stod(commands[4]) * unit, std::stod(commands[5]) * unit,
175 std::stod(commands[6]) * unit);
176
177 // Rotations are first X, then Y, Then Z
178 G4RotationMatrix rot;
179 rot.rotateX(std::stod(commands[8]) * pi / 180);
180 rot.rotateY(std::stod(commands[9]) * pi / 180);
181 rot.rotateZ(std::stod(commands[10]) * pi / 180);
182
183 chromosome = new EllipticalChromosome(name, center, sx, sy, sz, rot);
184 }
185 else {
186 G4cout << "The arguments for a ellipse are:" << G4endl;
187 G4cout << "1) name ellipse sx sy sz x y z unit" << G4endl;
188 G4cout << "2) name ellipse sx sy sz x y z unit rx ry rz" << G4endl;
189 G4cout << "Note that rotations are in degrees" << G4endl;
190 G4cout << "Note that dimensions (sx, sy, sz) are semi-major axes" << G4endl;
191 InvalidReading(chromosome_type);
192 }
193 }
194 // added by sara
195
196 else if(chromosome_type == BoxChromosome::fShape)
197 {
198 // interpret the command for Box
199 // expect ellipse xdim ydim zdim x y z unit rx ry rz
200 // rotations are in degrees
201 // xdim ydim zdim are edges of box
202 if(commands.size() == 8)
203 {
204 G4double unit = G4UnitDefinition::GetValueOf(commands[7]);
205 G4double xdim = std::stod(commands[1]) * unit;
206 G4double ydim = std::stod(commands[2]) * unit;
207 G4double zdim = std::stod(commands[3]) * unit;
208 G4ThreeVector center(std::stod(commands[4]) * unit,
209 std::stod(commands[5]) * unit,
210 std::stod(commands[6]) * unit);
211 chromosome = new BoxChromosome(name, center, xdim, ydim, zdim);
212 }
213 else if(commands.size() == 11) // euler angles given
214 {
215 G4double unit = G4UnitDefinition::GetValueOf(commands[7]);
216 G4double xdim = std::stod(commands[1]) * unit;
217 G4double ydim = std::stod(commands[2]) * unit;
218 G4double zdim = std::stod(commands[3]) * unit;
219 G4ThreeVector center(std::stod(commands[4]) * unit,
220 std::stod(commands[5]) * unit,
221 std::stod(commands[6]) * unit);
222
223 // Rotations are first X, then Y, Then Z
224 G4RotationMatrix rot;
225 rot.rotateX(std::stod(commands[8]) * pi / 180);
226 rot.rotateY(std::stod(commands[9]) * pi / 180);
227 rot.rotateZ(std::stod(commands[10]) * pi / 180);
228
229 chromosome = new BoxChromosome(name, center, xdim, ydim, zdim, rot);
230 }
231 else
232 {
233 G4cout << "The arguments for a box are:" << G4endl;
234 G4cout << "1) name box xdim ydim zdim x y z unit" << G4endl;
235 G4cout << "2) name box xdim ydim zdim x y z unit rx ry rz" << G4endl;
236 G4cout << "Note that rotations are in degrees" << G4endl;
237 G4cout << "Note that dimensions (xdim, ydim, zdim) are box edges"
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........oooOO0OOooo........oooOO0OOooo......
251 void ChromosomeFactory::InvalidReading(const G4String& chromosome_type)
252 {
253 G4ExceptionDescription errmsg;
254 errmsg << "Chromosome type: " << chromosome_type << " is not valid" << G4endl;
255 G4Exception("ChromosomeFactory::MakeChromosome", "", FatalException, errmsg);
256 }
257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 * G4UniformRand() + 2));
274 vec->push_back(std::to_string(5 * G4UniformRand() + 10));
275 }
276 else if (r < 0.5) {
277 vec->push_back("sphere");
278 vec->push_back(std::to_string(10 * G4UniformRand() + 2));
279 }
280 else if (r < 0.75) {
281 vec->push_back("cyl");
282 vec->push_back(std::to_string(10 * G4UniformRand() + 2));
283 vec->push_back(std::to_string(5 * G4UniformRand() + 5));
284 }
285 else {
286 vec->push_back("ellipse");
287 vec->push_back(std::to_string(10 * G4UniformRand() + 2));
288 vec->push_back(std::to_string(5 * G4UniformRand() + 5));
289 vec->push_back(std::to_string(20 * G4UniformRand() - 6));
290 }
291
292 vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
293 vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
294 vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
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 * G4UniformRand() + 2));
305 vec->push_back(std::to_string(5 * G4UniformRand() + 10));
306 }
307 else if (r < 0.5) {
308 vec->push_back("sphere");
309 vec->push_back(std::to_string(10 * G4UniformRand() + 2));
310 }
311 else if (r < 0.75) {
312 vec->push_back("cyl");
313 vec->push_back(std::to_string(10 * G4UniformRand() + 2));
314 vec->push_back(std::to_string(5 * G4UniformRand() + 5));
315 }
316 else {
317 vec->push_back("ellipse");
318 vec->push_back(std::to_string(10 * G4UniformRand() + 2));
319 vec->push_back(std::to_string(5 * G4UniformRand() + 5));
320 vec->push_back(std::to_string(20 * G4UniformRand() - 6));
321 }
322
323 vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
324 vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
325 vec->push_back(std::to_string(10 * (G4UniformRand() - 0.5)));
326
327 vec->push_back("um");
328
329 vec->push_back(std::to_string(720 * (G4UniformRand() - 0.5)));
330 vec->push_back(std::to_string(720 * (G4UniformRand() - 0.5)));
331 vec->push_back(std::to_string(720 * (G4UniformRand() - 0.5)));
332
333 tests.push_back(vec);
334 }
335
336 VirtualChromosome* chromo;
337 for (auto& test : tests) {
338 chromo = this->MakeChromosome("test", (*test));
339 G4int passes = 0;
340 G4int n = 1000;
341 for (G4int jj = 0; jj != n; jj++) {
342 G4ThreeVector point = chromo->RandomPointInChromosome();
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 " << chromo->GetShape() << " based on " << n
352 << " test points" << G4endl;
353 }
354 else {
355 G4cout << "Chromosome Test Failed for " << chromo->GetShape() << " based on " << n
356 << " test points (only " << passes << " pass)" << G4endl;
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........oooOO0OOooo........oooOO0OOooo......