X射线峰形描述技术及分支比约束在解谱中的应用

王思广 韦江波 张 欢

（北京大学物理学院 核物理及核技术国家重点实验室 北京 100871）

X射线峰形描述技术及分支比约束在解谱中的应用

王思广 韦江波 张 欢

（北京大学物理学院 核物理及核技术国家重点实验室 北京 100871）

尽管多道能谱中X射线与γ射线重叠区内峰下的计数常具有能够提供重要信息的可能，但有时因这种重叠能区难于解谱而被摒弃。能量依赖的探测效率会使X射线峰出现左右不对称。探测器对X及γ射线响应形状的准确描述，及将已知信息在拟合中作为约束项的应用是精确分解出多道能谱中重叠区各成份计数的保障。充分利用RooFit软件包，建议了一种X射线峰形描述方式：利用洛伦兹函数描述其本征形状并与探测效率相乘，之后与分辨函数卷积以描述探测后的形状。测试显示放入与拟合出的峰下计数在统计误差范围内一致。所建议的方法及提供的代码可作为X及γ射线重峰区的解谱参考。

解谱，拟合，多道能谱

利用多道能谱测量系统进行放射性核素的X及γ射线能谱的测量具有长期的历史，当前探测器探头有HPGe、NaI(Tl)、CdTe/CdZnTe等多种选择，能谱解谱算法及解谱分析通用软件也很多。但要精确给出能谱分析结果，响应函数的准确描述及已知信息的充分利用系根本保障。

X及γ射线的复合能区的解谱在实际工作中经常用到。前期工作中曾提出利用逐渐逼近的方法修正因不同能量所对应的探测效率的不同而导致的X射线形状的变形[1]。得益于对RooFit软件包[2]功能的理解，近期研究实现了直接利用该软件包提供的函数进行描述该形变的方法及描绘代码。鉴于RooFit软件包本身是开源代码，对外完全免费，利用我们的方法，可以很容易进行X及γ射线的复合能区的解谱工作，为国内外同行提供一些参考。

1 X射线形状

X射线的本征形状（探测器探测前）可用洛伦兹函数形式描述：

式中，x0为峰位中心位置；Γ对应半高处的全宽度。

对于每一次衰变，发出的X射线能量不定，其数值的分布可用该函数的形式随机取样模拟。

以高纯锗探测器为例，对一能量为E0的X射线，要被探测器探测到，需历经如下过程：(1) 不被放射源本身吸收；(2) 不被源与探测器之间的空气或其他介质吸收或散射掉；(3) 发出的角度正好在探测器的灵敏区对放射源所张的空间几何角度之内；(4) 顺利打入到探测器的灵敏区而不被探测器前窗（如有）及死层吸收；(5) 将全部能量转换为电子空穴对而被全部收集，转成对应信号幅度。所有这些因素都会对总的探测效率产生影响。并且(1)、(2)和(4)的反应截面与射线的能量有关，从而导致探测效率ε(x)为能量的函数。假设在X射线对应的能量区间探测效率随能量而增大，则X射线的左侧将较右侧压低得厉害，而右侧将较左侧高一些。考虑探测效率的影响后，探测器分辨无限好时的X射线的形状可用L(x)ε(x)表示。

因探测的有限分辨率，探测系统所给出的能量将是以入射能量E0为中心在一定范围内变化后的数值，作为例子，这里用高斯函数描述探测系统给出的能谱（实际能谱处理可根据具体情况引入合适的修正项）。故最终探测到的X射线的形状可用[L(x)ε(x)]⊗G(x)进行表述。其中⊗表示卷积，通常可用快速傅里叶变换处理实际数据[3]。

2 γ射线形状

γ射线所发的射线本征形状原则上也应该用洛伦兹函数描述，但因其能级的寿命较长，与探测系统的分辨相比，其宽度太小，故可简化成单一能量Eγ。同样因为有限的探测器分辨，该能量会被展宽，通常用高斯函数描述该展宽，故γ射线的形状在这里简化为高斯描述。对于实际的探测器，因电荷收集的不完全，有可能在低能端形成拖尾，有的仪器会因堆积效应，在高能端形成高能尾巴。实际使用时可以用较为复杂的函数及方法描述γ射线形状或仪器分辨[4−8]。

3 约束拟合信息

