杨振伟 清华大学 第十一讲：大作业介绍、计算机拟合方法 粒子物理与核物理实验中的数据分析 杨振伟 清华大学 第十一讲：大作业介绍、计算机拟合方法

Project 大作业 分为文字报告与口头报告两部分。 本学期有三个 projects 供大家选择： Project 1：光源位置重建定位 主要是利用最大似然法拟合发光点位置 Project 2：应用神经网络法甄别中子信号 主要侧重于如何处理较大批量数据，如何在 Root 平台应用神经网络法对信号与本底进 行甄别。 Project 3：应用解谱法还原真实电子能谱 主要侧重于如何应用 Geant4 进行简单探 测器模拟以及在 Root 平台进行解谱分析。

1：发光点位置重建 目的 要求 学会对数据作各种直方图与二维散点图 学会应用各种分布去研究各种图表和发现问题 学会如何应用MINUIT进行参数拟合研究 学会误差分析方法 学会如何做数据分析的最终研究报告 要求 报告需要有实验原理描述 各观察量分布与统计分析 对不同发光点数据进行位置拟合重建 系统误差分析 结论

超级神冈切仑科夫水探测器 在一只有效体积为高3620cm, 半径1690cm，重量约5万吨的纯净圆柱水桶中，一共安装了11146根光电倍增管，其中桶壁安装了安装了7650根，上下桶盖各安装了1748根。由于中微子与水中的核子发生相互作用可以产生带电 -子，而-子在水中飞行的速度大于光在水中的传播速度时，会发出切仑科夫光。根据光电倍增管所探测到的时间，我们可以确定中微子与水中的核子发生相互作用位置，也就可以确定中微子是否发生在给定的靶体积内。 42 m 40 m

激光装置 激光球发出各向同性的单色光，将激光球放置在不同位置并进行数据采集。

探测器坐标定义 中心 顶部 桶部 Z （cm） 底部

