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 "G4XnpElasticLowE.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 G4XnpElasticLowE::_lowLimit = 0.;
00058 const G4double G4XnpElasticLowE::_highLimit = 3.*GeV;
00059
00060
00061
00062 const G4double G4XnpElasticLowE::_eMinTable = 1.8964808;
00063 const G4double G4XnpElasticLowE::_eStepLog = 0.01;
00064
00065
00066
00067 const G4int G4XnpElasticLowE::_tableSize = 101;
00068 const G4double G4XnpElasticLowE::_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, 32.48, 30.76, 29.46, 28.53, 27.84, 27.20,
00073 26.53, 25.95, 25.59, 25.46, 25.00, 24.49, 24.08, 23.86,
00074 23.17, 22.70, 21.88, 21.48, 20.22, 19.75, 18.97, 18.39,
00075 17.98, 17.63, 17.21, 16.72, 16.68, 16.58, 16.42, 16.22,
00076 15.98, 15.71, 15.42, 15.14, 14.87, 14.65, 14.44, 14.26,
00077 14.10, 13.95, 13.80, 13.64, 13.47, 13.29, 13.09, 12.89,
00078 12.68, 12.47, 12.27, 12.06, 11.84, 11.76, 11.69, 11.60,
00079 11.50, 11.41, 11.29, 11.17, 11.06, 10.93, 10.81, 10.68,
00080 10.56, 10.44, 10.33, 10.21, 10.12, 10.03, 9.96, 9.89,
00081 9.83, 9.80, 9.77, 9.75, 9.74, 9.74, 9.74, 9.76,
00082 9.73, 9.70, 9.68, 9.65, 9.63, 9.60, 9.57, 9.55,
00083 9.52, 9.49, 9.46, 9.43
00084 };
00085
00086
00087 G4XnpElasticLowE::G4XnpElasticLowE()
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__, "G4XnpElasticLowE::G4XnpElasticLowE - Low energy limit not valid");
00099
00100 if (_highLimit > _eMax)
00101 throw G4HadronicException(__FILE__, __LINE__, "G4XnpElasticLowE::G4XnpElasticLowE - 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 G4XnpElasticLowE::~G4XnpElasticLowE()
00114 {
00115 delete _sigma;
00116 }
00117
00118
00119 G4bool G4XnpElasticLowE::operator==(const G4XnpElasticLowE &right) const
00120 {
00121 return (this == (G4XnpElasticLowE *) &right);
00122 }
00123
00124
00125 G4bool G4XnpElasticLowE::operator!=(const G4XnpElasticLowE &right) const
00126 {
00127 return (this != (G4XnpElasticLowE *) &right);
00128 }
00129
00130
00131 G4double G4XnpElasticLowE::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 G4XnpElasticLowE::Print() const
00158 {
00159
00160 G4cout << Name() << "Cross-section table: " << G4endl;
00161 G4bool dummy = false;
00162 G4int i;
00163
00164 for (i=0; i<_tableSize; i++)
00165 {
00166 G4double e = _sigma->GetLowEdgeEnergy(i) / GeV;
00167 G4double sigma = _sigma->GetValue(e,dummy) / millibarn;
00168 G4cout << i << ") e = " << e << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
00169 }
00170
00171 G4VCrossSectionSource::Print();
00172 }
00173
00174
00175 G4String G4XnpElasticLowE::Name() const
00176 {
00177 G4String name("npElasticLowE");
00178 return name;
00179 }
00180
00181
00182 G4bool G4XnpElasticLowE::IsValid(G4double e) const
00183 {
00184 G4bool answer = InLimits(e,_lowLimit,_highLimit);
00185
00186 return answer;
00187 }
00188
00189