include/Geant4/G4ParticleHPMadlandNixSpectrum.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
0028 //
0029 #ifndef G4ParticleHPMadlandNixSpectrum_h
0030 #define G4ParticleHPMadlandNixSpectrum_h 1
0031
0032 #include "G4Exp.hh"
0033 #include "G4Log.hh"
0034 #include "G4ParticleHPVector.hh"
0035 #include "G4Pow.hh"
0036 #include "G4VParticleHPEDis.hh"
0037 #include "G4ios.hh"
0038 #include "Randomize.hh"
0039 #include "globals.hh"
0040
0041 #include <CLHEP/Units/PhysicalConstants.h>
0042
0043 #include <cmath>
0044 #include <fstream>
0045
0046 //     #include <nag.h> @
0047 //     #include <nags.h> @
0048
0049 // we will need a List of these .... one per term.
0050
0051 class G4ParticleHPMadlandNixSpectrum : public G4VParticleHPEDis
0052 {
0053   public:
0054     G4ParticleHPMadlandNixSpectrum()
0055     {
0056       expm1 = G4Exp(-1.);
0057       theAvarageKineticPerNucleonForLightFragments = 0.0;
0058       theAvarageKineticPerNucleonForHeavyFragments = 0.0;
0059     }
0060     ~G4ParticleHPMadlandNixSpectrum() override = default;
0061
0062     inline void Init(std::istream& aDataFile) override
0063     {
0064       theFractionalProb.Init(aDataFile);
0065       aDataFile >> theAvarageKineticPerNucleonForLightFragments;
0066       theAvarageKineticPerNucleonForLightFragments *= CLHEP::eV;
0067       aDataFile >> theAvarageKineticPerNucleonForHeavyFragments;
0068       theAvarageKineticPerNucleonForHeavyFragments *= CLHEP::eV;
0069       theMaxTemp.Init(aDataFile);
0070     }
0071
0072     inline G4double GetFractionalProbability(G4double anEnergy) override
0073     {
0074       return theFractionalProb.GetY(anEnergy);
0075     }
0076
0077     G4double Sample(G4double anEnergy) override;
0078
0079   private:
0080     G4double Madland(G4double aSecEnergy, G4double tm);
0081
0082     inline G4double FissionIntegral(G4double tm, G4double anEnergy)
0083     {
0084       return 0.5
0085              * (GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForLightFragments)
0086                 + GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForHeavyFragments));
0087     }
0088
0089     G4double GIntegral(G4double tm, G4double anEnergy, G4double aMean);
0090
0091     inline G4double Gamma05(G4double aValue)
0092     {
0093       G4double result;
0094       // gamma(1.2,x*X) = std::sqrt(CLHEP::pi)*Erf(x)
0095       G4double x = std::sqrt(aValue);
0096       G4double t = 1. / (1 + 0.47047 * x);
0097       result =
0098         1
0099         - (0.3480242 * t - 0.0958798 * t * t + 0.7478556 * t * t * t) * G4Exp(-aValue);  // @ check
0100       result *= std::sqrt(CLHEP::pi);
0101       return result;
0102     }
0103
0104     inline G4double Gamma15(G4double aValue)
0105     {
0106       G4double result;
0107       // gamma(a+1, x) = a*gamma(a,x)-x**a*std::exp(-x)
0108       result = 0.5 * Gamma05(aValue) - std::sqrt(aValue) * G4Exp(-aValue);  // @ check
0109       return result;
0110     }
0111
0112     inline G4double Gamma25(G4double aValue)
0113     {
0114       G4double result;
0115       result =
0116         1.5 * Gamma15(aValue) - G4Pow::GetInstance()->powA(aValue, 1.5) * G4Exp(aValue);  // @ check
0117       return result;
0118     }
0119
0120     inline G4double E1(G4double aValue)
0121     {
0122       // good only for rather low aValue @@@ replace by the corresponding NAG function for the
0123       // exponential integral. (<5 seems ok.
0124       G4double gamma = 0.577216;
0125       G4double precision = 0.000001;
0126       G4double result = -gamma - G4Log(aValue);
0127       G4double term = -aValue;
0128       // 110527TKDB  Unnessary codes, Detected by gcc4.6 compiler
0129       // G4double last;
0130       G4int count = 1;
0131       result -= term;
0132       for (;;) {
0133         count++;
0134         // 110527TKDB  Unnessary codes, Detected by gcc4.6 compiler
0135         // last = result;
0136         term = -term * aValue * (count - 1) / (count * count);
0137         result -= term;
0138         if (std::fabs(term) / std::fabs(result) < precision) break;
0139       }
0140       //    NagError *fail; @
0141       //    result = nag_exp_integral(aValue, fail); @
0142       return result;
0143     }
0144
0145   private:
0146     G4double expm1;
0147
0148   private:
0149     G4ParticleHPVector theFractionalProb;
0150
0151     G4double theAvarageKineticPerNucleonForLightFragments;
0152     G4double theAvarageKineticPerNucleonForHeavyFragments;
0153
0154     G4ParticleHPVector theMaxTemp;
0155 };
0156
0157 #endif