---

MyLiSource.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 <complex>
#include <string>
#include "ReactorConstants.hh"
#include "ReactorFortran.hh"
#include "MyLiSource.hh"

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

  Isotope = newIsotope;
  Rmaxgen = newRmaxgen;

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

  // 9Li
  if(Isotope == 9){

    // 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[8] = {0.0238,0.3400,0.4400,0.4475,
                          0.4550,0.4894,0.4950,1.0000};

    // Q-value for beta decay
    float static betaQ[8] = {11.176,11.176,10.826,5.666,
                             5.666,2.323,2.323,0.000};

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

    // neutron energy
    float static neuE[8] = {0.764,0.856,1.115,3.235,
                            5.473,8.816,6.578,0.000};

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

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

      first = 0;
    }
    //
    // let the 9Li decay
    //
    int ntrial = 0;
    bool passed = false;
    while(!passed){
      const int len=1;
      float rv[len];
      ranlux_(rv,len);
      //
      for(int i=0;i<8; 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 << "   9Li neutron energy = " << neutronE << " MeV" << endl;
    //cout << "   9Li beta energy = " << betaE << " MeV" << endl;
    //cout << "   9Li weight = " << Weight << endl;

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

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

  } // done with 9Li
}

MyLiSource::~MyLiSource(){
  ;
}

MyLiSource::MyLiSource(const MyLiSource& LiSource){

  Isotope = LiSource.Isotope;
  Rmaxgen = LiSource.Rmaxgen;

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

  Weight = LiSource.Weight;
  DecayTime = LiSource.DecayTime;
} 
   

MyLiSource& MyLiSource::operator=(const MyLiSource& 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 MyLiSource::BetaHighestQ(float const& x){

  // beta decay function from Hannah and Jon
  static double me = MELECTRON;
  float beta = 0.;
  int A = 9;
  int Z = 3;
  double Q = 11.176;

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

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

  // beta decay function from Hannah and Jon
  static double me = MELECTRON;
  float beta = 0.;
  int A = 9;
  int Z = 3;
  double Q = 10.826;

  // 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 MyLiSource::BetaMiddleQ(float const& x){

  // beta decay function from Hannah and Jon
  static double me = MELECTRON;
  float beta = 0.;
  int A = 9;
  int Z = 3;
  double Q = 5.666;

  // 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 MyLiSource::BetaLowQ(float const& x){

  // beta decay function from Hannah and Jon
  static double me = MELECTRON;
  float beta = 0.;
  int A = 9;
  int Z = 3;
  double Q = 2.323;

  // 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 MyLiSource::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;
}

---