G4SPSEneDistribution.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// MODULE:        G4SPSEneDistribution.cc
//
// Version:      1.0
// Date:         5/02/04
// Author:       Fan Lei 
// Organisation: QinetiQ ltd.
// Customer:     ESA/ESTEC
//
//
// CHANGE HISTORY
// --------------
//
//
// Version 1.0, 05/02/2004, Fan Lei, Created.
//    Based on the G4GeneralParticleSource class in Geant4 v6.0
//
//
//  26/03/2014, Andrew Green.
//      Modification to used STL vectors instead of C-style arrays. This should save some space,
//      particularly when the blackbody function is not used. Also moved to dynamically allocated
//      memory in the LinearInterpolation, ExpInterpolation and LogInterpolation functions. Again,
//      this will save space if these functions are unused.
#include "G4SPSEneDistribution.hh"

#include "G4Exp.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"
#include "G4AutoLock.hh"
#include "G4Threading.hh"

G4SPSEneDistribution::G4SPSEneDistribution(): Epnflag(false),eneRndm(0),Splinetemp(0)
{
    G4MUTEXINIT(mutex);
    //
    // Initialise all variables
    particle_energy = 1.0 * MeV;
    EnergyDisType = "Mono";
    weight=1.;
    MonoEnergy = 1 * MeV;
    Emin = 0.;
    Emax = 1.e30;
    alpha = 0.;
    biasalpha = 0.;
    prob_norm = 1.0;
    Ezero = 0.;
    SE = 0.;
    Temp = 0.;
    grad = 0.;
    cept = 0.;
    Biased = false; // not biased
    EnergySpec = true; // true - energy spectra, false - momentum spectra
    DiffSpec = true; // true - differential spec, false integral spec
    IntType = "NULL"; // Interpolation type
    IPDFEnergyExist = false;
    IPDFArbExist = false;

    ArbEmin = 0.;
    ArbEmax = 1.e30;

    verbosityLevel = 0;
    
    //AG: Set these pointers to NULL so we know if they were used...
    Arb_grad = NULL;
    Arb_cept = NULL;
    Arb_alpha = NULL;
    Arb_Const = NULL;
    Arb_ezero = NULL;
    Arb_grad_cept_flag = false;
    Arb_alpha_Const_flag = false;
    Arb_ezero_flag = false;
    histInit = false;
    histCalcd = false;
    
    BBHist = NULL;
    Bbody_x = NULL;

    threadLocal_t& data = threadLocalData.Get();
    data.Emax = Emax;
    data.Emin = Emin;
    data.alpha =alpha;
    data.cept = cept;
    data.Ezero = Ezero;
    data.grad = grad;
    data.particle_energy = 0;
    data.particle_definition = NULL;
    data.weight = weight;
}

G4SPSEneDistribution::~G4SPSEneDistribution()
{
    G4MUTEXDESTROY(mutex);
    if(Arb_grad_cept_flag)
    {
        delete[] Arb_grad;
        delete[] Arb_cept;
    }
    
    if(Arb_alpha_Const_flag)
    {
        delete[] Arb_alpha;
        delete[] Arb_Const;
    }
    
    if(Arb_ezero_flag)
    {
        delete[] Arb_ezero;
    }
    delete Bbody_x;
    delete BBHist;
    for ( std::vector<G4DataInterpolation*>::iterator it = SplineInt.begin() ; it!=SplineInt.end() ; ++it)
    {
        delete *it;
        *it = 0;
    }
    SplineInt.clear();
}

void G4SPSEneDistribution::SetEnergyDisType(G4String DisType) {
    G4AutoLock l(&mutex);
    EnergyDisType = DisType;
    if (EnergyDisType == "User")
    {
        UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
        IPDFEnergyExist = false;
    }
    else if (EnergyDisType == "Arb")
    {
        ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
        IPDFArbExist = false;
    }
    else if (EnergyDisType == "Epn")
    {
        UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
        IPDFEnergyExist = false;
        EpnEnergyH = ZeroPhysVector;
    }
}

G4String G4SPSEneDistribution::GetEnergyDisType()
{
    G4AutoLock l(&mutex);
    return EnergyDisType;
}

void G4SPSEneDistribution::SetEmin(G4double emi) {
    G4AutoLock l(&mutex);
    Emin = emi;
    threadLocalData.Get().Emin = Emin;
}

G4double G4SPSEneDistribution::GetEmin()
{
    return threadLocalData.Get().Emin;
}

G4double G4SPSEneDistribution::GetArbEmin()
{
    G4AutoLock l(&mutex);
    return ArbEmin;
}

G4double G4SPSEneDistribution::GetArbEmax()
{
    G4AutoLock l(&mutex);
    return ArbEmax;
}

void G4SPSEneDistribution::SetEmax(G4double ema) {
    G4AutoLock l(&mutex);
    Emax = ema;
    threadLocalData.Get().Emax = Emax;
}

G4double G4SPSEneDistribution::GetEmax()
{
    return threadLocalData.Get().Emax;
}


void G4SPSEneDistribution::SetMonoEnergy(G4double menergy) {
    G4AutoLock l(&mutex);
    MonoEnergy = menergy;
}

void G4SPSEneDistribution::SetBeamSigmaInE(G4double e) {
    G4AutoLock l(&mutex);
    SE = e;
}
void G4SPSEneDistribution::SetAlpha(G4double alp) {
    G4AutoLock l(&mutex);
    alpha = alp;
    threadLocalData.Get().alpha = alpha;
}

void G4SPSEneDistribution::SetBiasAlpha(G4double alp) {
    G4AutoLock l(&mutex);
    biasalpha = alp;
    Biased = true;
}

void G4SPSEneDistribution::SetTemp(G4double tem) {
    G4AutoLock l(&mutex);
    Temp = tem;
}

void G4SPSEneDistribution::SetEzero(G4double eze) {
    G4AutoLock l(&mutex);
    Ezero = eze;
    threadLocalData.Get().Ezero = Ezero;
}

void G4SPSEneDistribution::SetGradient(G4double gr) {
    G4AutoLock l(&mutex);
    grad = gr;
    threadLocalData.Get().grad = grad;
}

void G4SPSEneDistribution::SetInterCept(G4double c) {
    G4AutoLock l(&mutex);
    cept = c;
    threadLocalData.Get().cept = cept;
}

G4String G4SPSEneDistribution::GetIntType()
{
    G4AutoLock l(&mutex);
    return IntType;
}

void G4SPSEneDistribution::SetBiasRndm(G4SPSRandomGenerator* a)
{
    G4AutoLock l(&mutex);
    eneRndm = a;
}

void G4SPSEneDistribution::SetVerbosity(G4int a)
{
    G4AutoLock l(&mutex);
    verbosityLevel = a;
}

G4double G4SPSEneDistribution::GetWeight()
{
    return threadLocalData.Get().weight;
}

G4double G4SPSEneDistribution::GetMonoEnergy()
{
    G4AutoLock l(&mutex);
    return MonoEnergy;
}

G4double G4SPSEneDistribution::GetSE()
{
    G4AutoLock l(&mutex);
    return SE;
}

