;+ ; NAME: ; FIT2WT ; ; PURPOSE: ; This routine fits the time interval distribution for the ; dead time and the average count rate of the observed source 1 ; and vetoed source 2. ; ; CATEGORY: ; Math ; ; CALLING SEQUENCE: ; ; FIT2WT, Ts, Minima ; ; INPUTS: ; Ts: An array of arrival times. ; ; OPTIONAL INPUT KEYWORD PARAMETERS: ; ; BINSIZE: The waiting time bin size to form the waiting time distribution, ; (Default = 1). ; ; MAX: Display the waiting time distribution up to this maximum waiting ; time, (Default = 1000*Binsize). If the RESTRICT keyword is NOT ; set then this value will also be the maximum value that the ; waiting time distribution will be fitted up to. ; ; RESTRICT: Set this keyword to fit a restricted range of data, (0=Default). ; ; VERBOSE: Set this keyword to display the current fit parameters before ; each iteration, (0=Default). ; ; OUTPUTS: ; Minima: Structure holding the results of the fit. Tags are: ; ; param : Structure holding the fitted parameters with tags: ; rate1 : Fitted count rate of source 1 without dead time. ; rate2 : Fitted count rate of source 2 without dead time. ; offset1: Fitted dead time of source 1. ; offset2: Fitted dead time of source 2. ; ; sigma : Structure holding the 1 sigma error of the fitted ; parameters with tags: ; rate1 : 1 sigma error on the fitted count rate1 parameter. ; rate2 : 1 sigma error on the fitted count rate2 parameter. ; offset1: 1 sigma error on the offset1 parameter. ; offset2: 1 sigma error on the offset2 parameter. ; ; wtbin: Waiting time bins of the fitted distribution. ; noccur: Number of occurrences at each waiting time bin. ; prob: Probability distribution of the fit. ; nwt: Total number of waiting times (normalization). ; chi2: Total chi-squared of the fit. ; dof: Number of degrees of freedom of the fit. ; CL: Associated confidence level of the chi-squared. ; ; ; COMMON BLOCKS: ; FIT2WT: Pass extra parameters to the USER function of CURVEFIT. ; ; MODIFICATION HISTORY: ; Written by: Han Wen, September 1996. (Adapted from FITWTIME) ;- pro NUM2EXP, n, A, Prob common FIT2WT, par2_ Mean1 = ABS(A(0)) Mean2 = ABS(A(1)) Offset1 = ABS(A(2)) Offset2 = ABS(A(3)) Norm = par2_.Norm Bin = par2_.Bin Verbose = par2_.Verbose if KEYWORD_SET(Verbose) then print,'params:',STRCOMPRESS(A) & wait,1 Prob = Norm*PROB2EXP(n, Mean1, Mean2, Offset1, Offset2, BIN=Bin, /SUM) end pro FIT2WT, Wts, Minima, MAX=Binmax, BINSIZE=Binsize, RESTRICT=Restrict, $ VERBOSE=Verbose common FIT2WT, par2_ NP = N_PARAMS() if (NP ne 2) then message, $ 'Must be called with 2 parameters: Dts, Minima' nwt = N_ELEMENTS(wts) if (N_ELEMENTS(Binsize) eq 0) then Binsize = 1 if (N_ELEMENTS(Binmax) eq 0) then Binmax = 1000L*Binsize if (N_ELEMENTS(Verbose) eq 0) then Verbose = 0 ; Histogram time intervals print,'Forming Time Interval Distribution...' & wait ,0.1 noccur = HIST1D(wts,MIN=0, MAX=Binmax, BINSIZE=Binsize, $ BINEDGE=-1, OBIN=wtbin) nbin = N_ELEMENTS(wtbin) i0 = -1 repeat i0=i0+1 until (noccur(i0) ne 0) tau1 = wtbin(i0) ; USER selects data region to fit print,'Applying USER-selected cuts...' ; Mark limit of the fit if KEYWORD_SET(Restrict) then begin x = wtbin(i0:nbin-1) plot,x,noccur(i0:nbin-1)>1,psym=10,/ytype oploterr,x,noccur(i0:nbin-1),SQRT(noccur(i0:nbin-1)),psym=3,/NOHAT wtmax= ROUND(MRK_LINE(/VERTICAL,/HELP)) wtmax= XQUERY('Selected cutoff:',wtmax) if (wtmax eq -1) then return nbin = LONG(wtmax/Binsize) + 1 endif noccur = noccur(0:nbin-1) ; Estimate initial parameters print,'Estimating Initial parameters...' nsig = 10 fit = FIT2EXP( wtbin(i0:nbin-1), noccur(i0:nbin-1), c2exp, sigc2, tau2 ) if (N_ELEMENTS(c2exp) ne 6) then return ; y = aL * exp( -bL*(x-cL) ), x < x0 ; = aR * exp( -bR*(x-cR) ), x > x0 aL = c2exp(0) & bL = c2exp(1) & cL = c2exp(2) aR = c2exp(3) & bR = c2exp(4) & cR = c2exp(5) saL = sigc2(0) & sbL = sigc2(1) & scL = sigc2(2) saR = sigc2(3) & sbR = sigc2(4) & scR = sigc2(5) ; Extract the "true" source1 and source2 parameters: ; ; p1 = count rate [cts/s] of source ; p2 = count rate [cts/s] of vetoed background ; pL = p1 + p2 = count rate [cts/s] for waiting times < Offset2 ; pR ~ p1/(1 + b*p2)= count rate [cts/s] for waiting times > Offset2 b = 1.4*(tau2 - tau1) ; Don't ask me why, but this seems to work nL = aL & pL = bL & TL = cL ; fairly well, phenomenologically. nR = aR & pR = bR & TR = cR p1 = pR*(1 + b*pL)/(1 + b*pR) & mean1 = 1/p1 p2 = (pL - pR)/(1 + b*pR) & mean2 = 1/p2 A = [Mean1, Mean2, tau1, tau2] & A0=A x = wtbin(i0:nbin-1) y = noccur(i0:nbin-1) w = 1./y par2_= { BIN:Binsize, NORM:nwt, VERBOSE:Verbose } yfit = CURVEFIT(x,y,w,A,sigA,FUNCTION_NAME='NUM2EXP', $ /NODERIVATIVE, CHI2=Chi2, ITER=Iter ) Prob = PROB2EXP(wtbin(0:nbin-1), A(0), A(1), A(2), A(3), $ BIN=Binsize, /SUM) if (A(2) lt wtbin(i0)) then message, /INFORMATIONAL, $ 'WARNING: Dead time from Source 1 < Smallest Waiting Time '+ $ 'with Counts.' param = { rate1 :1/A(0), $ rate2 :1/A(1), $ offset1 :A(2), $ offset2 :A(3) } sigma = { rate1 :sigA(0)/A(0)^2, $ rate2 :sigA(1)/A(1)^2, $ offset1 :sigA(2), $ offset2 :sigA(3) } dof = N_ELEMENTS(y) - N_ELEMENTS(A) CL = 1. - CHI_SQR1( Chi2*dof, dof ) minima = { param :param, $ sigma :sigma, $ wtbin :wtbin(0:nbin-1), $ noccur :noccur, $ prob :prob, $ nwt :nwt, $ chi2 :chi2*dof, $ dof :dof, $ CL :CL } end