source: trunk/source/event/src/G4SPSEneDistribution.cc@ 1354

//
// ********************************************************************
// * 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
//
///////////////////////////////////////////////////////////////////////////////
//
#include "Randomize.hh"
//#include <cmath>

#include "G4SPSEneDistribution.hh"

G4SPSEneDistribution::G4SPSEneDistribution() {
        //
        // 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;

}

G4SPSEneDistribution::~G4SPSEneDistribution() {
}

void G4SPSEneDistribution::SetEnergyDisType(G4String DisType) {
        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;
        }
}

void G4SPSEneDistribution::SetEmin(G4double emi) {
        Emin = emi;
}

void G4SPSEneDistribution::SetEmax(G4double ema) {
        Emax = ema;
}

void G4SPSEneDistribution::SetMonoEnergy(G4double menergy) {
        MonoEnergy = menergy;
}

void G4SPSEneDistribution::SetBeamSigmaInE(G4double e) {
        SE = e;
}
void G4SPSEneDistribution::SetAlpha(G4double alp) {
        alpha = alp;
}

void G4SPSEneDistribution::SetBiasAlpha(G4double alp) {
        biasalpha = alp;
        Biased = true;
}

void G4SPSEneDistribution::SetTemp(G4double tem) {
        Temp = tem;
}

void G4SPSEneDistribution::SetEzero(G4double eze) {
        Ezero = eze;
}

void G4SPSEneDistribution::SetGradient(G4double gr) {
        grad = gr;
}

void G4SPSEneDistribution::SetInterCept(G4double c) {
        cept = c;
}

void G4SPSEneDistribution::UserEnergyHisto(G4ThreeVector input) {
        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;
}

void G4SPSEneDistribution::ArbEnergyHisto(G4ThreeVector input) {
        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) {
        std::ifstream infile(filename, std::ios::in);
        if (!infile)
                G4Exception("Unable to open the histo ASCII file");
        G4double ehi, val;
        while (infile >> ehi >> val) {
                ArbEnergyH.InsertValues(ehi, val);
        }
}

void G4SPSEneDistribution::EpnEnergyHisto(G4ThreeVector input) {
        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;
        Epnflag = true;
}

void G4SPSEneDistribution::Calculate() {
        if (EnergyDisType == "Cdg")
                CalculateCdgSpectrum();
        else if (EnergyDisType == "Bbody")
                CalculateBbodySpectrum();
}

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] = Emin;
        if (Emin < 18 * keV) {
                n_par = 2;
                ene_line[2] = Emax;
                if (Emax < 18 * keV) {
                        n_par = 1;
                        ene_line[1] = Emax;
                }
        } else {
                n_par = 1;
                pfact[0] = 112.;
                spind[0] = 2.3;
                ene_line[1] = 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++;
        }
}

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 = Emax - Emin;
        G4double steps = erange / 10000.;
        G4double Bbody_y[10000];
        G4double k = 8.6181e-11; //Boltzmann const in MeV/K
        G4double h = 4.1362e-21; // Plancks const in MeV s
        G4double c = 3e8; // Speed of light
        G4double h2 = h * h;
        G4double c2 = c * c;
        G4int count = 0;
        G4double sum = 0.;
        BBHist[0] = 0.;
        while (count < 10000) {
                Bbody_x[count] = Emin + G4double(count * steps);
                Bbody_y[count] = (2. * std::pow(Bbody_x[count], 2.)) / (h2 * c2
                                * (std::exp(Bbody_x[count] / (k * Temp)) - 1.));
                sum = sum + Bbody_y[count];
                BBHist[count + 1] = BBHist[count] + Bbody_y[count];
                count++;
        }

        Bbody_x[10000] = Emax;
        // Normalise cumulative histo.
        count = 0;
        while (count < 10001) {
                BBHist[count] = BBHist[count] / sum;
                count++;
        }
}