实验数据的基本信息 光在水中的传播速度（群速）为：21.6834 cm/ns 数据文件位于/data/laser/idcal_lsr，ROOT格式。 数据内容位于skdata这个tree中，分为四个branches： skhead: 事例基本信息 skq: 与光电倍增管对应的电荷测量信息 (电荷单位 p.e.） skt: 与光电倍增管对应的时间测量信息 （时间单位 ns） skchnl: 击中信息 详细的每个leaf的说明参见/home/chensm/laser/ 另外，已知 光电倍增管位置：/home/chensm/geom/geom_sk3.dat 光电倍增管分辨率：/home/chensm/resolution/pmt_sig.dat 光在水中的传播速度（群速）为：21.6834 cm/ns 超级神冈IV共11146个电子学道，其中17道没有对应的光电倍增管，约50道被标识为坏的光电倍增管

各批次（run）数据的基本信息 F值越小，代表激光球强度越高，标有loose_c的情况下，激光球强度同比更弱

基本要求 根据所给的数据，作出 了解数据文件中每个leaf的意义 分析不同光强下的击中数，总电荷量，平均电荷量，给出相应的直方图 从桶部和上下桶盖分别选出一只光电倍增管，针对不同光强，作出对应的时间测量直方图，电荷测量直方图以及时间电荷二维散点图 比较不同数据的分布，对图形分布进行解释与说明，并尽可能作定量说明 利用光在水中传播的速度 21.6834 cm/ns，并根据设定的光源点位置与每个光电倍增管位置，计算光从设定位置达到每只光电倍增的飞行时间，作出时间随光电倍增管号码的变化图。 针对不同光强下的几个事例，对有击中的光电倍增管，做出时间测量值减去飞行时间随光电倍增管号码的变化图，随飞行距离的变化图

研究光源在中心的数据 利用扩散球在中心的数据，对不同光强，作出 利用激光球在中心, 发出的光强基本各向同性的特点, 估计桶部其中一圈光电倍增管各自对光的探测效率, 并作出探测效率与光电倍增管号码的关系分布图(提示, 每一个事例均是有激光从激光球中发出)。 利用各光电倍增管测量的时间, 减去光从激光球到各相应光电倍增管所用的时间, 可得到时间测量的残差t。作出桶部其中一圈光电倍增管t的分布图，并说明它们应该服从什么分布？各对应的平均值与方差开方值是多少？并作出平均值，方差开方值与光电倍增管号码的关系分布图。 尽量减少电子学噪音等因素的干扰，以及反射光和散射光的影响，重做2中要求的分布图。 由于光电倍增管时间测量精度与所接收的光电子数目是相关的。根据数据任选一只光电倍增管，分别作出光电子数大于5与小于5的时间测量的残差t，验证上述相关关系。

影响时间分布的因素研究 在研究光电倍增管对时间测量的分辨率时，我们需要比较实际观测的时间分布与理论假设是否吻合，从而找出是否存在影响时间测量精度的其它因素。 请任意挑选一只光电倍增管, 利用扩散球在中心的数据，针对不同电荷量，作出时间残差分布图，找出平均值。然后， 随机产生以该平均值为中心值的高斯分布，每个事例的时间分辨率大小估计可以利用 /home/chensm/resolution/ pmt_sig.dat文件中的电荷-时间分辨关系曲线来描述。作出蒙特卡罗模拟给出的时间残差分布。将实验数据的直方图与蒙特卡罗模拟的直方图画在一起，比较两者之间在分布宽度上有什么区别？能否说明时间测量的误差纯粹是由光电倍增管自身的原因造成的？

似然估计量的构造与参数拟合 利用扩散球在中心，中等光强(F ~ 1000)的数据。根据理想状况下，任何一个光电倍增管对时间的测量，在扣除飞行时间以及全局时间偏移以后，时间残差近似服从均值零的，标准误差为t(Q)假设，构造似然函数。以光源的位置(x0,y0,z0)为拟合参数，写出拟合程序，给出光源位置的估计值。 利用该拟合程序处理低光强和高光强数据，看拟合顶点的展宽和均值有何变化 考虑反射光，散射光以及电子学噪音等影响因素，适当加时间窗选择击中，再看顶点拟合结果有何变化。 考察每次顶点拟合后得到的 2，若不接近于1，说明什么？分别对大于1和小于1的情况做出解释。调整时间分辨率为原值的10倍，查看拟合结果的变化。 若拟合顶点值偏离真值，请分析其原因并给出相应解释。

置信区间与结果报告 根据 MINUIT 的说明给出不同方法误差的估计； 画出每两个参数的 68% 等高线； 分析实验中可能的系统误差来源； 书面报告

2：神经网络方法寻找中子信号 目的 要求 学会处理较大批量数据的方法 学会在 root 平台应用神经网络方法进行信号识别 学会作直方图、直方图归一化与它们之间的运算 掌握直方图运算中的误差传递原理 学会如何在 root 平台进行直方图拟合与结果评价 学会如何做数据分析的最终研究报告 要求 报告需要有实验原理描述 各观察量在信号与本底方面的分布与比较 找出凸显本底与信号分布的差异的神经元 利用模拟信号样本与本底数据样本训练神经网络 分析中子数据样本，直方图拟合给出中子在水中的寿命 误差分析 结论

关联中子在纯水中的探测原理 5 cm Am/Be 计时时间起点 计时时间终点 中微子反应截面太小，不能用于探测原理的检验，实验上采用放射源来模拟类似的反应，检验关联中子的探测方案。 Gamma 在本实验中只能通过康普顿散射打出电子，再由电子发出切伦科夫光后才能被光电倍增管观测到。 分析关键点：如何从数据分析中证明探测到了中子。

超级神冈实验探测器 41.4 m  39.3 m A 50k tons water Č detector located at 1k m underground 探测器分为内外层： 外层是直径为 2 米 的水，用来屏蔽环境 产生的中子，内层为 中微子探测器。 中子  Am/Be 1 km

实验数据 信号：通过蒙特卡罗模拟得到 本底：无放射源状态下，实验获取数据 含中子数据：探测器置入放射源，实验获取数据 运行号 放射源 063923 063926 放射源 有 无 装置的位置 (x,y,z) cm 35.3 -70.7 0 - - - 数据存放处：/data/AmBe/AmBe/06392*/ 不考虑噪音情况下的纯信号模拟数据： /data/AmBe/AmBe/mc/gamma 加入5 kHz噪音的信号模拟数据： /data/AmBe/AmBe/mc/gammanoise

数据分析基本原理（1） Because of time-of-flight difference to individual PMT, the PMT timings of 2.2MeV  can not form a peak against BG. 2.2MeV  # of hits Averaged BG PMT time # of hits Averaged BG PMT time 因此数据分析第一步是要做飞行时间修正以恢复信号峰。但是做飞行时间修正需要 顶点信息，而2.2 MeV gamma信号太弱而不能进行有效的顶点重建。解决方法是利 用已知的源的顶点代替，因为我们知道中子跑的离源不会太远（见问题1）。

数据分析基本原理（2） 分析步骤： 1. 利用源的顶点做飞行时间修正 2. 利用10 ns滑动的时间窗选出信号峰 10ns内击中数(N10) 大于某个值时即认为是信号。为什么用10ns窗呢？N10的阈值又怎么确定？见问题2和3。 3. 由于PMT噪音以及其它低能本底的偶然符合，利用步骤2中的方法得到“信号”中依然含有大量本底。此时可以利用更多信息以对信号和本底做进一步的分离。信号的特征是什么？本底的特征又是什么？找出尽可能多的刻画这些特征的变量。 做变量的分布图并研究变量之间的相关性，找出甄别能力大的几个变量来。 然后利用这些变量来构建神经网络。

基本问题 右表是中子在不同能量下弹性散射的截面，又知中子被水中 质子俘获的截面是0.332 barn。试根据这些数据估算一个3 MeV 中子在俘获前离源的平均距离，以证实利用源的顶点来做飞行 时间修正的合理性。 提示：中子首先要慢化变成热中子，热中子在被俘获前还要发生 多次弹性碰撞（这一点从俘获截面和弹性散射截面数值上可以看 出来）。因此中子离源的直线距离并不简单地等于平均自由程。 此时的热中子是做随机游动（random walk)，唯一可求的是方均 根距离，所得结果比平均自由程小得多，可参见费曼物理学讲义 第一卷41章。类似的推理可以估算慢化阶段中子跑的距离。 2. 利用纯信号数据研究2.2 MeV gamma信号的基本特征： 1）做击中数的分布 2）电荷值分布 3）T – TOF的分布， 注意其宽度，以验证10 ns窗的合理性， 并解释分布的左右不对称。 3. 研究N10的阈值。阈值越高，本底越少，但信号丢失的越多。 尝试4，5，6等不同的值，研究信噪比。假如仅根据N10，信噪比 能达到什么水平？试做出信号效率 VS本底概率的曲线。 进一步研究信号和本底的特征，尝试找出除N10之外的其它甄 别量并利用神经网络做多元分析。再做信号效率 VS本底概率的 曲线并同3中结果比较。

