#!/usr/bin/env python

class Bin:
    """
    Class representing a histogram bin.
    """

    def __init__(self, low, high, val, isvaldividedbybinwidth):

        self.l = low

        self.h = high

        self.val = val

        if not isvaldividedbybinwidth:
            self.val /= self.GetWidth()

    def GetLow(self):
        return self.l

    def GetHigh(self):
        return self.h

    def GetVal(self):
        return self.val

    def GetWidth(self):
        return self.h - self.l

    def GetMidPoint(self):
        return self.l + GetWidth()/2.

    def Print(self):
        print("* BIN: [%8.3f MeV, %8.3f MeV] -- %10.3f cm^{-2} s^{-1} MeV{-1}" % (self.l, self.h, self.val))

    def GetFlux(self):
        return self.GetVal() * self.GetWidth()

class FluxSpectrum:

    def __init__(self, particle, Es, IFSs, lowerbound, upperbound, dbg = False):
        """
        Constructor.
        """

        self.dbg = dbg

        ## GEANT4 string representing the particle type (e.g. e-, proton).
        self.particle = particle

        ## The lower bound of the spectrum.
        self.lowerbound = lowerbound

        ## The upper bound of the spectrum.
        self.upperbound = upperbound

        if self.dbg:
            print("* DEBUG:")
            print("* DEBUG: Initialising %s FluxSpectra object with lower bound %f MeV, upper bound %f MeV." % (self.particle, self.lowerbound, self.upperbound))
            print("* DEBUG:")

        ## The energy bin values.
        self.Es = Es

        ## The integrated flux spectra values (dictionary).
        self.IFSs = IFSs

        if len(self.Es) != len(self.IFSs):
            print("* ERROR: Mis-match in the energy values and integrated flux spectra values.")

        self.bins = {}

        # Populate the spectrum.
        for i in range(1, len(Es.values())):
            val = self.IFSs[i-1] - self.IFSs[i]
            width = self.Es[i] - self.Es[i-1]

            #if self.dbg:
            #    print("* Creating bin with E_low = %f, E_high = %f" % (self.Es[i-1], self.Es[i]))
            #    print("* Bin integrated     flux: %f - %f = %f cm-2 s-1" % (self.IFSs[i-1], self.IFSs[i], val))
            #    print("* Bin differentiated flux:         = %f cm-2 s-1 MeV-1" % (val/width))

            # Add the bin to the spectrum.
            self.bins[i] = Bin(self.Es[i-1], self.Es[i], val, False)

            if self.dbg:
                self.bins[i].Print()

        total_flux_within_bounds, total_flux = self.GetTotalFlux()

        if self.dbg:

            print("*")
            print("* Total %s flux [%f, %f) MeV = %f (%f) cm^{-2} s^{-1}" % (self.particle, self.lowerbound, self.upperbound, total_flux_within_bounds, total_flux))

    def GetNumBins(self):
        return len(self.bins)

    def GetTotalFlux(self):
        """
        Returns the total flux (cm-2 s-1) within the energy bounds and in total.
        """

        total_flux = 0.0
        total_flux_within_bounds = 0.0

        for i, mybin in self.bins.iteritems():
            total_flux += mybin.GetFlux()

            #print("*--> GetLow() = %f lowerbound = %f" % (mybin.GetLow(), self.lowerbound))

            if mybin.GetLow() >= self.lowerbound and mybin.GetHigh() <= self.upperbound:
                total_flux_within_bounds += mybin.GetFlux()

        return total_flux_within_bounds, total_flux

    def MakeSource(self):

        total_flux_within_bounds, total_flux = self.GetTotalFlux()

        s  = ""
        s += "# Spherical %s source for outer space \n" % (self.particle)
        s += "#----------------------------------\n"
        s += "/gps/source/intensity %f\n" % (total_flux_within_bounds)
        s += "/gps/particle %s\n" % (self.particle)
        s += "/gps/pos/type Surface\n"
        s += "/gps/pos/shape Sphere\n"
        s += "/gps/pos/centre 0. 0. 0. mm\n"
        s += "/gps/pos/radius 49.9 mm\n"
        s += "/gps/pos/confine PseudoDetector_phys\n"
        s += "/gps/ang/type cos\n"
        s += "/gps/ene/type User\n"
        s += "/gps/hist/type energy\n"
        s += "/gps/hist/point %f 0.0\n" % (self.bins[1].GetLow())

        # Now add the user-defined energy histogram bins.
        for i, mybin in self.bins.iteritems():
            if mybin.GetLow() >= self.lowerbound:
                s += "/gps/hist/point %f %f\n" % (mybin.GetHigh(), mybin.GetFlux()) # <- or should it be the int, i.e. multiplied by the bin width?
            else:
                s += "/gps/hist/point %f %f\n" % (mybin.GetHigh(), 0.0) # <- or should it be the int, i.e. multiplied by the bin width?

        s += "#\n"
        #s += "/allpix/lab/setSourceId Hemisphere, %s\n" % (self.particle)


        if self.dbg:
            print s

        return s

    def GetLowerBound(self):
        return self.lowerbound

    def GetUpperBound(self):
        return self.upperbound