#! /usr/bin/env python

#
import math
from random import gauss

index = 0
anum = 1
posn = [0.,0.,0.]
angle = [1.,0.,0.]

mass_mu = 105.6584            #MeV/c**2
#pmean = 28.654                  #MeV/c
#pfrwid = 0.020
pfrwid = 0.010                  #MeV/c, pmean=29.8/(1+2*pfrwid)
pmean = 29.8/(1. + 2.*pfrwid)
#pfrwid = 0.020
prms = pmean*pfrwid
print 'Mean momentum = ', pmean, ', rms width = ', prms

f = open('TRIM_29p922_0p010.DAT','w')

s = 'TRIM.DAT input file for TRIM\r\n'
f.write(s)

s = '\r\n'
f.write(s)
f.write(s)
f.write(s)
f.write(s)
f.write(s)

s = "Incident particle mass %.3f MeV/c**2\r\n" % (mass_mu,)
f.write(s)

s = "Mean momentum %.4f MeV/c, fractional width %.3f, rms width %.4f MeV/c\r\n" % \
    (pmean,pfrwid,prms)
f.write(s)

s = "Event Atom Energy    Depth     Position     -- Atom Direction --\r\n"
f.write(s)
s = "Name   No   (eV)     X (A)   Y (A)   Z (A)  Cos(X)  Cos(Y)  Cos(Z)\r\n"
f.write(s)



while index < 10000:
    sevent = str(index).zfill(5)
    pmom = 100.
    loops = 0
    while pmom > 29.79:
        loops = loops + 1
        pmom = gauss(pmean,prms)
        if (loops >= 100):
            print 'Maximum loop count', loops, 'exceeded'
            print 'at event', index
    energy = 1000000.*(math.sqrt(pmom**2 + mass_mu**2) - mass_mu)  #eV
    line = (sevent,anum,energy,posn[0],posn[1],posn[2],angle[0],angle[1],angle[2])
    s = "%s   %i   %.4e     %.1f     %.1f     %.1f     %.1f     %.1f     %.1f \r\n" % line
    #print s
    f.write(s)
    index = index + 1

#s = r'\004'
#f.write(s)
f.close()