In This Package:

beta.py

Go to the documentation of this file.

#!/usr/bin/env python



class Spectra(object):
    '''
    Look at Beta decay Spectra
    '''

    def __init__(self,z,a):
        '''
        Create Spectra for given ground state isotope of given Z and A.
        '''
        import random
        self.random = random.Random()

        from g4data import G4RadData
        rd = G4RadData()
        rd = G4RadData()
        self.nuc = rd.nuclideData(z,a)
        state = self.nuc.excitedState()

        betas = []
        for dk in state.decays:
            if dk.type[0:4] == 'Beta':
                betas.append(dk)
                continue
            continue

        if not betas:
            print 'Warning: no beta decays for: %s'%nuc
            return

        self.betas = betas
        return
    
    def modeledSpectra(self,nbins=100):
        total_fraction = sum([dk.fraction for dk in self.betas])
        maximum = max([dk.radiation.endpoint for dk in self.betas])

        ret = []
        from numpy import array
        total = array([0.0 for x in range(nbins)])
        bins = [(x+0.5)*maximum/nbins for x in range(nbins)]
        for dk in self.betas:
            ratio = dk.fraction/total_fraction
            Sm = []
            Sm = map(dk.radiation.dnde,bins)
            scaled = [ratio*x for x in Sm]
            total += scaled
            ret.append(scaled)
            continue
        return (ret,bins,total)

    def sampledSpectra(self,ndecays=1000,nbins=100):
        '''Return a list of tuples (decay.BetaDecay,numpy.histogram).
        Each histogram has the given number of bins and gives the
        sampled energy spectrum.  The total number of entries of all
        histograms is as given.  The branching fractions of the decay
        are used to divide up the total number of decays when sampling
        each branch.  Histograms are normalized so their sum is 1 and
        keeping the branching ratios intact'''
        from numpy import histogram, array

        total_fraction = float(sum([dk.fraction for dk in self.betas]))
        maximum = max([dk.radiation.endpoint for dk in self.betas])
        ret = []
        
        binwidth = maximum/nbins

        for dk in self.betas:
            ratio = dk.fraction/total_fraction
            energies = []
            dk_number = ndecays*ratio
            for count in range(0,int(dk_number)):
                energies.append(dk.radiation.sampleEnergy())
            h = histogram(energies,nbins,(0,maximum))
            norm = ratio/(binwidth*sum(h[0]))
            scaled = array([ bin*norm for bin in h[0] ])
            ret.append((dk,(scaled,h[1])))
            continue
        return ret

        
    pass

def Bi210():
    return plot(83,210)

def plot(z,a,nbins=100):
    spec = Spectra(z,a)

    dk_h = spec.sampledSpectra(ndecays=10000,nbins=nbins)
    
    from numpy import array

    firstDK,firstH = dk_h[0]
    total = array([0.0 for x in range(len(firstH[0]))])
    for dk,h in dk_h:
        total += h[0]
        continue

    mods,bins,modTotal = spec.modeledSpectra(nbins=nbins)

    import matplotlib.pyplot as plt
    for dk,h in dk_h:
        plt.plot(h[1],h[0])
        continue
    plt.plot(firstH[1],total)

    plt.plot(bins,modTotal)
    for mod in mods:
        plt.plot(bins,mod)
        continue

    import elements
    plt.title('%d%s'%(a,elements.symbol(z)))

    plt.show()

---