粒子物理與核物理試驗(yàn)中的數(shù)據(jù)分析課件

9/6/20231粒子物理與核物理實(shí)驗(yàn)中的數(shù)據(jù)分析楊振偉清華大學(xué)第四講:蒙特卡羅方法8/3/20231粒子物理與核物理實(shí)驗(yàn)中的數(shù)據(jù)分析楊振偉9/6/20232上一講回顧概率的基本概念隨機(jī)變量與概率密度函數(shù)隨機(jī)變量的平均值與方差能不通過(guò)實(shí)驗(yàn)對(duì)隨機(jī)變量進(jìn)行研究嗎？8/3/20232上一講回顧概率的基本概念隨機(jī)變量與概率密度9/6/20233本講要點(diǎn)蒙特卡羅方法隨機(jī)數(shù)產(chǎn)生子任意分布抽樣之函數(shù)變換法與舍選法蒙特卡羅方法中的精度問(wèn)題在粒子物理與核物理中的應(yīng)用8/3/20233本講要點(diǎn)蒙特卡羅方法9/6/20234蒙特卡羅方法簡(jiǎn)介蒙特卡羅方法就是利用一系列隨機(jī)數(shù)來(lái)計(jì)算各種概率大小和隨機(jī)變量均值等等的數(shù)值分析技術(shù)。通常的步驟為：產(chǎn)生一系列在[0，1]之間均勻分布的隨機(jī)數(shù)。利用這些隨機(jī)數(shù)按某些概率密度函數(shù)抽樣生成我們感興趣的另一隨機(jī)序列。利用這些值來(lái)估計(jì)的一些特性，例如：通過(guò)找到在區(qū)間的比例，給出積分值。第一層面上的應(yīng)用：蒙特卡羅計(jì)算=積分第二層面上的應(yīng)用：蒙特卡羅變量=“模擬的數(shù)據(jù)”8/3/20234蒙特卡羅方法簡(jiǎn)介蒙特卡羅方法就是利用一系列9/6/20235隨機(jī)數(shù)的產(chǎn)生用物理方法產(chǎn)生真正的隨機(jī)數(shù)不可重復(fù)產(chǎn)生速度慢用數(shù)學(xué)方法產(chǎn)生偽隨機(jī)數(shù)可以重復(fù)產(chǎn)生的速度快8/3/20235隨機(jī)數(shù)的產(chǎn)生用物理方法產(chǎn)生不可重復(fù)用數(shù)學(xué)方9/6/20236真隨機(jī)數(shù)與偽隨機(jī)數(shù)美國(guó)蘭德（RAND）公司在1950年代，利用真空管中產(chǎn)生的噪音制作了一個(gè)含十萬(wàn)個(gè)真正的隨機(jī)數(shù)表，并運(yùn)用于其開(kāi)展的所有模擬研究中。真正的隨機(jī)數(shù)與偽隨機(jī)數(shù)之間的區(qū)別在于：數(shù)據(jù)串是否具有可壓縮性，即能否用更短的形式來(lái)表示。真正的隨機(jī)數(shù)是不可壓縮的，非常不規(guī)則，以至于無(wú)法用更短的形式來(lái)表示它。在粒子物理與核物理研究中，隨機(jī)數(shù)的可重復(fù)性經(jīng)常也是非常有用的，尤其是程序的調(diào)試(debugging)。8/3/20236真隨機(jī)數(shù)與偽隨機(jī)數(shù)美國(guó)蘭德（R9/6/20237隨機(jī)數(shù)產(chǎn)生子目的是使在[0，1]范圍內(nèi)產(chǎn)生的偽隨機(jī)數(shù)滿足：均勻性；相互獨(dú)立性；長(zhǎng)周期性乘同余法友情推薦M=2K

=52q+1

0周期=2K-2232513123010923651312342·1010242517124010128/3/20237隨機(jī)數(shù)產(chǎn)生子目的是使在[0，1]范圍內(nèi)9/6/20238CERN庫(kù)的隨機(jī)數(shù)產(chǎn)生子PAW用戶…gRandom->SetSeed();…Float_t

random=gRandom->Rndm(1);……Real

random(1)Call

