---

MyHeSource.cpp

Go to the documentation of this file.

#include <math.h>
#include <iostream.h>
#include <fstream.h>   // file I/O
#include <iomanip.h>   // format manipulation
#include <string>
#include "ReactorConstants.hh"
#include "ReactorFortran.hh"
#include "MyHeSource.hh"

MyHeSource::MyHeSource(int newIsotope,
                       double newRmaxgen){

  Isotope = newIsotope;
  Rmaxgen = newRmaxgen;

  //cout << "MyHeSource::MyHeSource(" << Isotope << "," << Rmaxgen 
  //     << ")" << endl;

  // 8He
  if(Isotope == 8){

    // setup the decay parameters (from Jon Link)
    // Nneutron is always 1 and Ngamma always 0
    Nneutron = 1;
    Ngamma = 0;

    // branching fraction (the 1.0000 bin contains all decays 
    // which do nothing: neutron E = 0, beta Q = 0)
    float static br[4] = {0.080,0.160,0.169,1.0000};

    // Q-value for beta decay
    float static betaQ[4] = {7.442,5.252,1.852,0.000};

    // Q-value for the total decay chain
    float static totQ;
    totQ = 8.619;

    // neutron energy
    float static neuE[4] = {0.697,2.887,6.727,0.000};

    // prepare the beta decay functions
    bool static first=1;
    float static fspace[3][200];

    if(first){
      float xlo = 0.;
      funlxp_(BetaHighQ,fspace[0],xlo,betaQ[0]);
      funlxp_(BetaMiddleQ,fspace[1],xlo,betaQ[1]);
      funlxp_(BetaLowQ,fspace[2],xlo,betaQ[2]);

      first = 0;
    }
    //
    // let the 8He decay
    //
    int ntrial = 0;
    bool passed = false;
    while(!passed){
      const int len=1;
      float rv[len];
      ranlux_(rv,len);
      //
      for(int i=0;i<4; i++){
        if(rv[0] < br[i] && !passed){

          // pick an energy for each neutron
          neutronE = neuE[i];

          if(neutronE){

            // get the energy of the beta
            int len2 = 1;
            funlux_(fspace[i],betaE,len2);
            passed = true;
          }
        }
      }
      ntrial++;
    }
    Weight = 1./float(ntrial);

    //cout << "   8He neutron energy = " << neutronE << " MeV" << endl;
    //cout << "   8He beta energy = " << betaE << " MeV" << endl;
    //cout << "   8He weight = " << Weight << endl;

    // pick a decay time in ns
    double HalfLife = 1.2e8;

    const int len=1;
    float rv[len];
    ranlux_(rv,len);
    DecayTime = HalfLife*log(1/rv[0]);

  } // done with 8He
}

MyHeSource::~MyHeSource(){
  ;
}

MyHeSource::MyHeSource(const MyHeSource& HeSource){

  Isotope = HeSource.Isotope;
  Rmaxgen = HeSource.Rmaxgen;

  Nneutron = HeSource.Nneutron;
  Nbeta    = HeSource.Nbeta;
  Ngamma   = HeSource.Ngamma;
  neutronE = HeSource.neutronE;
  betaE    = HeSource.betaE;
  gammaE   = HeSource.gammaE;

  Weight = HeSource.Weight;
  DecayTime = HeSource.DecayTime;
}    

MyHeSource& MyHeSource::operator=(const MyHeSource& rhs){

  if(this == &rhs)return *this;

  Isotope = rhs.Isotope;  
  Rmaxgen = rhs.Rmaxgen;

  Nneutron = rhs.Nneutron;
  Nbeta    = rhs.Nbeta;
  Ngamma   = rhs.Ngamma;
  neutronE = rhs.neutronE;
  betaE    = rhs.betaE;
  gammaE   = rhs.gammaE;

  Weight = rhs.Weight;
  DecayTime = rhs.DecayTime;

  return *this;
}    

float MyHeSource::BetaHighQ(float const& x){

  // beta decay function from Hannah and Jon
  static double me = MELECTRON;
  float beta = 0.;
  int A = 8;
  int Z = 2;
  double Q = 7.442;

  // x -> 0, beta -> 4.222, but Fermi -> inf so don't call it
  if(x == 0.)
    return 4.222;
  else{
    beta = sqrt(x*x+2*x*me)*(Q-x)*(Q-x)*(x+me)*Fermi(Z,A,x);
    return beta;
  }
}

float MyHeSource::BetaMiddleQ(float const& x){

  // beta decay function from Hannah and Jon
  static double me = MELECTRON;
  float beta = 0.;
  int A = 8;
  int Z = 2;
  double Q = 5.252;

  // x -> 0, beta -> 1.1565, but Fermi -> inf so don't call it
  if(x == 0.)
    return 1.1565;
  else{
    beta = sqrt(x*x+2*x*me)*(Q-x)*(Q-x)*(x+me)*Fermi(Z,A,x);
    return beta;
  }
}

float MyHeSource::BetaLowQ(float const& x){

  // beta decay function from Hannah and Jon
  static double me = MELECTRON;
  float beta = 0.;
  int A = 8;
  int Z = 2;
  double Q = 1.852;

  // x -> 0, beta -> 0.1944, but Fermi -> inf so don't call it
  if(x == 0.)
    return 0.1944;
  else{
    beta = sqrt(x*x+2*x*me)*(Q-x)*(Q-x)*(x+me)*Fermi(Z,A,x);
    return beta;
  }
}

double MyHeSource::Fermi(int const& Z, int const& A, float const& x){

  // Fermi function for beta decay from Hannah and Jon
  static double me    = MELECTRON;
  static double alpha = 1./ALPHAINV;
  double fermi = 0.;

  double xx = x/me+1.;
  double pe = sqrt(xx*xx-1.);
  double g0 = 1.-(alpha*Z)*(alpha*Z);
  double rnu = alpha*Z*xx/pe;
  double r = pow(A,1./3.);

  complex<double> arg(g0,rnu);
  complex<double> cg = wgamma_(arg);
  double abscgsq = real(conj(cg)*cg);

  fermi = exp(PI*rnu)*abscgsq*pow((pe*r),(2.*(g0-1.)));
  return fermi;
}

---