| 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" 28 #include "G4Exp.hh"
29 #include "G4Log.hh" 29 #include "G4Log.hh"
30 #include "G4SystemOfUnits.hh" 30 #include "G4SystemOfUnits.hh"
31 #include "G4XNNElasticLowE.hh" 31 #include "G4XNNElasticLowE.hh"
32 #include "G4KineticTrack.hh" 32 #include "G4KineticTrack.hh"
33 #include "G4ParticleDefinition.hh" 33 #include "G4ParticleDefinition.hh"
34 #include "G4PhysicsVector.hh" 34 #include "G4PhysicsVector.hh"
35 #include "G4PhysicsLogVector.hh" << 35 #include "G4PhysicsLnVector.hh"
36 #include "G4Proton.hh" 36 #include "G4Proton.hh"
37 #include "G4Neutron.hh" 37 #include "G4Neutron.hh"
38 38
39 const G4double G4XNNElasticLowE::_lowLimit = 0 39 const G4double G4XNNElasticLowE::_lowLimit = 0.;
40 const G4double G4XNNElasticLowE::_highLimit = 40 const G4double G4XNNElasticLowE::_highLimit = 3.*GeV;
41 41
42 // Low energy limit of the cross-section table 42 // Low energy limit of the cross-section table (in GeV)
43 // Units are assigned when filling the Physics 43 // Units are assigned when filling the PhysicsVector
44 const G4double G4XNNElasticLowE::_eMinTable = 44 const G4double G4XNNElasticLowE::_eMinTable = 1.8964808;
45 const G4double G4XNNElasticLowE::_eStepLog = 0 45 const G4double G4XNNElasticLowE::_eStepLog = 0.01;
46 46
47 // Cross-sections in mb 47 // Cross-sections in mb
48 // Units are assigned when filling the Physics 48 // Units are assigned when filling the PhysicsVector
49 49
50 const G4int G4XNNElasticLowE::tableSize = 101; 50 const G4int G4XNNElasticLowE::tableSize = 101;
51 51
52 const G4double G4XNNElasticLowE::ppTable[101] 52 const G4double G4XNNElasticLowE::ppTable[101] =
53 { 53 {
54 60.00, //was 0. 54 60.00, //was 0.
55 33.48, 26.76, 25.26, 24.55, 23.94, 23.77, 23 55 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 56 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 57 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 58 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 59 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 60 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 61 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 62 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 63 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 64 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 65 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 66 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 67 9.90, 9.87, 9.84, 9.80
68 }; 68 };
69 69
70 const G4double G4XNNElasticLowE::npTable[101] 70 const G4double G4XNNElasticLowE::npTable[101] =
71 { 71 {
72 1500.00, // was 0. 72 1500.00, // was 0.
73 248.20, 93.38, 55.26, 44.50, 41.33, 38.48, 3 73 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 74 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 75 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 76 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 77 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 78 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 79 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 80 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 81 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 82 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. 83 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. 84 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 85 9.52, 9.49, 9.46, 9.43
86 }; 86 };
87 87
88 88
89 G4XNNElasticLowE::G4XNNElasticLowE() 89 G4XNNElasticLowE::G4XNNElasticLowE()
90 { 90 {
91 // Cross-sections are available in the range 91 // Cross-sections are available in the range (_eMin,_eMax)
92 92
93 _eMin = _eMinTable * GeV; 93 _eMin = _eMinTable * GeV;
94 _eMax = G4Exp(G4Log(_eMinTable) + tableSize 94 _eMax = G4Exp(G4Log(_eMinTable) + tableSize * _eStepLog) * GeV;
95 if (_eMin < _lowLimit) 95 if (_eMin < _lowLimit)
96 throw G4HadronicException(__FILE__, __LINE 96 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");
97 if (_highLimit > _eMax) 97 if (_highLimit > _eMax)
98 throw G4HadronicException(__FILE__, __LINE 98 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - High energy limit not valid");
99 G4PhysicsVector* pp = new G4PhysicsLogVector << 99 G4PhysicsVector* pp = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
100 100
101 _eMin = G4Exp(G4Log(_eMinTable)-_eStepLog)*G 101 _eMin = G4Exp(G4Log(_eMinTable)-_eStepLog)*GeV;
102 if (_eMin < _lowLimit) 102 if (_eMin < _lowLimit)
103 throw G4HadronicException(__FILE__, __LINE 103 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");
104 G4PhysicsVector* np = new G4PhysicsLogVector << 104 G4PhysicsVector* np = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
105 105
106 G4int i; 106 G4int i;
107 for (i=0; i<tableSize; i++) 107 for (i=0; i<tableSize; i++)
108 { 108 {
109 G4double value = ppTable[i] * millibarn; 109 G4double value = ppTable[i] * millibarn;
110 pp->PutValue(i,value); 110 pp->PutValue(i,value);
111 value = npTable[i] * millibarn; 111 value = npTable[i] * millibarn;
112 np->PutValue(i,value); 112 np->PutValue(i,value);
113 } 113 }
114 xMap[G4Proton::ProtonDefinition()] = pp; 114 xMap[G4Proton::ProtonDefinition()] = pp;
115 xMap[G4Neutron::NeutronDefinition()] = np; 115 xMap[G4Neutron::NeutronDefinition()] = np;
116 } 116 }
117 117
118 118
119 G4XNNElasticLowE::~G4XNNElasticLowE() 119 G4XNNElasticLowE::~G4XNNElasticLowE()
120 { 120 {
121 delete xMap[G4Proton::ProtonDefinition()]; 121 delete xMap[G4Proton::ProtonDefinition()];
122 delete xMap[G4Neutron::NeutronDefinition()]; 122 delete xMap[G4Neutron::NeutronDefinition()];
123 } 123 }
124 124
125 125
126 G4bool G4XNNElasticLowE::operator==(const G4XN 126 G4bool G4XNNElasticLowE::operator==(const G4XNNElasticLowE &right) const
127 { 127 {
128 return (this == (G4XNNElasticLowE *) &right) 128 return (this == (G4XNNElasticLowE *) &right);
129 } 129 }
130 130
131 131
132 G4bool G4XNNElasticLowE::operator!=(const G4XN 132 G4bool G4XNNElasticLowE::operator!=(const G4XNNElasticLowE &right) const
133 { 133 {
134 134
135 return (this != (G4XNNElasticLowE *) &right) 135 return (this != (G4XNNElasticLowE *) &right);
136 } 136 }
137 137
138 138
139 G4double G4XNNElasticLowE::CrossSection(const 139 G4double G4XNNElasticLowE::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const
140 { 140 {
141 G4double sigma = 0.; 141 G4double sigma = 0.;
142 G4double sqrtS = (trk1.Get4Momentum() + trk2 142 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
143 G4bool dummy = false; 143 G4bool dummy = false;
144 144
145 const G4ParticleDefinition * key = FindKeyPa 145 const G4ParticleDefinition * key = FindKeyParticle(trk1,trk2);
146 146
147 typedef std::map <const G4ParticleDefinition 147 typedef std::map <const G4ParticleDefinition *, G4PhysicsVector*, std::less<const G4ParticleDefinition *> > StringPhysMap;
148 148
149 if (xMap.find(key)!= xMap.end()) 149 if (xMap.find(key)!= xMap.end())
150 { 150 {
151 151
152 StringPhysMap::const_iterator iter; 152 StringPhysMap::const_iterator iter;
153 for (iter = xMap.begin(); iter != xMap.e 153 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
154 { 154 {
155 const G4ParticleDefinition * str = (*iter) 155 const G4ParticleDefinition * str = (*iter).first;
156 if (str == key) 156 if (str == key)
157 { 157 {
158 G4PhysicsVector* physVector = (*iter). 158 G4PhysicsVector* physVector = (*iter).second;
159 // G4PhysicsVector* physVector = x 159 // G4PhysicsVector* physVector = xMap[key];
160 if (sqrtS >= _eMin && sqrtS <= _eMax) 160 if (sqrtS >= _eMin && sqrtS <= _eMax)
161 { 161 {
162 sigma = physVector->GetValue(sqrtS,dummy 162 sigma = physVector->GetValue(sqrtS,dummy);
163 } else if ( sqrtS < _eMin ) 163 } else if ( sqrtS < _eMin )
164 { 164 {
165 sigma = physVector->GetValue 165 sigma = physVector->GetValue(_eMin,dummy);
166 } 166 }
167 //G4cout << " sqrtS / sigma " << sqrtS/GeV 167 //G4cout << " sqrtS / sigma " << sqrtS/GeV << " / " <<
168 // sigma/millibarn << G4endl; 168 // sigma/millibarn << G4endl;
169 } 169 }
170 } 170 }
171 } 171 }
172 return sigma; 172 return sigma;
173 } 173 }
174 174
175 175
176 void G4XNNElasticLowE::Print() const 176 void G4XNNElasticLowE::Print() const
177 { 177 {
178 // Dump the pp cross-section table 178 // Dump the pp cross-section table
179 179
180 G4cout << Name() << ", pp cross-section: " < 180 G4cout << Name() << ", pp cross-section: " << G4endl;
181 181
182 G4bool dummy = false; 182 G4bool dummy = false;
183 G4int i; 183 G4int i;
184 const G4ParticleDefinition * key = G4Proton: 184 const G4ParticleDefinition * key = G4Proton::ProtonDefinition();
185 G4PhysicsVector* pp = 0; 185 G4PhysicsVector* pp = 0;
186 186
187 typedef std::map <const G4ParticleDefinition 187 typedef std::map <const G4ParticleDefinition *, G4PhysicsVector*, std::less<const G4ParticleDefinition *> > StringPhysMap;
188 StringPhysMap::const_iterator iter; 188 StringPhysMap::const_iterator iter;
189 189
190 for (iter = xMap.begin(); iter != xMap.end() 190 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
191 { 191 {
192 const G4ParticleDefinition * str = (*ite 192 const G4ParticleDefinition * str = (*iter).first;
193 if (str == key) 193 if (str == key)
194 { 194 {
195 pp = (*iter).second; 195 pp = (*iter).second;
196 } 196 }
197 } 197 }
198 198
199 if (pp != 0) 199 if (pp != 0)
200 { 200 {
201 for (i=0; i<tableSize; i++) 201 for (i=0; i<tableSize; i++)
202 { 202 {
203 G4double e = pp->GetLowEdgeEnergy(i); 203 G4double e = pp->GetLowEdgeEnergy(i);
204 G4double sigma = pp->GetValue(e,dummy) / m 204 G4double sigma = pp->GetValue(e,dummy) / millibarn;
205 G4cout << i << ") e = " << e / GeV << " Ge 205 G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
206 } 206 }
207 } 207 }
208 208
209 // Dump the np cross-section table 209 // Dump the np cross-section table
210 210
211 G4cout << Name() << ", np cross-section: " < 211 G4cout << Name() << ", np cross-section: " << G4endl;
212 212
213 key = G4Neutron::NeutronDefinition(); 213 key = G4Neutron::NeutronDefinition();
214 G4PhysicsVector* np = 0; 214 G4PhysicsVector* np = 0;
215 for (iter = xMap.begin(); iter != xMap.end() 215 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
216 { 216 {
217 const G4ParticleDefinition * str = (*ite 217 const G4ParticleDefinition * str = (*iter).first;
218 if (str == key) 218 if (str == key)
219 { 219 {
220 np = (*iter).second; 220 np = (*iter).second;
221 } 221 }
222 } 222 }
223 223
224 // G4PhysicsVector* np = xMap[G4Neutron::Ne 224 // G4PhysicsVector* np = xMap[G4Neutron::NeutronDefinition()->GetParticleName()];
225 225
226 if (np != 0) 226 if (np != 0)
227 { 227 {
228 for (i=0; i<tableSize; i++) 228 for (i=0; i<tableSize; i++)
229 { 229 {
230 G4double e = np->GetLowEdgeEnergy(i); 230 G4double e = np->GetLowEdgeEnergy(i);
231 G4double sigma = np->GetValue(e,dummy) / m 231 G4double sigma = np->GetValue(e,dummy) / millibarn;
232 G4cout << i << ") e = " << e / GeV << " Ge 232 G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
233 } 233 }
234 } 234 }
235 G4VCrossSectionSource::Print(); 235 G4VCrossSectionSource::Print();
236 } 236 }
237 237
238 238
239 G4String G4XNNElasticLowE::Name() const 239 G4String G4XNNElasticLowE::Name() const
240 { 240 {
241 G4String name("NNElasticLowE"); 241 G4String name("NNElasticLowE");
242 return name; 242 return name;
243 } 243 }
244 244
245 245
246 246
247 G4bool G4XNNElasticLowE::IsValid(G4double e) c 247 G4bool G4XNNElasticLowE::IsValid(G4double e) const
248 { 248 {
249 G4bool answer = InLimits(e,_lowLimit,_highLi 249 G4bool answer = InLimits(e,_lowLimit,_highLimit);
250 250
251 return answer; 251 return answer;
252 } 252 }
253 253
254 254
255 255