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   double XCunits = XcMeVtoCmsqrd;
00018   double CosThetaC = (0.9741+0.9756)/2;
00019   //
00020   // Radiative correction constant
00021   //
00022   double RadCor = 0.024;
00023   //
00024   // check for threshold
00025   //
00026   double EminBeta = 
00027     ((MPROTON+DELTA+MELECTRON)*(MPROTON+MELECTRON+DELTA)
00028      -MPROTON*MPROTON)/2/MPROTON;
00029   if(Enu<EminBeta)return 0;
00030   //
00031   // overall scale
00032   //
00033   double Sigma0 = GFERMI*GFERMI*CosThetaC*CosThetaC/PI*(1+RadCor);
00034   //
00035   // couplings
00036   //
00037   double f = 1.00;
00038   double f2 = 3.706;
00039   double g = 1.26;
00040   //
00041   // order 0 terms
00042   //
00043   double E0 = Enu - DELTA;
00044   if(E0<MELECTRON)E0=MELECTRON;
00045   double p0 = sqrt(E0*E0-MELECTRON*MELECTRON);
00046   double v0 = p0/E0;
00047 
00048   if (order==0){
00049     double XC = ((f*f+3*g*g)+(f*f-g*g)*v0*CosThetaLab)*p0*E0;
00050     XC *= Sigma0/2;
00051     XC *= XCunits;
00052     return XC;
00053   }
00054   //
00055   //  order 1 terms
00056   //
00057   double Ysquared = (DELTA*DELTA-MELECTRON*MELECTRON)/2;
00058   double E1 = E0*(1-Enu/MPROTON*(1-v0*CosThetaLab))-Ysquared/MPROTON;
00059   if(E1<MELECTRON)E1=MELECTRON;
00060   double p1 = sqrt(E1*E1-MELECTRON*MELECTRON);
00061   double v1 = p1/E1;
00062 
00063   double Gamma = 
00064     2*(f+f2)*g*((2*E0+DELTA)*(1-v0*CosThetaLab)-
00065                 MELECTRON*MELECTRON/E0)
00066     +(f*f+g*g)*(DELTA*(1+v0*CosThetaLab)+MELECTRON*MELECTRON/E0)
00067     +(f*f+3*g*g)*((E0+DELTA)*(1-CosThetaLab/v0)-DELTA)
00068     +(f*f-g*g)*((E0+DELTA)*(1-CosThetaLab/v0)-DELTA)*v0*CosThetaLab;
00069 
00070   double XC = 
00071     ((f*f+3*g*g)+(f*f-g*g)*v1*CosThetaLab)*E1*p1
00072     -Gamma/MPROTON*E0*p0;
00073   XC *= Sigma0/2;
00074   XC *= XCunits;
00075   return XC;
00076   
00077 }