- Email: [email protected]
A glow curve analyzer (GCA) for routine analysis of personnel thermoluminescent dosemeter results
Radiation Measurements 43 (2008) 621 – 625 www.elsevier.com/locate/radmeas
A glow curve analyzer (GCA) for routine analysis of personnel thermoluminescent dosemeter results W.J. Chase a,∗ , M.D. Bezaire b , F.P. Vanderzwet c , C.E. Taylor a a Health Physics Department, Ontario Power Generation, Whitby, Ont., Canada b Department of Applied Health Sciences, University of Waterloo, Waterloo, Ont., Canada c Bruce Power, P.O. Box 1540, Building B13, Tiverton, Ont., Canada N0G 2T0
Abstract A glow curve analyzer (GCA) spreadsheet has been developed using Microsoft 䉸 Excel姠 to perform glow curve analysis on thermoluminescent dosimeter (TLD) data from a personnel dosimetry system. The TLD data come from cards with four LiF:Mg,Ti chips that have been annealed and therefore have a simple glow peak structure. GCA removes spikes in the glow curve data, and then smoothes it. After select start and end points for the glow peak, it fits a Boltzmann function to represent the glow curve signal background under the glow peak. The Boltzmann function is subtracted and two Weibull curves are fit to the remaining net signal between the start and end points. The first Weibull curve is fit to peak 5, and the second one to any small remaining contribution from peaks 3 and 4 or from contaminants. The sum of the two Weibull curves is the glow curve signal result. GCA provides rapid review and correction of all glow curves, improving the quality of the results and reducing the time required for complete processing of official dose results. © 2008 Elsevier Ltd. All rights reserved. Keywords: TLD; Glow curve analysis; Weibull distribution; Boltzmann function
1. Introduction Ontario Power Generation (OPG) and Bruce Power (BP) use the Harshaw 8828 badge case and a four-element thermoluminescent dosimeter (TLD) card for personnel dosimetry. Elements 1 and 2 of the badge are 100 mg cm−2 LiF:Mg,Ti chips enriched in 7 Li (Harshaw TLD-700). Element 1 is filtered by 1010 mg cm−2 of low Z material, Element 2 is filtered by 1010 mg cm−2 of tin and Teflon, and they are used to measure photon signals. Elements 3 and 4 are 40 mg cm−2 TLD-700 chips with 17 mg cm−2 of filtration for Element 3 and 74 mg cm−2 for Element 4. They are used to measure the beta signals (see Chase and Hirning, 2008 for further information). The cards are annealed in an oven at 80 ◦ C for 17 h before use to minimize fading, by significantly reducing the size of the low temperature TL peaks. As a result peak 5 is the main peak, with a small contribution from peaks 3 and 4.
∗ Corresponding author. Tel.: +1 905 430 2215x3242; fax: +1 905 430 8583.
E-mail address: [email protected] (W.J. Chase). 1350-4487/$ - see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.radmeas.2008.02.020
For each card element, the TL glow curve signal is integrated over a defined region to obtain the calibration region integral (CRI), corrected for the chip sensitivity. The system is calibrated by exposing cards free-in-air to a 137 Cs source so that TL signals are converted to units of exposure in mR by means of element-specific reader calibration factors (RCFs). Element-specific values for the background radiation exposure are subtracted from the field card element TLD readings. The net values are processed with an algorithm to obtain personal dose equivalent results for Hp (10) and Hp (0.07). The original operating mode used the vendor-supplied software to review all glow curves for cards with doses that would be assignable (Hp (10) or Hp (0.07) 0.095 mSv (9.5 mrem)) and determine an adjustment as necessary. Most adjustments were required for the glow curves from Elements 3 and 4. These glow curves are about 2.5 times smaller in absolute terms than the glow curves from Elements 1 and 2, and thus usually noisier. These glow curves often have a significant high-temperature tail under and after the main glow peak, which required manual adjustment to correct the value for the glow peak area. Fig. 1 shows a typical glow curve for Element 3 or 4, with spikes
622
W.J. Chase et al. / Radiation Measurements 43 (2008) 621 – 625
3. Program description
1.6 Calibration Region 1.4
The main functionality and display of glow curves in GCA is contained in one worksheet of GCA, although many other worksheets are used for processing data or displaying the results for a single element. The main functions of GCA are either initiated on the main worksheet or selected from an additional top-level “GCA” menu item. GCA also includes the algorithm used to process the four net element results and convert them into the personal dose equivalents, Hp (10) and Hp (0.07). Many of the details of the GCA functions are specific to the OPG and BP TLD systems, but the GCA glow curve analysis should be applicable to TLD glow curves with a single main peak.
Adjustments
TL Signal
1.2 1.0 0.8 0.6 0.4 0.2 0.0 1
21
41
61
81
101 121 Channel
141
161
181
Fig. 1. Element 4 glow curve requiring correction to remove spikes (left-hand set of bars) and tail after glow peak (right-hand set of bars). The calibration region integral (CRI) of the TL signal can be manually adjusted by the signal in the adjustment regions.
before the glow peak and a high-temperature tail after the peak. The left-hand set of vertical bars shows how the effect of the spikes could be determined, while the right-hand set of bars shows how the tail after the glow peak up to the end of the calibration region could be estimated. However, using this manual method it is not readily possible to determine the value of the glow curve background signal under the glow peak. Our interest was primarily in accurately determining the glow peak area after subtracting any spikes and any residual high-temperature signal under and after the glow peak. After reviewing a number of options, it was determined that glow curve analysis would provide the most reliable results. BP in collaboration with OPG developed a computer program called glow curve analyzer (GCA). The program has been evaluated (Bezaire, 2003), and accepted by the Canadian Nuclear Safety Commission (CNSC) to meet the Canadian regulatory requirements for dosimetry service providers (CNSC, 2006). 2. Choice of glow curve analysis functions and processing system OPG and BP each read more than 10,000 TLD cards per processing period (quarterly or monthly), with many cards having low-level signals. The analysis method and functions chosen must therefore not rely on user input of starting values or require numerous iterations to obtain the final values. Pagonis et al. (2001) demonstrated that the Weibull function could successfully model the TL intensity for peaks exhibiting first-order kinetics. After considering various options, we decided to develop our own GCA using Microsoft 䉸 Excel姠 . In addition to its programming capability through Visual Basic for Applications (VBA), Excel can drag and drop points on a graph by the use of a “Data Marker”, allowing ready changes to certain parameters. We derived parameters directly from smoothed glow peak data that allowed us to readily determine suitable parameters for calculating the Weibull distribution, and did not require subsequent iterations to improve the fit.
4. Glow curve analysis 4.1. Despiking The value for each channel of glow curve data is compared to the average for the four nearest channels, and replaced by the average if it exceeds a certain value. Because there are cases where spikes two channels wide are not removed, the despiking routine is run again on the despiked data. 4.2. Smoothing The despiked data must be smoothed to properly determine certain “critical points” that are key to determining certain points on the glow curve, such as the start and end of the glow peak. The smoothing function is applied to the despiked data points compressed by adding 1 and then taking the natural logarithm. The smoothing routine works by fitting a running 15-point quadratic least squares fit, centered on the channel in question, to each of the 200 points. The parabola is fit using the method of least squares and the center of the 15-point parabola is then plotted as the new point. A second pass of the previously smoothed data is done to further smooth the curve and ensure no discontinuities remain. This is crucial, as the first and second derivatives of the glow curve data at every channel are required to determine the critical points. The GCA smoothing routine is identical to the Savitsky–Golay smoothing routine applied to points equally spaced in x (Press et al., 1992). After the second pass, the data are transformed back to a linear scale by reversing the original transformation. 4.3. Determination of start and end points Two “critical points” GCA determines are the start and end points for the glow peak. The calculation of these is shown in Fig. 2. The start point (xsp , ysp ) is the point where the first derivative has first reached a certain fractional value (the start point factor) of the absolute maximum of the first derivative. In effect, the start point is the location where the tangent of the glow peak has risen above a certain angle. The calculation of the end point (xep , yep ) starts with the location of the channel (channel a) of the absolute minimum
W.J. Chase et al. / Radiation Measurements 43 (2008) 621 – 625
623
b
0.04
0.7
0.03
0.6
0.02
0.5
0.01
0.4
0
End End Point Point
Start StartPoint Point
0.3
-0.01 0.2
a
-0.02
0.1
-0.03
0
-0.04 1
25
50
75
100
125
150
175
200
1
25
50
75
100
125
150
175
200
Fig. 2. A typical example of the first and second derivatives. The two vertical lines represent the locations of the start and end points. In the lower left figure and the inset figure, “a” and “b” point to the first derivative and the other line is the second derivative. In the right-hand figure, note the start and end point positions with respect to the smoothed data.
E1
60
E2
E3
E4
E1/E2 Boltzmann
E3/E4 Boltzmann
50 Charge (pC/16)
of the first derivative. Next, the first channel to the right of this (channel b) where the first derivative equals 0 is located. Lastly, the first channel to the right of channel b where the second derivative passes through zero from positive to negative is located. This is an inflection point, and one channel to its left is taken as the end point of the glow peak, (i.e. In other words, the end point is the channel just before the first channel to the right of the peak where there is an inflection point and a positive slope).
40 30 20 10
4.4. Boltzmann function for background signal under glow peak A key reason for creating GCA was to eliminate the residual signal under and after the glow peak. A main reason for this residual signal is the emission of black body radiation by the chips. As the chip temperature increases, the black body output of light increases exponentially. Many glow curve analysis and deconvolution packages use a simple exponential fit to subtract the background TL signal, assuming constant linear heating throughout the read cycle. Our chips are heated using a linear increase in temperature from 50 to 300 ◦ C at a rate of 25 ◦ C s−1 for 10 s, followed by a steady temperature of 300 ◦ C for 6.67 s. This causes the background TL signals to level off as the chips reach thermal equilibrium. The rapid heating also means that the indicated temperature of the glow peak is not representative of the true glow peak temperature, and can shift several ◦ C because of variations in the readout (Stadtmann et al., 2006). OPG conducted a test of the background signal by repeatedly oven-annealing 67 cards at 80 ◦ C for 17 h and immediately reading them. Fig. 3 shows the results of the average signal from the fourth time this was done. It was found that the Boltzmann function fit the data well. Although Elements 3 and 4 were 2.5 times thinner than
0 0
50
100
150
200
Channel
Fig. 3. Residual signals on elements obtained on the fourth cycle of oven-annealing and reading.
Elements 1 and 2, the charge recorded in the high-temperature region is higher than for Elements 1 and 2. The Boltzmann function has the form shown in Eq. (1), where the constants A, B, and D are shown in Fig. 4. y=A−
B . 1 + exp[C(x − D)]
(1)
The constant C is proportional to the slope at the point (D, A − B/2) and inversely proportional to B. It can be shown that the first derivative of the Boltzmann function is a quadratic function with respect to y. This means that a Boltzmann function can be fit to a set of data once the plot of dy/dx vs. y for those points has been fitted with a parabola. Since a parabola can be fit to any set of data points using the quadratic least squares fitting method, a set of data points of dy/dx vs. y can be represented
624
W.J. Chase et al. / Radiation Measurements 43 (2008) 621 – 625
A
2nd Weibull
0.05
Residual
0.04 0.03 A-B/2
B
0.02 0.01 0.00
A-B
1
21
41
61
81
101 121 141 161 181
-0.01 D
-0.02 Fig. 4. The Boltzmann function with constants A, B, and D as shown. C depends on the slope at the point (D, A − B/2).
Boltzmann Background
Boltzmann Function
Raw
Fig. 6. Residual data and second Weibull fit to the residual data.
electronic noise, and are subtracted along with the Boltzmann background data.
Smooth
0.7
4.5. Fitting Weibull functions to glow peak data
0.6 0.5 0.4 0.3 0.2 0.1 0 1
21
41
61
81
101
121
141
161
181
Fig. 5. A typical Boltzmann fit. The function goes through and is tangential to the smoothed glow curve data at the end point. The Boltzmann background between the start and end points, and all data outside the start and end points, are subtracted.
by the following equation: jy = ay 2 + by + c. jx
(2)
Glow curve data points to the left of the start point and to the right of the end point are used to fit the parabola, but there are several restrictions on which points can be used. In particular, the points must create a parabola that is concave down. To do this, the chosen points to the left of the start point must be steadily increasing in height and slope, and the chosen points to the right of the end point must be increasing in height but decreasing in slope with a slope 0. After a number of algebraic manipulations, values for all four constants, A, B, C, and D, can be obtained. Fig. 5 shows a typical Boltzmann fit. Since the Boltzmann function does not necessarily represent the uneven portion of the curves often found to the right of the main glow peak, GCA only uses the Boltzmann function between the start and end points, xsp and xep . All data outside these points are considered
Pagonis et al. (2001) demonstrated that the Weibull function could successfully be used to model the TL intensity for peaks exhibiting first-order kinetics. In their study, the chips were heated linearly at a slow ramp rate. This allowed the Weibull function to have a constant shape factor for nearly every glow curve. Our heating cycle has a fast ramp rate and then flattens out. However, our observation is that the Weibull function is valid for nonlinear heating, as long as a new shape factor is calculated for each peak. The Weibull function fitted to the peak is of the form x − −1 x− W (x) = A , (3) exp − where A is the height adjustment factor, is the shape factor, is the scale factor, x is the channel number, and is the horizontal position factor. These values are referred to as the Weibull constants. Nine relationships were found relating the Weibull constants to several measurable quantities. Some of these include: the full width half maximum (FWHM); the right half-width (RHW); and yRHHM and yLHHM , where y signifies the first derivative of y with respect to x, LHHM is the left-hand half maximum, and RHHM is the right-hand half maximum. Each formula was calculated by comparing an ideal Weibull function to its measurable quantities. For example, can be expressed as a linear function of , which in turn can be expressed in terms of yRHHM and yLHHM . Once the four constants have been calculated, the Weibull function is known. GCA calculates the Weibull function and subtracts it from the background-subtracted data. The remaining data that remain after the Weibull are subtracted are called the residual, and are shown in Fig. 6. GCA fits a second Weibull function to the residual using the same method as was used to fit the first function. This ensures that any real dose due to the presence of peak 3 or 4 is accounted for. A flag GCA
W.J. Chase et al. / Radiation Measurements 43 (2008) 621 – 625 Boltzmann
Weibull
2nd Weibull
Sum
Raw
1.0
Start/End Points
Smooth
Sum of 2 Weibulls plus Boltzmann
0.9 0.8 0.7 1st Weibull
0.6
Background Noise
0.5
Smoothed Data
0.4 0.3 0.2
2nd Weibull
0.1 0.0 1
21
41
61
81
101
121
141
161
181
Fig. 7. Typical fitted glow curve showing comparison between the smoothed data and the sum of the Boltzmann function plus both Weibull functions.
uses during the glow curve analysis is a second Weibull (W2) result that is too large. If some contamination such as chalk or cement dust on the TLD card causes an erroneous glow signal, GCA will flag it for review. Once the Boltzmann and both Weibull functions have been calculated and plotted, they are summed. The summed curve is then used to compare the final fitted data to the original smoothed data, as shown in Fig. 7. The sum of the channel-bychannel differences between the summed curve and smoothed curve, divided by the sum of the smoothed curve, from channel 30 to the end point, is reported as the figure-of-merit (FOM), as shown in the following equation: xep |Sumi − Smoothi | FOM = i=30 . (4) xep j =30 Smoothj GCA flags all curves that have an FOM above a user-adjustable limit. The sum of the integrals of the two Weibull functions constitutes the element reading for a particular chip. Fig. 7 (without the text in the graph) is a typical graph from one of the element worksheets. The corresponding element graph on the main worksheet has data for the fit value (sum of two Weibull curves), smooth–Boltzmann (smoothed glow curve data minus the Boltzmann background between the start and end points), the CRI (total integral glow curve value in the calibration region), the Boltzmann background value, and the FOM. 5. Evaluation of program Bezaire (2003) conducted extensive testing of GCA. A subset of this testing will be described. The glow curve from one element of a low-dose (air kerma of 71 Gy, equivalent to exposure of 8 mR) card was used to test the despiking algorithm. Nineteen channels were chosen in places where it was felt the despiking algorithm would have the most difficulty or would have the most effect on dose, including some adjacent channels with spikes. At most about 3% of the dose was not removed.
625
Cards were exposed to various air kerma levels and the results compared to delivered air kermas ranging from 0 to 17.5 mGy (exposures ranging from 0 to 2000 mR). RCFs specific to the normal CRI method and to GCA were determined. The mean relative response for both methods was similar for all four elements and showed a small positive bias for low doses, perhaps because of some small inaccuracy in determining the radiation background signal to be subtracted. The coefficient of variation (CoV) in these results was similar for both methods for the 100 mg cm−2 chip, but for the 40 mg cm−2 chip the CoV for the GCA method was considerably smaller than for the CRI method. Other tests shows that GCA performs as well as or better than the CRI method in a wide variety of situations. The only area where GCA performs poorly compared to the CRI method is for very low dose cards (air kerma element result < 44 Gy, equivalent to exposure < 5 mR), where the poor quality of the glow curve data causes worse performance because of the extra data processing done by GCA. This is not a practical limitation for our personnel dosimetry system because all cards issued pick up sufficient signal from background exposure for processing by GCA. 6. Summary An Excel spreadsheet has been developed to do glow curve analysis of personnel dosimetry results. This has allowed review and correction of all glow curves, improving the quality of the results and reducing the total time required for complete processing of the official dosimetry results. Acknowledgments A.J. Pyatt did some of the initial work on inputting and plotting glow curve data. J. Aro performed testing which determined that certain methods, such as formulae using ROI values, would not be useful. B. Gaulke provided his spreadsheet using the Microsoft䉸 Excel姠 Excel Solver Add-in to analyze glow peaks using the Weibull function for our review. References Bezaire, M.D., 2003. Evaluation of the TLD glow curve analyzer computer program. Bruce Power Report B-REP-03146.14-00001, Rev 000, July 30. Chase, W.J., Hirning, C.R., 2008. Application of radiation physics in the design of the Harshaw 8828 beta-gamma TLD badge. Radiat. Meas., this issue, doi:10.1016/j.radmeas.2007.11.053. CNSC, 2006. Technical and quality assurance requirements for dosimetry services. Canadian Nuclear Safety Commission Regulatory Standard S-106 Revision 1, May. Pagonis, V., Mian, S.M., Kitis, G., 2001. Fit of first order thermoluminescence glow peaks using the Weibull distribution function. Radiat. Prot. Dosim. 93 (1), 11–17. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical recipes in FORTRAN. second ed.. Cambridge University Press, Cambridge. Stadtmann, H., Hranitzky, C., Brasik, N., 2006. Study of real time temperature profiles in routine TLD read out—influences of detector thickness and heating rate on glow curve shape. Radiat. Prot. Dosim. 119 (1–4), 310–313.
A glow curve analyzer (GCA) for routine analysis of personnel thermoluminescent dosemeter results W.J. Chase a,∗ , M.D. Bezaire b , F.P. Vanderzwet c , C.E. Taylor a a Health Physics Department, Ontario Power Generation, Whitby, Ont., Canada b Department of Applied Health Sciences, University of Waterloo, Waterloo, Ont., Canada c Bruce Power, P.O. Box 1540, Building B13, Tiverton, Ont., Canada N0G 2T0
Abstract A glow curve analyzer (GCA) spreadsheet has been developed using Microsoft 䉸 Excel姠 to perform glow curve analysis on thermoluminescent dosimeter (TLD) data from a personnel dosimetry system. The TLD data come from cards with four LiF:Mg,Ti chips that have been annealed and therefore have a simple glow peak structure. GCA removes spikes in the glow curve data, and then smoothes it. After select start and end points for the glow peak, it fits a Boltzmann function to represent the glow curve signal background under the glow peak. The Boltzmann function is subtracted and two Weibull curves are fit to the remaining net signal between the start and end points. The first Weibull curve is fit to peak 5, and the second one to any small remaining contribution from peaks 3 and 4 or from contaminants. The sum of the two Weibull curves is the glow curve signal result. GCA provides rapid review and correction of all glow curves, improving the quality of the results and reducing the time required for complete processing of official dose results. © 2008 Elsevier Ltd. All rights reserved. Keywords: TLD; Glow curve analysis; Weibull distribution; Boltzmann function
1. Introduction Ontario Power Generation (OPG) and Bruce Power (BP) use the Harshaw 8828 badge case and a four-element thermoluminescent dosimeter (TLD) card for personnel dosimetry. Elements 1 and 2 of the badge are 100 mg cm−2 LiF:Mg,Ti chips enriched in 7 Li (Harshaw TLD-700). Element 1 is filtered by 1010 mg cm−2 of low Z material, Element 2 is filtered by 1010 mg cm−2 of tin and Teflon, and they are used to measure photon signals. Elements 3 and 4 are 40 mg cm−2 TLD-700 chips with 17 mg cm−2 of filtration for Element 3 and 74 mg cm−2 for Element 4. They are used to measure the beta signals (see Chase and Hirning, 2008 for further information). The cards are annealed in an oven at 80 ◦ C for 17 h before use to minimize fading, by significantly reducing the size of the low temperature TL peaks. As a result peak 5 is the main peak, with a small contribution from peaks 3 and 4.
∗ Corresponding author. Tel.: +1 905 430 2215x3242; fax: +1 905 430 8583.
E-mail address: [email protected] (W.J. Chase). 1350-4487/$ - see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.radmeas.2008.02.020
For each card element, the TL glow curve signal is integrated over a defined region to obtain the calibration region integral (CRI), corrected for the chip sensitivity. The system is calibrated by exposing cards free-in-air to a 137 Cs source so that TL signals are converted to units of exposure in mR by means of element-specific reader calibration factors (RCFs). Element-specific values for the background radiation exposure are subtracted from the field card element TLD readings. The net values are processed with an algorithm to obtain personal dose equivalent results for Hp (10) and Hp (0.07). The original operating mode used the vendor-supplied software to review all glow curves for cards with doses that would be assignable (Hp (10) or Hp (0.07) 0.095 mSv (9.5 mrem)) and determine an adjustment as necessary. Most adjustments were required for the glow curves from Elements 3 and 4. These glow curves are about 2.5 times smaller in absolute terms than the glow curves from Elements 1 and 2, and thus usually noisier. These glow curves often have a significant high-temperature tail under and after the main glow peak, which required manual adjustment to correct the value for the glow peak area. Fig. 1 shows a typical glow curve for Element 3 or 4, with spikes
622
W.J. Chase et al. / Radiation Measurements 43 (2008) 621 – 625
3. Program description
1.6 Calibration Region 1.4
The main functionality and display of glow curves in GCA is contained in one worksheet of GCA, although many other worksheets are used for processing data or displaying the results for a single element. The main functions of GCA are either initiated on the main worksheet or selected from an additional top-level “GCA” menu item. GCA also includes the algorithm used to process the four net element results and convert them into the personal dose equivalents, Hp (10) and Hp (0.07). Many of the details of the GCA functions are specific to the OPG and BP TLD systems, but the GCA glow curve analysis should be applicable to TLD glow curves with a single main peak.
Adjustments
TL Signal
1.2 1.0 0.8 0.6 0.4 0.2 0.0 1
21
41
61
81
101 121 Channel
141
161
181
Fig. 1. Element 4 glow curve requiring correction to remove spikes (left-hand set of bars) and tail after glow peak (right-hand set of bars). The calibration region integral (CRI) of the TL signal can be manually adjusted by the signal in the adjustment regions.
before the glow peak and a high-temperature tail after the peak. The left-hand set of vertical bars shows how the effect of the spikes could be determined, while the right-hand set of bars shows how the tail after the glow peak up to the end of the calibration region could be estimated. However, using this manual method it is not readily possible to determine the value of the glow curve background signal under the glow peak. Our interest was primarily in accurately determining the glow peak area after subtracting any spikes and any residual high-temperature signal under and after the glow peak. After reviewing a number of options, it was determined that glow curve analysis would provide the most reliable results. BP in collaboration with OPG developed a computer program called glow curve analyzer (GCA). The program has been evaluated (Bezaire, 2003), and accepted by the Canadian Nuclear Safety Commission (CNSC) to meet the Canadian regulatory requirements for dosimetry service providers (CNSC, 2006). 2. Choice of glow curve analysis functions and processing system OPG and BP each read more than 10,000 TLD cards per processing period (quarterly or monthly), with many cards having low-level signals. The analysis method and functions chosen must therefore not rely on user input of starting values or require numerous iterations to obtain the final values. Pagonis et al. (2001) demonstrated that the Weibull function could successfully model the TL intensity for peaks exhibiting first-order kinetics. After considering various options, we decided to develop our own GCA using Microsoft 䉸 Excel姠 . In addition to its programming capability through Visual Basic for Applications (VBA), Excel can drag and drop points on a graph by the use of a “Data Marker”, allowing ready changes to certain parameters. We derived parameters directly from smoothed glow peak data that allowed us to readily determine suitable parameters for calculating the Weibull distribution, and did not require subsequent iterations to improve the fit.
4. Glow curve analysis 4.1. Despiking The value for each channel of glow curve data is compared to the average for the four nearest channels, and replaced by the average if it exceeds a certain value. Because there are cases where spikes two channels wide are not removed, the despiking routine is run again on the despiked data. 4.2. Smoothing The despiked data must be smoothed to properly determine certain “critical points” that are key to determining certain points on the glow curve, such as the start and end of the glow peak. The smoothing function is applied to the despiked data points compressed by adding 1 and then taking the natural logarithm. The smoothing routine works by fitting a running 15-point quadratic least squares fit, centered on the channel in question, to each of the 200 points. The parabola is fit using the method of least squares and the center of the 15-point parabola is then plotted as the new point. A second pass of the previously smoothed data is done to further smooth the curve and ensure no discontinuities remain. This is crucial, as the first and second derivatives of the glow curve data at every channel are required to determine the critical points. The GCA smoothing routine is identical to the Savitsky–Golay smoothing routine applied to points equally spaced in x (Press et al., 1992). After the second pass, the data are transformed back to a linear scale by reversing the original transformation. 4.3. Determination of start and end points Two “critical points” GCA determines are the start and end points for the glow peak. The calculation of these is shown in Fig. 2. The start point (xsp , ysp ) is the point where the first derivative has first reached a certain fractional value (the start point factor) of the absolute maximum of the first derivative. In effect, the start point is the location where the tangent of the glow peak has risen above a certain angle. The calculation of the end point (xep , yep ) starts with the location of the channel (channel a) of the absolute minimum
W.J. Chase et al. / Radiation Measurements 43 (2008) 621 – 625
623
b
0.04
0.7
0.03
0.6
0.02
0.5
0.01
0.4
0
End End Point Point
Start StartPoint Point
0.3
-0.01 0.2
a
-0.02
0.1
-0.03
0
-0.04 1
25
50
75
100
125
150
175
200
1
25
50
75
100
125
150
175
200
Fig. 2. A typical example of the first and second derivatives. The two vertical lines represent the locations of the start and end points. In the lower left figure and the inset figure, “a” and “b” point to the first derivative and the other line is the second derivative. In the right-hand figure, note the start and end point positions with respect to the smoothed data.
E1
60
E2
E3
E4
E1/E2 Boltzmann
E3/E4 Boltzmann
50 Charge (pC/16)
of the first derivative. Next, the first channel to the right of this (channel b) where the first derivative equals 0 is located. Lastly, the first channel to the right of channel b where the second derivative passes through zero from positive to negative is located. This is an inflection point, and one channel to its left is taken as the end point of the glow peak, (i.e. In other words, the end point is the channel just before the first channel to the right of the peak where there is an inflection point and a positive slope).
40 30 20 10
4.4. Boltzmann function for background signal under glow peak A key reason for creating GCA was to eliminate the residual signal under and after the glow peak. A main reason for this residual signal is the emission of black body radiation by the chips. As the chip temperature increases, the black body output of light increases exponentially. Many glow curve analysis and deconvolution packages use a simple exponential fit to subtract the background TL signal, assuming constant linear heating throughout the read cycle. Our chips are heated using a linear increase in temperature from 50 to 300 ◦ C at a rate of 25 ◦ C s−1 for 10 s, followed by a steady temperature of 300 ◦ C for 6.67 s. This causes the background TL signals to level off as the chips reach thermal equilibrium. The rapid heating also means that the indicated temperature of the glow peak is not representative of the true glow peak temperature, and can shift several ◦ C because of variations in the readout (Stadtmann et al., 2006). OPG conducted a test of the background signal by repeatedly oven-annealing 67 cards at 80 ◦ C for 17 h and immediately reading them. Fig. 3 shows the results of the average signal from the fourth time this was done. It was found that the Boltzmann function fit the data well. Although Elements 3 and 4 were 2.5 times thinner than
0 0
50
100
150
200
Channel
Fig. 3. Residual signals on elements obtained on the fourth cycle of oven-annealing and reading.
Elements 1 and 2, the charge recorded in the high-temperature region is higher than for Elements 1 and 2. The Boltzmann function has the form shown in Eq. (1), where the constants A, B, and D are shown in Fig. 4. y=A−
B . 1 + exp[C(x − D)]
(1)
The constant C is proportional to the slope at the point (D, A − B/2) and inversely proportional to B. It can be shown that the first derivative of the Boltzmann function is a quadratic function with respect to y. This means that a Boltzmann function can be fit to a set of data once the plot of dy/dx vs. y for those points has been fitted with a parabola. Since a parabola can be fit to any set of data points using the quadratic least squares fitting method, a set of data points of dy/dx vs. y can be represented
624
W.J. Chase et al. / Radiation Measurements 43 (2008) 621 – 625
A
2nd Weibull
0.05
Residual
0.04 0.03 A-B/2
B
0.02 0.01 0.00
A-B
1
21
41
61
81
101 121 141 161 181
-0.01 D
-0.02 Fig. 4. The Boltzmann function with constants A, B, and D as shown. C depends on the slope at the point (D, A − B/2).
Boltzmann Background
Boltzmann Function
Raw
Fig. 6. Residual data and second Weibull fit to the residual data.
electronic noise, and are subtracted along with the Boltzmann background data.
Smooth
0.7
4.5. Fitting Weibull functions to glow peak data
0.6 0.5 0.4 0.3 0.2 0.1 0 1
21
41
61
81
101
121
141
161
181
Fig. 5. A typical Boltzmann fit. The function goes through and is tangential to the smoothed glow curve data at the end point. The Boltzmann background between the start and end points, and all data outside the start and end points, are subtracted.
by the following equation: jy = ay 2 + by + c. jx
(2)
Glow curve data points to the left of the start point and to the right of the end point are used to fit the parabola, but there are several restrictions on which points can be used. In particular, the points must create a parabola that is concave down. To do this, the chosen points to the left of the start point must be steadily increasing in height and slope, and the chosen points to the right of the end point must be increasing in height but decreasing in slope with a slope 0. After a number of algebraic manipulations, values for all four constants, A, B, C, and D, can be obtained. Fig. 5 shows a typical Boltzmann fit. Since the Boltzmann function does not necessarily represent the uneven portion of the curves often found to the right of the main glow peak, GCA only uses the Boltzmann function between the start and end points, xsp and xep . All data outside these points are considered
Pagonis et al. (2001) demonstrated that the Weibull function could successfully be used to model the TL intensity for peaks exhibiting first-order kinetics. In their study, the chips were heated linearly at a slow ramp rate. This allowed the Weibull function to have a constant shape factor for nearly every glow curve. Our heating cycle has a fast ramp rate and then flattens out. However, our observation is that the Weibull function is valid for nonlinear heating, as long as a new shape factor is calculated for each peak. The Weibull function fitted to the peak is of the form x − −1 x− W (x) = A , (3) exp − where A is the height adjustment factor, is the shape factor, is the scale factor, x is the channel number, and is the horizontal position factor. These values are referred to as the Weibull constants. Nine relationships were found relating the Weibull constants to several measurable quantities. Some of these include: the full width half maximum (FWHM); the right half-width (RHW); and yRHHM and yLHHM , where y signifies the first derivative of y with respect to x, LHHM is the left-hand half maximum, and RHHM is the right-hand half maximum. Each formula was calculated by comparing an ideal Weibull function to its measurable quantities. For example, can be expressed as a linear function of , which in turn can be expressed in terms of yRHHM and yLHHM . Once the four constants have been calculated, the Weibull function is known. GCA calculates the Weibull function and subtracts it from the background-subtracted data. The remaining data that remain after the Weibull are subtracted are called the residual, and are shown in Fig. 6. GCA fits a second Weibull function to the residual using the same method as was used to fit the first function. This ensures that any real dose due to the presence of peak 3 or 4 is accounted for. A flag GCA
W.J. Chase et al. / Radiation Measurements 43 (2008) 621 – 625 Boltzmann
Weibull
2nd Weibull
Sum
Raw
1.0
Start/End Points
Smooth
Sum of 2 Weibulls plus Boltzmann
0.9 0.8 0.7 1st Weibull
0.6
Background Noise
0.5
Smoothed Data
0.4 0.3 0.2
2nd Weibull
0.1 0.0 1
21
41
61
81
101
121
141
161
181
Fig. 7. Typical fitted glow curve showing comparison between the smoothed data and the sum of the Boltzmann function plus both Weibull functions.
uses during the glow curve analysis is a second Weibull (W2) result that is too large. If some contamination such as chalk or cement dust on the TLD card causes an erroneous glow signal, GCA will flag it for review. Once the Boltzmann and both Weibull functions have been calculated and plotted, they are summed. The summed curve is then used to compare the final fitted data to the original smoothed data, as shown in Fig. 7. The sum of the channel-bychannel differences between the summed curve and smoothed curve, divided by the sum of the smoothed curve, from channel 30 to the end point, is reported as the figure-of-merit (FOM), as shown in the following equation: xep |Sumi − Smoothi | FOM = i=30 . (4) xep j =30 Smoothj GCA flags all curves that have an FOM above a user-adjustable limit. The sum of the integrals of the two Weibull functions constitutes the element reading for a particular chip. Fig. 7 (without the text in the graph) is a typical graph from one of the element worksheets. The corresponding element graph on the main worksheet has data for the fit value (sum of two Weibull curves), smooth–Boltzmann (smoothed glow curve data minus the Boltzmann background between the start and end points), the CRI (total integral glow curve value in the calibration region), the Boltzmann background value, and the FOM. 5. Evaluation of program Bezaire (2003) conducted extensive testing of GCA. A subset of this testing will be described. The glow curve from one element of a low-dose (air kerma of 71 Gy, equivalent to exposure of 8 mR) card was used to test the despiking algorithm. Nineteen channels were chosen in places where it was felt the despiking algorithm would have the most difficulty or would have the most effect on dose, including some adjacent channels with spikes. At most about 3% of the dose was not removed.
625
Cards were exposed to various air kerma levels and the results compared to delivered air kermas ranging from 0 to 17.5 mGy (exposures ranging from 0 to 2000 mR). RCFs specific to the normal CRI method and to GCA were determined. The mean relative response for both methods was similar for all four elements and showed a small positive bias for low doses, perhaps because of some small inaccuracy in determining the radiation background signal to be subtracted. The coefficient of variation (CoV) in these results was similar for both methods for the 100 mg cm−2 chip, but for the 40 mg cm−2 chip the CoV for the GCA method was considerably smaller than for the CRI method. Other tests shows that GCA performs as well as or better than the CRI method in a wide variety of situations. The only area where GCA performs poorly compared to the CRI method is for very low dose cards (air kerma element result < 44 Gy, equivalent to exposure < 5 mR), where the poor quality of the glow curve data causes worse performance because of the extra data processing done by GCA. This is not a practical limitation for our personnel dosimetry system because all cards issued pick up sufficient signal from background exposure for processing by GCA. 6. Summary An Excel spreadsheet has been developed to do glow curve analysis of personnel dosimetry results. This has allowed review and correction of all glow curves, improving the quality of the results and reducing the total time required for complete processing of the official dosimetry results. Acknowledgments A.J. Pyatt did some of the initial work on inputting and plotting glow curve data. J. Aro performed testing which determined that certain methods, such as formulae using ROI values, would not be useful. B. Gaulke provided his spreadsheet using the Microsoft䉸 Excel姠 Excel Solver Add-in to analyze glow peaks using the Weibull function for our review. References Bezaire, M.D., 2003. Evaluation of the TLD glow curve analyzer computer program. Bruce Power Report B-REP-03146.14-00001, Rev 000, July 30. Chase, W.J., Hirning, C.R., 2008. Application of radiation physics in the design of the Harshaw 8828 beta-gamma TLD badge. Radiat. Meas., this issue, doi:10.1016/j.radmeas.2007.11.053. CNSC, 2006. Technical and quality assurance requirements for dosimetry services. Canadian Nuclear Safety Commission Regulatory Standard S-106 Revision 1, May. Pagonis, V., Mian, S.M., Kitis, G., 2001. Fit of first order thermoluminescence glow peaks using the Weibull distribution function. Radiat. Prot. Dosim. 93 (1), 11–17. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical recipes in FORTRAN. second ed.. Cambridge University Press, Cambridge. Stadtmann, H., Hranitzky, C., Brasik, N., 2006. Study of real time temperature profiles in routine TLD read out—influences of detector thickness and heating rate on glow curve shape. Radiat. Prot. Dosim. 119 (1–4), 310–313.