In This Package:

decay.py

Go to the documentation of this file.

#!/usr/bin/env python

from data import units

class AlphaDecay(object):
    '''Nuclear decay by alpha emission'''

    def __init__(self,qvalue,parentZ=0):
        '''Create an alpha decay of the given qvalue, and optional
        parentZ (must be set eventually)'''
        self.qvalue = qvalue
        self.parentZ = parentZ
        return

    def __str__(self):
        return 'alpha: Z=%d Q=%.4f MeV'%(self.parentZ,self.qvalue/units.MeV)

    def energy(self):
        return (self.parentZ-4)/(self.parentZ*self.qvalue)

    pass


class BetaDecay(object):
    '''Nuclear decay by beta emission'''
    
    def __init__(self,qvalue,parentZ):
        'Create a beta decay of the given qvalue and parentZ'
        self.qvalue = float(qvalue)
        self.kine = None
        self.doff = True        # for testing w/out fermi function
        self.onlyff = False     # for testing w/only fermi function

        if parentZ < 0:         # beta+ decay
            self.parentZ = -parentZ
            self.daughterZ = self.parentZ - 1
            self.betaSign = +1
            self.endpoint = self.qvalue - 2.0*units.electron_mass_c2
        else:                   # beta- decay
            self.parentZ = parentZ
            self.daughterZ = self.parentZ + 1
            self.betaSign = -1
            self.endpoint = self.qvalue

        self._normalize()
        return
    def _normalize(self):
        steps = 1000
        dx = self.endpoint/steps
        lo = dx/2.0
        self.norm = 1.0
        self.norm = sum([self.dnde(i*dx+lo)*dx for i in range(steps)])
        self.maximum = max([self.dnde(i*dx+lo) for i in range(steps)])
        #print 'norm=%f max=%f %s'%(self.norm,self.maximum,self)
        return

    def __str__(self):
        pm = '+'
        if self.betaSign < 0: pm = '-'
        return 'beta%s: Z=%d(->%d) Q=%.3f MeV'%(pm,self.parentZ,self.daughterZ,self.qvalue/units.MeV)

    def energy(self):
        if self.kine is None:
            self.kine = self.sampleEnergy()
        return self.kine

    def sampleEnergy(self):
        from random import Random
        u = Random().uniform
        while True:
            T = u(0.0,self.endpoint)
            P = u(0.0,self.maximum)
            dnde = self.dnde(T)
            assert(type(T) is float)
            assert(type(P) is float)
            assert(type(dnde) is float)
            if P <= dnde: return T
            continue
        return None

    def dnde(self,T):
        if T > self.endpoint: return 0.0
        if T < 0.0: return 0.0

        if self.onlyff:
            ret = self.fermi_function(T)
        else:
            ret = self.dnde_noff(T)
            if self.doff: ret *= self.fermi_function(T)
        return ret/self.norm

    def dnde_noff(self,T):
        '''Return the unormalized dN/dE at the given beta kinetic
        energy and ignoring nuclear Coulomb field'''
        import math
        W = self.endpoint/units.electron_mass_c2 + 1.0
        E = T/units.electron_mass_c2 + 1.0
        return math.sqrt(E**2-1.0) * (W-E)**2 * E
        
    def fermi_function(self,T):
        "Return the unormalized Fermi function value for given kinetic energy"
        import math
        E = T/units.electron_mass_c2 + 1.0
        P = math.sqrt(E*E-1.0)
        U = -1*self.betaSign*self.daughterZ/137.0
        S = math.sqrt(1.0 - U*U) - 1.0
        Y = 2.0*math.pi*U*E/P
        A1 = U*U * E*E + P*P/4.0
        A2 = math.fabs(Y/(1.0-math.exp(-Y)))
        return math.pow(A1,S)*A2

class GammaDecay(object):
    '''A decay by radiating a photon.'''

    def __init__(self,energy):
        self._energy = energy
        return

    def energy(self): return self._energy

    def __str__(self):
        from data import units
        return 'gamma: %.4f MeV'%(self.energy()/units.MeV)

    pass


def plot_64Cu():
    # beta-
    bm = BetaDecay(0.5794*units.MeV,+29)
    print bm
    Tm = map(lambda x: bm.endpoint/100*(x+0.5),range(99))
    Sm = map(bm.dnde,Tm)

    # beta+
    bp = BetaDecay(1.6751*units.MeV,-29)
    print bp
    Tp = map(lambda x: bp.endpoint/100*(x+0.5),range(99))
    Sp = map(bp.dnde,Tp)


    import matplotlib.pyplot as plt
    plt.figure(1)
    plt.subplot(211)
    plt.plot(Tp,Sp,'r')
    plt.subplot(212)
    plt.plot(Tm,Sm,'b')
    plt.show()
    return (bp,bm)


if '__main__' == __name__:
    plot_64Cu()

---