PDF《粒子物理与核物理实验中的》

10 call ranmar(rvec,1) r=rvec(1)*2.0-1.0 ! from (0,1) to (-1,1) z=(1.+alpha*r+beta*r**2)/(2.+2.*beta/3.) if(z.gt.fmax)then fmax=z write(6,*)'

z greater than fmax'

end if call ranmar(rvec,1) u=rvec(1)*fmax if(u.le.z)then if(abs(r).lt.angle_cut)then ntot=ntot+1 x(ntot)=r if(ntot.ge.n)go to

20 end if end if go to

10 20 continue return end 对于采用蒙特卡罗研究参数估计值问题,可用下面的舍选法例子对任意 概率密度函数进行抽样,得到模拟的数据样本.

12 MINUIT数值求解极小值(4) 做500次模拟实验,重复最大似然拟合,每次都是模拟 n=2000 个事例: program MINUIT_FIT … double precision par(npar) real h(80000) common/pawc/h integer i call hlimit(80000) call hbook2(200,'

alpha vs beta'

,100,0.,1.,100,0.,1.,0.) do i=1,500 c Initialize MINUIT, set print level to C1 … call hfill(200,real(par(1)),real(par(2)),1.0) end do call hrput(0,'

minuit_fit.hbook'

,'

N'

) … mess $shell('

f77 minuit_fit.f -o minuit_fit `cernlib`'

) mess $shell('

minuit_fit'

) fun2

300 1./(1-0.42**2)*((x-0.5)**2/0.051**2+_ (y-0.5)**2/0.111**2-2*0.42*(x-0.5)/0.051*_ (y-0.5)/0.111)-1.

100 0. 1.

100 0. 1. h/fil

1 minuit_fit.hbook zone

2 2;

slix

200 1;

sliy

200 1;

h/proj

200 h/pl 200.sliy.1;

h/pl

200 Ve/cr tmp(50);

contour

300 1

1 tmp h/pl 200.slix.1 对前面的程序添加下列指令 minuit_fit.kumac 可以从 /home/chensm/examples/paw 拷贝 minuit_fit.f 和minuit_fit.kumac 运行 PAW>

exec minuit_fit

13 MINUIT数值求解极小值(5) const int npoints=2000;

Double_t x[npoints];

Double_t angle_cut=0.95;

void minuit_fit() { // Get data points get_input_data();

// Prepare for fit const int npar=2;

TMinuit *gMinuit = new TMinuit(npar);

gMinuit->

SetFCN(fcn);

Double_t arglist[10];

Int_t ierflg = 0;

arglist[0] = 1;

gMinuit->

mnexcm( SET ERR , arglist ,1,ierflg);

// Set starting values and step sizes for parameters static Double_t vstart[npar] = {0.5, 0.5 };

static Double_t step[npar] = {0.1 , 0.1 };

gMinuit->

mnparm(0, alpha , vstart[0], step[0], 0,0,ierflg);

gMinuit->

mnparm(1, beta , vstart[1], step[1], 0,0,ierflg);

// Now ready for minimization step arglist[0] = 500;

arglist[1] = 1.;

gMinuit->

mn........