粒子物理与核物理实验中的数据分析lecture-3-蒙特卡罗方法

粒子物理与核物理实验中的数据分析陈少敏清华大学第三讲:蒙特卡罗方法~chensm/lectures/lecture_3.pdf本讲要点„蒙特卡罗方法„随机数产生子„任意分布抽样之函数变换法与舍选法„蒙特卡罗方法中的精度问题„在粒子与核物理中的应用蒙特卡罗方法简介蒙特卡罗方法就是利用一系列随机数来计算各种概率和相关量的数值分析技术。通常的步骤为：1)产生一系列在[0，1]之间均匀分布的随机数。2)利用这些随机数按某些概率密度函数抽样生成我们感兴趣的另一随机序列。3)利用这些值来估计的一些特性，例如：通过找到在区间的比例，给出积分值。mrrr,...,,21)(xfnxxx,...,,21xx)(xfbxa≤≤∫badxxf)(第一层面上的应用：蒙特卡罗计算=积分第二层面上的应用：蒙特卡罗变量=“模拟的数据”随机数用物理方法产生真正的随机数不可重复产生速度慢用数学方法产生伪随机数可以重复产生的速度快随机数产生子目的是在[0，1]范围内产生的伪随机数满足：均匀性；相互独立性；长周期性乘同余法。,)/(011为随机数种子为选定常数与随机数余数ξλξλξξMMrMiiii−−=⇒≡友情推荐M=2Kλ=52q+1ξ0周期=2K-223251312030≈10923651312034≈2·101024251712040≈1012CERN库的随机数产生子PAW用户…gRandom-SetSeed();…Float_trandom=gRandom-Rndm(1);……Realrandom(1)CallRmarin(ISEED,0,0)…CallRanmar(random,1)…注意:用于产生子的随机数种子还可以用来保证后续进程的随机数不重复。Root用户粒子物理与核物理研究中，大都采用CERN程序库提供的随机数产生子。粒子物理与核物理研究中，大都采用CERN程序库提供的随机数产生子。随机数均匀性与相关性检验subroutinemcdoubleprecisionlamda,M,x,x0,ycallhbook1(10,'r',100,0.,1.,0.)callhbook2(20,'r(i+1)vs.r(i)',&100,0.,1.,100,0.,1.,0.)x0=1.lamda=1220703125!5**13M=4294967296.!2**32doi=1,10000x=Mod(lamda*x0,M)y=x/Mcallhfill(10,real(y),0.,1.0)if(i.gt.1)call&hfill(20,real(y_old),real(y),1.0)x0=xy_old=yenddoreturnend随机变量第I个随机变量第I+1个随机变量频数均匀性相关性20h/pl10;h/pl2;1zonePAWmc.fcallPAW函数变换法nxxxxfrx,...,,)()(,1][0,21分布的随机数找出服从通过适当的变换均匀分布的随机数从在∫∫∞−∞−===≤=≤')'())'((')'(')(，))'(()'(:rrxrxFdxxfrdrrgrxxPrrP即要求)()(1rFxrxF−==解出令均匀分布任意分布例子:指数分布抽样)1log()(),('1)0(1),(:)(0/'/rrxrxrdxexexfrxxx−−==≥=∫−−ξξξξξξ得到并解出令指数概率密度函数2=ξ呈一一对应关系与随机变量)(rxr抽样效率为100%。抽样效率为100%。舍选法)(:,,)(minmax1minxxrxxxxxf−+=即变量产生均匀分布的由第一个随机取值范围的自变量根据概率密度函数max2max,0,fruf=即范围内与均匀分布在变量产生第二个独立的随机。，，,),(新进行抽样从值拒绝该否则值则接受该如果xxxfu问题:如何找到函数的最大值?问题:如何找到函数的最大值?x的最大值)(xf)(xf舍选法举例subroutineacc_rejrealrvec(1)callhbook1(10,'x(r)',100,0.,10.,0.)callhbook1(20,'x(r)',100,0.,10.,0.)callhbook2(30,'f(x)vs.x(r)',100,0.,10.,100,0.,1.1,0.)fmax=-999.doi=1,100callranmar(rvec,1)r=0+rvec(1)*(10.-0.)f=0.5*exp(-r/2.)if(fmax.lt.f)fmax=fenddofmax=1.2*fmaxntot=0doi=1,10000callranmar(rvec,1)r=0+rvec(1)*(10.-0.)z=0.5*exp(-r/2.)if(z.gt.fmax)thenfmax=z*1.2write(6,*)'zgreaterthanfmax'endifcallhfill(10,r,0.,1.0)callranmar(rvec,1)u=rvec(1)*fmaxif(u.lt.z)thencallhfill(20,r,0.,1.0)callhfill(30,r,u,1.0)ntot=ntot+1endifenddowrite(6,*)'ntot=',ntotreturnend的最大值找)(xf30h/pl20;h/pl10;h/pl3;1zonePAWacc_rej.fcallPAW舍选法举例(续))(rx随机变量)(rx随机变量)(rx随机变量频数)(xf频数边缘分布舍选法存在效率问题。舍选法存在效率问题。二维均匀分布函数变换法与舍选法函数变换法优点：100%的抽样效率缺点：函数须解析可积舍选法优点：方法简单，可用于非常复杂的函数缺点：需要估计函数最大值，而且抽样效率低粒子物理与核物理中，对常用的概率密度函数有各种建议采用的方法(见)。除此之外，舍选法最为常用。粒子物理与核物理中，对常用的概率密度函数有各种建议采用的方法(见)。除此之外，舍选法最为常用。蒙特卡罗方法中的精度问题采用蒙特卡罗方法(MC)计算积分与传统的梯形法相比有如下特点一维积分:)(,1/:)(,1/:MC2为子区间的数目梯形法精度为产生的随机数精度nnnn∝∝9多维积分:)(,1/:)(,1/:MC/2为子区间的数目梯形法精度为产生的随机数且与维数无关精度维数nnnn∝∝9对于维数大于4的积分,用蒙特卡罗方计算积分总是最好。对于维数大于4的积分,用蒙特卡罗方计算积分总是最好。蒙特卡罗事例产生子目的：将理论用于某种物理过程的事例产生输出量：为对应某一物理过程的事例。对于每个事例，给出过程产生的末态粒子和对应的动量在粒子物理与核物理实验数据分析中，为了验证某一理论或模型，常常需要理论家提供蒙特卡罗产生子。在粒子物理与核物理实验数据分析中，为了验证某一理论或模型，常常需要理论家提供蒙特卡罗产生子。蒙特卡罗事例产生子(续)简单情形−+−+→μμee+e−e−μ+μ产生θ与φπφθθθ21)();coscos381();(cos2=++∝gAAfFBFB粒子物理与核物理中常用的较复杂产生子强子→−+ee强子→pp−+−+→WWeeJETSET(PYTHIA)HERWIGARIADNEISAJETPYTHIAHERWIGKORALWEXCALIBURERATO蒙特卡罗探测器模拟从产生子中输入粒子种类与动量模拟探测器响应多重库仑散射(产生散射角)粒子衰变(产生寿命)电离能损(产生能损)电磁与强子镞射产生信号，电子学响应…输出量=模拟的数据输入重建分析软件用途:预测“产生子层面上”给定的假设在“探测器层面上”应该观测到的响应。通用软件包：GEANT3(FORTRAN)，GEANT4(C++)CERN的蒙特卡罗模拟程序包GEANT4是模拟粒子经过物质时所发生的相互作用的一个软件包。它的应用范围包括:空间科学医学物理粒子物理，核物理和加速器物理蒙特卡罗方法应用举例)理论P()实验|理论()实验P()理论|实验P(=P如何确定在实验条件下，理论的概率密度函数例如：一质量为m共振宽度为Γ的共振态在实验上观察到的概率密度函数是什么形式？布莱格-魏格纳分布8探测器分辨率探测效率8),;(00ΓmMBW)|('MMR)|%100(Mp=ε贝叶斯定理：应用举例(续一)也就是说，对应于真实的M，实际的M’应该是怎样一个分布∫∫∫∞+∞−∞+∞−+∞∞−=Γ=Γ=Γ')|%100()|'(),;()|%100()|'(),;(),;'(000000dMdMMpMMRmMBWdMMpMMRmMBWmMBWεε如果假设22220.01GeV/),2)'(exp(21)|'(cMMMMR=−−=σσπσ)1.01(9.0)|%100(2MMp−==ε2020GeV/00426.0;GeV/19456.1ccm=Γ=应用举例(续二)BWRBW⊗pRBW⊗⊗真实物理的图像在实验观测中会发生变化。如果探测器的影响可以用函数来表达，有时积分可积。但大多数数情况下，不能用函数表示时，蒙特卡罗方法可以给出最好的近似。真实物理的图像在实验观测中会发生变化。如果探测器的影响可以用函数来表达，有时积分可积。但大多数数情况下，不能用函数表示时，蒙特卡罗方法可以给出最好的近似。应用举例(续三)应用蒙特卡罗方法的步骤：步骤一：写出布莱格-魏格纳产生子输出末态粒子的四动量步骤二：输入末态粒子四动量，模拟粒子在探测器的响应步骤三：输入各子探测器响应，重建探测粒子的四动量步骤四：输入探测粒子四动量，计算不变质量分布输出各子探测器响应输出探测粒子四动量实验条件下预期的布莱格-魏格纳分布小结„蒙特卡罗方法„随机数产生子„函数变换法„舍选法„蒙特卡罗方法的精度问题„在粒子与核物理中的应用利用随机数对概率或与概率有关的数值计算[0，1]均匀分布r，相互独立，长周期(伪随机数))()(')'(1rFxxFdxxfrx−∞−=→==∫的边缘分布xxfyxffryxr=→=→)()(;;max21n/1∝精度总是事例产生子与探测器模拟