;+ ; NAME: ; PROB2EXP ; ; PURPOSE: ; This function determines the probability function of the ; waiting times between successive events from source 1 in the ; presence of possible events from source 2 arriving between source 1 ; events. Both sources are assumed Poisson with a respective ; non-extended dead time. ; ; CATEGORY: ; Math. ; ; CALLING SEQUENCE: ; ; Result = FWT2EXP( Wtmax, Mean1, Mean2, Offset1, Offset2 ) ; ; INPUTS: ; Wt: An array of waiting times or the maximum waiting time to ; determine the probability density function at. ; ; Mean1: The average waiting time for source 1 with no dead time. ; ; Mean2: The average waiting time for source 2 with no dead time. ; ; Offset1: The dead time of source 1. ; ; Offset2: The dead time of source 2. ; ; OPTIONAL INPUT KEYWORD PARAMETERS: ; ; BIN: Time resolution of the waiting times, (1=Default). ; ; SUM: Set this keyword to return the total probability ; function, see OUTPUTS below, (0=Default). ; ; OUTPUTS: ; This function returns the probabilities at the waiting times ; specified by the Wt parameter if Wt is an array. If Wt is a scalar, ; the waiting times chosen are [0, BIN, 2*BIN, 3*BIN,... N*BIN < Wt]. ; Normally, this function will return a 2-D array, (nWT, nQ) where ; the first dimension corresponds to the waiting times and the second ; dimension specifies the number of events which occurred from ; source 2 in a given waiting time. If the SUM keyword is set, then ; the total probability is returned (i.e. TOTAL((nWT, nQ),2)). ; ; The probability density function follows the formula: ; ; psd = SUM(f_q) ; ; p1 = 1/Mean1 ; p2 = 1/Mean2 ; q = Number of events from source 2 occurring in the ; waiting time, T ; T' = T - (Offset1 + q*Offset2) ; ; f(p,T) = p*exp(-p*T), Waiting time distribution for Poisson source ; P(n,u) = u^n*exp(-u)/n!, Prob. of n counts for Poisson source ; with mean, u. ; ; f_q = f(p1,T')*P(q,p2*T'), T > Offset1 + q*Offset2 ; ; To calculate the probability, one must take into account that events ; can fall anywhere within the time resolution, BIN. ; ; EXAMPLE: ; Given data with 8 microsec resolution, with p1 = 500 Hz, p2 = 200 Hz, ; t1 = 20 microsec, and t2 = 1 ms, to determine the Probabilities from ; 0 to 1000 microsec, we do the following: ; ; Bin = 8e-6 ; p1 = 500. ; p2 = 200. ; t1 = 20e-6. ; t2 = 1e-3 ; wtmax= 1000.e-6/Bin ; nbin = wtmax/Bin ; wt = Bin*FINDGEN(nbin)/(nbin-1.) ; psd = PROB2EXP( wtmax, 1/p1, 1/p2, t1, t2, BIN=Bin ) ; ; MODIFICATION HISTORY: ; Written by: Han Wen, September 1996. ; 30-SEP-1996 Misc. bug fixes. Accept array for Wt parameter. ;- function POI_, Ncount, Mu lnP = 0d0*Mu h = WHERE(Mu gt 0,ngt0) if (ngt0 gt 0) then \$ lnP(h) = Ncount(h)*ALOG(DOUBLE(Mu(h)))-Mu(h)-LNGAMMA(Ncount(h)+1) return, EXP(lnP) end function PROB2EXP, Wt, Mean1, Mean2, Offset1, Offset2, BIN=Bin, SUM=Sum ; Check parameters NP = N_PARAMS() if (NP ne 5) then message, 'Must be called with 5 parameters: '+ \$ 'Wt, Mean1, Mean2, Offset1, Offset2' if (N_ELEMENTS(Bin) eq 0) then Bin=1 nWT = N_ELEMENTS(Wt) Wtmax= MAX(Wt) b = Bin t1 = Offset1 t2 = Offset2 p1 = 1./Mean1 p2 = 1./Mean2 p = p1 + p2 mu = p*b nmax = LONG(1.*WTmax/Bin) if (nmax lt 1 ) then message, 'BIN keyword set too large.' if (nWT eq 1) then Prob = DBLARR( nmax+1 ) else Prob = DBLARR( nWT ) if (WTmax lt t1) then return, Prob ; Determine max. number of events from Source 2 eps = 1e-25 qmax = LONG((WTmax - 1.*t1)/t2) nq = qmax + 1 if (KEYWORD_SET(Sum) eq 0) then Prob = Prob # REPLICATE(1.,nq) K1 = LONG(t1/b) if (nWT gt 1) then begin h1 = WHERE( Wt eq K1*b ) h1p1 = WHERE( Wt eq (K1+1)*b ) endif else begin h1 = K1 h1p1 = K1 + 1 endelse ; Probabilities of waiting times around the source 1 offset if (nmax lt K1) then goto, DONE if (h1(0) ne -1) then begin mu_k = p*(K1*b -(t1 + 0*t2)) p0q_k = (p1/(p*mu))*(poi_(0,mu_k+mu) + mu + mu_k - 1) Prob(h1,0)= p0q_k endif if (nmax lt (K1+1)) then goto, DONE if (h1p1(0) ne -1) then begin mu_kp1 = p*((K1+1)*b - (t1 + 0*t2)) p0q_kp1 = (p1/(p*mu)) \$ *(poi_(0,mu_kp1 + mu) -2*poi_(0,mu_kp1) \$ + mu - mu_kp1 + 1) Prob(h1p1,0)= p0q_kp1 endif if (nmax lt (K1+2)) then goto, DONE nbin = nmax - (K1+1) n0 = K1+2 n = (LINDGEN(nbin) + n0) ; Probabilities of waiting times from source 1 and q events from source 2 for q=0L,nq-1 do begin K2 = LONG((t1 + q*t2)/b) nmin = K2 + 2 if (nWT gt 1) then begin h = WHERE( (Wt ge nmin*b) and (Wt le nmax*b), nbin ) nm = LONG(Wt(h)/b) # REPLICATE(1,q+1) endif else begin nbin = nmax - nmin + 1 nm = n(nmin-n0:nmax-n0) # REPLICATE(1,q+1) endelse z = REPLICATE(1,nbin) # LINDGEN(q+1) mu_n = p*(nm*b -(t1 + q*t2)) pNq = (((p1/p)*(p2/p)^q)/mu) \$ *((z + 1)*( poi_(q-z,mu_n-mu) \$ +poi_(q-z,mu_n+mu) \$ -2*poi_(q-z,mu_n ))) if (q gt 0) then pNq = TOTAL(pNq,2) i = q*(1 - KEYWORD_SET(Sum)) if (nWT gt 1) then Prob(h,i) = Prob(h,i) + pNq \$ else Prob(nmin:nmax) = Prob(nmin:nmax) + pNq if (MAX(pNq) lt eps) then goto, DONE endfor DONE: return, Prob end