G4double G4SPSEneDistribution::Getalpha()
{
    return threadLocalData.Get().alpha;
}

G4double G4SPSEneDistribution::GetEzero()
{
    return threadLocalData.Get().Ezero;
}

G4double G4SPSEneDistribution::GetTemp()
{
    G4AutoLock l(&mutex);
    return Temp;
}

G4double G4SPSEneDistribution::Getgrad()
{
    return threadLocalData.Get().grad;
}

G4double G4SPSEneDistribution::Getcept()
{
    return threadLocalData.Get().cept;
}

G4PhysicsOrderedFreeVector G4SPSEneDistribution::GetUserDefinedEnergyHisto()
{
    G4AutoLock l(&mutex);
    return UDefEnergyH;
}

G4PhysicsOrderedFreeVector G4SPSEneDistribution::GetArbEnergyHisto()
{
    G4AutoLock l(&mutex);
    return ArbEnergyH;
}

void G4SPSEneDistribution::UserEnergyHisto(G4ThreeVector input) {
    G4AutoLock l(&mutex);
    G4double ehi, val;
    ehi = input.x();
    val = input.y();
    if (verbosityLevel > 1) {
        G4cout << "In UserEnergyHisto" << G4endl;
        G4cout << " " << ehi << " " << val << G4endl;
    }
    UDefEnergyH.InsertValues(ehi, val);
    Emax = ehi;
    threadLocalData.Get().Emax = Emax;
}

void G4SPSEneDistribution::ArbEnergyHisto(G4ThreeVector input)
{
    G4AutoLock l(&mutex);
    G4double ehi, val;
    ehi = input.x();
    val = input.y();
    if (verbosityLevel > 1)
    {
        G4cout << "In ArbEnergyHisto" << G4endl;
        G4cout << " " << ehi << " " << val << G4endl;
    }
    ArbEnergyH.InsertValues(ehi, val);
}

void G4SPSEneDistribution::ArbEnergyHistoFile(G4String filename)
{
    G4AutoLock l(&mutex);
    std::ifstream infile(filename, std::ios::in);
    if (!infile)
    {
        G4Exception("G4SPSEneDistribution::ArbEnergyHistoFile","Event0301",FatalException,"Unable to open the histo ASCII file");
    }
    G4double ehi, val;
    while (infile >> ehi >> val)
    {
        ArbEnergyH.InsertValues(ehi, val);
    }
}

void G4SPSEneDistribution::EpnEnergyHisto(G4ThreeVector input)
{
    G4AutoLock l(&mutex);
    G4double ehi, val;
    ehi = input.x();
    val = input.y();
    if (verbosityLevel > 1)
    {
        G4cout << "In EpnEnergyHisto" << G4endl;
        G4cout << " " << ehi << " " << val << G4endl;
    }
    EpnEnergyH.InsertValues(ehi, val);
    Emax = ehi;
    threadLocalData.Get().Emax = Emax;
    Epnflag = true;
}

void G4SPSEneDistribution::Calculate()
{
    G4AutoLock l(&mutex);
    if (EnergyDisType == "Cdg")
    {
        CalculateCdgSpectrum();
    }
    else if (EnergyDisType == "Bbody")
    {
        if(!histInit)
        {
            InitHists();
        }
        CalculateBbodySpectrum();
    }
}

//MT: Lock in caller
void G4SPSEneDistribution::InitHists()
{
    BBHist = new std::vector<G4double>(10001, 0.0);
    Bbody_x = new std::vector<G4double>(10001, 0.0);
    histInit = true;
}


//MT: Lock in caller
void G4SPSEneDistribution::CalculateCdgSpectrum()
{
    // This uses the spectrum from The INTEGRAL Mass Model (TIMM)
    // to generate a Cosmic Diffuse X/gamma ray spectrum.
    G4double pfact[2] = { 8.5, 112 };
    G4double spind[2] = { 1.4, 2.3 };
    G4double ene_line[3] = { 1. * keV, 18. * keV, 1E6 * keV };
    G4int n_par;

    ene_line[0] = threadLocalData.Get().Emin;
    if (threadLocalData.Get().Emin < 18 * keV)
    {
        n_par = 2;
        ene_line[2] = threadLocalData.Get().Emax;
        if (threadLocalData.Get().Emax < 18 * keV)
        {
            n_par = 1;
            ene_line[1] = threadLocalData.Get().Emax;
        }
    }
    else
    {
        n_par = 1;
        pfact[0] = 112.;
        spind[0] = 2.3;
        ene_line[1] = threadLocalData.Get().Emax;
    }

    // Create a cumulative histogram.
    CDGhist[0] = 0.;
    G4double omalpha;
    G4int i = 0;

    while (i < n_par)
    {
        omalpha = 1. - spind[i];
        CDGhist[i + 1] = CDGhist[i] + (pfact[i] / omalpha) * (std::pow(ene_line[i + 1] / keV, omalpha) - std::pow(ene_line[i] / keV,omalpha));
        i++;
    }

    // Normalise histo and divide by 1000 to make MeV.
    i = 0;
    while (i < n_par)
    {
        CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[n_par];
        //      G4cout << CDGhist[i] << CDGhist[n_par] << G4endl;
        i++;
    }
}

//MT: Lock in caller
void G4SPSEneDistribution::CalculateBbodySpectrum()
{
    // create bbody spectrum
    // Proved very hard to integrate indefinitely, so different
    // method. User inputs emin, emax and T. These are used to
    // create a 10,000 bin histogram.
    // Use photon density spectrum = 2 nu**2/c**2 * (std::exp(h nu/kT)-1)
    // = 2 E**2/h**2c**2 times the exponential
    
    G4double erange = threadLocalData.Get().Emax - threadLocalData.Get().Emin;
    G4double steps = erange / 10000.;
    
    const G4double k = 8.6181e-11; //Boltzmann const in MeV/K
    const G4double h = 4.1362e-21; // Plancks const in MeV s
    const G4double c = 3e8; // Speed of light
    const G4double h2 = h * h;
    const G4double c2 = c * c;
    G4int count = 0;
    G4double sum = 0.;
    BBHist->at(0) = 0.;
    
    while (count < 10000) {
        Bbody_x->at(count) = threadLocalData.Get().Emin + G4double(count * steps);
        G4double Bbody_y = (2. * std::pow(Bbody_x->at(count), 2.)) / (h2 * c2 * (std::exp(Bbody_x->at(count) / (k * Temp)) - 1.));
        sum = sum + Bbody_y;
        BBHist->at(count + 1) = BBHist->at(count) + Bbody_y;
        count++;
    }

    Bbody_x->at(10000) = threadLocalData.Get().Emax;
    // Normalise cumulative histo.
    count = 0;
    while (count < 10001) {
        BBHist->at(count) = BBHist->at(count) / sum;
        count++;
    }
}

void G4SPSEneDistribution::InputEnergySpectra(G4bool value) {
    G4AutoLock l(&mutex);
    // Allows user to specifiy spectrum is momentum
    EnergySpec = value; // false if momentum
    if (verbosityLevel > 1)
        G4cout << "EnergySpec has value " << EnergySpec << G4endl;
}

