計(jì)算材料學(xué).ppt

1、2020-07-12,1,粒子物理與核物理實(shí)驗(yàn)中的數(shù)據(jù)分析,楊振偉 清華大學(xué) 第四講:蒙特卡羅方法,2020-07-12,2,上一講回顧,概率的基本概念,隨機(jī)變量與概率密度函數(shù),隨機(jī)變量的平均值與方差,能不通過實(shí)驗(yàn)對(duì)隨機(jī)變量進(jìn)行研究嗎？,2020-07-12,3,本講要點(diǎn),蒙特卡羅方法 隨機(jī)數(shù)產(chǎn)生子 任意分布抽樣之函數(shù)變換法與舍選法 蒙特卡羅方法中的精度問題 在粒子物理與核物理中的應(yīng)用,2020-07-12,4,蒙特卡羅方法簡(jiǎn)介,蒙特卡羅方法就是利用一系列隨機(jī)數(shù)來計(jì)算各種概率大小和隨機(jī)變量均值等等的數(shù)值分析技術(shù)。通常的步驟為：,產(chǎn)生一系列在0，1之間均勻分布的隨機(jī)數(shù) 。 利用這些隨機(jī)數(shù)按某些概

2、率密度函數(shù) 抽樣生成我們感興趣的另一隨機(jī)序列 。 利用這些 值來估計(jì) 的一些特性，例如：通過找到在區(qū)間 的 比例，給出積分值 。,第一層面上的應(yīng)用：,蒙特卡羅計(jì)算 = 積分,第二層面上的應(yīng)用：,蒙特卡羅變量 = “模擬的數(shù)據(jù)”,2020-07-12,5,隨機(jī)數(shù)的產(chǎn)生,用物理方法產(chǎn)生 真正的隨機(jī)數(shù),不可重復(fù) 產(chǎn)生速度慢,用數(shù)學(xué)方法產(chǎn)生 偽隨機(jī)數(shù),可以重復(fù) 產(chǎn)生的速度快,2020-07-12,6,真隨機(jī)數(shù)與偽隨機(jī)數(shù),美國(guó)蘭德（RAND）公司在1950年代，利用真空管中產(chǎn)生的噪音制作了一個(gè)含十萬個(gè)真正的隨機(jī)數(shù)表，并運(yùn)用于其開展的所有模擬研究中。,真正的隨機(jī)數(shù)與偽隨機(jī)數(shù)之間的區(qū)別在于：數(shù)據(jù)串是否具有可

3、壓縮性，即能否用更短的形式來表示。,真正的隨機(jī)數(shù)是不可壓縮的，非常不規(guī)則，以至于無法用更短的形式來表示它。,在粒子物理與核物理研究中，隨機(jī)數(shù)的可重復(fù)性經(jīng)常也是非常有用的，尤其是程序的調(diào)試(debugging)。,2020-07-12,7,隨機(jī)數(shù)產(chǎn)生子,目的是使在 0，1 范圍內(nèi)產(chǎn)生的偽隨機(jī)數(shù)滿足：,均勻性；相互獨(dú)立性；長(zhǎng)周期性,乘同余法,友情推薦,2020-07-12,8,CERN庫(kù)的隨機(jī)數(shù)產(chǎn)生子,PAW用戶, gRandom-SetSeed(); Float_t random = gRandom-Rndm(1); , Real random(1) Call Rmarin(ISEED,0,0)

