import numpy as np
from scipy.special import factorial

def xenon_resolution(s1,nu):
'''See Eqs. (14) and (15) of 1103.0303'''
     deltan=20
     nmin=max(1,np.rint(s1)-deltan)
     nmax=np.rint(s1)+deltan
     sigma_pmt=0.5
     p=0.
     for n in np.arange(nmin,nmax+1,1.):
          sigma=np.sqrt(n)*sigma_pmt
          p+=gaussian(s1,n,sigma)*poisson(n,nu)
     return p

def gaussian(x,mu,sigma):
     return 1./(np.sqrt(2.*np.pi)*sigma)*np.exp(-(x-mu)**2/(2.*sigma**2))

def poisson(n,l):
     if l<100:
          return l**n/float(factorial(n,exact=False))*exp(-l)
     else:
          return gaussian(n,l,np.sqrt(l))