Changeset 1347 for trunk/source/event

Timestamp:

Dec 22, 2010, 3:52:27 PM (15 years ago)

Author:

garnier

Message:

geant4 tag 9.4

Location:

trunk/source/event

Files:

• GNUmakefile (modified) (1 diff)
• History (modified) (2 diffs)
• include/G4GeneralParticleSourceMessenger.hh (modified) (2 diffs)
• include/G4SPSEneDistribution.hh (modified) (1 diff)
• include/G4SPSRandomGenerator.hh (modified) (1 diff)
• include/G4SingleParticleSource.hh (modified) (1 diff)
• src/G4GeneralParticleSourceMessenger.cc (modified) (77 diffs)
• src/G4SPSEneDistribution.cc (modified) (1 diff)
• src/G4SPSRandomGenerator.cc (modified) (2 diffs)
• src/G4SingleParticleSource.cc (modified) (1 diff)

trunk/source/event/GNUmakefile

@@ -1 @@
-# $Id: GNUmakefile,v 1.8 2005/10/17 21:21:53 tinslay Exp $
+# $Id: GNUmakefile,v 1.10 2010/10/27 07:21:13 gcosmo Exp $
 # ------------------------------------------------------------
 # GNUmakefile for events library.  Makoto Asai, 5/9/95.

trunk/source/event/History

@@ -1 @@
-$Id: History,v 1.131 2010/06/12 04:07:45 asaim Exp $
+$Id: History,v 1.142 2010/12/15 22:15:07 asaim Exp $
 -------------------------------------------------------------------
 
@@ -17 +17 @@
      * Reverse chronological order (last date on top), please *
      ----------------------------------------------------------
+
+December 15th, 2010, M.Asai (event-V09-03-10)
+- Fixing incorrect destination stack for the case more than three track stacks are used.
+
+November 24th, 2010, M.Asai (event-V09-03-09)
+- Added protection against null pointer : G4TrackStack.cc, G4SmartTrackStack.cc
+
+Novemver 10th, 2010, F.Lei (event-V09-03-08)
+- User can now specify arbitrary energy distribution by supplying the data points
+  in a 2-column (energy<in MeV> and flux) ASCII file, via the UI command
+   /gps/hist/file your_file_name
+- Further update of the auto energy biasing implementation
+- Bug fix (new implementation) in the Spline interpolation implementation
+- Bug fix in the Logarithm interpolation implementation (alpha = -1 case).   
+
+November 8th, 2010, M.Asai (event-V09-03-07)
+- Setting touchable history of the origin to G4Track.
+
+November 2nd, 2010, F.Lei 
+- Implemented concrete destructor for G4SPSEneDistribution (bug #1149)
+- Added automatic energy biasing scheme. Users now can bias the energy distribution
+  sampling to a given power-law distribution (index alpha) via the UI command  
+    /gps/ene/biasAlpha alpha
+  Files affected are
+        G4SingleParticleSource.hh .cc
+        G4SPSEneDistribution.hh .cc  
+        G4SPSRandomGenerator.hh .cc  
+        G4GeneralParticleSourceMessenger.hh .cc
+
+October 27th, 2010, G.Cosmo (event-V09-03-06)
+- Restored DLL setup as originally. Withdrawn changes in last tag.
+
+October 19th, 2010, G.Cosmo (event-V09-03-05)
+- Replaced G4EVENT_ALLOC_EXPORT flag with G4ALLOC_EXPORT for DLL exported
+  symbols.
+
+August 9th, 2010, M.Asai (event-V09-03-04)
+- Fixing electron mass correction for ions in G4PrimaryTransformer.
 
 June 11th, 2010, M.Asai (event-V09-03-03)

trunk/source/event/include/G4GeneralParticleSourceMessenger.hh

@@ -212 +212 @@
   G4UIcmdWithADouble         *gradientCmd1;
   G4UIcmdWithADouble         *interceptCmd1;
+  G4UIcmdWithADouble         *arbeintCmd1;
   G4UIcmdWithoutParameter    *calculateCmd1;
   G4UIcmdWithABool           *energyspecCmd1;
@@ -237 +238 @@
   // old ones, will be removed soon
   G4UIcmdWith3Vector         *histpointCmd1;
+  G4UIcmdWithAString         *histfileCmd1;
   G4UIcmdWithAString         *histnameCmd1;
   G4UIcmdWithAString         *arbintCmd1;

trunk/source/event/include/G4SPSEneDistribution.hh

@@ -144 +144 @@
 #include "G4SPSRandomGenerator.hh"
 
-class G4SPSEneDistribution
-{
+class G4SPSEneDistribution {
 public:
-  G4SPSEneDistribution (); 
-  ~G4SPSEneDistribution ();
-
-  void SetEnergyDisType(G4String);
-  inline G4String GetEnergyDisType() {return EnergyDisType;};
-  void SetEmin(G4double);
-  inline G4double GetEmin() {return Emin;} ;
-  inline G4double GetArbEmin() {return ArbEmin;} ;
-  void SetEmax(G4double);
-  inline G4double GetEmax() {return Emax;} ;
-  inline G4double GetArbEmax() {return ArbEmax;};
-  void SetMonoEnergy(G4double);
-  void SetAlpha(G4double);
-  void SetTemp(G4double);
-  void SetBeamSigmaInE(G4double);
-  void SetEzero(G4double);
-  void SetGradient(G4double);
-  void SetInterCept(G4double);
-  void UserEnergyHisto(G4ThreeVector);
-  void ArbEnergyHisto(G4ThreeVector);
-  void EpnEnergyHisto(G4ThreeVector);
-
-  void InputEnergySpectra(G4bool);
-  void InputDifferentialSpectra(G4bool);
-  void ArbInterpolate(G4String);
-  inline G4String GetIntType() {return IntType;};
-  void Calculate();
-  //
-  void SetBiasRndm(G4SPSRandomGenerator* a) {eneRndm = a; };
-  // method to re-set the histograms
-  void ReSetHist(G4String);
-  // Set the verbosity level.
-  void SetVerbosity(G4int a) {verbosityLevel = a; } ;
-  //x
-
-  G4double GetMonoEnergy () {return MonoEnergy; }; //Mono-energteic energy
-  G4double GetSE() {return SE;}; // Standard deviation for Gaussion distrbution in energy
-  G4double Getalpha () {return alpha;}; // alpha (pow)
-  G4double GetEzero() {return Ezero;}; // E0 (exp)
-  G4double GetTemp() {return Temp;}; // Temp (bbody,brem)
-  G4double Getgrad() {return grad;}; // gradient and intercept for linear spectra
-  G4double Getcept() {return cept;};
-
-  inline G4PhysicsOrderedFreeVector  GetUserDefinedEnergyHisto(){return UDefEnergyH;};
-  inline G4PhysicsOrderedFreeVector  GetArbEnergyHisto(){return ArbEnergyH;};
-
-  G4double GenerateOne(G4ParticleDefinition*);
-  
+        G4SPSEneDistribution();
+        ~G4SPSEneDistribution();
+
+        void SetEnergyDisType(G4String);
+        inline G4String GetEnergyDisType() {
+                return EnergyDisType;
+        }
+        ;
+        void SetEmin(G4double);
+        inline G4double GetEmin() {
+                return Emin;
+        }
+        ;
+        inline G4double GetArbEmin() {
+                return ArbEmin;
+        }
+        ;
+        void SetEmax(G4double);
+        inline G4double GetEmax() {
+                return Emax;
+        }
+        ;
+        inline G4double GetArbEmax() {
+                return ArbEmax;
+        }
+        ;
+        void SetMonoEnergy(G4double);
+        void SetAlpha(G4double);
+        void SetBiasAlpha(G4double);
+        void SetTemp(G4double);
+        void SetBeamSigmaInE(G4double);
+        void SetEzero(G4double);
+        void SetGradient(G4double);
+        void SetInterCept(G4double);
+        void UserEnergyHisto(G4ThreeVector);
+        void ArbEnergyHisto(G4ThreeVector);
+        void ArbEnergyHistoFile(G4String);
+        void EpnEnergyHisto(G4ThreeVector);
+
+        void InputEnergySpectra(G4bool);
+        void InputDifferentialSpectra(G4bool);
+        void ArbInterpolate(G4String);
+        inline G4String GetIntType() {
+                return IntType;
+        }
+        ;
+        void Calculate();
+        //
+        void SetBiasRndm(G4SPSRandomGenerator* a) {
+                eneRndm = a;
+        }
+        ;
+        // method to re-set the histograms
+        void ReSetHist(G4String);
+        // Set the verbosity level.
+        void SetVerbosity(G4int a) {
+                verbosityLevel = a;
+        }
+        ;
+        //x
+        G4double GetWeight() {
+                return weight;
+        }
+
+        G4double GetMonoEnergy() {
+                return MonoEnergy;
+        }
+        ; //Mono-energteic energy
+        G4double GetSE() {
+                return SE;
+        }
+        ; // Standard deviation for Gaussion distrbution in energy
+        G4double Getalpha() {
+                return alpha;
+        }
+        ; // alpha (pow)
+        G4double GetEzero() {
+                return Ezero;
+        }
+        ; // E0 (exp)
+        G4double GetTemp() {
+                return Temp;
+        }
+        ; // Temp (bbody,brem)
+        G4double Getgrad() {
+                return grad;
+        }
+        ; // gradient and intercept for linear spectra
+        G4double Getcept() {
+                return cept;
+        }
+        ;
+
+        inline G4PhysicsOrderedFreeVector GetUserDefinedEnergyHisto() {
+                return UDefEnergyH;
+        }
+        ;
+        inline G4PhysicsOrderedFreeVector GetArbEnergyHisto() {
+                return ArbEnergyH;
+        }
+        ;
+
+        G4double GenerateOne(G4ParticleDefinition*);
+        G4double GetProbability (G4double);
+
+
 private:
-  void LinearInterpolation();
-  void LogInterpolation();
-  void ExpInterpolation();
-  void SplineInterpolation();
-  void CalculateCdgSpectrum(); 
-  void CalculateBbodySpectrum();
-
-  // The following methods generate energies according to the spectral
-  // parameters defined above.
-  void GenerateMonoEnergetic();
-  void GenerateLinearEnergies(G4bool);
-  void GeneratePowEnergies(G4bool);
-  void GenerateExpEnergies(G4bool );
-  void GenerateGaussEnergies();
-  void GenerateBremEnergies();
-  void GenerateBbodyEnergies();
-  void GenerateCdgEnergies();
-  void GenUserHistEnergies();
-  void GenEpnHistEnergies();
-  void GenArbPointEnergies();
-  // converts energy per nucleon to energy.
-  void ConvertEPNToEnergy();
+        void LinearInterpolation();
+        void LogInterpolation();
+        void ExpInterpolation();
+        void SplineInterpolation();
+        void CalculateCdgSpectrum();
+        void CalculateBbodySpectrum();
+
+        // The following methods generate energies according to the spectral
+        // parameters defined above.
+        void GenerateMonoEnergetic();
+        void GenerateLinearEnergies(G4bool);
+        void GeneratePowEnergies(G4bool);
+        void GenerateBiasPowEnergies();
+        void GenerateExpEnergies(G4bool);
+        void GenerateGaussEnergies();
+        void GenerateBremEnergies();
+        void GenerateBbodyEnergies();
+        void GenerateCdgEnergies();
+        void GenUserHistEnergies();
+        void GenEpnHistEnergies();
+        void GenArbPointEnergies();
+        // converts energy per nucleon to energy.
+        void ConvertEPNToEnergy();
+
 
 private:
 
-  G4String EnergyDisType; // energy dis type Variable  - Mono,Lin,Exp,etc
-  G4double MonoEnergy; //Mono-energteic energy
-  G4double SE ; // Standard deviation for Gaussion distrbution in energy
-  G4double Emin, Emax; // emin and emax
-  G4double alpha, Ezero, Temp; // alpha (pow), E0 (exp) and Temp (bbody,brem)
-  G4double grad, cept; // gradient and intercept for linear spectra
-  G4bool EnergySpec; // true - energy spectra, false - momentum spectra
-  G4bool DiffSpec;  // true - differential spec, false integral spec
-  G4bool ApplyRig; // false no rigidity cutoff, true then apply one
-  G4double ERig; // energy of rigidity cutoff
-  G4PhysicsOrderedFreeVector UDefEnergyH; // energy hist data
-  G4PhysicsOrderedFreeVector IPDFEnergyH;  
-  G4bool IPDFEnergyExist, IPDFArbExist, Epnflag;
-  G4PhysicsOrderedFreeVector ArbEnergyH; // Arb x,y histogram
-  G4PhysicsOrderedFreeVector IPDFArbEnergyH; // IPDF for Arb
-  G4PhysicsOrderedFreeVector EpnEnergyH; 
-  G4double CDGhist[3]; // cumulative histo for cdg
-  G4double BBHist[10001], Bbody_x[10001];
-  G4String IntType; // Interpolation type
-  G4double Arb_grad[1024], Arb_cept[1024]; // grad and cept for 1024 segments
-  G4double Arb_alpha[1024], Arb_Const[1024]; // alpha and constants
-  G4double Arb_ezero[1024]; // ezero
-  G4double ArbEmin, ArbEmax; // Emin and Emax for the whole arb distribution used primarily for debug.
-
-  G4double               particle_energy;
-  G4ParticleDefinition*  particle_definition;
-
-  G4SPSRandomGenerator* eneRndm;
-
-  // Verbosity
-  G4int verbosityLevel;
-
-  G4PhysicsOrderedFreeVector ZeroPhysVector ; // for re-set only 
-
-  G4DataInterpolation *SplineInt; // holds Spline stuff
+        G4String EnergyDisType; // energy dis type Variable  - Mono,Lin,Exp,etc
+        G4double weight; // particle weight
+        G4double MonoEnergy; //Mono-energteic energy
+        G4double SE; // Standard deviation for Gaussion distrbution in energy
+        G4double Emin, Emax; // emin and emax
+        G4double alpha, Ezero, Temp; // alpha (pow), E0 (exp) and Temp (bbody,brem)
+        G4double biasalpha; // biased power index
+        G4double grad, cept; // gradient and intercept for linear spectra
+        G4double prob_norm; // normalisation factor use in calculate the probability 
+        G4bool Biased; // true - biased to power-law
+        G4bool EnergySpec; // true - energy spectra, false - momentum spectra
+        G4bool DiffSpec; // true - differential spec, false integral spec
+        G4bool ApplyRig; // false no rigidity cutoff, true then apply one
+        G4double ERig; // energy of rigidity cutoff
+        G4PhysicsOrderedFreeVector UDefEnergyH; // energy hist data
+        G4PhysicsOrderedFreeVector IPDFEnergyH;
+        G4bool IPDFEnergyExist, IPDFArbExist, Epnflag;
+        G4PhysicsOrderedFreeVector ArbEnergyH; // Arb x,y histogram
+        G4PhysicsOrderedFreeVector IPDFArbEnergyH; // IPDF for Arb
+        G4PhysicsOrderedFreeVector EpnEnergyH;
+        G4double CDGhist[3]; // cumulative histo for cdg
+        G4double BBHist[10001], Bbody_x[10001];
+        G4String IntType; // Interpolation type
+        G4double Arb_grad[1024], Arb_cept[1024]; // grad and cept for 1024 segments
+        G4double Arb_alpha[1024], Arb_Const[1024]; // alpha and constants
+        G4double Arb_ezero[1024]; // ezero
+        G4double ArbEmin, ArbEmax; // Emin and Emax for the whole arb distribution used primarily for debug.
+
+        G4double particle_energy;
+        G4ParticleDefinition* particle_definition;
+
+        G4SPSRandomGenerator* eneRndm;
+
+        // Verbosity
+        G4int verbosityLevel;
+
+        G4PhysicsOrderedFreeVector ZeroPhysVector; // for re-set only
+
+        G4DataInterpolation *SplineInt[1024]; // holds Spline stuff required for sampling
+        G4DataInterpolation *Splinetemp; // holds a temp Spline used for calculating area
 
 };
 
-
 #endif
 
-
-
-

trunk/source/event/include/G4SPSRandomGenerator.hh

@@ -136 +136 @@
 #include "G4DataInterpolation.hh"
 
-class G4SPSRandomGenerator
-{
+class G4SPSRandomGenerator {
 public:
-  G4SPSRandomGenerator (); 
-  ~G4SPSRandomGenerator ();
-
-  //  static G4SPSRandomGenerator* getInstance ();
-
-  // Biasing Methods 
-  void SetXBias(G4ThreeVector);
-  void SetYBias(G4ThreeVector);
-  void SetZBias(G4ThreeVector);
-  void SetThetaBias(G4ThreeVector);
-  void SetPhiBias(G4ThreeVector);
-  void SetEnergyBias(G4ThreeVector);
-  void SetPosThetaBias(G4ThreeVector);
-  void SetPosPhiBias(G4ThreeVector);
-  G4double GenRandX();
-  G4double GenRandY();
-  G4double GenRandZ();
-  G4double GenRandTheta();
-  G4double GenRandPhi();
-  G4double GenRandEnergy();
-  G4double GenRandPosTheta();
-  G4double GenRandPosPhi();
-
-  inline void SetIntensityWeight(G4double weight) {bweights[8]=weight;};
-
-  inline G4double GetBiasWeight() 
-  { return bweights[0]*bweights[1]*bweights[2]*bweights[3]*bweights[4]*bweights[5]*bweights[6]*bweights[7]*bweights[8];};
- 
- // method to re-set the histograms
-  void ReSetHist(G4String);
-
-  // Set the verbosity level.
-   void SetVerbosity(G4int a) {verbosityLevel = a; } ;
+        G4SPSRandomGenerator();
+        ~G4SPSRandomGenerator();
+
+        //  static G4SPSRandomGenerator* getInstance ();
+
+        // Biasing Methods
+        void SetXBias(G4ThreeVector);
+        void SetYBias(G4ThreeVector);
+        void SetZBias(G4ThreeVector);
+        void SetThetaBias(G4ThreeVector);
+        void SetPhiBias(G4ThreeVector);
+        void SetEnergyBias(G4ThreeVector);
+        void SetPosThetaBias(G4ThreeVector);
+        void SetPosPhiBias(G4ThreeVector);
+        G4double GenRandX();
+        G4double GenRandY();
+        G4double GenRandZ();
+        G4double GenRandTheta();
+        G4double GenRandPhi();
+        G4double GenRandEnergy();
+        G4double GenRandPosTheta();
+        G4double GenRandPosPhi();
+
+        inline void SetIntensityWeight(G4double weight) {
+                bweights[8] = weight;
+        }
+        ;
+
+        inline G4double GetBiasWeight() {
+                return bweights[0] * bweights[1] * bweights[2] * bweights[3]
+                                * bweights[4] * bweights[5] * bweights[6] * bweights[7]
+                                * bweights[8];
+        }
+        ;
+
+        // method to re-set the histograms
+        void ReSetHist(G4String);
+
+        // Set the verbosity level.
+        void SetVerbosity(G4int a) {
+                verbosityLevel = a;
+        }
+        ;
 
 private:
 
-  //  static G4SPSRandomGenerator  *instance;
-
-  G4bool XBias, IPDFXBias;
-  G4PhysicsOrderedFreeVector XBiasH;
-  G4PhysicsOrderedFreeVector IPDFXBiasH;
-  G4bool YBias, IPDFYBias;
-  G4PhysicsOrderedFreeVector YBiasH;
-  G4PhysicsOrderedFreeVector IPDFYBiasH;
-  G4bool ZBias, IPDFZBias;
-  G4PhysicsOrderedFreeVector ZBiasH;
-  G4PhysicsOrderedFreeVector IPDFZBiasH;
-  G4bool ThetaBias, IPDFThetaBias;
-  G4PhysicsOrderedFreeVector ThetaBiasH;
-  G4PhysicsOrderedFreeVector IPDFThetaBiasH;
-  G4bool PhiBias, IPDFPhiBias;
-  G4PhysicsOrderedFreeVector PhiBiasH;
-  G4PhysicsOrderedFreeVector IPDFPhiBiasH;
-  G4bool EnergyBias, IPDFEnergyBias;
-  G4PhysicsOrderedFreeVector EnergyBiasH;
-  G4PhysicsOrderedFreeVector IPDFEnergyBiasH;
-  G4bool PosThetaBias, IPDFPosThetaBias;
-  G4PhysicsOrderedFreeVector PosThetaBiasH;
-  G4PhysicsOrderedFreeVector IPDFPosThetaBiasH;
-  G4bool PosPhiBias, IPDFPosPhiBias;
-  G4PhysicsOrderedFreeVector PosPhiBiasH;
-  G4PhysicsOrderedFreeVector IPDFPosPhiBiasH;
-
-  G4double bweights[9]; //record x,y,z,theta,phi,energy,posThet,posPhi,intensity weights
-
-  // Verbosity
-  G4int verbosityLevel;
-
-  G4PhysicsOrderedFreeVector ZeroPhysVector ; // for re-set only 
-  
+        //  static G4SPSRandomGenerator  *instance;
+
+        G4bool XBias, IPDFXBias;
+        G4PhysicsOrderedFreeVector XBiasH;
+        G4PhysicsOrderedFreeVector IPDFXBiasH;
+        G4bool YBias, IPDFYBias;
+        G4PhysicsOrderedFreeVector YBiasH;
+        G4PhysicsOrderedFreeVector IPDFYBiasH;
+        G4bool ZBias, IPDFZBias;
+        G4PhysicsOrderedFreeVector ZBiasH;
+        G4PhysicsOrderedFreeVector IPDFZBiasH;
+        G4bool ThetaBias, IPDFThetaBias;
+        G4PhysicsOrderedFreeVector ThetaBiasH;
+        G4PhysicsOrderedFreeVector IPDFThetaBiasH;
+        G4bool PhiBias, IPDFPhiBias;
+        G4PhysicsOrderedFreeVector PhiBiasH;
+        G4PhysicsOrderedFreeVector IPDFPhiBiasH;
+        G4bool EnergyBias, IPDFEnergyBias;
+        G4PhysicsOrderedFreeVector EnergyBiasH;
+        G4PhysicsOrderedFreeVector IPDFEnergyBiasH;
+        G4bool PosThetaBias, IPDFPosThetaBias;
+        G4PhysicsOrderedFreeVector PosThetaBiasH;
+        G4PhysicsOrderedFreeVector IPDFPosThetaBiasH;
+        G4bool PosPhiBias, IPDFPosPhiBias;
+        G4PhysicsOrderedFreeVector PosPhiBiasH;
+        G4PhysicsOrderedFreeVector IPDFPosPhiBiasH;
+
+        G4double alpha;   // for biasing energy
+
+        G4double bweights[9]; //record x,y,z,theta,phi,energy,posThet,posPhi,intensity weights
+
+        // Verbosity
+        G4int verbosityLevel;
+
+        G4PhysicsOrderedFreeVector ZeroPhysVector; // for re-set only
+
 };
 
-
 #endif
 
-
-
-

trunk/source/event/include/G4SingleParticleSource.hh

@@ -121 +121 @@
 #include "G4SPSRandomGenerator.hh"
 
-class G4SingleParticleSource : public G4VPrimaryGenerator
-{
+class G4SingleParticleSource: public G4VPrimaryGenerator {
 public:
-  G4SingleParticleSource (); 
-  ~G4SingleParticleSource ();
-  void GeneratePrimaryVertex(G4Event *evt);
-  //
-  G4SPSPosDistribution* GetPosDist() {return posGenerator;};
-  G4SPSAngDistribution* GetAngDist() {return angGenerator;};
-  G4SPSEneDistribution* GetEneDist() {return eneGenerator;};
-  G4SPSRandomGenerator* GetBiasRndm() {return biasRndm;};
-
-  // Set the verbosity level.
-  void SetVerbosity(G4int);
-
-  // Set the particle species
-  void SetParticleDefinition (G4ParticleDefinition * aParticleDefinition);
-  inline G4ParticleDefinition * GetParticleDefinition () { return particle_definition;} ;
-
-  inline void SetParticleCharge(G4double aCharge) { particle_charge = aCharge; } ;
-
-  // Set polarization
-  inline void SetParticlePolarization (G4ThreeVector aVal) {particle_polarization = aVal;};
-  inline G4ThreeVector GetParticlePolarization ()  {return particle_polarization;};
-
-  // Set Time.
-  inline void SetParticleTime(G4double aTime)  { particle_time = aTime; };
-  inline G4double GetParticleTime()  { return particle_time; };
-
-  inline void SetNumberOfParticles(G4int i)  { NumberOfParticlesToBeGenerated = i; };
-  //
-  inline G4int GetNumberOfParticles() { return NumberOfParticlesToBeGenerated; };
-  inline G4ThreeVector GetParticlePosition()  { return particle_position;};
-  inline G4ThreeVector GetParticleMomentumDirection()  { return particle_momentum_direction;};
-  inline G4double GetParticleEnergy()  {return particle_energy;};
+        G4SingleParticleSource();
+        ~G4SingleParticleSource();
+        void GeneratePrimaryVertex(G4Event *evt);
+        //
+        G4SPSPosDistribution* GetPosDist() {
+                return posGenerator;
+        }
+        ;
+        G4SPSAngDistribution* GetAngDist() {
+                return angGenerator;
+        }
+        ;
+        G4SPSEneDistribution* GetEneDist() {
+                return eneGenerator;
+        }
+        ;
+        G4SPSRandomGenerator* GetBiasRndm() {
+                return biasRndm;
+        }
+        ;
+
+        // Set the verbosity level.
+        void SetVerbosity(G4int);
+
+        // Set the particle species
+        void SetParticleDefinition(G4ParticleDefinition * aParticleDefinition);
+        inline G4ParticleDefinition * GetParticleDefinition() {
+                return particle_definition;
+        }
+        ;
+
+        inline void SetParticleCharge(G4double aCharge) {
+                particle_charge = aCharge;
+        }
+        ;
+
+        // Set polarization
+        inline void SetParticlePolarization(G4ThreeVector aVal) {
+                particle_polarization = aVal;
+        }
+        ;
+        inline G4ThreeVector GetParticlePolarization() {
+                return particle_polarization;
+        }
+        ;
+
+        // Set Time.
+        inline void SetParticleTime(G4double aTime) {
+                particle_time = aTime;
+        }
+        ;
+        inline G4double GetParticleTime() {
+                return particle_time;
+        }
+        ;
+
+        inline void SetNumberOfParticles(G4int i) {
+                NumberOfParticlesToBeGenerated = i;
+        }
+        ;
+        //
+        inline G4int GetNumberOfParticles() {
+                return NumberOfParticlesToBeGenerated;
+        }
+        ;
+        inline G4ThreeVector GetParticlePosition() {
+                return particle_position;
+        }
+        ;
+        inline G4ThreeVector GetParticleMomentumDirection() {
+                return particle_momentum_direction;
+        }
+        ;
+        inline G4double GetParticleEnergy() {
+                return particle_energy;
+        }
+        ;
 
 private:
 
-  G4SPSPosDistribution* posGenerator;
-  G4SPSAngDistribution* angGenerator;
-  G4SPSEneDistribution* eneGenerator;
-  G4SPSRandomGenerator* biasRndm;
-  //
-  // Other particle properties
-  G4int                  NumberOfParticlesToBeGenerated;
-  G4ParticleDefinition * particle_definition;
-  G4ParticleMomentum     particle_momentum_direction;
-  G4double               particle_energy;
-  G4double               particle_charge;
-  G4ThreeVector          particle_position;
-  G4double               particle_time;
-  G4ThreeVector          particle_polarization;
-  G4double               particle_weight;
-
-  // Verbosity
-  G4int verbosityLevel;
+        G4SPSPosDistribution* posGenerator;
+        G4SPSAngDistribution* angGenerator;
+        G4SPSEneDistribution* eneGenerator;
+        G4SPSRandomGenerator* biasRndm;
+        //
+        // Other particle properties
+        G4int NumberOfParticlesToBeGenerated;
+        G4ParticleDefinition * particle_definition;
+        G4ParticleMomentum particle_momentum_direction;
+        G4double particle_energy;
+        G4double particle_charge;
+        G4ThreeVector particle_position;
+        G4double particle_time;
+        G4ThreeVector particle_polarization;
+        G4double particle_weight;
+
+        // Verbosity
+        G4int verbosityLevel;
 
 };
 
-
 #endif
 
-
-
-

trunk/source/event/src/G4GeneralParticleSourceMessenger.cc

@@ -151 +151 @@
   directionCmd->SetGuidance("Set momentum direction.");
   directionCmd->SetGuidance("Direction needs not to be a unit vector.");
-  directionCmd->SetParameterName("Px","Py","Pz",true,true); 
+  directionCmd->SetParameterName("Px","Py","Pz",false,false); 
   directionCmd->SetRange("Px != 0 || Py != 0 || Pz != 0");
   
   energyCmd = new G4UIcmdWithADoubleAndUnit("/gps/energy",this);
   energyCmd->SetGuidance("Set kinetic energy.");
-  energyCmd->SetParameterName("Energy",true,true);
+  energyCmd->SetParameterName("Energy",false,false);
   energyCmd->SetDefaultUnit("GeV");
   //energyCmd->SetUnitCategory("Energy");
@@ -163 +163 @@
   positionCmd = new G4UIcmdWith3VectorAndUnit("/gps/position",this);
   positionCmd->SetGuidance("Set starting position of the particle.");
-  positionCmd->SetParameterName("X","Y","Z",true,true);
+  positionCmd->SetParameterName("X","Y","Z",false,false);
   positionCmd->SetDefaultUnit("cm");
   //positionCmd->SetUnitCategory("Length");
@@ -193 +193 @@
   timeCmd = new G4UIcmdWithADoubleAndUnit("/gps/time",this);
   timeCmd->SetGuidance("Set initial time of the particle.");
-  timeCmd->SetParameterName("t0",true,true);
+  timeCmd->SetParameterName("t0",false,false);
   timeCmd->SetDefaultUnit("ns");
   //timeCmd->SetUnitCategory("Time");
@@ -200 +200 @@
   polCmd = new G4UIcmdWith3Vector("/gps/polarization",this);
   polCmd->SetGuidance("Set polarization.");
-  polCmd->SetParameterName("Px","Py","Pz",true,true); 
+  polCmd->SetParameterName("Px","Py","Pz",false,false); 
   polCmd->SetRange("Px>=-1.&&Px<=1.&&Py>=-1.&&Py<=1.&&Pz>=-1.&&Pz<=1.");
 
   numberCmd = new G4UIcmdWithAnInteger("/gps/number",this);
   numberCmd->SetGuidance("Set number of particles to be generated per vertex.");
-  numberCmd->SetParameterName("N",true,true);
+  numberCmd->SetParameterName("N",false,false);
   numberCmd->SetRange("N>0");
 
@@ -225 +225 @@
   typeCmd1->SetGuidance("Sets source distribution type.");
   typeCmd1->SetGuidance("Either Point, Beam, Plane, Surface or Volume");
-  typeCmd1->SetParameterName("DisType",true,true);
+  typeCmd1->SetParameterName("DisType",false,false);
   typeCmd1->SetDefaultValue("Point");
   typeCmd1->SetCandidates("Point Beam Plane Surface Volume");
@@ -231 +231 @@
   shapeCmd1 = new G4UIcmdWithAString("/gps/pos/shape",this);
   shapeCmd1->SetGuidance("Sets source shape for Plan, Surface or Volume type source.");
-  shapeCmd1->SetParameterName("Shape",true,true);
+  shapeCmd1->SetParameterName("Shape",false,false);
   shapeCmd1->SetDefaultValue("NULL");
   shapeCmd1->SetCandidates("Circle Annulus Ellipse Square Rectangle Sphere Ellipsoid Cylinder Para");
@@ -238 +238 @@
   centreCmd1->SetGuidance("Set centre coordinates of source.");
   centreCmd1->SetGuidance("   same effect as the /gps/position command");
-  centreCmd1->SetParameterName("X","Y","Z",true,true);
+  centreCmd1->SetParameterName("X","Y","Z",false,false);
   centreCmd1->SetDefaultUnit("cm");
   //  centreCmd1->SetUnitCandidates("micron mm cm m km");
@@ -245 +245 @@
   posrot1Cmd1->SetGuidance("Set the 1st vector defining the rotation matrix'.");
   posrot1Cmd1->SetGuidance("It does not need to be a unit vector.");
-  posrot1Cmd1->SetParameterName("R1x","R1y","R1z",true,true); 
+  posrot1Cmd1->SetParameterName("R1x","R1y","R1z",false,false); 
   posrot1Cmd1->SetRange("R1x != 0 || R1y != 0 || R1z != 0");
 
@@ -251 +251 @@
   posrot2Cmd1->SetGuidance("Set the 2nd vector defining the rotation matrix'.");
   posrot2Cmd1->SetGuidance("It does not need to be a unit vector.");
-  posrot2Cmd1->SetParameterName("R2x","R2y","R2z",true,true); 
+  posrot2Cmd1->SetParameterName("R2x","R2y","R2z",false,false); 
   posrot2Cmd1->SetRange("R2x != 0 || R2y != 0 || R2z != 0");
 
   halfxCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/pos/halfx",this);
   halfxCmd1->SetGuidance("Set x half length of source.");
-  halfxCmd1->SetParameterName("Halfx",true,true);
+  halfxCmd1->SetParameterName("Halfx",false,false);
   halfxCmd1->SetDefaultUnit("cm");
   //  halfxCmd1->SetUnitCandidates("micron mm cm m km");
@@ -262 +262 @@
   halfyCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/pos/halfy",this);
   halfyCmd1->SetGuidance("Set y half length of source.");
-  halfyCmd1->SetParameterName("Halfy",true,true);
+  halfyCmd1->SetParameterName("Halfy",false,false);
   halfyCmd1->SetDefaultUnit("cm");
   //  halfyCmd1->SetUnitCandidates("micron mm cm m km");
@@ -268 +268 @@
   halfzCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/pos/halfz",this);
   halfzCmd1->SetGuidance("Set z half length of source.");
-  halfzCmd1->SetParameterName("Halfz",true,true);
+  halfzCmd1->SetParameterName("Halfz",false,false);
   halfzCmd1->SetDefaultUnit("cm");
   //  halfzCmd1->SetUnitCandidates("micron mm cm m km");
@@ -274 +274 @@
   radiusCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/pos/radius",this);
   radiusCmd1->SetGuidance("Set radius of source.");
-  radiusCmd1->SetParameterName("Radius",true,true);
+  radiusCmd1->SetParameterName("Radius",false,false);
   radiusCmd1->SetDefaultUnit("cm");
   //  radiusCmd1->SetUnitCandidates("micron mm cm m km");
@@ -280 +280 @@
   radius0Cmd1 = new G4UIcmdWithADoubleAndUnit("/gps/pos/inner_radius",this);
   radius0Cmd1->SetGuidance("Set inner radius of source when required.");
-  radius0Cmd1->SetParameterName("Radius0",true,true);
+  radius0Cmd1->SetParameterName("Radius0",false,false);
   radius0Cmd1->SetDefaultUnit("cm");
   //  radius0Cmd1->SetUnitCandidates("micron mm cm m km");
@@ -287 +287 @@
   possigmarCmd1->SetGuidance("Set standard deviation in radial of the beam positional profile");
   possigmarCmd1->SetGuidance(" applicable to Beam type source only");
-  possigmarCmd1->SetParameterName("Sigmar",true,true);
+  possigmarCmd1->SetParameterName("Sigmar",false,false);
   possigmarCmd1->SetDefaultUnit("cm");
   //  possigmarCmd1->SetUnitCandidates("micron mm cm m km");
@@ -294 +294 @@
   possigmaxCmd1->SetGuidance("Set standard deviation of beam positional profile in x-dir");
   possigmaxCmd1->SetGuidance(" applicable to Beam type source only");
-  possigmaxCmd1->SetParameterName("Sigmax",true,true);
+  possigmaxCmd1->SetParameterName("Sigmax",false,false);
   possigmaxCmd1->SetDefaultUnit("cm");
   //  possigmaxCmd1->SetUnitCandidates("micron mm cm m km");
@@ -301 +301 @@
   possigmayCmd1->SetGuidance("Set standard deviation of beam positional profile in y-dir");
   possigmayCmd1->SetGuidance(" applicable to Beam type source only");
-  possigmayCmd1->SetParameterName("Sigmay",true,true);
+  possigmayCmd1->SetParameterName("Sigmay",false,false);
   possigmayCmd1->SetDefaultUnit("cm");
   //  possigmayCmd1->SetUnitCandidates("micron mm cm m km");
@@ -307 +307 @@
   paralpCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/pos/paralp",this);
   paralpCmd1->SetGuidance("Angle from y-axis of y' in Para");
-  paralpCmd1->SetParameterName("paralp",true,true);
+  paralpCmd1->SetParameterName("paralp",false,false);
   paralpCmd1->SetDefaultUnit("rad");
   //  paralpCmd1->SetUnitCandidates("rad deg");
@@ -313 +313 @@
   partheCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/pos/parthe",this);
   partheCmd1->SetGuidance("Polar angle through centres of z faces");
-  partheCmd1->SetParameterName("parthe",true,true);
+  partheCmd1->SetParameterName("parthe",false,false);
   partheCmd1->SetDefaultUnit("rad");
   //  partheCmd1->SetUnitCandidates("rad deg");
@@ -319 +319 @@
   parphiCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/pos/parphi",this);
   parphiCmd1->SetGuidance("Azimuth angle through centres of z faces");
-  parphiCmd1->SetParameterName("parphi",true,true);
+  parphiCmd1->SetParameterName("parphi",false,false);
   parphiCmd1->SetDefaultUnit("rad");
   //  parphiCmd1->SetUnitCandidates("rad deg");
@@ -326 +326 @@
   confineCmd1->SetGuidance("Confine source to volume (NULL to unset).");
   confineCmd1->SetGuidance("usage: confine VolName");
-  confineCmd1->SetParameterName("VolName",true,true);
+  confineCmd1->SetParameterName("VolName",false,false);
   confineCmd1->SetDefaultValue("NULL");
 
@@ -333 +333 @@
   typeCmd->SetGuidance("Sets source distribution type. (obsolete!)");
   typeCmd->SetGuidance("Either Point, Beam, Plane, Surface or Volume");
-  typeCmd->SetParameterName("DisType",true,true);
+  typeCmd->SetParameterName("DisType",false,false);
   typeCmd->SetDefaultValue("Point");
   typeCmd->SetCandidates("Point Beam Plane Surface Volume");
@@ -339 +339 @@
   shapeCmd = new G4UIcmdWithAString("/gps/shape",this);
   shapeCmd->SetGuidance("Sets source shape type.(obsolete!)");
-  shapeCmd->SetParameterName("Shape",true,true);
+  shapeCmd->SetParameterName("Shape",false,false);
   shapeCmd->SetDefaultValue("NULL");
   shapeCmd->SetCandidates("Circle Annulus Ellipse Square Rectangle Sphere Ellipsoid Cylinder Para");
@@ -345 +345 @@
   centreCmd = new G4UIcmdWith3VectorAndUnit("/gps/centre",this);
   centreCmd->SetGuidance("Set centre coordinates of source.(obsolete!)");
-  centreCmd->SetParameterName("X","Y","Z",true,true);
+  centreCmd->SetParameterName("X","Y","Z",false,false);
   centreCmd->SetDefaultUnit("cm");
   //  centreCmd->SetUnitCandidates("micron mm cm m km");
@@ -352 +352 @@
   posrot1Cmd->SetGuidance("Set rotation matrix of x'.(obsolete!)");
   posrot1Cmd->SetGuidance("Posrot1 does not need to be a unit vector.");
-  posrot1Cmd->SetParameterName("R1x","R1y","R1z",true,true); 
+  posrot1Cmd->SetParameterName("R1x","R1y","R1z",false,false); 
   posrot1Cmd->SetRange("R1x != 0 || R1y != 0 || R1z != 0");
 
@@ -358 +358 @@
   posrot2Cmd->SetGuidance("Set rotation matrix of y'.(obsolete!)");
   posrot2Cmd->SetGuidance("Posrot2 does not need to be a unit vector.");
-  posrot2Cmd->SetParameterName("R2x","R2y","R2z",true,true); 
+  posrot2Cmd->SetParameterName("R2x","R2y","R2z",false,false); 
   posrot2Cmd->SetRange("R2x != 0 || R2y != 0 || R2z != 0");
 
   halfxCmd = new G4UIcmdWithADoubleAndUnit("/gps/halfx",this);
   halfxCmd->SetGuidance("Set x half length of source.(obsolete!)");
-  halfxCmd->SetParameterName("Halfx",true,true);
+  halfxCmd->SetParameterName("Halfx",false,false);
   halfxCmd->SetDefaultUnit("cm");
   //  halfxCmd->SetUnitCandidates("micron mm cm m km");
@@ -369 +369 @@
   halfyCmd = new G4UIcmdWithADoubleAndUnit("/gps/halfy",this);
   halfyCmd->SetGuidance("Set y half length of source.(obsolete!)");
-  halfyCmd->SetParameterName("Halfy",true,true);
+  halfyCmd->SetParameterName("Halfy",false,false);
   halfyCmd->SetDefaultUnit("cm");
   //  halfyCmd->SetUnitCandidates("micron mm cm m km");
@@ -375 +375 @@
   halfzCmd = new G4UIcmdWithADoubleAndUnit("/gps/halfz",this);
   halfzCmd->SetGuidance("Set z half length of source.(obsolete!)");
-  halfzCmd->SetParameterName("Halfz",true,true);
+  halfzCmd->SetParameterName("Halfz",false,false);
   halfzCmd->SetDefaultUnit("cm");
   //  halfzCmd->SetUnitCandidates("micron mm cm m km");
@@ -381 +381 @@
   radiusCmd = new G4UIcmdWithADoubleAndUnit("/gps/radius",this);
   radiusCmd->SetGuidance("Set radius of source.(obsolete!)");
-  radiusCmd->SetParameterName("Radius",true,true);
+  radiusCmd->SetParameterName("Radius",false,false);
   radiusCmd->SetDefaultUnit("cm");
   //  radiusCmd->SetUnitCandidates("micron mm cm m km");
@@ -387 +387 @@
   radius0Cmd = new G4UIcmdWithADoubleAndUnit("/gps/radius0",this);
   radius0Cmd->SetGuidance("Set inner radius of source.(obsolete!)");
-  radius0Cmd->SetParameterName("Radius0",true,true);
+  radius0Cmd->SetParameterName("Radius0",false,false);
   radius0Cmd->SetDefaultUnit("cm");
   //  radius0Cmd->SetUnitCandidates("micron mm cm m km");
@@ -393 +393 @@
   possigmarCmd = new G4UIcmdWithADoubleAndUnit("/gps/sigmaposr",this);
   possigmarCmd->SetGuidance("Set standard deviation of beam position in radial(obsolete!)");
-  possigmarCmd->SetParameterName("Sigmar",true,true);
+  possigmarCmd->SetParameterName("Sigmar",false,false);
   possigmarCmd->SetDefaultUnit("cm");
   // possigmarCmd->SetUnitCandidates("micron mm cm m km");
@@ -399 +399 @@
   possigmaxCmd = new G4UIcmdWithADoubleAndUnit("/gps/sigmaposx",this);
   possigmaxCmd->SetGuidance("Set standard deviation of beam position in x-dir(obsolete!)");
-  possigmaxCmd->SetParameterName("Sigmax",true,true);
+  possigmaxCmd->SetParameterName("Sigmax",false,false);
   possigmaxCmd->SetDefaultUnit("cm");
   //  possigmaxCmd->SetUnitCandidates("micron mm cm m km");
@@ -405 +405 @@
   possigmayCmd = new G4UIcmdWithADoubleAndUnit("/gps/sigmaposy",this);
   possigmayCmd->SetGuidance("Set standard deviation of beam position in y-dir(obsolete!)");
-  possigmayCmd->SetParameterName("Sigmay",true,true);
+  possigmayCmd->SetParameterName("Sigmay",false,false);
   possigmayCmd->SetDefaultUnit("cm");
   //  possigmayCmd->SetUnitCandidates("micron mm cm m km");
@@ -411 +411 @@
   paralpCmd = new G4UIcmdWithADoubleAndUnit("/gps/paralp",this);
   paralpCmd->SetGuidance("Angle from y-axis of y' in Para(obsolete!)");
-  paralpCmd->SetParameterName("paralp",true,true);
+  paralpCmd->SetParameterName("paralp",false,false);
   paralpCmd->SetDefaultUnit("rad");
   //  paralpCmd->SetUnitCandidates("rad deg");
@@ -417 +417 @@
   partheCmd = new G4UIcmdWithADoubleAndUnit("/gps/parthe",this);
   partheCmd->SetGuidance("Polar angle through centres of z faces(obsolete!)");
-  partheCmd->SetParameterName("parthe",true,true);
+  partheCmd->SetParameterName("parthe",false,false);
   partheCmd->SetDefaultUnit("rad");
   //  partheCmd->SetUnitCandidates("rad deg");
@@ -423 +423 @@
   parphiCmd = new G4UIcmdWithADoubleAndUnit("/gps/parphi",this);
   parphiCmd->SetGuidance("Azimuth angle through centres of z faces(obsolete!)");
-  parphiCmd->SetParameterName("parphi",true,true);
+  parphiCmd->SetParameterName("parphi",false,false);
   parphiCmd->SetDefaultUnit("rad");
   //  parphiCmd->SetUnitCandidates("rad deg");
@@ -430 +430 @@
   confineCmd->SetGuidance("Confine source to volume (NULL to unset)(obsolete!) .");
   confineCmd->SetGuidance("usage: confine VolName");
-  confineCmd->SetParameterName("VolName",true,true);
+  confineCmd->SetParameterName("VolName",false,false);
   confineCmd->SetDefaultValue("NULL");
 
@@ -440 +440 @@
   angtypeCmd1->SetGuidance("Sets angular source distribution type");
   angtypeCmd1->SetGuidance("Possible variables are: iso, cos, planar, beam1d, beam2d, focused or user");
-  angtypeCmd1->SetParameterName("AngDis",true,true);
+  angtypeCmd1->SetParameterName("AngDis",false,false);
   angtypeCmd1->SetDefaultValue("iso");
   angtypeCmd1->SetCandidates("iso cos planar beam1d beam2d focused user");
@@ -447 +447 @@
   angrot1Cmd1->SetGuidance("Sets the 1st vector for angular distribution rotation matrix");
   angrot1Cmd1->SetGuidance("Need not be a unit vector");
-  angrot1Cmd1->SetParameterName("AR1x","AR1y","AR1z",true,true);
+  angrot1Cmd1->SetParameterName("AR1x","AR1y","AR1z",false,false);
   angrot1Cmd1->SetRange("AR1x != 0 || AR1y != 0 || AR1z != 0");
 
@@ -453 +453 @@
   angrot2Cmd1->SetGuidance("Sets the 2nd vector for angular distribution rotation matrix");
   angrot2Cmd1->SetGuidance("Need not be a unit vector");
-  angrot2Cmd1->SetParameterName("AR2x","AR2y","AR2z",true,true);
+  angrot2Cmd1->SetParameterName("AR2x","AR2y","AR2z",false,false);
   angrot2Cmd1->SetRange("AR2x != 0 || AR2y != 0 || AR2z != 0");
 
   minthetaCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/ang/mintheta",this);
   minthetaCmd1->SetGuidance("Set minimum theta");
-  minthetaCmd1->SetParameterName("MinTheta",true,true);
+  minthetaCmd1->SetParameterName("MinTheta",false,false);
   minthetaCmd1->SetDefaultValue(0.);
   minthetaCmd1->SetDefaultUnit("rad");
@@ -465 +465 @@
   maxthetaCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/ang/maxtheta",this);
   maxthetaCmd1->SetGuidance("Set maximum theta");
-  maxthetaCmd1->SetParameterName("MaxTheta",true,true);
+  maxthetaCmd1->SetParameterName("MaxTheta",false,false);
   maxthetaCmd1->SetDefaultValue(pi);
   maxthetaCmd1->SetDefaultUnit("rad");
@@ -472 +472 @@
   minphiCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/ang/minphi",this);
   minphiCmd1->SetGuidance("Set minimum phi");
-  minphiCmd1->SetParameterName("MinPhi",true,true);
+  minphiCmd1->SetParameterName("MinPhi",false,false);
   minphiCmd1->SetDefaultUnit("rad");
   //  minphiCmd1->SetUnitCandidates("rad deg");
@@ -478 +478 @@
   maxphiCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/ang/maxphi",this);
   maxphiCmd1->SetGuidance("Set maximum phi");
-  maxphiCmd1->SetParameterName("MaxPhi",true,true);
+  maxphiCmd1->SetParameterName("MaxPhi",false,false);
   maxphiCmd1->SetDefaultValue(pi);
   maxphiCmd1->SetDefaultUnit("rad");
@@ -485 +485 @@
   angsigmarCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/ang/sigma_r",this);
   angsigmarCmd1->SetGuidance("Set standard deviation in direction for 1D beam.");
-  angsigmarCmd1->SetParameterName("Sigmara",true,true);
+  angsigmarCmd1->SetParameterName("Sigmara",false,false);
   angsigmarCmd1->SetDefaultUnit("rad");
   //  angsigmarCmd1->SetUnitCandidates("rad deg");
@@ -491 +491 @@
   angsigmaxCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/ang/sigma_x",this);
   angsigmaxCmd1->SetGuidance("Set standard deviation in direction in x-direc. for 2D beam");
-  angsigmaxCmd1->SetParameterName("Sigmaxa",true,true);
+  angsigmaxCmd1->SetParameterName("Sigmaxa",false,false);
   angsigmaxCmd1->SetDefaultUnit("rad");
   //  angsigmaxCmd1->SetUnitCandidates("rad deg");
@@ -497 +497 @@
   angsigmayCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/ang/sigma_y",this);
   angsigmayCmd1->SetGuidance("Set standard deviation in direction in y-direc. for 2D beam");
-  angsigmayCmd1->SetParameterName("Sigmaya",true,true);
+  angsigmayCmd1->SetParameterName("Sigmaya",false,false);
   angsigmayCmd1->SetDefaultUnit("rad");
   //  angsigmayCmd1->SetUnitCandidates("rad deg");
@@ -503 +503 @@
   angfocusCmd = new G4UIcmdWith3VectorAndUnit("/gps/ang/focuspoint",this);
   angfocusCmd->SetGuidance("Set the focusing point for the beam");
-  angfocusCmd->SetParameterName("x","y","z",true,true);
+  angfocusCmd->SetParameterName("x","y","z",false,false);
   angfocusCmd->SetDefaultUnit("cm");
   //  angfocusCmd->SetUnitCandidates("micron mm cm m km");
@@ -523 +523 @@
   angtypeCmd->SetGuidance("Sets angular source distribution type (obsolete!)");
   angtypeCmd->SetGuidance("Possible variables are: iso, cos planar beam1d beam2d or user");
-  angtypeCmd->SetParameterName("AngDis",true,true);
+  angtypeCmd->SetParameterName("AngDis",false,false);
   angtypeCmd->SetDefaultValue("iso");
   angtypeCmd->SetCandidates("iso cos planar beam1d beam2d user");
@@ -530 +530 @@
   angrot1Cmd->SetGuidance("Sets the x' vector for angular distribution(obsolete!) ");
   angrot1Cmd->SetGuidance("Need not be a unit vector");
-  angrot1Cmd->SetParameterName("AR1x","AR1y","AR1z",true,true);
+  angrot1Cmd->SetParameterName("AR1x","AR1y","AR1z",false,false);
   angrot1Cmd->SetRange("AR1x != 0 || AR1y != 0 || AR1z != 0");
 
@@ -536 +536 @@
   angrot2Cmd->SetGuidance("Sets the y' vector for angular distribution (obsolete!)");
   angrot2Cmd->SetGuidance("Need not be a unit vector");
-  angrot2Cmd->SetParameterName("AR2x","AR2y","AR2z",true,true);
+  angrot2Cmd->SetParameterName("AR2x","AR2y","AR2z",false,false);
   angrot2Cmd->SetRange("AR2x != 0 || AR2y != 0 || AR2z != 0");
 
   minthetaCmd = new G4UIcmdWithADoubleAndUnit("/gps/mintheta",this);
   minthetaCmd->SetGuidance("Set minimum theta (obsolete!)");
-  minthetaCmd->SetParameterName("MinTheta",true,true);
+  minthetaCmd->SetParameterName("MinTheta",false,false);
   minthetaCmd->SetDefaultUnit("rad");
   //  minthetaCmd->SetUnitCandidates("rad deg");
@@ -547 +547 @@
   maxthetaCmd = new G4UIcmdWithADoubleAndUnit("/gps/maxtheta",this);
   maxthetaCmd->SetGuidance("Set maximum theta (obsolete!)");
-  maxthetaCmd->SetParameterName("MaxTheta",true,true);
+  maxthetaCmd->SetParameterName("MaxTheta",false,false);
   maxthetaCmd->SetDefaultValue(3.1416);
   maxthetaCmd->SetDefaultUnit("rad");
@@ -554 +554 @@
   minphiCmd = new G4UIcmdWithADoubleAndUnit("/gps/minphi",this);
   minphiCmd->SetGuidance("Set minimum phi (obsolete!)");
-  minphiCmd->SetParameterName("MinPhi",true,true);
+  minphiCmd->SetParameterName("MinPhi",false,false);
   minphiCmd->SetDefaultUnit("rad");
   //  minphiCmd->SetUnitCandidates("rad deg");
@@ -560 +560 @@
   maxphiCmd = new G4UIcmdWithADoubleAndUnit("/gps/maxphi",this);
   maxphiCmd->SetGuidance("Set maximum phi(obsolete!)");
-  maxphiCmd->SetParameterName("MaxPhi",true,true);
+  maxphiCmd->SetParameterName("MaxPhi",false,false);
   maxphiCmd->SetDefaultUnit("rad");
   //  maxphiCmd->SetUnitCandidates("rad deg");
@@ -566 +566 @@
   angsigmarCmd = new G4UIcmdWithADoubleAndUnit("/gps/sigmaangr",this);
   angsigmarCmd->SetGuidance("Set standard deviation of beam direction in radial(obsolete!).");
-  angsigmarCmd->SetParameterName("Sigmara",true,true);
+  angsigmarCmd->SetParameterName("Sigmara",false,false);
   angsigmarCmd->SetDefaultUnit("rad");
   //  angsigmarCmd->SetUnitCandidates("rad deg");
@@ -572 +572 @@
   angsigmaxCmd = new G4UIcmdWithADoubleAndUnit("/gps/sigmaangx",this);
   angsigmaxCmd->SetGuidance("Set standard deviation of beam direction in x-direc(obsolete!).");
-  angsigmaxCmd->SetParameterName("Sigmaxa",true,true);
+  angsigmaxCmd->SetParameterName("Sigmaxa",false,false);
   angsigmaxCmd->SetDefaultUnit("rad");
   //  angsigmaxCmd->SetUnitCandidates("rad deg");
@@ -578 +578 @@
   angsigmayCmd = new G4UIcmdWithADoubleAndUnit("/gps/sigmaangy",this);
   angsigmayCmd->SetGuidance("Set standard deviation of beam direction in y-direc.(obsolete!)");
-  angsigmayCmd->SetParameterName("Sigmaya",true,true);
+  angsigmayCmd->SetParameterName("Sigmaya",false,false);
   angsigmayCmd->SetDefaultUnit("rad");
   //  angsigmayCmd->SetUnitCandidates("rad deg");
@@ -601 +601 @@
   energytypeCmd1 = new G4UIcmdWithAString("/gps/ene/type",this);
   energytypeCmd1->SetGuidance("Sets energy distribution type");
-  energytypeCmd1->SetParameterName("EnergyDis",true,true);
+  energytypeCmd1->SetParameterName("EnergyDis",false,false);
   energytypeCmd1->SetDefaultValue("Mono");
   energytypeCmd1->SetCandidates("Mono Lin Pow Exp Gauss Brem Bbody Cdg User Arb Epn");
@@ -607 +607 @@
   eminCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/ene/min",this);
   eminCmd1->SetGuidance("Sets minimum energy");
-  eminCmd1->SetParameterName("emin",true,true);
+  eminCmd1->SetParameterName("emin",false,false);
   eminCmd1->SetDefaultUnit("keV");
   //  eminCmd1->SetUnitCandidates("eV keV MeV GeV TeV PeV");
@@ -613 +613 @@
   emaxCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/ene/max",this);
   emaxCmd1->SetGuidance("Sets maximum energy");
-  emaxCmd1->SetParameterName("emax",true,true);
+  emaxCmd1->SetParameterName("emax",false,false);
   emaxCmd1->SetDefaultUnit("keV");
   //  emaxCmd1->SetUnitCandidates("eV keV MeV GeV TeV PeV");
@@ -619 +619 @@
   monoenergyCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/ene/mono",this);
   monoenergyCmd1->SetGuidance("Sets a monocromatic energy (same as  gps/energy)");
-  monoenergyCmd1->SetParameterName("monoenergy",true,true);
+  monoenergyCmd1->SetParameterName("monoenergy",false,false);
   monoenergyCmd1->SetDefaultUnit("keV");
   //  monoenergyCmd1->SetUnitCandidates("eV keV MeV GeV TeV PeV");
@@ -625 +625 @@
   engsigmaCmd1 = new G4UIcmdWithADoubleAndUnit("/gps/ene/sigma",this);
   engsigmaCmd1->SetGuidance("Sets the standard deviation for Gaussian energy dist.");
-  engsigmaCmd1->SetParameterName("Sigmae",true,true);
+  engsigmaCmd1->SetParameterName("Sigmae",false,false);
   engsigmaCmd1->SetDefaultUnit("keV");
   //  engsigmaCmd1->SetUnitCandidates("eV keV MeV GeV TeV PeV");
@@ -631 +631 @@
   alphaCmd1 = new G4UIcmdWithADouble("/gps/ene/alpha",this);
   alphaCmd1->SetGuidance("Sets Alpha (index) for power-law energy dist.");
-  alphaCmd1->SetParameterName("alpha",true,true);
+  alphaCmd1->SetParameterName("alpha",false,false);
   
   tempCmd1 = new G4UIcmdWithADouble("/gps/ene/temp",this);
   tempCmd1->SetGuidance("Sets the temperature for Brem and BBody distributions (in Kelvin)");
-  tempCmd1->SetParameterName("temp",true,true);
+  tempCmd1->SetParameterName("temp",false,false);
 
   ezeroCmd1 = new G4UIcmdWithADouble("/gps/ene/ezero",this);
   ezeroCmd1->SetGuidance("Sets E_0 for exponential distribution (in MeV)");
-  ezeroCmd1->SetParameterName("ezero",true,true);
+  ezeroCmd1->SetParameterName("ezero",false,false);
 
   gradientCmd1 = new G4UIcmdWithADouble("/gps/ene/gradient",this);
   gradientCmd1->SetGuidance("Sets the gradient for Lin distribution (in 1/MeV)");
-  gradientCmd1->SetParameterName("gradient",true,true);
+  gradientCmd1->SetParameterName("gradient",false,false);
 
   interceptCmd1 = new G4UIcmdWithADouble("/gps/ene/intercept",this);
   interceptCmd1->SetGuidance("Sets the intercept for Lin distributions (in MeV)");
-  interceptCmd1->SetParameterName("intercept",true,true);
+  interceptCmd1->SetParameterName("intercept",false,false);
+
+  arbeintCmd1 = new G4UIcmdWithADouble("/gps/ene/biasAlpha",this);
+  arbeintCmd1->SetGuidance("Set the power-law index for the energy sampling distri. )");
+  arbeintCmd1->SetParameterName("arbeint",false,false);
 
   calculateCmd1 = new G4UIcmdWithoutParameter("/gps/ene/calculate",this);
@@ -665 +669 @@
   energytypeCmd = new G4UIcmdWithAString("/gps/energytype",this);
   energytypeCmd->SetGuidance("Sets energy distribution type (obsolete!)");
-  energytypeCmd->SetParameterName("EnergyDis",true,true);
+  energytypeCmd->SetParameterName("EnergyDis",false,false);
   energytypeCmd->SetDefaultValue("Mono");
   energytypeCmd->SetCandidates("Mono Lin Pow Exp Gauss Brem Bbody Cdg User Arb Epn");
@@ -671 +675 @@
   eminCmd = new G4UIcmdWithADoubleAndUnit("/gps/emin",this);
   eminCmd->SetGuidance("Sets Emin (obsolete!)");
-  eminCmd->SetParameterName("emin",true,true);
+  eminCmd->SetParameterName("emin",false,false);
   eminCmd->SetDefaultUnit("keV");
   //  eminCmd->SetUnitCandidates("eV keV MeV GeV TeV PeV");
@@ -677 +681 @@
   emaxCmd = new G4UIcmdWithADoubleAndUnit("/gps/emax",this);
   emaxCmd->SetGuidance("Sets Emax (obsolete!)");
-  emaxCmd->SetParameterName("emax",true,true);
+  emaxCmd->SetParameterName("emax",false,false);
   emaxCmd->SetDefaultUnit("keV");
   //  emaxCmd->SetUnitCandidates("eV keV MeV GeV TeV PeV");
@@ -683 +687 @@
   monoenergyCmd = new G4UIcmdWithADoubleAndUnit("/gps/monoenergy",this);
   monoenergyCmd->SetGuidance("Sets Monoenergy (obsolete, use gps/energy instead!)");
-  monoenergyCmd->SetParameterName("monoenergy",true,true);
+  monoenergyCmd->SetParameterName("monoenergy",false,false);
   monoenergyCmd->SetDefaultUnit("keV");
   //  monoenergyCmd->SetUnitCandidates("eV keV MeV GeV TeV PeV");
@@ -689 +693 @@
   engsigmaCmd = new G4UIcmdWithADoubleAndUnit("/gps/sigmae",this);
   engsigmaCmd->SetGuidance("Sets the standard deviation for Gaussian energy dist.(obsolete!)");
-  engsigmaCmd->SetParameterName("Sigmae",true,true);
+  engsigmaCmd->SetParameterName("Sigmae",false,false);
   engsigmaCmd->SetDefaultUnit("keV");
   //  engsigmaCmd->SetUnitCandidates("eV keV MeV GeV TeV PeV");
@@ -695 +699 @@
   alphaCmd = new G4UIcmdWithADouble("/gps/alpha",this);
   alphaCmd->SetGuidance("Sets Alpha (index) for power-law energy dist(obsolete!).");
-  alphaCmd->SetParameterName("alpha",true,true);
+  alphaCmd->SetParameterName("alpha",false,false);
   
   tempCmd = new G4UIcmdWithADouble("/gps/temp",this);
   tempCmd->SetGuidance("Sets the temperature for Brem and BBody (in Kelvin)(obsolete!)");
-  tempCmd->SetParameterName("temp",true,true);
+  tempCmd->SetParameterName("temp",false,false);
 
   ezeroCmd = new G4UIcmdWithADouble("/gps/ezero",this);
   ezeroCmd->SetGuidance("Sets ezero exponential distributions (in MeV)(obsolete!)");
-  ezeroCmd->SetParameterName("ezero",true,true);
+  ezeroCmd->SetParameterName("ezero",false,false);
 
   gradientCmd = new G4UIcmdWithADouble("/gps/gradient",this);
   gradientCmd->SetGuidance("Sets the gradient for Lin distributions (in 1/MeV)(obsolete!)");
-  gradientCmd->SetParameterName("gradient",true,true);
+  gradientCmd->SetParameterName("gradient",false,false);
 
   interceptCmd = new G4UIcmdWithADouble("/gps/intercept",this);
   interceptCmd->SetGuidance("Sets the intercept for Lin distributions (in MeV)(obsolete!)");
-  interceptCmd->SetParameterName("intercept",true,true);
+  interceptCmd->SetParameterName("intercept",false,false);
 
   calculateCmd = new G4UIcmdWithoutParameter("/gps/calculate",this);
@@ -732 +736 @@
   histnameCmd1 = new G4UIcmdWithAString("/gps/hist/type",this);
   histnameCmd1->SetGuidance("Sets histogram type");
-  histnameCmd1->SetParameterName("HistType",true,true);
+  histnameCmd1->SetParameterName("HistType",false,false);
   histnameCmd1->SetDefaultValue("biasx");
   histnameCmd1->SetCandidates("biasx biasy biasz biast biasp biase biaspt biaspp theta phi energy arb epn");
@@ -738 +742 @@
   resethistCmd1 = new G4UIcmdWithAString("/gps/hist/reset",this);
   resethistCmd1->SetGuidance("Reset (clean) the histogram ");
-  resethistCmd1->SetParameterName("HistType",true,true);
+  resethistCmd1->SetParameterName("HistType",false,false);
   resethistCmd1->SetDefaultValue("energy");
   resethistCmd1->SetCandidates("biasx biasy biasz biast biasp biase biaspt biaspp theta phi energy arb epn");
@@ -748 +752 @@
   histpointCmd1->SetRange("Ehi >= 0. && Weight >= 0.");
 
+  histfileCmd1 = new G4UIcmdWithAString("/gps/hist/file",this);
+  histfileCmd1->SetGuidance("import the arb energy hist in an ASCII file");
+  histfileCmd1->SetParameterName("HistFile",false,false);
+
   arbintCmd1 = new G4UIcmdWithAString("/gps/hist/inter",this);
   arbintCmd1->SetGuidance("Sets the interpolation method for arbitrary distribution.");
-  arbintCmd1->SetParameterName("int",true,true);
+  arbintCmd1->SetParameterName("int",false,false);
   arbintCmd1->SetDefaultValue("Lin");
   arbintCmd1->SetCandidates("Lin Log Exp Spline");
@@ -757 +765 @@
   histnameCmd = new G4UIcmdWithAString("/gps/histname",this);
   histnameCmd->SetGuidance("Sets histogram type (obsolete!)");
-  histnameCmd->SetParameterName("HistType",true,true);
+  histnameCmd->SetParameterName("HistType",false,false);
   histnameCmd->SetDefaultValue("biasx");
   histnameCmd->SetCandidates("biasx biasy biasz biast biasp biase biaspt biaspp theta phi energy arb epn");
@@ -764 +772 @@
   resethistCmd = new G4UIcmdWithAString("/gps/resethist",this);
   resethistCmd->SetGuidance("Re-Set the histogram (obsolete!)");
-  resethistCmd->SetParameterName("HistType",true,true);
+  resethistCmd->SetParameterName("HistType",false,false);
   resethistCmd->SetDefaultValue("energy");
   resethistCmd->SetCandidates("biasx biasy biasz biast biasp biase biaspt biaspp theta phi energy arb epn");
@@ -771 +779 @@
   histpointCmd->SetGuidance("Allows user to define a histogram (obsolete!)");
   histpointCmd->SetGuidance("Enter: Ehi Weight");
-  histpointCmd->SetParameterName("Ehi","Weight","Junk",true,true);
+  histpointCmd->SetParameterName("Ehi","Weight","Junk",false,false);
   histpointCmd->SetRange("Ehi >= 0. && Weight >= 0.");
 
   arbintCmd = new G4UIcmdWithAString("/gps/arbint",this);
   arbintCmd->SetGuidance("Sets Arbitrary Interpolation type.(obsolete!) ");
-  arbintCmd->SetParameterName("int",true,true);
+  arbintCmd->SetParameterName("int",false,false);
   arbintCmd->SetDefaultValue("NULL");
   arbintCmd->SetCandidates("Lin Log Exp Spline");
@@ -870 +878 @@
   delete gradientCmd1;
   delete interceptCmd1;
+  delete arbeintCmd1;
   delete calculateCmd1;
   delete energyspecCmd1;
@@ -882 +891 @@
   delete resethistCmd1;
   delete histpointCmd1;
+  delete histfileCmd1;
   delete arbintCmd1;
 
@@ -1505 +1515 @@
       fParticleGun->GetEneDist()->SetInterCept(interceptCmd1->GetNewDoubleValue(newValues));
     }
+  else if(command == arbeintCmd1)
+    {
+      fParticleGun->GetEneDist()->SetBiasAlpha(arbeintCmd1->GetNewDoubleValue(newValues));
+    }
   else if(command == calculateCmd1)
     {
@@ -1520 +1534 @@
     {
       histtype = newValues;
+    }
+  else if(command == histfileCmd1)
+    {
+      histtype = "arb";
+      fParticleGun->GetEneDist()->ArbEnergyHistoFile(newValues);
     }
   else if(command == histpointCmd1)

trunk/source/event/src/G4SPSEneDistribution.cc

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

trunk/source/event/src/G4SPSRandomGenerator.cc

@@ -63 +63 @@
 //G4SPSRandomGenerator* G4SPSRandomGenerator::instance = 0;
 
-G4SPSRandomGenerator::G4SPSRandomGenerator()
-{
-  // Initialise all variables
-
-  // Bias variables
-  XBias = false;
-  IPDFXBias = false;
-  YBias = false;
-  IPDFYBias = false;
-  ZBias = false;
-  IPDFZBias = false;
-  ThetaBias = false;
-  IPDFThetaBias = false;
-  PhiBias = false;
-  IPDFPhiBias = false;
-  EnergyBias = false;
-  IPDFEnergyBias = false;
-  PosThetaBias = false;
-  IPDFPosThetaBias = false;
-  PosPhiBias = false;
-  IPDFPosPhiBias = false;
-  bweights[0] = bweights[1] = bweights[2] = bweights[3] = bweights[4] = bweights[5] = bweights[6] = bweights[7]= bweights[8]= 1. ;
-  verbosityLevel = 0 ;
-}
-
-G4SPSRandomGenerator::~G4SPSRandomGenerator()
-{}
+G4SPSRandomGenerator::G4SPSRandomGenerator() {
+        // Initialise all variables
+
+        // Bias variables
+        XBias = false;
+        IPDFXBias = false;
+        YBias = false;
+        IPDFYBias = false;
+        ZBias = false;
+        IPDFZBias = false;
+        ThetaBias = false;
+        IPDFThetaBias = false;
+        PhiBias = false;
+        IPDFPhiBias = false;
+        EnergyBias = false;
+        IPDFEnergyBias = false;
+        PosThetaBias = false;
+        IPDFPosThetaBias = false;
+        PosPhiBias = false;
+        IPDFPosPhiBias = false;
+        bweights[0] = bweights[1] = bweights[2] = bweights[3] = bweights[4]
+                        = bweights[5] = bweights[6] = bweights[7] = bweights[8] = 1.;
+        verbosityLevel = 0;
+}
+
+G4SPSRandomGenerator::~G4SPSRandomGenerator() {
+}
 
 //G4SPSRandomGenerator* G4SPSRandomGenerator::getInstance ()
@@ -99 +99 @@
 // Biasing methods
 
-void G4SPSRandomGenerator::SetXBias(G4ThreeVector input)
-{  
-  G4double ehi, val;
-  ehi = input.x();
-  val = input.y();
-  XBiasH.InsertValues(ehi, val);
-  XBias = true;
-}
-
-void G4SPSRandomGenerator::SetYBias(G4ThreeVector input)
-{
-  G4double ehi, val;
-  ehi = input.x();
-  val = input.y();
-  YBiasH.InsertValues(ehi, val);
-  YBias = true;
-}
-
-void G4SPSRandomGenerator::SetZBias(G4ThreeVector input)
-{
-  G4double ehi, val;
-  ehi = input.x();
-  val = input.y();
-  ZBiasH.InsertValues(ehi, val);
-  ZBias = true;
-}
-
-void G4SPSRandomGenerator::SetThetaBias(G4ThreeVector input)
-{
-  G4double ehi, val;
-  ehi = input.x();
-  val = input.y();
-  ThetaBiasH.InsertValues(ehi, val);
-  ThetaBias = true;
-}
-
-void G4SPSRandomGenerator::SetPhiBias(G4ThreeVector input)
-{
-  G4double ehi, val;
-  ehi = input.x();
-  val = input.y();
-  PhiBiasH.InsertValues(ehi, val);
-  PhiBias = true;
-}
-
-void G4SPSRandomGenerator::SetEnergyBias(G4ThreeVector input)
-{
-  G4double ehi, val;
-  ehi = input.x();
-  val = input.y();
-  EnergyBiasH.InsertValues(ehi, val);
-  EnergyBias = true;
-}
-
-void G4SPSRandomGenerator::SetPosThetaBias(G4ThreeVector input)
-{
-  G4double ehi, val;
-  ehi = input.x();
-  val = input.y();
-  PosThetaBiasH.InsertValues(ehi, val);
-  PosThetaBias = true;
-}
-
-void G4SPSRandomGenerator::SetPosPhiBias(G4ThreeVector input)
-{
-  G4double ehi, val;
-  ehi = input.x();
-  val = input.y();
-  PosPhiBiasH.InsertValues(ehi, val);
-  PosPhiBias = true;
-}
-
-void G4SPSRandomGenerator::ReSetHist(G4String atype)
-{
-  if ( atype == "biasx") {
-    XBias = false ;
-    IPDFXBias = false;
-    XBiasH = IPDFXBiasH = ZeroPhysVector ;}
-  else if ( atype == "biasy") {
-    YBias = false ;
-    IPDFYBias = false;
-    YBiasH = IPDFYBiasH = ZeroPhysVector ;} 
-  else if ( atype == "biasz") {
-    ZBias = false ;
-    IPDFZBias = false;
-    ZBiasH = IPDFZBiasH = ZeroPhysVector ;}
-  else if ( atype == "biast") {
-    ThetaBias = false ;
-    IPDFThetaBias = false;
-    ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector ;}
-  else if ( atype == "biasp") {
-    PhiBias = false ;
-    IPDFPhiBias = false;
-    PhiBiasH = IPDFPhiBiasH = ZeroPhysVector ;}
-  else if ( atype == "biase") {
-    EnergyBias = false ;
-    IPDFEnergyBias = false;
-    EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector ;}
-  else if ( atype == "biaspt") {
-    PosThetaBias = false ;
-    IPDFPosThetaBias = false;
-    PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector ;}
-  else if ( atype == "biaspp") {
-    PosPhiBias = false ;
-    IPDFPosPhiBias = false;
-    PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector ;}
-  else {
-    G4cout << "Error, histtype not accepted " << G4endl;
-  }
-}
-
-G4double G4SPSRandomGenerator::GenRandX()
-{
-  if(verbosityLevel >= 1)
-    G4cout << "In GenRandX" << G4endl;
-  if(XBias == false)
-    {
-      // X is not biased
-      G4double rndm = G4UniformRand();
-      return(rndm);
-    }
-  else
-    {
-      // X is biased
-      if(IPDFXBias == false)
-        {
-          // IPDF has not been created, so create it
-          G4double bins[1024],vals[1024], sum;
-          G4int ii;
-          G4int maxbin = G4int(XBiasH.GetVectorLength());
-          bins[0] = XBiasH.GetLowEdgeEnergy(size_t(0));
-          vals[0] = XBiasH(size_t(0));
-          sum = vals[0];
-          for(ii=1;ii<maxbin;ii++)
-            {
-              bins[ii] = XBiasH.GetLowEdgeEnergy(size_t(ii));
-              vals[ii] = XBiasH(size_t(ii)) + vals[ii-1];
-              sum = sum + XBiasH(size_t(ii));
-            }
-
-          for(ii=0;ii<maxbin;ii++)
-            {
-              vals[ii] = vals[ii]/sum;
-              IPDFXBiasH.InsertValues(bins[ii], vals[ii]);
-            }
-          // Make IPDFXBias = true
-          IPDFXBias = true;
-        }
-      // IPDF has been create so carry on
-      G4double rndm = G4UniformRand();
-
-      // Calculate the weighting: Find the bin that the determined
-      // rndm is in and the weigthing will be the difference in the
-      // natural probability (from the x-axis) divided by the 
-      // difference in the biased probability (the area).
-      size_t numberOfBin = IPDFXBiasH.GetVectorLength();
-      G4int biasn1 = 0;
-      G4int biasn2 = numberOfBin/2;
-      G4int biasn3 = numberOfBin - 1;
-      while (biasn1 != biasn3 - 1) {
-      if (rndm > IPDFXBiasH(biasn2))
-         biasn1 = biasn2;
-      else
-         biasn3 = biasn2;
-      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
-      }
-      // retrieve the areas and then the x-axis values
-      bweights[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1);
-      G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
-      G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2));
-      G4double NatProb = xaxisu - xaxisl;
-      //G4cout << "X Bin weight " << bweights[0] << " " << rndm << G4endl;
-      //G4cout << "lower and upper xaxis vals "<<xaxisl<<" "<<xaxisu<<G4endl;
-      bweights[0] = NatProb/bweights[0];
-      if(verbosityLevel >= 1)
-        G4cout << "X bin weight " << bweights[0] << " " << rndm << G4endl;
-      return(IPDFXBiasH.GetEnergy(rndm));
-    }
-}
-
-G4double G4SPSRandomGenerator::GenRandY()
-{
-  if(verbosityLevel >= 1)
-    G4cout << "In GenRandY" << G4endl;
-  if(YBias == false)
-    {
-      // Y is not biased
-      G4double rndm = G4UniformRand();
-      return(rndm);
-    }
-  else
-    {
-      // Y is biased
-      if(IPDFYBias == false)
-        {
-          // IPDF has not been created, so create it
-          G4double bins[1024],vals[1024], sum;
-          G4int ii;
-          G4int maxbin = G4int(YBiasH.GetVectorLength());
-          bins[0] = YBiasH.GetLowEdgeEnergy(size_t(0));
-          vals[0] = YBiasH(size_t(0));
-          sum = vals[0];
-          for(ii=1;ii<maxbin;ii++)
-            {
-              bins[ii] = YBiasH.GetLowEdgeEnergy(size_t(ii));
-              vals[ii] = YBiasH(size_t(ii)) + vals[ii-1];
-              sum = sum + YBiasH(size_t(ii));
-            }
-
-          for(ii=0;ii<maxbin;ii++)
-            {
-              vals[ii] = vals[ii]/sum;
-              IPDFYBiasH.InsertValues(bins[ii], vals[ii]);
-            }
-          // Make IPDFYBias = true
-          IPDFYBias = true;
-        }
-      // IPDF has been create so carry on
-      G4double rndm = G4UniformRand();
-      size_t numberOfBin = IPDFYBiasH.GetVectorLength();
-      G4int biasn1 = 0;
-      G4int biasn2 = numberOfBin/2;
-      G4int biasn3 = numberOfBin - 1;
-      while (biasn1 != biasn3 - 1) {
-      if (rndm > IPDFYBiasH(biasn2))
-         biasn1 = biasn2;
-      else
-         biasn3 = biasn2;
-      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
-      }
-      bweights[1] =  IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
-      G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
-      G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2));
-      G4double NatProb = xaxisu - xaxisl;
-      bweights[1] = NatProb/bweights[1];
-      if(verbosityLevel >= 1)
-        G4cout << "Y bin weight " << bweights[1] << " " << rndm << G4endl;
-      return(IPDFYBiasH.GetEnergy(rndm));
-    }
-}
-
-G4double G4SPSRandomGenerator::GenRandZ()
-{
-  if(verbosityLevel >= 1)
-    G4cout << "In GenRandZ" << G4endl;
-  if(ZBias == false)
-    {
-      // Z is not biased
-      G4double rndm = G4UniformRand();
-      return(rndm);
-    }
-  else
-    {
-      // Z is biased
-      if(IPDFZBias == false)
-        {
-          // IPDF has not been created, so create it
-          G4double bins[1024],vals[1024], sum;
-          G4int ii;
-          G4int maxbin = G4int(ZBiasH.GetVectorLength());
-          bins[0] = ZBiasH.GetLowEdgeEnergy(size_t(0));
-          vals[0] = ZBiasH(size_t(0));
-          sum = vals[0];
-          for(ii=1;ii<maxbin;ii++)
-            {
-              bins[ii] = ZBiasH.GetLowEdgeEnergy(size_t(ii));
-              vals[ii] = ZBiasH(size_t(ii)) + vals[ii-1];
-              sum = sum + ZBiasH(size_t(ii));
-            }
-
-          for(ii=0;ii<maxbin;ii++)
-            {
-              vals[ii] = vals[ii]/sum;
-              IPDFZBiasH.InsertValues(bins[ii], vals[ii]);
-            }
-          // Make IPDFZBias = true
-          IPDFZBias = true;
-        }
-      // IPDF has been create so carry on
-      G4double rndm = G4UniformRand();
-      //      size_t weight_bin_no = IPDFZBiasH.FindValueBinLocation(rndm);
-      size_t numberOfBin = IPDFZBiasH.GetVectorLength();
-      G4int biasn1 = 0;
-      G4int biasn2 = numberOfBin/2;
-      G4int biasn3 = numberOfBin - 1;
-      while (biasn1 != biasn3 - 1) {
-      if (rndm > IPDFZBiasH(biasn2))
-         biasn1 = biasn2;
-      else
-         biasn3 = biasn2;
-      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
-      }
-      bweights[2] =  IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
-      G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
-      G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2));
-      G4double NatProb = xaxisu - xaxisl;
-      bweights[2] = NatProb/bweights[2];
-      if(verbosityLevel >= 1)
-        G4cout << "Z bin weight " << bweights[2] << " " << rndm << G4endl;
-      return(IPDFZBiasH.GetEnergy(rndm));
-    }
-}
-
-G4double G4SPSRandomGenerator::GenRandTheta()
-{
-  if(verbosityLevel >= 1)
-    {
-      G4cout << "In GenRandTheta" << G4endl;
-      G4cout << "Verbosity " << verbosityLevel << G4endl;
-    }
-  if(ThetaBias == false)
-    {
-      // Theta is not biased
-      G4double rndm = G4UniformRand();
-      return(rndm);
-    }
-  else
-    {
-      // Theta is biased
-      if(IPDFThetaBias == false)
-        {
-          // IPDF has not been created, so create it
-          G4double bins[1024],vals[1024], sum;
-          G4int ii;
-          G4int maxbin = G4int(ThetaBiasH.GetVectorLength());
-          bins[0] = ThetaBiasH.GetLowEdgeEnergy(size_t(0));
-          vals[0] = ThetaBiasH(size_t(0));
-          sum = vals[0];
-          for(ii=1;ii<maxbin;ii++)
-            {
-              bins[ii] = ThetaBiasH.GetLowEdgeEnergy(size_t(ii));
-              vals[ii] = ThetaBiasH(size_t(ii)) + vals[ii-1];
-              sum = sum + ThetaBiasH(size_t(ii));
-            }
-
-          for(ii=0;ii<maxbin;ii++)
-            {
-              vals[ii] = vals[ii]/sum;
-              IPDFThetaBiasH.InsertValues(bins[ii], vals[ii]);
-            }
-          // Make IPDFThetaBias = true
-          IPDFThetaBias = true;
-        }
-      // IPDF has been create so carry on
-      G4double rndm = G4UniformRand();
-      //      size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
-      size_t numberOfBin = IPDFThetaBiasH.GetVectorLength();
-      G4int biasn1 = 0;
-      G4int biasn2 = numberOfBin/2;
-      G4int biasn3 = numberOfBin - 1;
-      while (biasn1 != biasn3 - 1) {
-      if (rndm > IPDFThetaBiasH(biasn2))
-         biasn1 = biasn2;
-      else
-         biasn3 = biasn2;
-      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
-      }
-      bweights[3] =  IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
-      G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
-      G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
-      G4double NatProb = xaxisu - xaxisl;
-      bweights[3] = NatProb/bweights[3];
-      if(verbosityLevel >= 1)
-        G4cout << "Theta bin weight " << bweights[3] << " " << rndm << G4endl;
-      return(IPDFThetaBiasH.GetEnergy(rndm));
-    }
-}
-
-G4double G4SPSRandomGenerator::GenRandPhi()
-{
-  if(verbosityLevel >= 1)
-    G4cout << "In GenRandPhi" << G4endl;
-  if(PhiBias == false)
-    {
-      // Phi is not biased
-      G4double rndm = G4UniformRand();
-      return(rndm);
-    }
-  else
-    {
-      // Phi is biased
-      if(IPDFPhiBias == false)
-        {
-          // IPDF has not been created, so create it
-          G4double bins[1024],vals[1024], sum;
-          G4int ii;
-          G4int maxbin = G4int(PhiBiasH.GetVectorLength());
-          bins[0] = PhiBiasH.GetLowEdgeEnergy(size_t(0));
-          vals[0] = PhiBiasH(size_t(0));
-          sum = vals[0];
-          for(ii=1;ii<maxbin;ii++)
-            {
-              bins[ii] = PhiBiasH.GetLowEdgeEnergy(size_t(ii));
-              vals[ii] = PhiBiasH(size_t(ii)) + vals[ii-1];
-              sum = sum + PhiBiasH(size_t(ii));
-            }
-
-          for(ii=0;ii<maxbin;ii++)
-            {
-              vals[ii] = vals[ii]/sum;
-              IPDFPhiBiasH.InsertValues(bins[ii], vals[ii]);
-            }
-          // Make IPDFPhiBias = true
-          IPDFPhiBias = true;
-        }
-      // IPDF has been create so carry on
-      G4double rndm = G4UniformRand();
-      //      size_t weight_bin_no = IPDFPhiBiasH.FindValueBinLocation(rndm);
-      size_t numberOfBin = IPDFPhiBiasH.GetVectorLength();
-      G4int biasn1 = 0;
-      G4int biasn2 = numberOfBin/2;
-      G4int biasn3 = numberOfBin - 1;
-      while (biasn1 != biasn3 - 1) {
-      if (rndm > IPDFPhiBiasH(biasn2))
-         biasn1 = biasn2;
-      else
-         biasn3 = biasn2;
-      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
-      }
-      bweights[4] =  IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
-      G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
-      G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
-      G4double NatProb = xaxisu - xaxisl;
-      bweights[4] = NatProb/bweights[4];
-      if(verbosityLevel >= 1)
-        G4cout << "Phi bin weight " << bweights[4] << " " << rndm << G4endl;
-      return(IPDFPhiBiasH.GetEnergy(rndm));
-    }
-}
-
-G4double G4SPSRandomGenerator::GenRandEnergy()
-{
-  if(verbosityLevel >= 1)
-    G4cout << "In GenRandEnergy" << G4endl;
-  if(EnergyBias == false)
-    {
-      // Energy is not biased
-      G4double rndm = G4UniformRand();
-      return(rndm);
-    }
-  else {
-    // ENERGY is biased
-    if(IPDFEnergyBias == false) {
-      // IPDF has not been created, so create it
-      G4double bins[1024],vals[1024], sum;
-      G4int ii;
-      G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
-      bins[0] = EnergyBiasH.GetLowEdgeEnergy(size_t(0));
-      vals[0] = EnergyBiasH(size_t(0));
-      sum = vals[0];
-      for(ii=1;ii<maxbin;ii++) {
-        bins[ii] = EnergyBiasH.GetLowEdgeEnergy(size_t(ii));
-        vals[ii] = EnergyBiasH(size_t(ii)) + vals[ii-1];
-        sum = sum + EnergyBiasH(size_t(ii));
-      }
-      
-      for(ii=0;ii<maxbin;ii++) {
-        vals[ii] = vals[ii]/sum;
-        IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
-      }
-      // Make IPDFEnergyBias = true
-      IPDFEnergyBias = true;
-    }
-    // IPDF has been create so carry on
-    G4double rndm = G4UniformRand();
-    //  size_t weight_bin_no = IPDFEnergyBiasH.FindValueBinLocation(rndm);
-    size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
-    G4int biasn1 = 0;
-    G4int biasn2 = numberOfBin/2;
-    G4int biasn3 = numberOfBin - 1;
-    while (biasn1 != biasn3 - 1) {
-      if (rndm > IPDFEnergyBiasH(biasn2))
-        biasn1 = biasn2;
-      else
-        biasn3 = biasn2;
-      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
-    }
-    bweights[5] =  IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
-    G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
-    G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2));
-    G4double NatProb = xaxisu - xaxisl;
-    bweights[5] = NatProb/bweights[5];
-    if(verbosityLevel >= 1)
-      G4cout << "Energy bin weight " << bweights[5] << " " << rndm << G4endl;
-    return(IPDFEnergyBiasH.GetEnergy(rndm));
-  }
-}
-
-
-G4double G4SPSRandomGenerator::GenRandPosTheta()
-{
-  if(verbosityLevel >= 1)
-    {
-      G4cout << "In GenRandPosTheta" << G4endl;
-      G4cout << "Verbosity " << verbosityLevel << G4endl;
-    }
-  if(PosThetaBias == false)
-    {
-      // Theta is not biased
-      G4double rndm = G4UniformRand();
-      return(rndm);
-    }
-  else
-    {
-      // Theta is biased
-      if(IPDFPosThetaBias == false)
-        {
-          // IPDF has not been created, so create it
-          G4double bins[1024],vals[1024], sum;
-          G4int ii;
-          G4int maxbin = G4int(PosThetaBiasH.GetVectorLength());
-          bins[0] = PosThetaBiasH.GetLowEdgeEnergy(size_t(0));
-          vals[0] = PosThetaBiasH(size_t(0));
-          sum = vals[0];
-          for(ii=1;ii<maxbin;ii++)
-            {
-              bins[ii] = PosThetaBiasH.GetLowEdgeEnergy(size_t(ii));
-              vals[ii] = PosThetaBiasH(size_t(ii)) + vals[ii-1];
-              sum = sum + PosThetaBiasH(size_t(ii));
-            }
-
-          for(ii=0;ii<maxbin;ii++)
-            {
-              vals[ii] = vals[ii]/sum;
-              IPDFPosThetaBiasH.InsertValues(bins[ii], vals[ii]);
-            }
-          // Make IPDFThetaBias = true
-          IPDFPosThetaBias = true;
-        }
-      // IPDF has been create so carry on
-      G4double rndm = G4UniformRand();
-      //      size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
-      size_t numberOfBin = IPDFPosThetaBiasH.GetVectorLength();
-      G4int biasn1 = 0;
-      G4int biasn2 = numberOfBin/2;
-      G4int biasn3 = numberOfBin - 1;
-      while (biasn1 != biasn3 - 1) {
-      if (rndm > IPDFPosThetaBiasH(biasn2))
-         biasn1 = biasn2;
-      else
-         biasn3 = biasn2;
-      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
-      }
-      bweights[6] =  IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
-      G4double xaxisl = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
-      G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
-      G4double NatProb = xaxisu - xaxisl;
-      bweights[6] = NatProb/bweights[6];
-      if(verbosityLevel >= 1)
-        G4cout << "PosTheta bin weight " << bweights[6] << " " << rndm << G4endl;
-      return(IPDFPosThetaBiasH.GetEnergy(rndm));
-    }
-}
-
-G4double G4SPSRandomGenerator::GenRandPosPhi()
-{
-  if(verbosityLevel >= 1)
-    G4cout << "In GenRandPosPhi" << G4endl;
-  if(PosPhiBias == false)
-    {
-      // PosPhi is not biased
-      G4double rndm = G4UniformRand();
-      return(rndm);
-    }
-  else
-    {
-      // PosPhi is biased
-      if(IPDFPosPhiBias == false)
-        {
-          // IPDF has not been created, so create it
-          G4double bins[1024],vals[1024], sum;
-          G4int ii;
-          G4int maxbin = G4int(PosPhiBiasH.GetVectorLength());
-          bins[0] = PosPhiBiasH.GetLowEdgeEnergy(size_t(0));
-          vals[0] = PosPhiBiasH(size_t(0));
-          sum = vals[0];
-          for(ii=1;ii<maxbin;ii++)
-            {
-              bins[ii] = PosPhiBiasH.GetLowEdgeEnergy(size_t(ii));
-              vals[ii] = PosPhiBiasH(size_t(ii)) + vals[ii-1];
-              sum = sum + PosPhiBiasH(size_t(ii));
-            }
-
-          for(ii=0;ii<maxbin;ii++)
-            {
-              vals[ii] = vals[ii]/sum;
-              IPDFPosPhiBiasH.InsertValues(bins[ii], vals[ii]);
-            }
-          // Make IPDFPosPhiBias = true
-          IPDFPosPhiBias = true;
-        }
-      // IPDF has been create so carry on
-      G4double rndm = G4UniformRand();
-      //      size_t weight_bin_no = IPDFPosPhiBiasH.FindValueBinLocation(rndm);
-      size_t numberOfBin = IPDFPosPhiBiasH.GetVectorLength();
-      G4int biasn1 = 0;
-      G4int biasn2 = numberOfBin/2;
-      G4int biasn3 = numberOfBin - 1;
-      while (biasn1 != biasn3 - 1) {
-      if (rndm > IPDFPosPhiBiasH(biasn2))
-         biasn1 = biasn2;
-      else
-         biasn3 = biasn2;
-      biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
-      }
-      bweights[7] =  IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
-      G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
-      G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
-      G4double NatProb = xaxisu - xaxisl;
-      bweights[7] = NatProb/bweights[7];
-      if(verbosityLevel >= 1)
-        G4cout << "PosPhi bin weight " << bweights[7] << " " << rndm << G4endl;
-      return(IPDFPosPhiBiasH.GetEnergy(rndm));
-    }
-}
-
-
-
-
-
+void G4SPSRandomGenerator::SetXBias(G4ThreeVector input) {
+        G4double ehi, val;
+        ehi = input.x();
+        val = input.y();
+        XBiasH.InsertValues(ehi, val);
+        XBias = true;
+}
+
+void G4SPSRandomGenerator::SetYBias(G4ThreeVector input) {
+        G4double ehi, val;
+        ehi = input.x();
+        val = input.y();
+        YBiasH.InsertValues(ehi, val);
+        YBias = true;
+}
+
+void G4SPSRandomGenerator::SetZBias(G4ThreeVector input) {
+        G4double ehi, val;
+        ehi = input.x();
+        val = input.y();
+        ZBiasH.InsertValues(ehi, val);
+        ZBias = true;
+}
+
+void G4SPSRandomGenerator::SetThetaBias(G4ThreeVector input) {
+        G4double ehi, val;
+        ehi = input.x();
+        val = input.y();
+        ThetaBiasH.InsertValues(ehi, val);
+        ThetaBias = true;
+}
+
+void G4SPSRandomGenerator::SetPhiBias(G4ThreeVector input) {
+        G4double ehi, val;
+        ehi = input.x();
+        val = input.y();
+        PhiBiasH.InsertValues(ehi, val);
+        PhiBias = true;
+}
+
+void G4SPSRandomGenerator::SetEnergyBias(G4ThreeVector input) {
+        G4double ehi, val;
+        ehi = input.x();
+        val = input.y();
+        EnergyBiasH.InsertValues(ehi, val);
+        EnergyBias = true;
+}
+
+void G4SPSRandomGenerator::SetPosThetaBias(G4ThreeVector input) {
+        G4double ehi, val;
+        ehi = input.x();
+        val = input.y();
+        PosThetaBiasH.InsertValues(ehi, val);
+        PosThetaBias = true;
+}
+
+void G4SPSRandomGenerator::SetPosPhiBias(G4ThreeVector input) {
+        G4double ehi, val;
+        ehi = input.x();
+        val = input.y();
+        PosPhiBiasH.InsertValues(ehi, val);
+        PosPhiBias = true;
+}
+
+void G4SPSRandomGenerator::ReSetHist(G4String atype) {
+        if (atype == "biasx") {
+                XBias = false;
+                IPDFXBias = false;
+                XBiasH = IPDFXBiasH = ZeroPhysVector;
+        } else if (atype == "biasy") {
+                YBias = false;
+                IPDFYBias = false;
+                YBiasH = IPDFYBiasH = ZeroPhysVector;
+        } else if (atype == "biasz") {
+                ZBias = false;
+                IPDFZBias = false;
+                ZBiasH = IPDFZBiasH = ZeroPhysVector;
+        } else if (atype == "biast") {
+                ThetaBias = false;
+                IPDFThetaBias = false;
+                ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector;
+        } else if (atype == "biasp") {
+                PhiBias = false;
+                IPDFPhiBias = false;
+                PhiBiasH = IPDFPhiBiasH = ZeroPhysVector;
+        } else if (atype == "biase") {
+                EnergyBias = false;
+                IPDFEnergyBias = false;
+                EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector;
+        } else if (atype == "biaspt") {
+                PosThetaBias = false;
+                IPDFPosThetaBias = false;
+                PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector;
+        } else if (atype == "biaspp") {
+                PosPhiBias = false;
+                IPDFPosPhiBias = false;
+                PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector;
+        } else {
+                G4cout << "Error, histtype not accepted " << G4endl;
+        }
+}
+
+G4double G4SPSRandomGenerator::GenRandX() {
+        if (verbosityLevel >= 1)
+                G4cout << "In GenRandX" << G4endl;
+        if (XBias == false) {
+                // X is not biased
+                G4double rndm = G4UniformRand();
+                return (rndm);
+        } else {
+                // X is biased
+                if (IPDFXBias == false) {
+                        // IPDF has not been created, so create it
+                        G4double bins[1024], vals[1024], sum;
+                        G4int ii;
+                        G4int maxbin = G4int(XBiasH.GetVectorLength());
+                        bins[0] = XBiasH.GetLowEdgeEnergy(size_t(0));
+                        vals[0] = XBiasH(size_t(0));
+                        sum = vals[0];
+                        for (ii = 1; ii < maxbin; ii++) {
+                                bins[ii] = XBiasH.GetLowEdgeEnergy(size_t(ii));
+                                vals[ii] = XBiasH(size_t(ii)) + vals[ii - 1];
+                                sum = sum + XBiasH(size_t(ii));
+                        }
+
+                        for (ii = 0; ii < maxbin; ii++) {
+                                vals[ii] = vals[ii] / sum;
+                                IPDFXBiasH.InsertValues(bins[ii], vals[ii]);
+                        }
+                        // Make IPDFXBias = true
+                        IPDFXBias = true;
+                }
+                // IPDF has been create so carry on
+                G4double rndm = G4UniformRand();
+
+                // Calculate the weighting: Find the bin that the determined
+                // rndm is in and the weigthing will be the difference in the
+                // natural probability (from the x-axis) divided by the
+                // difference in the biased probability (the area).
+                size_t numberOfBin = IPDFXBiasH.GetVectorLength();
+                G4int biasn1 = 0;
+                G4int biasn2 = numberOfBin / 2;
+                G4int biasn3 = numberOfBin - 1;
+                while (biasn1 != biasn3 - 1) {
+                        if (rndm > IPDFXBiasH(biasn2))
+                                biasn1 = biasn2;
+                        else
+                                biasn3 = biasn2;
+                        biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
+                }
+                // retrieve the areas and then the x-axis values
+                bweights[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1);
+                G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
+                G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2));
+                G4double NatProb = xaxisu - xaxisl;
+                //G4cout << "X Bin weight " << bweights[0] << " " << rndm << G4endl;
+                //G4cout << "lower and upper xaxis vals "<<xaxisl<<" "<<xaxisu<<G4endl;
+                bweights[0] = NatProb / bweights[0];
+                if (verbosityLevel >= 1)
+                        G4cout << "X bin weight " << bweights[0] << " " << rndm << G4endl;
+                return (IPDFXBiasH.GetEnergy(rndm));
+        }
+}
+
+G4double G4SPSRandomGenerator::GenRandY() {
+        if (verbosityLevel >= 1)
+                G4cout << "In GenRandY" << G4endl;
+        if (YBias == false) {
+                // Y is not biased
+                G4double rndm = G4UniformRand();
+                return (rndm);
+        } else {
+                // Y is biased
+                if (IPDFYBias == false) {
+                        // IPDF has not been created, so create it
+                        G4double bins[1024], vals[1024], sum;
+                        G4int ii;
+                        G4int maxbin = G4int(YBiasH.GetVectorLength());
+                        bins[0] = YBiasH.GetLowEdgeEnergy(size_t(0));
+                        vals[0] = YBiasH(size_t(0));
+                        sum = vals[0];
+                        for (ii = 1; ii < maxbin; ii++) {
+                                bins[ii] = YBiasH.GetLowEdgeEnergy(size_t(ii));
+                                vals[ii] = YBiasH(size_t(ii)) + vals[ii - 1];
+                                sum = sum + YBiasH(size_t(ii));
+                        }
+
+                        for (ii = 0; ii < maxbin; ii++) {
+                                vals[ii] = vals[ii] / sum;
+                                IPDFYBiasH.InsertValues(bins[ii], vals[ii]);
+                        }
+                        // Make IPDFYBias = true
+                        IPDFYBias = true;
+                }
+                // IPDF has been create so carry on
+                G4double rndm = G4UniformRand();
+                size_t numberOfBin = IPDFYBiasH.GetVectorLength();
+                G4int biasn1 = 0;
+                G4int biasn2 = numberOfBin / 2;
+                G4int biasn3 = numberOfBin - 1;
+                while (biasn1 != biasn3 - 1) {
+                        if (rndm > IPDFYBiasH(biasn2))
+                                biasn1 = biasn2;
+                        else
+                                biasn3 = biasn2;
+                        biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
+                }
+                bweights[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
+                G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
+                G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2));
+                G4double NatProb = xaxisu - xaxisl;
+                bweights[1] = NatProb / bweights[1];
+                if (verbosityLevel >= 1)
+                        G4cout << "Y bin weight " << bweights[1] << " " << rndm << G4endl;
+                return (IPDFYBiasH.GetEnergy(rndm));
+        }
+}
+
+G4double G4SPSRandomGenerator::GenRandZ() {
+        if (verbosityLevel >= 1)
+                G4cout << "In GenRandZ" << G4endl;
+        if (ZBias == false) {
+                // Z is not biased
+                G4double rndm = G4UniformRand();
+                return (rndm);
+        } else {
+                // Z is biased
+                if (IPDFZBias == false) {
+                        // IPDF has not been created, so create it
+                        G4double bins[1024], vals[1024], sum;
+                        G4int ii;
+                        G4int maxbin = G4int(ZBiasH.GetVectorLength());
+                        bins[0] = ZBiasH.GetLowEdgeEnergy(size_t(0));
+                        vals[0] = ZBiasH(size_t(0));
+                        sum = vals[0];
+                        for (ii = 1; ii < maxbin; ii++) {
+                                bins[ii] = ZBiasH.GetLowEdgeEnergy(size_t(ii));
+                                vals[ii] = ZBiasH(size_t(ii)) + vals[ii - 1];
+                                sum = sum + ZBiasH(size_t(ii));
+                        }
+
+                        for (ii = 0; ii < maxbin; ii++) {
+                                vals[ii] = vals[ii] / sum;
+                                IPDFZBiasH.InsertValues(bins[ii], vals[ii]);
+                        }
+                        // Make IPDFZBias = true
+                        IPDFZBias = true;
+                }
+                // IPDF has been create so carry on
+                G4double rndm = G4UniformRand();
+                //      size_t weight_bin_no = IPDFZBiasH.FindValueBinLocation(rndm);
+                size_t numberOfBin = IPDFZBiasH.GetVectorLength();
+                G4int biasn1 = 0;
+                G4int biasn2 = numberOfBin / 2;
+                G4int biasn3 = numberOfBin - 1;
+                while (biasn1 != biasn3 - 1) {
+                        if (rndm > IPDFZBiasH(biasn2))
+                                biasn1 = biasn2;
+                        else
+                                biasn3 = biasn2;
+                        biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
+                }
+                bweights[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
+                G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
+                G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2));
+                G4double NatProb = xaxisu - xaxisl;
+                bweights[2] = NatProb / bweights[2];
+                if (verbosityLevel >= 1)
+                        G4cout << "Z bin weight " << bweights[2] << " " << rndm << G4endl;
+                return (IPDFZBiasH.GetEnergy(rndm));
+        }
+}
+
+G4double G4SPSRandomGenerator::GenRandTheta() {
+        if (verbosityLevel >= 1) {
+                G4cout << "In GenRandTheta" << G4endl;
+                G4cout << "Verbosity " << verbosityLevel << G4endl;
+        }
+        if (ThetaBias == false) {
+                // Theta is not biased
+                G4double rndm = G4UniformRand();
+                return (rndm);
+        } else {
+                // Theta is biased
+                if (IPDFThetaBias == false) {
+                        // IPDF has not been created, so create it
+                        G4double bins[1024], vals[1024], sum;
+                        G4int ii;
+                        G4int maxbin = G4int(ThetaBiasH.GetVectorLength());
+                        bins[0] = ThetaBiasH.GetLowEdgeEnergy(size_t(0));
+                        vals[0] = ThetaBiasH(size_t(0));
+                        sum = vals[0];
+                        for (ii = 1; ii < maxbin; ii++) {
+                                bins[ii] = ThetaBiasH.GetLowEdgeEnergy(size_t(ii));
+                                vals[ii] = ThetaBiasH(size_t(ii)) + vals[ii - 1];
+                                sum = sum + ThetaBiasH(size_t(ii));
+                        }
+
+                        for (ii = 0; ii < maxbin; ii++) {
+                                vals[ii] = vals[ii] / sum;
+                                IPDFThetaBiasH.InsertValues(bins[ii], vals[ii]);
+                        }
+                        // Make IPDFThetaBias = true
+                        IPDFThetaBias = true;
+                }
+                // IPDF has been create so carry on
+                G4double rndm = G4UniformRand();
+                //      size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
+                size_t numberOfBin = IPDFThetaBiasH.GetVectorLength();
+                G4int biasn1 = 0;
+                G4int biasn2 = numberOfBin / 2;
+                G4int biasn3 = numberOfBin - 1;
+                while (biasn1 != biasn3 - 1) {
+                        if (rndm > IPDFThetaBiasH(biasn2))
+                                biasn1 = biasn2;
+                        else
+                                biasn3 = biasn2;
+                        biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
+                }
+                bweights[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
+                G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
+                G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
+                G4double NatProb = xaxisu - xaxisl;
+                bweights[3] = NatProb / bweights[3];
+                if (verbosityLevel >= 1)
+                        G4cout << "Theta bin weight " << bweights[3] << " " << rndm
+                                        << G4endl;
+                return (IPDFThetaBiasH.GetEnergy(rndm));
+        }
+}
+
+G4double G4SPSRandomGenerator::GenRandPhi() {
+        if (verbosityLevel >= 1)
+                G4cout << "In GenRandPhi" << G4endl;
+        if (PhiBias == false) {
+                // Phi is not biased
+                G4double rndm = G4UniformRand();
+                return (rndm);
+        } else {
+                // Phi is biased
+                if (IPDFPhiBias == false) {
+                        // IPDF has not been created, so create it
+                        G4double bins[1024], vals[1024], sum;
+                        G4int ii;
+                        G4int maxbin = G4int(PhiBiasH.GetVectorLength());
+                        bins[0] = PhiBiasH.GetLowEdgeEnergy(size_t(0));
+                        vals[0] = PhiBiasH(size_t(0));
+                        sum = vals[0];
+                        for (ii = 1; ii < maxbin; ii++) {
+                                bins[ii] = PhiBiasH.GetLowEdgeEnergy(size_t(ii));
+                                vals[ii] = PhiBiasH(size_t(ii)) + vals[ii - 1];
+                                sum = sum + PhiBiasH(size_t(ii));
+                        }
+
+                        for (ii = 0; ii < maxbin; ii++) {
+                                vals[ii] = vals[ii] / sum;
+                                IPDFPhiBiasH.InsertValues(bins[ii], vals[ii]);
+                        }
+                        // Make IPDFPhiBias = true
+                        IPDFPhiBias = true;
+                }
+                // IPDF has been create so carry on
+                G4double rndm = G4UniformRand();
+                //      size_t weight_bin_no = IPDFPhiBiasH.FindValueBinLocation(rndm);
+                size_t numberOfBin = IPDFPhiBiasH.GetVectorLength();
+                G4int biasn1 = 0;
+                G4int biasn2 = numberOfBin / 2;
+                G4int biasn3 = numberOfBin - 1;
+                while (biasn1 != biasn3 - 1) {
+                        if (rndm > IPDFPhiBiasH(biasn2))
+                                biasn1 = biasn2;
+                        else
+                                biasn3 = biasn2;
+                        biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
+                }
+                bweights[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
+                G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
+                G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
+                G4double NatProb = xaxisu - xaxisl;
+                bweights[4] = NatProb / bweights[4];
+                if (verbosityLevel >= 1)
+                        G4cout << "Phi bin weight " << bweights[4] << " " << rndm << G4endl;
+                return (IPDFPhiBiasH.GetEnergy(rndm));
+        }
+}
+
+G4double G4SPSRandomGenerator::GenRandEnergy() {
+        if (verbosityLevel >= 1)
+                G4cout << "In GenRandEnergy" << G4endl;
+        if (EnergyBias == false) {
+                // Energy is not biased
+                G4double rndm = G4UniformRand();
+                return (rndm);
+        } else {
+                // ENERGY is biased
+                if (IPDFEnergyBias == false) {
+                        // IPDF has not been created, so create it
+                        G4double bins[1024], vals[1024], sum;
+                        G4int ii;
+                        G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
+                        bins[0] = EnergyBiasH.GetLowEdgeEnergy(size_t(0));
+                        vals[0] = EnergyBiasH(size_t(0));
+                        sum = vals[0];
+                        for (ii = 1; ii < maxbin; ii++) {
+                                bins[ii] = EnergyBiasH.GetLowEdgeEnergy(size_t(ii));
+                                vals[ii] = EnergyBiasH(size_t(ii)) + vals[ii - 1];
+                                sum = sum + EnergyBiasH(size_t(ii));
+                        }
+                        IPDFEnergyBiasH = ZeroPhysVector;
+                        for (ii = 0; ii < maxbin; ii++) {
+                                vals[ii] = vals[ii] / sum;
+                                IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
+                        }
+                        // Make IPDFEnergyBias = true
+                        IPDFEnergyBias = true;
+                }
+                // IPDF has been create so carry on
+                G4double rndm = G4UniformRand();
+                //  size_t weight_bin_no = IPDFEnergyBiasH.FindValueBinLocation(rndm);
+                size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
+                G4int biasn1 = 0;
+                G4int biasn2 = numberOfBin / 2;
+                G4int biasn3 = numberOfBin - 1;
+                while (biasn1 != biasn3 - 1) {
+                        if (rndm > IPDFEnergyBiasH(biasn2))
+                                biasn1 = biasn2;
+                        else
+                                biasn3 = biasn2;
+                        biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
+                }
+                bweights[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
+                G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
+                G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2));
+                G4double NatProb = xaxisu - xaxisl;
+                bweights[5] = NatProb / bweights[5];
+                if (verbosityLevel >= 1)
+                        G4cout << "Energy bin weight " << bweights[5] << " " << rndm
+                                        << G4endl;
+                return (IPDFEnergyBiasH.GetEnergy(rndm));
+        }
+}
+
+G4double G4SPSRandomGenerator::GenRandPosTheta() {
+        if (verbosityLevel >= 1) {
+                G4cout << "In GenRandPosTheta" << G4endl;
+                G4cout << "Verbosity " << verbosityLevel << G4endl;
+        }
+        if (PosThetaBias == false) {
+                // Theta is not biased
+                G4double rndm = G4UniformRand();
+                return (rndm);
+        } else {
+                // Theta is biased
+                if (IPDFPosThetaBias == false) {
+                        // IPDF has not been created, so create it
+                        G4double bins[1024], vals[1024], sum;
+                        G4int ii;
+                        G4int maxbin = G4int(PosThetaBiasH.GetVectorLength());
+                        bins[0] = PosThetaBiasH.GetLowEdgeEnergy(size_t(0));
+                        vals[0] = PosThetaBiasH(size_t(0));
+                        sum = vals[0];
+                        for (ii = 1; ii < maxbin; ii++) {
+                                bins[ii] = PosThetaBiasH.GetLowEdgeEnergy(size_t(ii));
+                                vals[ii] = PosThetaBiasH(size_t(ii)) + vals[ii - 1];
+                                sum = sum + PosThetaBiasH(size_t(ii));
+                        }
+
+                        for (ii = 0; ii < maxbin; ii++) {
+                                vals[ii] = vals[ii] / sum;
+                                IPDFPosThetaBiasH.InsertValues(bins[ii], vals[ii]);
+                        }
+                        // Make IPDFThetaBias = true
+                        IPDFPosThetaBias = true;
+                }
+                // IPDF has been create so carry on
+                G4double rndm = G4UniformRand();
+                //      size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
+                size_t numberOfBin = IPDFPosThetaBiasH.GetVectorLength();
+                G4int biasn1 = 0;
+                G4int biasn2 = numberOfBin / 2;
+                G4int biasn3 = numberOfBin - 1;
+                while (biasn1 != biasn3 - 1) {
+                        if (rndm > IPDFPosThetaBiasH(biasn2))
+                                biasn1 = biasn2;
+                        else
+                                biasn3 = biasn2;
+                        biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
+                }
+                bweights[6] = IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
+                G4double xaxisl =
+                                IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
+                G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
+                G4double NatProb = xaxisu - xaxisl;
+                bweights[6] = NatProb / bweights[6];
+                if (verbosityLevel >= 1)
+                        G4cout << "PosTheta bin weight " << bweights[6] << " " << rndm
+                                        << G4endl;
+                return (IPDFPosThetaBiasH.GetEnergy(rndm));
+        }
+}
+
+G4double G4SPSRandomGenerator::GenRandPosPhi() {
+        if (verbosityLevel >= 1)
+                G4cout << "In GenRandPosPhi" << G4endl;
+        if (PosPhiBias == false) {
+                // PosPhi is not biased
+                G4double rndm = G4UniformRand();
+                return (rndm);
+        } else {
+                // PosPhi is biased
+                if (IPDFPosPhiBias == false) {
+                        // IPDF has not been created, so create it
+                        G4double bins[1024], vals[1024], sum;
+                        G4int ii;
+                        G4int maxbin = G4int(PosPhiBiasH.GetVectorLength());
+                        bins[0] = PosPhiBiasH.GetLowEdgeEnergy(size_t(0));
+                        vals[0] = PosPhiBiasH(size_t(0));
+                        sum = vals[0];
+                        for (ii = 1; ii < maxbin; ii++) {
+                                bins[ii] = PosPhiBiasH.GetLowEdgeEnergy(size_t(ii));
+                                vals[ii] = PosPhiBiasH(size_t(ii)) + vals[ii - 1];
+                                sum = sum + PosPhiBiasH(size_t(ii));
+                        }
+
+                        for (ii = 0; ii < maxbin; ii++) {
+                                vals[ii] = vals[ii] / sum;
+                                IPDFPosPhiBiasH.InsertValues(bins[ii], vals[ii]);
+                        }
+                        // Make IPDFPosPhiBias = true
+                        IPDFPosPhiBias = true;
+                }
+                // IPDF has been create so carry on
+                G4double rndm = G4UniformRand();
+                //      size_t weight_bin_no = IPDFPosPhiBiasH.FindValueBinLocation(rndm);
+                size_t numberOfBin = IPDFPosPhiBiasH.GetVectorLength();
+                G4int biasn1 = 0;
+                G4int biasn2 = numberOfBin / 2;
+                G4int biasn3 = numberOfBin - 1;
+                while (biasn1 != biasn3 - 1) {
+                        if (rndm > IPDFPosPhiBiasH(biasn2))
+                                biasn1 = biasn2;
+                        else
+                                biasn3 = biasn2;
+                        biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
+                }
+                bweights[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
+                G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
+                G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
+                G4double NatProb = xaxisu - xaxisl;
+                bweights[7] = NatProb / bweights[7];
+                if (verbosityLevel >= 1)
+                        G4cout << "PosPhi bin weight " << bweights[7] << " " << rndm
+                                        << G4endl;
+                return (IPDFPosPhiBiasH.GetEnergy(rndm));
+        }
+}
+

trunk/source/event/src/G4SingleParticleSource.cc

@@ -58 +58 @@
 #include "G4SingleParticleSource.hh"
 
-G4SingleParticleSource::G4SingleParticleSource()
-{
-  // Initialise all variables
-  // Position distribution Variables
+G4SingleParticleSource::G4SingleParticleSource() {
+        // Initialise all variables
+        // Position distribution Variables
 
-  NumberOfParticlesToBeGenerated = 1;
-  particle_definition = G4Geantino::GeantinoDefinition();
-  G4ThreeVector zero;
-  particle_momentum_direction = G4ParticleMomentum(1,0,0);
-  particle_energy = 1.0*MeV;
-  particle_position = zero;
-  particle_time = 0.0;
-  particle_polarization = zero;
-  particle_charge = 0.0;
-  particle_weight = 1.0;
+        NumberOfParticlesToBeGenerated = 1;
+        particle_definition = G4Geantino::GeantinoDefinition();
+        G4ThreeVector zero;
+        particle_momentum_direction = G4ParticleMomentum(1, 0, 0);
+        particle_energy = 1.0 * MeV;
+        particle_position = zero;
+        particle_time = 0.0;
+        particle_polarization = zero;
+        particle_charge = 0.0;
+        particle_weight = 1.0;
 
-  biasRndm = new G4SPSRandomGenerator();
-  posGenerator = new G4SPSPosDistribution();
-  posGenerator->SetBiasRndm(biasRndm);
-  angGenerator = new G4SPSAngDistribution();
-  angGenerator->SetPosDistribution(posGenerator);
-  angGenerator->SetBiasRndm(biasRndm);
-  eneGenerator = new G4SPSEneDistribution();
-  eneGenerator->SetBiasRndm(biasRndm);
+        biasRndm = new G4SPSRandomGenerator();
+        posGenerator = new G4SPSPosDistribution();
+        posGenerator->SetBiasRndm(biasRndm);
+        angGenerator = new G4SPSAngDistribution();
+        angGenerator->SetPosDistribution(posGenerator);
+        angGenerator->SetBiasRndm(biasRndm);
+        eneGenerator = new G4SPSEneDistribution();
+        eneGenerator->SetBiasRndm(biasRndm);
 
-  // verbosity
-  verbosityLevel = 0;
+        // verbosity
+        verbosityLevel = 0;
 
 }
 
-G4SingleParticleSource::~G4SingleParticleSource()
-{}
-
-void G4SingleParticleSource::SetVerbosity(int vL)
-{
-  verbosityLevel = vL;
-  posGenerator->SetVerbosity(vL);
-  angGenerator->SetVerbosity(vL);
-  eneGenerator->SetVerbosity(vL);
-  G4cout << "Verbosity Set to: " << verbosityLevel << G4endl;
+G4SingleParticleSource::~G4SingleParticleSource() {
+        delete biasRndm;
+        delete posGenerator;
+        delete angGenerator;
+        delete eneGenerator;
 }
 
-void G4SingleParticleSource::SetParticleDefinition
-  (G4ParticleDefinition* aParticleDefinition)
-{
-  particle_definition = aParticleDefinition;
-  particle_charge = particle_definition->GetPDGCharge();
+void G4SingleParticleSource::SetVerbosity(int vL) {
+        verbosityLevel = vL;
+        posGenerator->SetVerbosity(vL);
+        angGenerator->SetVerbosity(vL);
+        eneGenerator->SetVerbosity(vL);
+        G4cout << "Verbosity Set to: " << verbosityLevel << G4endl;
 }
 
-void G4SingleParticleSource::GeneratePrimaryVertex(G4Event *evt)
-{
-  if(particle_definition==NULL) return;
-
-  if(verbosityLevel > 1)
-    G4cout << " NumberOfParticlesToBeGenerated: "<<NumberOfParticlesToBeGenerated << G4endl;
-
-  // Position stuff
-  particle_position = posGenerator->GenerateOne();
-
-  // create a new vertex
-  G4PrimaryVertex* vertex = new G4PrimaryVertex(particle_position,particle_time);
-
-  for( G4int i=0; i<NumberOfParticlesToBeGenerated; i++ ) {
-    // Angular stuff
-    particle_momentum_direction = angGenerator->GenerateOne();
-    // Energy stuff
-    particle_energy = eneGenerator->GenerateOne(particle_definition);
-    
-    if(verbosityLevel >= 2)
-      G4cout << "Creating primaries and assigning to vertex" << G4endl;
-    // create new primaries and set them to the vertex
-    G4double mass =  particle_definition->GetPDGMass();
-    G4double energy = particle_energy + mass;
-    G4double pmom = std::sqrt(energy*energy-mass*mass);
-    G4double px = pmom*particle_momentum_direction.x();
-    G4double py = pmom*particle_momentum_direction.y();
-    G4double pz = pmom*particle_momentum_direction.z();
-
-    if(verbosityLevel > 1){
-      G4cout << "Particle name: "<<particle_definition->GetParticleName() << G4endl; 
-      G4cout << "       Energy: "<<particle_energy << G4endl;
-      G4cout << "     Position: "<<particle_position<< G4endl; 
-      G4cout << "    Direction: "<<particle_momentum_direction << G4endl;
-    }
-    G4PrimaryParticle* particle =
-      new G4PrimaryParticle(particle_definition,px,py,pz);
-    particle->SetMass( mass );
-    particle->SetCharge( particle_charge );
-    particle->SetPolarization(particle_polarization.x(),
-                              particle_polarization.y(),
-                              particle_polarization.z());
-    // Set bweight equal to the multiple of all non-zero weights
-    particle_weight = biasRndm->GetBiasWeight();
-    // pass it to primary particle
-    particle->SetWeight(particle_weight);
-
-    vertex->SetPrimary( particle );
-      
-  }
-  // now pass the weight to the primary vertex. CANNOT be used here!
-  //  vertex->SetWeight(particle_weight);
-  evt->AddPrimaryVertex( vertex );
-  if(verbosityLevel > 1)
-    G4cout << " Primary Vetex generated !"<< G4endl;   
+void G4SingleParticleSource::SetParticleDefinition(
+                G4ParticleDefinition* aParticleDefinition) {
+        particle_definition = aParticleDefinition;
+        particle_charge = particle_definition->GetPDGCharge();
 }
 
+void G4SingleParticleSource::GeneratePrimaryVertex(G4Event *evt) {
+        if (particle_definition == NULL)
+                return;
 
+        if (verbosityLevel > 1)
+                G4cout << " NumberOfParticlesToBeGenerated: "
+                                << NumberOfParticlesToBeGenerated << G4endl;
 
+        // Position stuff
+        particle_position = posGenerator->GenerateOne();
 
+        // create a new vertex
+        G4PrimaryVertex* vertex = new G4PrimaryVertex(particle_position,
+                        particle_time);
 
+        for (G4int i = 0; i < NumberOfParticlesToBeGenerated; i++) {
+                // Angular stuff
+                particle_momentum_direction = angGenerator->GenerateOne();
+                // Energy stuff
+                particle_energy = eneGenerator->GenerateOne(particle_definition);
 
+                if (verbosityLevel >= 2)
+                        G4cout << "Creating primaries and assigning to vertex" << G4endl;
+                // create new primaries and set them to the vertex
+                G4double mass = particle_definition->GetPDGMass();
+                G4double energy = particle_energy + mass;
+                G4double pmom = std::sqrt(energy * energy - mass * mass);
+                G4double px = pmom * particle_momentum_direction.x();
+                G4double py = pmom * particle_momentum_direction.y();
+                G4double pz = pmom * particle_momentum_direction.z();
 
+                if (verbosityLevel > 1) {
+                        G4cout << "Particle name: "
+                                        << particle_definition->GetParticleName() << G4endl;
+                        G4cout << "       Energy: " << particle_energy << G4endl;
+                        G4cout << "     Position: " << particle_position << G4endl;
+                        G4cout << "    Direction: " << particle_momentum_direction
+                                        << G4endl;
+                }
+                G4PrimaryParticle* particle = new G4PrimaryParticle(
+                                particle_definition, px, py, pz);
+                particle->SetMass(mass);
+                particle->SetCharge(particle_charge);
+                particle->SetPolarization(particle_polarization.x(),
+                                particle_polarization.y(), particle_polarization.z());
+                // Set bweight equal to the multiple of all non-zero weights
+                particle_weight = eneGenerator->GetWeight()*biasRndm->GetBiasWeight();
+                // pass it to primary particle
+                particle->SetWeight(particle_weight);
 
+                vertex->SetPrimary(particle);
 
+        }
+        // now pass the weight to the primary vertex. CANNOT be used here!
+        //  vertex->SetWeight(particle_weight);
+        evt->AddPrimaryVertex(vertex);
+        if (verbosityLevel > 1)
+                G4cout << " Primary Vetex generated !" << G4endl;
+}
+