📄 parainput.cpp

📁 用蒙特卡罗方法实现能谱测量的模拟,以及探测效率的计算.
💻 CPP
📖 第 1 页 / 共 2 页
字号:
12 下一页
// ParaInput.cpp : implementation file
//

#include "stdafx.h"
#include "moni.h"
#include "ParaInput.h"
#include<vector>
#include<time.h>
#include<math.h>
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
#define pi 3.1415926
double rand2(int &i);
void junyun(long n,int m) ;
void two_effects(float E0,float E1,float p,long n,float r,float h);//两种作用过程模拟函数
void three_effects(float E0,float E1,float p,long n,float r,float h);//两种作用过程模拟函数


/////////////////////////////////////////////////////////////////////////////
// CParaInput dialog


CParaInput::CParaInput(CWnd* pParent /*=NULL*/)
	: CDialog(CParaInput::IDD, pParent)
{
	//{{AFX_DATA_INIT(CParaInput)
	m_energy1 = 700.0f;
	m_energy2 = 800.0f;
	m_num = 1000;
	m_num1 = 0;
	m_num2 = 0;
	m_prob = 0.5f;
	m_qujian = 0;
	m_Nget = 0;
	m_Ntotal = 0;
	m_lamada1 = 0.0f;
	m_lamada2 = 0.0f;
	m_relative = 0.0f;
	m_height = 4.0f;
	m_radius = 2.0f;
	//}}AFX_DATA_INIT
}


void CParaInput::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CParaInput)
	DDX_Text(pDX, IDC_eng1, m_energy1);
	DDX_Text(pDX, IDC_eng2, m_energy2);
	DDX_Text(pDX, IDC_num, m_num);
	DDX_Text(pDX, IDC_num1, m_num1);
	DDX_Text(pDX, IDC_num2, m_num2);
	DDX_Text(pDX, IDC_prob, m_prob);
	DDX_Text(pDX, IDC_qujian, m_qujian);
	DDX_Text(pDX, IDC_Nget, m_Nget);
	DDX_Text(pDX, IDC_Ntotal, m_Ntotal);
	DDX_Text(pDX, IDC_lamada1, m_lamada1);
	DDX_Text(pDX, IDC_lamada2, m_lamada2);
	DDX_Text(pDX, IDC_relative, m_relative);
	DDX_Text(pDX, IDC_height, m_height);
	DDX_Text(pDX, IDC_radius, m_radius);
	//}}AFX_DATA_MAP
}


BEGIN_MESSAGE_MAP(CParaInput, CDialog)
	//{{AFX_MSG_MAP(CParaInput)
	ON_BN_CLICKED(IDC_runjun, Onrunjun)
	ON_BN_CLICKED(IDC_rundepen, Onrundepen)
	ON_BN_CLICKED(IDC_Runtwo, OnRuntwo)
	ON_BN_CLICKED(IDC_Runthree, OnRunthree)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CParaInput message handlers

void CParaInput::Onrunjun() 
{
	// TODO: Add your control notification handler code here
    UpdateData(TRUE);
	junyun(m_num1,m_qujian);
    UpdateData(FALSE);
	

}

