diff options
| -rwxr-xr-x | benchmark.py | 79 | ||||
| -rw-r--r-- | generator/g4gen.py | 2 | ||||
| -rw-r--r-- | generator/vertex.py | 21 |
3 files changed, 92 insertions, 10 deletions
diff --git a/benchmark.py b/benchmark.py index 5a344bf..8c75d35 100755 --- a/benchmark.py +++ b/benchmark.py @@ -148,6 +148,82 @@ def pdf(gpu_geometry, max_pmt_id, npdfs=10, nevents=100, nreps=16, ndaq=1, return nevents*nreps*ndaq/ufloat((np.mean(run_times),np.std(run_times))) +def pdf_eval(gpu_geometry, max_pmt_id, npdfs=10, nevents=25, nreps=16, ndaq=128, + nthreads_per_block=64, max_blocks=1024): + """ + Returns the average number of 100 MeV events that can be + histogrammed per second. + + Args: + - gpu_instance, chroma.gpu.GPU + The GPU instance passed to the GPUGeometry constructor. + - max_pmt_id, int + The channel number of the highest PMT + - npdfs, int + The number of pdf generations to average. + - nevents, int + The number of 100 MeV events to generate for each PDF. + - nreps, int + The number of times to propagate each event and add to PDF + - ndaq, int + The number of times to run the DAQ simulation on the propagated + event and add it to the PDF. + """ + rng_states = gpu.get_rng_states(nthreads_per_block*max_blocks) + + # Make data event + data_ev = g4generator.generate_events(itertools.islice(generator.vertex.constant_particle_gun('e-', (0,0,0), + (1,0,0), 100), + 1)).next() + gpu_photons = gpu.GPUPhotons(data_ev.photons_beg) + + gpu_photons.propagate(gpu_geometry, rng_states, + nthreads_per_block, max_blocks) + gpu_daq = gpu.GPUDaq(gpu_geometry, max_pmt_id) + data_ev_channels = gpu_daq.acquire(gpu_photons, rng_states, nthreads_per_block, max_blocks).get() + + # Setup PDF evaluation + gpu_pdf = gpu.GPUPDF() + gpu_pdf.setup_pdf_eval(data_ev_channels.hit, + data_ev_channels.t, + data_ev_channels.q, + 0.05, + (-0.5, 999.5), + 1.0, + (-0.5, 20), + min_bin_content=20, + time_only=True) + + run_times = [] + for i in tools.progress(range(npdfs)): + t0 = time.time() + gpu_pdf.clear_pdf_eval() + + vertex_gen = generator.vertex.constant_particle_gun('e-', (0,0,0), + (1,0,0), 100) + vertex_iter = itertools.islice(vertex_gen, nevents) + + for ev in g4generator.generate_events(vertex_iter): + gpu_photons = gpu.GPUPhotons(ev.photons_beg, ncopies=nreps) + + gpu_photons.propagate(gpu_geometry, rng_states, + nthreads_per_block, max_blocks) + for gpu_photon_slice in gpu_photons.iterate_copies(): + for idaq in xrange(ndaq): + gpu_channels = gpu_daq.acquire(gpu_photon_slice, rng_states, + nthreads_per_block, max_blocks) + gpu_pdf.accumulate_pdf_eval(gpu_channels, nthreads_per_block) + + cuda.Context.get_current().synchronize() + elapsed = time.time() - t0 + + if i > 0: + # first kernel call incurs some driver overhead + run_times.append(elapsed) + + return nevents*nreps*ndaq/ufloat((np.mean(run_times),np.std(run_times))) + + if __name__ == '__main__': from chroma import detectors import gc @@ -171,4 +247,7 @@ if __name__ == '__main__': print '%s 100 MeV events histogrammed/s' % \ tools.ufloat_to_str(pdf(gpu_geometry, max(lbne.pmtids))) + print '%s 100 MeV events/s accumulated in PDF evaluation data structure (100 GEANT4 x 16 Chroma x 128 DAQ)' % \ + tools.ufloat_to_str(pdf_eval(gpu_geometry, max(lbne.pmtids))) + context.pop() diff --git a/generator/g4gen.py b/generator/g4gen.py index 3836897..91ddc18 100644 --- a/generator/g4gen.py +++ b/generator/g4gen.py @@ -120,6 +120,8 @@ class G4Generator(object): self.particle_gun.SetParticleEnergy(total_energy) self.particle_gun.SetParticlePosition(G4ThreeVector(*vertex.pos)*m) self.particle_gun.SetParticleMomentumDirection(G4ThreeVector(*vertex.dir).unit()) + self.particle_gun.SetParticleTime(vertex.t0*s) + if vertex.pol is not None: self.particle_gun.SetParticlePolarization(G4ThreeVector(*vertex.pol).unit()) diff --git a/generator/vertex.py b/generator/vertex.py index b4e2d23..0ec4b31 100644 --- a/generator/vertex.py +++ b/generator/vertex.py @@ -40,17 +40,18 @@ def flat(e_lo, e_hi): # vertex generators -def particle_gun(particle_name_iter, pos_iter, dir_iter, ke_iter, start_id=0): - for i, particle_name, pos, dir, ke in izip(count(start_id), particle_name_iter, pos_iter, dir_iter, ke_iter): +def particle_gun(particle_name_iter, pos_iter, dir_iter, ke_iter, + t0_iter=constant(0.0), start_id=0): + for i, particle_name, pos, dir, ke, t0 in izip(count(start_id), particle_name_iter, pos_iter, dir_iter, ke_iter, t0_iter): dir /= np.linalg.norm(dir) - vertex = event.Vertex(particle_name, pos, dir, ke) + vertex = event.Vertex(particle_name, pos, dir, ke, t0=t0) ev_vertex = event.Event(i, vertex, [vertex]) yield ev_vertex -def pi0_gun(pos_iter, dir_iter, ke_iter, start_id=0): - for i, pos, dir, ke in izip(count(start_id), pos_iter, dir_iter, ke_iter): +def pi0_gun(pos_iter, dir_iter, ke_iter, t0_iter=constant(0.0), start_id=0): + for i, pos, dir, ke, t0 in izip(count(start_id), pos_iter, dir_iter, ke_iter, t0_iter): dir /= np.linalg.norm(dir) - primary_vertex = event.Vertex('pi0', pos, dir, ke) + primary_vertex = event.Vertex('pi0', pos, dir, ke, t0=t0) cos_theta_rest = np.random.random_sample() * 2 - 1 theta_rest = np.arccos(cos_theta_rest) @@ -59,15 +60,15 @@ def pi0_gun(pos_iter, dir_iter, ke_iter, start_id=0): (gamma1_e, gamma1_dir), (gamma2_e, gamma2_dir) = \ pi0_decay(ke+134.9766, dir, theta_rest, phi_rest) - gamma1_vertex = event.Vertex('gamma', pos, gamma1_dir, gamma1_e) - gamma2_vertex = event.Vertex('gamma', pos, gamma2_dir, gamma2_e) + gamma1_vertex = event.Vertex('gamma', pos, gamma1_dir, gamma1_e, t0=t0) + gamma2_vertex = event.Vertex('gamma', pos, gamma2_dir, gamma2_e, t0=t0) ev_vertex = event.Event(i, primary_vertex, [gamma1_vertex, gamma2_vertex]) yield ev_vertex -def constant_particle_gun(particle_name, pos, dir, ke, start_id=0): +def constant_particle_gun(particle_name, pos, dir, ke, t0=0.0, start_id=0): '''Convenience wrapper around particle gun that assumes all arguments are constants, rather than generators.''' - return particle_gun(constant(particle_name), constant(pos), constant(dir), constant(ke), start_id=start_id) + return particle_gun(constant(particle_name), constant(pos), constant(dir), constant(ke), constant(t0), start_id=start_id) |