function detect_limit,bkg,nsig,asrc=asrc,abkg=abkg,bgalt=bgalt,abgalt=abgalt,\$ nsim=nsim,ulsig=ulsig,ulsim=ulsim,gaussy=gaussy,nxbin=nxbin,\$ verbose=verbose, _extra=e ;+ ;function detect_limit ; compute and return the counts upper limit for ; detection of a source at a given significance, ; given the background counts. ; ;description ; compute the cumulative significance of obtaining a specified ; number of counts given the background, and assume that a source ; would be considered detected if the counts were to exceed the ; NSIG threshold. ; ; see Pease, Drake, & Kashyap (2006, ApJ 636, 426) for full description. ; briefly, computes the probability that as many as D counts can be ; observed for a given background b, p(=0, ULSIM[*,0] is the conservative limit that ; is derived from the coadded probability distributions ; that take into account the variations in the background ; gaussy [INPUT] if set, computes the limit corresponding to the ; significance matching the location of the NSIG-sigma ; *intercept* of a Gaussian, rather than matching the total ; area under the curve. ; nxbin [INPUT] number of bins to use in the integration ; * the integration is carrid out over a range of ; 0..5*E(bg) or 20, whichever is greater. ; by default, the number of bins is set by the step size, ; which is set to 1 count, ; -- unless E(bg) < 1, in which case a bin width of 0.05 ; is used by default ; * changing NXBIN does not change the range, but only ; changes the bin size. ; * a hard lower limit of 20 is set -- cannot use a bin ; width larger than 1 count ; verbose [INPUT] controls chatter ; ;subroutines ; LNPOISSON() ; KILROY ; ;history ; vinay kashyap (Apr2004) ; added keyword ULSIG; changed name from PUPLIM to DETECT_LIMIT ; (VK; May2004) ; added keyword NXBIN (VK; Sep2004) ; added keywords BGALT and ABGALT (VK; Dec2004) ; modified output behavior of ULSIM[*,0]; now BGALT can be 0 ; (VK; Mar2005) ; fixed bug for when NXBIN is set; made the numerical precision ; problem at high NSIG explicit by making the code return -1 as ; the UL (VK; Mar2009) ;- message,'Obsolete. Look at ULIM_SIMPLEPOI() instead.',/informational ; usage ok='ok' & np=n_params() & nb=n_elements(bkg) & ns=n_elements(nsig) if np eq 0 then ok='Insufficient parameters' else \$ if nb eq 0 then ok='BKG is undefined' if ok ne 'ok' then begin print,'Usage: ul=detect_limit(bkg,nsig,asrc=asrc,abkg=abkg,nsim=nsim,\$' print,' bgalt=bgalt,abgalt=abgalt,ulsig=ulsig,ulsim=ulsim,\$' print,' /gaussy,nxbin=nxbin,verbose=verbose)' print,' return counts limit for detection in the presence of background' if np ne 0 then message,ok,/informational return,-1L endif ; inputs vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1 ; bb=[bkg[*]] ; msig=fltarr(nb)+3. if ns gt 0 then begin msig[*]=nsig[ns-1L] if ns lt nb then msig[0L:ns-1L]=nsig[*] else msig[*]=nsig[0L:nb-1L] endif ; areas=fltarr(nb)+1. & nas=n_elements(asrc) if nas gt 0 then begin areas[*]=asrc[nas-1L] if nas lt nb then areas[0L:nas-1L]=asrc[*] else areas[*]=asrc[0L:nb-1L] endif ; areab=fltarr(nb)+1. & nab=n_elements(abkg) if nab gt 0 then begin areab[*]=abkg[nab-1L] if nab lt nb then areab[0L:nab-1L]=abkg[*] else areab[*]=abkg[0L:nb-1L] endif ; nalt1=n_elements(bgalt) & nalt2=n_elements(abgalt) if nalt1 gt 0 then begin altbg=[bgalt[*]] altareab=fltarr(nalt1)-1 if nalt2 gt 0 then begin altareab[*]=abgalt[0] altareab[0L:(nalt1 0 if numsim gt 0 then tmpsim=fltarr(numsim) cgauss=errorf(msig[i]/sqrt(2.D)) if keyword_set(gaussy) then begin cgauss=cgauss/2.D + 0.5D if vv gt 5 then print,\$ 'Computing the limit corresponding to a significance of ',cgauss endif for j=0L,numsim do begin ;{for each realization if vv gt 5 then kilroy if areab[i] gt 0 then area_ratio=areas[i]/areab[i] else area_ratio=1. bg=bb[i]*area_ratio if nalt1 gt 0 then begin ob=where(altareab gt 0,mob) & alt_area_ratio=fltarr(nalt1)+1. if mob gt 0 then alt_area_ratio[ob]=areas[i]/altareab[ob] alt_bg=altbg*alt_area_ratio bg=bg+total(alt_bg) endif if j eq 0 then begin nx=0 if keyword_set(nxbin) then begin ;(user defined number of bins if long(nxbin[0]) ge 20L then nx=long(nxbin[0]) endif ;NXBIN) if nx eq 0 then begin ;(default grid dx=1. & if bg lt 1 then dx=0.05 nmax=long(5*bg)>20L & nx=long(nmax/dx)+1 ;if bg lt 1 then nx=400L else nx=long(5*bg)>20L x=findgen(nx+1L)*dx endif else x=findgen(nx+1L) ;default grid) endif if j gt 0 then begin bg=randomu(seed,poisson=bb[i])*area_ratio if nalt1 gt 0 then begin alt_bg_sim=fltarr(nalt1) for k=0L,nalt1-1L do \$ if altbg[k] gt 0 then \$ alt_bg_sim[k]=randomu(seed,poisson=altbg[k])*alt_area_ratio[k] bg=bg+total(alt_bg_sim) endif endif dpr=exp(lnpoisson(x,bg)) if j eq 0 then mcdpr=dpr else mcdpr=mcdpr+dpr cpr=0.D*dpr+dpr[0] & for k=1L,nx do cpr[k]=cpr[k-1L]+dpr[k] ;cpr=0.D*dpr & for k=1L,nx do cpr[k]=cpr[k-1L]+dpr[k-1L] cpr=cpr/max(cpr) tmp=interpol(x,cpr,cgauss) ;Gaussian NSIG equivalent if finite(tmp) eq 0 then tmp=-1 ;a temporary hack if j eq 0 then ul[i]=tmp[0] else tmpsim[j-1L]=tmp[0] endfor ;J=0,NUMSIM} mccpr=0.*mcdpr+mcdpr[0] & for k=1L,nx do mccpr[k]=mccpr[k-1L]+mcdpr[k] mccpr=mccpr/max(mccpr) & tmp=interpol(x,mccpr,cgauss) if finite(tmp) eq 0 then tmp=0. ulsim[i,0]=tmp[0] if numsim gt 0 then ulsim[i,1L:numsim]=tmpsim[*] if numsim gt 1 then ulsig[i]=stddev(tmpsim) if vv gt 10 then plot,x,dpr,xtitle='ct',ytitle='p(ct|bkg)' if vv gt 10 then oplot,tmp*[1,1],[1e-30,1] if vv gt 100 then stop,'HALTing; type .CON to continue' endfor ;I=0,NB-1} return,ul end