//
// ********************************************************************
// * 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";
  MonoEnergy = 1*MeV;
  Emin = 0.;
  Emax = 1.e30;
  alpha = 0.;
  Ezero = 0.;
  SE = 0.;
  Temp = 0.;
  grad = 0.;
  cept = 0.;
  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::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::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 = Emax;
  ArbEmin = Emin;

  // 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++;
    }

  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;
      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++;
    }
  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++;
    }
  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)
  G4double  Arb_x[1024], Arb_y[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
      }
    }
  }
  //
  for(i=1;i<maxi;i++)
    Arb_y[i] += Arb_y[i-1];

  for(i=0;i<maxi;i++)
    Arb_y[i] /= Arb_y[maxi-1];
  // now Arb_y is accumulated normalised probabilities
  /*  for(i=0; i<maxi;i++) {
      if(verbosityLevel >1) 
       G4cout << i <<" "<< Arb_x[i] << " " << Arb_y[i] << G4endl;
       IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_y[i]);    
   }
   Emax = IPDFArbEnergyH.GetLowEdgeEnergy(IPDFArbEnergyH.GetVectorLength()-1);
   Emin = IPDFArbEnergyH.GetLowEdgeEnergy(0);
  */
  // Should now have normalised cumulative probabilities in Arb_y
  // and energy values in Arb_x.
  // maxi = maxi + 1;
  // Put y into x and x into y. The spline interpolation will then
  // go through x-axis to find where to interpolate (cum probability)
  // then generate a y (which will now be energy).
  SplineInt = new G4DataInterpolation(Arb_y,Arb_x,maxi,1e30,1e30);
  if(verbosityLevel >1 )
    {
      G4cout << SplineInt << G4endl;
      G4cout << SplineInt->LocateArgument(1.0) << G4endl;
    }
  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 powerlaw

  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 if(alpha == -1.)
    {
      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::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();
  if(IntType != "Spline")
    {
      //      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")
    {
      if(verbosityLevel > 1)
	G4cout << "IntType = Spline " << rndm << G4endl;
      // in SplineInterpolation created SplineInt
      // Now generate a random number put it into CubicSplineInterpolation
      // and you should get out an energy!?!
      particle_energy = -1e100;
      while (particle_energy < Emin || particle_energy > Emax ) {
	particle_energy = SplineInt->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(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;
}