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

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

=52q+1

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

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

random(1)Call

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

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隨機變量第I個隨機變量第I+1個隨機變量頻數(shù)均勻性相關性8/3/20239隨機數(shù)均勻性與相關性檢驗subrout9/6/202310隨機數(shù)均勻性與相關性檢驗隨機變量第I個隨機變量第I+1個隨機變量頻數(shù)均勻性相關性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隨機數(shù)均勻性與相關性檢驗隨機變量第I9/6/202311用蒙特卡羅法計算積分對于計算積分值解析解：數(shù)值解：蒙特卡羅方法：ABOx函數(shù)必須解析可積自變量不能太多對函數(shù)是否解析可積和是否太多自變量無要求在AB區(qū)間均勻投總數(shù)為N個點。8/3/202311用蒙特卡羅法計算積分對于計算積分值解析解9/6/202312蒙特卡羅方法中的精度問題采用蒙特卡羅方法(MC)計算積分與傳統(tǒng)的梯形法相比有如下特點一維積分:

多維積分:

對于維數(shù)大于4的積分,用蒙特卡羅方計算積分總是最好。8/3/202312蒙特卡羅方法中的精度問題采用蒙特卡羅方法9/6/202313從均勻分布到任意分布的隨機數(shù)函數(shù)變換法舍選法尋找某個函數(shù)，當函數(shù)的自變量取均勻分布值時，對應的函數(shù)值自動滿足給定分布。均勻分布給定分布從一個隨機變量與對應概率密度函數(shù)最大值構成的二維均勻分布中，按概率密度函數(shù)與自變量關系曲線切割得到。8/3/202313從均勻分布到任意分布的隨機數(shù)函數(shù)變換法舍9/6/202314函數(shù)變換法均勻分布任意分布8/3/202314函數(shù)變換法均勻分布任意分布9/6/202315例子:指數(shù)分布抽樣抽樣效率為100%。8/3/202315例子:指數(shù)分布抽樣抽樣效率為100%。可采用函數(shù)變換法抽樣的分布指數(shù)分布三維各向同性分布二維隨機角度的正、余弦分布高斯分布n

個自由度的

2

分布伽馬分布二項式分布泊松分布Student分布（/2008/reviews/monterpp.pdf）9/6/202316可采用函數(shù)變換法抽樣的分布指數(shù)分布8/3/2023169/6/202317舍選法問題:如何找到函數(shù)的最大值?8/3/202317舍選法問題:如何找到函數(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ù)舍選法存在效率問題。二維均勻分布8/3/202320舍選法舉例(續(xù))頻數(shù)頻數(shù)舍選法存在效率問9/6/202321函數(shù)變換法與舍選法函數(shù)變換法優(yōu)點：100%的抽樣效率缺點：函數(shù)須解析可積舍選法優(yōu)點：方法簡單，可用于非常復雜的函數(shù)缺點：需要估計函數(shù)最大值，而且抽樣效率低粒子物理與核物理中，對常用的概率密度函數(shù)有各種建議采用的方法(見/2008/reviews/monterpp.pdf)。除此之外，舍選法最為常用。初學者常犯的錯誤：對同一個過程做計算機模擬批處理，沒有考慮在批處理結果中存在隨機數(shù)的重復性。8/3/202321函數(shù)變換法與舍選法函數(shù)變換法優(yōu)點：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);}}產生平均值為mean標準偏差為sigma的高斯分布?？梢該Q為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)境下采用已有的分布，可以容易完成布置的練習。8/3/202322常用概率密度分布函數(shù)的抽樣高斯(正態(tài))分9/6/202323蒙特卡羅統(tǒng)計檢驗例如常用來檢驗理論與實驗符合好壞的

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

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

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

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

與

粒子物理與核物理中常用的產生子程序包JETSET(PYTHIA)HERWIGARIADNEISAJETPYTHIAHERWIGKORALWEXCALIBURERATO8/3/202326蒙特卡羅物理產生子(續(xù))簡單情形產生9/6/202327蒙特卡羅探測器模擬從產生子中輸入粒子種類與動量，然后模擬粒子的輸運過程模擬探測器響應多重散射(產生散射角)粒子衰變(產生壽命)電離能損(產生能損)電磁與強子簇射產生信號，電子學響應…輸出量=模擬的數(shù)據輸入重建分析軟件用途:預測“物理產生子層面上的”給定假設在“探測器層面上”應該觀測到的響應。通用軟件包：GEANT3(FORTRAN)，GEANT4(C++)8/3/202327蒙特卡羅探測器模擬從產生子中輸入粒子種類粒子與核物理中模擬的應用用于實驗初期的設計階段建模分析用于了解實驗可能遇到物理過程的基本特征用于了解實驗儀器自身所受到的各種影響因素與所影響的大小用于數(shù)據分析階段的系統(tǒng)分析…9/6/202328粒子與核物理中模擬的應用用于實驗初期的設計階段建模分析8/39/6/202329帶電粒子在水中的輸運過程模擬給定帶電粒子的四動量單位厘米產生多少光子？從均勻分布中產生滿足一定波長分布的光子沿期倫科夫光錐方向均勻給所有光子動量每個光子開始在水中傳播按光與水分子發(fā)生作用的概率抽樣該光子是否被吸收或散射…8/3/202329帶電粒子在水中的輸運過程模擬給定帶電粒子9/6/2023302MeV

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

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

是模