Hello everyone!
I’m to try integrate the expression of the photon’s cross section in the Compton scattering with ROOT via Monte Carlo. I’m using the rejection’s method but until now I didn’t have sucess. I have doubt about the units…
If somebody can help me with some observation, I thank so much!!
This is code that I’m working…
void Comptomc()
{
Int_t N=0,i=0;
const Int_t T=100000000;
Float_t pi=TMath::Pi(),cs[T],thetav[T],cst[T],theta;
Float_t ne,k,efic,sig,A,s=0,fmax=5.6e-17,a=0,b=pi,Iphi=2*pi,Itheta,It;
Float_t alpha=1./137,m=0.511e6,kk,phi,u;
printf("\n\n########################################\n");
printf( "### Compton Scattering ###\n");
printf( "########################################\n\n");
printf("What is the number of events with sucess?");
scanf("%f",&ne);
printf("\nWhat is the energy of photon?");
scanf("%f",&k);
r1=new TRandom(NULL);
r2=new TRandom(NULL+1); //O argumento "NULL+1" permite gerar séries diferentes para cada variável;
r3=new TRandom(NULL+2);
func=new TF1("csf","(((1./137)**2)/(2*0.511e6)**2)*((1./(1+((5.e3)/(0.511e6))*(1-cos(x))))**3+(1./(1+((5.e3)/(0.511e6))*(1-cos(x))))-(1./(1+((5.e3)/(0.511e6))*(1-cos(x))))**2*sin(x)**2)*sin(x)",0,pi);
A=fmax*(b-a);
do
{
i=i+1;
N=N+1;
cst[i]=fmax*(r2->Rndm());
theta=b*(r3->Rndm()); thetav[i]=theta; u=1-cos(theta);
kk=1./(1+(k/m)*u);
//cross section function;
cs[i]=((alpha**2)/(2*m**2))*(kk**3+kk-(kk**2)*sin(theta)**2)*sin(theta);
if (cst[i]<cs[i]) s=s+1;
}
while(s<ne);
efic=s/N;
Itheta=A*efic; //Integração em theta;
It=Itheta*Iphi; //Iphi é a integral em phi;
sig=(A/(TMath::Sqrt(N)))*(TMath::Sqrt(efic*(1-efic)));
printf("\nCross Section= %e +/- %e b\n\n",It,sig);
TGraph *dispersao=new TGraph(i,thetav,cst);
dispersao->SetTitle("Angular Distribution 5 keV");
c2=new TCanvas("c2","c2",700,10,400,400);
dispersao->Draw("AP");
func->Draw("same");
return ();
}