void G4SPSEneDistribution::InputDifferentialSpectra(G4bool value) {
    G4AutoLock l(&mutex);
    // Allows user to specify integral or differential spectra
    DiffSpec = value; // true = differential, false = integral
    if (verbosityLevel > 1)
        G4cout << "Diffspec has value " << DiffSpec << G4endl;
}

void G4SPSEneDistribution::ArbInterpolate(G4String IType) {
    G4AutoLock l(&mutex);
    if (EnergyDisType != "Arb")
        G4Exception("G4SPSEneDistribution::ArbInterpolate",
                    "Event0302",FatalException,
                    "Error: this is for arbitrary distributions");
    IntType = IType;
    ArbEmax = ArbEnergyH.GetMaxLowEdgeEnergy();
    ArbEmin = ArbEnergyH.GetMinLowEdgeEnergy();

    // Now interpolate points
    if (IntType == "Lin")
        LinearInterpolation();
    if (IntType == "Log")
        LogInterpolation();
    if (IntType == "Exp")
        ExpInterpolation();
    if (IntType == "Spline")
        SplineInterpolation();
}

//MT: Lock in caller
void G4SPSEneDistribution::LinearInterpolation() {
    // Method to do linear interpolation on the Arb points
    // Calculate equation of each line segment, max 1024.
    // Calculate Area under each segment
    // Create a cumulative array which is then normalised Arb_Cum_Area

    G4double Area_seg[1024]; // Stores area under each segment
    G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
    G4int i, count;
    G4int maxi = ArbEnergyH.GetVectorLength();
    for (i = 0; i < maxi; i++)
    {
        Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
        Arb_y[i] = ArbEnergyH(size_t(i));
    }
    // Points are now in x,y arrays. If the spectrum is integral it has to be
    // made differential and if momentum it has to be made energy.
    if (DiffSpec == false)
    {
        // Converts integral point-wise spectra to Differential
        for (count = 0; count < maxi - 1; count++)
        {
            Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
                    / (Arb_x[count + 1] - Arb_x[count]);
        }
        maxi--;
    }
    //
    if (EnergySpec == false)
    {
        // change currently stored values (emin etc) which are actually momenta
        // to energies.
        G4ParticleDefinition* pdef =threadLocalData.Get().particle_definition;
        if (pdef == NULL)
        {
            G4Exception("G4SPSEneDistribution::LinearInterpolation",
                        "Event0302",FatalException,
                        "Error: particle not defined");
        }
        else
        {
            // Apply Energy**2 = p**2c**2 + m0**2c**4
            // p should be entered as E/c i.e. without the division by c
            // being done - energy equivalent.
            G4double mass = pdef->GetPDGMass();
            // convert point to energy unit and its value to per energy unit
            G4double total_energy;
            for (count = 0; count < maxi; count++)
            {
                total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
                        * mass)); // total energy

                Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
                Arb_x[count] = total_energy - mass; // kinetic energy
            }
        }
    }
    //
    i = 1;
    // AG:
    // This *should* be the first use of Arb_grad and friends, so we *should* be able to
    // dynamically allocate here
    Arb_grad = new G4double [1024];
    Arb_cept = new G4double [1024];
    Arb_grad_cept_flag = true;
    
    Arb_grad[0] = 0.;
    Arb_cept[0] = 0.;
    Area_seg[0] = 0.;
    Arb_Cum_Area[0] = 0.;
    while (i < maxi)
    {
        // calc gradient and intercept for each segment
        Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[i - 1]);
        if (verbosityLevel == 2)
        {
            G4cout << Arb_grad[i] << G4endl;
        }
        if (Arb_grad[i] > 0.)
        {
            if (verbosityLevel == 2)
            {
                G4cout << "Arb_grad is positive" << G4endl;
            }
            Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]);
        }
        else if (Arb_grad[i] < 0.)
        {
            if (verbosityLevel == 2)
            {
                G4cout << "Arb_grad is negative" << G4endl;
            }
            Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]);
        }
        else
        {
            if (verbosityLevel == 2)
            {
                G4cout << "Arb_grad is 0." << G4endl;
            }
            Arb_cept[i] = Arb_y[i];
        }

        Area_seg[i] = ((Arb_grad[i] / 2) * (Arb_x[i] * Arb_x[i] - Arb_x[i - 1] * Arb_x[i - 1]) + Arb_cept[i] * (Arb_x[i] - Arb_x[i - 1]));
        Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
        sum = sum + Area_seg[i];
        if (verbosityLevel == 2)
        {
            G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i] << G4endl;
        }
        i++;
    }

    i = 0;
    while (i < maxi)
    {
        Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation
        IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
        i++;
    }

    // now scale the ArbEnergyH, needed by Probability()
    ArbEnergyH.ScaleVector(1., 1./sum);

    if (verbosityLevel >= 1)
    {
        G4cout << "Leaving LinearInterpolation" << G4endl;
        ArbEnergyH.DumpValues();
        IPDFArbEnergyH.DumpValues();
    }
}

//MT: Lock in caller
void G4SPSEneDistribution::LogInterpolation()
{
    // Interpolation based on Logarithmic equations
    // Generate equations of line segments
    // y = Ax**alpha => log y = alpha*logx + logA
    // Find area under line segments
    // create normalised, cumulative array Arb_Cum_Area
    G4double Area_seg[1024]; // Stores area under each segment
    G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
    G4int i, count;
    G4int maxi = ArbEnergyH.GetVectorLength();
    for (i = 0; i < maxi; i++)
    {
        Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
        Arb_y[i] = ArbEnergyH(size_t(i));
    }
    // Points are now in x,y arrays. If the spectrum is integral it has to be
    // made differential and if momentum it has to be made energy.
    if (DiffSpec == false)
    {
        // Converts integral point-wise spectra to Differential
        for (count = 0; count < maxi - 1; count++)
        {
            Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) / (Arb_x[count + 1] - Arb_x[count]);
        }
        maxi--;
    }
    //
    if (EnergySpec == false)
    {
        // change currently stored values (emin etc) which are actually momenta
        // to energies.
        G4ParticleDefinition* pdef = threadLocalData.Get().particle_definition;
        if (pdef == NULL)
        {
            G4Exception("G4SPSEneDistribution::LogInterpolation",
                        "Event0302",FatalException,
                        "Error: particle not defined");
        }
        else
        {
            // Apply Energy**2 = p**2c**2 + m0**2c**4
            // p should be entered as E/c i.e. without the division by c
            // being done - energy equivalent.
            G4double mass = pdef->GetPDGMass();
            // convert point to energy unit and its value to per energy unit
            G4double total_energy;
            for (count = 0; count < maxi; count++)
            {
                total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
                        * mass)); // total energy

                Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
                Arb_x[count] = total_energy - mass; // kinetic energy
            }
        }
    }
    //
    i = 1;
    
    //AG:  *Should* be the first use of these two arrays
    if ( Arb_ezero ) { delete[] Arb_ezero; Arb_ezero = 0; }
    if ( Arb_Const ) { delete[] Arb_Const; Arb_Const = 0; }
    Arb_alpha = new G4double [1024];
    Arb_Const = new G4double [1024];
    Arb_alpha_Const_flag = true;
    
    
    Arb_alpha[0] = 0.;
    Arb_Const[0] = 0.;
    Area_seg[0] = 0.;
    Arb_Cum_Area[0] = 0.;
    if (Arb_x[0] <= 0. || Arb_y[0] <= 0.)
    {
        G4cout << "You should not use log interpolation with points <= 0." << G4endl;
        G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl;
        if (Arb_x[0] <= 0.)
        {
            Arb_x[0] = 1e-20;
        }
        if (Arb_y[0] <= 0.)
        {
            Arb_y[0] = 1e-20;
        }
    }

    G4double alp;
    while (i < maxi)
    {
        // Incase points are negative or zero
        if (Arb_x[i] <= 0. || Arb_y[i] <= 0.)
        {
            G4cout << "You should not use log interpolation with points <= 0." << G4endl;
            G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl;
            if (Arb_x[i] <= 0.)
            {
                Arb_x[i] = 1e-20;
            }
            if (Arb_y[i] <= 0.)
            {
                Arb_y[i] = 1e-20;
            }
        }

        Arb_alpha[i] = (std::log10(Arb_y[i]) - std::log10(Arb_y[i - 1])) / (std::log10(Arb_x[i]) - std::log10(Arb_x[i - 1]));
        Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[i], Arb_alpha[i]));
        alp = Arb_alpha[i] + 1;
        if (alp == 0.)
        {
          Area_seg[i] = Arb_Const[i] * (std::log(Arb_x[i]) - std::log(Arb_x[i - 1])); 
        }
        else
        {
          Area_seg[i] = (Arb_Const[i] / alp) * (std::pow(Arb_x[i], alp) - std::pow(Arb_x[i - 1], alp));
        }
        sum = sum + Area_seg[i];
        Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
        if (verbosityLevel == 2)
        {
            G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl;
        }
        i++;
    }

    i = 0;
    while (i < maxi)
    {
        Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
        IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
        i++;
    }

    // now scale the ArbEnergyH, needed by Probability()
    ArbEnergyH.ScaleVector(1., 1./sum);

    if (verbosityLevel >= 1)
    {
        G4cout << "Leaving LogInterpolation " << G4endl;
    }
}

