00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include "globals.hh"
00046 #include "G4ios.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "G4XnpTotalLowE.hh"
00049 #include "G4KineticTrack.hh"
00050 #include "G4ParticleDefinition.hh"
00051 #include "G4PhysicsVector.hh"
00052 #include "G4PhysicsLnVector.hh"
00053 #include "G4ParticleDefinition.hh"
00054 #include "G4Proton.hh"
00055 #include "G4Neutron.hh"
00056
00057 const G4double G4XnpTotalLowE::_lowLimit = 0.;
00058 const G4double G4XnpTotalLowE::_highLimit = 3.*GeV;
00059
00060
00061
00062 const G4double G4XnpTotalLowE::_eMinTable = 1.8964808;
00063 const G4double G4XnpTotalLowE::_eStepLog = 0.01;
00064
00065
00066
00067 const G4int G4XnpTotalLowE::_tableSize = 101;
00068 const G4double G4XnpTotalLowE::_sigmaTable[101] =
00069 {
00070 1500.0,
00071 248.20, 93.38, 55.26, 44.50, 41.33, 38.48, 37.20, 35.98,
00072 35.02, 34.47, 34.37, 34.67, 35.23, 35.97, 36.75, 37.37,
00073 37.77, 38.03, 38.40, 38.83, 39.26, 39.67, 40.06, 40.45,
00074 40.79, 41.06, 41.31, 41.52, 41.70, 41.81, 41.87, 41.98,
00075 42.12, 42.29, 42.55, 42.82, 43.01, 43.12, 43.16, 43.14,
00076 43.06, 42.95, 42.81, 42.67, 42.54, 42.45, 42.38, 42.33,
00077 42.30, 42.29, 42.28, 42.26, 42.24, 42.21, 42.17, 42.14,
00078 42.10, 42.07, 42.06, 42.05, 42.04, 42.03, 42.02, 42.00,
00079 41.97, 41.94, 41.89, 41.84, 41.79, 41.73, 41.67, 41.61,
00080 41.55, 41.49, 41.44, 41.38, 41.34, 41.31, 41.29, 41.28,
00081 41.27, 41.28, 41.30, 41.33, 41.36, 41.40, 41.44, 41.49,
00082 41.50, 41.51, 41.51, 41.51, 41.52, 41.51, 41.51, 41.50,
00083 41.50, 41.49, 41.47, 41.46
00084 };
00085
00086
00087 G4XnpTotalLowE::G4XnpTotalLowE()
00088 {
00089
00090
00091 _eMin = _eMinTable * GeV;
00092 _eMin = std::exp(std::log(_eMinTable)-_eStepLog)*GeV;
00093 _eMax = std::exp(std::log(_eMinTable) + _tableSize * _eStepLog) * GeV;
00094
00095
00096
00097 if (_eMin < _lowLimit)
00098 throw G4HadronicException(__FILE__, __LINE__, "G4XnpTotalLowE::G4XnpTotalLowE - Low energy limit not valid");
00099
00100 if (_highLimit > _eMax)
00101 throw G4HadronicException(__FILE__, __LINE__, "G4XnpTotalLowE::G4XnpTotalLowE - High energy limit not valid");
00102
00103 _sigma = new G4PhysicsLnVector(_eMin,_eMax,_tableSize);
00104 G4int i;
00105 for (i=0; i<_tableSize; i++)
00106 {
00107 G4double value = _sigmaTable[i] * millibarn;
00108 _sigma->PutValue(i,value);
00109 }
00110 }
00111
00112
00113 G4XnpTotalLowE::~G4XnpTotalLowE()
00114 {
00115 delete _sigma;
00116 }
00117
00118
00119 G4bool G4XnpTotalLowE::operator==(const G4XnpTotalLowE &right) const
00120 {
00121 return (this == (G4XnpTotalLowE *) &right);
00122 }
00123
00124
00125 G4bool G4XnpTotalLowE::operator!=(const G4XnpTotalLowE &right) const
00126 {
00127 return (this != (G4XnpTotalLowE *) &right);
00128 }
00129
00130
00131 G4double G4XnpTotalLowE::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const
00132 {
00133 G4double sigma = 0.;
00134 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
00135 G4bool dummy = false;
00136
00137 G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
00138 G4ParticleDefinition* neutron = G4Neutron::NeutronDefinition();
00139
00140 G4ParticleDefinition* def1 = trk1.GetDefinition();
00141 G4ParticleDefinition* def2 = trk2.GetDefinition();
00142 if ( (def1 == proton && def2 == neutron) ||
00143 (def1 == neutron && def2 == proton) )
00144 {
00145 if (sqrtS >= _eMin && sqrtS <= _eMax)
00146 {
00147 sigma = _sigma->GetValue(sqrtS,dummy);
00148 } else if ( sqrtS < _eMin )
00149 {
00150 sigma = _sigma->GetValue(_eMin,dummy);
00151 }
00152 }
00153
00154 return sigma;
00155 }
00156
00157 void G4XnpTotalLowE::Print() const
00158 {
00159
00160
00161 G4cout << Name() << "Cross-section table: " << G4endl;
00162 G4bool dummy = false;
00163 G4int i;
00164
00165 for (i=0; i<_tableSize; i++)
00166 {
00167 G4double e = _sigma->GetLowEdgeEnergy(i) / GeV;
00168 G4double sigma = _sigma->GetValue(e,dummy) / millibarn;
00169 G4cout << i << ") e = " << e << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
00170 }
00171
00172 G4VCrossSectionSource::Print();
00173 }
00174
00175
00176 G4String G4XnpTotalLowE::Name() const
00177 {
00178 G4String name("NNTotalLowE");
00179 return name;
00180 }
00181
00182
00183 G4bool G4XnpTotalLowE::IsValid(G4double e) const
00184 {
00185 G4bool answer = InLimits(e,_lowLimit,_highLimit);
00186
00187 return answer;
00188 }
00189
00190