using namespace std;
double CParaInput::rand2(int &idum)//产生随机数的程序
{
    const unsigned long IM1=2147483563,IM2=2147483399;
    const unsigned long IA1=40014,IA2=40692,IQ1=53668,IQ2=52774;
    const unsigned long IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1;
    const unsigned long NDIV=1+IMM1/NTAB;
    const double EPS=3.0e-16,AM=1.0/IM1,RNMX=(1.0-EPS);
    static int iy=0,idum2=314159269;
    static vector<int> iv(NTAB);
    int j,k;
    double temp;
    
    if ( idum <=0 ){
        idum=(idum ==0 ? 1 : -idum);
        idum2=idum;
        for ( j=NTAB+7; j>=0; j--) {
            k=idum/IQ1;
            idum=IA1*(idum-k*IQ1)-k*IR1;
            if (idum < 0) idum+=IM1;
            if (j < NTAB) iv[j] = idum;
        }
        iy=iv[0];
    }
    k=idum/IQ1;
    idum=IA1*(idum-k*IQ1)-k*IR1;
    if (idum < 0) idum += IM1;
    k=idum2/IQ2;
    idum2=IA2*(idum2-k*IQ2)-k*IR2;
    if (idum2 < 0) idum2 +=IM2;
    j = iy/NDIV;
    iy = iv[j]-idum2;
    iv[j] = idum;
    if (iy < 1) iy += IMM1;
    if ((temp=AM*iy)>RNMX ) return RNMX;
    else return temp;
}
void CParaInput::junyun(long m_num1,int m_qujian)   // 均匀性检验程序
{
 	
   float x;
   long N;
   int kaka=12;
   int n[100];
  
   float fx[100];
   float m,lamada,lamada_a,t,sum=0,a=1;
   N=m_num1;
   m=m_qujian;

   for(int k=0;k<=m;k++)
	   fx[k]=0+k*a/m;    //creat the x belong to [0,1]
   for(k=0;k<m;k++)
	   n[k]=0;           //initiate the n[]
	   
   for(int i=1;i<=N;i++)
		
		{  
			x=rand2(kaka); 
			for(int j=0;j<m;j++)
		  {  
				if(x<fx[1+j]&&x>fx[0+j])
		        n[j]=n[j]+1;
		  }
	
		}
   

             for(i=0;i<m;i++)
		     sum=sum+(n[i]-N/m)*(n[i]-N/m);
		     lamada=m/N*sum;
			 t=(1.645+sqrt(2*(m-1)-1));
			 lamada_a=0.5*t*t;
			 m_lamada1=lamada_a;
			 m_lamada2=lamada;
				
}

void CParaInput::Onrundepen() 
{
	// TODO: Add your control notification handler code here
    UpdateData(TRUE);
	Duli(m_num2);
    UpdateData(FALSE);
}

void CParaInput::Duli(long m_num2)   //独立性检验程序
{
   float w1,w2,w3,p;
   long N;
   int kaka=12;
   double a,b=0,c=0,aver=0;
   double t[2]={0,0};
   N=m_num2;
    for(int i=1;i<=N;++i)
    {
	  
        t[i%2]=rand2(kaka);
        a=rand2(kaka);
        aver+=a;
		b+=a*a;
        c+=t[0]*t[1];
    }
    aver=aver/N;
	w1=aver*aver;
	w2=b/N;
	w3=c/N;
	p=(w3-w1)/(w2-w1);
    m_relative=p;
    }

void CParaInput::OnRuntwo() 
{
	// TODO: Add your control notification handler code here
	UpdateData(TRUE);
	two_effects(m_energy1,m_energy2,m_prob,m_num,m_radius,m_height);
	UpdateData(FALSE);
	
}

