針對中微子探測的多變量分析本底抑制方法

郭一, 張余煉,2, 魏志勇,2

(1.南京航空航天大學 航天學院,江蘇 南京 211106; 2.空間核技術應用與輻射防護工業和信息化部重點實驗室,江蘇 南京 211106)

中微子是一種具有極小質量且相互作用極其微弱的粒子,對其探測具有極大的挑戰性。然而,中微子在物理學、天體物理學和核能等領域中都扮演著重要角色。反應堆是地球上最大的中微子源之一,可以用來深入研究中微子的特性和行為。反貝塔衰變(inverse Beta decay, IBD)反應[1],通常被作為探測反應堆中微子的首選反應。盡管根據該反應產生的信號具有時間關聯的特性,可以排除大量本底,但仍存在快中子本底和宇宙線誘導的9Li/8He本底以及偶然符合本底等進入判選時間窗[2]。因此中微子探測技術的關鍵在于對環境本底的抑制。

在以往基于液體閃爍體的實驗中,探測器通常位于地下深處,利用天然巖石的屏蔽以抑制宇宙線本底[3]。與液體閃爍體相比,塑料閃爍體具有便于運輸、結構穩定性和無毒性等優點,并且可以制作成陣列結構。探測器的陣列結構能使信號和本底具有不同的拓撲結構(粒子在陣列探測器中所產生信號在各個單元上的分布),文獻[4-6]采用陣列結構的塑料閃爍體探測器進行中微子探測。Panda[5]和Ismran[6]通過判斷信號在每個單元中的能量沉積的范圍來區分本底,但中微子事件的重建效率相對較低。Sertac[7]提出了基于機器學習的多變量分析方法用于中微子事例重建,提高了變量分析的效率。隨后,Ismran[8]實驗也計劃使用基于多變量分析的本底抑制方案。

雖然多變量分析方法的概念已被提出,但相關的實驗研究仍然沒有進展。本文設計了一個由3×3的塑料閃爍體陣列組成的小型探測器并耦合硅光電倍增管(silicon photomultiplier,SiPM)進行信號讀出。設計SiPM的讀出電路,對探測器的能量進行了標定,并測試了探測器的脈沖形狀甄別性能。在此基礎上,考慮中子的能量校正以及探測器的能量分辨率對探測單元能量測量的影響,利用252Cf源的中子和伽馬射線仿真和實驗驗證了陣列探測器中不同粒子存在不同的拓撲結構。根據該能量校正方法,計算了中微子信號與本底在探測器中產生的不同拓撲結構,并應用多變量分析方法對其進行分析。

1 陣列探測器系統設計及搭建

為研究多變量分析方法在塑料閃爍體陣列探測器中的應用,本文設計并搭建了一個3×3的小型陣列探測器。針對探測器可讀出信號的特性,為探測器設計了讀出電路,期間對該讀出電路的參數進行優化選取,測試了探測器的能量分辨能力與脈沖形狀甄別性能。

1.1 探測器的結構與選型

探測器實物如圖1所示。陣列探測器由9條尺寸同為14.4 mm×14.4 mm×100 mm的EJ-276塑料閃爍體以3×3陣列的形式組成。每個EJ-276閃爍體被涂有4.8 mg/cm2Gd2O3的鋁箔包圍,以提高中子俘獲效率并增加閃爍光的收集效率。外殼是2 mm厚的鋁,以防止可見光的進入,同時對閃爍體陣列進行固定。每條塑料閃爍體的兩端各耦合一個型號為ArrayJ-60035-4P的SiPM陣列,利用LEMO接口進行信號輸出。探測器的上下兩端有2條RG-176同軸電纜作為兩端SiPM陣列的供電輸入。利用2塊CAEN公司生產的基于VME機框式的V1730B型數字化采集儀進行數據采集,其采樣頻率高達500 MHz,具有時間戳等功能,可以記錄每個事件的時間信息,最大通道數為16道。

圖1 探測器實物Fig.1 The physical object of the detector

1.2 放大及成型電路設計

反應堆附近地面的輻射水平更高。探測器的陣列結構使得通道數量大幅增加,而大面積的SiPM陣列也會產生較寬的脈沖信號,進一步增加了數據采集、處理和存儲的難度。因此在電流靈敏前置放大器后耦合了基于極零相消原理的整形電路。其電路模型如圖2所示,其中反饋電阻Rf決定了放大的倍數,反饋電容Cf用于提高反饋的穩定性。當輸入信號的衰減時間常數τi等于R1×C時,整形后輸出脈沖波形的衰減時間常數從τi縮短為τo:

