InverseBeta.cpp Source File

00001 #include <math.h>
00002 #include <iostream.h>
00003 #include "ReactorConstants.hh"
00004 //
00005 // Simple implementation of Beacom/Vogel inverse beat decac 
00006 // differential cross section.  order=0 means only use (1/M)**0,
00007 // order=1 adds (1/M)**1 corrections.  Notation follows published 
00008 // article
00009 //
00010 // TAB 11/20/03
00011 //
00012 double InverseBeta(double Enu,double CosThetaLab,int order){
00013   //
00014   // Cross section constants.  Some for overall scale are just
00015   // to allow absolute comparison to published articel.
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   // Radiative correction constant
00027   //
00028   double RadCor = 0.024;
00029   //
00030   // check for threshold
00031   //
00032   double Emin = 
00033     ((Mproton+Delta+Melectron)*(Mproton+Melectron+Delta)
00034      -Mproton*Mproton)/2/Mproton;
00035   if(Enu<Emin)return 0;
00036   //
00037   // overall scale
00038   //
00039   double Sigma0 = Gfermi*Gfermi*CosThetaC*CosThetaC/pi*(1+RadCor);
00040   //
00041   // couplings
00042   //
00043   double f = 1.00;
00044   double f2 = 3.706;
00045   double g = 1.26;
00046   //
00047   // order 0 terms
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   //  order 1 terms
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 }