Rmarin(ISEED,0,0)…CallRanmar(random,1)…注意:用于產(chǎn)生子的隨機(jī)數(shù)種子還可以用來(lái)保證后續(xù)進(jìn)程的隨機(jī)數(shù)不重復(fù)。Root用戶粒子物理與核物理研究中，大都采用CERN程序庫(kù)提供的隨機(jī)數(shù)產(chǎn)生子。8/3/20238CERN庫(kù)的隨機(jī)數(shù)產(chǎn)生子PAW用戶……注9/6/20239隨機(jī)數(shù)均勻性與相關(guān)性檢驗(yàn)

subroutine

mc

doubleprecision

lamda,M,x,x0,y

callhbook1(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**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=xy_old=y

enddo

return

end隨機(jī)變量第I個(gè)隨機(jī)變量第I+1個(gè)隨機(jī)變量頻數(shù)均勻性相關(guān)性8/3/20239隨機(jī)數(shù)均勻性與相關(guān)性檢驗(yàn)subrout9/6/202310隨機(jī)數(shù)均勻性與相關(guān)性檢驗(yàn)隨機(jī)變量第I個(gè)隨機(jī)變量第I+1個(gè)隨機(jī)變量頻數(shù)均勻性相關(guān)性voidrandom(){

UInt_tlambda,M,x0;

TH1F*h1=newTH1F("h1","",100,0,1);

TH2F*h2=newTH2F("h2","",100,0,1,100,0,1);lambda=1220703125;//5^13M=4294967296;//2^32x0=1;

doubley,y_old;

for(inti=0;i<10000;i++){x=(lambda*x0)%M;y=(double)x/M;h1->Fill(y);

if(i>1)h2->Fill(y_old,y);x0=x;y_old=y;}}8/3/202310隨機(jī)數(shù)均勻性與相關(guān)性檢驗(yàn)隨機(jī)變量第I9/6/202311用蒙特卡羅法計(jì)算積分對(duì)于計(jì)算積分值解析解：數(shù)值解：蒙特卡羅方法：ABOx函數(shù)必須解析可積自變量不能太多對(duì)函數(shù)是否解析可積和是否太多自變量無(wú)要求在AB區(qū)間均勻投總數(shù)為N個(gè)點(diǎn)。8/3/202311用蒙特卡羅法計(jì)算積分對(duì)于計(jì)算積分值解析解9/6/202312蒙特卡羅方法中的精度問(wèn)題采用蒙特卡羅方法(MC)計(jì)算積分與傳統(tǒng)的梯形法相比有如下特點(diǎn)一維積分:

多維積分:

對(duì)于維數(shù)大于4的積分,用蒙特卡羅方計(jì)算積分總是最好。8/3/202312蒙特卡羅方法中的精度問(wèn)題采用蒙特卡羅方法9/6/202313從均勻分布到任意分布的隨機(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)系曲線切割得到。8/3/202313從均勻分布到任意分布的隨機(jī)數(shù)函數(shù)變換法舍9/6/202314函數(shù)變換法均勻分布任意分布8/3/202314函數(shù)變換法均勻分布任意分布9/6/202315例子:指數(shù)分布抽樣抽樣效率為100%。8/3/202315例子:指數(shù)分布抽樣抽樣效率為100%?？刹捎煤瘮?shù)變換法抽樣的分布指數(shù)分布三維各向同性分布二維隨機(jī)角度的正、余弦分布高斯分布n

個(gè)自由度的

2

分布伽馬分布二項(xiàng)式分布泊松分布Student分布（/2008/reviews/monterpp.pdf）9/6/202316可采用函數(shù)變換法抽樣的分布指數(shù)分布8/3/2023169/6/202317舍選法問(wèn)題:如何找到函數(shù)的最大值?8/3/202317舍選法問(wèn)題:如何找到函數(shù)的最大值?9/6/202318舍選法舉例subroutineacc_rejrealrvec(1)call

hbook1(10,'x(r)',100,0.,10.,0.)call

