王思廣 韋江波 張 歡

（北京大學物理學院 核物理及核技術國家重點實驗室 北京 100871）

X射線峰形描述技術及分支比約束在解譜中的應用

王思廣 韋江波 張 歡

（北京大學物理學院 核物理及核技術國家重點實驗室 北京 100871）

盡管多道能譜中X射線與γ射線重疊區(qū)內峰下的計數(shù)常具有能夠提供重要信息的可能，但有時因這種重疊能區(qū)難于解譜而被摒棄。能量依賴的探測效率會使X射線峰出現(xiàn)左右不對稱。探測器對X及γ射線響應形狀的準確描述，及將已知信息在擬合中作為約束項的應用是精確分解出多道能譜中重疊區(qū)各成份計數(shù)的保障。充分利用RooFit軟件包，建議了一種X射線峰形描述方式：利用洛倫茲函數(shù)描述其本征形狀并與探測效率相乘，之后與分辨函數(shù)卷積以描述探測后的形狀。測試顯示放入與擬合出的峰下計數(shù)在統(tǒng)計誤差范圍內一致。所建議的方法及提供的代碼可作為X及γ射線重峰區(qū)的解譜參考。

解譜，擬合，多道能譜

利用多道能譜測量系統(tǒng)進行放射性核素的X及γ射線能譜的測量具有長期的歷史，當前探測器探頭有HPGe、NaI(Tl)、CdTe/CdZnTe等多種選擇，能譜解譜算法及解譜分析通用軟件也很多。但要精確給出能譜分析結果，響應函數(shù)的準確描述及已知信息的充分利用系根本保障。

X及γ射線的復合能區(qū)的解譜在實際工作中經(jīng)常用到。前期工作中曾提出利用逐漸逼近的方法修正因不同能量所對應的探測效率的不同而導致的X射線形狀的變形[1]。得益于對RooFit軟件包[2]功能的理解，近期研究實現(xiàn)了直接利用該軟件包提供的函數(shù)進行描述該形變的方法及描繪代碼。鑒于RooFit軟件包本身是開源代碼，對外完全免費，利用我們的方法，可以很容易進行X及γ射線的復合能區(qū)的解譜工作，為國內外同行提供一些參考。

1 X射線形狀

X射線的本征形狀（探測器探測前）可用洛倫茲函數(shù)形式描述：

式中，x0為峰位中心位置；Γ對應半高處的全寬度。

對于每一次衰變，發(fā)出的X射線能量不定，其數(shù)值的分布可用該函數(shù)的形式隨機取樣模擬。

以高純鍺探測器為例，對一能量為E0的X射線，要被探測器探測到，需歷經(jīng)如下過程：(1) 不被放射源本身吸收；(2) 不被源與探測器之間的空氣或其他介質吸收或散射掉；(3) 發(fā)出的角度正好在探測器的靈敏區(qū)對放射源所張的空間幾何角度之內；(4) 順利打入到探測器的靈敏區(qū)而不被探測器前窗（如有）及死層吸收；(5) 將全部能量轉換為電子空穴對而被全部收集，轉成對應信號幅度。所有這些因素都會對總的探測效率產生影響。并且(1)、(2)和(4)的反應截面與射線的能量有關，從而導致探測效率ε(x)為能量的函數(shù)。假設在X射線對應的能量區(qū)間探測效率隨能量而增大，則X射線的左側將較右側壓低得厲害，而右側將較左側高一些。考慮探測效率的影響后，探測器分辨無限好時的X射線的形狀可用L(x)ε(x)表示。

因探測的有限分辨率，探測系統(tǒng)所給出的能量將是以入射能量E0為中心在一定范圍內變化后的數(shù)值，作為例子，這里用高斯函數(shù)描述探測系統(tǒng)給出的能譜（實際能譜處理可根據(jù)具體情況引入合適的修正項）。故最終探測到的X射線的形狀可用[L(x)ε(x)]?G(x)進行表述。其中?表示卷積，通?？捎每焖俑道锶~變換處理實際數(shù)據(jù)[3]。

2 γ射線形狀

γ射線所發(fā)的射線本征形狀原則上也應該用洛倫茲函數(shù)描述，但因其能級的壽命較長，與探測系統(tǒng)的分辨相比，其寬度太小，故可簡化成單一能量Eγ。同樣因為有限的探測器分辨，該能量會被展寬，通常用高斯函數(shù)描述該展寬，故γ射線的形狀在這里簡化為高斯描述。對于實際的探測器，因電荷收集的不完全，有可能在低能端形成拖尾，有的儀器會因堆積效應，在高能端形成高能尾巴。實際使用時可以用較為復雜的函數(shù)及方法描述γ射線形狀或儀器分辨[4?8]。

