lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
216a75eeb1c9d83bd2f8d612fe761690fa561bba
geant4/source/event/src/G4SPSEneDistribution.cc
T
Gabriele Cosmo 216a75eeb1 Import Geant4 8.1.0 source tree
2016-06-09 14:44:26 +02:00

1227 lines
36 KiB
C++
Raw Blame History

//
// ********************************************************************
// * 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],omalpha)-std::pow(ene_line[i],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.;
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.;
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;
}
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];
Emax = ebins[maxcount-1];
// 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;
}
Reference in New Issue View Git Blame Copy Permalink