Analytical modeling of pulse-pileup distortion using the true pulse shape; applications to Fermi-GBM
Vandiver Chaplin, Narayana Bhat, Michael Briggs, Valerie Connaughton
AAnalytical modeling of pulse-pileup distortion using the true pulse shape;applications to
Fermi-GBM (cid:73)
Vandiver Chaplin a, ∗ , Narayana Bhat a , Michael S. Briggs a,b , Valerie Connaughton a,b a CSPAR, University of Alabama in Huntsville. 301 Sparkman Drive, Huntsville, AL 35899 USA b Department of Physics, University of Alabama in Huntsville. 301 Sparkman Drive, Huntsville, AL 35899 USA
Abstract
Pulse-pileup affects most photon counting systems and occurs when photon detections occur faster than thedetector’s shaping and recovery time. At high input rates, shaped pulses interfere and the source spectrum,as well as intensity information, get distorted. For instruments using bipolar pulse shaping there are twoaspects to consider: ‘peak’ and ‘tail’ pileup effects, which raise and lower the measured energy, respectively.Peak effects have been extensively modeled in the past. Tail effects have garnered less attention due toincreased complexity. We leverage previous work to derive an accurate, semi-analytical prediction for peakand tail pileup including high order effects. We use the pulse shape of the detectors of the
Fermi Gamma-ray Burst Monitor . The measured spectrum is calculated by expressing exposure time with a state-spaceexpansion of overlapping pileup states and is valid up to very high rates. The model correctly predictsdeadtime and pileup losses, and energy-dependent losses due to tail subtraction (sub-threshold) effects. Wediscuss total losses in terms of the true rate of photon detections versus the recorded count rate.
Keywords: pulse pileup, deadtime, GBM, TGF
1. Introduction
Pulse pileup affects X- and gamma-ray countingsystems in the presence of high-intensity sources.At high rates the detected events become increas-ingly different from those registered by the instru-ment. Detected events produce shaped electronicpulses that require some minimum time interval tomeasure the pulse height, and additional time to re-cover. The total time elapsed is the ‘pulse window’.Additional detections in this interval will modifythe expected pulse shape and distort the inferredpulse-height information and detection rate. In thispaper we elaborate on previous methods to developa correction for pileup using the true bipolar pulseshape, while also requiring fewer assumptions [1]. (cid:73) http://dx.doi.org/10.1016/j.nima.2013.03.067 ∗ Corresponding author. Electronic mail: [email protected]. Phone: 1-256-961-7514 Terms such as detected, input will refer to the photonsinteracting in the detector volume. Terms such as regis-tered, recorded, measured , observed will refer to the subsetof photons and pileup events which are finally counted bythe instrument We apply this method to the space-borne gamma-ray instrument
Fermi Gamma-ray Burst Monitor (GBM). GBM consists of a set of inorganic scintil-lation detectors and is one of two instruments on the
Fermi gamma-ray space telescope. An overview ofthe GBM detector hardware, science mission, andthe context for pulse-pileup studies is described insection 2.In bipolar-pulsed instruments pileup has two ef-fects on the measured pulse spectrum, depending onthe arrangement of pileup events.
Peak effects arisefrom the addition of positive pulse sections, causinga single apparent high-energy count.
Tail effects oc-cur when events are detected in the negative tail of aprevious event, causing a second count with reducedapparent energy. For pulse shapes with large neg-ative swings, modeling the tail pileup is especiallyimportant since losses (non-registered events) willoccur if the summed peak is below threshold. Suchlosses depend on the input spectrum, a dependencewhich gets stronger as the bipolar peak-to-peak ra-tio approaches unity.With some assumptions the distorted spectrumcan be predicted in terms of the true input rate and
Preprint submitted to Nuclear Instruments and Methods A November 11, 2018 a r X i v : . [ a s t r o - ph . I M ] A p r pectrum. Cano-Ott [2] and Danon [3] model peakeffects in the context of unipolar pulses. Taguchi[1] model both the peak and tail effects in an X-ray counter with a bipolar pulse. The approachgenerally taken is to model the “pileup response”of the detector, which is a conditional probability(a.k.a. “likelihood function”) of recording a certainenergy ε resulting from pileup of two events, withenergies E and E :Pr( ε | E , E )Calculating such functions requires basic informa-tion about the hardware response of the detector tophoton energy deposition, and a statistical modelfor the time-separation between detections.Assuming linearity in photo-current productionand pulse shaping, the detector’s observed analogsignal is prescribed using a superposition of individ-ual pulses. Then the likelihood functions are cal-culated assuming Poisson distributed events. Thepileup spectrum is written using “total probabil-ity”, an integral over measured energies (or sumover channels) of the pileup response times theprobability of input energies (i.e., the input spec-trum). Sections 6 and 7 derive the likelihood func-tions we use for modeling GBM pileup.A common simplification is to forego the truepulse shape when calculating the pileup-energy like-lihood, and instead use a mathematical approxima-tion. This has the advantage that likelihoods canhave closed form expressions. For example, [1] usea triangular or ‘delta’ approximation, such that thesummed peak-time (and amplitude) is given by sim-ple geometry, and its dependence on the event sep-aration is algebraically invertible. In section 6.2.3we demonstrate why invertibility is necessary for aclosed form likelihood expression. The disadvan-tage in approximating is that the true pulse shapegives significantly more accurate predictions thanthe approximate shapes. [2] demonstrate this byevaluating their model with an accurate Gaussianpulse and then with several approximate shapes.They compare results with Monte Carlo simulationand show that only the true pulse shape is accu-rate over the entire energy range. They includeonly first order peak effects and reach a maximumpileup fraction well below what is expected in ourapplications.The technique of [1] allows for modeling high or-der pileup effects using an iterative approximation.Tail modeling uses a convolution method rather than the pileup-likelihood technique used for thepeak effect. This admits a number of simplifyingassumptions about the timing and energy distribu-tion of tail events and avoids complications in cal-culating additional likelihood functions. However,when applied to the true pulse shape and specificinput spectral shapes in our application, we foundsuch assumptions did not yield sufficiently accurateresults. This is perhaps due to a larger negativeamplitude in our instrument, making the distortedspectrum more affected by tail pileup.The peak pileup method described in this paperis largely an application of the method of [1] to theGBM pulse shape. However, we present a differentmethod for tail modeling which extends the like-lihood function treatment to tail regions, allowingus to model the spectral and timing randomness oftail events. This method also accounts for eventsoccurring in a fixed deadtime that is separate fromthe peak pileup interval.A novelty of this approach is that the total spec-trum is obtained by partitioning time into overlap-ping pulse windows. Pulse windows are written asstates of a Poisson process, each one having an asso-ciated peak and tail spectrum. In section 4 a stateis defined by the number of events in the window,with pileup states having a non-zero number (the pileup order ). In sections 6 and 7, the recordedenergy distributions are derived based on the num-ber of events per state and how they are arrangedinto three sub-intervals of the pulse window. Thetotal time in which there is a non-zero analog sig-nal is then the union of all pulse windows, and thetotal spectrum is the corresponding union of peakand tail spectra per state. This is expressed as su-perposition of pulse states minus overlapping terms(section 9). We compare the model to Monte Carlosimulations in section 10.3.As a final motivation for using the true bipolarpulse shape, we note there is only a small increase incomputational overhead. Despite the apparent con-venience of closed form likelihood functions, the ac-tual difference in utility between them and a numer-ical function is small. Most spectra are sufficientlycomplicated that the total probability of likelihoodtimes spectrum over energy is not analytically inte-grable, thus numerical evaluation is required. Thisis especially the case for non-ideal detectors, whichhave a complicated instrument response R ( E, E γ ),and the deposited energy spectra are already ex-pressed numerically. Other than some storage andmemory access overhead, computation of the pileup2pectra requires about the same number of opera-tions in either approach. As will be shown, for theconstant Poisson process the likelihood functionsare independent of source intensity and thus canbe calculated a priori , stored in memory, and readby a separate program as required. In section 10we describe a simple numerical calculation of thelikelihood functions.
2. Fermi-GBM
GBM is a gamma-ray counting instrumentaboard the
Fermi Gamma-ray Space Telescope . Itwas launched in June 2008 and has been in near-continuous operation since the start of normal op-erations one month later. GBM consists of 12sodium-iodide (NaI) and two bismuth-germanate(BGO) detectors producing time and energy re-solved data sets. NaI detectors have an effectiveenergy range from approximately 8 keV to 1 MeV,and BGOs from 200 keV to 40 MeV. Pulse heightsare digitized into 128 pseudo-logarithmic spectralchannels per detector. Detector gains are calibratedsuch that channel boundaries lie between 0 to 5volts. More information on GBM hardware, detec-tors, and electronics is available in [4]. Instrumentcalibration (the relationship between peak voltageand input energy) studies are described in [5]. Each detector has its own shaping and digitiza-tion firmware for performing pulse height analysis(PHA) of detected events. Pulse-pileup, when it oc-curs, is on a detector-by-detector basis rather thanan integrated signal. The pulse shape has a finitewidth that requires 0 . µs to register a single eventand an additional 3 . − µs for baseline recovery(section 3). Generally speaking, pulse pileup oc-curs when the separation between detected eventsis smaller than the registration + recovery time.The negative peak amplitude is about 70% of thepositive peak amplitude. To help mitigate pileupeffects, a fixed deadtime of τ = 2 . µs is appliedonce a peak has been found, preventing pulse mea-surement during most of the recovery. Pulse pileuphowever extends the time for baseline recovery and For simplicity the figures in this paper are calculatedusing BGO channel energies that are approximated to belogarithmic: E m = (150keV)(1 . m gives the start ofthe m th energy channel in keV. Plots can be approximately“converted” to the NaI energy range for a channel m , using E m = (4keV) ∗ (1 . m . Actual channel energies are notprecisely logarithmic. τ becomes insufficient, adding tail distortions to themeasured spectrum.GBM is primarily an astrophysics experiment de-signed to detect transient cosmological sources suchas gamma-ray bursts (GRBs). Detectors are uncol-limated and typical source + background rates are1 − × counts per second (cps) per detector.In the majority of cases a simple non-paralyzabledeadtime correction is sufficient, but pileup is po-tentially problematic for several types of sources ob-served. Depending on the application, pileup effectswould become non-negligible at an input rate ( λ )around 10 cps in one detector, corresponding toan observed counting rate ( λ Rec ) around 60,000 -70,000 cps in that detector.In NaI detectors soft-gamma repeaters (SGRs)can have peak recorded rates near this maximum[6]. Solar flares are also routinely observed at theselevels. Terrestrial gamma-ray flashes (TGFs) typ-ically peak around λ Rec ∼ cps in BGO detec-tors [7]. In such cases the true intensity is muchgreater than what is inferred from a simple dead-time correction. At these levels a substantial frac-tion ( ∼
65% - 75%) of input photons pile-up (sec-tion 11), and the subset of resulting counts is highlydistorted. Modeling the spectrum and energetics ofsuch high-intensity transients thus requires under-standing pileup distortion in GBM, and our goal isto have a fast, flexible numerical prediction whichmight be employed in spectral analysis.The model we present gives accurate results upto very high values of the input rate λ . Demon-strations of spectral distortion versus rate are givenin section 10. The relation between λ Rec and λ isconsidered in more detail in section 11.
3. Pulse shape
GBM has well-studied electronic and detectorproperties that allow us to model its behavior withmathematical expressions. Shaping circuits makepulses whose first zero-crossing is fixed relative tothe pulse start, regardless of the total photocurrent.We assume the zeroth order shape can be writtenin a separable form where the shape f ( t ) is mod-ulated by a scalar that depends on the amount ofinput energy. Thus a given pulse can be written V ( t ) = v f ( t ), where v is an energy dependentcoefficient, and f ( t ) is the basic pulse shape. Us-ing data collected in pre-launch testing, we fit thefunction below to sampled analog signal forms, for3 range of input energies: f ( t ) = K ( c ∗ t α − c ∗ t β ) e − γt (1)where c , c , α, β, γ are fitted constants that param-eterize the pulse shape. K is a normalizing constantsuch that f ( t ) is unity at its peak (section 6.1).The best-fit shape parameters are c = 26 . , c =31 . , α = 1 . , β = 3 . , γ = 2 . , K = 0 . t in microseconds. The figures and specificresults in this paper refer to equation (1) with theseparameters. However, as evidence for the general-ity of the technique, we have tried a wide range ofparameter values and find the model, when com-pared to Monte Carlo simulations using the samepulse shape, to be equally accurate in all cases. Τ Τ A Τ B Τ C t p - - Μ s f H t L Figure 1: A single unit pulse showing the zeroth-order peaktime t p , the applied deadtime τ , and the partitions A, B,and C. The pulse “window” is ∆ τ = τ A + τ B + τ C = 4.5 µs . The pulse window initiated by detection of an ini-tial event will be partitioned into three regions, eachcorresponding to different distortion effects. We la-bel them ‘A’, ‘B’, and ‘C’, having widths τ A , τ B ,and τ C , respectively (figure 1). A and C will alsobe called the peak and tail regions. In our model,the particular distortion observed will be parame-terized in terms of the number of events in eachsub-interval, and pileup defined as the presence ofadditional photons in any of these sub-intervals.Choice of widths depends on the specifics of the pulse shape, the algorithm which measures pulseheight, and the deadtime implementation. How-ever, we can give the following generic definitionsfor each region, based on the case of first-orderpileup. Figure 2 depicts the three cases of first orderpileup, described below.Suppose two photons are detected within a pulse-width of each other. The initial pulse we call the‘zeroth’ event, and let it begin at time t = 0. Thenthe next event, beginning at t and having an as-sociated voltage v , can occur in one of the threeintervals A, B, or C. If t is in A, peaks will addand a ‘summed’ pulse height appears which is afunction of v , v , and t − t . If t is in B, v ismeasured correctly, and the second peak is shifteddown due to the negative baseline. In the case ofGBM this coincides with deadtime, so there is nosecond count. Finally if t is in C, there is a secondmeasurement, but the peak is shifted down from itsnominal value v . Thus the widths are defined asfollows. ‘ τ A ’ is the minimum time required to re-solve the pulse height of one event. For our model toapply it should be (approximately) ≤ the positivesignal width. ‘ τ B ’ is a dead interval: if the secondevent occurs in this region it is not measured, butwill influence the region C. ‘ τ C ’ is a live region, butlasts until the zeroth pulse is baseline-recovered.In the specific context of GBM, the width τ A ofthe first region is the time-to-max t p plus a bufferfollowing the peak: τ A = t p + τ Buff = 0 . µs This buffer is part of the instrument’s digital peakfinding algorithm, which requires a decreasing volt-age for a sequence of digital samples before regis-tering a pulse height. If additional events occurbefore the buffer expires (i.e., anywhere in regionA), they will cause peak-pileup. [8]Events in part B are lost. Applied deadtime τ contains the B interval. τ B smaller than τ becausea pulse occurring on its periphery might carry intothe live time and be measured. Likewise, a pulse oc-curring in the post-peak buffer would be measured.Thus τ B is the deadtime τ minus this buffer on the This scheme is designed to eliminate measurements dueto electrical noise, and to apply an energy-independent deadtime for each event. GBM samples analog pulses with aperiod of 0.104 µs . The buffer is programmed to be foursamples long, so τ Buff = 0 . µs . The deadtime results fromwaiting an additional 21 samples after a peak is found andbuffered, so τ = (4 + 21) ∗ .
104 = 2 . µs . τ B = τ − ( t p + τ Buff ) = τ − τ A = 1 . µs where τ = 2 . µs is the GBM deadtime.Region C is cutoff at an arbitrary point wherethe baseline (due to the initial pulse) can be con-sidered to be recovered. But the instrument is alsolive here, and peak pileup from additional eventsis possible. Thus we also require it to be shortenough that only a single peak can be recorded inthis region. For GBM modeling we use the value τ C = 1 . µs , to give a total pulse window of width τ A + τ B + τ C = ∆ τ = 4 . µs (2)Predictions from the model have some depen-dence on the value of τ C . We examined this bytrying both larger (up to 6 µs ) and smaller (downto 4 µs ) pulse window sizes, varying τ C only. Themodel’s accuracy was not appreciably different forwindow sizes between 4 . − . µs . For smaller win-dow size, the amount of tail distortion becomesunder-predicted when compared to MC simulations(i.e., not enough low-energy counts are predicted).For larger sizes, it becomes over-predicted. Thisis not surprising, since the probability of tail dis-tortion increases with increasing τ C (see equation(49)). For a successful model τ C must be a goodapproximation to the instrument’s actual recoverytime, such that the computed probability of taildistortions is equal (within statistics) to the frac-tion of tail counts from simulation. The value usedin this paper was not optimized for model accuracybut chosen a priori based on the known pulse shapeand the constraints given in the previous paragraph.
4. Pulse-pileup
Having partitioned the pulse shape into three ad-jacent regions, we can now describe pulse-pileup interms of the number of events in each one. Assum-ing a zeroth event to begin the window, one or moreadded events in region A constitutes peak pileup .One or more events occurring in region C cause tailpileup , or the tail effect. Additionally, pulses in Bcontribute to baseline reductions for the C events,and this is a more extreme category of tail effects.Together these produce artificial raising and lower-ing distortions in the measured spectrum, which itis our goal to model.Assuming a zeroth event with the A-B-C parti-tions, the state or configuration of pulse-pileup can be represented using three non-negative integers a, b, c which give the number of additional eventsin each interval (excluding event 0). We define thesymbol (cid:104) a, b, c (cid:105) to describe the pileup state in thepulse window ∆ τ . The total order of pileup is thenumber of additional events, such that order zerois a single count without pileup, first order is oneadded event in A, B, or C, etc. Thus the genericconfiguration (cid:104) a, b, c (cid:105) represents the case of an inci-dent event (the zeroth event) followed by ( a + b + c )events within the window. In later sections we ad-dress the probability and independence of pileupstates, but for now we need only the definition:Event (cid:104) a, b, c (cid:105)≡ a pulses in A, b in B, and c in CIn general this represents the event of ( a + b + c ) th order pileup, where ( a + b + c ) ≥
0, and (cid:104) , , (cid:105) is the event of a single recorded count, accuratelymeasured. In this scheme the maximum numberof measured counts is two (one from the peak, onefrom the tail).Figure 2 demonstrates the simplest case: first-order pileup. Figure 3 gives several examplesof higher order pileup, with the pulse partitionsshown. In the following sections we calculate pulseheight distributions corresponding to each state (cid:104) a, b, c (cid:105) , and an input energy spectrum. The peakspectrum will depend only on a , but the tail spec-trum depends on all three orders.One consequence of pulse pileup is that the dead-time is extended beyond its the nominal value τ .For example, pileup in region A results in a summedpeak at a time ≥ the single-pulse peak time t p . Andsince the deadtime is imposed from this registeredpeak, it effectively migrates forward, contributingto paralyzable deadtime [9]. Strictly speaking thepulse partitions should have a corresponding shift.Instead we use a semi-paralyzable model where thesummed peak time (section 10.1) can wander intoregion B, but the partitions and deadtime remainfixed. Figures 3(a-c) show how the registered peakshifts forward while ABC partitions are imposedfrom the zeroth event. The result is shown to beaccurate when compared to Monte Carlo simula-tions.
5. Interval Distribution
Within fixed intervals the interval distribution(i.e., the probability density (PDF) of separation5
Τ Τ Τ s @ Μ s D - - V H t L H a L Peak effect H st order L s @ Μ s D H b L Peak + deadtime H no distortion L s @ Μ s D H c L Tail effect H st order L Figure 2: The three cases of ‘first order’ pileup, (cid:104) (cid:105) , (cid:104) (cid:105) , and (cid:104) (cid:105) , showing the measured peak for two events of equalenergy, and the dead time τ as imposed by GBM hardware. (a) shows the peak effect, and (c) the tail effect. Panel (b) depictsa nominal case where one count is accurately measured and the next is lost. Typically this is not regarded as ‘pulse pileup’ asthere is no associated spectral distortion of the pulse height, only the count rate, which can be corrected by simpler means. Τ Τ A Τ B Τ C Τ Τ Τ A Τ B Τ C Τ Τ Τ A Τ B Τ C s @ Μ s D - V H t L H a L X \ s @ Μ s D H b L X \ s @ Μ s D H c L X \ Figure 3: Higher order pileup examples, with the A-B-C partitions shown. (a) second order peak pileup. (b) third orderpileup, with peak and and tail effects. (c) a third order case of the deadtime+tail effect. Recorded pulse height in C dependson the tail from pulses in A and B. time between events) can be calculated under theassumption of a fixed-rate Poisson process. Thisderivation is well-known, appearing in [1] amongother places, but because of its centrality in thecalculation of pileup likelihood we also present it.For the Poisson process of rate λ , the PDF of theseparation s between any event and the next is f s ( s ) = λe − λs (3)which is defined on s ∈ (0 , ∞ ). Define the separa-tion of the i th event from the ( i − th as s i ≡ t i − t i − (4)where t i is the time of event i . Now suppose a finiteinterval ranging from 0 to ∆ t , with the zeroth event outside the interval at t = 0. Then the probabilityof exactly one event in the interval is given by the(joint) probability that s ≤ ∆ t and s + s > ∆ t ,and the second condition is equivalent to s > ∆ t − s . Pr( s ≤ ∆ t and ( s > ∆ t − s ) )= s (cid:90) λe − λs ds ∞ (cid:90) ∆ t − s λe − λs ds (5)= λs e − λ ∆ t Differentiating in s gives the joint probability den-sity f s ( s , f s ( s , n ) forthe random variables s and n (number of events6 igure 4: 2 nd order pileup in the offset interval [ t, t + ∆ t ],showing events t , t , t , t . If we assume finite pulses occurat each time, events t , t are piled-up, with possible taileffects from t . s = s (cid:48) + t − t , and the random variables s (cid:48) , s are identically distributed. in ∆ t ) f s ( s ,
1) = dds Pr = λe − λ ∆ t (6)The conditional probability density, f s | n ( s ), is givenby dividing the joint density f s ( s ,
1) by the prob-ability that n = 1, a procedure equivalent to nor-malizing over the interval [0 , ∆ t ]: f s | n =1 ( s ) = f s ( s , (cid:82) ∆ t f s ( s (cid:48) , ds (cid:48) = λe − λ ∆ t λ ∆ te − λ ∆ t = 1∆ t (7)which is a constant, and is independent of λ . [10]For exactly two events in the interval, the condi-tion is s + s ≤ ∆ t and s + s + s > ∆ t (e.g.figure 4). This can also be written s ≤ ∆ t − s and s > ∆ t − ( s + s ). ThenPr( s ≤ ∆ t for 2 in ∆ t )= s (cid:90) λe − λs ds × ∆ t − s (cid:90) λe − λs ds ∞ (cid:90) ∆ t − ( s + s ) λe − λs ds (8)= λ e − λ ∆ t ( s ∆ t − s n = 2 in ∆ t givesthe second order conditional density: (9) f s | ( s ) = f s ( s , (cid:82) ∆ t f s ( s (cid:48) , ds (cid:48) = 2∆ t (∆ t − s )Because in this process event separations are inde-pendent, identically distributed random variables, s , s follow the same probability density f s | . Forthe general case of n events, the joint density func-tion f s ( s , n ) is calculated as in equation (8) and the conditional density f s | n ( s ) just like (9), whichcan be shown by induction to be, for all n ≥ f s | n ( s ) = n ∆ t n (∆ t − s ) n − (10) The above density was derived assuming t =0, but equation (10) turns out to have the sameform regardless of the position of the previous event t before the interval. Consider the second-orderintegral in equation (8), but now let the intervalbegin at some arbitrary time t ≥ t to t + ∆ t . Let the zeroth event occur at an arbitrarytime t ≤ t . Now the 2 nd order pileup conditions,for t , t in the interval, are: t − t < s < ( t + ∆ t ) − t s + s ≤ ( t + ∆ t ) − t s + s + s > ( t + ∆ t ) − t An example is shown in figure (4). Using the sym-bols T A ≡ t − t and T B ≡ ( t +∆ t ) − t for shorthand,the second-order joint probability isPr( s ,
2) = s (cid:90) T A λe − λs ds × T B − s (cid:90) λe − λs ds ∞ (cid:90) T B − ( s + s ) λe − λs ds (11)= λ e − λT B (cid:104) T A − s + 2 T B ( s − T A ) (cid:105) Differentiating and normalizing gives the condi-tional density: f s | ( s ) = 2( T B − s )( T B − T A ) = 2 ([ t + ∆ t − t ] − s )(∆ t ) (12)= 2(∆ t − [ s − ( t − t )])(∆ t ) Taking just the portion of s that is in the normal-ization interval, s (cid:48) = s − ( t − t ), and substitutingthis into the the above, f s | ( s (cid:48) ) = 2(∆ t − s (cid:48) )∆ t (13)where 0 < s (cid:48) ≤ ∆ t , which is exactly the resultwhen t is at the start of the interval, equation (9).The subtracted term ( t − t ) is a constant, so s and7 (cid:48) have the same distribution except for this offset.This generalizes to the result f s | n ( s ) = n ∆ t n (∆ t − s ) n − (14)for n events in an arbitrary interval [ t, t + ∆ t ],and 0 ≤ s ≤ ∆ t . Using this result we can writethe probability densities assuming a definite pileupstate (cid:104) a, b, c (cid:105) : f s | a ( s ) = a ( τ A ) a ( τ A − s ) a − (15) f s | b ( s ) = b ( τ B ) b ( τ B − s ) b − (16) f s | c ( s ) = c ( τ C ) c ( τ C − s ) c − (17)
6. Peak pileup effect
The peak modeling technique is an applicationof the one presented in [1]. The method is to firstderive expressions of the form Pr( ε | E (cid:48) , E (cid:48)(cid:48) ), whichgive the likelihood of recording energy ε when twoevents are detected with energies E (cid:48) , E (cid:48)(cid:48) , within atime τ A of each other. For first-order the modelis accurate without any approximation. At higherorders an iterative approximation is used since thedimensionality (i.e., number of random variables)becomes large. Then the total probability of ε isgiven in terms of the pileup likelihood and the inputspectrum. Without pileup, a single pulse at time 0 has avoltage signal that can be expressed as V ( t ) = v f ( t ) (18)where f ( t ) is the normalized pulse shape, such thatat its maximum, f = 1. Therefore the peak voltageis V = v . The peak-time is given by the firstsolution of df ( t ) dt (cid:12)(cid:12)(cid:12) t p = 0 (19)and we will denote it by the symbol t p . Thus, f (cid:48) ( t p ) = 0 f ( t p ) = 1 V ( t p ) = v t p s - - Μ s V H t L Figure 5: A case of 1 st order peak pileup, showing the mod-eled peak at t p . First order peak-pileup is depicted in figure 2(a)and figure 5. The apparent pulse is a linear combi-nation of unit pulses f ( t ), V ( t ) = v f ( t − t ) + v f ( t − t ) (20)where t , t are the incident event times, and v , v are their peak amplitudes. For the constant rateprocess, the total time offset of the pair is irrele-vant, so we can let t = 0. Then t is equal to theseparation of the two events, s . The maxima ofequation (20) are given by ddt V ( t ) (cid:12)(cid:12)(cid:12) t p = (cid:34) v df ( t ) dt + v df ( t − s ) dt (cid:35) t p = 0 (21)We define the symbol t p as the position of themeasured first-order peak, such that it satisfies theabove equation. As discussed in section 4 the ob-served peak is shifted forward, so t p ≥ t p , and itdepends on the random variables s , v , v : t p = t p ( s , v , v ) (22)Having a model for t p ( s , v , v ) is necessary so thatthe sum (20) can be evaluated for the observed first-order energy. Depending on the pulse shape f ( t ),equation (21) might have a closed form solution for t p . In previous work simplifications are employedsuch that t p is given by simple geometrical argu-ments. But in the present case, substitution of the8ulse shape, equation (1), into equation (21) resultsin an expression that cannot be readily inverted. Asa result we use an empirical model for the peak-timethat is presented in section 10.1.Note that for a given set of values ( s , v , v )there may be two solutions to equation (21) (twomaxima in the peak interval, for example in fig-ure 5). Selecting the correct peak requires knowl-edge of the peak-finding algorithm of the detectorsystem in question. The buffering scheme of theGBM pulse-height analyzers generally causes thelast maximum to be the measured one. However if v is much larger than v , and s is in the buffer, t p → t p because the falling derivative has a muchlarger (negative) value than the rising derivative inthis region. When the two add, the net change isstill negative, so the decreasing sample criterion issatisfied. Such discrete logic is expressed as piece-wise behavior in the model for t p . Let us assume that we have a sufficient modelfor t p ( s , v , v ). Unlike f ( t p ) = 1, f ( t p ) (cid:54) = 1 ingeneral. The recorded pulse height is V ( t p ) = v f ( t p ) + v f ( t p − s ) (23)and the recorded first-order pileup energy ε is ε = ξ [ v f ( t p ) + v f ( t p − s )] (24)where ξ [ v ] is the channel or voltage-to-energy con-version, determined by standard methods of instru-ment calibration. For modeling input spectra, weassume the inverse ξ − exists and the coefficientsof input pulses can be calculated as v i = ξ − [ E i ]where E i is energy deposited by a detected gamma-ray. E i is of course converted from incident photonenergy by the various physical processes of the de-tector, but this is extraneous to the current prob-lem. It is sufficient to say that E i is the recordedenergy of the i th count in the absence of pileup ef-fects.The recorded 1st-order energy, ε , for a given in-stance of the variables ( s , E , E ), is: (25) ε = ξ (cid:20) ξ − [ E ] f ( t p ) + ξ − [ E ] f ( t p − s ) (cid:21) ≡ ε ( s , E , E )Note that if ξ [ v ] is approximately of the form y = mx , then this simplifies to E f ( t p ) + E f ( t p − s ).In section 10.1 we give the expressions for t p and ε using the true GBM pulse shape. Figure 6 iscalculated with the above formula, and t p from thatsection. Discontinuity at the end of region A occursbecause the two pulses are sufficiently separated forthe first to be measured correctly. It is expressedas piecewise behavior in t p , given in section 10.1. ε What we need is a model giving the observeddistribution of the recorded energy ε in the eventof peak pileup. This means we need the dis-tribution of ε over all possible realizations of { ( s , E , E ) } . Let us examine the distribution ε due to s only , holding E , E fixed at some arbi-trary values. This gives a conditional form of equa-tion (25), ε ( s | E , E ).Given the interval PDF f s | n =1 ( s ) of equation(15), the PDF of the dependent variable ε can befound in the standard way [10]. For s in the peakregion, ε ( s | E , E ) is monotonic, and we canwrite its probability density as f (cid:104) (cid:105) peak ( ε ( s ) ) = f s | n =1 ( s ) | ∂ε ( s ) ∂s | s = s (26)where ε ( s | E , E ) is written ε ( s ) for brevity.It is understood that t p and ε are functions of thethree random variables.We can write the probability of pulse pileup intoan energy interval ε ( s ) to ε ( s ) + dε as d [Pr (cid:104) (cid:105) peak ( ε ( s ) | E , E )] = f (cid:104) (cid:105) peak ( ε ) dε (27)Then we can write the following for the probabilitythat ε is in the discrete channel ε to ε + ∆ ε , in theevent that it is due to E , E pileup:Pr (cid:104) (cid:105) peak ( ε | E , E ) = ε +∆ ε (cid:90) ε f (cid:104) (cid:105) peak ( ε (cid:48) ) dε (cid:48) (28)= ε +∆ ε (cid:90) ε f s | n =1 ( s (cid:48) ) | ∂ε ( s ) ∂s | s = s (cid:48) dε (cid:48) (29)This can be converted to an integral over s usingthe positive-definite Jacobian determinant (proba-bility must be positive) (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( ε (cid:48) ) ∂ ( s ) (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) ∂ε ( s ) ∂s (cid:12)(cid:12)(cid:12)(cid:12) s = s (30)9 Ç Ç Ç Ç Ç Ç Ç ÇÅ Å Å Å Å Å Å Å Å Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç ÇÅ Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å ÅÇ Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç Ç ÇÅ Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å Å ÅÅ ÇÇ
A B C
Model for peak energyModel for tail energyMonte carlo measuremnt with phase uncertainty0 1 2 3 4 5 601000200030004000500060007000 s @ Μ s D ¶ @ k e V D ¶ H s L with E = E = Figure 6: Plot of 1st-order recorded energy ε ( s , E , E ) for E , E at 3500 keV, if the second event occurs at s . The curvesrepresent the model for measured energy using the function t p , which is different in the peak than in the tail. Shaded areasbetween points give the range of recorded energy due to the phasing of digital samples. Discontinuity in the peak curve (BLUE)occurs when the separation is large enough for the initial pulse height to be measured correctly. The divisions A, B, C referto the pulse partition of section 3.1. The curve in A refers to a single measurement of two events that have pulse pileup. Thecurves in B-C refer to measurement of a second event (tail effect), which is addressed in section 7. In the latter case, two measurements are made, with peak zero measured at the input value of E . giving Pr (cid:104) (cid:105) peak ( ε | E , E )= s ( ε +∆ ε ) (cid:90) s ( ε ) f s | ( s (cid:48) ) | ∂ε ( s ) ∂s | s (cid:48) (cid:12)(cid:12)(cid:12)(cid:12) ∂ε ( s ) ∂s (cid:12)(cid:12)(cid:12)(cid:12) s (cid:48) ds (cid:48) (31)= s ( ε +∆ ε ) (cid:90) s ( ε ) f s | ( s (cid:48) ) ds (cid:48) Examples of this likelihood function are shown infigure 7. Evaluation of the above integration limitsevidently requires knowledge of an inverse functiongiving the separation in terms of the summed peakand two energies. Therefore, if a closed form ex-pression for equation (31) is desired, the recordedenergy ε ( s | E , E ) must be an analytical expres-sion and have an inverse, s = ε − ( ε, E , E ).The simplification of [1] using triangular pulsesallows such an inversion. In the present case, thetrue pulse shape results in an expression for ε which is too complicated to invert. Thus we havethe pileup likelihood expressed in implicit form, andmust evaluate it numerically. In section 10 we de- scribe a simple algorithm for doing the implicit eval-uation.Regardless, we can write the first-order peakpileup contribution to the spectrum as (32) P (cid:104) (cid:105) ( ε ) = ∞ (cid:90) ∞ (cid:90) Pr (cid:104) (cid:105) peak ( ε | E , E ) × S ( E ) S ( E ) dE dE S ( E ) is the PDF of the detected energy spectrum,and is normalized to unity. S ( E ) is typically mod-eled in terms of an externally incident photon spec-trum, S γ ( E γ ), and the detector response function R ( E | E γ ), for example as described in [11]: S ( E ) = (cid:90) R ( E | E γ ) S γ ( E γ ) dE γ (33)where S γ is the PDF of the photon spectrum. Sincereal detector systems are sensitive over a finite en-ergy range, we note that the limits of integration inequation (32) can be finite values E min → E max aslong as S ( E ) is normalized over this interval. ThePDF of the spectral contribution is p (cid:104) (cid:105) ( ε ) = ∂∂ε P (cid:104) (cid:105) ( ε ) (34)10 (cid:61)
500 keVE (cid:61)
500 keVE (cid:61)
500 keVE (cid:61) (cid:61)
500 keVE (cid:61) (cid:61)
500 keVE (cid:61) (cid:61)
700 keVE (cid:61)
500 keVE (cid:61) (cid:61)
500 keVE (cid:61) (cid:61)
500 keVE (cid:61) (cid:61)
500 keV P r (cid:72) (cid:182) (cid:200) E , E (cid:76) (cid:72) a (cid:76) (cid:72) b (cid:76) (cid:72) c (cid:76) (cid:182) (cid:64) keV (cid:68) P r (cid:72) (cid:182) (cid:200) E , E (cid:76) (cid:72) d (cid:76) (cid:72) e (cid:76) (cid:72) f (cid:76) (cid:72) g (cid:76) (cid:182) (cid:64) keV (cid:68) (cid:72) h (cid:76) Figure 7: Pr ( ε | E , E ), for several combinations of E , E . Assuming two peak pileup events, these recorded energydistributions result from varying their separation. When E ≥ E (panels (a) - (d)), peak probability occurs at ε = E + E ,indicating this is the most frequent case. The width of each distribution demonstrates that E + E is not a good approximationfor ε in general. A second peak, appearing in panels (e) - (h), indicates that E has a higher probability of being distinguishedas it becomes larger than E . For n th order pileup we require an expression giv-ing the recorded energy within the peak interval.The n +1 pulses will be superimposed and recordedas a single pulse height ε n . This energy would gen-erally be dependent on the energies and separationsof the piled-up events: ε n = ξ [ n (cid:88) i =0 v i f ( t − t i )] ⇒ ε n = ε n ( s , . . . , s n | E , E , E , . . . , E n )for s i = t i − t i − The other expressions of the previous section wouldgeneralize in a similar way, leading to the n th order peak correction P (cid:104) n (cid:105) ( ε n )= (cid:90) dE (cid:90) dE · · · (cid:90) dE n × Pr (cid:104) n (cid:105) peak ( ε n | E , . . . , E n ) S ( E ) · · · S ( E n ) (35)These expressions are complicated and unwieldy, soinstead we use the iterative approximation of [1].This scheme approximates higher order variationsdue to the added random variables by using resultsfrom the previous order.To calculate the n th order correction, the previ-ous order term p (cid:104) n − (cid:105) is used:11 (cid:104) n (cid:105) ( ε n ) = ∞ (cid:90) ∞ (cid:90) Pr (cid:104) n (cid:105) peak ( ε n | E n − , E n ) × p (cid:104) n − (cid:105) ( E n − ) S ( E n ) dE n − dE n (36)with the PDF p (cid:104) n (cid:105) ( ε n ) = ∂∂ε n P (cid:104) n (cid:105) ( ε n ) (37)and the kernel Pr (cid:104) n (cid:105) peak is evaluated using the n th order interval statistics:Pr (cid:104) n (cid:105) peak ( ε n | E n − , E n ) = s ( ε +∆ ε ) (cid:90) s ( ε ) f s | n ( s (cid:48) ) ds (cid:48) (38)We make the approximation that ε n ≈ ε n ( s n | E n − , E n ) and has the same form as ε ( s n | E n − , E n ). f s | n ( s ) is the interval dis-tribution for n events in the peak, equation(15).We can unify the 1 st order correction with (36)by defining p (cid:104) (cid:105) ( E ) = S ( E ). Then equation (36)reduces to (32) for n = 1.
7. Tail effect
The “tail effect” is spectral distortion caused bythe bipolar pulse tail, which reduces pulse heightsoccurring before baseline recovery. In addition toenergy-lowering distortion, energy-dependent lossescan occur if the reduced peak is below a lower-levelthreshold. The model of [1] presents a clear and rel-atively accurate model for tail effects in the outputspectrum. However, in their method there is notthe intervening deadtime region (interval ‘B’). Formodeling purposes this means tail subtraction ef-fects are taken to depend only on the peak state. Inthe present context we must consider the events inthe peak and deadtime, since pulses in either regioncontribute to a negative tail when the instrumentagain becomes live. The method of [1] uses furthersimplifications by reducing the number of randomvariables: input tail events are assumed to be uni-formly distributed in time and have the same en-ergy, which is taken to be the mean energy from themodeled input spectrum; i.e., S ( E ) → δ ( E − (cid:104) E (cid:105) ). We present an alternative tail technique whichdeals with the intervening deadtime interval (region‘B’ in figures 1 and 3), and models energy and tim-ing variations of tail pileup events. This method issimilar to that used for peak pileup in that pileuplikelihood functions of the form of equation (31) arederived for measurement in the tail region (region‘C’). Measurement in this region is sensitive to theprevious two, so the likelihood scheme is more com-plex than in the peak case. The number of eventsin A and B affect the C measurement because theirnegative tails combine (figure 3(c)). Additionally,peak pileup in C can occur (figure 3(b)). The re-sulting pulse heights in the latter case are modifiedby both peak and tail effects. We account for thevarious possibilities by isolating an ‘A+C’ effect,which is the case of pulses in A or C only (or both),and a ‘B+C’ effect, which assumes no zeroth eventand pulses only in B or C (or both). For each casea likelihood function for the tail pileup energy iscalculated. The spectrum of the total pulse config-uration, with ‘A+B+C’ dependence, is calculatedby convolving the corresponding ‘A+C’ componentwith the ‘B+C’ likelihood. In region A we employed the function t p ( s , v , v ) to give the time of the first-orderpeak. For tail pileup we make the simplifyingassumption that the location of the summed peak,were it to be measured in C, is approximately atthe peak of the second input pulse, i.e., t p ≈ t p + s .Thus equation (25) becomes (39) ε ≈ ξ (cid:20) ξ − [ E ] f ( t p + s ) + ξ − [ E ] f ( t p ) (cid:21) ≡ ε ( s , E , E ), for tail energiesFor modeling the second tail effect, we also needan expression for peaks measured in B. This is givenin equation (71), which is a smooth joining of thecomplicated function for region A peaks, and thesimplified expression for region C. Figure 6 showsthat this method is reasonably accurate. The recorded tail spectrum depends on the totalpulse the state (cid:104) a, b, c (cid:105) . We define the followingsymbols to represent the spectrum of measured tailenergies associated with a window state:The effect of events only in A upon measurementsin C is first expressed in the (cid:104) a, , c (cid:105) components.12 (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) during deadtime E in the peak E (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) (cid:61) P r (cid:72) (cid:182) (cid:200) E , E (cid:76) (cid:72) a (cid:76) (cid:72) b (cid:76) (cid:72) c (cid:76) (cid:182) (cid:64) keV (cid:68) P r (cid:72) (cid:182) (cid:200) E , E (cid:76) (cid:72) d (cid:76) (cid:72) e (cid:76) (cid:72) f (cid:76) (cid:72) g (cid:76) (cid:182) (cid:64) keV (cid:68) (cid:72) h (cid:76) Figure 8: Recorded energy distributions for the two tail effects. In each plot the lighter histogram plots Pr A + C ( ε | E , E ),which assumes E is the peak energy, and the separation s ranges from τ A + τ B to τ A + τ B + τ C . The darker histogram plotsPr B + C ( ε | E , E ), which assumes E is in the dead region (B), and s varies from τ B to τ B + τ C . Q (cid:104) a,b,c (cid:105) ( ε ) ≡ probability of recordingenergy ε in the tail q (cid:104) a,b,c (cid:105) ( ε ) ≡ PDF of recorded energy, ∂∂ε Q (cid:104) a,b,c (cid:105) Table 1: Symbols used for the distribution and PDF ofrecorded energy of the tail events in a state
For c = 0 there is no tail count and thus no tailmeasurement, so such terms are zero. For a + c > Q (cid:104) , , (cid:105) → Q (cid:104) , ,c (cid:105) → Q (cid:104) a, ,c (cid:105) → Q (cid:104) a,b,c (cid:105) . t p s s' Τ A Τ B Τ C - - Μ s V H t L Figure 9: (cid:104) , , (cid:105) pileup event. t p ≈ t p + s for the tail pulse. s (cid:48) is distributed according to equation (17). .2.1. First tail effect, Q (cid:104) a, ,c (cid:105) ( ε )We begin by calculating the likelihood of record-ing energy ε in region C, in the event of tail pileup.Here we assume a zeroth event to define the pulsewindow, and additional events in A and C. For firstorder, a = 0 and c = 1. Figure 9 depicts this con-figuration.Recall that the interval PDF for a single tail eventis f s | c =1 ( s (cid:48) ) = τ C , where s (cid:48) = s − ( τ A + τ B ) (equa-tion (17)). However the full separation s is used tocalculate the measured energy using equation (39).Assuming E in τ A , E in τ C , the probability oftail pileup into the discrete energy bin [ ε, ε + ∆ ε ] isthen:Pr (cid:104) (cid:105) A + C ( ε | E , E ) = s ( ε +∆ ε ) (cid:90) s ( ε ) f s | c =1 ( s (cid:48) ) ds (cid:48) (40)and the 1 st order spectral component is (41) Q (cid:104) , , (cid:105) ( ε ) = ∞ (cid:90) ∞ (cid:90) Pr (cid:104) (cid:105) A + C ( ε | E , E ) × S ( E ) S ( E ) dE dE For (cid:104) a, , (cid:105) states we approximate the effect of themultiple orders in A by using p (cid:104) a (cid:105) in the integral,and using the Pr (cid:104) (cid:105) A + C kernel, since the interval dis-tribution in C does not depend on the number ofevents in A (section 5.1). The spectrum measuredin the tail for these configurations isFor a ≥ , b = 0 , c = 1, (42) Q (cid:104) a, , (cid:105) ( ε ) = ∞ (cid:90) ∞ (cid:90) Pr (cid:104) (cid:105) A + C ( ε | E (cid:48) , E (cid:48)(cid:48) ) × p (cid:104) a (cid:105) ( E (cid:48) ) S ( E (cid:48)(cid:48) ) dE (cid:48) dE (cid:48)(cid:48) For a = 0 this reduces to equation (41) since p (cid:104) (cid:105) ( E ) = S ( E ). For zero events in the tail, therecan be no recorded tail energy, so q (cid:104) a,b, (cid:105) ( E ) = 0for all a, b .For (cid:104) , , c (cid:105) states and c >
1, we do an iterativeapproximation using the previous order. However,there are now multiple pulses in the region C, andwe approximate their measured energy as a case of( c − th order peak pileup. But now the primary and additional event have a negative baseline, andso we approximate them using q (cid:104) , ,c − (cid:105) as the pri-mary peak distribution , and q (cid:104) , , (cid:105) as the spec-trum of additional events. Then we use the kernelPr (cid:104) c − (cid:105) peak from equation (38). This is an approxima-tion since Pr (cid:104) n (cid:105) peak models pileup in the interval τ A instead of τ C . However multiple events in τ C tendto be clustered due to the distribution I (cid:104) n (cid:105) , and sothe small difference between interval sizes τ A and τ C is not a major source of error. Then the (cid:104) , , c (cid:105) tail spectrum isFor a = 0 , b = 0 , c >
1, (43) Q (cid:104) , ,c (cid:105) ( ε )= ∞ (cid:90) ∞ (cid:90) Pr (cid:104) c − (cid:105) peak ( ε | E (cid:48) , E (cid:48)(cid:48) ) × q (cid:104) , ,c − (cid:105) ( E (cid:48) ) q (cid:104) , , (cid:105) ( E (cid:48)(cid:48) ) dE (cid:48) dE (cid:48)(cid:48) Finally, for (cid:104) a, , c (cid:105) states, we use the a th orderpeak spectrum with the q (cid:104) , ,c (cid:105) ( ε ) just calculated:For a > , b = 0 , c > Q (cid:104) a, ,c (cid:105) ( ε ) = ∞ (cid:90) ∞ (cid:90) Pr (cid:104) c (cid:105) A + C ( ε | E (cid:48) , E (cid:48)(cid:48) ) × p (cid:104) a (cid:105) ( E (cid:48) ) q (cid:104) , ,c − (cid:105) ( E (cid:48)(cid:48) ) dE (cid:48) dE (cid:48)(cid:48) (44)where the c th order kernel is approximated asPr (cid:104) c (cid:105) A + C ( ε | E a , E c ) = s ( ε +∆ ε ) (cid:90) s ( ε ) f s | c ( s (cid:48) ) ds (cid:48) (45)and the recorded energy is equation (39) with E a , E c replacing E , E and f s | c is the tail inter-val density, equation (17). Q (cid:104) a,b,c (cid:105) ( ε )The second effect models the case when events inB pileup and create a large negative peak in C. It Much like the a th order peak iteration uses p a − ( ε ) inthe primary position. p s s' Τ A Τ B Τ C - - Μ s V H t L Figure 10: Configuration used to calculate Pr (cid:104) (cid:105) B + C ( ε | E , E ),which assumes no zeroth event. t p ≈ t p + s for the tail pulse. s (cid:48) is distributed according to equation (17). is only present at orders ≥ b random events in B using the peak spectrum termscalculated for order b − (cid:104) c =1 (cid:105) B + C by varying s (cid:48) .Equation (39) is used for the recorded energy, withthe separation given as s = s (cid:48) + τ B , with s (cid:48) dis-tributed according to equation (17). The probabil-ity of the pulse height lowering into [ ε, ε + ∆ ε ], duejust to an event in B, isPr (cid:104) c (cid:105) B + C ( ε | E b − , E c ) = s ( ε +∆ ε ) (cid:90) s ( ε ) f s | c ( s (cid:48) ) ds (cid:48) (46)for 0 ≤ s ≤ τ C .Finally, the additional event in C is an event fromthe input spectrum, shifted down due to tails fromA, but also shifted up due to pulses in C. In otherwords, it is approximately distributed as q (cid:104) a, ,c (cid:105) ( ε ).Thus we use these two terms with the kernel aboveto calculate the second componentFor a ≥ , b > , c > Q (cid:104) a,b,c (cid:105) ( ε ) = ∞ (cid:90) ∞ (cid:90) Pr (cid:104) c (cid:105) B + C ( ε | E (cid:48) , E (cid:48)(cid:48) ) × p (cid:104) b − (cid:105) ( E (cid:48) ) q (cid:104) a, ,c (cid:105) ( E (cid:48)(cid:48) ) dE (cid:48) dE (cid:48)(cid:48) (47) This spectral component clearly depends on thetotal positive signal from both A and B. Thus wehave an approximation for the highest order statesusing the succession of terms, Q (cid:104) , , (cid:105) → Q (cid:104) , ,c (cid:105) → Q (cid:104) a, ,c (cid:105) → Q (cid:104) a,b,c (cid:105) .
8. Probability of pileup events
Since the time intervals A, B, and C do not over-lap, the number of events in each interval is inde-pendent of the other two. Furthermore, they arePoisson distributed random variables. We assume λ is constant throughout the detection of the inputspectrum S ( E ). The probability of each state isthen Pr(state | rate = λ )= Pr( a ∈ A and b ∈ B and c ∈ C | λ )= ( λτ A ) a ( λτ B ) b ( λτ C ) c a ! ∗ b ! ∗ c ! e − λ ( τ A + τ B + τ C ) (48) ≡ Pr( (cid:104) a, b, c (cid:105)| λ )The rate λ is separately modeled so the conditioncan be ignored and the expression for state proba-bility isPr( (cid:104) a, b, c (cid:105) ) = ( λτ A ) a ( λτ B ) b ( λτ C ) c a ! ∗ b ! ∗ c ! e − λ ∆ τ for τ A + τ B + τ C = ∆ τ (49)This gives the probability of having non-overlapping pulse intervals of width ∆ τ with thestate (cid:104) a, b, c (cid:105) . These states are distinct , meaninga single pulse window ∆ τ can only be ‘in’ or ‘de-scribed by’ a single state. Such propositions definea state space , which unifies the total time processwith the total set of configurations. In this schemethe total exposure time is decomposed into pulsewindows of width ∆ τ , and since every window hasa state, the temporal composition is equivalent toa superposition of independent states. The fractionof windows with order k pileup is (50)( λ ∆ τ ) k k ! e − λ ∆ τ = k (cid:88) i =0 k − i (cid:88) j =0 Pr( (cid:104) i, k − ( i + j ) , j (cid:105) )where the right-hand-side is the sum over the var-ious combinations with a + b + c = k . In the next15ection it will become clear why each order is sub-divided. The terms above are treated as expan-sion coefficients with the spectral components ap-pended.The total time process is just the superpositionof independent states of order k , and thus ∞ (cid:88) k =0 ( λ ∆ τ ) k k ! e − λ ∆ τ = ∞ (cid:88) a,b,c Pr (cid:104) a, b, c (cid:105) = 1 (51)For example, if there are N events incident dur-ing a long exposure time, the average number ofnon-overlapping windows ∆ τ with configuration (cid:104) n A , n B , n C (cid:105) is equal to N ∗ Pr (cid:104) n A , n B , n C (cid:105) .This section has assumed that the pulse windows∆ τ are non-overlapping. In the next section, weconsider why the total spectrum is not correctlyreconstructed under this assumption. The error re-sults from the fact that tail events require addi-tional recovery time extending past the initial pulsewindow. We will develop a semi-empirical proba-bility expression that models pulse extension due totail events by overlapping two windows.
9. Full correction expansion
The final correction is a combination of the P and Q terms and their associated state probabili-ties. We first describe the total observed spectrum(i.e., the total process), as a superposition of in-dependent (non-overlapping) pulse measurements.We then devise a correction term based on the factthat measurement states overlap, due to the possi-bility of recording a tail event. Adopting the assumptions that pileup states areindependent, the total spectrum can be written asa superposition of the peak and tail componentsderived given in sections 6.3 and 7.2. For k th orderpileup, a + b + c = k , and the k events can becombined into A, B, C to give the set of possiblepileup states Ω k = { (cid:104) a, b, c (cid:105) : a + b + c = k } . Thetotal number of states per order k is ( k +1)( k +2)2 .The k th order spectrum is an expansion into thespectral components with the state probability ascoefficients. However, since there are k + 1 eventsand a maximum of two can be recorded (one peak,one tail), each p, q term is multiplied by k +1 . Now the k th order term is written (52) f ( k ) ( ε ) = 1 k + 1 k (cid:88) i =0 k − i (cid:88) j =0 Pr (cid:104) i, ∆ kij , j (cid:105)× (cid:110) p i ( ε ) + q i, ∆ kij ,j ( ε ) (cid:111) (53)= 1 k + 1 (cid:110) f ( k ) p ( ε ) + f ( k ) q ( ε ) (cid:111) where ∆ kij ≡ k − ( i + j ) p (cid:104) (cid:105) ( ε ) = S ( ε ) q (cid:104) a,b, (cid:105) ( ε ) = 0and the p , q terms have been collected: f ( k ) p ( ε ) = k (cid:88) i =0 k − i (cid:88) j =0 Pr (cid:104) i, ∆ kij , j (cid:105) p i ( ε ) (54) f ( k ) q ( ε ) = k (cid:88) i =0 k − i (cid:88) j =0 Pr (cid:104) i, ∆ kij , j (cid:105) q i, ∆ kij ,j ( ε ) (55)The total spectrum (PDF), assuming indepen-dent pulse states, up to order n , is f ( ε ) = n (cid:88) k =0 f ( k ) ( ε ) (56)= n (cid:88) k =0 k + 1 (cid:110) f ( k ) p ( ε ) + f ( k ) q ( ε ) (cid:111) (57)The approximation order n is the value at which ( λ ∆ τ ) n n ! e − λ ∆ τ becomes negligible. This expansion is accurate at low rates, when theprobability of having two adjacent pileup states issmall. But we must account for the fact that a tailmeasurement causes a new pulse interval to ‘begin’before the end of the first. A theoretically correctconstruction from the partitioned light curve (i.e.,the total process) becomes rather complex due tothis fact. Each state with a tail count, by defini-tion, also contains the ‘beginning’ of the next pulseinterval. Events occurring within ∆ τ of the firsttail count must also be modeled. In terms of thestate spectral components p and q , the total mea-surement from this joined-pulse interval contains apeak count, a tail count, and a possible third tailcount (from the adjoined pulse).16 (cid:68)Τ Τ(cid:68)Τ Τ(cid:68)Τ Τ(cid:68)Τ Τ(cid:68)Τ Τ(cid:68)Τ (cid:45) (cid:45) (cid:64) Μ s (cid:68) V Figure 11: Sample output V ( t ) showing how exposure timeis partitioned into pulse states, with measured energies dis-tributed as p a ( ε ) & q abc ( ε ). Subsequent pulse heights within∆ τ of a tail measurements are distributed as q abc ( ε ). A com-bination of independent states must be adjusted for the factof overlapping regions. monte carloinput spectrumno overlapw (cid:144) overlap10 - - - - - P D F @ k e V - D ΛΤ » Λ = ´ cps Energy @ keV D - - E rr o r @ % D Figure 12: Comparison with Monte Carlo simulation ofthe model spectrum assuming non-overalapping pulse states(BLUE, dashed) vs. the spectrum corrected for overlap(BLACK), for GBM. The spectrum is a cutoff power-law, S ( E ) ∼ E . exp( − E/ One way to model this is to consider all possi-ble combinations of states. If Ω is the set of pos-sible states, then we would have to partition theprocess in terms of the product states Ω ⊗ Ω. Atvery high rates it may even be necessary to considerΩ ⊗ Ω ⊗ Ω, since additional pulse can occur with ∆ τ of the third tail measurement. Due to the obviouscomplications, we have developed a semi-empiricalapproximation which adjust the weighting of peak and tail components.In the approximation two overlapping intervalscan be represented as dependent random events M and M (cid:48) , whose probability is generally constructedasPr( M ∪ M (cid:48) ) = Pr( M )+Pr( M (cid:48) ) − Pr( M ∩ M (cid:48) ) (58)In general the adjoined state, which we may iden-tify as M (cid:48) = (cid:104) a (cid:48) , b (cid:48) , c (cid:48) (cid:105) can be begin at any randomtime in the interval τ C , if M contains a tail pulse.Thus M = (cid:104) a, b, c (cid:105) with c ≥ a, b ). Figure11 shows several pulse windows with tail measure-ments. For a given state (cid:104) a, b, c (cid:105) with c ≥
1, thecounts c constitute the peak of a subsequent pileupevent, denoted by the symbol (cid:104) c − , b (cid:48) , c (cid:48) (cid:105) . Thenwe approximate the probability Pr( M ∩ M (cid:48) ) as theprobability of having c counts in the tail of a state M and c − M (cid:48) (with-out regard to the a, b or b (cid:48) , c (cid:48) ) (59)Pr (cid:18) (cid:104) abc (cid:105) ∩ (cid:104) ( c − b (cid:48) c (cid:48) (cid:105) (cid:19) ≈ ( λτ C ) c e − λτ C c ! × ( λτ A ) c − e − λτ A ( c − ≡ Pr( c, a (cid:48) ) δ c − ,a (cid:48) The corrected peak contribution is reduced byone count per configuration of this tail-peak over-lap, so we have a single term subtracted at eachorder, for k >
0, from f ( k ) p : f ( k ) p,Corr ( ε ) = (61)1 k + 1 (cid:110) f ( k ) p ( ε ) − Pr( k + 1 , k ) p k ( ε ) (cid:111) withPr( k +1 , k ) = ( λτ C ) k +1 e − λτ C ( k + 1)! × ( λτ A ) k e − λτ A k ! (62)For k = 0 there is no tail count, thus no need toworry about carry-over effects. Therefore we cansay f ( k =0) p,Corr ( ε ) = f (0) p ( ε ) = e − λ ∆ τ S ( ε )which is the weighted zeroth order term giving thefraction of counts without peak pileup or a trailingdeadtime / tail event.The tail correction adds weight to the contribu-tion of f ( k ) q in the total spectrum, since the adjoinedpulse has its own tail region. The energy of thesecounts has a similar distribution to tail counts of17he initial pulse. Therefore we use the tail terms q , but adjust the associated probability. We modelit as the complement of the overlap probability allconfigurations; i.e., 1 − (cid:80) nk =0 Pr( k + 1 , k ). The cor-rected tail contribution is: f ( k ) q,Corr ( ε ) = (63)1 k + 1 (cid:110)(cid:104) − ∞ (cid:88) k (cid:48) =0 Pr( k (cid:48) + 1 , k (cid:48) ) (cid:105) f ( k ) q + f ( k ) q (cid:111) (64)= 1 k + 1 (cid:104) − ∞ (cid:88) k (cid:48) =0 Pr( k (cid:48) + 1 , k (cid:48) ) (cid:105) f ( k ) q Again for k = 0 there are no tail counts, so f (0) q,Corr ( ε ) = f (0) q ( ε ) = 0. Now the total measuredspectrum is written f ( ε ) = f (0) p ( ε ) + n (cid:88) k =1 k + 1 × (65) (cid:26) f ( k ) p ( ε ) − Pr( k + 1 , k ) p k ( ε )+ (cid:104) − ∞ (cid:88) k (cid:48) =0 Pr( k (cid:48) + 1 , k (cid:48) ) (cid:105) f ( k ) q ( ε ) (cid:27) or simply f ( ε ) = n (cid:88) k =0 (cid:110) f ( k ) p,Corr ( ε ) + f ( k ) q,Corr ( ε ) (cid:111) (66)This method is evaluated by comparing com-puted spectra with Monte Carlo simulations. Fig-ure 12 demonstrates a spectrum computed underthe assumption of independent pulse states, ver-sus the overlap correction. The Monte Carlo isdescribed with more comparisons in section 10.3.Losses due to pileup and subtraction effects can besummarized as follows. If the true detection rate is λ , the predicted counting rate is λ Rec = λ (cid:90) f ( ε ) dε (67)= λ n (cid:88) k =0 (cid:90) f ( k ) ( ε ) dε (68)The relationship between λ Rec and λ is discussedfurther in section 11.
10. Numerical evaluation
In this section we apply the model to GBM andcompare it with Monte Carlo simulations. The firststep is numerical evaluation of the probability like-lihoods, which are used to derive pileup spectralcomponents.
The single-pulse peak time is given by f (cid:48) ( t p ) = 0.Recall that the first-order peak time, defined in gen-eral by equation (21), cannot be found in closedform for the true pulse shape f ( t ). Instead we maketwo approximations of the function t p ( s , v , v ) forregion A and region C, and join them smoothly inregion B. For peak pileup (A), we use an empiri-cal expression with two constant parameters Λ s , Λ v .Events are sampled on a ( s , v , v ) grid, and PHAmeasurement is simulated with the same routinesused in the Monte Carlo. Resulting peak times arefit using a non-linear least-squares technique andthe functional form below: t p ( s , v , v ) = ¯ t p + F ( e Λ s s − e − Λ v ( v − v v v ) (69)where ¯ t p is the weighted average of the zeroth-orderpeak time of each pulse,¯ t p = v t p + v ( t p + s ) v + v = t p + v v + v s (70)and F = 1 µs . The second term adjusts ¯ t p . Thebest-fit parameters are Λ v = 2 . , Λ s = 0 .
40, for s in µs . Of course these parameters, and thefunctional form itself, are dependent on the spe-cific pulse shape. The ones given here correspondto the best-fit GBM pulse shape from section 3.In region C we use the approximation t p ≈ t p + s .For region B these two expressions are smoothlyjoined using a logistic function about s = 1 µs ,which is approximately the positive pulse width.The joining scale is an arbitrary parameter, whichwe set to Λ = 30[ µs ] − . Thus the full peak timeexpression is: t p ( s , v , v ) = (cid:104) ¯ t p + ( e Λ s s − e − Λ v ( v − v v v ) (cid:105) × (cid:16) − g ( s ) (cid:17) + (cid:104) t p + s (cid:105) × g ( s )(71)where g ( s ) = 1 / [1 + exp( − Λ ( s − only (i.e., v , v in A), under certain conditions the initial pulse isaccurately distinguished from the summed signal.This results in a piecewise step in the peak-time for-mula such that t p → t p , observable as discontinuityin region A of figure 6. This generally occurs when v > v and the separations are sufficiently large.An inspection of the simulated peak-time data re-veals that the precise condition is a complicatedfunction of all three variables, however its strongestdependence is on s . By manual estimation we havedetermined the following additional criteria for thepeak-pileup time: if t p > t p ∗ exp(1 . s ), where t p isfrom equation (69), or s > τ A , then set t p → t p .This is approximates the effect of digital bufferingdescribed in section 3, and is only for the case ofpeak-pileup. (cid:104) n (cid:105) peak , Pr (cid:104) n (cid:105) A + C , Pr (cid:104) n (cid:105) B + C Because we are using the more accurate pulseshape, the conditional probability kernels cannotbe calculated in closed form. We calculate themnumerically by first first defining a discrete set ofchannel energies in which probabilities are calcu-lated, E = { E i } . Since pileup events have beendetected by the instrument it’s sufficient to use thechannel definitions of actual data.Defining a discrete time step ∆ s , integrals likePr( ε | E, E (cid:48) ) = (cid:90) f s | n ( s (cid:48) ) ds (cid:48) (72)are evaluated discretely over the finite ∆ s sam-ple ( s = { , ∆ s, s, . . . , l ∆ s } ). At each step ofthe integrand, the function ε ( s | E, E (cid:48) ) is evaluated.Then ε corresponds to one of the channels in E ,and an appropriate lookup algorithm returns its in-dex (channel) i . For E, E (cid:48) also discretely sampledfrom E and corresponding to channels j, k , the con-ditional probability is a triply-indexed discrete ob-ject Pr( ε | E, E (cid:48) ) → Pr[ i, j, k ].This probability ‘array’ is initialized to 0 for all i, j, k , and incremented by a value ∆ P l for each stepin the RHS numerical integration, where l specifiesthe integration step s l to s l + ∆ s , and i is the chan-nel corresponding to ε ( s l | E j , E k ). The intervaldistributions are integrable, so ∆ P l is exact for the l th time step:∆ P l = s l +1 (cid:90) s l nτ n ( τ − s (cid:48) ) n − ds (cid:48) = ( τ − s l ) n − ( τ − s l +1 ) n τ n where τ is the normalization interval (A, B, or C)and s l ∈ [0 , τ ]. The calculation and storage of thePr( ε | E, E (cid:48) ) likelihoods introduces some computa-tional overhead, since we must have one per orderin each interval. But since no information about thesource spectrum or intensity is required, they canbe calculated a priori and stored. A separate codecomputing spectrum corrections can read each datablock as necessary. In our implementation we use128 energy channels and store each function in itsown 3D data set in an HDF5 file [12]. This librarywas chosen for its stability and ease-of-use in storinga set of multi-dimensional arrays and keyword pa-rameters in a single file. Using 4-byte floating pointdata, the total uncompressed requirement for all or-ders up to 5 is (3 ∗ ) × (4 bytes) × We have implemented a Monte Carlo simulationthat includes both the arrival of random eventsand the discrete pulse-height measurements madeby the instrument. We simulate a sequence of ex-ponentially distributed event times over an arbi-trary exposure interval, during which the processintensity is constant. Their energies are sampledfrom an input spectrum S ( E ), which represents therecorded spectrum in the absence of pileup distor-tion. Our simulation includes the processing of sig-nals by the GBM pulse-height electronics, but notphysical interaction of gamma-rays in the detectors.Because the detector response is neglected, the ex-ample input spectra shown in the next section areidealized, e.g., monochromatic gamma-rays wouldnot generate a pure Gaussian input spectrum be-cause, for some of the incident photons, only a por-tion of the energy is deposited in the crystal (Comp-ton scattering or pair production), and the remain-ing energy might escape. This information is sum-marized in the detector response information thatconverts an external gamma-ray spectrum into theinput spectrum S ( E ) (equation (33)), and is not afocus of this paper.19he analytical pileup model is calculated up toorder 5. For k >
5, an approximation is used bysubstituting fifth-order likelihoods in for higher or-ders in the iterative calculation. At high orders thisresults in a slight overestimation for states having a + b > (cid:29) c . Monte Carlo simulations are exe-cuted for 200,000 input events, at several differentrates. This simulation size is large enough to limitstatistical fluctuations in the output spectra, allow-ing a more systematic comparison with the analyt-ical model.An instructive case is that of a narrow line spec-trum. We simulate a pure Gaussian line with µ = 2 . M eV and σ = 0 . µ , and compare it withthe model prediction (figure 13 ). At low rates theamount of spectral distortion is quite small, but asthe rate increases the line becomes distorted. Be-cause p a ( ε ) and q abc ( ε ) components are indepen-dent of the rate (only their expansion coefficientsvary), a range of input rates can be tested for agiven S ( E ) without much difficulty.Plot residuals (errors) are calculated relative tothe Monte Carlo output. If n i is the number ofcounts in the channel i of the output MC spec-trum, and m i is the analytical prediction (i.e., themodel), residuals are calculated r i = ( m i − n i ) /n i (times 100). The model is reasonably accurate forthe Gaussian spectrum, and improves for the otherspectra shown. For spectral shapes S ( E ) domi-nated by narrow-band features, such as figure 13where the line occupies about 1.6% of the total volt-age range, the model error tends to be higher inchannels away from the line, where the only countsare due to high order pileup. To the left of theline are primarily tail distorted counts, and to theright are mainly counts suffering from peak pileup.However for broad-band models, such as the one infigure 14, the error is much smaller and more uni-form. The notable difference is that more channelshave both peak and tail pileup counts when theinput spectrum is broad, suggesting that the peakand tail modeling assumptions have complementaryerrors, which tend to cancel out when combined.Figures 15 and 16 demonstrate the effects of pulsepileup when the zeroth order count spectrum is asmooth, cut-off power-law. Low-energy spectral in-dices become flattened as lower energy pulses areshifted up through peak-pileup or removed throughtail subtraction. The consequence is that pileupcan significantly affect inferences about the power-law index, peak or cut-off energy parameters, andmeasures of spectral hardness unless corrected. In general count spectra will become harder as the in-put rate increases.Figure 17 is a fictitious spectrum demonstrat-ing four emission lines without a continuum. Evenwithout a continuum lines become distorted at highrates. Additional false lines can appear from peakpileup of each line source.Finally, the energy binning used to generate plotsin this section are approximately those of the GBMBGO detectors. However, the pileup process occursin voltage space. The model MC Model Rate (cid:72) cps (cid:76) P D F (cid:64) k e V (cid:45) (cid:68) Energy (cid:64) keV (cid:68) (cid:45) (cid:45) E rr o r (cid:64) (cid:37) (cid:68) Figure 13: Model vs. simulation comparison for several inputrates. The input spectrum is a single Gaussian-shaped lineat 2.2 MeV, and is shown by the black dashed line. At highrates the line is completely distorted.
11. Pileup losses
As already mentioned in the opening sectionspulse-pileup introduces counting losses which ex-ceed the expectation from the conventional non-paralyzed deadtime correction. By simple argu-ments this function is λ Rec = λ / (1 + λ τ ), andis algebraically invertible to give λ = λ Rec / (1 − λ Rec τ ) (73)At high rates it is insufficient both because ofpileup, and because the derivation assumes λ Rec < /τ . [1], [9]Pileup changes the picture for two reasons. Thefirst is that it randomly extends instrumental dead-time. In GBM this occurs only when events are20 C Model Rate H cps L P D F @ k e V - D Energy @ keV D - - E rr o r @ % D Figure 14: A spectrum of two lines (4 and 6 MeV) on a cut-off power-law continuum. The model becomes more accuratefor spectra with a continuum, suggesting that energy raisingand lowering associated with peak and tail modeling havecanceling errors.
MC Model Rate H cps L - - - - - - P D F @ k e V - D Energy @ keV D - - E rr o r @ % D Figure 15: A cut-off power law model, whose PDF is ∼ E − . exp( − E/ within τ A of each preceding event, since the addi-tional detections delay peak measurement. The sec-ond is that the tail effect causes energy-dependentbaseline subtraction losses, when a tail shifts smallpulses below the recording threshold. Such lossesare more significant for ‘flat’ or broadband ener-getic spectra because they entail a mixture of large MC Model Rate H cps L - - - - - - P D F @ k e V - D Energy @ keV D - - E rr o r @ % D Figure 16: A cut-off power law model, whose PDF is ∼ E . exp( − E/ MC Model Rate H cps L - - - - - - P D F @ k e V - D Energy @ keV D - - E rr o r @ % D Figure 17: A spectrum of four fictitious emission lines, withno continuum. and small pulses, and the smaller ones are morelikely to be lost. In GBM these additional lossesare mitigated by the imposed 2.6 µs deadtime, butat high rates they must be accounted for to fullycorrect the relation λ Rec ( λ ).The first effect (paralyzable deadtime) by itselfcan be predicted, to first order, using Poisson statis-tics and employing a root-finding algorithm. The21rocedure is described in [9], and yields a numericalsolution. We do not explicitly use this treatment inour model. Rather, because the model’s predictedlosses agree well with simulation, we conclude thatthe method of overlapping windows, equation (65),is an effective proxy for paralyzable deadtime. Thisseems reasonable since the subtracted overlap termsin equation (61) reduce the proportion of peak mea-surements, which is roughly the same effect as par-alyzable deadtime. An alternative interpretation isthat severe paralyzable deadtime is improbable inGBM even at the high rates tested, due to the pulseshape and peak finding algorithm. At rates beyond10 cps, it is likely that the instrument becomesnon-linear and our assumptions would not hold.The second effect, tail subtraction loss, is some-what more serious than the first in our case, be-cause the pulse tail is longer than the peak and thenegative amplitude is large. Because of spectraldependence, tail losses can only be predicted byassuming a zeroth order pulse-height distribution S ( E ). These losses can be investigated by plottingmarginal distributions from the tail likelihoods cal-culated in section 7. Figures 20 and 21 show thefirst-order case of the ‘A+C’ and ‘B+C’ pileup sce-narios. When tail spectral components are calcu-lated using probability integrals like equation (47),they turn out to have total probability less thanone ( (cid:80) ε Q ( ε ) < λ Rec ( λ ), demonstrating tail losses due torate and spectral shape.At present an analytical inversion giving λ ( λ Rec )with deadtime and pileup has not been found,though in principle one exists until the the turnoverin figure 19. Future efforts using a technique suchas series reversion of equation (65) may be fruit-ful. However we note that spectral fitting with thepileup correction is itself an inversion process andresults in a possible solution for λ given λ Rec fromreal data.
12. Conclusions
We have derived a model which accurately pre-dicts the recorded number and spectrum of a con-stant intensity Poisson process with pulse-pileup,
MC simulationmodel Λ @ Μ s - D F r ac ti on r ec o r d e d Λ Τ Figure 18: Fraction recorded vs. detection rate. when compared to Monte Carlo simulations. Wehave used the peak modeling technique and iter-ation method of [1], and extended the treatmentto a three-region bipolar pulse. The novelty in thismethod is that it provides a way to model input en-ergy and timing statistics of tail pileup events. Thetotal spectrum is written as a state-space expansionof overlapping pulses. The technique generally ap-plies to bipolar shaping instruments, and has beendemonstrated using the true pulse shape for GBM.
13. Acknowledgements
This work is supported primarily by funds fromthe Fermi-GBM project, with additional supportfrom the Fermi Guest Investigation (GI) programon Terrestrial Gamma-ray Flashes. The authorswould also like to acknowledge the GBM instrumentteam and builders for their detailed specification ofthe detectors and electronics.
14. ReferencesReferences [1] K. Taguchi, E. C. Frey, X. Wang, J. S. Iwanczyk,W. C. Barber, An analytical model of the effects ofpulse pileup on the energy spectrum recorded by en-ergy resolved photon counting x-ray detectors, MedicalPhysics 37 (2010) 3957–3969. (cid:72) E (cid:76) (cid:181) E (cid:45) (cid:227) (cid:45) (cid:73) E (cid:77) S (cid:72) E (cid:76) (cid:181) E (cid:227) (cid:45) (cid:73) E (cid:77) S (cid:72) E (cid:76) (cid:181) E (cid:45) S (cid:72) E (cid:76) (cid:181) E (cid:45) Λ (cid:43) Λ Τ Λ (cid:64) Μ s (cid:45) (cid:68) Λ R ec (cid:64) Μ s (cid:45) (cid:68) Λ Τ Figure 19: Relation between the detection rate λ and theobserved rate λ Rec , in millions of counts per second. As thedetection rate increases, a simple deadtime approximationbecomes insufficient to explain the recorded rate. Moreover,energy dependent losses occur due to the tail effect, whichare worse for broadband, energetic (“flat”) spectra becausethey generate small and large pulses. Pulse pileup greatlyincreases uncertainty about the true rate given the observedrate. As the input rate exceeds 2 ∗ cps ( λ τ ≈ . λ given λ Rec becomes difficult (only a lower limitcan be determined). Note, since GBM detectors have theirown counting electronics, λ is the rate in a single detector.Cut-off energies of 1 MeV and 7 MeV are roughly analagousto 20 keV and 150 keV in NaI detectors, respectively.[2] D. Cano-Ott, J. Tain, A. Gadea, B. Rubio, L. Batist,M. Karny, E. Roeckl, Pulse pileup correction of largenai(tl) total absorption spectra using the true pulseshape, Nuclear Instruments and Methods in PhysicsResearch Section A: Accelerators, Spectrometers, De-tectors and Associated Equipment 430 (1999) 488 – 497.[3] Y. Danon, B. Sones, R. Block, Dead time and pileupin pulsed parametric x-ray spectroscopy, Nuclear In-struments and Methods in Physics Research Section A:Accelerators, Spectrometers, Detectors and AssociatedEquipment 524 (2004) 287 – 294.[4] C. Meegan, G. Lichti, P. N. Bhat, E. Bissaldi,M. S. Briggs, V. Connaughton, R. Diehl, G. Fishman,J. Greiner, A. S. Hoover, A. J. van der Horst, A. vonKienlin, R. M. Kippen, C. Kouveliotou, S. McBreen,W. S. Paciesas, R. Preece, H. Steinle, M. S. Wallace,R. B. Wilson, C. Wilson-Hodge, The fermi gamma-rayburst monitor, The Astrophysical Journal 702 (2009)791.[5] E. Bissaldi, A. von Kienlin, G. Lichti, H. Steinle, P. N.Bhat, M. S. Briggs, G. J. Fishman, A. S. Hoover, R. M.Kippen, M. Krumrey, M. Gerlach, V. Connaughton,R. Diehl, J. Greiner, A. J. van der Horst, C. Kou-veliotou, S. McBreen, C. A. Meegan, W. S. Paciesas, E > Ch.20 E > Ch.60 E > Ch.80 E > Ch.100 E > Ch.128 E channel F r ac ti on l o s t Channel energy @ keV D Figure 20: Baseline subtraction losses due to the (first or-der) A+C effect. E is energy of the primary count. E isthe input energy of the second count, which occurs in thetail of E (i.e., region C). This is the same configuration asdepicted in figures 2(c) and 9. As the peak amplitude E increases, its tail becomes more negative. If E is too small,the summed peak in C is below registration threshold and E is lost. 100% loss occurs when E (cid:29) E . Each curve is cal-culated directly from the tail likelihood by integrating overthe recorded pulse height (sum over channels in the discretecase): P r loss ( E , E ) = (1 − (cid:80) ε P r A + C ( ε | E , E )) E > Ch.20 E > Ch.60 E > Ch.80 E > Ch.100 E > Ch.128 E channel F r ac ti on l o s t Channel energy @ keV D Figure 21: Losses in the B+C effect (figure 10). Eachcurve is calculated by integrating the measurement prob-ability in the event of E , E : P r loss ( E , E ) = (1 − (cid:80) ε P r B + C ( ε | E , E )).R. D. Preece, C. A. Wilson-Hodge, Ground-based cal-ibration and characterization of the Fermi gamma-rayburst monitor detectors, Experimental Astronomy 24(2009) 47–88.[6] A. J. van der Horst, C. Kouveliotou, N. M. Gorgone,Y. Kaneko, M. G. Baring, S. Guiriec, E. G¨o˘g¨u¸s, J. Gra-not, A. L. Watts, L. Lin, P. N. Bhat, E. Bissaldi,V. L. Chaplin, M. H. Finger, N. Gehrels, M. H. Gibby,M. M. Giles, A. Goldstein, D. Gruber, A. K. Harding,L. Kaper, A. von Kienlin, M. van der Klis, S. McBreen, . Mcenery, C. A. Meegan, W. S. Paciesas, A. Pe’er,R. D. Preece, E. Ramirez-Ruiz, A. Rau, S. Wachter,C. Wilson-Hodge, P. M. Woods, R. A. M. J. Wijers,Sgr j1550–5418 bursts detected with the fermi gamma-ray burst monitor during its most prolific activity, TheAstrophysical Journal 749 (2012) 122.[7] M. S. Briggs, G. J. Fishman, V. Connaughton, P. N.Bhat, W. S. Paciesas, R. D. Preece, C. Wilson-Hodge,V. L. Chaplin, R. M. Kippen, A. von Kienlin, C. A.Meegan, E. Bissaldi, J. R. Dwyer, D. M. Smith, R. H.Holzworth, J. E. Grove, A. Chekhtman, First results onterrestrial gamma ray flashes from the fermi gamma-rayburst monitor, J. Geophys. Res. 115 (2010) A07323.[8] L. Theis, S. Persyn, M. Johnson, K. Smith, B. Walls,M. Epperly, Development of a high-speed multi-channelanalog data acquisitioning architecture, in: AerospaceConference, 2006 IEEE, 2006, p. 9 pp. doi: .[9] W. Leo, Techniques for Nuclear and Particle PhysicsExperiments, revised 2nd ed., Springer-Verlag, 1994.[10] G. Grimmet, D. Stirzaker, Probability and RandomProcesses, Clarendon Press, Oxford, 1992.[11] G. Knoll, Radiation Detection and Measurement, 2nded., Wiley, New York, 1989.[12] .-. The HDF Group. Hierarchical data format version 5,2012, Hdf5, URL: ..