3 約束擬合信息

對于重疊較嚴重的復合峰區(qū)，要準確給出各成份的貢獻，除提供正確的各射線峰的形狀外，還需要進一步提供約束信息。例如：峰的位置、X射線的本征寬度、假定探測效率已知的情況下的峰下計數(shù)之比。前兩者只需要能量刻度信息及查閱并放入對應的核數(shù)據(jù)即可，后一約束條件的使用除分支比信息外還需巧妙處理：對于來自同一核素的兩條γ射線，其能量分別用E1和E2表示，所查出的射線的分支比分別為Br1和Br2，如果待擬合的效率用ε(x)表示，則兩射線所對應的能譜上的貢獻N1與N2之比為：

在進行解譜時，假定N2為待擬合的未知數(shù)，則N1可利用式(2)表達出來參加擬合，減少未知量。

4 測試及結果

4.1 數(shù)據(jù)的產生

數(shù)據(jù)的產生使用Root軟件進行編程[9]。該軟件系開源軟件，其所提供的TRandom3偽隨機數(shù)產生子[10]的周期長度達106000。測試工作中所產生的能譜有兩條γ射線和一條X射線，其中X射線位于兩條γ射線之間。假定探測效率ε(x)及能量分辨的σ(x)值隨能量而增大（或減?。?/p>

對來自γ射線的兩信號，按照高斯函數(shù)隨機抽樣產生。信號的分辨率σ值按照既定的隨能量變化函數(shù)σ(x)計算給出，峰下計數(shù)按照分支比及對應的探測效率根據(jù)式(2)定出的比例給出。

X射線信號的產生按照以下步驟進行：首先用式(1)隨機抽樣出能量E0，然后根據(jù)ε(E0)進行抽樣決定是否被探測到（如果在[0,1]區(qū)間抽取的隨機數(shù)R＜ε(E0)則保留，否則重新抽取E0），對于探測到的事件，用該能量及既定的分辨率隨能量的變化，計算出σ(E0)值，之后利用以σ(E0)為寬度、以0為平均值的高斯函數(shù)隨機抽取出ΔE，E0+ΔE即為被探測器測量的X射線的能量。將所產生的γ射線及X射線的數(shù)據(jù)合在一起即為MC產生的能譜數(shù)據(jù)。之后將擬合該數(shù)據(jù)，以測試擬合出的量是否在誤差允許范圍內以及與放入的各成份的數(shù)目相符。

4.2 數(shù)據(jù)的擬合

數(shù)據(jù)的擬合是用Root的可選軟件包RooFit實現(xiàn)。我們的測試結果如圖1所示，圖1(a)與(b)系同一擬合結果，不同點在于圖1(b)用對數(shù)坐標。兩條γ射線分別位于86.0 keV(γ1)和93.0 keV(γ2)；X射線峰位于90.0 keV。

圖1 擬合測試結果Fig.1 Test results of fit.

為便于同行參考，本文給出完整的代碼：

using namespace RooFit;

void SetSgStyle();

void txt(Double_t x0, Double_t y0,

TString str, Int_t iColor=kBlack,Double_t fsize = 0.06);

//main function