//MT: Lock in caller
void G4SPSEneDistribution::ExpInterpolation()
{
    // Interpolation based on Exponential equations
    // Generate equations of line segments
    // y = Ae**-(x/e0) => ln y = -x/e0 + lnA
    // Find area under line segments
    // create normalised, cumulative array Arb_Cum_Area
    G4double Area_seg[1024]; // Stores area under each segment
    G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
    G4int i, count;
    G4int maxi = ArbEnergyH.GetVectorLength();
    for (i = 0; i < maxi; i++)
    {
        Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
        Arb_y[i] = ArbEnergyH(size_t(i));
    }
    // Points are now in x,y arrays. If the spectrum is integral it has to be
    // made differential and if momentum it has to be made energy.
    if (DiffSpec == false)
    {
        // Converts integral point-wise spectra to Differential
        for (count = 0; count < maxi - 1; count++)
        {
            Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
                    / (Arb_x[count + 1] - Arb_x[count]);
        }
        maxi--;
    }
    //
    if (EnergySpec == false)
    {
        // change currently stored values (emin etc) which are actually momenta
        // to energies.
        G4ParticleDefinition* pdef = threadLocalData.Get().particle_definition;
        if (pdef == NULL)
        {
            G4Exception("G4SPSEneDistribution::ExpInterpolation",
                        "Event0302",FatalException,
                        "Error: particle not defined");
        }
        else
        {
            // Apply Energy**2 = p**2c**2 + m0**2c**4
            // p should be entered as E/c i.e. without the division by c
            // being done - energy equivalent.
            G4double mass = pdef->GetPDGMass();
            // convert point to energy unit and its value to per energy unit
            G4double total_energy;
            for (count = 0; count < maxi; count++)
            {
                total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass * mass)); // total energy
                Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
                Arb_x[count] = total_energy - mass; // kinetic energy
            }
        }
    }
    //
    i = 1;
    
    //AG: Should be first use...
    if ( Arb_ezero ) { delete[] Arb_ezero; Arb_ezero = 0; }
    if ( Arb_Const ) { delete[] Arb_Const; Arb_Const = 0; }
    Arb_ezero = new G4double [1024];
    Arb_Const = new G4double [1024];
    Arb_ezero_flag = true;
    
    Arb_ezero[0] = 0.;
    Arb_Const[0] = 0.;
    Area_seg[0] = 0.;
    Arb_Cum_Area[0] = 0.;
    while (i < maxi)
    {
        G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i - 1]);
        if (test > 0. || test < 0.)
        {
            Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1]) / (std::log(Arb_y[i]) - std::log(Arb_y[i - 1]));
            Arb_Const[i] = Arb_y[i] / (std::exp(-Arb_x[i] / Arb_ezero[i]));
            Area_seg[i] = -(Arb_Const[i] * Arb_ezero[i]) * (std::exp(-Arb_x[i] / Arb_ezero[i]) - std::exp(-Arb_x[i - 1] / Arb_ezero[i]));
        }
        else
        {
            G4Exception("G4SPSEneDistribution::ExpInterpolation",
                        "Event0302",JustWarning,
                        "Flat line segment: problem, setting to zero parameters.");
            G4cout << "Flat line segment: problem" << G4endl;
            Arb_ezero[i] = 0.;
            Arb_Const[i] = 0.;
            Area_seg[i] = 0.;
        }
        sum = sum + Area_seg[i];
        Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
        if (verbosityLevel == 2)
        {
            G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl;
        }
        i++;
    }

    i = 0;
    while (i < maxi)
    {
        Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
        IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
        i++;
    }

    // now scale the ArbEnergyH, needed by Probability()
    ArbEnergyH.ScaleVector(1., 1./sum);

    if (verbosityLevel >= 1)
    {
        G4cout << "Leaving ExpInterpolation " << G4endl;
    }
}