数据中的基本信息 与本 project 有关的探测器主要信息为： 已知：光在水中的传播速度（群速）为：21.66 cm/ns class Header { int nrunsk; // run # int nevsk; // event # int idtgsk；// trigger id … }； class TQReal { Int_t nhits; // # of hit PMT vector<int> cables; // cable # vector<float> T; // timing vector<float> Q; // charge … }； 数据详细定义参见：/data/AmBe/AmBe/definition/* 已知：光在水中的传播速度（群速）为：21.66 cm/ns 一共有 11,129 根光电倍增管 光电倍增管位置：/home/chensm/geom/geom_sk3.dat 光电倍增管分辨率：/home/chensm/resolution/pmt_sig.dat

数据预处理 导出 ROOT 数据分析程序模板 .C 和.h文件,可以参考 /home/chensm/neutron/makeSelector.C 利用 >root –l makeSelector.C 参看生成的文件里的信息。 比较 /home/chensm/neutron/ref 中对应的程序，理解如何在 ROOT 中分析一个事例的流程，以及如何添加和删减分支 branch，在正式分析数据前去掉无用信息，减少以后分析所占用的 CPU 时间与内存。 在参考该目录下阅读 runAnalysis.C 文件，理解如何处理多个 root 数据文件的方法。 根据要分析的 run 号对应得装置位置，修改 histAnalysis.C 中 VX，VY 与 VZ。 用 root –b runAnalysis.C 提交作业生成简化了的数据文件（在 /home/chensm/neutron/ 中查找）。 用同样程序处理模拟的信号样本，注意改变输入输出文件名。 比较本底与模拟的 nhits,ani,dks,tdiv,dmean,ddiff 分布，在每个量的分布图时，需要用不同颜色和虚实线在同一画板中同时画出信号与本底的分布。

训练神经网络 参考 root 说明书给出的例子，例如 /projects/chensm/chensm/nn/mlpHiggs.C /projects/chensm/chensm/nn/mlpHiggs.root 先运行 >root –l mlpHiggs.C 观察训练过程中 root 给出的各种分布，理解神经网络训练的过程。 对例子中的程序加以改造，替换神经元名字，利用 TChain *signal = new TChain(“t_nn_sam”); signal->Add(“./signal.root”); signal->SetBranchAddress(“nhits”,&nhits); … 改写程序并运行，参看输出的神经网络函数.cxx和.h文件，理解神经网络函数。

从信号样本中选择事例 把神经网络程序放置在预处理数据过程采用的程序 histAnalysis.C 中，在完成计算 nhits、 ani, dks, tdiv, dmean, ddiff 处理后，输入神经网络中进行甄别，将神经网络输出量以及事例与原初事例的时间差存入新建的分支 branch 中，例如 t_nn_val=new TTree(“t_nn_val”,”NN value output”); t_nn_val->Branch(“nn_val”,&g_nn_val,”g_nn_val/D”); t_nn_val->Branch(“Td_sci”,&g_Td_val,”g_Td_val/F”); 运行程序产生只有神经网络输出量和中子寿命两个变量的 root 数据文件。 作出不对神经网络输出量大小进行要求的时间差分布图； 分别要求神经网络输出量大于 0.95 与 0.99 是对应的时间差分布，比较分布形状的变化。

研究不同放射源位置数据与本底 对不同放射源位置的数据重复前面关于预处理数据、神经网络训练与信号选择，观察时间差分布的变化。 用同样的神经网络函数处理有同样装置但无放射源的本底数据，观察要求不同神经网络输出量大小时，时间差分布于信号数据样本的区别，对事例进行归一化后把信号与本底分布画在同一个画板上，并加以说明。然后用信号分布减去本底分布得到纯中子信号时间分布 采用分区最大似然法拟合减去本底的时间分布，给出不同放射源位置测量的中子平均寿命。 尽可能研究测量中可能包含的系统误差，给出结果。

3：解谱法还原真实电子能谱 目的 要求 学会利用Geant4进行简单的探测器模拟与数据读出 学会对数据作各种直方图与二维散点图分析 学会如何在 root 平台应用解谱法还原真实谱分布 学会如何比较直方图并定量给出差异的结果 学会如何做数据分析的最终研究报告 要求 报告需要有实验原理描述 各观察量分布与统计分析 做出事例显示图并检验模拟程序的正误 采用射程与能量直接转换方法作出电子能谱分布图 应用解谱法还原电子能谱并与真实能谱比较 定量给出采用和不采用解谱法时，能谱测量的精度 结论

实验背景 拟议实验可在散列中子源实验（美国橡树岭实验室，英国 卢瑟福实验室、中国广东东莞中国散裂中子源）进行。

物理背景 中微子在物质中可能与中子反应生成质子和电子(反β衰变): 中微子与物质相互作用的微分截面对超新星爆发等物理非常重要。为了测量反应的微分截面，我们需要知道参与反应的中微子的能谱，其前提是获得反应产物之一——电子的能谱。 从实验上，如果中微子与铅或者铁相互作用，我们无法直接测量电子能谱，相对来测量电子在靶物质中的射程是可行的。 于是问题变成：如果实验中测量到电子的射程，如何得到该电子的能量？ 假如我们可以很好地通过测量电子的射程得到电子的能量，那么我们就可以建造相应的中微子探测器测量微分散射截面。

关键问题和目的 由于电子在铁/铅等物质中的多重散射现象，给定能量的电子的射程具有很大的涨落，一般而言，其射程与入射能量没有很好的线性对应关系。也就是说，通过测量到电子的射程为 L 时，通常不能精确地得到电子的初始能量 Einit （即使从统计意义上）。 尝试采用解谱法解决这个问题。

问题的简化处理 该Project只考虑电子能谱还原部分，不考虑如何由电子能谱还原中微子的能谱。Geant4模拟的起点为已知能量分布的电子。 不考虑探测器几何的具体细节，简化探测器的几何设置，并认为探测器对粒子的射程是非常理想的。

主要步骤 定义理想探测器 用单能粒子进行模拟，并通过敏感探测器记录有用的信息，检验探测器几何、物质以及数据读出是否正确 用单能电子进行模拟，给出平均射程以及方差，从而验证无法由射程直接反推能量。 产生指定能谱的电子并进行模拟，尝试用解谱法还原电子能谱。

步骤 1 定义理想探测器 探测器可以直接定义为立方体，材料为铅/铁/铝。可以将整个探测器设为敏感探测器，从中记录需要的信息，主要包括： 射程 能量沉积 粒子的初始信息(能量、动量等)

步骤 2 验证探测器定义以及数据读出是否正确 1) 将探测器几何画出来 2) 将入射粒子设为 +介子，入射动量为205 MeV，将探测器材料设为塑料闪烁体。模拟大约 2000 事例，分别做出射程与能量沉积的直方图，用高斯拟合，看是否为 30.5 cm 和 108.4 MeV。 3）将入射粒子设为 +子，入射动量为 235 MeV,分析模拟分布是具有均值否为 54 cm 和 152 MeV。