void G4SPSEneDistribution::InputEnergySpectra(G4bool value) {
        // 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) {
        // 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) {
        if (EnergyDisType != "Arb")
                G4cout << "Error: this is for arbitrary distributions" << G4endl;
        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();
}

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.
                if (particle_definition == NULL)
                        G4cout << "Error: particle not defined" << G4endl;
                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 = particle_definition->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_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();
        }
}

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.
                if (particle_definition == NULL)
                        G4cout << "Error: particle not defined" << G4endl;
                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 = particle_definition->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_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;
}

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.
                if (particle_definition == NULL)
                        G4cout << "Error: particle not defined" << G4endl;
                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 = particle_definition->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_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 {
                        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;
}

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.
                if (particle_definition == NULL)
                        G4cout << "Error: particle not defined" << G4endl;
                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 = particle_definition->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];
        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.) { 
              G4cout <<   "Warning: G4DataInterpolation returns value < 0  " << prob[count] <<" "<<ei[count]<< G4endl;
              G4Exception("         Please use an alternative method, e.g. Lin, for interpolation");
            }
            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.
        particle_energy = MonoEnergy;
}

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

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

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

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

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

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

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

        if (particle_energy < 0.)
                particle_energy = -particle_energy;

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

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

        G4double rndm;
        G4double emina, emaxa;

        emina = std::pow(Emin, alpha + 1);
        emaxa = std::pow(Emax, alpha + 1);

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

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

void G4SPSEneDistribution::GenerateBiasPowEnergies() {
        // Method to generate particle energies distributed as
        // in biased power-law and calculate its weight

        G4double rndm;
        G4double emina, emaxa, emin, emax;

        G4double normal = 1. ;

        emin = Emin;
        emax = 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);
                particle_energy = ((rndm * (emaxa - emina)) + emina);
                particle_energy = std::pow(particle_energy, (1. / (biasalpha + 1.)));
                normal = 1./(1+biasalpha) * (emaxa - emina);
        } else {
                particle_energy = (std::log(emin) + rndm * (std::log(emax) - std::log(
                                emin)));
                particle_energy = std::exp(particle_energy);
                normal = std::log(emax) - std::log(emin) ;
        }
        weight = GetProbability(particle_energy) / (std::pow(particle_energy,biasalpha)/normal);

        if (verbosityLevel >= 1)
                G4cout << "Energy is " << 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();

        particle_energy = -Ezero * (std::log(rndm * (std::exp(-Emax / Ezero)
                        - std::exp(-Emin / Ezero)) + std::exp(-Emin / Ezero)));
        if (verbosityLevel >= 1)
                G4cout << "Energy is " << 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

        expmax = std::exp(-Emax / (k * Temp));
        expmin = std::exp(-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.)
                G4cout << "*****EXPMAX=0. Choose different E's or Temp" << G4endl;
        if (expmin == 0.)
                G4cout << "*****EXPMIN=0. Choose different E's or Temp" << G4endl;

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

        G4double bigc = (tempvar - k * Temp * 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 = Emax - Emin;
        G4double steps = erange / 1000.;
        G4int i;
        G4double etest, diff, err;

        err = 100000.;

        for (i = 1; i < 1000; i++) {
                etest = 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;
                        particle_energy = etest;
                }
        }
        if (verbosityLevel >= 1)
                G4cout << "Energy is " << 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();

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

        // Now interpolate in that bin to find the correct output value.
        G4double x1, x2, y1, y2, m, q;
        x1 = Bbody_x[nbelow];
        x2 = Bbody_x[nbelow + 1];
        y1 = BBHist[nbelow];
        y2 = BBHist[nbelow + 1];
        m = (y2 - y1) / (x2 - x1);
        q = y1 - m * x1;

        particle_energy = (rndm - q) / m;

        if (verbosityLevel >= 1) {
                G4cout << "Energy is " << 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];
        G4double omalpha[2];
        if (Emin < 18 * keV && Emax < 18 * keV) {
                omalpha[0] = 1. - 1.4;
                ene_line[0] = Emin;
                ene_line[1] = Emax;
        }
        if (Emin < 18 * keV && Emax > 18 * keV) {
                omalpha[0] = 1. - 1.4;
                omalpha[1] = 1. - 2.3;
                ene_line[0] = Emin;
                ene_line[1] = 18. * keV;
                ene_line[2] = Emax;
        }
        if (Emin > 18 * keV) {
                omalpha[0] = 1. - 2.3;
                ene_line[0] = Emin;
                ene_line[1] = Emax;
        }
        rndm = eneRndm->GenRandEnergy();
        rndm2 = eneRndm->GenRandEnergy();

        G4int i = 0;
        while (rndm >= CDGhist[i]) {
                i++;
        }
        // Generate final energy.
        particle_energy = (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);
        particle_energy = std::pow(particle_energy, (1. / omalpha[i - 1]));

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

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

                if ((EnergySpec == false) && (particle_definition == NULL))
                        G4cout << "Error: particle definition is NULL" << G4endl;

                if (maxbin > 1024) {
                        G4cout << "Maxbin > 1024" << G4endl;
                        G4cout << "Setting maxbin to 1024, other bins are lost" << G4endl;
                }

                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 = 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();
        }

        // IPDF has been create so carry on
        G4double rndm = eneRndm->GenRandEnergy();
        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;
        }
        if (IntType == "Lin") {
          Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
          Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
          grad = Arb_grad[nbelow + 1];
          cept = Arb_cept[nbelow + 1];
          //      G4cout << rndm << " " << Emax << " " << Emin << " " << grad << " " << cept << G4endl;
          GenerateLinearEnergies(true);
        } else if (IntType == "Log") {
          Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
          Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
          alpha = Arb_alpha[nbelow + 1];
          //      G4cout << rndm << " " << Emax << " " << Emin << " " << alpha << G4endl;
          GeneratePowEnergies(true);
        } else if (IntType == "Exp") {
          Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
          Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
          Ezero = Arb_ezero[nbelow + 1];
          //      G4cout << rndm << " " << Emax << " " << Emin << " " << Ezero << G4endl;
          GenerateExpEnergies(true);
        } else if (IntType == "Spline") {
          Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
          Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
          particle_energy = -1e100;
          rndm = eneRndm->GenRandEnergy();
          while (particle_energy < Emin || particle_energy > Emax) {
            particle_energy = SplineInt[nbelow+1]->CubicSplineInterpolation(rndm);
            rndm = eneRndm->GenRandEnergy();
          }
          if (verbosityLevel >= 1)
            G4cout << "Energy is " << particle_energy << G4endl;
        } else
                G4cout << "Error: IntType unknown type" << G4endl;
}

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

        // Firstly convert to energy if not already done.
        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();
        }

        //  G4cout << "Creating IPDFEnergy if not already done so" << G4endl;
        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));
                }

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

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

