Uses lmfit to fit a "Double LorGauss" distribution to a provided histogram.

    Uses a grid of starting guesses to try to avoid local minima.

    Parameters
    ----------
    bins, h : array-like
        Bins and heights of a histogram, as returned by, e.g., `np.histogram`.

    Ntry : int, optional
        Spacing of grid for starting guesses.  Will try `Ntry**2` different
        initial values of the "Gaussian strength" parameters `G1` and `G2`.

    Returns
    -------
    parameters : tuple
        Parameters of best-fit "double LorGauss" distribution.

    Raises
    ------
    ImportError
        If the lmfit module is not available.

def fit_double_lorgauss(bins,h,Ntry=5):
    """Uses lmfit to fit a "Double LorGauss" distribution to a provided histogram.

    Uses a grid of starting guesses to try to avoid local minima.

    Parameters
    ----------
    bins, h : array-like
        Bins and heights of a histogram, as returned by, e.g., `np.histogram`.

    Ntry : int, optional
        Spacing of grid for starting guesses.  Will try `Ntry**2` different
        initial values of the "Gaussian strength" parameters `G1` and `G2`.

    Returns
    -------
    parameters : tuple
        Parameters of best-fit "double LorGauss" distribution.

    Raises
    ------
    ImportError
        If the lmfit module is not available.
    """
    try:
        from lmfit import minimize, Parameters, Parameter, report_fit
    except ImportError:
        raise ImportError('you need lmfit to use this function.')
        
    #make sure histogram is normalized
    h /= np.trapz(h,bins)

    #zero-pad the ends of the distribution to keep fits positive
    N = len(bins)
    dbin = (bins[1:]-bins[:-1]).mean()
    newbins = np.concatenate((np.linspace(bins.min() - N/10*dbin,bins.min(),N/10),
                             bins,
                             np.linspace(bins.max(),bins.max() + N/10*dbin,N/10)))
    newh = np.concatenate((np.zeros(N/10),h,np.zeros(N/10)))

    
    mu0 = bins[np.argmax(newh)]
    sig0 = abs(mu0 - newbins[np.argmin(np.absolute(newh - 0.5*newh.max()))])

    def set_params(G1,G2):
        params = Parameters()
        params.add('mu',value=mu0)
        params.add('sig1',value=sig0)
        params.add('sig2',value=sig0)
        params.add('gam1',value=sig0/10)
        params.add('gam2',value=sig0/10)
        params.add('G1',value=G1)
        params.add('G2',value=G2)
        return params

    sum_devsq_best = np.inf
    outkeep = None
    for G1 in np.linspace(0.1,1.9,Ntry):
        for G2 in np.linspace(0.1,1.9,Ntry):
            params = set_params(G1,G2)
        
            def residual(ps):
                pars = (params['mu'].value,
                        params['sig1'].value,
                        params['sig2'].value,
                        params['gam1'].value,
                        params['gam2'].value,
                        params['G1'].value,
                        params['G2'].value)
                hmodel = double_lorgauss(newbins,pars)
                return newh-hmodel

            out = minimize(residual,params)
            pars = (out.params['mu'].value,out.params['sig1'].value,
                    out.params['sig2'].value,out.params['gam1'].value,
                    out.params['gam2'].value,out.params['G1'].value,
                    out.params['G2'].value)

            sum_devsq = ((newh - double_lorgauss(newbins,pars))**2).sum()
            #print 'devs = %.1f; initial guesses for G1, G2; %.1f, %.1f' % (sum_devsq,G1, G2)
            if sum_devsq < sum_devsq_best:
                sum_devsq_best = sum_devsq
                outkeep = out
            

    return (outkeep.params['mu'].value,abs(outkeep.params['sig1'].value),
            abs(outkeep.params['sig2'].value),abs(outkeep.params['gam1'].value),
            abs(outkeep.params['gam2'].value),abs(outkeep.params['G1'].value),
            abs(outkeep.params['G2'].value))