include/Geant4/G4SPSEneDistribution.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4SPSEneDistribution
0027 //
0028 // Class Description:
0029 //
0030 // To generate the energy of a primary vertex according to the
0031 // defined distribution. This is a shared class between threads.
0032 // Only one thread should use the set-methods here.
0033 // Note that this is exactly what is achieved using UI commands.
0034 // If you use the set methods to set defaults in your application take
0035 // care that only one thread is executing them.
0036 // In addition take care of calling these methods before the run is
0037 // started. Do not use the setters during the event loop
0038
0039 // Author:       Fan Lei, QinetiQ ltd.
0040 // Customer:     ESA/ESTEC
0041 // History:
0042 // - 05/02/2004, Fan Lei - Created.
0043 //     Based on the G4GeneralParticleSource class.
0044 // - 26/03/2014, Andrew Green.
0045 //     Modification to use STL vectors instead of C-style arrays.
0046 //     Also moved to dynamically allocated memory in the LinearInterpolation(),
0047 //     ExpInterpolation() and LogInterpolation() functions.
0048 // - 06/06/2014, Andrea Dotti.
0049 //     For thread safety: this is a shared object.
0050 //     Added mutex to control access to shared resources (data members).
0051 //     in Getters and Setters, mutex is NOT used in GenerateOne() because it
0052 //     is assumed that properties are not changed during event loop.
0053 // - 24/11/2017, Fan Lei
0054 //    Added cutoff power-law distribution option. Implementation is similar
0055 //    to that of the BlackBody one.
0056 // --------------------------------------------------------------------
0057 #ifndef G4SPSEneDistribution_hh
0058 #define G4SPSEneDistribution_hh 1
0059
0060 #include "G4PhysicsFreeVector.hh"
0061 #include "G4ParticleMomentum.hh"
0062 #include "G4ParticleDefinition.hh"
0063 #include "G4DataInterpolation.hh"
0064 #include "G4Threading.hh"
0065 #include "G4Cache.hh"
0066 #include <vector>
0067
0068 #include "G4SPSRandomGenerator.hh"
0069
0070 class G4SPSEneDistribution
0071 {
0072   public:
0073
0074     G4SPSEneDistribution();
0075       // Constructor: initializes variables
0076    ~G4SPSEneDistribution();
0077       // Destructor
0078
0079     void SetEnergyDisType(const G4String&);
0080       // Allows the user to choose the energy distribution type.
0081       // The arguments are: Mono (mono-energetic), Lin (linear),
0082       // Pow (power-law), Exp (exponential), Gauss (gaussian),
0083       // Brem (bremsstrahlung), BBody (black-body),
0084       // Cdg (cosmic diffuse gamma-ray), User (user-defined),
0085       // Arb (arbitrary point-wise), Epn (energy per nucleon)
0086
0087     const G4String& GetEnergyDisType();
0088
0089     void SetEmin(G4double);
0090       // Sets the minimum energy
0091
0092     G4double GetEmin() const;
0093     G4double GetArbEmin();
0094
0095     void SetEmax(G4double);
0096       // Sets the maximum energy
0097
0098     G4double GetEmax() const;
0099     G4double GetArbEmax();
0100
0101     void SetMonoEnergy(G4double);
0102       // Sets energy for mono-energetic distribution
0103
0104     void SetAlpha(G4double);
0105       // Sets alpha for a power-law distribution
0106
0107     void SetBiasAlpha(G4double);
0108
0109     void SetTemp(G4double);
0110       // Sets Temperature for a Brem or BBody distributions
0111
0112     void SetBeamSigmaInE(G4double);
0113
0114     void SetEzero(G4double);
0115       // Sets Ezero for an exponential distribution
0116
0117     void SetGradient(G4double);
0118       // Sets gradient for a linear distribution
0119
0120     void SetInterCept(G4double);
0121       // Sets intercept for a linear distribution
0122
0123     void UserEnergyHisto(const G4ThreeVector&);
0124       // Allows user to defined a histogram for the energy distribution
0125
0126     void ArbEnergyHisto(const G4ThreeVector&);
0127       // Allows the user to define an Arbitrary set of points for the
0128       // energy distribution
0129
0130     void ArbEnergyHistoFile(const G4String&);
0131
0132     void EpnEnergyHisto(const G4ThreeVector&);
0133       // Allows the user to define an Energy per nucleon histogram
0134
0135     void InputEnergySpectra(G4bool);
0136       // Allows the user to choose between momentum and energy histograms
0137       // for user-defined histograms and arbitrary point-wise spectra.
0138       // The default is true (energy)
0139
0140     void InputDifferentialSpectra(G4bool);
0141       // Allows the user to choose between integral and differential
0142       // distributions when using the arbitrary point-wise option
0143
0144     void ArbInterpolate(const G4String&);
0145       // Allows the user to specify the type of function to
0146       // interpolate the Arbitrary points spectrum with
0147
0148     const G4String& GetIntType();
0149
0150     void Calculate();
0151       // Controls the calculation of Integral PDF for the Cdg and BBody
0152       // distributions
0153
0154     void SetBiasRndm(G4SPSRandomGenerator* a);
0155       // Sets the biased random number generator
0156
0157     void ReSetHist(const G4String&);
0158       // Resets the histogram for user defined distribution
0159
0160     void SetVerbosity(G4int a);
0161       // Sets the verbosity level
0162
0163     G4double GetWeight() const;
0164
0165     G4double GetMonoEnergy();
0166       // Mono-energetic energy
0167
0168     G4double GetSE();
0169       // Standard deviation for Gaussian distribution in energy
0170
0171     G4double Getalpha() const;
0172       // Alpha (pow)
0173
0174     G4double GetEzero() const;
0175       // E0 (exp)
0176
0177     G4double GetTemp();
0178       // Temp (bbody,brem)
0179
0180     G4double Getgrad() const;
0181       // Gradient and intercept for linear spectra
0182
0183     G4double Getcept() const;
0184
0185     G4PhysicsFreeVector GetUserDefinedEnergyHisto();
0186
0187     G4PhysicsFreeVector GetArbEnergyHisto();
0188
0189     G4double GenerateOne(G4ParticleDefinition*);
0190        // Generate one random energy for the specified particle
0191
0192     G4double GetProbability (G4double);
0193
0194     G4double GetArbEneWeight(G4double);
0195
0196     inline void ApplyEnergyWeight(G4bool val) { applyEvergyWeight = val; }
0197     inline G4bool IfApplyEnergyWeight() const { return applyEvergyWeight; }
0198
0199   private:
0200
0201     void LinearInterpolation();
0202     void LogInterpolation();
0203     void ExpInterpolation();
0204     void SplineInterpolation();
0205     void CalculateCdgSpectrum();
0206     void CalculateBbodySpectrum();
0207     void CalculateCPowSpectrum();
0208
0209     // The following methods generate energies according
0210     // to the spectral parameters defined above
0211
0212     void GenerateMonoEnergetic();
0213     void GenerateBiasPowEnergies();
0214     void GenerateGaussEnergies();
0215     void GenerateBremEnergies();
0216     void GenerateBbodyEnergies();
0217     void GenerateCdgEnergies();
0218     void GenUserHistEnergies();
0219     void GenEpnHistEnergies();
0220     void GenArbPointEnergies(); // NOTE: REQUIRES UPDATE OF DATA MEMBERS
0221     void GenerateExpEnergies(G4bool);
0222     void GenerateLinearEnergies(G4bool);
0223     void GeneratePowEnergies(G4bool);
0224     void GenerateCPowEnergies();
0225
0226     void ConvertEPNToEnergy();
0227       // Converts energy per nucleon to energy
0228
0229     void BBInitHists();
0230     void CPInitHists();
0231
0232   private:  // Non invariant data members become G4Cache
0233
0234     G4String EnergyDisType; // energy dis type Variable  - Mono,Lin,Exp,etc
0235     G4double weight; // particle weight //// NOT INVARIANT
0236     G4double MonoEnergy; //Mono-energteic energy
0237     G4double SE; // Standard deviation for Gaussian distribution in energy
0238
0239     G4double Emin, Emax; // emin and emax //// NOT INVARIANT
0240     G4double alpha, Ezero;// alpha (pow), E0 (exp) //// NOT INVARIANT
0241     G4double Temp; // Temp (bbody,brem)
0242     G4double biasalpha; // biased power index
0243     G4double grad, cept; // gradient and intercept for linear spectra //// NOT INVARIANT
0244     G4double prob_norm; // normalisation factor use in calculate the probability
0245     G4bool Biased = false; // biased to power-law
0246     G4bool EnergySpec = true; // energy spectra, false - momentum spectra
0247     G4bool DiffSpec = true; // differential spec, false integral spec
0248
0249     G4PhysicsFreeVector UDefEnergyH; // energy hist data
0250     G4PhysicsFreeVector IPDFEnergyH;
0251     G4bool IPDFEnergyExist = false, IPDFArbExist = false, Epnflag = false;
0252     G4PhysicsFreeVector ArbEnergyH; // Arb x,y histogram
0253     G4PhysicsFreeVector IPDFArbEnergyH; // IPDF for Arb
0254     G4PhysicsFreeVector EpnEnergyH;
0255     G4double CDGhist[3]; // cumulative histo for cdg
0256
0257     std::vector<G4double>* BBHist = nullptr;
0258     std::vector<G4double>* Bbody_x = nullptr;
0259     G4bool BBhistInit = false;
0260     G4bool BBhistCalcd = false;
0261
0262     //  For cutoff power-law
0263     //
0264     std::vector<G4double>* CPHist = nullptr;
0265     std::vector<G4double>* CP_x = nullptr;
0266     G4bool CPhistInit = false;
0267     G4bool CPhistCalcd = false;
0268
0269     G4String IntType; // Interpolation type
0270     G4double* Arb_grad = nullptr;
0271     G4double* Arb_cept = nullptr;
0272     G4bool Arb_grad_cept_flag = false;
0273     G4double* Arb_alpha = nullptr;
0274     G4double* Arb_Const = nullptr;
0275     G4bool Arb_alpha_Const_flag = false;
0276     G4double* Arb_ezero = nullptr;
0277     G4bool Arb_ezero_flag = false;
0278
0279     G4bool applyEvergyWeight = false;
0280
0281     G4double ArbEmin, ArbEmax;
0282       // Emin and Emax for the whole arb distribution used primarily for debug.
0283
0284     G4double particle_energy;
0285
0286     G4SPSRandomGenerator* eneRndm = nullptr;
0287
0288     G4int verbosityLevel;
0289
0290     G4PhysicsFreeVector ZeroPhysVector; // for re-set only
0291
0292     std::vector<G4DataInterpolation*> SplineInt;
0293       // Holds Spline stuff required for sampling
0294     G4DataInterpolation* Splinetemp = nullptr;
0295       // Holds a temp Spline used for calculating area
0296
0297     G4Mutex mutex; // protect access to shared resources
0298
0299     // Thread local data (non-invariant during event loop).
0300     // These are copied from master one at the beginning of
0301     // generation of each event
0302     //
0303     struct threadLocal_t
0304     {
0305       G4double Emin;
0306       G4double Emax;
0307       G4double alpha;
0308       G4double Ezero;
0309       G4double grad;
0310       G4double cept;
0311       G4ParticleDefinition* particle_definition;
0312       G4double weight;
0313       G4double particle_energy;
0314     };
0315     G4Cache<threadLocal_t> threadLocalData;
0316 };
0317
0318 #endif