hbook1(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,100

call

ranmar(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,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,*)'zgreaterthanfmax'

endif

call

hfill(10,r,0.,1.0)

callranmar(rvec,1)

u=rvec(1)*fmax

if(u.lt.z)then

callhfill(20,r,0.,1.0)

callhfill(30,r,u,1.0)

ntot=ntot+1

endifenddowrite(6,*)'ntot=',ntotreturnend8/3/202318舍選法舉例subroutineacc_9/6/202319舍選法舉例voidacc_rej(){TH1F*h11=newTH1F("h11","",100,0,10);TH1F*h12=newTH1F("h12","",100,0,10);TH2F*h2=newTH2F("h2","",100,0,10,100,0,1);doublefmax=-999.;for(inti=0;i<100;i++){doubler=gRandom->Uniform(0,10);doublef=0.5*exp(-r/2.);if(fmax<f)fmax=f;}fmax*=1.2;cout<<"fmax="<<fmax<<

endl;intntot=0;for(inti=0;i<10000;i++){doubler=gRandom->Uniform(0,10);doublez=0.5*exp(-r/2.);if(z>fmax)fmax=1.2*z;h11->Fill(r);doubleu=gRandom->Uniform(0,fmax);if(u<z){h12->Fill(r);h2->Fill(r,u);ntot+=1;}}cout<<"ntot="<<ntot<<endl;}ROOT腳本8/3/202319舍選法舉例voidacc_rej()9/6/202320舍選法舉例(續(xù))頻數(shù)頻數(shù)舍選法存在效率問(wèn)題。二維均勻分布8/3/202320舍選法舉例(續(xù))頻數(shù)頻數(shù)舍選法存在效率問(wèn)9/6/202321函數(shù)變換法與舍選法函數(shù)變換法優(yōu)點(diǎn)：100%的抽樣效率缺點(diǎn)：函數(shù)須解析可積舍選法優(yōu)點(diǎn)：方法簡(jiǎn)單，可用于非常復(fù)雜的函數(shù)缺點(diǎn)：需要估計(jì)函數(shù)最大值，而且抽樣效率低粒子物理與核物理中，對(duì)常用的概率密度函數(shù)有各種建議采用的方法(見(jiàn)/2008/reviews/monterpp.pdf)。除此之外，舍選法最為常用。初學(xué)者常犯的錯(cuò)誤：對(duì)同一個(gè)過(guò)程做計(jì)算機(jī)模擬批處理，沒(méi)有考慮在批處理結(jié)果中存在隨機(jī)數(shù)的重復(fù)性。8/3/202321函數(shù)變換法與舍選法函數(shù)變換法優(yōu)點(diǎn)：1009/6/202322常用概率密度分布函數(shù)的抽樣高斯(正態(tài))分布{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;i<kUPDATE;i++){

x=gRandom->Gaus(mean,sigma);hx->Fill(x);}}產(chǎn)生平均值為mean標(biāo)準(zhǔn)偏差為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環(huán)境下采用已有的分布，可以容易完成布置的練習(xí)。8/3/202322常用概率密度分布函數(shù)的抽樣高斯(正態(tài))分9/6/202323蒙特卡羅統(tǒng)計(jì)檢驗(yàn)例如常用來(lái)檢驗(yàn)理論與實(shí)驗(yàn)符合好壞的

2分布。四個(gè)服從N(0,1)正態(tài)分布的且相互獨(dú)立的隨機(jī)變量平方和一定符合自由度為4的

2分布思考：如果出現(xiàn)不符合的情況，該如何解釋?zhuān)?/3/202323蒙特卡羅統(tǒng)計(jì)檢驗(yàn)例如常用來(lái)檢驗(yàn)理論與實(shí)驗(yàn)9/6/202324Toy

蒙特卡羅方法粒子物理與核物理在實(shí)驗(yàn)的早期設(shè)計(jì)階段，通常利用Toy蒙特卡羅來(lái)估計(jì)可達(dá)到的測(cè)量精度（也稱黑盒子方法）。AB+CCBD+EF+GHIJK末態(tài)有D，F(xiàn)，H，I，J，K。研究測(cè)量E質(zhì)量時(shí)實(shí)驗(yàn)可以達(dá)到的分辨率。在不做探測(cè)器模擬的情況下，可以對(duì)穩(wěn)定的末態(tài)粒子動(dòng)量各分量進(jìn)行含高斯分辨率的抽樣，能損大小進(jìn)行朗道分布抽樣，壽命進(jìn)行指數(shù)分布抽樣，等等，然后在所有末態(tài)中尋找中間不穩(wěn)定態(tài)E，根據(jù)能動(dòng)量關(guān)系計(jì)算其對(duì)應(yīng)的質(zhì)量，得到的質(zhì)量分布稱為T(mén)oy

