王思廣 韋江波 張 歡
(北京大學(xué)物理學(xué)院 核物理及核技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室 北京 100871)
X射線峰形描述技術(shù)及分支比約束在解譜中的應(yīng)用
王思廣 韋江波 張 歡
(北京大學(xué)物理學(xué)院 核物理及核技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室 北京 100871)
盡管多道能譜中X射線與γ射線重疊區(qū)內(nèi)峰下的計(jì)數(shù)常具有能夠提供重要信息的可能,但有時(shí)因這種重疊能區(qū)難于解譜而被摒棄。能量依賴的探測(cè)效率會(huì)使X射線峰出現(xiàn)左右不對(duì)稱。探測(cè)器對(duì)X及γ射線響應(yīng)形狀的準(zhǔn)確描述,及將已知信息在擬合中作為約束項(xiàng)的應(yīng)用是精確分解出多道能譜中重疊區(qū)各成份計(jì)數(shù)的保障。充分利用RooFit軟件包,建議了一種X射線峰形描述方式:利用洛倫茲函數(shù)描述其本征形狀并與探測(cè)效率相乘,之后與分辨函數(shù)卷積以描述探測(cè)后的形狀。測(cè)試顯示放入與擬合出的峰下計(jì)數(shù)在統(tǒng)計(jì)誤差范圍內(nèi)一致。所建議的方法及提供的代碼可作為X及γ射線重峰區(qū)的解譜參考。
解譜,擬合,多道能譜
利用多道能譜測(cè)量系統(tǒng)進(jìn)行放射性核素的X及γ射線能譜的測(cè)量具有長(zhǎng)期的歷史,當(dāng)前探測(cè)器探頭有HPGe、NaI(Tl)、CdTe/CdZnTe等多種選擇,能譜解譜算法及解譜分析通用軟件也很多。但要精確給出能譜分析結(jié)果,響應(yīng)函數(shù)的準(zhǔn)確描述及已知信息的充分利用系根本保障。
X及γ射線的復(fù)合能區(qū)的解譜在實(shí)際工作中經(jīng)常用到。前期工作中曾提出利用逐漸逼近的方法修正因不同能量所對(duì)應(yīng)的探測(cè)效率的不同而導(dǎo)致的X射線形狀的變形[1]。得益于對(duì)RooFit軟件包[2]功能的理解,近期研究實(shí)現(xiàn)了直接利用該軟件包提供的函數(shù)進(jìn)行描述該形變的方法及描繪代碼。鑒于RooFit軟件包本身是開(kāi)源代碼,對(duì)外完全免費(fèi),利用我們的方法,可以很容易進(jìn)行X及γ射線的復(fù)合能區(qū)的解譜工作,為國(guó)內(nèi)外同行提供一些參考。
X射線的本征形狀(探測(cè)器探測(cè)前)可用洛倫茲函數(shù)形式描述:
式中,x0為峰位中心位置;Γ對(duì)應(yīng)半高處的全寬度。
對(duì)于每一次衰變,發(fā)出的X射線能量不定,其數(shù)值的分布可用該函數(shù)的形式隨機(jī)取樣模擬。
以高純鍺探測(cè)器為例,對(duì)一能量為E0的X射線,要被探測(cè)器探測(cè)到,需歷經(jīng)如下過(guò)程:(1) 不被放射源本身吸收;(2) 不被源與探測(cè)器之間的空氣或其他介質(zhì)吸收或散射掉;(3) 發(fā)出的角度正好在探測(cè)器的靈敏區(qū)對(duì)放射源所張的空間幾何角度之內(nèi);(4) 順利打入到探測(cè)器的靈敏區(qū)而不被探測(cè)器前窗(如有)及死層吸收;(5) 將全部能量轉(zhuǎn)換為電子空穴對(duì)而被全部收集,轉(zhuǎn)成對(duì)應(yīng)信號(hào)幅度。所有這些因素都會(huì)對(duì)總的探測(cè)效率產(chǎn)生影響。并且(1)、(2)和(4)的反應(yīng)截面與射線的能量有關(guān),從而導(dǎo)致探測(cè)效率ε(x)為能量的函數(shù)。假設(shè)在X射線對(duì)應(yīng)的能量區(qū)間探測(cè)效率隨能量而增大,則X射線的左側(cè)將較右側(cè)壓低得厲害,而右側(cè)將較左側(cè)高一些。考慮探測(cè)效率的影響后,探測(cè)器分辨無(wú)限好時(shí)的X射線的形狀可用L(x)ε(x)表示。
因探測(cè)的有限分辨率,探測(cè)系統(tǒng)所給出的能量將是以入射能量E0為中心在一定范圍內(nèi)變化后的數(shù)值,作為例子,這里用高斯函數(shù)描述探測(cè)系統(tǒng)給出的能譜(實(shí)際能譜處理可根據(jù)具體情況引入合適的修正項(xiàng))。故最終探測(cè)到的X射線的形狀可用[L(x)ε(x)]?G(x)進(jìn)行表述。其中?表示卷積,通??捎每焖俑道锶~變換處理實(shí)際數(shù)據(jù)[3]。
γ射線所發(fā)的射線本征形狀原則上也應(yīng)該用洛倫茲函數(shù)描述,但因其能級(jí)的壽命較長(zhǎng),與探測(cè)系統(tǒng)的分辨相比,其寬度太小,故可簡(jiǎn)化成單一能量Eγ。同樣因?yàn)橛邢薜奶綔y(cè)器分辨,該能量會(huì)被展寬,通常用高斯函數(shù)描述該展寬,故γ射線的形狀在這里簡(jiǎn)化為高斯描述。對(duì)于實(shí)際的探測(cè)器,因電荷收集的不完全,有可能在低能端形成拖尾,有的儀器會(huì)因堆積效應(yīng),在高能端形成高能尾巴。實(shí)際使用時(shí)可以用較為復(fù)雜的函數(shù)及方法描述γ射線形狀或儀器分辨[4?8]。
對(duì)于重疊較嚴(yán)重的復(fù)合峰區(qū),要準(zhǔn)確給出各成份的貢獻(xiàn),除提供正確的各射線峰的形狀外,還需要進(jìn)一步提供約束信息。例如:峰的位置、X射線的本征寬度、假定探測(cè)效率已知的情況下的峰下計(jì)數(shù)之比。前兩者只需要能量刻度信息及查閱并放入對(duì)應(yīng)的核數(shù)據(jù)即可,后一約束條件的使用除分支比信息外還需巧妙處理:對(duì)于來(lái)自同一核素的兩條γ射線,其能量分別用E1和E2表示,所查出的射線的分支比分別為Br1和Br2,如果待擬合的效率用ε(x)表示,則兩射線所對(duì)應(yīng)的能譜上的貢獻(xiàn)N1與N2之比為:
在進(jìn)行解譜時(shí),假定N2為待擬合的未知數(shù),則N1可利用式(2)表達(dá)出來(lái)參加擬合,減少未知量。
4.1 數(shù)據(jù)的產(chǎn)生
數(shù)據(jù)的產(chǎn)生使用Root軟件進(jìn)行編程[9]。該軟件系開(kāi)源軟件,其所提供的TRandom3偽隨機(jī)數(shù)產(chǎn)生子[10]的周期長(zhǎng)度達(dá)106000。測(cè)試工作中所產(chǎn)生的能譜有兩條γ射線和一條X射線,其中X射線位于兩條γ射線之間。假定探測(cè)效率ε(x)及能量分辨的σ(x)值隨能量而增大(或減?。?/p>
對(duì)來(lái)自γ射線的兩信號(hào),按照高斯函數(shù)隨機(jī)抽樣產(chǎn)生。信號(hào)的分辨率σ值按照既定的隨能量變化函數(shù)σ(x)計(jì)算給出,峰下計(jì)數(shù)按照分支比及對(duì)應(yīng)的探測(cè)效率根據(jù)式(2)定出的比例給出。
X射線信號(hào)的產(chǎn)生按照以下步驟進(jìn)行:首先用式(1)隨機(jī)抽樣出能量E0,然后根據(jù)ε(E0)進(jìn)行抽樣決定是否被探測(cè)到(如果在[0,1]區(qū)間抽取的隨機(jī)數(shù)R<ε(E0)則保留,否則重新抽取E0),對(duì)于探測(cè)到的事件,用該能量及既定的分辨率隨能量的變化,計(jì)算出σ(E0)值,之后利用以σ(E0)為寬度、以0為平均值的高斯函數(shù)隨機(jī)抽取出ΔE,E0+ΔE即為被探測(cè)器測(cè)量的X射線的能量。將所產(chǎn)生的γ射線及X射線的數(shù)據(jù)合在一起即為MC產(chǎn)生的能譜數(shù)據(jù)。之后將擬合該數(shù)據(jù),以測(cè)試擬合出的量是否在誤差允許范圍內(nèi)以及與放入的各成份的數(shù)目相符。
4.2 數(shù)據(jù)的擬合
數(shù)據(jù)的擬合是用Root的可選軟件包RooFit實(shí)現(xiàn)。我們的測(cè)試結(jié)果如圖1所示,圖1(a)與(b)系同一擬合結(jié)果,不同點(diǎn)在于圖1(b)用對(duì)數(shù)坐標(biāo)。兩條γ射線分別位于86.0 keV(γ1)和93.0 keV(γ2);X射線峰位于90.0 keV。
圖1 擬合測(cè)試結(jié)果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();
}
擬合能譜中的三種成份共用同一個(gè)效率修正函數(shù)和同一表征分辨隨能量變化的σ(x)函數(shù)。產(chǎn)生和擬合時(shí)兩條γ射線來(lái)自同一核素且其分支比為1:2的信息已知并在擬合過(guò)程中運(yùn)用約束技術(shù)。并且假定X射線的寬度也為已知量(當(dāng)寬度未知時(shí),僅需在定義WidthX時(shí)給定允許變化的范圍即可)。
圖1是在兩條γ射線距離X射線峰位比較遠(yuǎn),峰與峰之間能明顯區(qū)分開(kāi)時(shí)的擬合結(jié)果。擬合中峰的位置系待擬合量。如果將兩條γ射線峰向X峰靠攏并重新產(chǎn)生數(shù)據(jù),如圖2所示。圖2(a)與(b)系同一擬合結(jié)果,不同點(diǎn)在于圖2(b)用對(duì)數(shù)坐標(biāo)。兩條γ射線分別位于88.5 keV(γ1)和92.0 keV(γ2);X射線峰位于90.0 keV。這時(shí)候如果僅靠直接觀察,并不易感知在90.0 keV處有X射線峰的存在。即便隨意放上X射線峰形也很難保障被準(zhǔn)確擬合出。但如果充分利用已知各成份的能量及刻度信息,在擬合的時(shí)候固定下各成份的峰位,擬合結(jié)果將基本準(zhǔn)確。
圖2 擬合測(cè)試結(jié)果Fig.2 Test results of fit.
測(cè)試結(jié)果:圖1的X射線峰下的原始計(jì)數(shù)為60005,擬合值為(6.002±0.030)×104;86.0 keV和93.0keV處的γ射線峰的原始計(jì)數(shù)分別為43970和143868,擬合值分別為(4.383±0.023)×104及(1.4400±0.0041)×105,擬合時(shí)峰的位置系浮動(dòng)待擬合參數(shù)。對(duì)于圖2,X射線峰下的原始計(jì)數(shù)為60039,擬合的結(jié)果為(6.041±0.043)×104;88.5 keV處的γ射線下原始計(jì)數(shù)為54062,擬合值為(5.405± 0.025)×104;92.0 keV處的γ射線的原始計(jì)數(shù)為135948,擬合值為(1.3560± 0.0040)×105。故兩種情況下分別對(duì)應(yīng)的圖1、圖2所擬合的結(jié)果與能譜產(chǎn)生時(shí)放入的各成份的計(jì)數(shù)在誤差允許范圍內(nèi)一致。
從X射線的產(chǎn)生、傳播、被探測(cè)到的過(guò)程完整機(jī)制出發(fā),用Root軟件所帶的RooFit軟件包精細(xì)描述各因素對(duì)其形狀的影響,并提供一完整擬合復(fù)合峰區(qū)的解譜代碼。該代碼還示出利用來(lái)自同一核素不同能量的γ射線的分支比的比值在解譜中的約束技術(shù)的應(yīng)用方法。測(cè)試顯示,采用我們的X射線描述方法可以有效實(shí)現(xiàn)因探測(cè)效率的不同而導(dǎo)致的峰形左右不對(duì)稱的描述,測(cè)試結(jié)果顯示輸入、輸出值在誤差允許范圍內(nèi)相符。該方法及完整測(cè)試代碼可為同類研究提供借鑒。
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 王思廣, 冒亞軍, 唐培家, 等. 多道γ能譜的精細(xì)分析[J].核技術(shù), 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ì)數(shù)器γ能譜全能峰函數(shù)的探討[J]. 核技術(shù), 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
項(xiàng)目(No.10775006、No.10979010)資助
王思廣,男,1971年出生,2006年于北京大學(xué)博士后研究工作出站并留校任教,研究領(lǐng)域?yàn)榱W游锢砼c原子核物理
2014-11-17,
2014-12-09