SUBROUTINE pileup (ear, ne, param, ifl, photar, photer) INTEGER ne, ifl REAL ear(0:ne), param(6), photar(ne), photer(ne) c John Davis' forward-folding algorithm for dealing with pile-up. c This assumes that the input photar array has been multiplied by c the effective areas. The model parameters are the frame time (sec) c and the number of pile-ups to consider. c Model parameters are : c 1 frame time c 2 max number of piled up photons c 3 g0 - good event fraction c 4 alpha - grade migration factor c 5 psffrac - fraction of the counts to pile-up c 6 nregions - number of independent regions c Arguments : c ear r i: Energy ranges c ne i i: Number of elements in photar array c param r i: Model parameters c ifl i i: Data set c photar r i/r: Model flux c photer r i/r: Model flux errors c Calculates the expression c c exp(-tau/g0 Int dE' S'(E') ) [(exp(tau F)-1)/(tau F)] o S'(E) c c where tau is the frame time, S' is the input, and F is the functional c defined by c c F o f(E) = alpha Int_0^E dE' S'(E') f(E-E') c c The expression [(exp(tau F)-1)/(tau F)] is expanded in an infinite c series c c Sum_p (tau F)^(p-1)/p! INCLUDE '../../inc/xspec.inc' INTEGER ioutar, itmpar, iconv, ipfrac INTEGER status CHARACTER contxt*72 status = 0 c Grab the memory we need CALL udmget(ne, 7, ioutar, status) contxt = 'Failed to get memory for outar' IF ( status .NE. 0 ) GOTO 999 CALL udmget(ne, 7, itmpar, status) contxt = 'Failed to get memory for tmpar' IF ( status .NE. 0 ) GOTO 999 CALL udmget(ne, 7, iconv, status) contxt = 'Failed to get memory for conv' IF ( status .NE. 0 ) GOTO 999 CALL udmget(NINT(param(2))+1, 7, ipfrac, status) contxt = 'Failed to get memory for pfrac' IF ( status .NE. 0 ) GOTO 999 c Run the pile-up calculation CALL dopile(param(1), NINT(param(2)), param(3), param(4), & param(5), NINT(param(6)), ne, ear, ifl, photar, & MEMD(ioutar), MEMD(itmpar), MEMD(iconv), & MEMD(ipfrac)) c Free the temporary memory CALL udmfre(ioutar, 7, status) contxt = 'Failed to free memory for outar' IF ( status .NE. 0 ) GOTO 999 CALL udmfre(itmpar, 7, status) contxt = 'Failed to free memory for tmpar' IF ( status .NE. 0 ) GOTO 999 CALL udmfre(iconv, 7, status) contxt = 'Failed to free memory for conv' IF ( status .NE. 0 ) GOTO 999 CALL udmfre(ipfrac, 7, status) contxt = 'Failed to free memory for pfrac' IF ( status .NE. 0 ) GOTO 999 999 CONTINUE IF ( status .NE. 0 ) THEN CALL xwrite(contxt, 10) ENDIF RETURN END c ********************************************************************** SUBROUTINE dopile(frame, npiled, g0, alpha, psffrac, nregions, & ne, ear, ifl, photar, outar, tmpar, conv, pfrac) INTEGER ne, ifl, npiled, nregions REAL ear(0:ne), photar(ne), frame, g0, alpha, psffrac DOUBLE PRECISION outar(ne), tmpar(ne), conv(ne), pfrac(0:npiled) c Subroutine to do the actual pile-up calculation. Split out from the c main routine to make dynamic memory easier to handle. c c Arguments : c frame R i: Frame time (in seconds) c npiled I i: Max. number of photons to pile up c g0 R i: Grade 0 fraction c alpha R i: Index for parametrization of grade morphing c psffrac R i: Fraction of the flux to pile-up c nregions I i: Number of independent regions c ne I i: Number of energy bins c ifl I i: Dataset number (not used) c photar R i/r: Predicted counts/sec in each energy range c outar D w: Workspace array c tmpar D w: Workspace array c conv D w: Workspace array c pfrac D w: Workspace array DOUBLE PRECISION factor, pnorm, piledflux, sum, incounts DOUBLE PRECISION expfac, fract INTEGER ie, je, ipile, ioff CHARACTER wrtstr*255 c This is the fraction of the input counts that are used to calculate the c pile-up fract = psffrac/nregions c This is an offset we will require in the convolution in case ear(0) is c non-zero. Note that we assume the energy array is linear and starts at c a multiple of the energy bin size. ioff = -ear(0)/(ear(1)-ear(0)) c Set the output array and the temporary arrays for the single photon case c and find the sum of the input array pnorm = 0.d0 DO ie = 1, ne outar(ie) = DBLE(photar(ie)*fract) tmpar(ie) = outar(ie) pnorm = pnorm + outar(ie) ENDDO incounts = pnorm * frame WRITE(wrtstr,*) 'Pileup: input model count rate/frame/region = ', & incounts CALL xwrite(wrtstr, 15) c Calculate the exponential factor. expfac = EXP(-1 * frame * pnorm / g0) pfrac(0) = 0.0d0 pfrac(1) = pnorm * frame * expfac sum = pfrac(1) c Loop round pile-ups factor = 1.0d0 DO ipile = 2, npiled c Do the convolution (should probably do an fft here but use BFI for now). DO ie = 1, ne conv(ie) = 0. DO je = 1, ie+ioff conv(ie) = conv(ie) + DBLE(photar(je)*fract) & *tmpar(ie-je+1+ioff) ENDDO ENDDO c Add the latest iteration of the temporary array into the summation factor = factor * DBLE(alpha * frame/float(ipile)) pfrac(ipile) = 0.d0 DO ie = 1, ne tmpar(ie) = conv(ie) piledflux = tmpar(ie)*factor outar(ie) = outar(ie) + piledflux pfrac(ipile) = pfrac(ipile) + piledflux ENDDO pfrac(ipile) = pfrac(ipile) * frame * expfac sum = sum + pfrac(ipile) ENDDO c Now multiply by the exponential term and add in the unpiled fraction DO ie = 1, ne piledflux = piledflux + outar(ie)*expfac photar(ie) = nregions*SNGL(outar(ie)*expfac) & + (1.-psffrac)*photar(ie) ENDDO WRITE(wrtstr,*) 'Pileup: output model count rate/frame/region = ', & piledflux * frame CALL xwrite(wrtstr, 15) c Write out diagnostic information if required. The second number is the c probability of seeing that number of photons piled up in a frame and the c third number is the probability of an observed event being due to that c number of piled up photons. IF ( sum .GT. 0.0 ) THEN WRITE(wrtstr,'(i4,3(1x,1pg12.6))') 0, expfac, pfrac(0), & pfrac(0)/sum CALL xwrite(wrtstr, 25) pnorm = expfac DO ipile = 1, npiled pnorm = pnorm * incounts / ipile WRITE(wrtstr,'(i4,3(1x,1pg12.6))') ipile, pnorm, & pfrac(ipile), pfrac(ipile)/sum CALL xwrite(wrtstr, 25) ENDDO ENDIF RETURN END