蒙特卡羅結(jié)果。8/3/202324Toy蒙特卡羅方法粒子物理與核物理在實(shí)9/6/202325蒙特卡羅物理產(chǎn)生子目的：將理論用于某種物理過(guò)程的事例產(chǎn)生輸出量：為對(duì)應(yīng)某一物理過(guò)程的事例。對(duì)于每個(gè)事例，給出過(guò)程產(chǎn)生的末態(tài)粒子和對(duì)應(yīng)的動(dòng)量在粒子物理與核物理實(shí)驗(yàn)數(shù)據(jù)分析中，為了驗(yàn)證某一理論或模型，常常需要理論家提供蒙特卡羅物理產(chǎn)生子。8/3/202325蒙特卡羅物理產(chǎn)生子目的：輸出量：在粒子物9/6/202326蒙特卡羅物理產(chǎn)生子(續(xù))簡(jiǎn)單情形產(chǎn)生

與

粒子物理與核物理中常用的產(chǎn)生子程序包JETSET(PYTHIA)HERWIGARIADNEISAJETPYTHIAHERWIGKORALWEXCALIBURERATO8/3/202326蒙特卡羅物理產(chǎn)生子(續(xù))簡(jiǎn)單情形產(chǎn)生9/6/202327蒙特卡羅探測(cè)器模擬從產(chǎn)生子中輸入粒子種類(lèi)與動(dòng)量，然后模擬粒子的輸運(yùn)過(guò)程模擬探測(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++)8/3/202327蒙特卡羅探測(cè)器模擬從產(chǎn)生子中輸入粒子種類(lèi)粒子與核物理中模擬的應(yīng)用用于實(shí)驗(yàn)初期的設(shè)計(jì)階段建模分析用于了解實(shí)驗(yàn)可能遇到物理過(guò)程的基本特征用于了解實(shí)驗(yàn)儀器自身所受到的各種影響因素與所影響的大小用于數(shù)據(jù)分析階段的系統(tǒng)分析…9/6/202328粒子與核物理中模擬的應(yīng)用用于實(shí)驗(yàn)初期的設(shè)計(jì)階段建模分析8/39/6/202329帶電粒子在水中的輸運(yùn)過(guò)程模擬給定帶電粒子的四動(dòng)量單位厘米產(chǎn)生多少光子？從均勻分布中產(chǎn)生滿足一定波長(zhǎng)分布的光子沿期倫科夫光錐方向均勻給所有光子動(dòng)量每個(gè)光子開(kāi)始在水中傳播按光與水分子發(fā)生作用的概率抽樣該光子是否被吸收或散射…8/3/202329帶電粒子在水中的輸運(yùn)過(guò)程模擬給定帶電粒子9/6/2023302MeV

電子在水中的輸運(yùn)過(guò)程模擬結(jié)果顯示了電子在水中發(fā)出期倫科夫光，損失能量直至被停止在水中的過(guò)程。入射電子期倫科夫光子期倫科夫光子被水吸收2米長(zhǎng)2米寬2米高水立方空氣水8/3/2023302MeV電子在水中的輸運(yùn)過(guò)程模擬結(jié)果9/6/202331200MeV

電子在水中的輸運(yùn)過(guò)程入射電子2米長(zhǎng)2米寬2米高水立方空氣水圖中只顯示能量大于1MeV的粒子原初電子在水中的軌跡電子韌致輻射產(chǎn)生的光子光子在水中散射發(fā)生了康普頓效應(yīng)打出了電子8/3/202331200MeV電子在水中的輸運(yùn)過(guò)程入射探測(cè)器模擬（幾何設(shè)置）9/6/202332探測(cè)器模擬（幾何設(shè)置）8/3/202332探測(cè)器模擬（物理過(guò)程）9/6/202333這種模擬可以提供對(duì)探測(cè)器效率與預(yù)期性能的很好估計(jì)。探測(cè)器模擬（物理過(guò)程）8/3/202333這種模擬可以提供對(duì)9/6/202334CERN的蒙特卡羅模擬程序包GEANT4

是模