针对中微子探测的多变量分析本底抑制方法

郭一, 张余炼,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)本文通过使用小型阵列探测器进行反应堆中微子探测,虽然需要长时间的测量,但证实了多变量分析方法可以有效抑制中微子本底。未来的研究可以考虑采用更大结构、更多单元的探测器进行进一步研究,以期望获得更多实验结果。