粒子物理与核物理试验中的数据分析课件.ppt

1、6/23/20221粒子物理与核物理实验中的数粒子物理与核物理实验中的数据分析据分析杨振伟清华大学第第四四讲讲:蒙特卡罗方法蒙特卡罗方法6/23/20222上一讲回顾上一讲回顾概率的基本概念随机变量与概率密度函数随机变量的平均值与方差能不通过实验对随机变量进行研究吗？6/23/20223本讲要点本讲要点n蒙特卡罗方法n随机数产生子n任意分布抽样之函数变换法与舍选法n蒙特卡罗方法中的精度问题n在粒子物理与核物理中的应用6/23/20224蒙特卡罗方法简介蒙特卡罗方法简介蒙特卡罗方法蒙特卡罗方法就是利用一系列随机数来计算各种概率大小和就是利用一系列随机数来计算各种概率大小和随机变量均值等等的数值分

2、析技术随机变量均值等等的数值分析技术。通常的步骤为：。通常的步骤为：1)产生一系列在0，1之间均匀分布的随机数 。 2)利用这些随机数按某些概率密度函数 抽样生成我们感兴趣的另一随机序列 。3)利用这些 值来估计 的一些特性，例如：通过找到在区间 的 比例，给出积分值 。mrrr,.,21)(xfnxxx,.,21xx)(xfbxabadxxf)(第一层面上的应用第一层面上的应用：蒙特卡罗计算蒙特卡罗计算 = 积分积分 第二层面上的应用第二层面上的应用：蒙特卡罗蒙特卡罗变量变量 = “模拟的数据模拟的数据” 6/23/20225随机数的产生随机数的产生用物理方法产生用物理方法产生真正的真正的随

3、机数随机数不可重复不可重复产生速度慢产生速度慢用数学方法产生用数学方法产生伪伪随机数随机数可以重复可以重复产生的速度快产生的速度快6/23/20226真随机数与伪随机数真随机数与伪随机数 美国兰德（RAND）公司在1950年代，利用真空管中产生的噪音制作了一个含十万个真正的随机数表，并运用于其开展的所有模拟研究中。 真正的随机数与伪随机数之间的区别在于：数据串是否具有可压缩性，即能否用更短的形式来表示。 真正的随机数是不可压缩的，非常不规则，以至于无法用更短的形式来表示它。 在粒子物理与核物理研究中，随机数的可重复性经常也是非常有用的，尤其是程序的调试(debugging)。6/23/2022

4、7随机数产生子随机数产生子目的目的是使在是使在 0，1 范围内产生的范围内产生的伪随机数伪随机数满足：满足：均匀性；相互独立性；长周期性乘同余法乘同余法110() , / iiiiMrMM余数随机数与为选定常数为随机数种子。友情推荐友情推荐M=2K=52q+10周期=2K-22325131230 1092365131234 210102425171240 10126/23/20228CERN库的随机数产生子库的随机数产生子PAWPAW用户用户 gRandom-SetSeed(); Float_t random = gRandom-Rndm(1); Real random(1)Call Rmar

5、in(ISEED,0,0)Call Ranmar(random,1)注意注意: 用于产生子的用于产生子的随机数种子还可以用随机数种子还可以用来保证后续进程的随来保证后续进程的随机数不重复。机数不重复。Root Root 用户用户 粒子物理与核物理研究中粒子物理与核物理研究中，大都采用，大都采用CERN程序库提供的随机数产生子。程序库提供的随机数产生子。6/23/20229随机数均匀性与相关性检验随机数均匀性与相关性检验 subroutine mc double precision lamda,M,x,x0,y call hbook1(10,r,100,0.,1.,0.) call hbook2

6、(20,r(i+1) vs. r(i),&100,0.,1.,100,0.,1.,0.) x0=1. lamda=1220703125 ! 5*13 M=4294967296. ! 2*32 do i=1,10000 x=Mod(lamda*x0,M) y=x/M call hfill(10,real(y),0.,1.0) if(i.gt.1)call & hfill(20,real(y_old),real(y),1.0) x0=x y_old=y end do return end随机变量随机变量第第 I 个个随机变量随机变量第第 I+1 个个随机变量随机变量频数频数均匀性均匀性相关性相关性

7、20 h/pl 10; h/pl 2; 1 zonePAWmc.f callPAW6/23/202210随机数均匀性与相关性检验随机数均匀性与相关性检验随机变量随机变量第第 I 个个随机变量随机变量第第 I+1 个个随机变量随机变量频数频数均匀性均匀性相关性相关性Draw()h2- root2Draw()h1- root1random.C.x root0 void random() UInt_t lambda,M,x0; TH1F *h1 = new TH1F(h1,100,0,1); TH2F *h2 = new TH2F(h2,100,0,1,100,0,1); lambda=122070

8、3125; /513 M =4294967296; /232 x0 =1; double y, y_old; for (int i=0; iFill(y); if (i1) h2-Fill(y_old,y); x0=x; y_old=y; 6/23/202211用蒙特卡罗法计算积分用蒙特卡罗法计算积分对于计算积分值( )baf x dx解析解：( )( )( )bx bx aaf x dxF xF x数值解：1( )( )nbiaibaf x dxf xn蒙特卡罗方法：abABOx( )f x函数必须解析可积自变量不能太多对函数是否解析可积和是否太多自变量无要求 in ( )( )() ()b

9、aNf xf x dxAOBON在AB区间均匀投总数为N个点。6/23/202212蒙特卡罗方法中的精度问题蒙特卡罗方法中的精度问题采用蒙特卡罗方法采用蒙特卡罗方法(MC)计算积分计算积分与传统的梯形法相比有如下特点与传统的梯形法相比有如下特点一维积分一维积分:) ( ,1/ :) ( ,1/ : MC2为子区间的数目梯形法精度为产生的随机数精度nnnn多维积分多维积分:) ( ,1/ :) ( , 1/ : MC/2为子区间的数目梯形法精度为产生的随机数且与维数无关精度维数nnnn对于维数大于对于维数大于4的积分的积分, 用蒙特卡罗方计算积分总是最好。用蒙特卡罗方计算积分总是最好。 6/23

10、/202213从均匀分布到任意分布的随机数从均匀分布到任意分布的随机数函数变换法舍选法寻找某个函数，当函数的自变量取均匀分布值时，对应的函数值自动满足给定分布。均匀分布给定分布从一个随机变量与对应概率密度函数最大值构成的二维均匀分布中，按概率密度函数与自变量关系曲线切割得到。6/23/202214函数变换法函数变换法nxxxxfrx,., )( )( , 1 0, 21分布的随机数找出服从通过适当的变换均匀分布的随机数从在) ()() ()( ，)() ( :rrxrxFdxxfrdrrgrxxPrrP即要求)( )(1rFxrxF解出令均匀分布均匀分布任意分布任意分布6/23/202215例

11、子例子:指数分布抽样指数分布抽样)1log()( ),( 1 )0( 1),( :)(0/ /rrxrxrdxexexfrxxx得到并解出令指数概率密度函数2呈一一对应关系与随机变量)( rxr抽样效率抽样效率为为100%。可采用函数变换法抽样的分布可采用函数变换法抽样的分布n指数分布指数分布n三维各向同性分布三维各向同性分布n二维随机角度的正、余弦分布二维随机角度的正、余弦分布n高斯分布高斯分布nn 个自由度的个自由度的 2 分布分布n伽马分布伽马分布n二项式分布二项式分布n泊松分布泊松分布nStudent 分布分布（http:/pdg.lbl.gov/2008/reviews/monter

12、pp.pdf）6/23/2022166/23/202217舍选法舍选法)(: , , )( minmax1minxxrxxxxxf即变量产生均匀分布的由第一个随机取值范围的自变量根据概率密度函数max2max , 0 ,fruf即范围内与均匀分布在变量产生第二个独立的随机 。， ，, ),( 新进行抽样从值拒绝该否则值则接受该如果xxxfu 问题问题: 如何找到函数的最大值如何找到函数的最大值?x的最大值)(xf)(xf6/23/202218舍选法举例舍选法举例subroutine acc_rejreal rvec(1)call hbook1(10,x(r),100,0.,10.,0.)cal

13、l hbook1(20,x(r),100,0.,10.,0.)call hbook2(30,f(x) vs. x(r),100,0.,10.,100,0.,1.1,0.)fmax=-999.do i=1,100 call ranmar(rvec,1) r=0+rvec(1)*(10.-0.) f=0.5*exp(-r/2.) if(fmax.lt.f)fmax=fend dofmax=1.2*fmaxntot=0do i=1,10000 call ranmar(rvec,1) r=0+rvec(1)*(10.-0.) z=0.5*exp(-r/2.) if(z.gt.fmax)then fma

14、x=z*1.2 write(6,*)z greater than fmax end if call hfill(10,r,0.,1.0) call ranmar(rvec,1) u=rvec(1)*fmax if(u.lt.z)then call hfill(20,r,0.,1.0) call hfill(30,r,u,1.0) ntot=ntot+1 end ifend dowrite(6,*)ntot=,ntotreturnend的最大值找 )( xf30 h/pl 20; h/pl 10; h/pl 3; 1 zonePAWacc_rej.f callPAW6/23/202219舍选法举

15、例舍选法举例void acc_rej() TH1F *h11 = new TH1F(h11,100,0,10); TH1F *h12 = new TH1F(h12,100,0,10); TH2F *h2 = new TH2F(h2,100,0,10,100,0,1); double fmax=-999.; for (int i=0;iUniform(0,10); double f = 0.5*exp(-r/2.); if (fmax f) fmax=f; fmax *= 1.2; cout fmax= fmax endl; int ntot=0; for (int i=0;iUniform(0

16、,10); double z=0.5*exp(-r/2.); if (zfmax) fmax=1.2*z; h11-Fill(r); double u=gRandom-Uniform(0,fmax); if (uFill(r); h2-Fill(r,u); ntot += 1; cout ntot= ntot Reset(); hx = new TH1F(hx,“x dis., 100,-10,10); gRandom-SetSeed(); Double_t x; const Double_t sigma=2.0; const Double_t mean=1.0; const Int_t kU

17、PDATE = 1000; for ( Int_t i=0; iGaus(mean,sigma); hx-Fill(x); 产生平均值为mean标准偏差为sigma的高斯分布。可以换为x = gRandom-Rndm(i); x = gRandom-Uniform(xup); x = gRandom-Integer(Imax); x = gRandom-Landau(mean,sigma); x = gRandom-Binomial(ntot,prob); x = gRandom-Poisson(mean); x = gRandom-PoissonD(mean); x = gRandom-Ex

18、p(tau); x = gRandom-BreitWigner (me,sig); 在ROOT环境下采用已有的分布，可以容易完成布置的练习。6/23/202223蒙特卡罗统计检验蒙特卡罗统计检验例如常用来检验理论与实验符合好坏的例如常用来检验理论与实验符合好坏的 2 分布。分布。四个服从 N(0,1) 正态分布的且相互独立的随机变量平方和一定符合自由度为 4 的 2 分布011( )exp()2 22(2 1)! ( )4xxf xE xxf x dx思考：如果出现不符合的情况，该如何解释？思考：如果出现不符合的情况，该如何解释？6/23/202224Toy 蒙特卡罗方法蒙特卡罗方法粒子物理与

19、核物理在实验的早期设计阶段，通常利用粒子物理与核物理在实验的早期设计阶段，通常利用Toy 蒙特卡罗来估计可达到的测量精度（也称蒙特卡罗来估计可达到的测量精度（也称黑盒子方法黑盒子方法）。）。AB+CCBD+EF+GHIJK末态有D，F，H，I，J，K。研究测量E质量时实验可以达到的分辨率。在不做探测器模拟的情况下，可以对稳定的末态粒子动量在不做探测器模拟的情况下，可以对稳定的末态粒子动量各分量进行含高斯分辨率的抽样，能损大小进行朗道分布各分量进行含高斯分辨率的抽样，能损大小进行朗道分布抽样，寿命进行指数分布抽样，等等，然后在所有末态中抽样，寿命进行指数分布抽样，等等，然后在所有末态中寻找中间不

20、稳定态寻找中间不稳定态E，根据能动量关系计算其对应的质量，根据能动量关系计算其对应的质量，得到的质量分布称为得到的质量分布称为Toy 蒙特卡罗结果蒙特卡罗结果。6/23/202225蒙特卡罗物理产生子蒙特卡罗物理产生子目的目的：将理论用于某种物理过程的事例产生将理论用于某种物理过程的事例产生输出量输出量：为对应某一物理过程的事例。对于每个事例，给出过程为对应某一物理过程的事例。对于每个事例，给出过程产生的末态粒子和对应的动量产生的末态粒子和对应的动量在粒子物理与核物理实验数据分析中在粒子物理与核物理实验数据分析中，为了验证某一理，为了验证某一理论或模型，常常需要理论家提供蒙特卡罗物理产生子。论

21、或模型，常常需要理论家提供蒙特卡罗物理产生子。6/23/202226蒙特卡罗物理产生子蒙特卡罗物理产生子(续续)简单情形简单情形eeee产生产生 与与 21)( );coscos381 ();(cos2gAAfFBFB粒子物理与核物理中常用的产生子程序包粒子物理与核物理中常用的产生子程序包强子ee强子ppWWeeJETSET(PYTHIA)HERWIGARIADNEISAJETPYTHIAHERWIGKORALWEXCALIBURERATO6/23/202227蒙特卡罗探测器模拟蒙特卡罗探测器模拟从产生子中输入粒子种类与动量，然后模拟粒子的输运过程从产生子中输入粒子种类与动量，然后模拟粒子的输

22、运过程模拟探测器响应模拟探测器响应多重散射(产生散射角)粒子衰变(产生寿命)电离能损(产生能损)电磁与强子簇射产生信号，电子学响应输出量输出量 = 模拟的数据模拟的数据输入重建分析软件输入重建分析软件用途用途: 预测预测“物理产生子层面上的物理产生子层面上的”给定假设在给定假设在“探测器层探测器层面上面上”应该观测到的响应。应该观测到的响应。通用软件包通用软件包：GEANT3(FORTRAN)，GEANT4(C+)粒子与核物理中模拟的应用n用于实验初期的设计阶段建模分析n用于了解实验可能遇到物理过程的基本特征n用于了解实验仪器自身所受到的各种影响因素与所影响的大小n用于数据分析阶段的系统分析n

23、6/23/2022286/23/202229带电粒子在水中的输运过程模拟带电粒子在水中的输运过程模拟给定带电粒子的四动量单位厘米产生多少光子？从均匀分布中产生满足一定波长分布的光子沿期伦科夫光锥方向均匀给所有光子动量每个光子开始在水中传播按光与水分子发生作用的概率抽样该光子是否被吸收或散射6/23/2022302 MeV 电子在水中的输运过程电子在水中的输运过程模拟结果显示了电子在水中发出期伦科夫光，损失能量直至被停止在水中的过程。入射电子期伦科期伦科夫光子夫光子期伦科夫期伦科夫光子被水光子被水吸收吸收2 米长2 米宽2 米高水立方空气水6/23/202231200 MeV 电子在水中的输运过

24、程电子在水中的输运过程入射电子入射电子2 米长2 米宽2 米高水立方空气空气水水图中只显示能量大于图中只显示能量大于 1MeV 的粒子的粒子原初电子在原初电子在水中的轨迹水中的轨迹电子韧致辐射电子韧致辐射产生的光子产生的光子光子在水光子在水中散射中散射发生了康普发生了康普顿效应打出顿效应打出了电子了电子探测器模拟（几何设置）6/23/202232探测器模拟（物理过程）6/23/202233这种模拟可以提供对探测器效率与预期性能的很好估计。6/23/202234CERN的蒙特卡罗模拟程序包的蒙特卡罗模拟程序包G GEANTEANT4 4 是模拟粒子经过物质时所发生的相互作用的一个软件是模拟粒子经

25、过物质时所发生的相互作用的一个软件包。包。 它的应用范围包括它的应用范围包括: : 空间科学空间科学医学物理医学物理粒子物理，核物理和加速器物理粒子物理，核物理和加速器物理http:/geant4.web.cern.ch/geant46/23/202235蒙特卡罗方法应用举例蒙特卡罗方法应用举例)理论P()实验| 理论()实验P()理论 | 实验P(P如何确定在实验条件下如何确定在实验条件下，理论的概率密度函数，理论的概率密度函数例如例如：一质量为一质量为 m 共振宽度为共振宽度为 的共振态在实验上观察的共振态在实验上观察到的到的概率密度函数概率密度函数是什么形式？是什么形式？布莱布莱特特-魏

26、格纳魏格纳分布分布探测器探测器分辨率分辨率探测探测效率效率),;(00mMBW)|(MMR)|%100(Mp贝叶斯定理贝叶斯定理：6/23/202236应用举例应用举例(续一续一)也就是说也就是说，对应于真实的，对应于真实的 M，实际的实际的 M 应该是怎样一个应该是怎样一个分布分布 )|%100()| (),;()|%100()| (),;(),; (000000dMdMMpMMRmMBWdMMpMMRmMBWmMBW如果假设如果假设22220.01GeV/ ),2)(exp(21)| (cMMMMR)1 . 01 (9 . 0)|%100(2MMp2020GeV/00426. 0 ;GeV

27、/19456. 1ccm6/23/202237应用举例应用举例(续二续二)BWRBW pRBW真实物理的图真实物理的图像在实验观测像在实验观测中会发生变化中会发生变化。如果探测器的如果探测器的影响可以用函影响可以用函数来表达，有数来表达，有时积分可积。时积分可积。但大多数数情但大多数数情况下，不能用况下，不能用函数表示时，函数表示时，蒙特卡罗方法蒙特卡罗方法可以给出最好可以给出最好的近似。的近似。6/23/202238应用举例应用举例(续三续三)应用蒙特卡罗方法的步骤应用蒙特卡罗方法的步骤：步骤一步骤一：写出布莱特：写出布莱特-魏格纳产生子魏格纳产生子输出末态粒子的四动量输出末态粒子的四动量步

28、骤二步骤二：输入末态粒子四动量输入末态粒子四动量，模拟粒子在探测器的响应，模拟粒子在探测器的响应步骤三步骤三：输入各子探测器响应输入各子探测器响应，重建探测粒子的四动量重建探测粒子的四动量步骤四步骤四：输入探测粒子四动量输入探测粒子四动量，计算，计算不变质量分布不变质量分布输出各子探测器响应输出各子探测器响应输出探测粒子四动量输出探测粒子四动量实验条件下预期的布莱实验条件下预期的布莱特特-魏格纳分布魏格纳分布6/23/202239小结小结n蒙特卡罗方法蒙特卡罗方法n随机数产生子随机数产生子n函数变换法函数变换法n舍选法舍选法n蒙特卡罗方法的精度问题蒙特卡罗方法的精度问题n在粒子与核物理中的应用在粒子与核物理中的应用利用随机数对概率或与概率有关的数值计算0，1均匀分布 r，相互独立，长周期(伪随机数)()() (1rFxxFdxxfrx的边缘分布xxfyxffryxr)()( ; ;max21n/1精度总是事例产生子与探测器模拟