//MT: Lock in caller
void G4SPSEneDistribution::SplineInterpolation()
{
    // Interpolation using Splines.
    // Create Normalised arrays, make x 0->1 and y hold
    // the function (Energy)
        // 
        // Current method based on the above will not work in all cases. 
        // new method is implemented below.
  
    G4double sum, Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
    G4int i, count;

    G4int maxi = ArbEnergyH.GetVectorLength();
    for (i = 0; i < maxi; i++) {
        Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
        Arb_y[i] = ArbEnergyH(size_t(i));
    }
    // Points are now in x,y arrays. If the spectrum is integral it has to be
    // made differential and if momentum it has to be made energy.
    if (DiffSpec == false) {
        // Converts integral point-wise spectra to Differential
        for (count = 0; count < maxi - 1; count++)
        {
            Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) / (Arb_x[count + 1] - Arb_x[count]);
        }
        maxi--;
    }
    //
    if (EnergySpec == false) {
        // change currently stored values (emin etc) which are actually momenta
        // to energies.
        G4ParticleDefinition* pdef = threadLocalData.Get().particle_definition;
        if (pdef == NULL)
                    G4Exception("G4SPSEneDistribution::SplineInterpolation",
                                "Event0302",FatalException,
                    "Error: particle not defined");
        else {
            // Apply Energy**2 = p**2c**2 + m0**2c**4
            // p should be entered as E/c i.e. without the division by c
            // being done - energy equivalent.
            G4double mass = pdef->GetPDGMass();
            // convert point to energy unit and its value to per energy unit
            G4double total_energy;
            for (count = 0; count < maxi; count++) {
                total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
                        * mass)); // total energy

                Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
                Arb_x[count] = total_energy - mass; // kinetic energy
            }
        }
    }

    //
    i = 1;
    Arb_Cum_Area[0] = 0.;
    sum = 0.;
    Splinetemp = new G4DataInterpolation(Arb_x, Arb_y, maxi, 0., 0.);
    G4double ei[101],prob[101];
    for ( std::vector<G4DataInterpolation*>::iterator it = SplineInt.begin() ; it!=SplineInt.end() ; ++it)
    {
        delete *it;
        *it = 0;
    }
    SplineInt.clear();
    SplineInt.resize(1024,NULL);
    while (i < maxi) {
      // 100 step per segment for the integration of area
      G4double de = (Arb_x[i] - Arb_x[i - 1])/100.;
      G4double area = 0.;

      for (count = 0; count < 101; count++) {
        ei[count] = Arb_x[i - 1] + de*count ;
        prob[count] =  Splinetemp->CubicSplineInterpolation(ei[count]);
        if (prob[count] < 0.) { 
              G4ExceptionDescription ED;
          ED << "Warning: G4DataInterpolation returns value < 0  " << prob[count] <<" "<<ei[count]<< G4endl;
              G4Exception("G4SPSEneDistribution::SplineInterpolation","Event0303",
              FatalException,ED);
        }
        area += prob[count]*de;
      }
      Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + area;
      sum += area; 

      prob[0] = prob[0]/(area/de);
      for (count = 1; count < 100; count++)
        prob[count] = prob[count-1] + prob[count]/(area/de);

      SplineInt[i]=new G4DataInterpolation(prob, ei, 101, 0., 0.);
      // note i start from 1!
      i++;
    }
    i = 0;
    while (i < maxi) {
        Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation
        IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
        i++;
    }
    // now scale the ArbEnergyH, needed by Probability()
    ArbEnergyH.ScaleVector(1., 1./sum);

    if (verbosityLevel > 0)
      G4cout << "Leaving SplineInterpolation " << G4endl;
}

void G4SPSEneDistribution::GenerateMonoEnergetic() {
    // Method to generate MonoEnergetic particles.
    threadLocalData.Get().particle_energy = MonoEnergy;
}

void G4SPSEneDistribution::GenerateGaussEnergies() {
    // Method to generate Gaussian particles.
    G4double ene = G4RandGauss::shoot(MonoEnergy,SE);
    if (ene < 0) ene = 0.;
    threadLocalData.Get().particle_energy = ene;
}

void G4SPSEneDistribution::GenerateLinearEnergies(G4bool bArb = false) {
    G4double rndm;
    threadLocal_t& params = threadLocalData.Get();
    G4double emaxsq = std::pow(params.Emax, 2.); //Emax squared
    G4double eminsq = std::pow(params.Emin, 2.); //Emin squared
    G4double intersq = std::pow(params.cept, 2.); //cept squared

    if (bArb)
        rndm = G4UniformRand();
    else
        rndm = eneRndm->GenRandEnergy();

    G4double bracket = ((params.grad / 2.) * (emaxsq - eminsq) + params.cept * (params.Emax - params.Emin));
    bracket = bracket * rndm;
    bracket = bracket + (params.grad / 2.) * eminsq + params.cept * params.Emin;
    // Now have a quad of form m/2 E**2 + cE - bracket = 0
    bracket = -bracket;
    //  G4cout << "BRACKET" << bracket << G4endl;
    if (params.grad != 0.)
    {
        G4double sqbrack = (intersq - 4 * (params.grad / 2.) * (bracket));
        //      G4cout << "SQBRACK" << sqbrack << G4endl;
        sqbrack = std::sqrt(sqbrack);
        G4double root1 = -params.cept + sqbrack;
        root1 = root1 / (2. * (params.grad / 2.));

        G4double root2 = -params.cept - sqbrack;
        root2 = root2 / (2. * (params.grad / 2.));

        //      G4cout << root1 << " roots " << root2 << G4endl;

        if (root1 > params.Emin && root1 < params.Emax)
        {
            params.particle_energy = root1;
        }
        if (root2 > params.Emin && root2 < params.Emax)
        {
            params.particle_energy = root2;
        }
    }
    else if (params.grad == 0.)
    {
        // have equation of form cE - bracket =0
        params.particle_energy = bracket / params.cept;
    }

    if (params.particle_energy < 0.)
    {
        params.particle_energy = -params.particle_energy;
    }

    if (verbosityLevel >= 1)
    {
        G4cout << "Energy is " << params.particle_energy << G4endl;
    }
}

void G4SPSEneDistribution::GeneratePowEnergies(G4bool bArb = false)
{
    // Method to generate particle energies distributed as
    // a power-law

    G4double rndm;
    G4double emina, emaxa;
    
    threadLocal_t& params = threadLocalData.Get();
    
    emina = std::pow(params.Emin, params.alpha + 1);
    emaxa = std::pow(params.Emax, params.alpha + 1);

    if (bArb)
    {
        rndm = G4UniformRand();
    }
    else
    {
        rndm = eneRndm->GenRandEnergy();
    }

    if (params.alpha != -1.)
    {
        G4double ene = ((rndm * (emaxa - emina)) + emina);
        ene = std::pow(ene, (1. / (params.alpha + 1.)));
        params.particle_energy = ene;
    }
    else
    {
        G4double ene = (std::log(params.Emin) + rndm * (std::log(params.Emax) - std::log(params.Emin)));
        params.particle_energy = std::exp(ene);
    }
    if (verbosityLevel >= 1)
    {
        G4cout << "Energy is " << params.particle_energy << G4endl;
    }
}

void G4SPSEneDistribution::GenerateBiasPowEnergies()//G4double& ene,G4double& wweight)
{
  // Method to generate particle energies distributed as
  // in biased power-law and calculate its weight
    
    threadLocal_t& params = threadLocalData.Get();
    
    G4double rndm;
    G4double emina, emaxa, emin, emax;

    G4double normal = 1.;

    emin = params.Emin;
    emax = params.Emax;
    //  if (EnergyDisType == "Arb") {
    //  emin = ArbEmin;
    //  emax = ArbEmax;
    //}

    rndm = eneRndm->GenRandEnergy();

    if (biasalpha != -1.)
    {
        emina = std::pow(emin, biasalpha + 1);
        emaxa = std::pow(emax, biasalpha + 1);
        G4double ee = ((rndm * (emaxa - emina)) + emina);
        params.particle_energy = std::pow(ee, (1. / (biasalpha + 1.)));
        normal = 1./(1+biasalpha) * (emaxa - emina);
    }
    else
    {
        G4double ee = (std::log(emin) + rndm * (std::log(emax) - std::log(emin)));
        params.particle_energy = std::exp(ee);
        normal = std::log(emax) - std::log(emin) ;
    }
    params.weight = GetProbability(params.particle_energy)
                      / (std::pow(params.particle_energy,biasalpha)/normal);

    if (verbosityLevel >= 1)
    {
        G4cout << "Energy is " << params.particle_energy << G4endl;
    }
}

