Feature #1872 ยป time_structure.patch
src/common_py/ConfigurationDefaults.py 2016-09-16 10:32:44 +0000 | ||
---|---|---|
90 | 90 |
# Dictionary of variable to be set when using G4BL to generate particles. If |
91 | 91 |
# "get_magnet_currents_pa_cdb" is set to True magnet currents & proton absorber |
92 | 92 |
# thickness will be retrieved from the CDB for the run_number specified |
93 |
g4bl = {"run_number":2873,"q_1":1.066,"q_2":-1.332,"q_3":0.927,"d_1":-1.302,"d_2":-0.396,\ |
|
94 |
"d_s":3.837,"particles_per_spill":0,"rotation_angle":0,"translation_z":1000.0,\ |
|
95 |
"protonabsorberin":1,"proton_absorber_thickness":93,"proton_number":1E9,"proton_weight":1,\ |
|
96 |
"particle_charge":'all',"file_path":'MAUS_ROOT_DIR/src/map/MapPyBeamlineSimulation/G4bl',\ |
|
97 |
"get_magnet_currents_pa_cdb":False,"random_seed":1, |
|
93 |
g4bl = { |
|
94 |
"run_number":2873, |
|
95 |
"q_1":1.066, |
|
96 |
"q_2":-1.332, |
|
97 |
"q_3":0.927, |
|
98 |
"d_1":-1.302, |
|
99 |
"d_2":-0.396, |
|
100 |
"d_s":3.837, |
|
101 |
"particles_per_spill":0, |
|
102 |
"rotation_angle":0, |
|
103 |
"translation_z":1000.0, |
|
104 |
"protonabsorberin":1, |
|
105 |
"proton_absorber_thickness":93, |
|
106 |
"proton_number":1E9, |
|
107 |
"proton_weight":1, |
|
108 |
"particle_charge":'all', |
|
109 |
"file_path":'MAUS_ROOT_DIR/src/map/MapPyBeamlineSimulation/G4bl', |
|
110 |
"get_magnet_currents_pa_cdb":False, |
|
111 |
"random_seed":1, |
|
98 | 112 |
"seed_algorithm":"random_seed_and_spill_number", |
113 |
"time_structure":"uniform", |
|
114 |
"spill_gate":3e6, # 3e6 ns = 3 ms |
|
99 | 115 |
} |
100 | 116 | |
101 | 117 |
# Used by MapPyRemoveTracks. |
src/map/MapPyBeamlineSimulation/MapPyBeamlineSimulation.py 2016-09-16 11:11:47 +0000 | ||
---|---|---|
89 | 89 |
self.seed_algorithm = "use_random_seed" |
90 | 90 |
self.master_seed = 0 |
91 | 91 |
self.random_seed = 0 |
92 |
self.time_structure = "none" |
|
93 |
self.spill_gate = 0. |
|
92 | 94 | |
93 | 95 |
def birth(self, json_configuration): #pylint: disable=R0912, R0915 |
94 | 96 |
"birth doc string" |
... | ... | |
238 | 240 |
print('Error: Number of protons must be greater than zero!') |
239 | 241 |
good_birth = False |
240 | 242 |
except Exception: #pylint: disable=W0703 |
241 |
print("Error: proton_numer is not found in the config file!") |
|
243 |
print("Error: proton_number is not found in the config file!")
|
|
242 | 244 |
good_birth = False |
243 | 245 | |
244 | 246 |
try: |
... | ... | |
251 | 253 |
self.seed_algorithm = "use_random_seed" |
252 | 254 | |
253 | 255 |
try: |
256 |
self.time_structure = config_doc["g4bl"]["time_structure"] |
|
257 |
if self.time_structure not in ["uniform", "none"]: |
|
258 |
raise ValueError("time structure "+self.time_structure+\ |
|
259 |
" should be 'uniform' or 'none'") |
|
260 |
if self.time_structure == "uniform": |
|
261 |
self.spill_gate = float(config_doc["g4bl"]["spill_gate"]) |
|
262 |
if self.spill_gate < 0.: |
|
263 |
raise ValueError("spill gate should be > 0.") |
|
264 |
except KeyError: |
|
265 |
print config_doc["g4bl"] |
|
266 |
print("Error: time_structure set up error in the config file!") |
|
267 |
good_birth = False |
|
268 | ||
269 |
try: |
|
254 | 270 |
self.random_seed = config_doc["g4bl"]["random_seed"] |
255 | 271 |
if self.random_seed == -1: |
256 | 272 |
self.random_seed = random.randint(1, 1000000) |
... | ... | |
357 | 373 |
for i in range(self.queue.qsize()): |
358 | 374 |
primary = (self.queue.get_nowait()) |
359 | 375 |
spill['mc_events'].append({"primary":primary}) |
376 |
self.apply_time_structure(spill) |
|
360 | 377 |
with open('G4BLoutput.json', 'w') as outfile: |
361 | 378 |
json.dump(spill, outfile) |
362 |
#self.queue = 0
|
|
379 |
#self.queue = 0 |
|
363 | 380 |
return json.dumps(spill) |
364 | 381 |
else: |
365 | 382 |
spill = {} |
... | ... | |
400 | 417 |
for i in range(self.particles_per_spill): |
401 | 418 |
primary = (self.queue.get_nowait()) |
402 | 419 |
spill['mc_events'].append({"primary":primary}) |
420 |
self.apply_time_structure(spill) |
|
403 | 421 |
with open('G4BLoutputsimulation'+str(i)+'.json', 'w') \ |
404 | 422 |
as outfile: |
405 | 423 |
json.dump(spill, outfile) |
... | ... | |
408 | 426 |
for i in range(self.particles_per_spill): |
409 | 427 |
primary = (self.queue.get_nowait()) |
410 | 428 |
spill['mc_events'].append({"primary":primary}) |
429 |
self.apply_time_structure(spill) |
|
411 | 430 |
with open('G4BLoutputsimulation'+str(i)+'.json', 'w') \ |
412 | 431 |
as outfile: |
413 | 432 |
json.dump(spill, outfile) |
... | ... | |
415 | 434 |
else: |
416 | 435 |
print("Warning: G4BL simulated zero output particles!") |
417 | 436 |
raise IndexError("G4BL simulated zero output particles") |
437 | ||
438 |
def apply_time_structure(self, spill): |
|
439 |
""" |
|
440 |
apply time structure corresponding to the MICE spill gate |
|
441 |
- if self.time_structure is uniform, add time uniformly between 0. and |
|
442 |
spill gate to the G4BL output time structure |
|
443 |
""" |
|
444 |
if self.time_structure == "uniform": |
|
445 |
for event in spill["mc_events"]: |
|
446 |
numpy.random.seed(event["primary"]["random_seed"]) |
|
447 |
event["primary"]["time"] += \ |
|
448 |
numpy.random.uniform(0., self.spill_gate) |
|
418 | 449 |
|
419 | 450 |
def death(self): #pylint: disable = R0201 |
420 | 451 |
""" |
src/map/MapPyBeamlineSimulation/test_mappybeamlinesimulation.py 2016-09-16 11:10:19 +0000 | ||
---|---|---|
27 | 27 |
adjusting number of protons from target. |
28 | 28 |
""" #pylint: disable = W0105 |
29 | 29 |
TEST_1 = { |
30 |
"g4bl":{"q_1":1.16, "q_2":-1.45, "q_3":1.006, "d_1":-1.41,\ |
|
30 |
"g4bl":{"q_1":1.16, "q_2":-1.45, "q_3":1.006, "d_1":-1.41, "d_2":-0.7,\
|
|
31 | 31 |
"proton_absorber_thickness":11,'d_s':200,\ |
32 | 32 |
"proton_number":500000000,'particles_per_spill':0,\ |
33 | 33 |
"particle_charge":"all", 'rotation_angle':0, 'translation_z':20000.0,\ |
34 |
"random_seed":0,"run_number":0,"get_magnet_currents_pa_cdb":False,\ |
|
35 |
"proton_weight":1 |
|
34 |
"random_seed":0,"run_number":0,"get_magnet_currents_pa_cdb":False, |
|
35 |
"proton_weight":1, |
|
36 |
"protonabsorberin":0, |
|
37 |
"time_structure":"none", |
|
38 |
"spill_gate":3e6, |
|
36 | 39 | |
37 | 40 |
} |
38 | 41 |
} |
... | ... | |
47 | 50 |
"proton_number":500000000,'particles_per_spill':200,\ |
48 | 51 |
"particle_charge":"all", 'rotation_angle':0, 'translation_z':20000.0,\ |
49 | 52 |
"random_seed":0,"run_number":0,"get_magnet_currents_pa_cdb":False,\ |
50 |
"proton_weight":1 |
|
53 |
"proton_weight":1, |
|
54 |
"time_structure":"none", |
|
55 |
"spill_gate":3e6, |
|
51 | 56 | |
52 | 57 |
} |
53 | 58 |
} |
... | ... | |
62 | 67 |
"proton_number":500000000,'particles_per_spill':500,\ |
63 | 68 |
"particle_charge":"all", 'rotation_angle':0, 'translation_z':20000.0,\ |
64 | 69 |
"random_seed":0,"run_number":0,"get_magnet_currents_pa_cdb":False,\ |
65 |
"proton_weight":1 |
|
70 |
"proton_weight":1, |
|
71 |
"time_structure":"none", |
|
72 |
"spill_gate":3e6, |
|
66 | 73 | |
67 | 74 |
} |
68 | 75 |
} |
... | ... | |
77 | 84 |
"proton_number":10000,'particles_per_spill':200,\ |
78 | 85 |
"particle_charge":"all", 'rotation_angle':0, 'translation_z':20000.0,\ |
79 | 86 |
"random_seed":0,"run_number":0,"get_magnet_currents_pa_cdb":False,\ |
80 |
"proton_weight":1 |
|
87 |
"proton_weight":1, |
|
88 |
"time_structure":"none", |
|
89 |
"spill_gate":3e6, |
|
81 | 90 | |
82 | 91 |
} |
83 | 92 |
} |
... | ... | |
92 | 101 |
"proton_number":-10000,'particles_per_spill':200,\ |
93 | 102 |
"particle_charge":"all", 'rotation_angle':0, 'translation_z':20000.0,\ |
94 | 103 |
"random_seed":0,"run_number":0,"get_magnet_currents_pa_cdb":False,\ |
95 |
"proton_weight":1 |
|
104 |
"proton_weight":1, |
|
105 |
"time_structure":"none", |
|
106 |
"spill_gate":3e6, |
|
96 | 107 | |
97 | 108 |
} |
98 | 109 |
} |
... | ... | |
108 | 119 |
"particle_charge":"1.21 jiggawatts", 'rotation_angle':0,\ |
109 | 120 |
'translation_z':20000.0,\ |
110 | 121 |
"random_seed":0,"run_number":0,"get_magnet_currents_pa_cdb":False,\ |
111 |
"proton_weight":1 |
|
122 |
"proton_weight":1,
|
|
112 | 123 | |
113 | 124 |
} |
114 | 125 |
} |
... | ... | |
123 | 134 |
"particle_charge":"negative", 'rotation_angle':0,\ |
124 | 135 |
'translation_z':20000.0,\ |
125 | 136 |
"random_seed":0,"run_number":0,"get_magnet_currents_pa_cdb":False,\ |
126 |
"proton_weight":1 |
|
137 |
"proton_weight":1, |
|
138 |
"time_structure":"none", |
|
139 |
"spill_gate":3e6, |
|
127 | 140 | |
128 | 141 |
} |
129 | 142 |
} |
130 | 143 | |
131 | 144 |
""" |
132 | 145 |
Test random seed - seed from spill |
133 | ||
134 | 146 |
""" #pylint: disable = W0105 |
135 | 147 |
TEST_8 = { |
136 | 148 |
"g4bl":{"q_1":1.16, "q_2":-1.45, "q_3":1.006, "d_1":-1.41, "d_2":1., |
... | ... | |
141 | 153 |
"random_seed":-2,"run_number":0,"get_magnet_currents_pa_cdb":False, |
142 | 154 |
"proton_weight":1, |
143 | 155 |
"seed_algorithm":"random_seed_and_spill_number", |
156 |
"time_structure":"none", |
|
157 |
"spill_gate":3e6, |
|
158 |
} |
|
159 |
} |
|
160 | ||
161 |
""" |
|
162 |
Test time structure - uniform |
|
163 |
""" #pylint: disable = W0105 |
|
164 |
TEST_9 = { |
|
165 |
"g4bl":{"q_1":1.16, "q_2":-1.45, "q_3":1.006, "d_1":-1.41, "d_2":1., |
|
166 |
"proton_absorber_thickness":11,'d_s':200, |
|
167 |
"proton_number":100,'particles_per_spill':1000, |
|
168 |
"particle_charge":"positive", 'rotation_angle':0, |
|
169 |
'translation_z':20000.0, |
|
170 |
"random_seed":-2,"run_number":0,"get_magnet_currents_pa_cdb":False, |
|
171 |
"proton_weight":1, |
|
172 |
"seed_algorithm":"random_seed_and_spill_number", |
|
173 |
"time_structure":"uniform", |
|
174 |
"spill_gate":3e6, |
|
144 | 175 |
} |
145 | 176 |
} |
146 | 177 | |
... | ... | |
154 | 185 |
""" Initialise the class """ |
155 | 186 |
self.particles = {} |
156 | 187 |
for i in range(self.number_of_entries): |
157 |
self.particles['entry'+str(i)] = {} |
|
188 |
self.particles['entry'+str(i)] = { |
|
189 |
"position":{"x":0.0, "y":-0.0, "z":-0.0}, |
|
190 |
"momentum":{"x":0.0, "y":0.0, "z":1.0}, |
|
191 |
"particle_id":-13, |
|
192 |
"energy":226.0, "time":0.0, |
|
193 |
"random_seed":i, |
|
194 |
"spin":{"x":0.0, "y":0.0, "z":1.0} |
|
195 |
} |
|
196 | ||
158 | 197 | |
159 | 198 |
number_of_entries = 10 |
160 | 199 | |
... | ... | |
274 | 313 |
self.assertNotEqual(self.g4bl.random_seed, random_seed) |
275 | 314 |
MapPyBeamlineSimulation.CallG4bl = CallG4bl |
276 | 315 | |
316 |
def test_mappybeamline_time_structure(self): |
|
317 |
""" Test mappybeamline time structure""" |
|
318 |
MapPyBeamlineSimulation.CallG4bl = MockCallG4bl |
|
319 |
self.g4bl = MapPyBeamlineSimulation.MapPyBeamlineSimulation() |
|
320 | ||
321 |
self.g4bl.birth(json.dumps(TEST_8)) |
|
322 |
spill = self.g4bl.process(json.dumps({"spill_number":0})) |
|
323 |
spill = json.loads(spill) |
|
324 |
for event in spill["mc_events"]: |
|
325 |
self.assertAlmostEqual(event["primary"]["time"], 0.) |
|
326 |
self.g4bl.death() |
|
327 | ||
328 |
self.g4bl.birth(json.dumps(TEST_9)) |
|
329 |
spill = self.g4bl.process(json.dumps({"spill_number":0})) |
|
330 |
spill = json.loads(spill) |
|
331 |
time_list = [event["primary"]["time"] for event in spill["mc_events"]] |
|
332 |
bin_1, bin_2 = 0, 0 |
|
333 |
for time in time_list: |
|
334 |
if time > 0.-1e-9 and time < 1.5e6: |
|
335 |
bin_1 += 1 |
|
336 |
elif time >= 1.5e6 and time < 3.0e6+1e-9: |
|
337 |
bin_2 += 1 |
|
338 |
self.assertEqual(bin_1+bin_2, 1000) |
|
339 |
# should do KS test or something; this will do |
|
340 |
self.assertGreater(bin_1, 100) |
|
341 |
self.assertGreater(bin_2, 100) |
|
342 |
self.g4bl.death() |
|
343 |
|
|
344 | ||
277 | 345 | |
278 | 346 |
if __name__ == "__main__": |
279 | 347 |
unittest.main() |