4、 Call Ranmar(random,1) ,注意: 用于產(chǎn)生子的 隨機(jī)數(shù)種子還可以用 來保證后續(xù)進(jìn)程的隨 機(jī)數(shù)不重復(fù)。,Root 用戶,粒子物理與核物理研究中，大都采用CERN程序庫(kù)提供的隨機(jī)數(shù)產(chǎn)生子。,2020-07-12,9,隨機(jī)數(shù)均勻性與相關(guān)性檢驗(yàn),subroutine mc double precision lamda,M,x,x0,y call hbook1(10,r,100,0.,1.,0.) call hbook2(20,r(i+1) vs. r(i), TH1F *h1 = new TH1F(h1,100,0,1); TH2F *h2 = new TH2F(h2,100,0

5、,1,100,0,1); lambda=1220703125; /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; ,2020-07-12,11,用蒙特卡羅法計(jì)算積分,對(duì)于計(jì)算積分值,解析解：,數(shù)值解：,函數(shù)必須解析可積,自變量不能太多,對(duì)函數(shù)是否解析可積和 是否太多自變量無要求,在AB區(qū)間均勻投總數(shù)為N個(gè)點(diǎn)。,2020-07-12,12,蒙特卡羅方法中的精度問題,采用蒙特卡羅方法(MC)計(jì)算積分 與傳統(tǒng)的梯形法相比有如

6、下特點(diǎn),一維積分:,多維積分:,對(duì)于維數(shù)大于4的積分, 用蒙特卡羅方計(jì)算積分總是最好。,2020-07-12,13,從均勻分布到任意分布的隨機(jī)數(shù),函數(shù)變換法,舍選法,尋找某個(gè)函數(shù)，當(dāng)函數(shù)的自變量取均勻分布值時(shí)，對(duì)應(yīng)的函數(shù)值自動(dòng)滿足給定分布。,從一個(gè)隨機(jī)變量與對(duì)應(yīng)概率密度函數(shù)最大值構(gòu)成的二維均 勻分布中，按概率密度函數(shù)與自變量關(guān)系曲線切割得到。,2020-07-12,14,函數(shù)變換法,均勻分布,任意分布,2020-07-12,15,例子:指數(shù)分布抽樣,抽樣效率為100%。,可采用函數(shù)變換法抽樣的分布,指數(shù)分布 三維各向同性分布 二維隨機(jī)角度的正、余弦分布 高斯分布 n 個(gè)自由度的 2 分布 伽馬

7、分布 二項(xiàng)式分布 泊松分布 Student 分布 （/2008/reviews/monterpp.pdf）,2020-07-12,16,2020-07-12,17,舍選法,問題: 如何找到函數(shù)的最大值?,2020-07-12,18,舍選法舉例,subroutine acc_rej real rvec(1) call hbook1(10,x(r),100,0.,10.,0.) call 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

8、=-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=f end do fmax=1.2*fmax ntot=0 do i=1,10000 call ranmar(rvec,1) r=0+rvec(1)*(10.-0.) z=0.5*exp(-r/2.),if(z.gt.fmax)then fmax=z*1.2 write(6,*)z greater than fmax end if call hfill(10,r,0.,1.0) call ranmar(rve

9、c,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 if end do write(6,*)ntot=,ntot return end,2020-07-12,19,舍選法舉例,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

10、=-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,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 n

11、tot= ntot endl; ,ROOT腳本,2020-07-12,20,舍選法舉例(續(xù)),頻數(shù),頻數(shù),舍選法存在效率問題。,二維均勻分布,2020-07-12,21,函數(shù)變換法與舍選法,函數(shù)變換法,優(yōu)點(diǎn)：100%的抽樣效率,缺點(diǎn)：函數(shù)須解析可積,舍選法,優(yōu)點(diǎn)：方法簡(jiǎn)單，可用于非常復(fù)雜的函數(shù),缺點(diǎn)：需要估計(jì)函數(shù)最大值，而且抽樣效率低,粒子物理與核物理中，對(duì)常用的概率密度函數(shù)有各種建議采用的方法(見/2008/reviews/monterpp.pdf)。除此之外，舍選法最為常用。,初學(xué)者常犯的錯(cuò)誤：對(duì)同一個(gè)過程做計(jì)算機(jī)模擬批處理， 沒有考慮在批處理結(jié)果中存在隨

12、機(jī)數(shù)的重復(fù)性。,2020-07-12,22,常用概率密度分布函數(shù)的抽樣,高斯(正態(tài))分布, gROOT-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 kUPDATE = 1000; for ( Int_t i=0; iGaus(mean,sigma); hx-Fill(x); ,產(chǎn)生平均值為mean 標(biāo)準(zhǔn)偏差為sigma的 高斯分布。,在ROOT環(huán)境下采用已有

13、的分布，可以容易完成布置的練習(xí)。,2020-07-12,23,蒙特卡羅統(tǒng)計(jì)檢驗(yàn),例如常用來檢驗(yàn)理論與實(shí)驗(yàn)符合好壞的2 分布。,四個(gè)服從 N(0,1) 正態(tài) 分布的且相互獨(dú)立的 隨機(jī)變量平方和,一定符合自由度 為 4 的 2 分布,思考：如果出現(xiàn)不符合的情況，該如何解釋？,2020-07-12,24,Toy 蒙特卡羅方法,粒子物理與核物理在實(shí)驗(yàn)的早期設(shè)計(jì)階段，通常利用Toy 蒙特卡羅來估計(jì)可達(dá)到的測(cè)量精度（也稱黑盒子方法）。,在不做探測(cè)器模擬的情況下，可以對(duì)穩(wěn)定的末態(tài)粒子動(dòng)量 各分量進(jìn)行含高斯分辨率的抽樣，能損大小進(jìn)行朗道分布 抽樣，壽命進(jìn)行指數(shù)分布抽樣，等等，然后在所有末態(tài)中 尋找中間不穩(wěn)定態(tài)

14、E，根據(jù)能動(dòng)量關(guān)系計(jì)算其對(duì)應(yīng)的質(zhì)量， 得到的質(zhì)量分布稱為Toy 蒙特卡羅結(jié)果。,2020-07-12,25,蒙特卡羅物理產(chǎn)生子,目的： 將理論用于某種物理過程的事例產(chǎn)生,輸出量： 為對(duì)應(yīng)某一物理過程的事例。對(duì)于每個(gè)事例，給出過程 產(chǎn)生的末態(tài)粒子和對(duì)應(yīng)的動(dòng)量,在粒子物理與核物理實(shí)驗(yàn)數(shù)據(jù)分析中，為了驗(yàn)證某一理論或模型，常常需要理論家提供蒙特卡羅物理產(chǎn)生子。,2020-07-12,26,蒙特卡羅物理產(chǎn)生子(續(xù)),簡(jiǎn)單情形,產(chǎn)生 與 ,粒子物理與核物理中常用的產(chǎn)生子程序包,JETSET(PYTHIA) HERWIG ARIADNE,ISAJET PYTHIA HERWIG,KORALW EXCALIB

15、UR ERATO,2020-07-12,27,蒙特卡羅探測(cè)器模擬,從產(chǎn)生子中輸入粒子種類與動(dòng)量，然后模擬粒子的輸運(yùn)過程,模擬探測(cè)器響應(yīng),多重散射(產(chǎn)生散射角) 粒子衰變(產(chǎn)生壽命) 電離能損(產(chǎn)生能損) 電磁與強(qiáng)子簇射 產(chǎn)生信號(hào)，電子學(xué)響應(yīng) ,輸出量 = 模擬的數(shù)據(jù),輸入重建分析軟件,用途: 預(yù)測(cè)“物理產(chǎn)生子層面上的”給定假設(shè)在“探測(cè)器層面上”應(yīng)該觀測(cè)到的響應(yīng)。,通用軟件包：GEANT3(FORTRAN)，GEANT4(C+),粒子與核物理中模擬的應(yīng)用,用于實(shí)驗(yàn)初期的設(shè)計(jì)階段建模分析 用于了解實(shí)驗(yàn)可能遇到物理過程的基本特征 用于了解實(shí)驗(yàn)儀器自身所受到的各種影響因素與所影響的大小 用于數(shù)據(jù)分析階

16、段的系統(tǒng)分析 ,2020-07-12,28,2020-07-12,29,帶電粒子在水中的輸運(yùn)過程模擬,給定帶電粒子的四動(dòng)量,單位厘米產(chǎn)生多少光子？,從均勻分布中產(chǎn)生滿足一定波長(zhǎng)分布的光子,沿期倫科夫光錐方向均勻給所有光子動(dòng)量,每個(gè)光子開始在水中傳播,按光與水分子發(fā)生作用的概率抽樣該光子是否被吸收或散射,2020-07-12,30,2 MeV 電子在水中的輸運(yùn)過程,模擬結(jié)果顯示了電子在水中發(fā)出期倫科夫光，損失能量直至 被停止在水中的過程。,入射電子,期倫科夫光子,期倫科夫 光子被水 吸收,2 米長(zhǎng) 2 米寬 2 米高 水立方,空氣,水,2020-07-12,31,200 MeV 電子在水中的輸運(yùn)

17、過程,入射電子,2 米長(zhǎng) 2 米寬 2 米高 水立方,空氣,水,圖中只顯示能量大于 1MeV 的粒子,原初電子在 水中的軌跡,電子韌致輻射 產(chǎn)生的光子,光子在水 中散射,發(fā)生了康普 頓效應(yīng)打出 了電子,探測(cè)器模擬（幾何設(shè)置）,2020-07-12,32,探測(cè)器模擬（物理過程）,2020-07-12,33,這種模擬可以提供對(duì)探測(cè)器效率與預(yù)期性能的很好估計(jì)。,2020-07-12,34,CERN的蒙特卡羅模擬程序包,GEANT4 是模擬粒子經(jīng)過物質(zhì)時(shí)所發(fā)生的相互作用的一個(gè)軟件 包。 它的應(yīng)用范圍包括:,空間科學(xué),醫(yī)學(xué)物理,粒子物理，核物理和加速器物理,http:/geant4.web.cern.c

18、h/geant4,2020-07-12,35,蒙特卡羅方法應(yīng)用舉例,如何確定在實(shí)驗(yàn)條件下，理論的概率密度函數(shù),例如：一質(zhì)量為 m 共振寬度為 的共振態(tài)在實(shí)驗(yàn)上觀察 到的概率密度函數(shù)是什么形式？,布萊特-魏格納 分布,探測(cè)器 分辨率,探測(cè) 效率,貝葉斯定理：,2020-07-12,36,應(yīng)用舉例(續(xù)一),也就是說，對(duì)應(yīng)于真實(shí)的 M，實(shí)際的 M 應(yīng)該是怎樣一個(gè) 分布,如果假設(shè),2020-07-12,37,應(yīng)用舉例(續(xù)二),真實(shí)物理的圖像在實(shí)驗(yàn)觀測(cè)中會(huì)發(fā)生變化。如果探測(cè)器的影響可以用函數(shù)來表達(dá)，有時(shí)積分可積。但大多數(shù)數(shù)情況下，不能用函數(shù)表示時(shí)，蒙特卡羅方法可以給出最好的近似。,2020-07-12,38,