对于重叠较严重的复合峰区，要准确给出各成份的贡献，除提供正确的各射线峰的形状外，还需要进一步提供约束信息。例如：峰的位置、X射线的本征宽度、假定探测效率已知的情况下的峰下计数之比。前两者只需要能量刻度信息及查阅并放入对应的核数据即可，后一约束条件的使用除分支比信息外还需巧妙处理：对于来自同一核素的两条γ射线，其能量分别用E1和E2表示，所查出的射线的分支比分别为Br1和Br2，如果待拟合的效率用ε(x)表示，则两射线所对应的能谱上的贡献N1与N2之比为：

在进行解谱时，假定N2为待拟合的未知数，则N1可利用式(2)表达出来参加拟合，减少未知量。

4 测试及结果

4.1 数据的产生

数据的产生使用Root软件进行编程[9]。该软件系开源软件，其所提供的TRandom3伪随机数产生子[10]的周期长度达106000。测试工作中所产生的能谱有两条γ射线和一条X射线，其中X射线位于两条γ射线之间。假定探测效率ε(x)及能量分辨的σ(x)值随能量而增大（或减小）。

对来自γ射线的两信号，按照高斯函数随机抽样产生。信号的分辨率σ值按照既定的随能量变化函数σ(x)计算给出，峰下计数按照分支比及对应的探测效率根据式(2)定出的比例给出。

X射线信号的产生按照以下步骤进行：首先用式(1)随机抽样出能量E0，然后根据ε(E0)进行抽样决定是否被探测到（如果在[0,1]区间抽取的随机数R＜ε(E0)则保留，否则重新抽取E0），对于探测到的事件，用该能量及既定的分辨率随能量的变化，计算出σ(E0)值，之后利用以σ(E0)为宽度、以0为平均值的高斯函数随机抽取出ΔE，E0+ΔE即为被探测器测量的X射线的能量。将所产生的γ射线及X射线的数据合在一起即为MC产生的能谱数据。之后将拟合该数据，以测试拟合出的量是否在误差允许范围内以及与放入的各成份的数目相符。

4.2 数据的拟合

数据的拟合是用Root的可选软件包RooFit实现。我们的测试结果如图1所示，图1(a)与(b)系同一拟合结果，不同点在于图1(b)用对数坐标。两条γ射线分别位于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();

}

拟合能谱中的三种成份共用同一个效率修正函数和同一表征分辨随能量变化的σ(x)函数。产生和拟合时两条γ射线来自同一核素且其分支比为1:2的信息已知并在拟合过程中运用约束技术。并且假定X射线的宽度也为已知量（当宽度未知时，仅需在定义WidthX时给定允许变化的范围即可）。

图1是在两条γ射线距离X射线峰位比较远，峰与峰之间能明显区分开时的拟合结果。拟合中峰的位置系待拟合量。如果将两条γ射线峰向X峰靠拢并重新产生数据，如图2所示。图2(a)与(b)系同一拟合结果，不同点在于图2(b)用对数坐标。两条γ射线分别位于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射线峰下的原始计数为60005，拟合值为(6.002±0.030)×104；86.0 keV和93.0keV处的γ射线峰的原始计数分别为43970和143868，拟合值分别为(4.383±0.023)×104及(1.4400±0.0041)×105，拟合时峰的位置系浮动待拟合参数。对于图2，X射线峰下的原始计数为60039，拟合的结果为(6.041±0.043)×104；88.5 keV处的γ射线下原始计数为54062，拟合值为(5.405± 0.025)×104；92.0 keV处的γ射线的原始计数为135948，拟合值为(1.3560± 0.0040)×105。故两种情况下分别对应的图1、图2所拟合的结果与能谱产生时放入的各成份的计数在误差允许范围内一致。

5 结语

从X射线的产生、传播、被探测到的过程完整机制出发，用Root软件所带的RooFit软件包精细描述各因素对其形状的影响，并提供一完整拟合复合峰区的解谱代码。该代码还示出利用来自同一核素不同能量的γ射线的分支比的比值在解谱中的约束技术的应用方法。测试显示，采用我们的X射线描述方法可以有效实现因探测效率的不同而导致的峰形左右不对称的描述，测试结果显示输入、输出值在误差允许范围内相符。该方法及完整测试代码可为同类研究提供借鉴。

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全身计数器γ能谱全能峰函数的探讨[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