from Geant4 import *
import g4py.ezgeom
import g4py.NISTmaterials
import g4py.ParticleGun
import pyublas
import numpy as np
import pi0

try:
    import G4chroma
except:
    # Try building the module
    import subprocess
    import sys, os
    module_dir = os.path.split(os.path.realpath(__file__))[0]
    print >>sys.stderr, 'Compiling G4chroma.so...'
    retcode = subprocess.call('g++ -o \'%s/G4chroma.so\' -shared \'%s/G4chroma.cc\' -fPIC `geant4-config --cflags --libs` `python-config --cflags --libs --ldflags` -lboost_python' % (module_dir, module_dir), shell=True)
    assert retcode == 0
    import G4chroma

class G4Generator(object):
    def __init__(self, material, seed=None):
        '''Create generator to produce photons inside the specified material.

           material: chroma.geometry.Material object with density, 
                     composition dict and refractive_index.

                     composition dictionary should be 
                        { element_symbol : fraction_by_weight, ... }.
           seed: Random number generator seed for HepRandom. If None,
                 generator is not seeded.
        '''
        if seed is not None:
            HepRandom.setTheSeed(seed)

        g4py.NISTmaterials.Construct()
        g4py.ezgeom.Construct()
        self.physics_list = G4chroma.ChromaPhysicsList()
        gRunManager.SetUserInitialization(self.physics_list)
        self.particle_gun = g4py.ParticleGun.Construct()
        
        self.world_material = self.create_g4material(material)
        g4py.ezgeom.SetWorldMaterial(self.world_material)

        self.world = g4py.ezgeom.G4EzVolume('world')
        self.world.CreateBoxVolume(self.world_material, 100*m, 100*m, 100*m)
        self.world.PlaceIt(G4ThreeVector(0,0,0))

        self.tracking_action = G4chroma.PhotonTrackingAction()
        gRunManager.SetUserAction(self.tracking_action)
        gRunManager.Initialize()

    def create_g4material(self, material):
        g4material = G4Material('world_material', material.density * g / cm3,
                                len(material.composition))

        # Add elements
        for element_name, element_frac_by_weight in material.composition.items():
            g4material.AddElement(G4Element.GetElement(element_name, True),
                                  element_frac_by_weight)

        # Set index of refraction
        prop_table = G4MaterialPropertiesTable()
        # Reverse entries so they are in ascending energy order rather
        # than wavelength
        energy = list((2*pi*hbarc / (material.refractive_index[::-1,0] * nanometer)).astype(float))
        values = list(material.refractive_index[::-1, 1].astype(float))
        prop_table.AddProperty('RINDEX', energy, values)

        # Load properties
        g4material.SetMaterialPropertiesTable(prop_table)
        return g4material

    def generate_pi0(self, total_energy, position, direction):
        '''Use GEANT4 to generate photons produced by a pi0 decay.

        pi0 event will be generated by running two gammas with the appropriate
        energy and angular distribution.

               total_energy: Total energy of pi0 (incl rest mass) in MeV
               position: 3-tuple of position of particle in global coordinates
               direction: 3-tuple direction vector.  
                          Does not have to be normalized.
        '''
        direction = direction / np.linalg.norm(direction)
        
        cos_theta_rest = np.random.random_sample() * 2 - 1
        theta_rest = np.arccos(cos_theta_rest)
        phi_rest = np.random.random_sample() * 2 * np.pi

        (gamma1_e, gamma1_dir), (gamma2_e, gamma2_dir) = pi0.pi0_decay(energy=total_energy,
                                                                       direction=direction,
                                                                       theta=theta_rest,
                                                                       phi=phi_rest)

        gamma1 = self.generate_photons('gamma', gamma1_e, position, gamma1_dir)
        gamma2 = self.generate_photons('gamma', gamma2_e, position, gamma2_dir)

        photons = {'subtracks' : [ {'name':'gamma', 'pos': position, 't0' : 0.0,
                                    'dir' : gamma1_dir, 'total_e' : gamma1_e },
                                   {'name':'gamma', 'pos': position, 't0' : 0.0,
                                    'dir' : gamma2_dir, 'total_e' : gamma2_e } ]
                    }
        for key in gamma1.keys():
            photons[key] = np.concatenate((gamma1[key], gamma2[key]))
                                                                       
        return photons

    def generate_photons(self, particle_name, total_energy, position, direction):
        '''Use GEANT4 to generate photons produced by the given particle.
           
               particle_name: GEANT4 name of particle.  'e-', 'mu-', etc
               total_energy: Total energy of particle (incl rest mass) in MeV
               position: 3-tuple of position of particle in global coordinates
               direction: 3-tuple direction vector.  
                          Does not have to be normalized.
        '''
        self.particle_gun.SetParticleByName(particle_name)
        self.particle_gun.SetParticleEnergy(total_energy * MeV)
        self.particle_gun.SetParticlePosition(G4ThreeVector(*position))
        self.particle_gun.SetParticleMomentumDirection(G4ThreeVector(*direction).unit())

        self.tracking_action.Clear()

        gRunManager.BeamOn(1)
        n = self.tracking_action.GetNumPhotons()        
        pos = np.zeros(shape=(n,3), dtype=np.float32)
        pos[:,0] = self.tracking_action.GetX()
        pos[:,1] = self.tracking_action.GetY()
        pos[:,2] = self.tracking_action.GetZ()

        dir = np.zeros(shape=(n,3), dtype=np.float32)
        dir[:,0] = self.tracking_action.GetDirX()
        dir[:,1] = self.tracking_action.GetDirY()
        dir[:,2] = self.tracking_action.GetDirZ()

        pol = np.zeros(shape=(n,3), dtype=np.float32)
        pol[:,0] = self.tracking_action.GetPolX()
        pol[:,1] = self.tracking_action.GetPolY()
        pol[:,2] = self.tracking_action.GetPolZ()
        
        wavelength = self.tracking_action.GetWavelength().astype(np.float32)
        t0 = self.tracking_action.GetT0().astype(np.float32)
        return { 'pos' : pos, 
                 'dir' : dir, 
                 'pol' : pol,
                 'wavelength' : wavelength, 
                 't0' : t0 }

if __name__ == '__main__':
    import time
    import optics
    gen = G4Generator(optics.water)
    # prime things
    gen.generate_photons('e-', 1, (0,0,0), (1,0,0))

    start = time.time()   
    n = 0
    for i in xrange(100):
        photons = gen.generate_photons('mu-', 700, (0,0,0), (1,0,0))
        n += len(photons['t0'])
        print photons['pos'][0].min(), photons['pos'][0].max()
    stop = time.time()
    print stop - start, 'sec'
    print n / (stop-start), 'photons/sec'