00001 #include <math.h>
00002 #include <iostream.h>
00003 #include "ReactorConstants.hh"
00004
00005
00006
00007
00008
00009
00010
00011
00012 double InverseBeta(double Enu,double CosThetaLab,int order){
00013
00014
00015
00016
00017 static ReactorConstants k;
00018 double pi = k.pi;
00019 double Mproton = k.Mproton;
00020 double Delta = k.Delta;
00021 double Melectron = k.Melectron;
00022 double Gfermi = k.Gfermi;
00023 double XCunits = k.XcMeVtoCmsqrd;
00024 double CosThetaC = (0.9741+0.9756)/2;
00025
00026
00027
00028 double RadCor = 0.024;
00029
00030
00031
00032 double Emin =
00033 ((Mproton+Delta+Melectron)*(Mproton+Melectron+Delta)
00034 -Mproton*Mproton)/2/Mproton;
00035 if(Enu<Emin)return 0;
00036
00037
00038
00039 double Sigma0 = Gfermi*Gfermi*CosThetaC*CosThetaC/pi*(1+RadCor);
00040
00041
00042
00043 double f = 1.00;
00044 double f2 = 3.706;
00045 double g = 1.26;
00046
00047
00048
00049 double E0 = Enu - Delta;
00050 if(E0<Melectron)E0=Melectron;
00051 double p0 = sqrt(E0*E0-Melectron*Melectron);
00052 double v0 = p0/E0;
00053
00054 if (order==0){
00055 double XC = ((f*f+3*g*g)+(f*f-g*g)*v0*CosThetaLab)*p0*E0;
00056 XC *= Sigma0/2;
00057 XC *= XCunits;
00058 return XC;
00059 }
00060
00061
00062
00063 double Ysquared = (Delta*Delta-Melectron*Melectron)/2;
00064 double E1 = E0*(1-Enu/Mproton*(1-v0*CosThetaLab))-Ysquared/Mproton;
00065 if(E1<Melectron)E1=Melectron;
00066 double p1 = sqrt(E1*E1-Melectron*Melectron);
00067 double v1 = p1/E1;
00068
00069 double Gamma =
00070 2*(f+f2)*g*((2*E0+Delta)*(1-v0*CosThetaLab)-
00071 Melectron*Melectron/E0)
00072 +(f*f+g*g)*(Delta*(1+v0*CosThetaLab)+Melectron*Melectron/E0)
00073 +(f*f+3*g*g)*((E0+Delta)*(1-CosThetaLab/v0)-Delta)
00074 +(f*f-g*g)*((E0+Delta)*(1-CosThetaLab/v0)-Delta)*v0*CosThetaLab;
00075
00076 double XC =
00077 ((f*f+3*g*g)+(f*f-g*g)*v1*CosThetaLab)*E1*p1
00078 -Gamma/Mproton*E0*p0;
00079 XC *= Sigma0/2;
00080 XC *= XCunits;
00081 return XC;
00082
00083 }