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

18/03/20091粒子物理与核物理实验中的数据分析陈少敏清华大学第四讲:蒙特卡罗方法18/03/20092上一讲回顾概率的基本概念随机变量与概率密度函数随机变量的平均值与方差能不通过实验对随机变量进行研究吗？18/03/20093本讲要点蒙特卡罗方法随机数产生子任意分布抽样之函数变换法与舍选法蒙特卡罗方法中的精度问题在粒子与核物理中的应用18/03/20094蒙特卡罗方法简介蒙特卡罗方法就是利用一系列随机数来计算各种概率大小和随机变量均值等等的数值分析技术。通常的步骤为：1)产生一系列在[0，1]之间均匀分布的随机数。2)利用这些随机数按某些概率密度函数抽样生成我们感兴趣的另一随机序列。3)利用这些值来估计的一些特性，例如：通过找到在区间的比例，给出积分值。mrrr,...,,21)(xfnxxx,...,,21xx)(xfbxabadxxf)(第一层面上的应用：蒙特卡罗计算=积分第二层面上的应用：蒙特卡罗变量=“模拟的数据”18/03/20095随机数的产生用物理方法产生真正的随机数•不可重复•产生速度慢用数学方法产生伪随机数•可以重复•产生的速度快18/03/20096真随机数与伪随机数美国兰德（RAND）公司在1950年代，利用真空管中产生的噪音制作了一个含十万个真正的随机数表，并运用于其开展的所有模拟研究中。真正的随机数与伪随机数之间的区别在于：数据串是否具有可压缩性，即能否可以用更短的形式来表示。真正的随机数是不可压缩的，是非常地不规则，以至于无法用更短的形式来表示它。在粒子物理与核物理研究中，随机数的可以重复性时常也是非常有用的，尤其是程序的debugging。18/03/20097随机数产生子目的是使在[0，1]范围内产生的伪随机数满足：均匀性；相互独立性；长周期性乘同余法110(),/iiiiMrMM余数随机数与为选定常数为随机数种子。友情推荐M=2K=52q+10周期=2K-2232513123010923651312342·10102425171240101218/03/20098CERN库的随机数产生子PAW用户…gRandom-SetSeed();…Float_trandom=gRandom-Rndm(1);……Realrandom(1)CallRmarin(ISEED,0,0)…CallRanmar(random,1)…注意:用于产生子的随机数种子还可以用来保证后续进程的随机数不重复。Root用户粒子物理与核物理研究中，大都采用CERN程序库提供的随机数产生子。18/03/20099随机数均匀性与相关性检验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大家可以在用C++语言进行上述验证。18/03/200910用蒙特卡罗法计算积分对于计算积分值()bafxdx解析解：()()()bxbxaafxdxFxFx数值解：1()()nbiaibafxdxfxn蒙特卡罗方法：abABOx()fx函数必须解析可积自变量不能太多对函数是否解析可积和是否太多自变量无要求in()()()()baNfxfxdxAOBON在AB区间均匀投总数为N个点。18/03/200911蒙特卡罗方法中的精度问题采用蒙特卡罗方法(MC)计算积分与传统的梯形法相比有如下特点一维积分:)(,1/:)(,1/:MC2为子区间的数目梯形法精度为产生的随机数精度nnnn多维积分:)(,1/:)(,1/:MC/2为子区间的数目梯形法精度为产生的随机数且与维数无关精度维数nnnn对于维数大于4的积分,用蒙特卡罗方计算积分总是最好。18/03/200912从均匀分布到任意分布的随机数函数变换法舍选法从服从一给定分布函数中找出到任一随机变量取值的累积函数，然后求出累积函数的逆。当逆函数的自变量取均匀分布值时，对应的逆函数值自动满足给定分布。均匀分布给定分布从一个随机变量与对应概率密度函数最大值构成的二维均匀分布中，按概率密度函数与自变量关系曲线切割得到。18/03/200913函数变换法nxxxxfrx,...,,)()(,1][0,21分布的随机数找出服从通过适当的变换均匀分布的随机数从在')'())'((')'(')(，))'(()'(:rrxrxFdxxfrdrrgrxxPrrP即要求)()(1rFxrxF解出令均匀分布任意分布18/03/200914例子:指数分布抽样)1log()(),('1)0(1),(:)(0/'/rrxrxrdxexexfrxxx得到并解出令指数概率密度函数2呈一一对应关系与随机变量)(rxr抽样效率为100%。可采用函数变换法抽样的分布指数分布三维各向同性分布二维随机角度的正、余弦分布高斯分布n个自由度的2分布伽马分布二项式分布泊松分布Student分布（）18/03/20091518/03/200916舍选法)(:,,)(minmax1minxxrxxxxxf即变量产生均匀分布的由第一个随机取值范围的自变量根据概率密度函数max2max,0,fruf即范围内与均匀分布在变量产生第二个独立的随机。，，,),(新进行抽样从值拒绝该否则值则接受该如果xxxfu问题:如何找到函数的最大值?x的最大值)(xf)(xf18/03/200917舍选法举例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18/03/200918舍选法举例(续))(rx随机变量)(rx随机变量)(rx随机变量频数)(xf频数边缘分布舍选法存在效率问题。二维均匀分布18/03/200919函数变换法与舍选法函数变换法优点：100%的抽样效率缺点：函数须解析可积舍选法优点：方法简单，可用于非常复杂的函数缺点：需要估计函数最大值，而且抽样效率低粒子物理与核物理中，对常用的概率密度函数有各种建议采用的方法(见)。除此之外，舍选法最为常用。初学者常犯的错误：对同一个过程做计算机模拟批处理，没有考虑在批处理结果中存在随机数的重复性。18/03/200920常用概率密度分布函数的抽样高斯(正态)分布{gROOT-Reset();hx=newTH1F(hx,“xdis.,100,-10,10);gRandom-SetSeed();Double_tx;constDouble_tsigma=2.0;constDouble_tmean=1.0;constInt_tkUPDATE=1000;for(Int_ti=0;i1000;i++){x=gRandom-Gaus(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-Exp(tau);x=gRandom-BreitWigner(me,sig);在ROOT环境下采用已有的分布，可以容易完成布置的练习。18/03/200921蒙特卡罗统计检验例如常用来检验理论与实验符合好坏的2分布。四个服从N(0,1)正态分布的且相互独立的随机变量平方和一定符合自由度为4的2分布011()exp()222(21)![]()4xxfxExxfxdx思考：如果出现不符合的情况，该如何解释？18/03/200922Toy蒙特卡罗方法粒子物理与核物理在实验的早期设计阶段，通常利用Toy蒙特卡罗来估计可达到的测量精度（也称黑盒子方法）。AB+CCBD+EF+GHIJK末态有D，F，H，I，J，K。研究测量E质量时实验可以达到的分辨率。在不做探测器模拟的情况下，可以对稳定的末态粒子动量各分量进行含高斯分辨率的抽样，能损大小进行朗道分布抽样，寿命进行指数分布抽样，等等，然后在所有末态中寻找中间不稳定态E，根据能动量关系计算其对应的质量，得到的质量分布称为Toy蒙特卡罗结果。18/03/200923蒙特卡罗物理产生子目的：将理论用于某种物理过程的事例产生输出量：为对应某一物理过程的事例。对于每个事例，给出过程产生的末态粒子和对应的动量在粒子物理与核物理实验数据分析中，为了验证某一理论或模型，常常需要理论家提供蒙特卡罗物理产生子。18/03/200924蒙特卡罗物理产生子(续)简单情形eeee产生与21)();coscos381();(cos2gAAfFBFB粒子物理与核物理中常用的产生子程序包强子ee强子ppWWeeJETSET(PYTHIA)HERWIGARIADNEISAJETPYTHIAHERWIGKORALWEXCALIBURERATO18/03/200925蒙特卡罗探测器模拟从产生子中输入粒子种类与动量，然后模拟粒子的输运过程模拟探测器响应多重散射(产生散射角)粒子衰变(产生寿命)电离能损(产生能损)电磁与强子镞射产生信号，电子学响应…输出量=模拟的数据输入重建分析软件用途:预测“物理产生子层面上的”给定假设在“探测器层面上”应该观测到的响应。通用软件包：GEANT3(FORTRAN)，GEANT4(C