粒子物理与核物理实验中的数据分析
合集下载
粒子物理与核物理实验中的数据分析剖析
任何估计量(不仅仅是最大似然法)的方差下界为
V
[ˆ]
1
b
2
E
2 log
2
L
( b为偏置)
也称为 Rao-Cramér-Frechet 不等式(信息不等式)。
如果等式满足，就可以说 ˆ 是有效的。
最大似然估计量对大的样本统计量 n 几乎总是有效的。 通常假设上述结论为真，利用RCF边界估计 V[ˆ]
6
估计量的方差：蒙特卡罗方法
通常情况下，ˆ 的具体形式 g(ˆ; , n) 并不知道。对于此类情况，
可采用蒙特卡罗方法得到 g(ˆ; , n)
例如，对指数的pdf，我们有 ˆ=1.062。在蒙特卡罗中，把其作 为 的真值。产生 n 50 的样本，并重复1000次实验。计算 每次实验的 ˆ,并填入直方图。
✓ 如果已经知道概率密度函数 pdf 的具体形式，但是函 数中包含有未知参数，如何从有限的样本中估计未知 参数的期待值、方差与相关系数…
✓ 基于假设为真时，会使观测结果的概率最大，构造似 然函数的方法
指数参数确定举例
10/16/2020
2
高斯概率密度函数中的参数
考虑一个样本服从高斯概率密度函数，其参数 ， 2 未知。
粒子物理与核物理实验中的 数据分析
杨振伟 清华大学
第八讲： 最大似然法（续）
上一章回顾
估计量，均值，方差与协方差的估计量
✓ 给出了在不知道概率密度函数 pdf 的情况下，如果没 有未知参数，如何从有限的数据样本中，估计出随机 变量的期待值、方差与相关系数。
✓ 讨论了偏置问题
似然函数，最大似然估计量
指数函数例子
也就是
log
L(ˆ
ˆˆ
粒子物理与核物理实验中的
= E[r2 ] − E2[r]
=λ
举例：光电倍增管暗电流影响
在有11146根PMT的探测 器中，已知每根PMT暗电 流产生的误击中为3.5kHz。 求探测器在任意总长度为 500μs时间段观察到每隔 10ns PMT误击中数目分 别为5和6的总次数
在10ns间隔观测到PMT误 击中的平均数目为
…
得到(n1,n2,…,nm)概率为
f
G (n)
=
N! n1!n2!...nm!
p n1 1
p n2 2
...
p n3 m
平均值 : E[ni ] = Npi
方差 :V[ni ] = Npi (1− pi ) 协方差 :Vij = −Npi p j (i ≠ j)
适适用用于于直直方方图图 频频数数误误差差估估计计。。
MRPC记录的击中数目N’
MRPC探测效率 测量值及其误差
从二项式到多项式分布
类似于二项式分布,但允许结果的可能性m大于两种,概率为
G p = ( p1, p2..., pm )
m
∑ pi = 1
尝试N次，结果为
可能性1：n1 可能性2：n2
G n = (n1, n2 ,..., nm )
i =1
举例：角分布中的前后不对称
e+
e+
θ e-
+ e− → J /ψ → e+ + e− -1 B：后向计数；F：前向计数；N=B+F
0
1
cos θe+e+
若上述过程平均事例数为ν，则观测到N个事例的 概率服从泊松分布
PP
=
e−νν N
N!
在这N个事例中，如果单个事例为前向的概率为f， 则观测到F个前向事例的概率满足二项式分布
=λ
举例：光电倍增管暗电流影响
在有11146根PMT的探测 器中，已知每根PMT暗电 流产生的误击中为3.5kHz。 求探测器在任意总长度为 500μs时间段观察到每隔 10ns PMT误击中数目分 别为5和6的总次数
在10ns间隔观测到PMT误 击中的平均数目为
…
得到(n1,n2,…,nm)概率为
f
G (n)
=
N! n1!n2!...nm!
p n1 1
p n2 2
...
p n3 m
平均值 : E[ni ] = Npi
方差 :V[ni ] = Npi (1− pi ) 协方差 :Vij = −Npi p j (i ≠ j)
适适用用于于直直方方图图 频频数数误误差差估估计计。。
MRPC记录的击中数目N’
MRPC探测效率 测量值及其误差
从二项式到多项式分布
类似于二项式分布,但允许结果的可能性m大于两种,概率为
G p = ( p1, p2..., pm )
m
∑ pi = 1
尝试N次，结果为
可能性1：n1 可能性2：n2
G n = (n1, n2 ,..., nm )
i =1
举例：角分布中的前后不对称
e+
e+
θ e-
+ e− → J /ψ → e+ + e− -1 B：后向计数；F：前向计数；N=B+F
0
1
cos θe+e+
若上述过程平均事例数为ν，则观测到N个事例的 概率服从泊松分布
PP
=
e−νν N
N!
在这N个事例中，如果单个事例为前向的概率为f， 则观测到F个前向事例的概率满足二项式分布
粒子物理与核物理实验中的
这个例子被广泛任为是“无辜偏向性”的受害者
在粒子物理与核物理实验研究中，历史上类似的“无辜偏向性”受害者 很多。为了避免这种情况的发生，近几年来盲分析方法渐渐成为潮流。
8
盲分析方法分类
信号区隐藏法
•最适合于对稀有或禁戒物理过程的实验测量； •本底必须从非信号区、或模拟样本、或额外的子实验来估计。
•对信号与本底事例数的估计既可以采用蒙特卡罗模拟，也可以直接利用 与信号有很多共同特征的本底数据样本。 •选择条件的调整与优化一定要尽量避免受到统计涨落的影响。
•为了验证选择条件，可将样本随机分为1/3与2/3两部分。 •判断选择条件和本底估计有否偏向性，可以通过比较1/3与2/3样本的 结果是否在统计上相符来进行。
•导致结果因此包含不可定量估计的系统误差。
4
实验结果的偏向性
中子寿命测量
K0S寿命测量
年代
在粒子物理与核物理实验中有不少实验结果显示了非常明显的时间相关 性，前一次实验与后一次实验要么在误差范围内，要么在几倍误差之外
5
并合不等精度实验结果时的困惑
实验结果与理论符合得太好。 但是χ2/dof 远小于1
13
例子：稀有K衰变中的本底
理论预言：Br(K + → π +νν ) = (0.79 ± 0.12) ×10−10
过程 所有Κ衰变 K+ → μ+νμ K+ → π+π0 K+ → μ+νμ γ 束流本底
电荷交换本底K+ n→ K0p, K0 → p+l-n
信号
事例数 1010
0.6343×1010 0.2113 ×1010 0.0055 ×1010
π − p → p + MM −
在粒子物理与核物理实验研究中，历史上类似的“无辜偏向性”受害者 很多。为了避免这种情况的发生，近几年来盲分析方法渐渐成为潮流。
8
盲分析方法分类
信号区隐藏法
•最适合于对稀有或禁戒物理过程的实验测量； •本底必须从非信号区、或模拟样本、或额外的子实验来估计。
•对信号与本底事例数的估计既可以采用蒙特卡罗模拟，也可以直接利用 与信号有很多共同特征的本底数据样本。 •选择条件的调整与优化一定要尽量避免受到统计涨落的影响。
•为了验证选择条件，可将样本随机分为1/3与2/3两部分。 •判断选择条件和本底估计有否偏向性，可以通过比较1/3与2/3样本的 结果是否在统计上相符来进行。
•导致结果因此包含不可定量估计的系统误差。
4
实验结果的偏向性
中子寿命测量
K0S寿命测量
年代
在粒子物理与核物理实验中有不少实验结果显示了非常明显的时间相关 性，前一次实验与后一次实验要么在误差范围内，要么在几倍误差之外
5
并合不等精度实验结果时的困惑
实验结果与理论符合得太好。 但是χ2/dof 远小于1
13
例子：稀有K衰变中的本底
理论预言：Br(K + → π +νν ) = (0.79 ± 0.12) ×10−10
过程 所有Κ衰变 K+ → μ+νμ K+ → π+π0 K+ → μ+νμ γ 束流本底
电荷交换本底K+ n→ K0p, K0 → p+l-n
信号
事例数 1010
0.6343×1010 0.2113 ×1010 0.0055 ×1010
π − p → p + MM −
粒子物理与核物理实验中的数据分析
10/04/2021
14
例子：对长寿命 K 介子的鉴别
强子量能器
h–
K
0 L
利用KL0粒子 不受磁场影 响而且较少
发生电磁簇
射的特点把
它和带电强
子区分开来。
电磁量能器
为常数，其余为实验观测量
10/04/2021
Eur.Phys.J.C10，1（1999）
把一个2-维甄别问题 简化为一维甄别问题。
通常情况下很难处理多维的
x
问题，
因此， 常常构造低维的统计检验,在
不失去甄别各种假设能力的条件下， 使得 t(x)成为精简后的数据样本。
那么此时的统计量 t 具有概率密度函数 g(t | H0 ), g(t | H1),...
10/04/2021
6
拒绝域、第一与第二类误差
考虑统计检验量t 服从 g(t | H0 ), g(t | H1),... g(t)
上一层节点函数可写为
n
hi (x) s(wi0 wij x j )
j 1
ai , wij为权重或者联结强度。
t(x) 输出定义为
n
t(x) s[a0 aihi (x)]
i 1
越多节点
神经网络越接
近优化的 t(x)
但需要定更多的参数!
10/04/2021
22
神经网络中的误差函数最小化
参数取值通常根据误差函数的最小化结果来决定
单元数为M n。
f (x | H0 ) f (x | H1)
但是如果 n 太大时，实 际运用会很 困难。
10/04/2021
12
例子:蒙特卡罗近似求二维p.d.f.
M.C.
M.C.
粒子物理与核物理实验中的数据分析
有A21的概率被测量为b2
真值为x2，有A12的概率被测量为b1
有A22的概率被测量为b2
归一化：A11 A21 1，A12 A22 1
(不考虑测量效率和本底)
实际的问题经常是：已知测量值b，希望求出真值x
如果可以知道A，则问题可以解决。
问题转化为：如何得到A？
9
响应矩阵与二维散点图
14
练习
1. 用舍选法，产生f(x)=3x2+2x的分布 (0<x<1)，并用自己写似然函数进行拟合。 注意归一化问题。
2. 实验测到数据存放在expdata.root中。试 读取之，画出直方图，估计可能的函数形 式，用最大似然法拟合之。
/~yangzw/Cour seDataAna/expdata.root
15
基于ROOT的解谱法(Unfolding) ROOT手册263页 /~adye/software/unf old/RooUnfold.html
ROOT中各种拟合
(前2个简单介绍，在project练习的过程中自己学习 使用，最后一个主要是本节练习)
3
12
矩阵求逆的解谱法
MC得到响应矩阵R 填入TMatrixD中，调用Invert()函数求逆
矩阵 利用逆矩阵以及测量值还原真实值
(得到的结果振荡很大，方差很大，需要引入 正规化函数进行平滑处理，减少方差。
但要注意的是：矩阵求逆方法是无偏和有效的， 平滑处理引入了偏置性，即系统误差)
13
1 1
,
〈0 〈1
1: 完美探测器
0 : 分辨率很差
Aˆ 1
1
2
粒子物理与核物理实验中数据分析
G4Element* H = new G4Element(name="Hydrogen",symbol="H" , z= 1., a); G4Element* O = new G4Element(name="Oxygen" ,symbol="O" , z= 8., a); density = 1.000*g/cm3; G4Material* H2O = new G4Material(name="Water", density, ncomponents=2); H2O->AddElement(H, natoms=2); H2O->AddElement(O, natoms=1); //定义水，给定密度、元素种类数目、添加元素
• 参考 资料 1) http://geant4.cern.ch 2)Nuclear Instruments and Methods in Physics Research A 506 (2003) 250-303, and IEEE Transactions on Nuclear Science 53 No. 1 (2006) 270-278. 最新版为9.1版，于2008年2月5日发布
2)下载安装CLHEP程序包(这是唯一需要预安装的程 序)
3)下载Geant4软件包以及相应的数据文件(用于各 种物理模型)，按照安装手册进行编译安装
如果系统版本相同(内核版本和g++版本)，把已经 编译好的程序直接复制到其它机器上即可使用。
比如，对SLC3系统，直接复制training服务器 /projects/soft/ext/clhep.tgz和g4.tgz到 本地机器，解压缩到相应目录即可。
粒子物理与核物理实验中的数据分析-第1讲-基本概念
24/02/2009
A
9
文恩图（Venn diagram）检验
A B
A B A B A ( A B) A ( A B) ( A B) A A B ( A B) ( A B) ( A B) A ( B C ) ( A B) ( A C )
2 Q1Q2 0.25Q0 exp( L / L0 ), L0 2 z ln(Q1 / Q2 )
24/02/2009 6
举例：测量闪烁体衰减长度（续）
2 Q1Q2 0.25Q0 exp( L / L0 ),
L0 2 z ln(Q1 / Q2 )
实验采用恒定光源，因此 Q0 为常数，对待测闪烁体 L0 也 为常数。理论上只要在给定一个位置 z，测量闪烁体两端的 电荷输出量即可。但在实际中，往往需要做多点测量。
也就是说，你可能没什么问题！? 从你的观点上看：对自己染上AIDS结果的可信度为3.2%。 从医生角度上看：象你这样的人有3.2%感染上了AIDS。 涉及到如何诠释结果（概率）的问题！
24/02/2009 16
概率含义的诠释
相对频率（频率论者）
假设A，B，…是一可重复实验的结果，则概率就是
结果为A P ( A) lim n n 次 实验
考虑任何一次AIDS检查的结果只有阴性(-)或阳性(+)两种
P( | AIDS) 0.98 P( | AIDS) 0.02 P( | no AIDS) 0.03 P( | no AIDS) 0.97 AIDS感 染患者阳性的 概率 AIDS感 染患者阴性的 概率 AIDS未 感 染 者阳性的概 率 AIDS未 感 染 者阴性的概 率
需要解决好 •A 的定义 •适当的误差
A
9
文恩图（Venn diagram）检验
A B
A B A B A ( A B) A ( A B) ( A B) A A B ( A B) ( A B) ( A B) A ( B C ) ( A B) ( A C )
2 Q1Q2 0.25Q0 exp( L / L0 ), L0 2 z ln(Q1 / Q2 )
24/02/2009 6
举例：测量闪烁体衰减长度（续）
2 Q1Q2 0.25Q0 exp( L / L0 ),
L0 2 z ln(Q1 / Q2 )
实验采用恒定光源，因此 Q0 为常数，对待测闪烁体 L0 也 为常数。理论上只要在给定一个位置 z，测量闪烁体两端的 电荷输出量即可。但在实际中，往往需要做多点测量。
也就是说，你可能没什么问题！? 从你的观点上看：对自己染上AIDS结果的可信度为3.2%。 从医生角度上看：象你这样的人有3.2%感染上了AIDS。 涉及到如何诠释结果（概率）的问题！
24/02/2009 16
概率含义的诠释
相对频率（频率论者）
假设A，B，…是一可重复实验的结果，则概率就是
结果为A P ( A) lim n n 次 实验
考虑任何一次AIDS检查的结果只有阴性(-)或阳性(+)两种
P( | AIDS) 0.98 P( | AIDS) 0.02 P( | no AIDS) 0.03 P( | no AIDS) 0.97 AIDS感 染患者阳性的 概率 AIDS感 染患者阴性的 概率 AIDS未 感 染 者阳性的概 率 AIDS未 感 染 者阴性的概 率
需要解决好 •A 的定义 •适当的误差
粒子物理与核物理实验中的数据分析
对于二维情形，其概率密度函数可表示为
f
( x1 ,
x2; 1,
2 ,1, 2 ,
)
2
1
1 2
1
2
cov[x1, x2 ] /(1 2 )
exp
2(1
1
2
)
x1 1 1
2
x2 2 2
2
我们为大亚湾实验研制生产触发电子学板。按设计在一年内
需要修理的电路板为10%。如果在实验所需的20块板中有 5块在第一年使用时需要进行维修，那么这种故障率是否可 以接受？
解答：首先找出在一年内20块板中有5块或更多出现问题需 要进行维修的概率
20
4
f (n; 20, 0.1) 1 f (n; 20, 0.1)
07/11/2019
8
泊松分布式是二项分布的近似
概率的第三公理：如果 A1，A2，A3，… 是在空间 S 中互斥 事例的有限或无限序列，则
P( A1 A2 A3 ...) P( A1 ) P( A2 ) P( A3 ) ...
e n e n e e 1 (函数为 ex 的麦克劳林公式 )
i , V[ y]
2 i
i 1
i 1
i 1
07/11/2019
15
高斯分布与泊松分布
Probability
=2
泊松分布 高斯分布 (=,=)
=5
=10
r successes (or failures)
07/11/2019
泊松分布只有非负 整数定义。 高斯分布是连续且 可延伸到正负无穷。 当泊松分布的平均 值越大，与高斯分布 的区别就越小。 实际应用时，当计 数或事例数大于5 时， 可认为误差满足高斯 分布。
粒子物理与核物理试验中的数据分析
们感兴趣的另一随机序列 x1 , x2 ,..., xn。 3) 利用这些 x 值来估计 f ( x) 的一些特性，例如：通过找 b 到在区间 a x b 的 x 比例，给出积分值 a f ( x)dx 。 第一层面上的应用： 蒙特卡罗计算 = 积分
第二层面上的应用：
12/11/2018
蒙特卡罗变量 = “模拟的数据”
4
随机数的产生
用物理方法产生 真正的随机数 用数学方法产生 伪随机数
•不可重复 •产生速度慢
•可以重复 •产生的速度快
12/11/2018
5
真随机数与伪随机数
美国兰德（RAND）公司在1950年代，利用真空管中 产生的噪音制作了一个含十万个真正的随机数表，并运用于 其开展的所有模拟研究中。
真正的随机数与伪随机数之间的区别在于：数据串是否 具有可压缩性，即能否用更短的形式来表示。 真正的随机数是不可压缩的，非常不规则，以至于无法 用更短的形式来表示它。
粒子物理与核物理实验中的数 据分析
杨振伟 清华大学
第四讲:蒙特卡罗方法
12/11/2018 1
上一讲回顾
概率的基本概念 随机变量与概率密度函数 随机变量的平均值与方差
能不通过实验对随机变量进行研究吗？
12/11/2018
2
本讲要点
蒙特卡罗方法 随机数产生子
任意分布抽样之函数变换法与舍选法
数值解： a
b
函数必须解析可积 自变量不能太多
f ( x)dx f ( xi )
i 1
n
ba n
蒙特卡罗方法： A
在AB区间均匀投总数为N个点。
f ( x)
b
粒子物理与核物理实验中的数据分析
1
t
1
L
•实验中观测量的分布
F ( L′; < l >, k ) = α i f ( L′; < l >, k ) + (1 − α ) BG本底 ( L′) f ( L′; < l >, k ) =
+∞
−∞
∫ P( L; < l >, k )iGauss( L − L′, σ
L
)dL
拟合参数： l >, k < α = 信号所占的比率(由蒙特卡罗确定) σ L = k iσ (蒙特卡罗确定σ ，修正量k由拟合确定) 18
一般情况下检查区别
如果两种分析样本不完全重叠 引入加权平均 方差
2 2 σ Δ = σ 12 + σ 2 − 2 ρσ 1σ 2
a ( w) = wa1 + (1 − w)a2
2 ∂σ a ( w)
2 2 σ a ( w) = w2σ 12 + (1 − w) 2 σ 2 + 2w(1 − w) ρσ 1σ 2
不确定性 错误
•统计分析还提供了如何鉴别一个错误的方法，但它不能告诉我们下一步 该如何做，因为它无法告诉我们错误的根源在哪里。
5
从语义学上定义系统误差
•物理学家通常将随机(统计)误差定义为随机不确定性而不是随机的错误 •为了与上述定义保持一致，应该将系统误差定义为系统不确定性而不是 系统的错误
systematic error = systematic uncertainty ≠ systematic mistake
N ee = ∫ Liσ ee dt = σ ee i ∫ Ldt
如果理论计算精度只到第三阶 亮度计算结果总是给出同样的不准确性。
粒子物理与核物理实验中的数据分析
0
2 2 T 0
设经过 n 次迭代以后，找到一组解 x(n) , (n) ，得到函数
2 ( x(n) , (n) ) 的值。在 x(n)上对 (n) 进行线性展开，并略
去高阶项，得到
(n) Fx(n) ( x(n1) x(n) ) 0
而 n+1 次迭代后
ˆ 没有偏置，而且得到的方差
最小(高斯-马尔可夫定理)。
这里 aj(x) 是 x 的任意线性独立函数。 用矩阵来表示时，令 Aij=aj(xi)，有
2 ( ) ( y )TV 1( y ) ( y A )TV 1( y A )
对 i 求偏微分，并令结果等于零，有
E[ yi ] i (xi; )
这里，x1 , ..., xN
与
V
[
yi
]
2 i
已知。
为了估计参数 ，可以用曲线
拟合所有的测量点(见右图)。
对于独立的高斯变量 yi，联合概率密度函数为
N
g( y; , 2 )
i 1
1
2
2 i
exp
( yi i )2
如果 yi 是多维高斯变量，协方差矩阵为V ，满足
g(
y;
,V
)
(2
)N
1
/2
|V
|1/ 2
exp
1 2
y T V 1
y
那么其对数似然函数为
log
L(
)
1 2
N i 1
[ yi
( xi ;
)](V
粒子物理与核物理实验中的数据分析00002
根据概率的相对频率定义，在 n 次测量中出现 ti 频率为一次
P(ti
)
1 n
因此，期待值（或平均寿命）为
E[t]
n
ti
i 1
1 n
1 n
n
ti
i 1
思考：如果频率为 mi 次，结果会不同吗？
26/08/2020
21
误差传递
假设
x
( x1 ,...,
xn
) 服从某一联合
p.d.f.
f
( x ) ，我们也许并不
结论：只有1）与4）是合理的。
评论：作为一个合格的实验研究人员，一定要具备判断 结果是否合理的能力！
26/08/2020
11
举例：检查经验概率密度函数
实验上经常经验性地从直方图中给出概率密度函数（例如 通过拟合直方图分布等等），但是需要确定得到的函数是否 满足概率密度函数的定义，例如
1) f (x) x 2 2
j 1
y xi
y x j
Vij
x
26/08/2020
23
误差传递(续二)
两项合起来给出 y(x) 的方差
2 y
V[y]
n y
i, j1 xi
y x j
Vij x
如果 xi 之间是无关的，则Vij i2ij ，那么上式变为
2 y
V[y]
n i 1
2
y xi
x
2 i
证明：由于 A 与 A 根据定义是互斥的，并且从文恩图得到
A A S
因此可以写出
P( A) P( A) P( A A) P(S) 1
P( A) 1 P( A)
26/08/2020
粒子物理和核物理实验中的数据分析
＝
➢将观测直方图乘以修正因子直方图得到理论 (真实)的直方图
＝
20
正规化的解谱法
考虑“合理的”估计量，对选定的logL 满足
r
l o g L () l o g L m a x l o g L
logL描 述 了 数 据 n r与 期
待 值 r之 间 的 “ 距 离 ” 。
r 估计量满足该不等式并且最光滑，等价于 将下式求最大值
这个随机 1/数 放被 大后，结果信 中息 有完 用全 的被非湮 物没 理。 振
所以，通常情接 况求 下反 ，应 直矩阵x的 获方 得法 真， 值尽管理 是严格的，是的 无， 偏但 有结 效果并意 没义 有。 物理
解决办法是进行平滑处理，消除无意义的统计涨落。 但平滑会带来偏向性，需要在涨落与偏向性之间找到平衡。 13
通常取 k=2，使得 S 约等于曲率平方的平均 值。对直方图而言，也就是
r M 2
S( ) (i
2i1i2)2
Sov.
Math.5(1963)1035
i1
注意：2 阶导数对直方图的第一和最后的区 间没有很好的定义。
24
Tikhonov 规则（续）
如果在
log L 1 2
2
下，采用Tikhonov(k=2)规则
注意：,是常数，n会受
到统计涨落的影响。
真实直方图 离散化的p.d.f.
观测数据 和数据的期待值
7
效率、本底
有时，事例可能会不被探测到: 效率
N
N
Rij P(观测值在第i 区|真实值在第j 区)
i1
i1
真实直方图第
P(观测值在全范围|真实值在第j区) j 区探测效率
j(效率)
CERN ROOT-粒子物理与核物理实验中的数据分析-第二讲.
如果还有其它额外的信息，应该给出不同的先
验概率。这种贝叶斯统计的特点必定是主观的。例 如，受检者有过吸毒历史。一旦验前概率改变，贝 叶斯定理就会告诉患病的可能性。对阳性结果的诠 释就会改变。
问题：能否构造含自变量的概率？
26/09/2020
2
随机变量与概率密度函数
假设实验结果为 x (记作样本空间中元素)的概率为
cx , o y ] E [ v x ( x ) [ y ( y ) E ] [ x ] y xy
相关系数定义为
xy
covx,[y],
xy
1xy1
如果 x，y 独立，即
f(x,y)fx(x)fy(y)
则
covx,[y]0
26/09/2020
23
举例：样本平均值
假设实验上研究一核素衰变寿命，在探测效率为100%的情况 下，每次探测到的寿命为 ti，一共测量了 n 次，求平均寿命 （也就是寿命的期待值）。
g(a)da f(x)dx dS dS a在[a,ada]内的x空间范围
x(ada)
g(a)da
f (x)dx
x(a)
x(a)dx da
da f(x)dx
x(a)
g(a) f(x(a)) dx da
26/09/2020
17
函数的逆不唯一情况
假如 a(x) 的逆不唯一，则函数的 p.d.f. 应将 dS 中对应于 da 的所有 dx 的区间包括进来
下列各种情况给出的概率值是否是合理的：
1 ) P (A )1/3 ,P (B )1/3 ,P (C )1/3 2) P (A )0.64 ,P (B )0.38 ,P (C ) 0.02 3 ) P (A )0.35 ,P (B )0.52 ,P (C )0.26 4) P (A )0.57,P (B )0.24 ,P (C )0.19
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
G4HCofThisEvent* hc = evt->GetHCofThisEvent(); G4int NbOfColl = hc->GetNumberOfCollections(); //获得第0个HitsCollection，即ExN02TrackerHitsCollection //也可以通过CollectionID获得 ExN02TrackerHitsCollection *hitsC = hc->GetHC(0); //该Collection中Hits数目 G4int sizehits = hitsC->entries(); .......
产生事例：G4ParticleGun 物理过程：电磁、强作用、衰变、光轻子-强 子作用、光学、参数化、输运(必要过程)
2015-7-8 2
本讲要点
产生主事例 G4HEPEvtInterface 灵敏探测器 取出灵敏探测器的数据，并存入ROOT 格式文件
2015-7-8
3
事例产生子接口
G4HEPEvtInterface 很多时候，事例产生子已经存在，而且是Fortran 语言。Geant4并不直接链接这些Fortran程序， 而是提供了一个接口： G4HEPEvtInterface读取事例产生子生成的 ASCII文件中的信息，重新生成 G4PrimaryParticle对象，并关联到对应的 G4PrimaryVertex 也就是说，G4HEPEvtInterface将/HEPEVT/公共 块的信息转换为一个O-O数据结构。这个公共块 在高能物理中被广泛使用。
粒子物理与核物理实验中的数 据分析
杨振伟 清华大学 第八讲:Geant4 的探测 器模拟介绍(3)
2015-7-8 1
上讲回顾
粒子定义
G4ParticleDefinition 6大类粒子：G4LeptonConstructor G4BosonConstructor G4MesonConstructor G4BaryonConstructor G4IonConstructor G4ShortlivedConstructor
G4SDManager* SDman = G4SDManager::GetSDMpointer(); G4String trackerChamberSDname = "ExN02/TrackerChamberSD"; ExN02TrackerSD* aTrackerSD = new ExN02TrackerSD( trackerChamberSDname );
5
以HEPEVT格式输出的ASCII文件
比如：下面这个事例表示该事例共102个粒子(包括中间态)， 随后的102行分别为这102个粒子的具体信息： 第一列为粒子状态(3：对撞入射粒子或其它；2：衰变了；1： 存在的粒子；0：空)， 第2列为粒子PDG号， 最后4列分别为粒子的x，y，z方向动量和质量。
用户灵敏探测器继承自抽象基类G4VSensitiveDetector，用户需要完成 3个主要函数： ProcessHits(G4Step* aStep, G4TouchableHistory*) 构造“击中”，被G4SteppingManager调用 Initialize(G4HCofThisEvent* HCE) 初始化，事例开始时调用，指定构造的“击中”与当前事例关联起来 EndOfEvent(G4HCofThisEvent*) 事例结束时调用
参见例子N02/src/
2015-7-8 8
定义和添加灵敏探测器(1)
1.定义Hits，如 2.定义SD，如 3.在DetectorConstruction()中添加SD 在探测器构造中添加敏感探测器，比如： //SDManager //创建敏感探测器
175
201Байду номын сангаас-7-8
6
使用HEPEvtInterface的例子
参见例子N04，在中：
ExN04PrimaryGeneratorAction::ExN04PrimaryGeneratorAction() { const char* filename = "pythia_event.data"; //读取pythia_event.data HEPEvt = new G4HEPEvtInterface(filename); } void ExN04PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) { //设定主顶点位置，产生主顶点 HEPEvt->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,0.*cm); HEPEvt->GeneratePrimaryVertex(anEvent); } 其中HEPEvt在头文件中定义： G4VPrimaryGenerator* HEPEvt;
2015-7-8 9
定义和添加灵敏探测器(2)
将多个logical体积添加为 灵敏探测器时： V1 假设有3个体积V1，V2， V3 如果定义这3个体积的时 候，先定义V1，再定义V2， 最后定义V3，则V1，V2， 各自被覆盖掉一部分。
V2
V3
如果希望蓝色区域为SD，则需要 V2->SetSensitiveDetector(....)
2015-7-8 13
练习
在例子N02的基础上，将模拟的信息存储 到root文件中。这些信息包括：粒子的 PDG号、质量、能量沉积、径迹长度。生 成root文件后画出这些信息的直方图，并 进行分析 修改探测器物质和入射粒子，重新运行， 得到新的root文件，并画出储存信息的直 方图。 在N03的基础上，加入敏感探测器。
2015-7-8 4
用/HEPEVT/公共块生成ASCII文件
common block
将以下量写入文件中
2015-7-8
第一行：NHEP，当前事例粒子数(包括中间态) 随后的NHEP行：每个粒子的ISTHEP,IDHEP,JDAHEP,PHEP信息 ISTHEP：粒子状态；IDHEP：粒子PDG号；JDAHEP：粒子衰变产物 位置的指针；PHEP(1-3,5)：粒子x,y,z动量，能量，质量
2015-7-8 14
参考资料
1. 2. 3.
Geant4应用开发手册3.6节 Geant4应用开发手册4.4节 Geant4例子novice/N02,N04
2015-7-8
15
2015-7-8 10
读取敏感探测器的信息
在EventAction类的EndOfEventAction()函数中，可以读取 该事例中存储的Hits。比如可以在中 加入下面代码，查看每个事例中的Hits数目： //获得该事例的HitsCollection(可能不止一个)
参见
/~yangzw/CourseDataAna/examples/Lec8.tgz
2015-7-8 12
小结
G4HEPEvtInterface 主产生子(PrimaryGenerator)的一种，直接 读取ASCII文件中以HEPEVT格式存储的事例。 敏感探测器的添加和定义 在DetectorConstruction中，不但要将SD 添加给SDManager，还要指定相应的logical 体积。 将结果存储到root文件中 在EventAction中收集数据，或者在SD中直接 收集。
//添加到SDManager
SDman->AddNewDetector( aTrackerSD );
//为logical体积设定敏感探测器!!!
logicChamber->SetSensitiveDetector( aTrackerSD ); 参见例子N02/src/
2015-7-8 注：main 函数或者mac文件中设定beamOn事例数不能超过ASCII中事例数。 7
灵敏探测器(Sensitive Detector)
灵敏探测器(SD)的首要任务是通过粒子“迹”(track)上的 “步”(step)的信息，构造“击中”(hit)。 这些击中经过数字化，被读出模块读出的信息是真正的模拟结果。(当 然在模拟中我们也可以忽略数字化而直接读出hit的信息或者其它信息， 这些信息实际上是所谓的"Monte Carlo Truth")
102 3 3 3 3 3 3 3 3 3 2 1 1 ...... ...... 11 -11 11 -11 11 -11 23 2 -2 23 22 22 0 0 0 0 0 0 0 0 0 16 0 0 0 0.00000000E+00 0.00000000E+00 0.25000000E+03 0.51000000E-03 0 0.00000000E+00 0.00000000E+00 -0.25000000E+03 0.51000000E-03 0 0.00000000E+00 0.00000000E+00 0.24999999E+03 0.00000000E+00 0 0.00000000E+00 0.00000000E+00 -0.25000000E+03 0.00000000E+00 0 0.37396914E-02 0.15234913E-02 0.24138585E+03 0.00000000E+00 0 -0.93164320E-02 0.27396574E-01 -0.24687934E+03 0.00000000E+00 0 -0.55767406E-02 0.28920065E-01 -0.54934906E+01 0.48823428E+03 0 0.19070032E+02 0.24337596E+03 -0.48627266E+01 0.33000000E+00 0 -0.19075609E+02 -0.24334704E+03 -0.63076405E+00 0.33000000E+00 26 -0.55767406E-02 0.28920065E-01 -0.54934906E+01 0.48823428E+03 0 0.93164331E-02 -0.27396573E-01 -0.31205891E+01 0.00000000E+00 0 -0.81046576E-03 -0.82301151E-04 0.14162632E+00 0.00000000E+00
产生事例：G4ParticleGun 物理过程：电磁、强作用、衰变、光轻子-强 子作用、光学、参数化、输运(必要过程)
2015-7-8 2
本讲要点
产生主事例 G4HEPEvtInterface 灵敏探测器 取出灵敏探测器的数据，并存入ROOT 格式文件
2015-7-8
3
事例产生子接口
G4HEPEvtInterface 很多时候，事例产生子已经存在，而且是Fortran 语言。Geant4并不直接链接这些Fortran程序， 而是提供了一个接口： G4HEPEvtInterface读取事例产生子生成的 ASCII文件中的信息，重新生成 G4PrimaryParticle对象，并关联到对应的 G4PrimaryVertex 也就是说，G4HEPEvtInterface将/HEPEVT/公共 块的信息转换为一个O-O数据结构。这个公共块 在高能物理中被广泛使用。
粒子物理与核物理实验中的数 据分析
杨振伟 清华大学 第八讲:Geant4 的探测 器模拟介绍(3)
2015-7-8 1
上讲回顾
粒子定义
G4ParticleDefinition 6大类粒子：G4LeptonConstructor G4BosonConstructor G4MesonConstructor G4BaryonConstructor G4IonConstructor G4ShortlivedConstructor
G4SDManager* SDman = G4SDManager::GetSDMpointer(); G4String trackerChamberSDname = "ExN02/TrackerChamberSD"; ExN02TrackerSD* aTrackerSD = new ExN02TrackerSD( trackerChamberSDname );
5
以HEPEVT格式输出的ASCII文件
比如：下面这个事例表示该事例共102个粒子(包括中间态)， 随后的102行分别为这102个粒子的具体信息： 第一列为粒子状态(3：对撞入射粒子或其它；2：衰变了；1： 存在的粒子；0：空)， 第2列为粒子PDG号， 最后4列分别为粒子的x，y，z方向动量和质量。
用户灵敏探测器继承自抽象基类G4VSensitiveDetector，用户需要完成 3个主要函数： ProcessHits(G4Step* aStep, G4TouchableHistory*) 构造“击中”，被G4SteppingManager调用 Initialize(G4HCofThisEvent* HCE) 初始化，事例开始时调用，指定构造的“击中”与当前事例关联起来 EndOfEvent(G4HCofThisEvent*) 事例结束时调用
参见例子N02/src/
2015-7-8 8
定义和添加灵敏探测器(1)
1.定义Hits，如 2.定义SD，如 3.在DetectorConstruction()中添加SD 在探测器构造中添加敏感探测器，比如： //SDManager //创建敏感探测器
175
201Байду номын сангаас-7-8
6
使用HEPEvtInterface的例子
参见例子N04，在中：
ExN04PrimaryGeneratorAction::ExN04PrimaryGeneratorAction() { const char* filename = "pythia_event.data"; //读取pythia_event.data HEPEvt = new G4HEPEvtInterface(filename); } void ExN04PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) { //设定主顶点位置，产生主顶点 HEPEvt->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,0.*cm); HEPEvt->GeneratePrimaryVertex(anEvent); } 其中HEPEvt在头文件中定义： G4VPrimaryGenerator* HEPEvt;
2015-7-8 9
定义和添加灵敏探测器(2)
将多个logical体积添加为 灵敏探测器时： V1 假设有3个体积V1，V2， V3 如果定义这3个体积的时 候，先定义V1，再定义V2， 最后定义V3，则V1，V2， 各自被覆盖掉一部分。
V2
V3
如果希望蓝色区域为SD，则需要 V2->SetSensitiveDetector(....)
2015-7-8 13
练习
在例子N02的基础上，将模拟的信息存储 到root文件中。这些信息包括：粒子的 PDG号、质量、能量沉积、径迹长度。生 成root文件后画出这些信息的直方图，并 进行分析 修改探测器物质和入射粒子，重新运行， 得到新的root文件，并画出储存信息的直 方图。 在N03的基础上，加入敏感探测器。
2015-7-8 4
用/HEPEVT/公共块生成ASCII文件
common block
将以下量写入文件中
2015-7-8
第一行：NHEP，当前事例粒子数(包括中间态) 随后的NHEP行：每个粒子的ISTHEP,IDHEP,JDAHEP,PHEP信息 ISTHEP：粒子状态；IDHEP：粒子PDG号；JDAHEP：粒子衰变产物 位置的指针；PHEP(1-3,5)：粒子x,y,z动量，能量，质量
2015-7-8 14
参考资料
1. 2. 3.
Geant4应用开发手册3.6节 Geant4应用开发手册4.4节 Geant4例子novice/N02,N04
2015-7-8
15
2015-7-8 10
读取敏感探测器的信息
在EventAction类的EndOfEventAction()函数中，可以读取 该事例中存储的Hits。比如可以在中 加入下面代码，查看每个事例中的Hits数目： //获得该事例的HitsCollection(可能不止一个)
参见
/~yangzw/CourseDataAna/examples/Lec8.tgz
2015-7-8 12
小结
G4HEPEvtInterface 主产生子(PrimaryGenerator)的一种，直接 读取ASCII文件中以HEPEVT格式存储的事例。 敏感探测器的添加和定义 在DetectorConstruction中，不但要将SD 添加给SDManager，还要指定相应的logical 体积。 将结果存储到root文件中 在EventAction中收集数据，或者在SD中直接 收集。
//添加到SDManager
SDman->AddNewDetector( aTrackerSD );
//为logical体积设定敏感探测器!!!
logicChamber->SetSensitiveDetector( aTrackerSD ); 参见例子N02/src/
2015-7-8 注：main 函数或者mac文件中设定beamOn事例数不能超过ASCII中事例数。 7
灵敏探测器(Sensitive Detector)
灵敏探测器(SD)的首要任务是通过粒子“迹”(track)上的 “步”(step)的信息，构造“击中”(hit)。 这些击中经过数字化，被读出模块读出的信息是真正的模拟结果。(当 然在模拟中我们也可以忽略数字化而直接读出hit的信息或者其它信息， 这些信息实际上是所谓的"Monte Carlo Truth")
102 3 3 3 3 3 3 3 3 3 2 1 1 ...... ...... 11 -11 11 -11 11 -11 23 2 -2 23 22 22 0 0 0 0 0 0 0 0 0 16 0 0 0 0.00000000E+00 0.00000000E+00 0.25000000E+03 0.51000000E-03 0 0.00000000E+00 0.00000000E+00 -0.25000000E+03 0.51000000E-03 0 0.00000000E+00 0.00000000E+00 0.24999999E+03 0.00000000E+00 0 0.00000000E+00 0.00000000E+00 -0.25000000E+03 0.00000000E+00 0 0.37396914E-02 0.15234913E-02 0.24138585E+03 0.00000000E+00 0 -0.93164320E-02 0.27396574E-01 -0.24687934E+03 0.00000000E+00 0 -0.55767406E-02 0.28920065E-01 -0.54934906E+01 0.48823428E+03 0 0.19070032E+02 0.24337596E+03 -0.48627266E+01 0.33000000E+00 0 -0.19075609E+02 -0.24334704E+03 -0.63076405E+00 0.33000000E+00 26 -0.55767406E-02 0.28920065E-01 -0.54934906E+01 0.48823428E+03 0 0.93164331E-02 -0.27396573E-01 -0.31205891E+01 0.00000000E+00 0 -0.81046576E-03 -0.82301151E-04 0.14162632E+00 0.00000000E+00