(1)

圖2 耦合極零相消電路的跨阻抗放大器模型Fig.2 Model of transimpedance amplifiers coupled with pole-zero cancellation circuit

式中:R1為極零相消電路中的串聯電阻阻值;C為極零相消電路中的串聯電容容值;R2為極零相消電路中的接地電阻阻值。

實驗中需要對脈沖信號進行積分以得到對應的粒子能量,過窄的脈沖波形將失去粒子的能量信息。電路參數在保持C為10 nF和R1為82 Ω不變,通過改變R2的阻值在不同程度上縮短脈沖寬度。其中將不加入脈沖整形模塊的測試結果作為對照組。圖3顯示了不同R2阻值時的脈沖形狀,與沒有成形電路的結果相比,具有更窄的脈沖。

圖3 不同R2阻值下的波形形狀對比Fig.3 Comparison of the pulse shape with different resistance values of R2

2 陣列探測器的優化與性能測試

探測器的脈沖形狀甄別性能是后續分析中子與伽馬信號的基礎。由于脈沖的整形會改變脈沖的形狀,可能會對探測器的脈沖形狀甄別性能產生影響。因此對以上5組電路參數下的脈沖形狀甄別性能進行了測試。

2.1 探測器的能量標定

探測器的脈沖形狀甄別性能與探測到的粒子能量相關,在進行脈沖形狀甄別能力測試前需進行探測器的能量標定。塑料閃爍體的密度和原子序數較低,光電峰很難出現在能譜中,無法通過擬合光電峰來獲得能量與道數的關系。由于康普頓邊緣并不等同于光電峰,通過高斯擬合康普頓邊緣難以獲得準確的能量標定結果,另外利用該方法得到的能量分辨率與實際值也有較大偏差。本文采用了高斯能量拓寬(Gaussian energy broadening, GEB)的方法逼近實驗值[9],可更準確地得到探測器的能量分辨率與康普頓邊緣能量對應的道數。

實驗中利用3種137Cs、60Co和22Na伽馬射線源進行探測器的能量標定。首先在GEANT4[10]中模擬了3種放射源在EJ-276閃爍體產生的理論能譜。利用GEB方法對理論能譜按能量分辨率以及康普頓邊緣的位置進行展寬,然后將其與實驗結果進行比較。利用最小二乘法計算擬合的相似程度:

(2)

式中:Aexp(Ch)和Asim(Ch)是實驗和模擬譜中通道Ch中的事件數。以R2為200 Ω時的能量標定結果為例。表1為當阻值為200 Ω時,探測器單個單元的137Cs能譜的標定結果。χ2的計算范圍為65 000～76 000個ADC道數,數據的選擇間隔為200個ADC道數。當能量分辨率(energy resolution, ER)EER為16%,康普頓邊緣(Compton edge,CE)ECE為66 000 ADC道數時,χ2達到最小值。圖4比較了最小χ2擬合時對137Cs源的實驗和模擬譜的比較。根據137Cs、60Co和22Na源的標定結果獲得能量與ADC道數的線性曲線,如圖5所示。由于閃爍體在低能量下的光輸出是非線性的,因此線性擬合結果存在一個較小的截距,其值為25 keV。

表1 不同康普頓邊緣能量與通道值和不同能量分辨率的χ2結果

圖4 最小χ2擬合時的模擬譜與實驗譜比較Fig.4 Comparison between simulated spectrum and experimental spectrum during minimum χ2 fitting

圖5 137Cs、 60Co和22Na 3種放射源利用展寬擬合法的標定結果Fig.5 The calibration results using the broadening fitting method for 137Cs,60Co and 22Na radioactive sources

2.2 探測器的脈沖形狀甄別

進行能量標定后,利用252Cf中子源分別測試了每組電路參數下探測器的脈沖形狀甄別性能,測試時長統一為30 min。使用電荷比較法(charge comparison method, CCM),利用入射粒子脈沖波形在不同時間區間積分電荷比值的差異進行甄別。對2個具有不同形狀的波形取相同的時間起點t0作積分,至某時刻ts短門積分為Qshort,至時刻tL長門積分為Qlong。2個不同的時間寬度對應長門和短門。在脈沖下降的起始部分以及脈沖的后尾沿部分往往區分程度較小。需通過長門與短門選取脈沖差異最大的區域進行中子與伽馬的區分。由此可計算變量RPSD為:

(3)

本文以R2為200 Ω時的測試結果為例。R2=200 Ω時探測器的RPSD值分布如圖6所示。粒子沉積的能量越大,中子與伽馬的區分越明顯。由于閃爍體和反射層對閃爍光子的吸收,使得低能事件被收集到的光子數少,導致其脈沖形狀的失真。另外低能的脈沖波形由于脈沖幅度低,其受后端電子學噪聲的影響較大。基于以上因素,能量較大的脈沖信號往往能保持較好的脈沖形狀甄別能力。塑料閃爍體對于不同種類的帶電粒子的光產額不同,而實驗中在無法首先判斷粒子種類的情況下,在測得能量后,一般將等效能量統計為等效電子能量,單位為MeVee。統計能量大于1.5 MeVee事件的RPSD值分布,得到RPSD值的分布如圖7所示。伽馬峰與中子峰,利用高斯分布擬合2個峰位分別得到對應的均值μn、μγ和標準差σn、σγ。通過品質因子(figure of merit,FOM)RFOM進行脈沖形狀甄別效果評估:

(4)

圖6 R2=200 Ω時探測器的RPSD值與能量的關系Fig.6 The relationship between RPSD value and energy of the detector at R2=200 Ω

圖7 中子和伽馬信號的RPSD值分布Fig.7 The distribution of RPSD corresponds to the neutron and Gamma signal

式中:Гn、Гγ分別為是中子峰與伽馬峰的半高寬,半高寬Г可表示為2.355倍的標準差。RFOM值越大,表明中子和伽馬越容易區分。表2顯示了不同整形電路下的RFOM值。結果顯示當R2為200 Ω時,探測器保持了最佳的脈沖形狀甄別能力,其RFOM值達到1.28。后續實驗均在R2為200 Ω時進行。

表2 252Cf中子源測試下脈沖形狀甄別的最佳RFOM值及對應高斯擬合結果

3 多變量分析方法

3.1 探測器中不同粒子的拓撲結構

3.1.1 GEANT4仿真

探測器的分段結構能夠根據事件的拓撲結構來分離信號和本底。粒子在探測器中的拓撲結構可以體現為粒子在陣列探測器各個單元產生能量的分布,即最大能量沉積單元Efirst、總能量沉積Etotal以及總沉積能量在0.2 MeV以上的響應單元數Nunit。在本節中,測量了252Cf 源產生的中子和伽馬射線在探測器中的事件拓撲結構和并驗證了能量校正方法。

在仿真中考慮到質子與電子在塑料閃爍體中的光產額不同,需要將對中子對探測器的能量響應進行校正,與實驗相匹配。對質子的能量校正函數為:

L(Ep)=A×Ep-B×(1-e(-C×Ep))

(5)

式中:L(Ep)為質子的等效電子能量;Ep為質子能量;A、B、C為實驗常數。EJ-276塑料閃爍體對應的A=0.75,B=3.2,C=0.22[11]。此外受探測器具有一定能量分辨的影響,實際探測到的能譜是理論能譜按能量分辨率進行展寬后的結果。探測器的能量分辨率根據2.1節的標定結果得到。在仿真中采用能量校正前后中子在Efirst、Etotal、Nunit3種變量的分布情況如圖8所示。

圖8 校正前后的中子在各變量上的分布對比Fig.8 Comparison of neutron distribution in various variables before and after correction

從中子各個變量在校正前后的對比可以看出,校正前后各變量有較大變化:Efirst的最高相對計數由1.35降低到0.7,沉積總能量分布范圍由0～5 MeVee變為0～4 MeVee(圖8(a))。Etotal也出現最高相對計數降低和沉積能量分布范圍變窄的情況(圖8(b)),變化幅度與Efirst基本相同。另外響應單元數為1的相對計數也從7下降到6.6(圖8(c))。經過校正后中子的能量相較于校正前有所減小,使得Efirst與Etotal2個與能量相關的變量的縱坐標壓縮。由于中子在探測器中產生的能量主要集中在低能部分,且在低能段校正函數的影響較低,因此響應單元數在2和3時占比的變化程度不明顯。

采用能量校正方法后中子與伽馬在變量上的分布結果如圖9所示。

圖9 仿真得到的中子伽馬事件在各變量上的分布Fig.9 The distribution of simulated neutron Gamma events on various variables

