In [1]:
# import some useful utils
%matplotlib notebook
from matplotlib import pyplot as plt
import numpy as np
import histogram.hdf as hh, histogram as H


/home/lj7/anaconda2/envs/dev-mph/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

In [3]:
def getDOS(eventnxs, Emin=-100, Emax=100, dE=1.,
           Qmin=0, Qmax=15., dQ=0.1, T=300, Ecutoff=50., 
           elastic_E_cutoff=(-20., 7), M=50.94,
           C_ms=0.3, Ei=116.446, workdir='work'):
    import os
    # reduce
    cmd = "mcvine instruments arcs nxs reduce %(eventnxs)s --out=iqe.nxs"
    cmd += " --eaxis %(Emin)s %(Emax)s %(dE)s --qaxis %(Qmin)s %(Qmax)s %(dQ)s"
    cmd = cmd % locals()
    if os.system(cmd):
        raise RuntimeError("%s failed" % cmd)
    # to histogram
    cmd = "mcvine mantid extract_iqe iqe.nxs iqe.h5"
    if os.system(cmd):
        raise RuntimeError("%s failed" % cmd)
    # to DOS
    iqehist = hh.load("./v30585-iqe.h5") 
    # interpolate data
    from multiphonon.sqe import interp
    # probably don't need this line
    newiqe = interp(iqehist, newE = np.arange(Emin, Emax, dE))
    # save interpolated data
    hh.dump(newiqe, 'iqe-interped.h5')
    # create processing engine
    from multiphonon.backward import sqe2dos
    iterdos = sqe2dos.sqe2dos(
      newiqe, T=T, Ecutoff=Ecutoff, 
      elastic_E_cutoff=elastic_E_cutoff, M=M,
      C_ms=C_ms, Ei=Ei, workdir='work')
    doslist = list(iterdos)
    return doslist

In [4]:
eventnxs="/SNS/ARCS/2012_2_18_CAL/0/30585/NeXus/ARCS_30585_event.nxs"

In [5]:
doslist = getDOS(eventnxs, Emin=-100, Emax=100, dE=1.,
       Qmin=0, Qmax=15., dQ=0.1, T=300, Ecutoff=50., 
       elastic_E_cutoff=(-20., 7), M=50.94,
       C_ms=0.3, Ei=116.446, workdir='work')


/home/lj7/anaconda2/envs/dev-mph/lib/python2.7/site-packages/histogram/hdf/Loader.py:129: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if 'storage' in list(dataGroup): # this uses the 'storage' convention
/home/lj7/dv/sns-chops/multiphonon/multiphonon/forward/phonon.py:227: RuntimeWarning: divide by zero encountered in divide
  return np.cosh(x)/np.sinh(x)
/home/lj7/dv/sns-chops/multiphonon/multiphonon/forward/phonon.py:242: RuntimeWarning: invalid value encountered in multiply
  f = coth(beta * E/2.) * g/E
/home/lj7/dv/sns-chops/multiphonon/multiphonon/forward/phonon.py:200: RuntimeWarning: divide by zero encountered in divide
  t = 1./(1-np.exp(-E*beta)) # XXX
/home/lj7/dv/sns-chops/multiphonon/multiphonon/forward/phonon.py:201: RuntimeWarning: invalid value encountered in divide
  t = g/(E*g0)*t
/home/lj7/dv/sns-chops/multiphonon/multiphonon/backward/singlephonon_sqe2dos.py:63: RuntimeWarning: divide by zero encountered in divide
  dos = initdos * (expse/simse)
/home/lj7/dv/sns-chops/multiphonon/multiphonon/backward/singlephonon_sqe2dos.py:63: RuntimeWarning: invalid value encountered in multiply
  dos = initdos * (expse/simse)

In [ ]: