Energy resolution (dual phase xenon detectors)import numpy as npfrom 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))*np.exp(-l) else: return gaussian(n,l,np.sqrt(l)) |