由于252Cf產生的中子和伽馬主要通過與探測器介質撞擊產生的帶電粒子在探測器中產生能量,其能量很難達到中微子反貝塔衰變反應在探測器中產生的能量范圍。從結果中可以看出252Cf產生的中子與伽馬的總能量范圍主要集中在0～3 MeVee。Efirst能量分布與Etotal類似,粒子大部分能量沉積在一個單元中。在1.5～3 MeVee能量區間內,中子在Efirst與Etotal上的相對計數略高于伽馬,如中子在Efirst為1.6 MeVee時的相對計數為0.25,伽馬則為0.15。另外中子在響應單元數為2的相對計數為1.8高于伽馬的1.3。這是由于中子在探測器內的徑跡較長并可能被Gd、H等元素俘獲,會在多個單元內產生能量;伽馬基本上只發生一次康普頓散射,所產生反沖電子的徑跡決定了響應單元的數量與能量沉積分布。

3.1.2 實驗測試

由于脈沖形狀甄別方法能夠進行中子與伽馬的區分,可以從實驗上獲得相應的中子和伽馬信號。利用探測器對252Cf放射源產生的中子與伽馬進行20 h探測,得到了符合能量范圍的70萬個伽馬事件與6萬個中子事件。圖10所示為中子與伽馬在各個變量上的分布,變量與模擬中的變量相吻合。在1.5～3 MeVee能量區間內,中子在Efirst與Etotal上的相對計數略高于伽馬,如中子在Efirst為1.6 MeVee時的相對計數為0.5,伽馬則為0.3。中子的響應單元數量總體上大于伽馬,如中子在響應單元數為2的相對計數為2高于伽馬的1.2。

圖10 實驗得到的中子伽馬事件在各變量上的分布Fig.10 The distribution of neutron Gamma events obtained from experiments on various variables

利用該修正方法得到的模擬結果與實驗結果基本一致。然而仿真與實驗結果仍然存在些許差異,這是因為利用單一的能量分辨率來模擬具有一定的誤差,另外利用其他尺寸探測器得到的中子校正結果,并不能完全精確地描述中子能量。為達到精確的能量信息,后續也需要對探測器進行中子的能量標定。

3.2 中微子信號與本底的多變量分析

本文考慮了難以屏蔽完全的快中子與伽馬2種本底,以3.1.1節采用的能量校正方法為基礎,模擬中微子信號與本底在探測器內產生的拓撲結構,并分析多變量分析方法對本底的抑制效果。將1～50 MeV均勻分布的伽馬作為本底,中子本底能譜則采用文獻[12]中地面實測的1～103MeV中子數據。中微子反貝塔衰變反應會產生一個次級正電子(快信號)與一個次級中子(慢信號),可以根據2種信號具有時間關聯的特性進行挑選[1]。因此在分析中微子信號的拓撲結構時,2種信號是相互獨立的。中微子反貝塔衰變反應產生的次級中子能量范圍在1～100 keV,次級正電子能量范圍在0～8 MeV[13]。

由于快中子會在探測器內部通過產生的反沖質子形成假快信號,中子被逐漸慢化后被Gd、H等元素捕獲將形成假慢信號,具有時間關聯特性。需要將快中子事件拆分為假快信號與假慢信號,分別與IBD快信號與慢信號進行分析。在仿真中統計反沖質子信號即可得到快中子產生的假快信號,統計中子俘獲反應發生后產生的電子信號即可得到快中子產生的假慢信號。而伽馬本底沒有時間關聯,其隨機出現在IBD判選的時間窗內,因此將伽馬本底與IBD快信號與慢信號分別進行分析。

圖11顯示了各變量下探測器的信號和本底分布。信號和本底在Efirst、Etotal和Nunit的分布中有所不同:IBD次級正電子、IBD次級中子、伽馬本底、中子假快信號和中子假慢信號在Efirst中最大相對計數對應的能量分別為1.75、0.95、0.85、0.55和0.65 MeV(圖11(a)),Etotal中最大相對計數對應的能量分別為2.05、1.83、0.85、0.66和1.95 MeV(圖11(b)),響應單元數為2的相對計數分別為0.82、0.57、0.53、0.58和0.60(圖11(c))。伽馬本底在探測器內的總能量沉積基本均勻分布在能量閾值內,其響應單元數為1的概率相比IBD信號與快中子本底更高,與3.1節的規律一致。IBD次級中子與快中子本底在Efirst與Etotal上的分布類似,在Nunit上的分布有所差別。這是由于IBD次級中子能量較低,信號主要來自于中子俘獲產生的次級伽馬,其過程與快中子引起的慢信號類似。兩者的主要區別在于中子發射的位置不同:IBD次級中子在探測器有效介質內均勻產生,快中子則從外部入射到探測器。IBD正電子由于在探測器內部均勻產生,主要參與電離相互作用,其能量基本沉積在探測器內,正電子在探測內的徑跡以及正電子湮滅產生的2條伽馬射線決定了響應單元數的分布,相較于IBD次級中子以及2種本底,IBD次級正電子更容易在多個單元內沉積能量。