步骤 3 用单能电子模拟得到平均射程和方差 将探测器材料分别设为铅/铁/铝。在1-50MeV范围内均匀分10-20个区间，用单能电子多次模拟，得到不同能量下的平均射程和方差，做出直方图。 对三种不同材质都进行模拟，对三种不同材质的结果进行比较。给出为什么无法直接由射程给出能量

步骤 4-1 产生指定能谱的电子进行模拟，解谱法还原电子能量 将探测器材料分别设为铅/铁。 产生1-50MeV范围内的均匀分布，进行模拟，从而得到响应矩阵R。画出电子能量与射程的分布以及它们之间的二维散点图。

步骤 4-2 产生指定能谱的电子进行模拟，解谱法还原电子能量 求反应矩阵 R 的逆矩阵 R-1。 重新产生一组电子能量为 1-50 MeV 范围内的均匀分布，模拟得到电子射程，利用前面得到 R-1 还原电子能谱，画出还原能谱的直方图以及原始的均匀分布，并进行比较，说明引入光滑函数进行正规化处理的必要性。

步骤 4-3 产生指定能谱的电子进行模拟，解谱法还原电子能量 利用 RooUnfold 软件包或者自己写 Unfold 程序，对前面均匀分布的模拟结果进行解谱，研究说明正规化参数的选择。做出还原后能谱与原始能谱的对比图，并将两个能谱相除，给出还原的精度。 RooUnfold是基于ROOT平台的解谱包，见下面网页： http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html

