parainput.cpp

用蒙特卡罗方法实现能谱测量的模拟,以及探测效率的计算.

		    FWHM=0.01+0.05*sqrt(Edown+0.4*Edown*Edown);
		    delta=0.427*FWHM;
		    Eget=Edown+delta*rand2(kaka);
			for(int t=0;t<1024;t++)
			{
				if(ceil(Eget)==channel_energy[t])
					num[t]=num[t]+1;
			}
		

		}//以上是能量沉积记录
	}//end for：结束一个N个光子模拟过程

m_Nget=nget;
m_Ntotal=ntotal;
	
}//end two_effects:结束两种作用模拟函数


//下面是插值函数

float CParaInput::chazhi(float x[],float y[],float e)
{
		float x1[3];
		float l;
		int t=0;
		int ii=0;
		float a,b;
		for(int i=0;i<=23;i++)
		{
			if(e==x[i])
			{	
				t=1;
				ii=i;
				break;
			}
		
	
	    	if(e>x[i]&&e<x[i+1])
			{
				for(int j=0;j<2;j++)
				x1[j]=x[i+j];
				ii=i;
			}
			   continue;

		}

		
		if(t==1)
		{
			l=y[ii];
          return l;
		}
    	
	else
	{	
		a=((e-x1[0])/(x1[1]-x1[0]))*y[ii+1];
	    b=((e-x1[1])/(x1[0]-x1[1]))*y[ii];
	}
	return (a+b);
		



}





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



void CParaInput::three_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,1000.0f,1500.0f,2000.0f};
	float sigmaedata[]={46476.0f,15765.0f,8154.0f,2899.2f,1510.0f,4063.3f,3065.6f,1433.5f,796.21f,259.49f,119.22f,38.009f,111.06f,61.075f,37.201f,16.589f,8.7618f,2.704f,1.1687f,0.3582f,0.1622f,0.090189f,0.056554f,0.02891f,0.01794f,0.0082092f,0.005043f};
	float sigmacdata[]={0.024588f,0.039799f,0.059925f,0.093192f,0.13189f,0.18037f,0.21117f,0.26029f,0.30539f,0.38311f,0.44303f,0.50444f,0.55244f,0.57217f,0.58218f,0.58042f,0.57058f,0.53287f,0.48766f,0.41889f,0.36962f,0.32616f,0.30918f,0.26963f,0.2402f,0.19566f,0.16762f};
	float sigmapdata[]={0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.003712f,0.012546f};
    long N,nget=0,ntotal=0,collidetime;
	float E,Eget;
	float r,rnew,z,theta,miu,fai,sigmae,sigmac,sigmap,sigmat,l,cdth,sdth,dtheta,tt;
	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 R,H,p;
	float E0,E1;
	R=m_radius;
	H=m_height;
	E0=m_energy1;
	E1=m_energy2;
	p=m_prob;
	N=m_num;
	signn=1;

	int channel_energy[2049];
	for(int k=0;k<2048;k++)
	{
		channel_energy[k]=k;
		num1[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(pi/4))
			continue;
		else
		{
			z=0;
		    r=2*sqrt(1-miu*miu)/miu;
	    	theta=fai;
		}
	
		while(E>1)
		{
			sigmae=chazhi3(Edata,sigmaedata,E);
		    sigmac=chazhi3(Edata,sigmacdata,E);
			sigmap=chazhi3(Edata,sigmapdata,E);//插值函数调用获得截面数据
			sigmat=sigmae+sigmac+sigmap;
			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((tt=rand2(kaka))&&tt<=sigmae/sigmat)//光电效应抽样
				E=0;
			else if((tt>sigmae/sigmat)&&(tt<=(sigmae+sigmac)/sigmat))//康普顿效应抽样
		
			{
				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 else :结束康普顿效应模拟过程

			else//电子对效应抽样
			{
			


					if((tt=rand2(kaka))&&tt<=0.333)//电子对效应获得全能峰抽样
					{
						if(sign==0)
				    	    Edown=E0;
						else
							Edown=E1;
		                FWHM=0.01+0.05*sqrt(Edown+0.4*Edown*Edown);
		                 delta=0.427*FWHM;
		                 Eget=Edown+delta*rand2(kaka);
			            for(int t=0;t<2048;t++)
						{
				          if(ceil(Eget)==channel_energy[t])
				          num1[t]=num1[t]+1;
						}
					}


					else if(tt>0.333&&tt<=0.666)//电子对效应获得单逃逸峰抽样
					{
						if(sign==0)
					     	Edown=E0-511;
						else
							Edown=E1-511;
						FWHM=0.01+0.05*sqrt(Edown+0.4*Edown*Edown);
		                 delta=0.427*FWHM;
		                 Eget=Edown+delta*rand2(kaka);
			            for(int t=0;t<2048;t++)
						{
				          if(ceil(Eget)==channel_energy[t])
				          num1[t]=num1[t]+1;
						}
						
					}
					else//电子对效应获得双逃逸峰抽样
					{
						if(sign==0)
                         	Edown=E0-1020;
						else
							Edown=E1-1020;
						FWHM=0.01+0.05*sqrt(Edown+0.4*Edown*Edown);
		                 delta=0.427*FWHM;
		                 Eget=Edown+delta*rand2(kaka);
			            for(int t=0;t<2048;t++)
						{
				          if(ceil(Eget)==channel_energy[t])
				          num1[t]=num1[t]+1;
						}
					}	
				 
			}//end else：结束电子对效应模拟过程
	
		}//end while:结束三种作用跟踪过程

	
	
		ntotal=ntotal+1;
		if(collidetime>0)
		{
			nget=nget+1;
			if(sign==0)
		        Edown=E0-E;
			else
				Edown=E1-E;
		    FWHM=0.01+0.05*sqrt(Edown+0.4*Edown*Edown);
		    delta=0.427*FWHM;
		    Eget=Edown+delta*rand2(kaka);
			for(int t=0;t<2048;t++)
			{
				if(ceil(Eget)==channel_energy[t])
					num1[t]=num1[t]+1;
			}
		

		}
	}//end for:结束一个N个光子模拟过程

m_Nget=nget;
m_Ntotal=ntotal;
	
	
}//end two_effect :结束三种效应模拟函数


float CParaInput::chazhi3(float x[],float y[],float e)
{
		float x1[3];
		float l;
		int t=0;
		int ii=0;
		float a,b;
		for(int i=0;i<=26;i++)
		{
			if(e==x[i])
			{	
				t=1;
				ii=i;
				break;
			}
		
	
	    	if(e>x[i]&&e<x[i+1])
			{
				for(int j=0;j<2;j++)
				x1[j]=x[i+j];
				ii=i;
			}
			   continue;

		}

		
		if(t==1)
		{
			l=y[ii];
            return l;
		}
    	
	else
	{	
		a=((e-x1[0])/(x1[1]-x1[0]))*y[ii+1];
	    b=((e-x1[1])/(x1[0]-x1[1]))*y[ii];
	}
	return (a+b);
		



}