From 945179f7c1e763412b4d2d5a32f191c94342da6f Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Mon, 8 Aug 2011 11:51:48 -0400 Subject: Switch sim.py to spawn a separate process for GEANT4 and start it generating photons while the detector geometry is being built. --- sim.py | 88 ++++++++++++++++++++++++++++++++++-------------------------------- 1 file changed, 46 insertions(+), 42 deletions(-) diff --git a/sim.py b/sim.py index b8deef6..44fb1dd 100755 --- a/sim.py +++ b/sim.py @@ -2,6 +2,7 @@ import sys import optparse import time +import multiprocessing import detectors import optics @@ -27,8 +28,31 @@ def info(type, value, tb): sys.excepthook = info -if __name__ == '__main__': +class GeneratorThread(multiprocessing.Process): + def __init__(self, particle, energy, position, direction, nevents, material): + multiprocessing.Process.__init__(self) + + self.particle = particle + self.energy = energy + self.position = position + self.direction = direction + self.nevents = nevents + self.material = material + self.queue = multiprocessing.Queue() + + def run(self): + print >>sys.stderr, 'Starting generator thread...' + generator = g4gen.G4Generator(self.material) + + for i in xrange(self.nevents): + photons = generator.generate_photons(particle_name=self.particle, + total_energy=self.energy, + position=self.position, + direction=self.direction) + self.queue.put(photons) + +def main(): parser = optparse.OptionParser('%prog') parser.add_option('-b', type='int', dest='nbits', default=10) parser.add_option('-j', type='int', dest='device', default=None) @@ -59,11 +83,20 @@ if __name__ == '__main__': print >>sys.stderr, 'Creating detector...' detector.build(bits=options.nbits) - print >>sys.stderr, 'Initializing generator...' + print >>sys.stderr, 'Creating generator...' detector_material = optics.water - generator = g4gen.G4Generator(detector_material) + generator_thread = GeneratorThread(particle=options.particle, + energy=options.energy, + position=position, + direction=direction, + nevents=options.nevents, + material=detector_material) print >>sys.stderr, 'WARNING: ASSUMING DETECTOR IS WATER!!' + # Do this now so we can get ahead of the photon propagation + print >>sys.stderr, 'Starting GEANT4 generator...' + generator_thread.start() + print >>sys.stderr, 'Initializing GPU...' gpu_worker = gpu.GPU(options.device) @@ -74,61 +107,29 @@ if __name__ == '__main__': gpu_worker.setup_propagate() gpu_worker.setup_daq(max(detector.pmtids)) - print >>sys.stderr, 'Loading tables in GEANT4...' - # Do easy event to force tables to load, throw away photons - generator.generate_photons(particle_name='e-', - total_energy=1.5, - position=(0,0,0), - direction=(1,0,0)) - - # Create output file f = ROOT.TFile(output_filename, 'RECREATE') ev, T = root.make_tree('T') print >>sys.stderr, 'Starting simulation...' start_sim = time.time() - #### Do the generation and writing in this offset order to ensure - #### that propagation on the GPU overlaps with CPU work to create - #### and save photons. - #### WARNING: THIS MAKES THE CODE LOOK A LITTLE CRAZY. I AM SORRY - - # Do first event - photons = generator.generate_photons(particle_name=options.particle, - total_energy=options.energy, - position=position, - direction=direction) - nphotons = len(photons['pos']) - gpu_worker.load_photons(**photons) - gpu_worker.propagate() - gpu_worker.run_daq() - - for i in xrange(1, options.nevents-1): - # photons for next event while previous event propagates on GPU - photons = generator.generate_photons(particle_name=options.particle, - total_energy=options.energy, - position=position, - direction=direction) + nphotons = 0 + for i in xrange(options.nevents): + photons = generator_thread.queue.get() + + assert len(photons['pos']) > 0, 'GEANT4 generated event with no photons!' + nphotons += len(photons['pos']) - - # this will stop and wait for event to finish - hit_times = gpu_worker.get_hits() - # turn around next event gpu_worker.load_photons(**photons) gpu_worker.propagate() gpu_worker.run_daq() - - # write out this event + hit_times = gpu_worker.get_hits() root.fill_event(T, ev, i-1, position, len(hit_times), hit_times) if i % 10 == 0: print >>sys.stderr, "\rEvent:", i, - # Get results for last event and write it out - hit_times = gpu_worker.get_hits() - root.fill_event(T, ev, options.nevents - 1, position, len(hit_times), hit_times) - end_sim = time.time() print >>sys.stderr, "\rEvent:", options.nevents - 1 @@ -138,3 +139,6 @@ if __name__ == '__main__': print >>sys.stderr, 'Done. %1.1f events/sec, %1.0f photons/sec.' % (options.nevents/(end_sim - start_sim), nphotons/(end_sim - start_sim)) gpu_worker.shutdown() + +if __name__ == '__main__': + main() -- cgit From 09e042b8888342ed8fc7a8c5b05ea1caa47a3842 Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Mon, 8 Aug 2011 15:29:57 -0400 Subject: Modify data structure to hold photon start and stop vertices separately. Add --save-photon-start and --save-photon-stop options to sim.py to save these vertices. Otherwise, only hit information is recorded. --- io/root.C | 58 ++++++++++++++++++++++++++++++++++++---------------------- io/root.py | 3 ++- sim.py | 43 +++++++++++++++++++++++++++++++++++++++++-- 3 files changed, 79 insertions(+), 25 deletions(-) diff --git a/io/root.C b/io/root.C index dfcb3c0..43f4e73 100644 --- a/io/root.C +++ b/io/root.C @@ -9,23 +9,19 @@ struct Photon { TVector3 dir; TVector3 pol; double wavelength; // nm - int state; - int last_triangle; - int last_solid; + int history; + int last_hit_triangle; }; struct MC { std::string particle; TVector3 gen_pos; TVector3 gen_dir; - double gen_ke; + double gen_total_e; - std::vector photon; + std::vector photon_start; + std::vector photon_stop; - void clear() { - particle = "none"; - photon.clear(); - } }; struct Channel { @@ -41,23 +37,42 @@ struct Event { int nhit; std::vector channel; - void clear() { - mc.clear(); - channel.clear(); - } }; -void fill_event(TTree *T, - Event *ev, int event_id, double *gen_pos, int nchannels, - float *channel_times) +void fill_photons(Event *ev, bool start, + unsigned int nphotons, float *pos, float *dir, + float *pol, float *wavelength, float *t0, + int *histories=0, int *last_hit_triangle=0) +{ + std::vector &photons = start ? ev->mc.photon_start : ev->mc.photon_stop; + photons.resize(nphotons); + + for (unsigned int i=0; i < nphotons; i++) { + Photon &photon = photons[i]; + photon.t = t0[i]; + photon.pos.SetXYZ(pos[3*i], pos[3*i + 1], pos[3*i + 2]); + photon.dir.SetXYZ(dir[3*i], dir[3*i + 1], dir[3*i + 2]); + photon.pol.SetXYZ(pol[3*i], pol[3*i + 1], pol[3*i + 2]); + photon.wavelength = wavelength[i]; + if (histories) + photon.history = histories[i]; + else + photon.history = 0; + + if (last_hit_triangle) + photon.last_hit_triangle = last_hit_triangle[i]; + else + photon.last_hit_triangle = -1; + } +} + + +void fill_hits(Event *ev, unsigned int nchannels, float *channel_times) { - ev->clear(); - ev->event_id = event_id; - MC &mc = ev->mc; - mc.gen_pos.SetXYZ(gen_pos[0], gen_pos[1], gen_pos[2]); + ev->channel.resize(0); ev->nhit = 0; Channel ch; - for (int i=0; i < nchannels; i++) { + for (unsigned int i=0; i < nchannels; i++) { if (channel_times[i] < 1e8) { ev->nhit++; ch.channel_id = i; @@ -66,7 +81,6 @@ void fill_event(TTree *T, ev->channel.push_back(ch); } } - T->Fill(); } diff --git a/io/root.py b/io/root.py index d970d1d..e43f5d4 100644 --- a/io/root.py +++ b/io/root.py @@ -5,7 +5,8 @@ ROOT.gROOT.ProcessLine('.L '+os.path.join(os.path.dirname(__file__), 'root.C+g') from ROOT import Event -fill_event = ROOT.fill_event +fill_photons = ROOT.fill_photons +fill_hits = ROOT.fill_hits def make_tree(name, desc=''): '''Create a ROOT tree for holding event information. diff --git a/sim.py b/sim.py index 44fb1dd..0a1f9e6 100755 --- a/sim.py +++ b/sim.py @@ -52,6 +52,29 @@ class GeneratorThread(multiprocessing.Process): self.queue.put(photons) +def write_event(T, ev, event_id, hits, photon_start=None, photon_stop=None): + ev.event_id = event_id + if photon_start is not None: + photons = photon_start + root.fill_photons(ev, True, + len(photons['pos']), + np.ravel(photons['pos']), + np.ravel(photons['dir']), + np.ravel(photons['pol']), + photons['wavelength'], photons['t0']) + if photon_stop is not None: + photons = photon_stop + root.fill_photons(ev, False, + len(photons['pos']), + np.ravel(photons['pos']), + np.ravel(photons['dir']), + np.ravel(photons['pol']), + photons['wavelength'], photons['t0'], + photons['histories'], photons['last_hit_triangles']) + + root.fill_hits(ev, len(hits), hits) + T.Fill() + def main(): parser = optparse.OptionParser('%prog') parser.add_option('-b', type='int', dest='nbits', default=10) @@ -64,6 +87,12 @@ def main(): parser.add_option('--energy', type='float', dest='energy', default=100.0) parser.add_option('--pos', type='string', dest='pos', default='(0,0,0)') parser.add_option('--dir', type='string', dest='dir', default='(1,0,0)') + parser.add_option('--save-photon-start', action='store_true', + dest='save_photon_start', default=False, + help='Save initial photon vertices to disk') + parser.add_option('--save-photon-stop', action='store_true', + dest='save_photon_stop', default=False, + help='Save final photon vertices to disk') options, args = parser.parse_args() if len(args) != 1: @@ -125,8 +154,18 @@ def main(): gpu_worker.propagate() gpu_worker.run_daq() hit_times = gpu_worker.get_hits() - root.fill_event(T, ev, i-1, position, len(hit_times), hit_times) - + + if options.save_photon_start: + photon_start = photons + else: + photon_start = None + if options.save_photon_stop: + photon_stop = gpu_worker.get_photons() + else: + photon_stop = None + + write_event(T, ev, i, hit_times, + photon_start=photon_start, photon_stop=photon_stop) if i % 10 == 0: print >>sys.stderr, "\rEvent:", i, -- cgit From f407488241846ec4c566131a0954689f75b8d1cf Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Mon, 8 Aug 2011 15:50:15 -0400 Subject: Oops. The generator is a process, not a thread. --- sim.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/sim.py b/sim.py index 9e5f9bd..13e53c6 100755 --- a/sim.py +++ b/sim.py @@ -28,7 +28,7 @@ def info(type, value, tb): sys.excepthook = info -class GeneratorThread(multiprocessing.Process): +class GeneratorProcess(multiprocessing.Process): def __init__(self, particle, energy, position, direction, nevents, material): multiprocessing.Process.__init__(self) @@ -114,12 +114,12 @@ def main(): print >>sys.stderr, 'Creating generator...' detector_material = optics.water - generator_thread = GeneratorThread(particle=options.particle, - energy=options.energy, - position=position, - direction=direction, - nevents=options.nevents, - material=detector_material) + generator_thread = GeneratorProcess(particle=options.particle, + energy=options.energy, + position=position, + direction=direction, + nevents=options.nevents, + material=detector_material) print >>sys.stderr, 'WARNING: ASSUMING DETECTOR IS WATER!!' # Do this now so we can get ahead of the photon propagation -- cgit