步骤 4-4 产生指定能谱的电子进行模拟，解谱法还原电子能量 产生 Michel能谱的电子(1-50MeV)，进行模拟。 利用 RooUnfold 软件包或者自己写 Unfold 程序，对前面 Michel电子能谱的模拟结果进行解谱，研究说明正规化参数的选择。并画出还原后能谱与原始能谱的对比，给出两个直方图的相除结果以及还原精度。 定量给出解谱还原电子能谱的精度，并给出结论。

一个实用的求函数最小值程序 核与粒子物理研究中，大家普遍使用的求函数最小值程序是 MINUIT 软件包 http://hep.tsinghua.edu.cn/~chensm/lectures/minuit.pdf 求 log L 最大值等效于求 –log L 的最小值。 它可以对目标函数为似然函数、最小二乘函数或用户自定义函数求极值，并提供几种最小化过程和相应的误差分析。最初为Fortran程序，后改为 C++，用于当时西欧核子中心(CERN)的实验数据分析，现也被应用于粒子物理以外的领域。主要有三种使用方式： 在 MINUIT 框架下单独使用(适于data-driven 模式)； 在物理分析工作站 (PAW) 环境下互动调用 MINUIT (基于Fortran)； 在 PAW 升级软件包 ROOT 环境下互动调用 MINUIT (基于C++)。 39