void drawFig1(){

//For plot style, can be commented out

SetSgStyle();

Int_t Nbin = 200;

Double_t xmin = 80;

Double_t xmax = 100;

TH1D *hSpec = new TH1D("hSpec","",Nbin,xmin,xmax);

//Initial Peak Position

Double_t fMeanG1 = 86.0;

Double_t fMeanG2 = 93.0;

Double_t fMeanX = 90.0;

//read the data of spectrum

FILE *pf = fopen("Spectrum.txt","rt");

Int_t n = 1;

while(n==1){

Float_teng;

n = fscanf(pf,"%f ",&eng);

if(n==1){

if(xmin＜=eng&&eng＜xmax) hSpec-＞Fill(eng);

}

}

fclose(pf);

RooRealVar x("x","x",xmin,xmax);

//X-ray peak

//peak position

RooRealVar MeanX("MeanX", "Mean of X-rayPeak", fMeanX,xmin,xmax);

//peak width.fix to 1-keV assumed we know

RooRealVar WidthX("WidthX", "Width of X-Ray Peak", 1.);

RooBreitWigner peakX0("peakX0", "X-ray Peak BeforeDetected",x,MeanX,WidthX);

//Efficiency

RooRealVar a("a","1st coefficient of efficiency",0.2, 0.01, 0.4);

RooRealVar b("b","2nd coefficient of efficiency",0.8, 0.0, 2.0);

RooFormulaVar effFun("effFun","a+b*(x-80.)/ (100.-80.)", RooArgList(a,b,x))；

//Resolution

RooRealVar mg("mg","mg",0);

RooRealVar c1("c1","1st coefficient of sigma",0.1, 0.0001, 0.5);

RooRealVar c2("c2","2nd coefficient of sigma",0.2, 0.0001, 1.);

RooFormulaVar

sg("sg","c1+c2*sqrt(MeanX-80.0)",RooArgList(c1,c2,MeanX)) ;

RooGaussian R("R","resolution",x,mg,sg);

//X-ray Peak * effFun

RooEffProd peakX0eff("peakX0eff","peakX0 * effFun", peakX0,effFun);

//Set #bins to be used for FFT sampling to 10000

x.setBins(10000,"cache");

//peakX0eff convoluted with Resolution

RooFFTConvPdf peakX("peakX","(peakX0*Eff) (X) gauss", x, peakX0eff, R);

//1st gamma peak

RooRealVar Mean1("Mean1","Mean of 1st GammaPeak", fMeanG1,fMeanG1-5,fMeanG1+5);

//Share the same resolution function with X-ray by

//using same coefficients c1 and c2

RooFormulaVar

Sigma1("Sigma1","(c1+c2*sqrt(Mean1-80.))",

RooArgList(c1,c2,Mean1));

RooGaussian peak1("peak1","1st Peak", x, Mean1, Sigma1);

//2nd gamma peak

RooRealVar Mean2("Mean2","Mean of 2nd gamma Peak",fMeanG2, fMeanG2-5, fMeanG2+5);

RooFormulaVar

Sigma2("Sigma2","(c1+c2*sqrt(Mean2-80.0))",RooArgList(c1, c2,Mean2));

RooGaussian peak2("peak2","2nd Peak", x, Mean2, Sigma2);

Double_tfN1,fNX,fN2;

//initial value of counts under each peak

fN1 = hSpec-＞GetEntries()*.33;

fNX = fN1; fN2 = fN1;

RooRealVar NX("NX","Counts under X Peak", fNX, 0, 3*fNX);

RooRealVar N2("N2","Counts under 2ndGamma Peak", fN2, 0, 3*fN2);

//Br1:Br2=1:2 N1:N2=(Br1*Eff1): (Br2*Eff2)

RooRealVar Br1vsBr2("Br1vsBr2","ratio between Br1 and Br2",0.5);

//Using the branch ratio as constrain RooFormulaVar

N1("N1","(a+b*(Mean1-80.)/(100.-80.))*

Br1vsBr2/(a+b*(Mean2-80.)/(100.-80.))*N2",

RooArgList(a,b,Mean1,Mean2,N2,Br1vsBr2));

//total functions

RooAddPdf model("model","model",RooArgList(peak1, peakX,peak2),RooArgList(N1,NX,N2));

RooDataHistd h("dh","data histogram to hold spec.", x,Import(*hSpec));

//fit the function to the data

RooFitResult *frlt = model.fitTo(dh,Save());

//draw the fit spectrum

TCanvas *cPlot = new TCanvas("cPlot","Fig.1",550,800); cPlot-＞Divide(1,2);

cPlot-＞cd(1);

RooPlot* frame = x.frame();

dh.plotOn(frame);

model.plotOn(frame, Components(peak1),

LineStyle(2),LineColor(kViolet));

model.plotOn(frame, Components(peakX),

LineStyle(3),LineColor(kRed));

model.plotOn(frame, Components(peak2),

LineStyle(4),LineColor(kBlue));

model.plotOn(frame,LineStyle(kSolid),

LineColor(kBlack));

frame-＞SetTitle(";Energy [keV]; Counts/(0.1keV)");

frame-＞SetMaximum(8000); //to fix Y axis in HighEdge

frame-＞SetMinimum(1); //to fix Y axis in LowEdge

frame-＞DrawClone();

//print the fitted results

Int_t nParsToFit = (frlt-＞floatParsFinal()).getSize();

//reduced chi-squared = chi2/ndof

Double_t chi2_red = frame-＞chiSquare(nParsToFit);

txt(0.22,0.88,Form("#chi^{2}/ndof=%.2f",chi2_red)); txt(0.22,0.82,Form("N_{X}=%.0f#pm%.0f",

NX.getVal(),NX.getError()),kRed);

txt(0.22,0.76, Form("N_{#gamma2}=%.0f#pm%.0f",

N2.getVal(),N2.getError()),kBlue);

Double_t N1Error = sqrt(N1.getVal()/N2.getVal()) *N2.getError();

txt(0.22,0.70, Form("N_{#gamma1}=%.0f#pm%.0f",

N1.getVal(),N1Error),kViolet);

txt(0.85,0.82,"a)");

cPlot-＞cd(2);

gPad-＞SetLogy(1);

frame-＞SetMaximum(10000); //to make sure displayed out 10^4

frame-＞SetMinimum(1); //logY from 1

frame-＞Draw();

txt(0.85,0.82,"b)");

cPlot-＞cd();

cPlot-＞Modified();

cPlot-＞Update();

cPlot-＞SaveAs("Fig1.eps");

}

