lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
2cce0e189df02c66f78a459a9a5a5d477f738c02
geant4/source/event/src/G4SPSEneDistribution.cc
T
Gabriele Cosmo d4af681f38 Import Geant4 10.2.0 source tree
2016-06-10 14:11:04 +02:00

1859 lines
51 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
//
///////////////////////////////////////////////////////////////////////////////
//
// 26/03/2014, Andrew Green.
// Modification to used STL vectors instead of C-style arrays. This should save some space,
// particularly when the blackbody function is not used. Also moved to dynamically allocated
// memory in the LinearInterpolation, ExpInterpolation and LogInterpolation functions. Again,
// this will save space if these functions are unused.
///////////////////////////////////////////////////////////////////////////////
#include "G4SPSEneDistribution.hh"
#include "G4Exp.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"
#include "G4AutoLock.hh"
#include "G4Threading.hh"
G4SPSEneDistribution::G4SPSEneDistribution(): Epnflag(false),eneRndm(0),Splinetemp(0)
{
G4MUTEXINIT(mutex);
//
// Initialise all variables
particle_energy = 1.0 * MeV;
EnergyDisType = "Mono";
weight=1.;
MonoEnergy = 1 * MeV;
Emin = 0.;
Emax = 1.e30;
alpha = 0.;
biasalpha = 0.;
prob_norm = 1.0;
Ezero = 0.;
SE = 0.;
Temp = 0.;
grad = 0.;
cept = 0.;
Biased = false; // not biased
EnergySpec = true; // true - energy spectra, false - momentum spectra
DiffSpec = true; // true - differential spec, false integral spec
IntType = "NULL"; // Interpolation type
IPDFEnergyExist = false;
IPDFArbExist = false;
ArbEmin = 0.;
ArbEmax = 1.e30;
verbosityLevel = 0;
//AG: Set these pointers to NULL so we know if they were used...
Arb_grad = NULL;
Arb_cept = NULL;
Arb_alpha = NULL;
Arb_Const = NULL;
Arb_ezero = NULL;
Arb_grad_cept_flag = false;
Arb_alpha_Const_flag = false;
Arb_ezero_flag = false;
histInit = false;
histCalcd = false;
BBHist = NULL;
Bbody_x = NULL;
threadLocal_t& data = threadLocalData.Get();
data.Emax = Emax;
data.Emin = Emin;
data.alpha =alpha;
data.cept = cept;
data.Ezero = Ezero;
data.grad = grad;
data.particle_energy = 0;
data.particle_definition = NULL;
data.weight = weight;
}
G4SPSEneDistribution::~G4SPSEneDistribution()
{
G4MUTEXDESTROY(mutex);
if(Arb_grad_cept_flag)
{
delete[] Arb_grad;
delete[] Arb_cept;
}
if(Arb_alpha_Const_flag)
{
delete[] Arb_alpha;
delete[] Arb_Const;
}
if(Arb_ezero_flag)
{
delete[] Arb_ezero;
}
delete Bbody_x;
delete BBHist;
for ( std::vector<G4DataInterpolation*>::iterator it = SplineInt.begin() ; it!=SplineInt.end() ; ++it)
{
delete *it;
*it = 0;
}
SplineInt.clear();
}
void G4SPSEneDistribution::SetEnergyDisType(G4String DisType) {
G4AutoLock l(&mutex);
EnergyDisType = DisType;
if (EnergyDisType == "User")
{
UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
IPDFEnergyExist = false;
}
else if (EnergyDisType == "Arb")
{
ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
IPDFArbExist = false;
}
else if (EnergyDisType == "Epn")
{
UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
IPDFEnergyExist = false;
EpnEnergyH = ZeroPhysVector;
}
}
G4String G4SPSEneDistribution::GetEnergyDisType()
{
G4AutoLock l(&mutex);
return EnergyDisType;
}
void G4SPSEneDistribution::SetEmin(G4double emi) {
G4AutoLock l(&mutex);
Emin = emi;
threadLocalData.Get().Emin = Emin;
}
G4double G4SPSEneDistribution::GetEmin()
{
return threadLocalData.Get().Emin;
}
G4double G4SPSEneDistribution::GetArbEmin()
{
G4AutoLock l(&mutex);
return ArbEmin;
}
G4double G4SPSEneDistribution::GetArbEmax()
{
G4AutoLock l(&mutex);
return ArbEmax;
}
void G4SPSEneDistribution::SetEmax(G4double ema) {
G4AutoLock l(&mutex);
Emax = ema;
threadLocalData.Get().Emax = Emax;
}
G4double G4SPSEneDistribution::GetEmax()
{
return threadLocalData.Get().Emax;
}
void G4SPSEneDistribution::SetMonoEnergy(G4double menergy) {
G4AutoLock l(&mutex);
MonoEnergy = menergy;
}
void G4SPSEneDistribution::SetBeamSigmaInE(G4double e) {
G4AutoLock l(&mutex);
SE = e;
}
void G4SPSEneDistribution::SetAlpha(G4double alp) {
G4AutoLock l(&mutex);
alpha = alp;
threadLocalData.Get().alpha = alpha;
}
void G4SPSEneDistribution::SetBiasAlpha(G4double alp) {
G4AutoLock l(&mutex);
biasalpha = alp;
Biased = true;
}
void G4SPSEneDistribution::SetTemp(G4double tem) {
G4AutoLock l(&mutex);
Temp = tem;
}
void G4SPSEneDistribution::SetEzero(G4double eze) {
G4AutoLock l(&mutex);
Ezero = eze;
threadLocalData.Get().Ezero = Ezero;
}
void G4SPSEneDistribution::SetGradient(G4double gr) {
G4AutoLock l(&mutex);
grad = gr;
threadLocalData.Get().grad = grad;
}
void G4SPSEneDistribution::SetInterCept(G4double c) {
G4AutoLock l(&mutex);
cept = c;
threadLocalData.Get().cept = cept;
}
G4String G4SPSEneDistribution::GetIntType()
{
G4AutoLock l(&mutex);
return IntType;
}
void G4SPSEneDistribution::SetBiasRndm(G4SPSRandomGenerator* a)
{
G4AutoLock l(&mutex);
eneRndm = a;
}
void G4SPSEneDistribution::SetVerbosity(G4int a)
{
G4AutoLock l(&mutex);
verbosityLevel = a;
}
G4double G4SPSEneDistribution::GetWeight()
{
return threadLocalData.Get().weight;
}
G4double G4SPSEneDistribution::GetMonoEnergy()
{
G4AutoLock l(&mutex);
return MonoEnergy;
}
G4double G4SPSEneDistribution::GetSE()
{
G4AutoLock l(&mutex);
return SE;
}
G4double G4SPSEneDistribution::Getalpha()
{
return threadLocalData.Get().alpha;
}
G4double G4SPSEneDistribution::GetEzero()
{
return threadLocalData.Get().Ezero;
}
G4double G4SPSEneDistribution::GetTemp()
{
G4AutoLock l(&mutex);
return Temp;
}
G4double G4SPSEneDistribution::Getgrad()
{
return threadLocalData.Get().grad;
}
G4double G4SPSEneDistribution::Getcept()
{
return threadLocalData.Get().cept;
}
G4PhysicsOrderedFreeVector G4SPSEneDistribution::GetUserDefinedEnergyHisto()
{
G4AutoLock l(&mutex);
return UDefEnergyH;
}
G4PhysicsOrderedFreeVector G4SPSEneDistribution::GetArbEnergyHisto()
{
G4AutoLock l(&mutex);
return ArbEnergyH;
}
void G4SPSEneDistribution::UserEnergyHisto(G4ThreeVector input) {
G4AutoLock l(&mutex);
G4double ehi, val;
ehi = input.x();
val = input.y();
if (verbosityLevel > 1) {
G4cout << "In UserEnergyHisto" << G4endl;
G4cout << " " << ehi << " " << val << G4endl;
}
UDefEnergyH.InsertValues(ehi, val);
Emax = ehi;
threadLocalData.Get().Emax = Emax;
}
void G4SPSEneDistribution::ArbEnergyHisto(G4ThreeVector input)
{
G4AutoLock l(&mutex);
G4double ehi, val;
ehi = input.x();
val = input.y();
if (verbosityLevel > 1)
{
G4cout << "In ArbEnergyHisto" << G4endl;
G4cout << " " << ehi << " " << val << G4endl;
}
ArbEnergyH.InsertValues(ehi, val);
}
void G4SPSEneDistribution::ArbEnergyHistoFile(G4String filename)
{
G4AutoLock l(&mutex);
std::ifstream infile(filename, std::ios::in);
if (!infile)
{
G4Exception("G4SPSEneDistribution::ArbEnergyHistoFile","Event0301",FatalException,"Unable to open the histo ASCII file");
}
G4double ehi, val;
while (infile >> ehi >> val)
{
ArbEnergyH.InsertValues(ehi, val);
}
}
void G4SPSEneDistribution::EpnEnergyHisto(G4ThreeVector input)
{
G4AutoLock l(&mutex);
G4double ehi, val;
ehi = input.x();
val = input.y();
if (verbosityLevel > 1)
{
G4cout << "In EpnEnergyHisto" << G4endl;
G4cout << " " << ehi << " " << val << G4endl;
}
EpnEnergyH.InsertValues(ehi, val);
Emax = ehi;
threadLocalData.Get().Emax = Emax;
Epnflag = true;
}
void G4SPSEneDistribution::Calculate()
{
G4AutoLock l(&mutex);
if (EnergyDisType == "Cdg")
{
CalculateCdgSpectrum();
}
else if (EnergyDisType == "Bbody")
{
if(!histInit)
{
InitHists();
}
CalculateBbodySpectrum();
}
}
//MT: Lock in caller
void G4SPSEneDistribution::InitHists()
{
BBHist = new std::vector<G4double>(10001, 0.0);
Bbody_x = new std::vector<G4double>(10001, 0.0);
histInit = true;
}
//MT: Lock in caller
void G4SPSEneDistribution::CalculateCdgSpectrum()
{
// This uses the spectrum from The INTEGRAL Mass Model (TIMM)
// to generate a Cosmic Diffuse X/gamma ray spectrum.
G4double pfact[2] = { 8.5, 112 };
G4double spind[2] = { 1.4, 2.3 };
G4double ene_line[3] = { 1. * keV, 18. * keV, 1E6 * keV };
G4int n_par;
ene_line[0] = threadLocalData.Get().Emin;
if (threadLocalData.Get().Emin < 18 * keV)
{
n_par = 2;
ene_line[2] = threadLocalData.Get().Emax;
if (threadLocalData.Get().Emax < 18 * keV)
{
n_par = 1;
ene_line[1] = threadLocalData.Get().Emax;
}
}
else
{
n_par = 1;
pfact[0] = 112.;
spind[0] = 2.3;
ene_line[1] = threadLocalData.Get().Emax;
}
// Create a cumulative histogram.
CDGhist[0] = 0.;
G4double omalpha;
G4int i = 0;
while (i < n_par)
{
omalpha = 1. - spind[i];
CDGhist[i + 1] = CDGhist[i] + (pfact[i] / omalpha) * (std::pow(ene_line[i + 1] / keV, omalpha) - std::pow(ene_line[i] / keV,omalpha));
i++;
}
// Normalise histo and divide by 1000 to make MeV.
i = 0;
while (i < n_par)
{
CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[n_par];
// G4cout << CDGhist[i] << CDGhist[n_par] << G4endl;
i++;
}
}
//MT: Lock in caller
void G4SPSEneDistribution::CalculateBbodySpectrum()
{
// create bbody spectrum
// Proved very hard to integrate indefinitely, so different
// method. User inputs emin, emax and T. These are used to
// create a 10,000 bin histogram.
// Use photon density spectrum = 2 nu**2/c**2 * (std::exp(h nu/kT)-1)
// = 2 E**2/h**2c**2 times the exponential
G4double erange = threadLocalData.Get().Emax - threadLocalData.Get().Emin;
G4double steps = erange / 10000.;
const G4double k = 8.6181e-11; //Boltzmann const in MeV/K
const G4double h = 4.1362e-21; // Plancks const in MeV s
const G4double c = 3e8; // Speed of light
const G4double h2 = h * h;
const G4double c2 = c * c;
G4int count = 0;
G4double sum = 0.;
BBHist->at(0) = 0.;
while (count < 10000) {
Bbody_x->at(count) = threadLocalData.Get().Emin + G4double(count * steps);
G4double Bbody_y = (2. * std::pow(Bbody_x->at(count), 2.)) / (h2 * c2 * (std::exp(Bbody_x->at(count) / (k * Temp)) - 1.));
sum = sum + Bbody_y;
BBHist->at(count + 1) = BBHist->at(count) + Bbody_y;
count++;
}
Bbody_x->at(10000) = threadLocalData.Get().Emax;
// Normalise cumulative histo.
count = 0;
while (count < 10001) {
BBHist->at(count) = BBHist->at(count) / sum;
count++;
}
}
void G4SPSEneDistribution::InputEnergySpectra(G4bool value) {
G4AutoLock l(&mutex);
// Allows user to specifiy spectrum is momentum
EnergySpec = value; // false if momentum
if (verbosityLevel > 1)
G4cout << "EnergySpec has value " << EnergySpec << G4endl;
}
void G4SPSEneDistribution::InputDifferentialSpectra(G4bool value) {
G4AutoLock l(&mutex);
// Allows user to specify integral or differential spectra
DiffSpec = value; // true = differential, false = integral
if (verbosityLevel > 1)
G4cout << "Diffspec has value " << DiffSpec << G4endl;
}
void G4SPSEneDistribution::ArbInterpolate(G4String IType) {
G4AutoLock l(&mutex);
if (EnergyDisType != "Arb")
G4Exception("G4SPSEneDistribution::ArbInterpolate",
"Event0302",FatalException,
"Error: this is for arbitrary distributions");
IntType = IType;
ArbEmax = ArbEnergyH.GetMaxLowEdgeEnergy();
ArbEmin = ArbEnergyH.GetMinLowEdgeEnergy();
// Now interpolate points
if (IntType == "Lin")
LinearInterpolation();
if (IntType == "Log")
LogInterpolation();
if (IntType == "Exp")
ExpInterpolation();
if (IntType == "Spline")
SplineInterpolation();
}
//MT: Lock in caller
void G4SPSEneDistribution::LinearInterpolation() {
// Method to do linear interpolation on the Arb points
// Calculate equation of each line segment, max 1024.
// Calculate Area under each segment
// Create a cumulative array which is then normalised Arb_Cum_Area
G4double Area_seg[1024]; // Stores area under each segment
G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
G4int i, count;
G4int maxi = ArbEnergyH.GetVectorLength();
for (i = 0; i < maxi; i++)
{
Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
Arb_y[i] = ArbEnergyH(size_t(i));
}
// Points are now in x,y arrays. If the spectrum is integral it has to be
// made differential and if momentum it has to be made energy.
if (DiffSpec == false)
{
// Converts integral point-wise spectra to Differential
for (count = 0; count < maxi - 1; count++)
{
Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
/ (Arb_x[count + 1] - Arb_x[count]);
}
maxi--;
}
//
if (EnergySpec == false)
{
// change currently stored values (emin etc) which are actually momenta
// to energies.
G4ParticleDefinition* pdef =threadLocalData.Get().particle_definition;
if (pdef == NULL)
{
G4Exception("G4SPSEneDistribution::LinearInterpolation",
"Event0302",FatalException,
"Error: particle not defined");
}
else
{
// Apply Energy**2 = p**2c**2 + m0**2c**4
// p should be entered as E/c i.e. without the division by c
// being done - energy equivalent.
G4double mass = pdef->GetPDGMass();
// convert point to energy unit and its value to per energy unit
G4double total_energy;
for (count = 0; count < maxi; count++)
{
total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
* mass)); // total energy
Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
Arb_x[count] = total_energy - mass; // kinetic energy
}
}
}
//
i = 1;
// AG:
// This *should* be the first use of Arb_grad and friends, so we *should* be able to
// dynamically allocate here
Arb_grad = new G4double [1024];
Arb_cept = new G4double [1024];
Arb_grad_cept_flag = true;
Arb_grad[0] = 0.;
Arb_cept[0] = 0.;
Area_seg[0] = 0.;
Arb_Cum_Area[0] = 0.;
while (i < maxi)
{
// calc gradient and intercept for each segment
Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[i - 1]);
if (verbosityLevel == 2)
{
G4cout << Arb_grad[i] << G4endl;
}
if (Arb_grad[i] > 0.)
{
if (verbosityLevel == 2)
{
G4cout << "Arb_grad is positive" << G4endl;
}
Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]);
}
else if (Arb_grad[i] < 0.)
{
if (verbosityLevel == 2)
{
G4cout << "Arb_grad is negative" << G4endl;
}
Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]);
}
else
{
if (verbosityLevel == 2)
{
G4cout << "Arb_grad is 0." << G4endl;
}
Arb_cept[i] = Arb_y[i];
}
Area_seg[i] = ((Arb_grad[i] / 2) * (Arb_x[i] * Arb_x[i] - Arb_x[i - 1] * Arb_x[i - 1]) + Arb_cept[i] * (Arb_x[i] - Arb_x[i - 1]));
Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
sum = sum + Area_seg[i];
if (verbosityLevel == 2)
{
G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i] << G4endl;
}
i++;
}
i = 0;
while (i < maxi)
{
Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation
IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
i++;
}
// now scale the ArbEnergyH, needed by Probability()
ArbEnergyH.ScaleVector(1., 1./sum);
if (verbosityLevel >= 1)
{
G4cout << "Leaving LinearInterpolation" << G4endl;
ArbEnergyH.DumpValues();
IPDFArbEnergyH.DumpValues();
}
}
//MT: Lock in caller
void G4SPSEneDistribution::LogInterpolation()
{
// Interpolation based on Logarithmic equations
// Generate equations of line segments
// y = Ax**alpha => log y = alpha*logx + logA
// Find area under line segments
// create normalised, cumulative array Arb_Cum_Area
G4double Area_seg[1024]; // Stores area under each segment
G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
G4int i, count;
G4int maxi = ArbEnergyH.GetVectorLength();
for (i = 0; i < maxi; i++)
{
Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
Arb_y[i] = ArbEnergyH(size_t(i));
}
// Points are now in x,y arrays. If the spectrum is integral it has to be
// made differential and if momentum it has to be made energy.
if (DiffSpec == false)
{
// Converts integral point-wise spectra to Differential
for (count = 0; count < maxi - 1; count++)
{
Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) / (Arb_x[count + 1] - Arb_x[count]);
}
maxi--;
}
//
if (EnergySpec == false)
{
// change currently stored values (emin etc) which are actually momenta
// to energies.
G4ParticleDefinition* pdef = threadLocalData.Get().particle_definition;
if (pdef == NULL)
{
G4Exception("G4SPSEneDistribution::LogInterpolation",
"Event0302",FatalException,
"Error: particle not defined");
}
else
{
// Apply Energy**2 = p**2c**2 + m0**2c**4
// p should be entered as E/c i.e. without the division by c
// being done - energy equivalent.
G4double mass = pdef->GetPDGMass();
// convert point to energy unit and its value to per energy unit
G4double total_energy;
for (count = 0; count < maxi; count++)
{
total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
* mass)); // total energy
Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
Arb_x[count] = total_energy - mass; // kinetic energy
}
}
}
//
i = 1;
//AG: *Should* be the first use of these two arrays
if ( Arb_ezero ) { delete[] Arb_ezero; Arb_ezero = 0; }
if ( Arb_Const ) { delete[] Arb_Const; Arb_Const = 0; }
Arb_alpha = new G4double [1024];
Arb_Const = new G4double [1024];
Arb_alpha_Const_flag = true;
Arb_alpha[0] = 0.;
Arb_Const[0] = 0.;
Area_seg[0] = 0.;
Arb_Cum_Area[0] = 0.;
if (Arb_x[0] <= 0. || Arb_y[0] <= 0.)
{
G4cout << "You should not use log interpolation with points <= 0." << G4endl;
G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl;
if (Arb_x[0] <= 0.)
{
Arb_x[0] = 1e-20;
}
if (Arb_y[0] <= 0.)
{
Arb_y[0] = 1e-20;
}
}
G4double alp;
while (i < maxi)
{
// Incase points are negative or zero
if (Arb_x[i] <= 0. || Arb_y[i] <= 0.)
{
G4cout << "You should not use log interpolation with points <= 0." << G4endl;
G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl;
if (Arb_x[i] <= 0.)
{
Arb_x[i] = 1e-20;
}
if (Arb_y[i] <= 0.)
{
Arb_y[i] = 1e-20;
}
}
Arb_alpha[i] = (std::log10(Arb_y[i]) - std::log10(Arb_y[i - 1])) / (std::log10(Arb_x[i]) - std::log10(Arb_x[i - 1]));
Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[i], Arb_alpha[i]));
alp = Arb_alpha[i] + 1;
if (alp == 0.)
{
Area_seg[i] = Arb_Const[i] * (std::log(Arb_x[i]) - std::log(Arb_x[i - 1]));
}
else
{
Area_seg[i] = (Arb_Const[i] / alp) * (std::pow(Arb_x[i], alp) - std::pow(Arb_x[i - 1], alp));
}
sum = sum + Area_seg[i];
Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
if (verbosityLevel == 2)
{
G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl;
}
i++;
}
i = 0;
while (i < maxi)
{
Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
i++;
}
// now scale the ArbEnergyH, needed by Probability()
ArbEnergyH.ScaleVector(1., 1./sum);
if (verbosityLevel >= 1)
{
G4cout << "Leaving LogInterpolation " << G4endl;
}
}
//MT: Lock in caller
void G4SPSEneDistribution::ExpInterpolation()
{
// Interpolation based on Exponential equations
// Generate equations of line segments
// y = Ae**-(x/e0) => ln y = -x/e0 + lnA
// Find area under line segments
// create normalised, cumulative array Arb_Cum_Area
G4double Area_seg[1024]; // Stores area under each segment
G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
G4int i, count;
G4int maxi = ArbEnergyH.GetVectorLength();
for (i = 0; i < maxi; i++)
{
Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
Arb_y[i] = ArbEnergyH(size_t(i));
}
// Points are now in x,y arrays. If the spectrum is integral it has to be
// made differential and if momentum it has to be made energy.
if (DiffSpec == false)
{
// Converts integral point-wise spectra to Differential
for (count = 0; count < maxi - 1; count++)
{
Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
/ (Arb_x[count + 1] - Arb_x[count]);
}
maxi--;
}
//
if (EnergySpec == false)
{
// change currently stored values (emin etc) which are actually momenta
// to energies.
G4ParticleDefinition* pdef = threadLocalData.Get().particle_definition;
if (pdef == NULL)
{
G4Exception("G4SPSEneDistribution::ExpInterpolation",
"Event0302",FatalException,
"Error: particle not defined");
}
else
{
// Apply Energy**2 = p**2c**2 + m0**2c**4
// p should be entered as E/c i.e. without the division by c
// being done - energy equivalent.
G4double mass = pdef->GetPDGMass();
// convert point to energy unit and its value to per energy unit
G4double total_energy;
for (count = 0; count < maxi; count++)
{
total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass * mass)); // total energy
Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
Arb_x[count] = total_energy - mass; // kinetic energy
}
}
}
//
i = 1;
//AG: Should be first use...
if ( Arb_ezero ) { delete[] Arb_ezero; Arb_ezero = 0; }
if ( Arb_Const ) { delete[] Arb_Const; Arb_Const = 0; }
Arb_ezero = new G4double [1024];
Arb_Const = new G4double [1024];
Arb_ezero_flag = true;
Arb_ezero[0] = 0.;
Arb_Const[0] = 0.;
Area_seg[0] = 0.;
Arb_Cum_Area[0] = 0.;
while (i < maxi)
{
G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i - 1]);
if (test > 0. || test < 0.)
{
Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1]) / (std::log(Arb_y[i]) - std::log(Arb_y[i - 1]));
Arb_Const[i] = Arb_y[i] / (std::exp(-Arb_x[i] / Arb_ezero[i]));
Area_seg[i] = -(Arb_Const[i] * Arb_ezero[i]) * (std::exp(-Arb_x[i] / Arb_ezero[i]) - std::exp(-Arb_x[i - 1] / Arb_ezero[i]));
}
else
{
G4Exception("G4SPSEneDistribution::ExpInterpolation",
"Event0302",JustWarning,
"Flat line segment: problem, setting to zero parameters.");
G4cout << "Flat line segment: problem" << G4endl;
Arb_ezero[i] = 0.;
Arb_Const[i] = 0.;
Area_seg[i] = 0.;
}
sum = sum + Area_seg[i];
Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
if (verbosityLevel == 2)
{
G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl;
}
i++;
}
i = 0;
while (i < maxi)
{
Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
i++;
}
// now scale the ArbEnergyH, needed by Probability()
ArbEnergyH.ScaleVector(1., 1./sum);
if (verbosityLevel >= 1)
{
G4cout << "Leaving ExpInterpolation " << G4endl;
}
}
//MT: Lock in caller
void G4SPSEneDistribution::SplineInterpolation()
{
// Interpolation using Splines.
// Create Normalised arrays, make x 0->1 and y hold
// the function (Energy)
//
// Current method based on the above will not work in all cases.
// new method is implemented below.
G4double sum, Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
G4int i, count;
G4int maxi = ArbEnergyH.GetVectorLength();
for (i = 0; i < maxi; i++) {
Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
Arb_y[i] = ArbEnergyH(size_t(i));
}
// Points are now in x,y arrays. If the spectrum is integral it has to be
// made differential and if momentum it has to be made energy.
if (DiffSpec == false) {
// Converts integral point-wise spectra to Differential
for (count = 0; count < maxi - 1; count++)
{
Arb_y[count] = (Arb_y[count] - Arb_y[count + 1]) / (Arb_x[count + 1] - Arb_x[count]);
}
maxi--;
}
//
if (EnergySpec == false) {
// change currently stored values (emin etc) which are actually momenta
// to energies.
G4ParticleDefinition* pdef = threadLocalData.Get().particle_definition;
if (pdef == NULL)
G4Exception("G4SPSEneDistribution::SplineInterpolation",
"Event0302",FatalException,
"Error: particle not defined");
else {
// Apply Energy**2 = p**2c**2 + m0**2c**4
// p should be entered as E/c i.e. without the division by c
// being done - energy equivalent.
G4double mass = pdef->GetPDGMass();
// convert point to energy unit and its value to per energy unit
G4double total_energy;
for (count = 0; count < maxi; count++) {
total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
* mass)); // total energy
Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
Arb_x[count] = total_energy - mass; // kinetic energy
}
}
}
//
i = 1;
Arb_Cum_Area[0] = 0.;
sum = 0.;
Splinetemp = new G4DataInterpolation(Arb_x, Arb_y, maxi, 0., 0.);
G4double ei[101],prob[101];
for ( std::vector<G4DataInterpolation*>::iterator it = SplineInt.begin() ; it!=SplineInt.end() ; ++it)
{
delete *it;
*it = 0;
}
SplineInt.clear();
SplineInt.resize(1024,NULL);
while (i < maxi) {
// 100 step per segment for the integration of area
G4double de = (Arb_x[i] - Arb_x[i - 1])/100.;
G4double area = 0.;
for (count = 0; count < 101; count++) {
ei[count] = Arb_x[i - 1] + de*count ;
prob[count] = Splinetemp->CubicSplineInterpolation(ei[count]);
if (prob[count] < 0.) {
G4ExceptionDescription ED;
ED << "Warning: G4DataInterpolation returns value < 0 " << prob[count] <<" "<<ei[count]<< G4endl;
G4Exception("G4SPSEneDistribution::SplineInterpolation","Event0303",
FatalException,ED);
}
area += prob[count]*de;
}
Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + area;
sum += area;
prob[0] = prob[0]/(area/de);
for (count = 1; count < 100; count++)
prob[count] = prob[count-1] + prob[count]/(area/de);
SplineInt[i]=new G4DataInterpolation(prob, ei, 101, 0., 0.);
// note i start from 1!
i++;
}
i = 0;
while (i < maxi) {
Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation
IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
i++;
}
// now scale the ArbEnergyH, needed by Probability()
ArbEnergyH.ScaleVector(1., 1./sum);
if (verbosityLevel > 0)
G4cout << "Leaving SplineInterpolation " << G4endl;
}
void G4SPSEneDistribution::GenerateMonoEnergetic() {
// Method to generate MonoEnergetic particles.
threadLocalData.Get().particle_energy = MonoEnergy;
}
void G4SPSEneDistribution::GenerateGaussEnergies() {
// Method to generate Gaussian particles.
G4double ene = G4RandGauss::shoot(MonoEnergy,SE);
if (ene < 0) ene = 0.;
threadLocalData.Get().particle_energy = ene;
}
void G4SPSEneDistribution::GenerateLinearEnergies(G4bool bArb = false) {
G4double rndm;
threadLocal_t& params = threadLocalData.Get();
G4double emaxsq = std::pow(params.Emax, 2.); //Emax squared
G4double eminsq = std::pow(params.Emin, 2.); //Emin squared
G4double intersq = std::pow(params.cept, 2.); //cept squared
if (bArb)
rndm = G4UniformRand();
else
rndm = eneRndm->GenRandEnergy();
G4double bracket = ((params.grad / 2.) * (emaxsq - eminsq) + params.cept * (params.Emax - params.Emin));
bracket = bracket * rndm;
bracket = bracket + (params.grad / 2.) * eminsq + params.cept * params.Emin;
// Now have a quad of form m/2 E**2 + cE - bracket = 0
bracket = -bracket;
// G4cout << "BRACKET" << bracket << G4endl;
if (params.grad != 0.)
{
G4double sqbrack = (intersq - 4 * (params.grad / 2.) * (bracket));
// G4cout << "SQBRACK" << sqbrack << G4endl;
sqbrack = std::sqrt(sqbrack);
G4double root1 = -params.cept + sqbrack;
root1 = root1 / (2. * (params.grad / 2.));
G4double root2 = -params.cept - sqbrack;
root2 = root2 / (2. * (params.grad / 2.));
// G4cout << root1 << " roots " << root2 << G4endl;
if (root1 > params.Emin && root1 < params.Emax)
{
params.particle_energy = root1;
}
if (root2 > params.Emin && root2 < params.Emax)
{
params.particle_energy = root2;
}
}
else if (params.grad == 0.)
{
// have equation of form cE - bracket =0
params.particle_energy = bracket / params.cept;
}
if (params.particle_energy < 0.)
{
params.particle_energy = -params.particle_energy;
}
if (verbosityLevel >= 1)
{
G4cout << "Energy is " << params.particle_energy << G4endl;
}
}
void G4SPSEneDistribution::GeneratePowEnergies(G4bool bArb = false)
{
// Method to generate particle energies distributed as
// a power-law
G4double rndm;
G4double emina, emaxa;
threadLocal_t& params = threadLocalData.Get();
emina = std::pow(params.Emin, params.alpha + 1);
emaxa = std::pow(params.Emax, params.alpha + 1);
if (bArb)
{
rndm = G4UniformRand();
}
else
{
rndm = eneRndm->GenRandEnergy();
}
if (params.alpha != -1.)
{
G4double ene = ((rndm * (emaxa - emina)) + emina);
ene = std::pow(ene, (1. / (params.alpha + 1.)));
params.particle_energy = ene;
}
else
{
G4double ene = (std::log(params.Emin) + rndm * (std::log(params.Emax) - std::log(params.Emin)));
params.particle_energy = std::exp(ene);
}
if (verbosityLevel >= 1)
{
G4cout << "Energy is " << params.particle_energy << G4endl;
}
}
void G4SPSEneDistribution::GenerateBiasPowEnergies()//G4double& ene,G4double& wweight)
{
// Method to generate particle energies distributed as
// in biased power-law and calculate its weight
threadLocal_t& params = threadLocalData.Get();
G4double rndm;
G4double emina, emaxa, emin, emax;
G4double normal = 1.;
emin = params.Emin;
emax = params.Emax;
// if (EnergyDisType == "Arb") {
// emin = ArbEmin;
// emax = ArbEmax;
//}
rndm = eneRndm->GenRandEnergy();
if (biasalpha != -1.)
{
emina = std::pow(emin, biasalpha + 1);
emaxa = std::pow(emax, biasalpha + 1);
G4double ee = ((rndm * (emaxa - emina)) + emina);
params.particle_energy = std::pow(ee, (1. / (biasalpha + 1.)));
normal = 1./(1+biasalpha) * (emaxa - emina);
}
else
{
G4double ee = (std::log(emin) + rndm * (std::log(emax) - std::log(emin)));
params.particle_energy = std::exp(ee);
normal = std::log(emax) - std::log(emin) ;
}
params.weight = GetProbability(params.particle_energy)
/ (std::pow(params.particle_energy,biasalpha)/normal);
if (verbosityLevel >= 1)
{
G4cout << "Energy is " << params.particle_energy << G4endl;
}
}
void G4SPSEneDistribution::GenerateExpEnergies(G4bool bArb = false)
{
// Method to generate particle energies distributed according
// to an exponential curve.
G4double rndm;
if (bArb)
{
rndm = G4UniformRand();
}
else
{
rndm = eneRndm->GenRandEnergy();
}
threadLocal_t& params = threadLocalData.Get();
params.particle_energy = -params.Ezero * (std::log(rndm * (std::exp(-params.Emax / params.Ezero) - std::exp(-params.Emin / params.Ezero)) + std::exp(-params.Emin / params.Ezero)));
if (verbosityLevel >= 1)
{
G4cout << "Energy is " << params.particle_energy << G4endl;
}
}
void G4SPSEneDistribution::GenerateBremEnergies()
{
// Method to generate particle energies distributed according
// to a Bremstrahlung equation of
// form I = const*((kT)**1/2)*E*(e**(-E/kT))
G4double rndm;
rndm = eneRndm->GenRandEnergy();
G4double expmax, expmin, k;
k = 8.6181e-11; // Boltzmann's const in MeV/K
G4double ksq = std::pow(k, 2.); // k squared
G4double Tsq = std::pow(Temp, 2.); // Temp squared
threadLocal_t& params = threadLocalData.Get();
expmax = std::exp(-params.Emax / (k * Temp));
expmin = std::exp(-params.Emin / (k * Temp));
// If either expmax or expmin are zero then this will cause problems
// Most probably this will be because T is too low or E is too high
if (expmax == 0.)
{
G4Exception("G4SPSEneDistribution::GenerateBremEnergies",
"Event0302",FatalException,
"*****EXPMAX=0. Choose different E's or Temp");
}
if (expmin == 0.)
{
G4Exception("G4SPSEneDistribution::GenerateBremEnergies",
"Event0302",FatalException,
"*****EXPMIN=0. Choose different E's or Temp");
}
G4double tempvar = rndm * ((-k) * Temp * (params.Emax * expmax - params.Emin * expmin) - (ksq * Tsq * (expmax - expmin)));
G4double bigc = (tempvar - k * Temp * params.Emin * expmin - ksq * Tsq * expmin) / (-k * Temp);
// This gives an equation of form: Ee(-E/kT) + kTe(-E/kT) - C =0
// Solve this iteratively, step from Emin to Emax in 1000 steps
// and take the best solution.
G4double erange = params.Emax - params.Emin;
G4double steps = erange / 1000.;
G4int i;
G4double etest, diff, err;
err = 100000.;
for (i = 1; i < 1000; i++)
{
etest = params.Emin + (i - 1) * steps;
diff = etest * (std::exp(-etest / (k * Temp))) + k * Temp * (std::exp(-etest / (k * Temp))) - bigc;
if (diff < 0.)
{
diff = -diff;
}
if (diff < err)
{
err = diff;
params.particle_energy = etest;
}
}
if (verbosityLevel >= 1)
{
G4cout << "Energy is " << params.particle_energy << G4endl;
}
}
void G4SPSEneDistribution::GenerateBbodyEnergies()
{
// BBody_x holds Energies, and BBHist holds the cumulative histo.
// binary search to find correct bin then lin interpolation.
// Use the earlier defined histogram + RandGeneral method to generate
// random numbers following the histos distribution.
G4double rndm;
G4int nabove, nbelow = 0, middle;
nabove = 10001;
rndm = eneRndm->GenRandEnergy();
G4AutoLock l(&mutex);
G4bool done = histCalcd;
l.unlock();
if(!done)
{
Calculate(); //This is has a lock inside, risk is to do it twice
l.lock();
histCalcd = true;
l.unlock();
}
// Binary search to find bin that rndm is in
while (nabove - nbelow > 1)
{
middle = (nabove + nbelow) / 2;
if (rndm == BBHist->at(middle))
{
break;
}
if (rndm < BBHist->at(middle))
{
nabove = middle;
}
else
{
nbelow = middle;
}
}
// Now interpolate in that bin to find the correct output value.
G4double x1, x2, y1, y2, t, q;
x1 = Bbody_x->at(nbelow);
if(nbelow+1 == static_cast<G4int>(Bbody_x->size()))
{
x2 = Bbody_x->back();
}
else
{
x2 = Bbody_x->at(nbelow + 1);
}
y1 = BBHist->at(nbelow);
if(nbelow+1 == static_cast<G4int>(BBHist->size()))
{
G4cout << BBHist->back() << G4endl;
y2 = BBHist->back();
}
else
{
y2 = BBHist->at(nbelow + 1);
}
t = (y2 - y1) / (x2 - x1);
q = y1 - t * x1;
threadLocalData.Get().particle_energy = (rndm - q) / t;
if (verbosityLevel >= 1) {
G4cout << "Energy is " << threadLocalData.Get().particle_energy << G4endl;
}
}
void G4SPSEneDistribution::GenerateCdgEnergies() {
// Gen random numbers, compare with values in cumhist
// to find appropriate part of spectrum and then
// generate energy in the usual inversion way.
// G4double pfact[2] = {8.5, 112};
// G4double spind[2] = {1.4, 2.3};
// G4double ene_line[3] = {1., 18., 1E6};
G4double rndm, rndm2;
G4double ene_line[3]={0,0,0};
G4double omalpha[2]={0,0};
threadLocal_t& params = threadLocalData.Get();
if (params.Emin < 18 * keV && params.Emax < 18 * keV)
{
omalpha[0] = 1. - 1.4;
ene_line[0] = params.Emin;
ene_line[1] = params.Emax;
}
if (params.Emin < 18 * keV && params.Emax > 18 * keV)
{
omalpha[0] = 1. - 1.4;
omalpha[1] = 1. - 2.3;
ene_line[0] = params.Emin;
ene_line[1] = 18. * keV;
ene_line[2] = params.Emax;
}
if (params.Emin > 18 * keV)
{
omalpha[0] = 1. - 2.3;
ene_line[0] = params.Emin;
ene_line[1] = params.Emax;
}
rndm = eneRndm->GenRandEnergy();
rndm2 = eneRndm->GenRandEnergy();
G4int i = 0;
while (rndm >= CDGhist[i])
{
i++;
}
// Generate final energy.
G4double ene = (std::pow(ene_line[i - 1], omalpha[i - 1]) + (std::pow(ene_line[i], omalpha[i - 1]) - std::pow(ene_line[i - 1], omalpha[i- 1])) * rndm2);
params.particle_energy = std::pow(ene, (1. / omalpha[i - 1]));
if (verbosityLevel >= 1)
{
G4cout << "Energy is " << params.particle_energy << G4endl;
}
}
void G4SPSEneDistribution::GenUserHistEnergies()
{
// Histograms are DIFFERENTIAL.
// G4cout << "In GenUserHistEnergies " << G4endl;
G4AutoLock l(&mutex);
if (IPDFEnergyExist == false)
{
G4int ii;
G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
G4double bins[1024], vals[1024], sum;
for ( ii = 0 ; ii<1024 ; ++ii ) { bins[ii]=0; vals[ii]=0; }
sum = 0.;
if ((EnergySpec == false) && (threadLocalData.Get().particle_definition == NULL))
{
G4Exception("G4SPSEneDistribution::GenUserHistEnergies",
"Event0302",FatalException,
"Error: particle definition is NULL");
}
if (maxbin > 1024)
{
G4Exception("G4SPSEneDistribution::GenUserHistEnergies",
"Event0302",JustWarning,"Maxbin>1024\n Setting maxbin to 1024, other bins are lost");
maxbin = 1024;
}
if (DiffSpec == false)
{
G4cout << "Histograms are Differential!!! " << G4endl;
}
else
{
bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
vals[0] = UDefEnergyH(size_t(0));
sum = vals[0];
for (ii = 1; ii < maxbin; ii++)
{
bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1];
sum = sum + UDefEnergyH(size_t(ii));
}
}
if (EnergySpec == false)
{
G4double mass = threadLocalData.Get().particle_definition->GetPDGMass();
// multiply the function (vals) up by the bin width
// to make the function counts/s (i.e. get rid of momentum
// dependence).
for (ii = 1; ii < maxbin; ii++)
{
vals[ii] = vals[ii] * (bins[ii] - bins[ii - 1]);
}
// Put energy bins into new histo, plus divide by energy bin width
// to make evals counts/s/energy
for (ii = 0; ii < maxbin; ii++)
{
bins[ii] = std::sqrt((bins[ii] * bins[ii]) + (mass * mass))- mass; //kinetic energy
}
for (ii = 1; ii < maxbin; ii++)
{
vals[ii] = vals[ii] / (bins[ii] - bins[ii - 1]);
}
sum = vals[maxbin - 1];
vals[0] = 0.;
}
for (ii = 0; ii < maxbin; ii++)
{
vals[ii] = vals[ii] / sum;
IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
}
// Make IPDFEnergyExist = true
IPDFEnergyExist = true;
if (verbosityLevel > 1)
{
IPDFEnergyH.DumpValues();
}
}
l.unlock();
// IPDF has been create so carry on
G4double rndm = eneRndm->GenRandEnergy();
threadLocalData.Get().particle_energy= IPDFEnergyH.GetEnergy(rndm);
if (verbosityLevel >= 1)
{
G4cout << "Energy is " << particle_energy << G4endl;
}
}
void G4SPSEneDistribution::GenArbPointEnergies()
{
if (verbosityLevel > 0)
{
G4cout << "In GenArbPointEnergies" << G4endl;
}
G4double rndm;
rndm = eneRndm->GenRandEnergy();
// IPDFArbEnergyH.DumpValues();
// Find the Bin
// have x, y, no of points, and cumulative area distribution
G4int nabove, nbelow = 0, middle;
nabove = IPDFArbEnergyH.GetVectorLength();
// G4cout << nabove << G4endl;
// Binary search to find bin that rndm is in
while (nabove - nbelow > 1)
{
middle = (nabove + nbelow) / 2;
if (rndm == IPDFArbEnergyH(size_t(middle)))
{
break;
}
if (rndm < IPDFArbEnergyH(size_t(middle)))
{
nabove = middle;
}
else
{
nbelow = middle;
}
}
threadLocal_t& params = threadLocalData.Get();
if (IntType == "Lin")
{
//Update thread-local copy of parameters
params.Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
params.grad = Arb_grad[nbelow + 1];
params.cept = Arb_cept[nbelow + 1];
// G4cout << rndm << " " << Emax << " " << Emin << " " << grad << " " << cept << G4endl;
GenerateLinearEnergies(true);
}
else if (IntType == "Log")
{
params.Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
params.alpha = Arb_alpha[nbelow + 1];
// G4cout << rndm << " " << Emax << " " << Emin << " " << alpha << G4endl;
GeneratePowEnergies(true);
}
else if (IntType == "Exp")
{
params.Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
params.Ezero = Arb_ezero[nbelow + 1];
// G4cout << rndm << " " << Emax << " " << Emin << " " << Ezero << G4endl;
GenerateExpEnergies(true);
}
else if (IntType == "Spline")
{
params.Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
params.Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
params.particle_energy = -1e100;
rndm = eneRndm->GenRandEnergy();
while (params.particle_energy < params.Emin || params.particle_energy > params.Emax)
{
params.particle_energy = SplineInt[nbelow+1]->CubicSplineInterpolation(rndm);
rndm = eneRndm->GenRandEnergy();
}
if (verbosityLevel >= 1)
{
G4cout << "Energy is " << params.particle_energy << G4endl;
}
}
else
{
G4Exception("G4SPSEneDistribution::GenArbPointEnergies","Event0302",
FatalException,"Error: IntType unknown type");
}
}
void G4SPSEneDistribution::GenEpnHistEnergies() {
// G4cout << "In GenEpnHistEnergies " << Epnflag << G4endl;
// Firstly convert to energy if not already done.
G4AutoLock l(&mutex);
if (Epnflag == true)
// epnflag = true means spectrum is epn, false means e.
{
// convert to energy by multiplying by A number
ConvertEPNToEnergy();
// EpnEnergyH will be replace by UDefEnergyH.
// UDefEnergyH.DumpValues();
}
if (IPDFEnergyExist == false)
{
// IPDF has not been created, so create it
G4double bins[1024], vals[1024], sum;
G4int ii;
G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
vals[0] = UDefEnergyH(size_t(0));
sum = vals[0];
for (ii = 1; ii < maxbin; ii++)
{
bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1];
sum = sum + UDefEnergyH(size_t(ii));
}
l.lock();
for (ii = 0; ii < maxbin; ii++)
{
vals[ii] = vals[ii] / sum;
IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
}
// Make IPDFEpnExist = true
IPDFEnergyExist = true;
}
l.unlock();
// IPDFEnergyH.DumpValues();
// IPDF has been create so carry on
G4double rndm = eneRndm->GenRandEnergy();
threadLocalData.Get().particle_energy = IPDFEnergyH.GetEnergy(rndm);
if (verbosityLevel >= 1)
{
G4cout << "Energy is " << threadLocalData.Get().particle_energy << G4endl;
}
}
//MT: lock in caller
void G4SPSEneDistribution::ConvertEPNToEnergy()
{
// Use this before particle generation to convert the
// currently stored histogram from energy/nucleon
// to energy.
// G4cout << "In ConvertEpntoEnergy " << G4endl;
threadLocal_t& params = threadLocalData.Get();
if (params.particle_definition == NULL)
{
G4cout << "Error: particle not defined" << G4endl;
}
else
{
// Need to multiply histogram by the number of nucleons.
// Baryon Number looks to hold the no. of nucleons.
G4int Bary = params.particle_definition->GetBaryonNumber();
// G4cout << "Baryon No. " << Bary << G4endl;
// Change values in histogram, Read it out, delete it, re-create it
G4int count, maxcount;
maxcount = G4int(EpnEnergyH.GetVectorLength());
// G4cout << maxcount << G4endl;
G4double ebins[1024], evals[1024];
if (maxcount > 1024)
{
G4Exception("G4SPSEneDistribution::ConvertEPNToEnergy()","gps001", JustWarning,
"Histogram contains more than 1024 bins!\nThose above 1024 will be ignored");
maxcount = 1024;
}
if (maxcount < 1)
{
G4Exception("G4SPSEneDistribution::ConvertEPNToEnergy()","gps001", FatalException,"Histogram contains less than 1 bin!\nRedefine the histogram");
return;
}
for (count = 0; count < maxcount; count++)
{
// Read out
ebins[count] = EpnEnergyH.GetLowEdgeEnergy(size_t(count));
evals[count] = EpnEnergyH(size_t(count));
}
// Multiply the channels by the nucleon number to give energies
for (count = 0; count < maxcount; count++)
{
ebins[count] = ebins[count] * Bary;
}
// Set Emin and Emax
params.Emin = ebins[0];
if (maxcount > 1)
{
params.Emax = ebins[maxcount - 1];
}
else
{
params.Emax = ebins[0];
}
// Put energy bins into new histogram - UDefEnergyH.
for (count = 0; count < maxcount; count++)
{
UDefEnergyH.InsertValues(ebins[count], evals[count]);
}
Epnflag = false; //so that you dont repeat this method.
}
}
//
void G4SPSEneDistribution::ReSetHist(G4String atype)
{
G4AutoLock l(&mutex);
if (atype == "energy")
{
UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
IPDFEnergyExist = false;
Emin = 0.;
Emax = 1e30;
}
else if (atype == "arb")
{
ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
IPDFArbExist = false;
}
else if (atype == "epn")
{
UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
IPDFEnergyExist = false;
EpnEnergyH = ZeroPhysVector;
}
else
{
G4cout << "Error, histtype not accepted " << G4endl;
}
}
G4double G4SPSEneDistribution::GenerateOne(G4ParticleDefinition* a)
{
//Copy global shared status to thread-local one
threadLocal_t& params = threadLocalData.Get();
params.particle_definition=a;
params.particle_energy=-1;
params.Emax = Emax;
params.Emin = Emin;
params.alpha = alpha;
params.Ezero = Ezero;
params.grad = grad;
params.cept = cept;
params.weight = weight;
//particle_energy = -1.;
while ((EnergyDisType == "Arb") ? (params.particle_energy < ArbEmin || params.particle_energy > ArbEmax) : (params.particle_energy < params.Emin|| params.particle_energy > params.Emax))
{
if (Biased)
{
GenerateBiasPowEnergies();
}
else
{
if (EnergyDisType == "Mono")
{
GenerateMonoEnergetic();
}
else if (EnergyDisType == "Lin")
{
GenerateLinearEnergies(false);
}
else if (EnergyDisType == "Pow")
{
GeneratePowEnergies(false);
}
else if (EnergyDisType == "Exp")
{
GenerateExpEnergies(false);
}
else if (EnergyDisType == "Gauss")
{
GenerateGaussEnergies();
}
else if (EnergyDisType == "Brem")
{
GenerateBremEnergies();
}
else if (EnergyDisType == "Bbody")
{
GenerateBbodyEnergies();
}
else if (EnergyDisType == "Cdg")
{
GenerateCdgEnergies();
}
else if (EnergyDisType == "User")
{
GenUserHistEnergies();
}
else if (EnergyDisType == "Arb")
{
GenArbPointEnergies();
}
else if (EnergyDisType == "Epn")
{
GenEpnHistEnergies();
}
else
{
G4cout << "Error: EnergyDisType has unusual value" << G4endl;
}
}
}
return params.particle_energy;
}
G4double G4SPSEneDistribution::GetProbability(G4double ene)
{
G4double prob = 1.;
threadLocal_t& params = threadLocalData.Get();
if (EnergyDisType == "Lin")
{
if (prob_norm == 1.)
{
prob_norm = 0.5*params.grad*params.Emax*params.Emax + params.cept*params.Emax - 0.5*params.grad*params.Emin*params.Emin - params.cept*params.Emin;
}
prob = params.cept + params.grad * ene;
prob /= prob_norm;
}
else if (EnergyDisType == "Pow")
{
if (prob_norm == 1.)
{
if (alpha != -1.)
{
G4double emina = std::pow(params.Emin, params.alpha + 1);
G4double emaxa = std::pow(params.Emax, params.alpha + 1);
prob_norm = 1./(1.+alpha) * (emaxa - emina);
}
else
{
prob_norm = std::log(params.Emax) - std::log(params.Emin) ;
}
}
prob = std::pow(ene, params.alpha)/prob_norm;
}
else if (EnergyDisType == "Exp")
{
if (prob_norm == 1.)
{
prob_norm = -params.Ezero*(std::exp(-params.Emax/params.Ezero) - std::exp(params.Emin/params.Ezero));
}
prob = std::exp(-ene / params.Ezero);
prob /= prob_norm;
}
else if (EnergyDisType == "Arb")
{
prob = ArbEnergyH.Value(ene);
// prob = ArbEInt->CubicSplineInterpolation(ene);
//G4double deltaY;
//prob = ArbEInt->PolynomInterpolation(ene, deltaY);
if (prob <= 0.)
{
//G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << " " <<deltaY<< G4endl;
G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << G4endl;
prob = 1e-30;
}
// already normalised
}
else
{
G4cout << "Error: EnergyDisType not supported" << G4endl;
}
return prob;
}
Reference in New Issue View Git Blame Copy Permalink