在 PAW 平台调用 MINUIT program MINUIT_FIT implicit NONE integer npar parameter (npar=2) character*10 chnam(npar) integer ierr, ird, isav, istat, ivarbl, iwr integer npari, nparx double precision arglis(10), bnd1, bnd2, deriv(npar), dpar(npar) double precision fmin, fedm, errdef, covmat(npar,npar), log_l external FCN double precision par(npar) c Initialize MINUIT, set print level to -1 ird = 5 ! unit number for input to Minuit (keyboard) iwr = 6 ! unit number for output from Minuit (screen) isav= 7 ! unit number for use of the SAVE command call MNINIT(ird,iwr,isav) arglis(1)=-1.d0 call MNEXCM(FCN,'SET PRIN',arglis,1,ierr,0) c Define parameters alpha and beta, give initial values and step sizes call MNPARM(1,'alpha',0.5d0,0.1d0,0.d0,0.d0,ierr) call MNPARM(2,'beta ',0.5d0,0.1d0,0.d0,0.d0,ierr) c Get input data by calling FCN with iflag=1 arglis(1)=1.d0 call MNEXCM(FCN,'CALL',arglis,1,ierr,0) c Minimize using SIMPLEX and MIGRAD, get covariance c matrix with HESSE call MNEXCM(FCN,'MIGRAD',arglis,0,ierr,0) call MNEXCM(FCN,'HESSE',arglis,0,ierr,0) call MNEXCM(FCN,‘MINOS',arglis,0,ierr,0) c Get result of fit (for least squares, fmin is chi2) call MNSTAT(fmin,fedm,errdef,npari,nparx,istat) call MNPOUT(1,chnam(1),par(1),dpar(1),bnd1,bnd2,ivarbl) call MNPOUT(2,chnam(2),par(2),dpar(2),bnd1,bnd2,ivarbl) call MNEMAT(covmat,npar) log_l=-0.5*fmin write(6,*)'Fit results:' write(6,*) write(6,*)'alpha =',par(1),'+/-',dpar(1) write(6,*)'beta =',par(2),'+/-',dpar(2) write(6,*)'cov[alpha,beta]=',covmat(1,2) write(6,*)'rho[alpha,beta]=',covmat(1,2)/(dpar(1)*dpar(2)) write(6,*)'log_l =',log_l stop end 编译:f77 minuit_fit.f -o minuit_fit `cernlib` <return> 运行:minuit_fit <return> 40

在 PAW平台调用 MINUIT 用户应提供似然函数 FCN 注意：概率密度函数需要在有 效测量范围进行归一化处理。 如果计算有困难，可以将归一 化因子作为参数来拟合。这种 折衷方法虽简便，但会给其它 待定参数的确定带来误差。 subroutine FCN(npar,grad,chi2,par,iflag,futil) c Input: integer npar number of parameters to fit c double precision par(npar) parameter vector c integer iflag select what to do c double precision futil optional external function c Output: double precision grad(npar) gradient vector (not filled) c double precision chi2 function to be minimized implicit NONE integer npar,iflag double precision futil,chi2,par(*),grad(*) integer n_max parameter (n_max=10000) integer i,n double precision alpha,beta,f,log_l,x(n_max),angle_cut,f_nor c Begin n=2000 angle_cut=0.95 if(iflag.eq.1)then ! get n, array x call GET_INPUT_DATA(x,n,n_max,angle_cut) end if c Calculate log-likelihood alpha = par(1) beta = par(2) log_l = 0. c Normalization factor for [-angle_cut,angle_cut] f_nor = 2*angle_cut+2*beta/3.*angle_cut**3 do i=1,n f=(1.+alpha*x(i)+beta*x(i)**2)/f_nor log_l=log_l+DLOG(f) end do chi2=-2.*log_l ! 2 gets errors right return end 41

ROOT 直方图与 MINUIT 调用 root>.x minuit_fit.C 在 ROOT 环境下的运行程序 // Now ready for minimization step arglist[0] = 500; arglist[1] = 1.; gMinuit->mnexcm("MIGRAD", arglist ,0,ierflg); gMinuit->mnexcm("HESSE", arglist ,0,ierflg); gMinuit->mnexcm(“MINOS", arglist ,0,ierflg); // Print results Double_t fmin,fedm,errdef,covmat[npar][npar]; Double_t alpha,alpha_err,beta,beta_err; Int_t nvpar,nparx,icstat; gMinuit->mnstat(fmin,fedm,errdef,nvpar,nparx,icstat); gMinuit->GetParameter(0, alpha, alpha_err ); gMinuit->GetParameter(1, beta, beta_err ); gMinuit->mnemat(&covmat[0][0],npar); Double_t rho=covmat[0][1]/(alpha_err*beta_err); cout <<"alpha = "<<alpha<<" +/- "<<alpha_err<<endl; cout <<"beta = "<<beta <<" +/- "<<beta_err<<endl; cout <<"cov[alpha][beta]= "<<covmat[0][1]<<endl; cout <<"rho[alpha][beta]= "<<rho<<endl; cout <<"log_l = "<<-0.5*fmin<<endl; } const int npoints=2000; Double_t x[npoints]; Double_t angle_cut=0.95; void minuit_fit() { // Get data points get_input_data(); // Prepare for fit const int npar=2; TMinuit *gMinuit = new TMinuit(npar); gMinuit->SetFCN(fcn); Double_t arglist[10]; Int_t ierflg = 0; arglist[0] = 1; gMinuit->mnexcm("SET ERR", arglist ,1,ierflg); // Set starting values and step sizes for parameters Double_t vstart[npar] = {0.5, 0.5 }; Double_t step[npar] = {0.1 , 0.1 }; gMinuit->mnparm(0, "alpha", vstart[0], step[0], 0,0,ierflg); gMinuit->mnparm(1, "beta", vstart[1], step[1], 0,0,ierflg); root>.x minuit_fit.C 42

MINUIT中的似然函数 似然函数 void fcn(Int_t &npar, Double_t *gin, Double_t &chi2, Double_t *par, Int_t iflag) { // Calculate log-likelihood Double_t log_l = 0; Double_t alpha,beta,f_nor,f; alpha = par[0]; beta = par[1]; f_nor = 2*angle_cut+2*beta/3.*angle_cut**3; for (Int_t i=0;i<npoints; i++) { f=(1.+alpha*x[i]+beta*x[i]**2)/f_nor; log_l += TMath::Log(f); } // 2 gets errors right chi2=-2.*log_l; 43

如何确认拟合结果的合理性 在粒子物理与核物理实验的参数估计中，利用计算机进行数值计算求极值越来越普遍。 好处是可以较快地得到参数的估计值与误差。 但是，随着越来越多实验研究采用类似的参数拟合程序，对拟合结果的合理性要求越来越高，尤其是在相同数据条件下，因分析不同而造成结果不同的情况屡见不鲜，对确认拟合结果合理性的研究要求在不断地提高。 目前在如何正确使用MINUIT程序包时已经有一定的共识。

使用 MINUIT 应注意的问题 求极值问题时，不可能自动找到合理的初值。 原因：似然函数 存在多个极值。 -log L 区域最小值 参数值 区域最小值 真实最小值 原因：似然函数 存在多个极值。 提供较好的参数不确定范围可以有助于找到真值。 太大的不确定范围会使MINUIT碰巧在远离真值处找到 一个区域最小值。

MINUIT中的解决方案 利用MINUIT函数 MIGRAD 找 –log L 的最小值 利用MINUIT函数 HESSE 计算参数的误差 还可利用MINUIT函数 MINOS 做更好的误差估计

MINUIT函数 MIGRAD

函数 MIGRAD（续一）

函数 MIGRAD（续二）

MINUIT函数 HESSE

函数 HESSE （续一）

函数 HESSE （续二）

函数 HESSE （续三）

MINUIT函数 MINOS 用途：更好的误差计算方法 注意：有太多参数时可能会耗费非常多的CPU时间 调用：在ROOT 或 PAW 使用“E” 选项

经常会碰到的问题：不收敛 时常会出现参数拟合过程不收敛，例如 很多情况下是 HESSE 相关系数矩阵是很好的调查起点 MIGRAD 不能找到最小值 HESSE 得到负的二阶导数 很多情况下是 参数拟合函数有误（程序书写错误或理论模型不正确） 数值计算中出现精度问题 数值计算中出现稳定性问题 HESSE 相关系数矩阵是很好的调查起点 表明两个参数几乎100%关联，求最小值过程无法建立。

消除影响稳定性的因素：相关性 策略一：更恰当地选取拟合参数 例如，对于采用类似的双高斯宽度拟合 HESSE 相关系数矩阵 事例数 HESSE 相关系数矩阵 双高斯宽度的大小与它们 所占的份额有很强的关联

影响稳定性因素：相关性（续） 策略二：固定所有高度相关的参数而只保留一个 可以采用不同参数化过程来解决 第二个高斯分布宽度与份额比的相关度从0.92减小至0.68 完全是参数化过程自身的问题 策略二：固定所有高度相关的参数而只保留一个 如果参数高度相关表明它们中有些是多余的，不会对 拟合模型的自由度带来贡献。

消除影响稳定性的因素：多项式 注意：普通多项式 a0+a1x+a2x2+a3x3 的参数化过程 例如，导致拟合稳定性问题，通常不能准确地找到高阶 系数的解。 解决方案：采用数学上的正交多项式， 例如，勒让德多项式，第一类契贝谢夫多项式等等

参数拟合问题：参数的限界 有时候需要对拟合参数进行区间的限制 但是，可能会导致 MINUIT 不能正确估计误差。 例如，份额比参量只能在【0，1】区间，三角函数值 只能在【-1，+1】之间，等等。 但是，可能会导致 MINUIT 不能正确估计误差。 因此，必须小心使用参数的拟合区间限制。 参考的解决方案： 如果区间限制的引入可使拟合稳定，最好考虑寻找有否其它可替代的无区间限制的参数化方案。 如果区间限制的引入可以避免出现“非物理”但仍为合理统计解时，应考虑不要强加区间限制，并对结果采用统计意义上的解释。一般是同时给出物理解与“非物理”解。前者是使结果易于物理解释，后者是使结果可与其它实验结果进行统计意义上的正确合并。

举例：加速器中微子振荡实验 Phys.Rev.D74:072003,2006

如何确认拟合结果的有效性 所谓的结果有效性是指 如果对结果的有效性有疑虑，可以考虑 拟合结果无偏 拟合误差符合统计不确定性 进行P-值计算，估计结果是否为极端情况。 采用蒙特卡罗样本，重复同样统计量的参数拟合过程至少100次，确认拟合的参数值是否有偏向性。如果有偏向性，例如，发现有偏差量 a±b，应考虑对结果进行 a 因子修正，并把 ±b 当作系统误差。 检查蒙特卡罗研究中均值的妳散是否与拟合参数的误差符合，以确定误差的合理性。

举例：拟合信号数目 一维直方图信号的拟合 研究拟合参数 Nsig 信号服从高斯分布 本底分布已知 产生 Nsig 蒙特卡罗模拟研究1000次 产生的信号和本底数目分别为 1000次结果分布如右图 结果表明拟合没有偏向性 Nsig=100， Nbg=200 Nsig(fit)

举例：拟合结果的“Pull”分析 如何判断误差是否合理？ 解决方案是检查“Pull”分布 从蒙特卡罗研究中很难进行诠释 (Nsig) 从蒙特卡罗研究中很难进行诠释 没有与 Nsig 等效的分布研究误差 解决方案是检查“Pull”分布 定义 Pull 的特点 本例是一个在统计精度内无偏而且误差正确的拟合。 Pull(Nsig) 如果拟合结果无偏则均值为0； 如果拟合误差正确则宽度为1。

低统计下的拟合结果检查 对低统计样本或者大样本小信号情况需要特别注意 可能出现的问题 一般地，偏置的出现使得拟合误差正确性会有问题 似然评估量不再有效，导致从二阶导数中估计误差不再准确 最小二乘法中误差不再服从高斯分布而是泊松分布 偏置项在评估量中受统计量影响1/N部分不再可以忽略 一般地，偏置的出现使得拟合误差正确性会有问题 采用点估计而不是分区的最大似然法拟合 尽最大可能检查拟合的有效性

举例：小样本情况的偏向性 低统计量情况举例 蒙特卡罗模拟研究 必须对结果进行修正并引入恰当的误差。 Nsig=10020， Nbg=200 Pull 的均值已经不再是 0，有2.3 的偏离 分布在低统计下已 经不再是对称分布 Nsig(fit) (Nsig) Pull(Nsig) 必须对结果进行修正并引入恰当的误差。