#!/usr/bin/env python
import sys
import optparse
import time
import multiprocessing

import detectors
import optics
import gpu
import g4gen
from io import root
import numpy as np
import math
import ROOT

def info(type, value, tb):
    if hasattr(sys, 'ps1') or not sys.stderr.isatty():
       # we are in interactive mode or we don't have a tty-like
       # device, so we call the default hook
       sys.__excepthook__(type, value, tb)
    else:
       import traceback, pdb
       # we are NOT in interactive mode, print the exception...
       traceback.print_exception(type, value, tb)
       print
       # ...then start the debugger in post-mortem mode.
       pdb.pm()

sys.excepthook = info

class GeneratorProcess(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 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'], hits['q'], hits['history'])
    T.Fill()

#@profile
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)
    parser.add_option('-n', type='int', dest='nblocks', default=64)

    parser.add_option('--detector', type='string', dest='detector', default='microlbne')
    parser.add_option('--nevents', type='int', dest='nevents', default=100)    
    parser.add_option('--particle', type='string', dest='particle', default='e-')
    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:
        print 'Must specify output filename!'
        sys.exit(1)
    else:
        output_filename = args[0]

    if options.nevents <= 0:
        print '--nevents must be greater than 0!'
        sys.exit(1)

    position = np.array(eval(options.pos), dtype=float)
    direction = np.array(eval(options.dir), dtype=float)
    detector = detectors.find(options.detector)

    print >>sys.stderr, 'Creating BVH for detector "%s" with %d bits...' % (options.detector, options.nbits)
    detector.build(bits=options.nbits)

    print >>sys.stderr, 'Creating generator...'
    detector_material = optics.water
    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
    print >>sys.stderr, 'Starting GEANT4 generator...'
    generator_thread.start()

    print >>sys.stderr, 'Initializing GPU...'
    gpu_worker = gpu.GPU(options.device)
    
    print >>sys.stderr, 'Loading detector onto GPU...'
    gpu_worker.load_geometry(detector)
    
    print >>sys.stderr, 'Initializing random numbers generators...'
    gpu_worker.setup_propagate()
    gpu_worker.setup_daq(max(detector.pmtids))

    # Create output file
    f = ROOT.TFile(output_filename, 'RECREATE')
    ev, T = root.make_tree('T')

    # Set generator info
    ev.mc.particle = options.particle
    ev.mc.gen_pos.SetXYZ(*position)
    ev.mc.gen_dir.SetXYZ(*direction)
    ev.mc.gen_total_e = options.energy

    print >>sys.stderr, 'Starting simulation...'
    start_sim = time.time()
    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'])

        gpu_worker.load_photons(**photons)
        gpu_worker.propagate()
        gpu_worker.run_daq()
        hits = gpu_worker.get_hits()

        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, hits,
                    photon_start=photon_start, photon_stop=photon_stop)
        if i % 10 == 0:
            print >>sys.stderr, "\rEvent:", i,

    end_sim = time.time()
    print >>sys.stderr, "\rEvent:", options.nevents - 1

    T.Write()
    f.Close()

    print >>sys.stderr, 'Done. %1.1f events/sec, %1.0f photons/sec.' % (options.nevents/(end_sim - start_sim), nphotons/(end_sim - start_sim))

if __name__ == '__main__':
    main()