void G4SPSEneDistribution::ConvertEPNToEnergy() {
        // Use this before particle generation to convert  the
        // currently stored histogram from energy/nucleon
        // to energy.
        //  G4cout << "In ConvertEpntoEnergy " << G4endl;
        if (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 = 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) {
                        G4cout << "Histogram contains more than 1024 bins!" << G4endl;
                        G4cout << "Those above 1024 will be ignored" << G4endl;
                        maxcount = 1024;
                }
                if (maxcount < 1) {
                        G4cout << "Histogram contains less than 1 bin!" << G4endl;
                        G4cout << "Redefine the histogram" << G4endl;
                        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
                Emin = ebins[0];
                if (maxcount > 1)
                        Emax = ebins[maxcount - 1];
                else
                        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) {
        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) {
        particle_definition = a;
        particle_energy = -1.;

        while ((EnergyDisType == "Arb") ? (particle_energy < ArbEmin
                        || particle_energy > ArbEmax) : (particle_energy < Emin
                        || particle_energy > Emax)) {
                if (Biased) {
                        GenerateBiasPowEnergies();
                } else {
                        if (EnergyDisType == "Mono")
                                GenerateMonoEnergetic();
                        else if (EnergyDisType == "Lin")
                                GenerateLinearEnergies();
                        else if (EnergyDisType == "Pow")
                                GeneratePowEnergies();
                        else if (EnergyDisType == "Exp")
                                GenerateExpEnergies();
                        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 particle_energy;
}

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

        if (EnergyDisType == "Lin") {
          if (prob_norm == 1.) {
            prob_norm = 0.5*grad*Emax*Emax + cept*Emax - 0.5*grad*Emin*Emin - cept*Emin;
          }
          prob = cept + grad * ene;
          prob /= prob_norm;
        }
        else if (EnergyDisType == "Pow") {
          if (prob_norm == 1.) {
            if (alpha != -1.) {
              G4double emina = std::pow(Emin, alpha + 1);
              G4double emaxa = std::pow(Emax, alpha + 1);
              prob_norm = 1./(1.+alpha) * (emaxa - emina);
            } else {
              prob_norm = std::log(Emax) - std::log(Emin) ;
            }
          }
          prob = std::pow(ene, alpha)/prob_norm;
        }
        else if (EnergyDisType == "Exp"){
          if (prob_norm == 1.) {
            prob_norm = -Ezero*(std::exp(-Emax/Ezero) - std::exp(Emin/Ezero));
          }  
          prob = std::exp(-ene / 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;
}