;+ ; NAME: ; POIDEAD ; ; PURPOSE: ; This function calculates the probability of detecting N counts in a ; given measuring time bin for a non-extended dead-time-distorted ; poisson process. The algorithm uses formulae from Muller, ; "Some Formulae for a Dead-Time-Distorted Poisson Process", NIM 117 ; (1974) 401. ; ; CATEGORY: ; Math. ; ; CALLING SEQUENCE: ; ; Result = POIDEAD( Ncount, Mean [, Tinterval, Deadtime ]) ; ; INPUTS: ; Ncount: A scalar or array of number of counts DETECTED in the given ; measuring time bin. ; ; Mean: A scalar or array of TRUE mean number of counts in the given ; measuring time bin due to the ORIGINAL Poissonian process. ; ; OPTIONAL INPUTS: ; ; Tinterval:The length of the measuring time bin, (Default=1). ; ; Deadtime: The length of the non-extended dead time in the SAME units ; as Tinterval, (Default=0). ; ; OUTPUTS: ; An scalar or array is returned giving the probability of detecting ; the counts specified in each element of Ncount and/or Mean. ; ; RESTRICTIONS: ; If Ncount and Mean are both arrays, they must have the same number ; of elements. ; Negative values in the input array will be returned as zeros. ; ; EXAMPLE: ; To determine the probability of detecting 0..100 counts in a 10 msec ; time bin if the true mean is 50 and the dead time is 20 usec, do the ; following: ; ; N = lindgen(101) ; mu = 50.0 ; bin = 10.0 ; tau = 20e-3 ; Prob = POIDEAD( N, mu, bin, tau ) ; plot,N,Prob,psym=10 ; ; MODIFICATION HISTORY: ; Written by: Han Wen, March 1996. ; 15-JUL-1996 Suppress underflow error messages. ; Bugfix: corrected roundoff error in Kappa ; 17-JUL-1996 Accept an array argument for the Mean parameter. ;- function POI_, Ncount, Mu if (Mu eq 0) then return, FLTARR(N_ELEMENTS(Ncount)) lnP = Ncount*ALOG(double(Mu))-Mu-LNGAMMA(Ncount+1) return, EXP(lnP) end function POIDEAD, Ncount, Mean, Tinterval, Deadtime, NSIGMA_CUT=Nsigma_cut ON_ERROR, 2 ; on error, return control to caller NP = N_PARAMS() if (NP lt 2) then message,'Must be called with 3-4 parameters: '+$ 'Ncount, Mean [, Tinterval, Deadtime]' if N_ELEMENTS(Tinterval)eq 0 then Tinterval= 1.0 if N_ELEMENTS(Deadtime) eq 0 then Deadtime = 0.0 Tbin = Tinterval Tau = Deadtime ; Kappa = long(Tbin/Tau) ; Bugfix, 07/15/96 Kappa = CEIL(Tbin/Tau-1) ; for cases when Tbin/Tau = integer Nc = N_ELEMENTS(Ncount) Nm = N_ELEMENTS(Mean) npts = Nc > Nm ct = Ncount mu = Mean > 0.0 if (npts gt 1) and (Nc ne Nm) then begin case 1 of (Nc eq 1) : ct = REPLICATE( Ncount(0), npts ) (Nm eq 1) : mu = REPLICATE( Mean(0)>0.0, npts ) else : message, 'Invalid array dimensions: Ncount, Mean.' endcase endif ; If no dead time, then return normal Poisson distribution ; if (Tau eq 0) then begin if (Nm eq 1) then return, POI_(ct,mu(0)) Prob = FLTARR(npts,/NOZERO) for i=0,npts-1 do Prob(i) = POI_(ct(i),mu(i)) return, Prob endif junk = CHECK_MATH(0,1) ; Suppress math error messages Rho = mu/Tbin Prob = FLTARR(npts) for i=0,npts-1 do begin k = ct(i) s = 0.0 x = Rho(i)*Tau ; Delta_k(t) term ; if (k gt (Kappa+1)) then goto, FILL $ else if (k eq Kappa) then s = (Kappa+1)*(1+x)-mu(i) $ else if (k eq (Kappa+1)) then s = mu(i) - Kappa*(1+x) ; NOTE: We need double-precision in TOTAL and POI_() to avoid ; numerical round-off errors. ; ; U(T_k-1) term ; if (k ge 2) and (k le (Kappa+1)) then begin j = LINDGEN(k-1) mu_eff = Rho(i)*(Tbin-(k-1)*tau) s = s + TOTAL( (k-1-j)*POI_(j,mu_eff),/DOUBLE) endif ; U(T_k) term ; if (k ge 1) and (k le Kappa) then begin j = LINDGEN(k) mu_eff = Rho(i)*(Tbin-k*tau) s = s - 2*TOTAL( (k-j)*POI_(j,mu_eff),/DOUBLE) endif ; U(T_k+1) term ; if (k ge 0) and (k le (Kappa-1)) then begin j = LINDGEN(k+1) mu_eff = Rho(i)*(Tbin-(k+1)*tau) s = s + TOTAL( (k+1-j)*POI_(j,mu_eff),/DOUBLE) endif ; Scale by lambda = 1/(1+x) factor ; s = s/(1+x) > 0.0 FILL: Prob(i) = s endfor if (npts eq 1) then Prob = Prob(0) junk = CHECK_MATH(0,0) ; Turn err msgs back on return, Prob end