void G4SPSEneDistribution::GenerateExpEnergies(G4bool bArb = false)
{
    // Method to generate particle energies distributed according
    // to an exponential curve.
    G4double rndm;

    if (bArb)
    {
        rndm = G4UniformRand();
    }
    else
    {
        rndm = eneRndm->GenRandEnergy();
    }

    threadLocal_t& params = threadLocalData.Get();
    params.particle_energy = -params.Ezero * (std::log(rndm * (std::exp(-params.Emax / params.Ezero) - std::exp(-params.Emin / params.Ezero)) + std::exp(-params.Emin / params.Ezero)));
    if (verbosityLevel >= 1)
    {
        G4cout << "Energy is " << params.particle_energy  << G4endl;
    }
}

void G4SPSEneDistribution::GenerateBremEnergies()
{
    // Method to generate particle energies distributed according
    // to a Bremstrahlung equation of
    // form I = const*((kT)**1/2)*E*(e**(-E/kT))

    G4double rndm;
    rndm = eneRndm->GenRandEnergy();
    G4double expmax, expmin, k;
    
    k = 8.6181e-11; // Boltzmann's const in MeV/K
    G4double ksq = std::pow(k, 2.); // k squared
    G4double Tsq = std::pow(Temp, 2.); // Temp squared

    threadLocal_t& params = threadLocalData.Get();

    expmax = std::exp(-params.Emax / (k * Temp));
    expmin = std::exp(-params.Emin / (k * Temp));

    // If either expmax or expmin are zero then this will cause problems
    // Most probably this will be because T is too low or E is too high

    if (expmax == 0.)
    {
        G4Exception("G4SPSEneDistribution::GenerateBremEnergies",
                    "Event0302",FatalException,
                    "*****EXPMAX=0. Choose different E's or Temp");
    }
    if (expmin == 0.)
    {
        G4Exception("G4SPSEneDistribution::GenerateBremEnergies",
                    "Event0302",FatalException,
                    "*****EXPMIN=0. Choose different E's or Temp");
    }

    G4double tempvar = rndm * ((-k) * Temp * (params.Emax * expmax - params.Emin * expmin) - (ksq * Tsq * (expmax - expmin)));

    G4double bigc = (tempvar - k * Temp * params.Emin * expmin - ksq * Tsq * expmin) / (-k * Temp);

    // This gives an equation of form: Ee(-E/kT) + kTe(-E/kT) - C =0
    // Solve this iteratively, step from Emin to Emax in 1000 steps
    // and take the best solution.

    G4double erange = params.Emax - params.Emin;
    G4double steps = erange / 1000.;
    G4int i;
    G4double etest, diff, err;

    err = 100000.;

    for (i = 1; i < 1000; i++)
    {
        etest = params.Emin + (i - 1) * steps;

        diff = etest * (std::exp(-etest / (k * Temp))) + k * Temp * (std::exp(-etest / (k * Temp))) - bigc;

        if (diff < 0.)
        {
            diff = -diff;
        }

        if (diff < err)
        {
            err = diff;
            params.particle_energy = etest;
        }
    }
    if (verbosityLevel >= 1)
    {
        G4cout << "Energy is " << params.particle_energy << G4endl;
    }
}

void G4SPSEneDistribution::GenerateBbodyEnergies()
{
    // BBody_x holds Energies, and BBHist holds the cumulative histo.
    // binary search to find correct bin then lin interpolation.
    // Use the earlier defined histogram + RandGeneral method to generate
    // random numbers following the histos distribution.
    G4double rndm;
    G4int nabove, nbelow = 0, middle;
    nabove = 10001;
    rndm = eneRndm->GenRandEnergy();
    G4AutoLock l(&mutex);
    G4bool done = histCalcd;
    l.unlock();
    if(!done)
    {
        Calculate(); //This is has a lock inside, risk is to do it twice
        l.lock();
        histCalcd = true;
        l.unlock();
    }

    // Binary search to find bin that rndm is in
    while (nabove - nbelow > 1)
    {
        middle = (nabove + nbelow) / 2;
        if (rndm == BBHist->at(middle))
        {
            break;
        }
        if (rndm < BBHist->at(middle))
        {
            nabove = middle;
        }
        else
        {
            nbelow = middle;
        }
    }

    // Now interpolate in that bin to find the correct output value.
    G4double x1, x2, y1, y2, t, q;
    x1 = Bbody_x->at(nbelow);
    if(nbelow+1 == static_cast<G4int>(Bbody_x->size()))
    {
        x2 = Bbody_x->back();
    }
    else
    {
        x2 = Bbody_x->at(nbelow + 1);
    }
    y1 = BBHist->at(nbelow);
    if(nbelow+1 == static_cast<G4int>(BBHist->size()))
    {
        G4cout << BBHist->back() << G4endl;
        y2 = BBHist->back();
    }
    else
    {
        y2 = BBHist->at(nbelow + 1);
    }
    t = (y2 - y1) / (x2 - x1);
    q = y1 - t * x1;

    threadLocalData.Get().particle_energy = (rndm - q) / t;

    if (verbosityLevel >= 1) {
        G4cout << "Energy is " << threadLocalData.Get().particle_energy << G4endl;
    }
}

void G4SPSEneDistribution::GenerateCdgEnergies() {
    // Gen random numbers, compare with values in cumhist
    // to find appropriate part of spectrum and then
    // generate energy in the usual inversion way.
    //  G4double pfact[2] = {8.5, 112};
    // G4double spind[2] = {1.4, 2.3};
    // G4double ene_line[3] = {1., 18., 1E6};
    G4double rndm, rndm2;
    G4double ene_line[3]={0,0,0};
    G4double omalpha[2]={0,0};
    threadLocal_t& params = threadLocalData.Get();
    if (params.Emin < 18 * keV && params.Emax < 18 * keV)
    {
        omalpha[0] = 1. - 1.4;
        ene_line[0] = params.Emin;
        ene_line[1] = params.Emax;
    }
    if (params.Emin < 18 * keV && params.Emax > 18 * keV)
    {
        omalpha[0] = 1. - 1.4;
        omalpha[1] = 1. - 2.3;
        ene_line[0] = params.Emin;
        ene_line[1] = 18. * keV;
        ene_line[2] = params.Emax;
    }
    if (params.Emin > 18 * keV)
    {
        omalpha[0] = 1. - 2.3;
        ene_line[0] = params.Emin;
        ene_line[1] = params.Emax;
    }
    rndm = eneRndm->GenRandEnergy();
    rndm2 = eneRndm->GenRandEnergy();

    G4int i = 0;
    while (rndm >= CDGhist[i])
    {
        i++;
    }
    // Generate final energy.
    G4double ene = (std::pow(ene_line[i - 1], omalpha[i - 1]) + (std::pow(ene_line[i], omalpha[i - 1]) - std::pow(ene_line[i - 1], omalpha[i- 1])) * rndm2);
    params.particle_energy = std::pow(ene, (1. / omalpha[i - 1]));

    if (verbosityLevel >= 1)
    {
        G4cout << "Energy is " << params.particle_energy << G4endl;
    }
}

