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 #include "globals.hh"
00027 #include "G4ios.hh"
00028 #include "G4SystemOfUnits.hh"
00029 #include "G4XNNElasticLowE.hh"
00030 #include "G4KineticTrack.hh"
00031 #include "G4ParticleDefinition.hh"
00032 #include "G4PhysicsVector.hh"
00033 #include "G4PhysicsLnVector.hh"
00034 #include "G4ParticleDefinition.hh"
00035 #include "G4Proton.hh"
00036 #include "G4Neutron.hh"
00037
00038 const G4double G4XNNElasticLowE::_lowLimit = 0.;
00039 const G4double G4XNNElasticLowE::_highLimit = 3.*GeV;
00040
00041
00042
00043 const G4double G4XNNElasticLowE::_eMinTable = 1.8964808;
00044 const G4double G4XNNElasticLowE::_eStepLog = 0.01;
00045
00046
00047
00048
00049 const G4int G4XNNElasticLowE::tableSize = 101;
00050
00051 const G4double G4XNNElasticLowE::ppTable[101] =
00052 {
00053 60.00,
00054 33.48, 26.76, 25.26, 24.55, 23.94, 23.77, 23.72, 23.98,
00055 25.48, 27.52, 27.72, 27.21, 25.80, 26.00, 24.32, 23.81,
00056 24.37, 24.36, 23.13, 22.43, 21.71, 21.01, 20.83, 20.74,
00057 20.25, 20.10, 20.59, 20.04, 20.83, 20.84, 21.07, 20.83,
00058 20.79, 21.88, 21.15, 20.92, 19.00, 18.60, 17.30, 17.00,
00059 16.70, 16.50, 16.20, 15.80, 15.57, 15.20, 15.00, 14.60,
00060 14.20, 14.00, 13.80, 13.60, 13.40, 13.20, 13.00, 12.85,
00061 12.70, 12.60, 12.50, 12.40, 12.30, 12.20, 12.10, 12.00,
00062 11.90, 11.80, 11.75, 11.70, 11.64, 11.53, 11.41, 11.31,
00063 11.22, 11.13, 11.05, 10.97, 10.89, 10.82, 10.75, 10.68,
00064 10.61, 10.54, 10.48, 10.41, 10.35, 10.28, 10.22, 10.16,
00065 10.13, 10.10, 10.08, 10.05, 10.02, 9.99, 9.96, 9.93,
00066 9.90, 9.87, 9.84, 9.80
00067 };
00068
00069 const G4double G4XNNElasticLowE::npTable[101] =
00070 {
00071 1500.00,
00072 248.20, 93.38, 55.26, 44.50, 41.33, 38.48, 37.20, 35.98,
00073 35.02, 34.47, 32.48, 30.76, 29.46, 28.53, 27.84, 27.20,
00074 26.53, 25.95, 25.59, 25.46, 25.00, 24.49, 24.08, 23.86,
00075 23.17, 22.70, 21.88, 21.48, 20.22, 19.75, 18.97, 18.39,
00076 17.98, 17.63, 17.21, 16.72, 16.68, 16.58, 16.42, 16.22,
00077 15.98, 15.71, 15.42, 15.14, 14.87, 14.65, 14.44, 14.26,
00078 14.10, 13.95, 13.80, 13.64, 13.47, 13.29, 13.09, 12.89,
00079 12.68, 12.47, 12.27, 12.06, 11.84, 11.76, 11.69, 11.60,
00080 11.50, 11.41, 11.29, 11.17, 11.06, 10.93, 10.81, 10.68,
00081 10.56, 10.44, 10.33, 10.21, 10.12, 10.03, 9.96, 9.89,
00082 9.83, 9.80, 9.77, 9.75, 9.74, 9.74, 9.74, 9.76,
00083 9.73, 9.70, 9.68, 9.65, 9.63, 9.60, 9.57, 9.55,
00084 9.52, 9.49, 9.46, 9.43
00085 };
00086
00087
00088 G4XNNElasticLowE::G4XNNElasticLowE()
00089 {
00090
00091
00092 _eMin = _eMinTable * GeV;
00093 _eMax = std::exp(std::log(_eMinTable) + tableSize * _eStepLog) * GeV;
00094 if (_eMin < _lowLimit)
00095 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");
00096 if (_highLimit > _eMax)
00097 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - High energy limit not valid");
00098 G4PhysicsVector* pp = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
00099
00100 _eMin = std::exp(std::log(_eMinTable)-_eStepLog)*GeV;
00101 if (_eMin < _lowLimit)
00102 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");
00103 G4PhysicsVector* np = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
00104
00105 G4int i;
00106 for (i=0; i<tableSize; i++)
00107 {
00108 G4double value = ppTable[i] * millibarn;
00109 pp->PutValue(i,value);
00110 value = npTable[i] * millibarn;
00111 np->PutValue(i,value);
00112 }
00113 xMap[G4Proton::ProtonDefinition()] = pp;
00114 xMap[G4Neutron::NeutronDefinition()] = np;
00115 }
00116
00117
00118 G4XNNElasticLowE::~G4XNNElasticLowE()
00119 {
00120 delete xMap[G4Proton::ProtonDefinition()];
00121 delete xMap[G4Neutron::NeutronDefinition()];
00122 }
00123
00124
00125 G4bool G4XNNElasticLowE::operator==(const G4XNNElasticLowE &right) const
00126 {
00127 return (this == (G4XNNElasticLowE *) &right);
00128 }
00129
00130
00131 G4bool G4XNNElasticLowE::operator!=(const G4XNNElasticLowE &right) const
00132 {
00133
00134 return (this != (G4XNNElasticLowE *) &right);
00135 }
00136
00137
00138 G4double G4XNNElasticLowE::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const
00139 {
00140 G4double sigma = 0.;
00141 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
00142 G4bool dummy = false;
00143
00144 G4ParticleDefinition * key = FindKeyParticle(trk1,trk2);
00145
00146 typedef std::map <G4ParticleDefinition *, G4PhysicsVector*, std::less<G4ParticleDefinition *> > StringPhysMap;
00147
00148 if (xMap.find(key)!= xMap.end())
00149 {
00150
00151 StringPhysMap::const_iterator iter;
00152 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
00153 {
00154 G4ParticleDefinition * str = (*iter).first;
00155 if (str == key)
00156 {
00157 G4PhysicsVector* physVector = (*iter).second;
00158
00159 if (sqrtS >= _eMin && sqrtS <= _eMax)
00160 {
00161 sigma = physVector->GetValue(sqrtS,dummy);
00162 } else if ( sqrtS < _eMin )
00163 {
00164 sigma = physVector->GetValue(_eMin,dummy);
00165 }
00166
00167
00168 }
00169 }
00170 }
00171 return sigma;
00172 }
00173
00174
00175 void G4XNNElasticLowE::Print() const
00176 {
00177
00178
00179 G4cout << Name() << ", pp cross-section: " << G4endl;
00180
00181 G4bool dummy = false;
00182 G4int i;
00183 G4ParticleDefinition * key = G4Proton::ProtonDefinition();
00184 G4PhysicsVector* pp = 0;
00185
00186 typedef std::map <G4ParticleDefinition *, G4PhysicsVector*, std::less<G4ParticleDefinition *> > StringPhysMap;
00187 StringPhysMap::const_iterator iter;
00188
00189 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
00190 {
00191 G4ParticleDefinition * str = (*iter).first;
00192 if (str == key)
00193 {
00194 pp = (*iter).second;
00195 }
00196 }
00197
00198 if (pp != 0)
00199 {
00200 for (i=0; i<tableSize; i++)
00201 {
00202 G4double e = pp->GetLowEdgeEnergy(i);
00203 G4double sigma = pp->GetValue(e,dummy) / millibarn;
00204 G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
00205 }
00206 }
00207
00208
00209
00210 G4cout << Name() << ", np cross-section: " << G4endl;
00211
00212 key = G4Neutron::NeutronDefinition();
00213 G4PhysicsVector* np = 0;
00214 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
00215 {
00216 G4ParticleDefinition * str = (*iter).first;
00217 if (str == key)
00218 {
00219 np = (*iter).second;
00220 }
00221 }
00222
00223
00224
00225 if (np != 0)
00226 {
00227 for (i=0; i<tableSize; i++)
00228 {
00229 G4double e = np->GetLowEdgeEnergy(i);
00230 G4double sigma = np->GetValue(e,dummy) / millibarn;
00231 G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
00232 }
00233 }
00234 G4VCrossSectionSource::Print();
00235 }
00236
00237
00238 G4String G4XNNElasticLowE::Name() const
00239 {
00240 G4String name("NNElasticLowE");
00241 return name;
00242 }
00243
00244
00245
00246 G4bool G4XNNElasticLowE::IsValid(G4double e) const
00247 {
00248 G4bool answer = InLimits(e,_lowLimit,_highLimit);
00249
00250 return answer;
00251 }
00252
00253