//for the special style of the plot, can be commented out void SetSgStyle(){

gStyle-＞SetCanvasBorderMode(0);

gStyle-＞SetCanvasBorderSize(0);

gStyle-＞SetCanvasColor(10);

//Format for axes

gStyle-＞SetLabelFont(22,"xyz");

gStyle-＞SetLabelSize(0.06,"xyz");

gStyle-＞SetLabelOffset(0.01,"xyz");

gStyle-＞SetNdivisions(510,"xyz");

gStyle-＞SetTitleFont(22,"xyz");

gStyle-＞SetTitleColor(1,"xyz");

gStyle-＞SetTitleSize(0.06,"xyz");

gStyle-＞SetTitleOffset(0.91);

gStyle-＞SetTitleYOffset(1.1);

//No pad borders

gStyle-＞SetPadBorderMode(0);

gStyle-＞SetPadBorderSize(0);

gStyle-＞SetPadColor(10);

//Margins for labels etc.

gStyle-＞SetPadLeftMargin(0.15);

gStyle-＞SetPadBottomMargin(0.15);

gStyle-＞SetPadRightMargin(0.05);

gStyle-＞SetPadTopMargin(0.06);

//No error bars in x direction

gStyle-＞SetErrorX(0);

//Format legend

gStyle-＞SetLegendBorderSize(0);

gStyle-＞SetLegendFont(22);

gStyle-＞SetFillStyle(0);

}

//for label text on the plot

void txt(Double_t x0, Double_t y0,

TString str, Int_t iColor,Double_t fsize){

TLatex *ltx = new TLatex();

ltx-＞SetNDC(kTRUE);

ltx-＞SetTextColor(iColor);

ltx -＞SetTextFont(22);

ltx-＞SetTextSize(fsize);

ltx-＞DrawLatex(x0,y0,str.Data());

gPad-＞Modified();

gPad-＞Update();

}

擬合能譜中的三種成份共用同一個效率修正函數(shù)和同一表征分辨隨能量變化的σ(x)函數(shù)。產生和擬合時兩條γ射線來自同一核素且其分支比為1:2的信息已知并在擬合過程中運用約束技術。并且假定X射線的寬度也為已知量（當寬度未知時，僅需在定義WidthX時給定允許變化的范圍即可）。

圖1是在兩條γ射線距離X射線峰位比較遠，峰與峰之間能明顯區(qū)分開時的擬合結果。擬合中峰的位置系待擬合量。如果將兩條γ射線峰向X峰靠攏并重新產生數(shù)據(jù)，如圖2所示。圖2(a)與(b)系同一擬合結果，不同點在于圖2(b)用對數(shù)坐標。兩條γ射線分別位于88.5 keV(γ1)和92.0 keV(γ2)；X射線峰位于90.0 keV。這時候如果僅靠直接觀察，并不易感知在90.0 keV處有X射線峰的存在。即便隨意放上X射線峰形也很難保障被準確擬合出。但如果充分利用已知各成份的能量及刻度信息，在擬合的時候固定下各成份的峰位，擬合結果將基本準確。

圖2 擬合測試結果Fig.2 Test results of fit.

測試結果：圖1的X射線峰下的原始計數(shù)為60005，擬合值為(6.002±0.030)×104；86.0 keV和93.0keV處的γ射線峰的原始計數(shù)分別為43970和143868，擬合值分別為(4.383±0.023)×104及(1.4400±0.0041)×105，擬合時峰的位置系浮動待擬合參數(shù)。對于圖2，X射線峰下的原始計數(shù)為60039，擬合的結果為(6.041±0.043)×104；88.5 keV處的γ射線下原始計數(shù)為54062，擬合值為(5.405± 0.025)×104；92.0 keV處的γ射線的原始計數(shù)為135948，擬合值為(1.3560± 0.0040)×105。故兩種情況下分別對應的圖1、圖2所擬合的結果與能譜產生時放入的各成份的計數(shù)在誤差允許范圍內一致。