void G4SPSEneDistribution::GenUserHistEnergies()
{
    // Histograms are DIFFERENTIAL.
    //  G4cout << "In GenUserHistEnergies " << G4endl;
    G4AutoLock l(&mutex);
    if (IPDFEnergyExist == false)
    {
        G4int ii;
        G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
        G4double bins[1024], vals[1024], sum;
        for ( ii = 0 ; ii<1024 ; ++ii ) { bins[ii]=0; vals[ii]=0; }
        sum = 0.;

        if ((EnergySpec == false) && (threadLocalData.Get().particle_definition == NULL))
        {
            G4Exception("G4SPSEneDistribution::GenUserHistEnergies",
                        "Event0302",FatalException,
                        "Error: particle definition is NULL");
        }

        if (maxbin > 1024)
        {
            G4Exception("G4SPSEneDistribution::GenUserHistEnergies",
                        "Event0302",JustWarning,"Maxbin>1024\n Setting maxbin to 1024, other bins are lost");
            maxbin = 1024;
        }

        if (DiffSpec == false)
        {
            G4cout << "Histograms are Differential!!! " << G4endl;
        }
        else
        {
            bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
            vals[0] = UDefEnergyH(size_t(0));
            sum = vals[0];
            for (ii = 1; ii < maxbin; ii++)
            {
                bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
                vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1];
                sum = sum + UDefEnergyH(size_t(ii));
            }
        }

        if (EnergySpec == false)
        {
            G4double mass = threadLocalData.Get().particle_definition->GetPDGMass();
            // multiply the function (vals) up by the bin width
            // to make the function counts/s (i.e. get rid of momentum
            // dependence).
            for (ii = 1; ii < maxbin; ii++)
            {
                vals[ii] = vals[ii] * (bins[ii] - bins[ii - 1]);
            }
            // Put energy bins into new histo, plus divide by energy bin width
            // to make evals counts/s/energy
            for (ii = 0; ii < maxbin; ii++)
            {
                bins[ii] = std::sqrt((bins[ii] * bins[ii]) + (mass * mass))- mass; //kinetic energy
            }
            for (ii = 1; ii < maxbin; ii++)
            {
                vals[ii] = vals[ii] / (bins[ii] - bins[ii - 1]);
            }
            sum = vals[maxbin - 1];
            vals[0] = 0.;
        }
        for (ii = 0; ii < maxbin; ii++)
        {
            vals[ii] = vals[ii] / sum;
            IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
        }

        // Make IPDFEnergyExist = true
        IPDFEnergyExist = true;
        if (verbosityLevel > 1)
        {
            IPDFEnergyH.DumpValues();
        }
    }
    l.unlock();
    
    // IPDF has been create so carry on
    G4double rndm = eneRndm->GenRandEnergy();
    threadLocalData.Get().particle_energy= IPDFEnergyH.GetEnergy(rndm);

    if (verbosityLevel >= 1)
    {
        G4cout << "Energy is " << particle_energy << G4endl;
    }
}

void G4SPSEneDistribution::GenArbPointEnergies()
{
    if (verbosityLevel > 0)
    {
        G4cout << "In GenArbPointEnergies" << G4endl;
    }
    G4double rndm;
    rndm = eneRndm->GenRandEnergy();
    //      IPDFArbEnergyH.DumpValues();
    // Find the Bin
    // have x, y, no of points, and cumulative area distribution
    G4int nabove, nbelow = 0, middle;
    nabove = IPDFArbEnergyH.GetVectorLength();
    //      G4cout << nabove << G4endl;
    // Binary search to find bin that rndm is in
    while (nabove - nbelow > 1)
    {
        middle = (nabove + nbelow) / 2;
        if (rndm == IPDFArbEnergyH(size_t(middle)))
        {
            break;
        }
        if (rndm < IPDFArbEnergyH(size_t(middle)))
        {
            nabove = middle;
        }
        else
        {
            nbelow = middle;
        }
    }
    threadLocal_t& params = threadLocalData.Get();
    if (IntType == "Lin")
    {
        //Update thread-local copy of parameters
        params.Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
        params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
        params.grad = Arb_grad[nbelow + 1];
        params.cept = Arb_cept[nbelow + 1];
        //    G4cout << rndm << " " << Emax << " " << Emin << " " << grad << " " << cept << G4endl;
        GenerateLinearEnergies(true);
    }
    else if (IntType == "Log")
    {
        params.Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
        params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
        params.alpha = Arb_alpha[nbelow + 1];
        //    G4cout << rndm << " " << Emax << " " << Emin << " " << alpha << G4endl;
        GeneratePowEnergies(true);
    }
    else if (IntType == "Exp")
    {
        params.Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
        params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
        params.Ezero = Arb_ezero[nbelow + 1];
        //    G4cout << rndm << " " << Emax << " " << Emin << " " << Ezero << G4endl;
        GenerateExpEnergies(true);
    }
    else if (IntType == "Spline")
    {
        params.Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
        params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
        params.particle_energy = -1e100;
        rndm = eneRndm->GenRandEnergy();
        while (params.particle_energy < params.Emin || params.particle_energy > params.Emax)
        {
            params.particle_energy = SplineInt[nbelow+1]->CubicSplineInterpolation(rndm);
            rndm = eneRndm->GenRandEnergy();
        }
        if (verbosityLevel >= 1)
        {
            G4cout << "Energy is " << params.particle_energy << G4endl;
        }
    }
    else
    {
        G4Exception("G4SPSEneDistribution::GenArbPointEnergies","Event0302",
                    FatalException,"Error: IntType unknown type");
    }
}

void G4SPSEneDistribution::GenEpnHistEnergies() {
    //  G4cout << "In GenEpnHistEnergies " << Epnflag << G4endl;

    // Firstly convert to energy if not already done.
    G4AutoLock l(&mutex);
    if (Epnflag == true)
    // epnflag = true means spectrum is epn, false means e.
    {
        // convert to energy by multiplying by A number
        ConvertEPNToEnergy();
        // EpnEnergyH will be replace by UDefEnergyH.
        //      UDefEnergyH.DumpValues();
    }
    if (IPDFEnergyExist == false)
    {
        // IPDF has not been created, so create it
        G4double bins[1024], vals[1024], sum;
        G4int ii;
        G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
        bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
        vals[0] = UDefEnergyH(size_t(0));
        sum = vals[0];
        for (ii = 1; ii < maxbin; ii++)
        {
            bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
            vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1];
            sum = sum + UDefEnergyH(size_t(ii));
        }

        l.lock();
        for (ii = 0; ii < maxbin; ii++)
        {
            vals[ii] = vals[ii] / sum;
            IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
        }
        // Make IPDFEpnExist = true
        IPDFEnergyExist = true;
       
    }
    l.unlock();
    //  IPDFEnergyH.DumpValues();
    // IPDF has been create so carry on
    G4double rndm = eneRndm->GenRandEnergy();
    threadLocalData.Get().particle_energy = IPDFEnergyH.GetEnergy(rndm);

    if (verbosityLevel >= 1)
    {
        G4cout << "Energy is " << threadLocalData.Get().particle_energy << G4endl;
    }
}

