ARTICLE IN PRESS
Nuclear Instruments and Methods in Physics Research A 526 (2004) 439–446
A novel iterative energy calibration method for composite germanium detectors N.S. Pattabiraman, S.N. Chintalapudi, S.S. Ghugre* Inter University Consortium For DAE Facilities, Calcutta Center, Sector III, Block LB, Plot 8, Bidhan Nagar, Kolkata 700 098, India Received 24 December 2001; received in revised form 2 February 2004; accepted 24 February 2004
Abstract An automatic method for energy calibration of the observed experimental spectrum has been developed. The method presented is based on an iterative algorithm and presents an efficient way to perform energy calibrations after establishing the weights of the calibration data. An application of this novel technique for data acquired using composite detectors in an in-beam g-ray spectroscopy experiment is presented. r 2004 Elsevier B.V. All rights reserved. PACS: 23.90.+w; 02.60.Pn; 02.60.Ed; 29.85.+c Keywords: (HI, xng) coincidence data; Clover detectors; Energy calibration
1. Introduction Advances in nuclear spectroscopy studies have been feasible due to the innovations in the detector and the associated pulse processing techniques. The enhancement in the photo-peak efficiency which results in a higher detection limit and improved sensitivity of modern gamma detector arrays are due to the use of composite detectors. Thus, these modern state of the art g-ray spectroscopy arrays provide a growing amount of data making it imperative to automate most of the pre-sorting and sorting procedures. The use of composite detectors presents the user with a major challenge, primarily due to the large number of detector channels, which have to be calibrated individually *Corresponding author. E-mail address:
[email protected] (S.S. Ghugre).
prior to sorting of the data into the conventional gN histograms [1,2]. An automated calibration procedure would minimize to a great extent the human intervention and book keeping, resulting in an efficient pre-sorting of the data. However, due to the large number of detectors used, both germanium as well as the ancillary detectors, the pulse processing and the digitizing techniques employed are different. Hence, the use of single mathematical function for energy calibration (conversion from the observed channel number to the corresponding energy), will have certain limitations. A new iterative method for energy calibration of such inbeam g-ray spectroscopy data is presented in this paper. This iterative algorithm has been incorporated into a PC-based g analysis software ‘‘IUCSORT’’ [3]. An application highlighting the power of this procedure to a nonlinear analog-to-digital converter is also discussed.
0168-9002/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.nima.2004.02.033
ARTICLE IN PRESS 440
N.S. Pattabiraman et al. / Nuclear Instruments and Methods in Physics Research A 526 (2004) 439–446
2. Energy calibration and gain matching The acquired coincidence list mode data is sorted into the conventional gN histogram. Prior to this sorting, for simplification it is desirable that the data obtained from each detector should have a constant energy dispersion. This process is termed as gain matching and could be defined as, 1 gain matched channel x0 ¼ epc Eg ; where Eg is the energy corresponding to an ADC channel x: The typical value of epc (energy per channel) for the gain matched data is 0:5 keV/channel for 4096 4096 Eg –Eg matrix and 2 keV/channel for 1024 1024 1024 Eg –Eg –Eg cube. The gain matching is achieved only after performing accurate energy calibration for all detectors. The relation between the ADC channel (x) and the gamma ray energy Eg is parameterized using a polynomial equation of order n; Eg ¼ a10 þ a11 xb11 þ a12 xb12 þ ? þ a1n xb1n
ð1Þ
where ai are the coefficients, and bi are natural numbers. The order of the polynomial is decided depending on the availability of data points. These coefficients have to be computed and are difficult to derive manually. One of the major problems encountered during the derivation of the calibration curve (calibration parameters) is the nonlinearity in the observed channel-energy relation. This nonlinearity originates from the electronics used and cannot be identical for all the detectors. Hence, it is not advisable to use an uniform form of the above equation for all the detectors. The calibration should be derived after establishing the weights for the calibration points. An iterative procedure to obtain the appropriate form for this relationship is described below. Standard radioactive g calibration sources are used to obtain the energy-channel calibration curve [4–6]. The distribution of data points plays an important role in the convergence of the solutions for obtaining the values of the coefficients in Eq. (1). The beam-off radioactivity data provides unique calibration data points at high energies even beyond 2 MeV: Hence, it has been used to obtain the energy calibration curve. Traditionally, a polynomial equation of order 2 was used to describe the relation between the Eg
and the recorded pffiffiffi channel x: Higher-order polynomial and or x are essential to account for the observed nonlinearity in the electronics [4]. Palacz and co-workers [4] suggested the following relation for the energy calibration curve, pffiffiffi Eg ¼ a0 þ a1 x þ a2 x2 þ a3 x3 þ a4 x: ð2Þ The cubic term in the above equation results in a better fit atphigher energies (above 1:5 MeV), ffiffiffi whereas the x is useful at lower energies (below 200 keV). When the energy calibration was initially performed using a quadratic polynomial, it is observed that the difference ðEactual Efit Þ is typically of the order of 200 eV; except at low energy, where it is about 400 eV; which is not desirable. The calibration curve was then obtained using the relation given by Palacz et al. [4]. This narrowed down the energy difference between the fitted and actual energy to less than 100 eV: Liu et al. [7], have employed the following equation to derive calibration curve between energy range 0.3 and 0:6 MeV Eg ¼ a0 þ a1 x þ a2 x2 þ a3 x3 þ eðxÞ
ð3Þ
where eðxÞ is a function that takes care of the residuals of remaining higher order terms. With larger number of data points, this equation narrows down the energy difference between the fitted and actual value to about 3 eV: The success of this procedure, within the small energy range, prompted us to explore of the possibility of dividing the entire energy range into two smaller ranges and to perform an independent calibration for each region. A similar procedure has been reported by Hu and co-workers in Ref. [5]. The calibrations could be performed using either the same functional form, or a different form for the two regions, i.e., ðEg Þ ¼ dðx; limitÞða10 þ a11 xb11 þ a12 xb12 þ a13 xb13 Þ þd% ðx; limitÞða20 þ a21 xb21 þ a22 xb22 þ a23 xb23 Þ ð4Þ where dðx; limitÞ ¼ 1; d% ðx; limitÞ ¼ 0 when xolimit; dðx; limitÞ ¼ 0; d% ðx; limitÞ ¼ 1 when xXlimit
ARTICLE IN PRESS N.S. Pattabiraman et al. / Nuclear Instruments and Methods in Physics Research A 526 (2004) 439–446
N X ððEactuali Efiti Þ2 wli wsi Þ
ð5Þ
i¼1
where N is the total number of data points used in calibration and n df1 ; n df2 correspond to number of degrees of freedom available in first and second data set and Eactual and Efit denote the actual and the fitted energy of the gamma ray. Weight ws is equal to variance of number of counts in the photo peak. Weight wl takes care of energy gap Eg between the data points and it is defined as x x 2 1 0 wli ¼ if i ¼ 1 2 x x 2 iþ1 i1 if i > 1 and io N wli ¼ 2 x x 2 i i1 wli ¼ if i ¼ N: 2 The value of wl is interchanged between nearby points, in such a way that the peak which had a greater value of ws should have a greater value of wl ; provided the energy difference between these points did not exceed a certain DE value. For the present case DE corresponds to 400 keV: This ensures that a strong and well separated peak from its nearest group of g-rays, is fitted more precisely. Thus appropriate weight during calibration is given to the data point which is relatively well separated from it’s neighbors. The method described above was used to obtain the energy calibration curve for the data acquired
wsm owsn > wso > wsp owsq owsr owss wlm > wln > wlo > wlp > wlq owlr owls :
500 200 100 50
ws
(a) n
s r q
20 o 10 m p 5 2 106
(b)
105 m 104
wl
1 w2 ¼ P2 ð j¼1 n dfj ÞðxN x1 Þ
using a 8 Clover detector array [8,9]. The experiment was performed at the 14 UD BARC-TIFR Pelletron facility at Mumbai, India [10]. The experiment was aimed to observe the high spin states in AB90 region, using the 31 P þ63 Cu reaction at an incident beam energy of 120 MeV: The calibration points (B18 points) are obtained from a combination of radioactive sources as well as the beam-off radioactivity data. Fig. 1 illustrates, the arrangements of wl in accordance with ws : In the figure plots (a) and (b) represent ws and wl of individual g-rays which are used for the calibration. In these plots few sample peaks are marked with letters m, n, o, p and q. The comparison of the weights ws and wl for these peaks are
Before adjustment
nop qr s
103 102 101 106
(c)
105
n
104
wl
and limit corresponds to the channel number which separates the data points into two subgroups, for independent calibration. The parameters a11 ; a12 ; a13 ; a21 ; a22 ; a23 are coefficients, and a10 ; a20 are constant terms. The parameters bij are natural numbers, except the parameter b11 ; where the b11 can be either 1 or 0.5. The assignment of b11 ¼ 0:5; usually results in better fit to the lower energy calibration data points [4]. The sum of squares of difference between the actual and fitted value (a figure of merit for goodness of the fit) is calculated as
441
After adjustment
m o rs pq
103 102 101 200
600
1000
1400
1800
2200
E (keV) Fig. 1. The weights ws (plot (a)) and wl (plot (b)) of individual g-rays (53 keVoEg o2320 keV). Plot (c) illustrates, the arrangement of wl in plot (b) in accordance with ws in plot (a). Refer text for details.
ARTICLE IN PRESS 442
N.S. Pattabiraman et al. / Nuclear Instruments and Methods in Physics Research A 526 (2004) 439–446
Once the wl is arranged in accordance with ws ; the above comparison is transformed to wlm owln > wlo > wlp owlq owlr owls : This arrangement is illustrated in Fig. 1(c). The values of bij are set to b11 ¼ 1; b12 ¼ 2; b13 ¼ 3; b21 ¼ 1; b22 ¼ 2; b23 ¼ 3 resulting in Eg ¼ dðx; limitÞða10 þ a11 x1 þ a12 x2 þ a13 x3 Þ þ d% ðx; limitÞða20 þ a21 x1 þ a22 x2 þ a23 x3 Þ: ð6Þ Similar to the procedure described in SAMPO80 [11], the data point after the limit and a data point before the limit has been included in the first and second sets, respectively. This ensures a smooth interpolation in the region between the two data sets. These 18 data points are initially sub-divided into 2 groups with channel number corresponding
to the 5th data point as the boundary between the two sets. Thus, there are 6 points in the first set and 15 points in the second. A linear calibration is performed for both the sets. The value of w2 as given by Eq. (5) is obtained. The w2 value along with the calibration parameters which include both the coefficients and the channel limit between the two regions are stored. The number of data points is then increased by 1 in the first set. This point is accordingly reduced from the second set. Linear calibration is then repeated, for both these sets. If the obtained w2 is less than the value from the earlier iteration, the new w2 and the calibration parameters are stored. The procedure is repeated till the first data set encompasses 15 points and the remaining 6 points belong to the second set. In each case, parameters obtained from the earlier iteration are replaced by the new
Fig. 2. A representative pffiffifficalibration fit using the present procedure. The calibration curve in first data set is given by the equation, Eg ¼ ð0:866 þ 0:170 x þ 0:397xÞ and the equation, Eg ¼ ð1:915 þ 0:399x 3:419 107 x2 þ 3:421 1011 x3 Þ represents the calibration curve in the second data set. These equations result in the best fit to the data and are derived from the iterative process. The open and filled diamonds represents the deviation of the fitted value from the actual value by using the first and second equations, respectively. The points ‘f’ and ‘s’ indicate the last calibration point in the first set and first calibration point in the second set, respectively. The limit is indicated by the mark ‘l’. The dotted and dashed lines are to guide an eye and indicate the calibration fit using the equations for first and second data sets, respectively.
ARTICLE IN PRESS N.S. Pattabiraman et al. / Nuclear Instruments and Methods in Physics Research A 526 (2004) 439–446
parameters only if these result in a better value for the obtained w2 : The above procedure is then repeated for a linear-quadratic combination where the first data set is subjected to a linear and the second set to a quadratic parameterization. This procedure is then extended to the cyclic linear-quadratic-cubic combinations, i.e., linear-cubic, quadratic-linear, quadratic-cubic, cubic-linear..., where there are 9 such combinations in all, for the present calibration data. Then the values of bij are set to b11 ¼ 0:5; b12 ¼ 1; b13 ¼ 2; b21 ¼ 1; b22 ¼ 2; b23 ¼ 3 resulting in pffiffiffi Eg ¼ dðx; limitÞða10 þ a11 x þ a12 x1 þ a13 x2 Þ þ d% ðx; limitÞða20 þ a21 x1 þ a22 x2 þ a23 x3 Þ: ð7Þ
443
to these the calibration parameters from linear, quadratic and cubic fits to the entire data are also recorded.
3. Applications This procedure, is referred to as an iterative multi-functional procedure. A representative energy calibration using the above procedure is shown in the Fig. 2. The plotted errors include contribution from statistical and fitting errors. Fig. 3 illustrates the corresponding values of w2 for all the 32 detectors obtained from the above iterative process. The w2 values obtained using Eq. (2) are also indicated for comparison. As seen from the figure this iterative procedure resulted in a better fit to the known energy values. For the 32 detectors it was found that the different combinations of the above generalized equation resulted in a better fit to the calibration data. This emphasizes the fact that use of a fixed pffiffi equation like ( ðxÞ + cubic, or quadratic) to derive the energy-channel relation (calibration curve) may be inappropriate in some cases. In
The iterations pare ffiffiffi performed using the above equation. The x term is used only for the first data set, since this functional form results in a better fit to the data at low energies, as mentioned earlier. At the end of this iterative procedure, the calibration parameters aij ; bij and limiting channel number which resulted in the least w2 value among the above 18 combinations are stored. In addition
- Cubic + Channel 45
- Itterative
2
35
25
15
5 5
15
25
Detector Number 2 Fig. 3. The results of the present iterative calibration process, illustrated pffiffiffi by the resultant w (refer text for more details), for all the 32 detectors. The values obtained by using the conventional cubic + x calibration equation are indicated for comparison.
ARTICLE IN PRESS 444
N.S. Pattabiraman et al. / Nuclear Instruments and Methods in Physics Research A 526 (2004) 439–446
Fig. 4. Plot of the observed nonlinearity as a function of the g-ray energy. The nonlinearity is obtained by subtracting out the contribution from the linear response. The points ‘f’ and ‘s’ indicate the last calibration point in the first set and first calibration point in the second set, respectively. The limit is indicated by the mark ‘l’. The dash+dotted and dashed lines indicates the calibration fit pffiffiffi using the first and second set of equations, respectively. Similar results for the quadratic and cubic + x calibration curve are illustrated for comparison. The errors include the fitting and statistical errors.
general, one expects that a equation with a larger number of free parameters would result in a better fit to the data. However, one would also have to include the appropriate weight for the points based on their nearest neighbors and sub-divide the data into groups to account for any localized nonlinearity in the data. As seen from the present method, the inclusion of these indicates that a higher order equation need not always result in a better fit to the data. The distinct advantage of the present method is that the functional form for the parameterization of the calibration curve does not require an individual treatment. The user only has to provide the calibration points (channel and energy values). The parameters for the 2 regions are automatically obtained from these data without any manual intervention.
To test the energy dependence of the observed nonlinearity in the calibration curve, a linear fit was performed for the entire data (shown in Fig. 2). Following this the calibration was performed using the above procedure. Then the contribution from the linear fit was subtracted. This is plotted in Fig. 4. As seen from the figure the nonlinearity has a distinctly different functional form in the 2 regions, which were automatically identified by the present process. For comparison similar results for the calibration curve obtained pffiffiffifrom the conventional quadratic and cubic + x fitting are included. As seen from the figure,pcurve obtained using the conventional ffiffiffi cubic + x equation does not fit the observed data reasonably, in a region, which was identified by the present method as the limit between the two
ARTICLE IN PRESS N.S. Pattabiraman et al. / Nuclear Instruments and Methods in Physics Research A 526 (2004) 439–446
(a)
(b)
445
between the predicted and actual energy value is reasonable. This method has also been tested when analogto-digital converter used had a nonlinear response especially at higher channels. The calibration data was analyzed using the present method. Fig. 5(b) illustrates the deviations between the predicted and actual energy values. For comparison we have also included the results using the established cubic + pffiffiffi x fit. As seen from the figure the agreement between the predicted and the actual energy values using the present procedure is remarkable. This establishes the reliability of the procedure which does not require any additional input from the user except the conventional calibration data points.
4. Conclusions Fig. 5. The difference between the actual and predicted energy values for the points which are not included in the calibration pffiffiffi fit. The open circles represent the conventional cubic + x calibration and the filled circles represent the present method. Plot (a) illustrates the comparison for a representative sample case. Plot (b) illustrates the comparison for a special case wherein the ADC has a nonlinear performance.
data sets, which require a different functional form for the calibration curve as mentioned earlier. In other regions, the two methods give almost identical results. This validates the assumption in the present method to parameterize the entire energy range into 2 different regions. The energy calibration is supposed to predict reliably the Eg ; for a given channel x: This could be ascertained by performing the above procedure and obtaining the Eg for a point which is not included in the parametrization. The above procedure was done by excluding one data point at a time from the fitting procedure. The predicted energy for this point is obtained from the derived calibration curve. This is repeated for all the data points. The deviations form the actual energy values using the above procedure is illustrated in Fig. 5(a). For comparison we have pffiffiffiincluded the results obtained from the cubic + x fit. As seen from the figure the agreement
An automatic method has been developed for the energy calibration using a multi-functional equation, which also incorporates the weight for the data points during the calibration procedure. The true power of the presented method is the ability to deal with a nonlinearity which is manifested only within small part of the spectrum. The method has been successfully applied to inbeam data acquired using a multi-Clover array.
Acknowledgements We are indebted to Dr. S.K. Basu whose comments and suggestions were helpful throughout this work. We would like to express our sincere gratitude to Dr. A.K. Sinha, Dr. R.K. Bhowmik, Dr. P.K. Joshi, Mr. R.P. Singh for their suggestions and constant encouragement during this work. The help provided by Ms. Sudatta Ray, Mr. Angha Chakraborty and Mr. Krishichayan during the data analysis is gratefully acknowledged.
References [1] S.L. Tabor, Nucl. Instr. and Meth. A265 (1988) 495.
ARTICLE IN PRESS 446
N.S. Pattabiraman et al. / Nuclear Instruments and Methods in Physics Research A 526 (2004) 439–446
. [2] C.W. Beausang, D. Pr!evost, M.H. Bergstrom, G. de France, B. Haas, J.C. Lisle, Ch. Theisen, J. Tim!ar, P.J. Twin, J.N. Wilson, Nucl. Instr. and Meth. A364 (1995) 560. [3] N.S. Pattabiraman, S.N. Chintalapudi, S.S. Ghugre, Nucl. Instr. and Meth. A, accompanying paper. [4] M. Palacz, J. Cederk.all, M. Lipoglav$sek, L.-O. Norlin, J. Nyberg, J. Persson, Nucl. Instr. and Meth. A383 (1996) 473. [5] Z. Hu, T. Glasmacher, W.F. Mueller, I. Wiedenhover, Nucl. Instr. and Meth. A482 (2002) 715. [6] R.G. Helmer, C. van der Leun, Nucl. Instr. and Meth. A450 (2000) 35.
[7] Zhigang Liu, Liangqiang Peng, Long Wei, Tianbao Chang, Nucl. Instr. and Meth. A432 (1999) 122. [8] P.K. Joshi, H.C. Jain, A.S. Medhi, S. Chattopadhyay, S. Bhattacharya, A. Goswami, Nucl. Instr. and Meth. A399 (1997) 51. [9] G. Duch#ene, F.A. Beck, P.J. Twin, G. de France, D. Curien, L. Han, C.W. Beausang, M.A. Bentley, P.J. Nolan, J. Simpson, Nucl. Instr. and Meth. A432 (1999) 90. [10] S. Lakshmi, H.C. Jain, P.K. Joshi, Amita, P. Agarwal, A.K. Jain, S.S. Malik, Phys. Rev. C 66 (2002) 043103(R). [11] M.J. Koskelo, P.A. Aarnio, J.T. Routi, Comput. Phys. Commun. 24 (1981) 11.