| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 #include "globals.hh" 26 #include "globals.hh"
27 #include "G4ios.hh" 27 #include "G4ios.hh"
28 #include "G4Exp.hh" <<
29 #include "G4Log.hh" <<
30 #include "G4SystemOfUnits.hh" 28 #include "G4SystemOfUnits.hh"
31 #include "G4XNNElasticLowE.hh" 29 #include "G4XNNElasticLowE.hh"
32 #include "G4KineticTrack.hh" 30 #include "G4KineticTrack.hh"
33 #include "G4ParticleDefinition.hh" 31 #include "G4ParticleDefinition.hh"
34 #include "G4PhysicsVector.hh" 32 #include "G4PhysicsVector.hh"
35 #include "G4PhysicsLogVector.hh" << 33 #include "G4PhysicsLnVector.hh"
36 #include "G4Proton.hh" 34 #include "G4Proton.hh"
37 #include "G4Neutron.hh" 35 #include "G4Neutron.hh"
38 36
39 const G4double G4XNNElasticLowE::_lowLimit = 0 37 const G4double G4XNNElasticLowE::_lowLimit = 0.;
40 const G4double G4XNNElasticLowE::_highLimit = 38 const G4double G4XNNElasticLowE::_highLimit = 3.*GeV;
41 39
42 // Low energy limit of the cross-section table 40 // Low energy limit of the cross-section table (in GeV)
43 // Units are assigned when filling the Physics << 41 // Units are assigned while filling the PhysicsVector
44 const G4double G4XNNElasticLowE::_eMinTable = 42 const G4double G4XNNElasticLowE::_eMinTable = 1.8964808;
45 const G4double G4XNNElasticLowE::_eStepLog = 0 43 const G4double G4XNNElasticLowE::_eStepLog = 0.01;
46 44
47 // Cross-sections in mb 45 // Cross-sections in mb
48 // Units are assigned when filling the Physics << 46 // Units are assigned while filling the PhysicsVector
49 47
50 const G4int G4XNNElasticLowE::tableSize = 101; 48 const G4int G4XNNElasticLowE::tableSize = 101;
51 49
52 const G4double G4XNNElasticLowE::ppTable[101] 50 const G4double G4XNNElasticLowE::ppTable[101] =
53 { 51 {
54 60.00, //was 0. 52 60.00, //was 0.
55 33.48, 26.76, 25.26, 24.55, 23.94, 23.77, 23 53 33.48, 26.76, 25.26, 24.55, 23.94, 23.77, 23.72, 23.98,
56 25.48, 27.52, 27.72, 27.21, 25.80, 26.00, 24 54 25.48, 27.52, 27.72, 27.21, 25.80, 26.00, 24.32, 23.81,
57 24.37, 24.36, 23.13, 22.43, 21.71, 21.01, 20 55 24.37, 24.36, 23.13, 22.43, 21.71, 21.01, 20.83, 20.74,
58 20.25, 20.10, 20.59, 20.04, 20.83, 20.84, 21 56 20.25, 20.10, 20.59, 20.04, 20.83, 20.84, 21.07, 20.83,
59 20.79, 21.88, 21.15, 20.92, 19.00, 18.60, 17 57 20.79, 21.88, 21.15, 20.92, 19.00, 18.60, 17.30, 17.00,
60 16.70, 16.50, 16.20, 15.80, 15.57, 15.20, 15 58 16.70, 16.50, 16.20, 15.80, 15.57, 15.20, 15.00, 14.60,
61 14.20, 14.00, 13.80, 13.60, 13.40, 13.20, 13 59 14.20, 14.00, 13.80, 13.60, 13.40, 13.20, 13.00, 12.85,
62 12.70, 12.60, 12.50, 12.40, 12.30, 12.20, 12 60 12.70, 12.60, 12.50, 12.40, 12.30, 12.20, 12.10, 12.00,
63 11.90, 11.80, 11.75, 11.70, 11.64, 11.53, 11 61 11.90, 11.80, 11.75, 11.70, 11.64, 11.53, 11.41, 11.31,
64 11.22, 11.13, 11.05, 10.97, 10.89, 10.82, 10 62 11.22, 11.13, 11.05, 10.97, 10.89, 10.82, 10.75, 10.68,
65 10.61, 10.54, 10.48, 10.41, 10.35, 10.28, 10 63 10.61, 10.54, 10.48, 10.41, 10.35, 10.28, 10.22, 10.16,
66 10.13, 10.10, 10.08, 10.05, 10.02, 9.99, 9 64 10.13, 10.10, 10.08, 10.05, 10.02, 9.99, 9.96, 9.93,
67 9.90, 9.87, 9.84, 9.80 65 9.90, 9.87, 9.84, 9.80
68 }; 66 };
69 67
70 const G4double G4XNNElasticLowE::npTable[101] 68 const G4double G4XNNElasticLowE::npTable[101] =
71 { 69 {
72 1500.00, // was 0. 70 1500.00, // was 0.
73 248.20, 93.38, 55.26, 44.50, 41.33, 38.48, 3 71 248.20, 93.38, 55.26, 44.50, 41.33, 38.48, 37.20, 35.98,
74 35.02, 34.47, 32.48, 30.76, 29.46, 28.53, 27 72 35.02, 34.47, 32.48, 30.76, 29.46, 28.53, 27.84, 27.20,
75 26.53, 25.95, 25.59, 25.46, 25.00, 24.49, 24 73 26.53, 25.95, 25.59, 25.46, 25.00, 24.49, 24.08, 23.86,
76 23.17, 22.70, 21.88, 21.48, 20.22, 19.75, 18 74 23.17, 22.70, 21.88, 21.48, 20.22, 19.75, 18.97, 18.39,
77 17.98, 17.63, 17.21, 16.72, 16.68, 16.58, 16 75 17.98, 17.63, 17.21, 16.72, 16.68, 16.58, 16.42, 16.22,
78 15.98, 15.71, 15.42, 15.14, 14.87, 14.65, 14 76 15.98, 15.71, 15.42, 15.14, 14.87, 14.65, 14.44, 14.26,
79 14.10, 13.95, 13.80, 13.64, 13.47, 13.29, 13 77 14.10, 13.95, 13.80, 13.64, 13.47, 13.29, 13.09, 12.89,
80 12.68, 12.47, 12.27, 12.06, 11.84, 11.76, 11 78 12.68, 12.47, 12.27, 12.06, 11.84, 11.76, 11.69, 11.60,
81 11.50, 11.41, 11.29, 11.17, 11.06, 10.93, 10 79 11.50, 11.41, 11.29, 11.17, 11.06, 10.93, 10.81, 10.68,
82 10.56, 10.44, 10.33, 10.21, 10.12, 10.03, 9 80 10.56, 10.44, 10.33, 10.21, 10.12, 10.03, 9.96, 9.89,
83 9.83, 9.80, 9.77, 9.75, 9.74, 9.74, 9. 81 9.83, 9.80, 9.77, 9.75, 9.74, 9.74, 9.74, 9.76,
84 9.73, 9.70, 9.68, 9.65, 9.63, 9.60, 9. 82 9.73, 9.70, 9.68, 9.65, 9.63, 9.60, 9.57, 9.55,
85 9.52, 9.49, 9.46, 9.43 83 9.52, 9.49, 9.46, 9.43
86 }; 84 };
87 85
88 86
89 G4XNNElasticLowE::G4XNNElasticLowE() 87 G4XNNElasticLowE::G4XNNElasticLowE()
90 { 88 {
91 // Cross-sections are available in the range 89 // Cross-sections are available in the range (_eMin,_eMax)
92 90
93 _eMin = _eMinTable * GeV; 91 _eMin = _eMinTable * GeV;
94 _eMax = G4Exp(G4Log(_eMinTable) + tableSize << 92 _eMax = std::exp(std::log(_eMinTable) + tableSize * _eStepLog) * GeV;
95 if (_eMin < _lowLimit) 93 if (_eMin < _lowLimit)
96 throw G4HadronicException(__FILE__, __LINE 94 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");
97 if (_highLimit > _eMax) 95 if (_highLimit > _eMax)
98 throw G4HadronicException(__FILE__, __LINE 96 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - High energy limit not valid");
99 G4PhysicsVector* pp = new G4PhysicsLogVector << 97 G4PhysicsVector* pp = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
100 98
101 _eMin = G4Exp(G4Log(_eMinTable)-_eStepLog)*G << 99 _eMin = std::exp(std::log(_eMinTable)-_eStepLog)*GeV;
102 if (_eMin < _lowLimit) 100 if (_eMin < _lowLimit)
103 throw G4HadronicException(__FILE__, __LINE 101 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");
104 G4PhysicsVector* np = new G4PhysicsLogVector << 102 G4PhysicsVector* np = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
105 103
106 G4int i; 104 G4int i;
107 for (i=0; i<tableSize; i++) 105 for (i=0; i<tableSize; i++)
108 { 106 {
109 G4double value = ppTable[i] * millibarn; 107 G4double value = ppTable[i] * millibarn;
110 pp->PutValue(i,value); 108 pp->PutValue(i,value);
111 value = npTable[i] * millibarn; 109 value = npTable[i] * millibarn;
112 np->PutValue(i,value); 110 np->PutValue(i,value);
113 } 111 }
114 xMap[G4Proton::ProtonDefinition()] = pp; 112 xMap[G4Proton::ProtonDefinition()] = pp;
115 xMap[G4Neutron::NeutronDefinition()] = np; 113 xMap[G4Neutron::NeutronDefinition()] = np;
116 } 114 }
117 115
118 116
119 G4XNNElasticLowE::~G4XNNElasticLowE() 117 G4XNNElasticLowE::~G4XNNElasticLowE()
120 { 118 {
121 delete xMap[G4Proton::ProtonDefinition()]; 119 delete xMap[G4Proton::ProtonDefinition()];
122 delete xMap[G4Neutron::NeutronDefinition()]; 120 delete xMap[G4Neutron::NeutronDefinition()];
123 } 121 }
124 122
125 123
126 G4bool G4XNNElasticLowE::operator==(const G4XN 124 G4bool G4XNNElasticLowE::operator==(const G4XNNElasticLowE &right) const
127 { 125 {
128 return (this == (G4XNNElasticLowE *) &right) 126 return (this == (G4XNNElasticLowE *) &right);
129 } 127 }
130 128
131 129
132 G4bool G4XNNElasticLowE::operator!=(const G4XN 130 G4bool G4XNNElasticLowE::operator!=(const G4XNNElasticLowE &right) const
133 { 131 {
134 132
135 return (this != (G4XNNElasticLowE *) &right) 133 return (this != (G4XNNElasticLowE *) &right);
136 } 134 }
137 135
138 136
139 G4double G4XNNElasticLowE::CrossSection(const 137 G4double G4XNNElasticLowE::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const
140 { 138 {
141 G4double sigma = 0.; 139 G4double sigma = 0.;
142 G4double sqrtS = (trk1.Get4Momentum() + trk2 140 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
143 G4bool dummy = false; 141 G4bool dummy = false;
144 142
145 const G4ParticleDefinition * key = FindKeyPa 143 const G4ParticleDefinition * key = FindKeyParticle(trk1,trk2);
146 144
147 typedef std::map <const G4ParticleDefinition 145 typedef std::map <const G4ParticleDefinition *, G4PhysicsVector*, std::less<const G4ParticleDefinition *> > StringPhysMap;
148 146
149 if (xMap.find(key)!= xMap.end()) 147 if (xMap.find(key)!= xMap.end())
150 { 148 {
151 149
152 StringPhysMap::const_iterator iter; 150 StringPhysMap::const_iterator iter;
153 for (iter = xMap.begin(); iter != xMap.e 151 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
154 { 152 {
155 const G4ParticleDefinition * str = (*iter) 153 const G4ParticleDefinition * str = (*iter).first;
156 if (str == key) 154 if (str == key)
157 { 155 {
158 G4PhysicsVector* physVector = (*iter). 156 G4PhysicsVector* physVector = (*iter).second;
159 // G4PhysicsVector* physVector = x 157 // G4PhysicsVector* physVector = xMap[key];
160 if (sqrtS >= _eMin && sqrtS <= _eMax) 158 if (sqrtS >= _eMin && sqrtS <= _eMax)
161 { 159 {
162 sigma = physVector->GetValue(sqrtS,dummy 160 sigma = physVector->GetValue(sqrtS,dummy);
163 } else if ( sqrtS < _eMin ) 161 } else if ( sqrtS < _eMin )
164 { 162 {
165 sigma = physVector->GetValue 163 sigma = physVector->GetValue(_eMin,dummy);
166 } 164 }
167 //G4cout << " sqrtS / sigma " << sqrtS/GeV 165 //G4cout << " sqrtS / sigma " << sqrtS/GeV << " / " <<
168 // sigma/millibarn << G4endl; 166 // sigma/millibarn << G4endl;
169 } 167 }
170 } 168 }
171 } 169 }
172 return sigma; 170 return sigma;
173 } 171 }
174 172
175 173
176 void G4XNNElasticLowE::Print() const 174 void G4XNNElasticLowE::Print() const
177 { 175 {
178 // Dump the pp cross-section table 176 // Dump the pp cross-section table
179 177
180 G4cout << Name() << ", pp cross-section: " < 178 G4cout << Name() << ", pp cross-section: " << G4endl;
181 179
182 G4bool dummy = false; 180 G4bool dummy = false;
183 G4int i; 181 G4int i;
184 const G4ParticleDefinition * key = G4Proton: 182 const G4ParticleDefinition * key = G4Proton::ProtonDefinition();
185 G4PhysicsVector* pp = 0; 183 G4PhysicsVector* pp = 0;
186 184
187 typedef std::map <const G4ParticleDefinition 185 typedef std::map <const G4ParticleDefinition *, G4PhysicsVector*, std::less<const G4ParticleDefinition *> > StringPhysMap;
188 StringPhysMap::const_iterator iter; 186 StringPhysMap::const_iterator iter;
189 187
190 for (iter = xMap.begin(); iter != xMap.end() 188 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
191 { 189 {
192 const G4ParticleDefinition * str = (*ite 190 const G4ParticleDefinition * str = (*iter).first;
193 if (str == key) 191 if (str == key)
194 { 192 {
195 pp = (*iter).second; 193 pp = (*iter).second;
196 } 194 }
197 } 195 }
198 196
199 if (pp != 0) 197 if (pp != 0)
200 { 198 {
201 for (i=0; i<tableSize; i++) 199 for (i=0; i<tableSize; i++)
202 { 200 {
203 G4double e = pp->GetLowEdgeEnergy(i); 201 G4double e = pp->GetLowEdgeEnergy(i);
204 G4double sigma = pp->GetValue(e,dummy) / m 202 G4double sigma = pp->GetValue(e,dummy) / millibarn;
205 G4cout << i << ") e = " << e / GeV << " Ge 203 G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
206 } 204 }
207 } 205 }
208 206
209 // Dump the np cross-section table 207 // Dump the np cross-section table
210 208
211 G4cout << Name() << ", np cross-section: " < 209 G4cout << Name() << ", np cross-section: " << G4endl;
212 210
213 key = G4Neutron::NeutronDefinition(); 211 key = G4Neutron::NeutronDefinition();
214 G4PhysicsVector* np = 0; 212 G4PhysicsVector* np = 0;
215 for (iter = xMap.begin(); iter != xMap.end() 213 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
216 { 214 {
217 const G4ParticleDefinition * str = (*ite 215 const G4ParticleDefinition * str = (*iter).first;
218 if (str == key) 216 if (str == key)
219 { 217 {
220 np = (*iter).second; 218 np = (*iter).second;
221 } 219 }
222 } 220 }
223 221
224 // G4PhysicsVector* np = xMap[G4Neutron::Ne 222 // G4PhysicsVector* np = xMap[G4Neutron::NeutronDefinition()->GetParticleName()];
225 223
226 if (np != 0) 224 if (np != 0)
227 { 225 {
228 for (i=0; i<tableSize; i++) 226 for (i=0; i<tableSize; i++)
229 { 227 {
230 G4double e = np->GetLowEdgeEnergy(i); 228 G4double e = np->GetLowEdgeEnergy(i);
231 G4double sigma = np->GetValue(e,dummy) / m 229 G4double sigma = np->GetValue(e,dummy) / millibarn;
232 G4cout << i << ") e = " << e / GeV << " Ge 230 G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
233 } 231 }
234 } 232 }
235 G4VCrossSectionSource::Print(); 233 G4VCrossSectionSource::Print();
236 } 234 }
237 235
238 236
239 G4String G4XNNElasticLowE::Name() const 237 G4String G4XNNElasticLowE::Name() const
240 { 238 {
241 G4String name("NNElasticLowE"); 239 G4String name("NNElasticLowE");
242 return name; 240 return name;
243 } 241 }
244 242
245 243
246 244
247 G4bool G4XNNElasticLowE::IsValid(G4double e) c 245 G4bool G4XNNElasticLowE::IsValid(G4double e) const
248 { 246 {
249 G4bool answer = InLimits(e,_lowLimit,_highLi 247 G4bool answer = InLimits(e,_lowLimit,_highLimit);
250 248
251 return answer; 249 return answer;
252 } 250 }
253 251
254 252
255 253