//MT: lock in caller
void G4SPSEneDistribution::ConvertEPNToEnergy()
{
    // Use this before particle generation to convert  the
    // currently stored histogram from energy/nucleon
    // to energy.
    //  G4cout << "In ConvertEpntoEnergy " << G4endl;
    threadLocal_t& params = threadLocalData.Get();
    if (params.particle_definition == NULL)
    {
        G4cout << "Error: particle not defined" << G4endl;
    }
    else
    {
        // Need to multiply histogram by the number of nucleons.
        // Baryon Number looks to hold the no. of nucleons.
        G4int Bary = params.particle_definition->GetBaryonNumber();
        //      G4cout << "Baryon No. " << Bary << G4endl;
        // Change values in histogram, Read it out, delete it, re-create it
        G4int count, maxcount;
        maxcount = G4int(EpnEnergyH.GetVectorLength());
        //      G4cout << maxcount << G4endl;
        G4double ebins[1024], evals[1024];
        if (maxcount > 1024)
        {
            G4Exception("G4SPSEneDistribution::ConvertEPNToEnergy()","gps001", JustWarning,
                        "Histogram contains more than 1024 bins!\nThose above 1024 will be ignored");
            maxcount = 1024;
        }
        if (maxcount < 1)
        {
            G4Exception("G4SPSEneDistribution::ConvertEPNToEnergy()","gps001", FatalException,"Histogram contains less than 1 bin!\nRedefine the histogram");
            return;
        }
        for (count = 0; count < maxcount; count++)
        {
            // Read out
            ebins[count] = EpnEnergyH.GetLowEdgeEnergy(size_t(count));
            evals[count] = EpnEnergyH(size_t(count));
        }

        // Multiply the channels by the nucleon number to give energies
        for (count = 0; count < maxcount; count++)
        {
            ebins[count] = ebins[count] * Bary;
        }

        // Set Emin and Emax
        params.Emin = ebins[0];
        if (maxcount > 1)
        {
            params.Emax = ebins[maxcount - 1];
        }
        else
        {
            params.Emax = ebins[0];
        }
        // Put energy bins into new histogram - UDefEnergyH.
        for (count = 0; count < maxcount; count++)
        {
            UDefEnergyH.InsertValues(ebins[count], evals[count]);
        }
        Epnflag = false; //so that you dont repeat this method.
    }
}

//
void G4SPSEneDistribution::ReSetHist(G4String atype)
{
    G4AutoLock l(&mutex);
    if (atype == "energy")
    {
        UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
        IPDFEnergyExist = false;
        Emin = 0.;
        Emax = 1e30;
    }
    else if (atype == "arb")
    {
        ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
        IPDFArbExist = false;
    }
    else if (atype == "epn")
    {
        UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
        IPDFEnergyExist = false;
        EpnEnergyH = ZeroPhysVector;
    }
    else
    {
        G4cout << "Error, histtype not accepted " << G4endl;
    }
}

G4double G4SPSEneDistribution::GenerateOne(G4ParticleDefinition* a)
{
    //Copy global shared status to thread-local one
    threadLocal_t& params = threadLocalData.Get();
    params.particle_definition=a;
    params.particle_energy=-1;
    params.Emax = Emax;
    params.Emin = Emin;
    params.alpha = alpha;
    params.Ezero = Ezero;
    params.grad = grad;
    params.cept = cept;
    params.weight = weight;
    //particle_energy = -1.;
    while ((EnergyDisType == "Arb") ? (params.particle_energy < ArbEmin || params.particle_energy > ArbEmax) : (params.particle_energy < params.Emin|| params.particle_energy > params.Emax))
    {
        if (Biased)
        {
            GenerateBiasPowEnergies();
        }
        else
        {
            if (EnergyDisType == "Mono")
            {
                GenerateMonoEnergetic();
            }
            else if (EnergyDisType == "Lin")
            {
                GenerateLinearEnergies(false);
            }
            else if (EnergyDisType == "Pow")
            {
                GeneratePowEnergies(false);
            }
            else if (EnergyDisType == "Exp")
            {
                GenerateExpEnergies(false);
            }
            else if (EnergyDisType == "Gauss")
            {
                GenerateGaussEnergies();
            }
            else if (EnergyDisType == "Brem")
            {
                GenerateBremEnergies();
            }
            else if (EnergyDisType == "Bbody")
            {
                GenerateBbodyEnergies();
            }
            else if (EnergyDisType == "Cdg")
            {
                GenerateCdgEnergies();
            }
            else if (EnergyDisType == "User")
            {
                GenUserHistEnergies();
            }
            else if (EnergyDisType == "Arb")
            {
                GenArbPointEnergies();
            }
            else if (EnergyDisType == "Epn")
            {
                GenEpnHistEnergies();
            }
            else
            {
                G4cout << "Error: EnergyDisType has unusual value" << G4endl;
            }
        }
    }
    return params.particle_energy;
}

G4double G4SPSEneDistribution::GetProbability(G4double ene)
{
    G4double prob = 1.;

    threadLocal_t& params = threadLocalData.Get();
    if (EnergyDisType == "Lin")
    {
        if (prob_norm == 1.)
        {
            prob_norm = 0.5*params.grad*params.Emax*params.Emax + params.cept*params.Emax - 0.5*params.grad*params.Emin*params.Emin - params.cept*params.Emin;
        }
        prob = params.cept + params.grad * ene;
        prob /= prob_norm;
    }
    else if (EnergyDisType == "Pow")
    {
        if (prob_norm == 1.)
        {
            if (alpha != -1.)
            {
                G4double emina = std::pow(params.Emin, params.alpha + 1);
                G4double emaxa = std::pow(params.Emax, params.alpha + 1);
                prob_norm = 1./(1.+alpha) * (emaxa - emina);
            }
            else
            {
                prob_norm = std::log(params.Emax) - std::log(params.Emin) ;
            }
        }
        prob = std::pow(ene, params.alpha)/prob_norm;
    }
    else if (EnergyDisType == "Exp")
    {
        if (prob_norm == 1.)
        {
            prob_norm = -params.Ezero*(std::exp(-params.Emax/params.Ezero) - std::exp(params.Emin/params.Ezero));
        }  
        prob = std::exp(-ene / params.Ezero);
        prob /= prob_norm;
    }
    else if (EnergyDisType == "Arb")
    {
        prob = ArbEnergyH.Value(ene);
        //  prob = ArbEInt->CubicSplineInterpolation(ene);
        //G4double deltaY;
        //prob = ArbEInt->PolynomInterpolation(ene, deltaY);
        if (prob <= 0.)
        {
            //G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << " " <<deltaY<< G4endl;
            G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << G4endl;
            prob = 1e-30;
        }
        // already normalised
    }
    else
    {
        G4cout << "Error: EnergyDisType not supported" << G4endl;
    }
       
    return prob;
}