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