5 結語

從X射線的產生、傳播、被探測到的過程完整機制出發(fā)，用Root軟件所帶的RooFit軟件包精細描述各因素對其形狀的影響，并提供一完整擬合復合峰區(qū)的解譜代碼。該代碼還示出利用來自同一核素不同能量的γ射線的分支比的比值在解譜中的約束技術的應用方法。測試顯示，采用我們的X射線描述方法可以有效實現(xiàn)因探測效率的不同而導致的峰形左右不對稱的描述，測試結果顯示輸入、輸出值在誤差允許范圍內相符。該方法及完整測試代碼可為同類研究提供借鑒。

1 Wang S G, Mao Y J, Tang P J. Method for unfolding multiple regions of X- and γ-ray spectra with a detection efficiency constraint avoiding inflecting the peak shapes for correction results[J]. Nuclear Instruments and Methods in Physical Research Section A, 2009, 600: 445?452

2 Verkerke W, Kirkby D. RooFituser's manual, v2.91[CP]. :ftp://root.cern.ch/root/doc/RooFit_Users_Ma nual_2.91-33.pdf, 2008-10-14

3 Frigo M, Johnson S G. FFTW version 3.3.3[CP]. http://fftw.org/fftw3.pdf, 2012-11

4 Wang S G, Mao Y J, Tang P J. A method for interpolating asymmetric peak shapes in multiple γ-ray spectra[J]. Chinese Physics C, 2009, 33(3): 383?386

5 Siegert H, Janssen H. Precise determination of gamma-ray peak areas[J]. Nuclear Instruments and Methods in Physical Research Section A, 1990, 286(3): 415?420

6 Ebeid M R, Kaid M A, Ali M G S. Intrinsic shape analysis of X-ray diffraction using Stokes deconvolution[J]. Canadian Journal of Physics, 2014, 92(4): 316?320

7 王思廣, 冒亞軍, 唐培家, 等. 多道γ能譜的精細分析[J].核技術, 2006, 29(7): 495?498

WANG Siguang, MAO Yajun, TANG Peijia, et al. A precise analysis of multi-channel γ-ray spectrum[J]. Nuclear Techniques, 2006, 29(7): 495?498

8 齊榮, 毛永, 陳熙萌. 快速NaI全身計數(shù)器γ能譜全能峰函數(shù)的探討[J]. 核技術, 2008, 31(5): 330?334

QI Rong, MAO Yong, CHEN Ximeng. Full peak fitting function to γ-ray spectra from a quick NaI whole body counter[J]. Nuclear Techniques, 2008, 31(5): 330?334

9 Brun R, Rademakers F. ROOT user's guide[CP]. http://root.cern.ch/root/html534/guides/users-guide/ROO TUsersGuideA4.pdf , 2013-5

10 Matsumoto M M, Nishimura T. Mersenne twister: a 623-dimensionally equidistributed uniform pseudorandom number generator[J]. ACM Transactions on Modeling and Computer Simulations, 1998, 8(1): 3?20

CLC TL99

X-ray peak shape description method and branch ratio constraint applied in spectrum unfolding

WANG Siguang WEI Jiangbo ZHANG Huan
(School of Physics and State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, China)

Background: Although counts under peaks in overlap region of multi-channel spectrum from X- and γ-rays can provide useful information, they are usually discarded because of much difficulty in being unfolded. The line-shape of X-ray peak could be asymmetric in left and right side caused by an energy-dependent detector efficiency. Purpose: To accurately fit out the shapes for the counts under the peaks, X- and γ-ray shapes should be described accurately, and all known information should be used as constraints for reliable results. Methods: With RooFit package, a method to describe the X-ray shape is proposed: take its intrinsic shape of X-ray described by a Lorentz function and multiplied by a detector efficiency function, then convoluted with a resolution function. Results: The results tested with the codes included in this paper show: the events number of input and output agrees within their statistical errors checked even with very heavily overlapped peaks unfolded with branch ratio constraint method. Conclusion: The method and the code introduced in this paper are feasible, and they can be used as an example to unfold overlap regions in X- and γ-ray multi-channel spectrum.

Spectrum unfolding, Fit, Multi-channel spectrum

TL99

10.11889/j.0253-3219.2015.hjs.38.020502

項目(No.10775006、No.10979010)資助

王思廣，男，1971年出生，2006年于北京大學博士后研究工作出站并留校任教，研究領域為粒子物理與原子核物理

2014-11-17，

2014-12-09