In This Package:

AmC_Co60_Composite.py

Go to the documentation of this file.

#!/usr/bin/env python
#
# Simple composite HepEVT generator
#  Combines the standard Co60.exe generator with a periodic Am-13C neutron
#

def main():
    import sys, getopt
    opts,args = getopt.getopt(sys.argv[1:],"f:s:n:")
    neutronFraction = 0.005
    seed = "0"
    nEvents = 10000
    for opt,arg in opts:
        if opt == "-f":
            neutronFraction = float(arg)
        if opt == "-s":
            seed = arg
        if opt == "-n":
            nEvents = int(arg)

    import os
    co60HepEvts = os.popen("Co60.exe -n %d -seed %s" % (nEvents, seed))
    import random, math
    random.seed(int(seed))
    MeV = 1.0
    GeV = 1000.0
    minNeutronKE = 0.1 * MeV
    for evtIdx in range(nEvents):
        if random.random() < neutronFraction:
            # Throw a neutron
            mass = 939.565 * MeV
            meanKE = 4.25 * MeV
            sigmaKE = 0.6 * MeV
            kinEnergy = random.gauss(meanKE, sigmaKE)
            while kinEnergy < minNeutronKE:
                # Avoid Geant failure for very low energy neutrons
                kinEnergy = random.gauss(meanKE, sigmaKE)
            totEnergy = kinEnergy + mass
            momentum = math.sqrt( (totEnergy+mass) * (totEnergy-mass) )
            cosTheta = random.random()
            sinTheta = math.sqrt( 1 - cosTheta*cosTheta )
            phi = 2*math.pi*random.random()
            mom = [ momentum*sinTheta*math.cos(phi),
                    momentum*sinTheta*math.sin(phi),
                    momentum*cosTheta ]
            print "1"
            print "1     2112 0 0 %f %f %f %f 0.0" % (mom[0]/GeV,
                                                      mom[1]/GeV,
                                                      mom[2]/GeV,
                                                      mass/GeV)
        else:
            # Throw a Co60 decay
            headLine = co60HepEvts.readline().strip()
            while headLine[0] == '#':
                headLine = co60HepEvts.readline().strip()
            nParticles = int(headLine)
            event = [headLine]
            for particle in range(nParticles):
                event.append(co60HepEvts.readline().strip())
            for line in event:
                print line

if __name__ == "__main__":
    main()

---