void CParaInput::two_effects(float m_energy1,float m_energy2,float m_prob,long m_num,float m_radius,float m_height)//两种作用过程模拟函数
{

	float Edata[]={1.0f,1.5f,2.0f,3.0f,4.0f,5.0f,6.0f,8.0f,10.0f,15.0f,20.0f,30.0f,40.0f,50.0f,60.0f,80.0f,100.0f,150.0f,200.0f,300.0f,400.0f,500.0f,600.0f,800.0f};
	float sigmaedata[]={28370.0f,13845.0f,6908.0f,2555.0f,1223.0f,2602.0f,1925.0f,905.3f,479.7f,164.5f,74.24f,23.86f,66.60f,36.62f,22.29f,9.978f,5.298f,1.668f,0.7378f,0.2361f,0.1099f,0.06211f,0.03939f,0.02030f};
	float sigmacdata[]={0.0220f,0.0393f,0.0568f,0.0904f,0.1209f,0.1479f,0.1722f,0.2136f,0.2480f,0.3092f,0.3486f,0.3932f,0.4153f,0.4268f,0.4319f,0.4291f,0.4215f,0.3969f,0.3691f,0.3269f,0.2944f,0.2709f,0.2512f,0.2209f};
    long N,nget=0,ntotal=0,collidetime;
	float E,Eget;
	float R,H,r,rnew,z,theta,miu,fai,sigmae,sigmac,sigmat,l,cdth,sdth,dtheta;
	float alpha,alphat,x,miul,a,b,randangle,miunew,sdf,cdf,cfn,sfn,fainew,Edown;
	float FWHM,delta;
	int kaka=12;
	int flag,sign;
	float E0,E1,p;
	E0=m_energy1;
	E1=m_energy2;
	p=m_prob;
	N=m_num;
	R=m_radius;
	H=m_height;
	int channel_energy[1025];
	for(int k=0;k<=1024;k++)
	{
		channel_energy[k]=k;
		num[k]=0;
	}//初始化多道数组

	//以上是设定初始值


	//以上是参数输入
	for(int i=0;i<N;i++)
	{
		collidetime=0;
		if(rand2(kaka)<=p)
		{
			E=E0;
			sign=0;
		}
		else
		{
			E=E1;
			sign=1;
		}
		z=-2;
		r=0;
		theta=2*pi*rand2(kaka);
		miu=2*rand2(kaka)-1;
		fai=2*pi*rand2(kaka);
		if(miu<cos(atan(R/2)))
			continue;
		else
		{
			z=0;
		    r=2*sqrt(1-miu*miu)/miu;
	    	theta=fai;
		}

		while(E>1)
		{
			sigmae=chazhi(Edata,sigmaedata,E);
		    sigmac=chazhi(Edata,sigmacdata,E);
			sigmat=sigmae+sigmac;  //插值函数调用获得截面数据
			l=-log(rand2(kaka))/sigmat;
			rnew=sqrt(r*r+l*l*(1-miu*miu)+2*r*l*sqrt(1-miu*miu)*cos(fai-theta));
			z=z+l*miu;
			cdth=(rnew*rnew+r*r-l*l*(1-miu*miu))/2/r/rnew;
			sdth=l*sqrt(1-miu*miu)*sin(fai-theta)/rnew;
			dtheta=asin(sdth);
			if(cdth<0)
				dtheta=pi-dtheta;
			theta=theta+dtheta;
			r=rnew;//以上是跟踪粒子过程坐标变换
			if((r>R)||(z>=H)||(z<0))
				break;
			else
				collidetime=collidetime+1;
			if(rand2(kaka)<sigmae/sigmat)//光电效应抽样
				E=0;
			else//康普顿效应抽样
		
			{
				alpha=E/511;
		    	flag=0;
			   while(flag==0)
			   {
				   if(rand2(kaka)<=27/(4*alpha+29))
				   {
					  x=(1+2*alpha)/(1+2*alpha*rand2(kaka));
				      if(rand2(kaka)<=0.5*(((alpha+1-x)/alpha)*((alpha+1-x)/alpha)+1))
				   	  flag=1;
				   }
				  else

				  {
					x=1+2*alpha*rand2(kaka);
				    if(rand2(kaka)<=(27/4)*((x-1)*(x-1))/(x*x*x))
					flag=1;
				  }


			   }//以上是康普顿散射后能量抽样（乘加抽样）
			  E=E/x;
			  alphat=alpha/x;
			  miul=1-1/alphat+1/alpha;
			  a=miul;
			  b=sqrt(1-a*a);
			  randangle=2*pi*rand2(kaka);
			  miunew=a*miu+b*sqrt(1-miu*miu)*cos(randangle);
			  sdf=b*sin(randangle)/sqrt(1-miu*miu);
			  cdf=(a-miu*miunew)/sqrt(1-miu*miu)/sqrt(1-miunew*miunew);
			  sfn=sdf*cos(fai)+cdf*sin(fai);
			  cfn=cdf*cos(fai)-sdf*sin(fai);
			  if(cfn<0)
			  fainew=pi-fainew;
			  fai=fainew;
			  miu=miunew;//以上是康普顿散射后坐标变换
			}
	
		}//end while：结束粒子跟踪过程
		ntotal=ntotal+1;
		if(collidetime>0)
		{
			nget=nget+1;
			if(sign==0)
				Edown=E0-E;
			else
				Edown=E1-E;
12 下一页
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -