G4NeutronHPMadlandNixSpectrum.hh

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 #ifndef G4NeutronHPMadlandNixSpectrum_h
00030 #define G4NeutronHPMadlandNixSpectrum_h 1
00031 
00032 #include <fstream>
00033 #include <cmath>
00034 #include <CLHEP/Units/PhysicalConstants.h>
00035 
00036 #include "globals.hh"
00037 #include "G4ios.hh"
00038 #include "Randomize.hh"
00039 #include "G4NeutronHPVector.hh"
00040 #include "G4VNeutronHPEDis.hh"
00041 
00042 //     #include <nag.h> @
00043 //     #include <nags.h> @
00044 
00045 
00046 // we will need a List of these .... one per term.
00047 
00048 class G4NeutronHPMadlandNixSpectrum : public G4VNeutronHPEDis
00049 {
00050   public:
00051   G4NeutronHPMadlandNixSpectrum()
00052   {
00053     expm1 = std::exp(-1.);
00054   }
00055   ~G4NeutronHPMadlandNixSpectrum()
00056   {
00057   }
00058   
00059   inline void Init(std::ifstream & aDataFile)
00060   {
00061     theFractionalProb.Init(aDataFile);
00062     aDataFile>> theAvarageKineticPerNucleonForLightFragments;
00063     theAvarageKineticPerNucleonForLightFragments*=CLHEP::eV;
00064     aDataFile>> theAvarageKineticPerNucleonForHeavyFragments;
00065     theAvarageKineticPerNucleonForHeavyFragments*=CLHEP::eV;
00066     theMaxTemp.Init(aDataFile);
00067   }
00068   
00069   inline G4double GetFractionalProbability(G4double anEnergy)
00070   {
00071     return theFractionalProb.GetY(anEnergy);
00072   }
00073   
00074   G4double Sample(G4double anEnergy);
00075     
00076   private:
00077   
00078   G4double Madland(G4double aSecEnergy, G4double tm);
00079   
00080   inline G4double FissionIntegral(G4double tm, G4double anEnergy)
00081   {
00082     return 0.5*(  GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForLightFragments) 
00083                        +GIntegral(tm, anEnergy, theAvarageKineticPerNucleonForHeavyFragments) );
00084   }
00085   
00086   G4double GIntegral(G4double tm, G4double anEnergy, G4double aMean);
00087   
00088   inline G4double Gamma05(G4double aValue)
00089   {
00090     G4double result;
00091     // gamma(1.2,x*X) = std::sqrt(CLHEP::pi)*Erf(x)
00092     G4double x = std::sqrt(aValue);
00093     G4double t = 1./(1+0.47047*x);
00094     result = 1- (0.3480242*t - 0.0958798*t*t + 0.7478556*t*t*t)*std::exp(-aValue); // @ check
00095     result *= std::sqrt(CLHEP::pi);
00096     return result;
00097   }
00098   
00099   inline G4double Gamma15(G4double aValue)
00100   {
00101     G4double result;
00102     // gamma(a+1, x) = a*gamma(a,x)-x**a*std::exp(-x)
00103     result = 0.5*Gamma05(aValue) - std::sqrt(aValue)*std::exp(-aValue); // @ check
00104     return result;
00105   }
00106   
00107   inline G4double Gamma25(G4double aValue)
00108   {
00109     G4double result;
00110     result = 1.5*Gamma15(aValue) - std::pow(aValue,1.5)*std::exp(aValue); // @ check
00111     return result;
00112   }
00113   
00114   inline G4double E1(G4double aValue)
00115   {
00116   // good only for rather low aValue @@@ replace by the corresponding NAG function for the
00117   // exponential integral. (<5 seems ok.
00118     G4double gamma = 0.577216;
00119     G4double precision = 0.000001;
00120     G4double result =-gamma - std::log(aValue);
00121     G4double term = -aValue;
00122     //110527TKDB  Unnessary codes, Detected by gcc4.6 compiler 
00123     //G4double last;
00124     G4int count = 1;
00125     result -= term;
00126     for(;;)
00127     {
00128       count++;
00129       //110527TKDB  Unnessary codes, Detected by gcc4.6 compiler 
00130       //last = result;
00131       term = -term*aValue*(count-1)/(count*count);
00132       result -=term;
00133       if(std::fabs(term)/std::fabs(result)<precision) break;
00134     }
00135 //    NagError *fail; @
00136 //    result = nag_exp_integral(aValue, fail); @
00137     return result;
00138   }
00139   
00140   private:
00141   
00142   G4double expm1;
00143   
00144   private: 
00145   
00146   G4NeutronHPVector theFractionalProb;
00147   
00148   G4double theAvarageKineticPerNucleonForLightFragments;
00149   G4double theAvarageKineticPerNucleonForHeavyFragments;
00150 
00151   G4NeutronHPVector theMaxTemp;
00152   
00153 };
00154 
00155 #endif

Generated on Mon May 27 17:49:02 2013 for Geant4 by  doxygen 1.4.7