圖11 IBD信號與本底在3種變量中的分布Fig.11 The distribution of IBD signal and background in three variables

決策樹(boost decision tree,BDT)是實驗物理分析中都有著廣泛應用的多變量分析算法,BDT將數據分割成節點,其中最終節點將事件分類為信號或本底。由于中微子信號為2種次級粒子,與本文考慮的2種本底,對應4種決策樹。在多變量分析中采用BDT算法綜合考慮各個因素,將中微子信號與本底進行分析,得到了信號與本底的BDT響應分布。其中次級正電子信號與2種本底的BDT響應如圖12 (a)、(b)所示。次級中子與2種本底的BDT響應如圖12 (c)、(d)所示。BDT算法通過整合計算,將信號的BDT響應值集中在正數部分,本底集中在負數部分:對于次級正電子信號與伽馬本底,信號的BDT響應值小于0的概率為35.72%,本底則為83.67%(圖12(a));對于次級正電子信號與快中子本底,信號的BDT響應值小于0的概率為17.35%,本底則為64.10%(圖12(b));對于次級中子信號與伽馬底,信號的BDT響應值小于0的概率為13.54%,本底則為61.57%(圖12(c));對于次級中子信號與快中子本底,信號的BDT響應值小于0的概率為60.09%,本底則為73.62%(圖12(d))。對比圖11信號與本底分布差異性更明顯。多變量分析方法將變量結合起來,相較于單獨判斷單個變量可以更好地區分信號與本底。但由于IBD次級中子與中子本底的物理過程類似,多變量分析方法對其區分效率較低,IBD次級中子信號與中子本底的區分程度并不顯著。

圖12 各組信號與本底的BDT響應分布Fig.12 BDT response distribution of each group of signals and background

圖13 不同信號與本底的信號效率和本底抑制關系Fig.13 Signal efficiency and background suppression relationships for different signals and backgrounds

根據 BDT 的響應分布結果,可以通過設置閾值,以最大化信號的識別效率并最小化本底的誤判率。探測器對中微子事例重建效率和本底拒絕效率的關系可以通過接收機工作特性( receiver operating characteristic,ROC) 曲線表征。ROC曲線通常用于展示在不同的判別閾值下,系統的真陽性率(true positive rate, TPR)和真陰性率(ture negative rate, TNR),即信號效率與本底拒絕效率之間的關系。在物理實驗中,理想的情況是保持高信號效率的同時,盡可能多地排除本底。因此,ROC曲線越接近右上角,表示信號與本底區的區分效果越好。4 種信號與本底組合的 ROC 曲線如圖 13 所示。

在應用4種決策樹時統一保持信號效率為90%,分析次級正電子信號與本底時,可以排除48%的伽馬本底和54%的中子本底;分析次級中子信號與本底時,可以排除57%的伽馬本底和16%的中子本底。即可以保留65%的中微子信號,排除77%的伽馬本底,與84%的中子本底。

4 結論

1)本文優化了探測器的信號讀出電子學,提高了其脈沖形狀的甄別性能,縮短了脈沖寬度。

2)實驗測試了252Cf源產生的中子與伽馬在探測器單元內產生的能量分布情況并與蒙特卡羅計算對比分析。考慮使用能量校正方法后,仿真結果與實驗結果達到較好的一致性,另外也證實了不同種類粒子在探測器內具有不同的拓撲結構。

3)本文利用多變量分析方法對信號與本底進行了區分。多變量分析方法可以保留65%的中微子信號,排除66%的伽馬本底與84%的快中子本底。

4)本文通過使用小型陣列探測器進行反應堆中微子探測,雖然需要長時間的測量,但證實了多變量分析方法可以有效抑制中微子本底。未來的研究可以考慮采用更大結構、更多單元的探測器進行進一步研究,以期望獲得更多實驗結果。