| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 #include "globals.hh" 23 #include "globals.hh"
27 #include "G4ios.hh" 24 #include "G4ios.hh"
28 #include "G4Exp.hh" <<
29 #include "G4Log.hh" <<
30 #include "G4SystemOfUnits.hh" <<
31 #include "G4XNNElasticLowE.hh" 25 #include "G4XNNElasticLowE.hh"
32 #include "G4KineticTrack.hh" 26 #include "G4KineticTrack.hh"
33 #include "G4ParticleDefinition.hh" 27 #include "G4ParticleDefinition.hh"
34 #include "G4PhysicsVector.hh" 28 #include "G4PhysicsVector.hh"
35 #include "G4PhysicsLogVector.hh" << 29 #include "G4PhysicsLnVector.hh"
>> 30 #include "G4ParticleDefinition.hh"
36 #include "G4Proton.hh" 31 #include "G4Proton.hh"
37 #include "G4Neutron.hh" 32 #include "G4Neutron.hh"
38 33
>> 34
39 const G4double G4XNNElasticLowE::_lowLimit = 0 35 const G4double G4XNNElasticLowE::_lowLimit = 0.;
40 const G4double G4XNNElasticLowE::_highLimit = 36 const G4double G4XNNElasticLowE::_highLimit = 3.*GeV;
41 37
42 // Low energy limit of the cross-section table 38 // Low energy limit of the cross-section table (in GeV)
43 // Units are assigned when filling the Physics << 39 // Units are assigned while filling the PhysicsVector
44 const G4double G4XNNElasticLowE::_eMinTable = 40 const G4double G4XNNElasticLowE::_eMinTable = 1.8964808;
45 const G4double G4XNNElasticLowE::_eStepLog = 0 41 const G4double G4XNNElasticLowE::_eStepLog = 0.01;
46 42
47 // Cross-sections in mb 43 // Cross-sections in mb
48 // Units are assigned when filling the Physics << 44 // Units are assigned while filling the PhysicsVector
49 45
50 const G4int G4XNNElasticLowE::tableSize = 101; 46 const G4int G4XNNElasticLowE::tableSize = 101;
51 47
52 const G4double G4XNNElasticLowE::ppTable[101] 48 const G4double G4XNNElasticLowE::ppTable[101] =
53 { 49 {
54 60.00, //was 0. 50 60.00, //was 0.
55 33.48, 26.76, 25.26, 24.55, 23.94, 23.77, 23 51 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 52 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 53 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 54 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 55 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 56 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 57 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 58 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 59 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 60 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 61 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 62 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 63 9.90, 9.87, 9.84, 9.80
68 }; 64 };
69 65
70 const G4double G4XNNElasticLowE::npTable[101] 66 const G4double G4XNNElasticLowE::npTable[101] =
71 { 67 {
72 1500.00, // was 0. 68 1500.00, // was 0.
73 248.20, 93.38, 55.26, 44.50, 41.33, 38.48, 3 69 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 70 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 71 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 72 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 73 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 74 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 75 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 76 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 77 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 78 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. 79 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. 80 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 81 9.52, 9.49, 9.46, 9.43
86 }; 82 };
87 83
88 84
89 G4XNNElasticLowE::G4XNNElasticLowE() 85 G4XNNElasticLowE::G4XNNElasticLowE()
90 { 86 {
91 // Cross-sections are available in the range 87 // Cross-sections are available in the range (_eMin,_eMax)
92 88
93 _eMin = _eMinTable * GeV; 89 _eMin = _eMinTable * GeV;
94 _eMax = G4Exp(G4Log(_eMinTable) + tableSize << 90 _eMax = exp(log(_eMinTable) + tableSize * _eStepLog) * GeV;
95 if (_eMin < _lowLimit) 91 if (_eMin < _lowLimit)
96 throw G4HadronicException(__FILE__, __LINE 92 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");
97 if (_highLimit > _eMax) 93 if (_highLimit > _eMax)
98 throw G4HadronicException(__FILE__, __LINE 94 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - High energy limit not valid");
99 G4PhysicsVector* pp = new G4PhysicsLogVector << 95 G4PhysicsVector* pp = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
100 96
101 _eMin = G4Exp(G4Log(_eMinTable)-_eStepLog)*G << 97 _eMin = exp(log(_eMinTable)-_eStepLog)*GeV;
102 if (_eMin < _lowLimit) 98 if (_eMin < _lowLimit)
103 throw G4HadronicException(__FILE__, __LINE 99 throw G4HadronicException(__FILE__, __LINE__, "G4XNNElasticLowE::G4XNNElasticLowE - Low energy limit not valid");
104 G4PhysicsVector* np = new G4PhysicsLogVector << 100 G4PhysicsVector* np = new G4PhysicsLnVector(_eMin,_eMax,tableSize);
105 101
106 G4int i; 102 G4int i;
107 for (i=0; i<tableSize; i++) 103 for (i=0; i<tableSize; i++)
108 { 104 {
109 G4double value = ppTable[i] * millibarn; 105 G4double value = ppTable[i] * millibarn;
110 pp->PutValue(i,value); 106 pp->PutValue(i,value);
111 value = npTable[i] * millibarn; 107 value = npTable[i] * millibarn;
112 np->PutValue(i,value); 108 np->PutValue(i,value);
113 } 109 }
114 xMap[G4Proton::ProtonDefinition()] = pp; << 110 xMap[G4Proton::ProtonDefinition()->GetParticleName()] = pp;
115 xMap[G4Neutron::NeutronDefinition()] = np; << 111 xMap[G4Neutron::NeutronDefinition()->GetParticleName()] = np;
116 } 112 }
117 113
118 114
119 G4XNNElasticLowE::~G4XNNElasticLowE() 115 G4XNNElasticLowE::~G4XNNElasticLowE()
120 { 116 {
121 delete xMap[G4Proton::ProtonDefinition()]; << 117 delete xMap[G4Proton::ProtonDefinition()->GetParticleName()];
122 delete xMap[G4Neutron::NeutronDefinition()]; << 118 delete xMap[G4Neutron::NeutronDefinition()->GetParticleName()];
123 } 119 }
124 120
125 121
126 G4bool G4XNNElasticLowE::operator==(const G4XN 122 G4bool G4XNNElasticLowE::operator==(const G4XNNElasticLowE &right) const
127 { 123 {
128 return (this == (G4XNNElasticLowE *) &right) 124 return (this == (G4XNNElasticLowE *) &right);
129 } 125 }
130 126
131 127
132 G4bool G4XNNElasticLowE::operator!=(const G4XN 128 G4bool G4XNNElasticLowE::operator!=(const G4XNNElasticLowE &right) const
133 { 129 {
134 130
135 return (this != (G4XNNElasticLowE *) &right) 131 return (this != (G4XNNElasticLowE *) &right);
136 } 132 }
137 133
138 134
139 G4double G4XNNElasticLowE::CrossSection(const 135 G4double G4XNNElasticLowE::CrossSection(const G4KineticTrack& trk1, const G4KineticTrack& trk2) const
140 { 136 {
141 G4double sigma = 0.; 137 G4double sigma = 0.;
142 G4double sqrtS = (trk1.Get4Momentum() + trk2 138 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
143 G4bool dummy = false; 139 G4bool dummy = false;
144 140
145 const G4ParticleDefinition * key = FindKeyPa << 141 G4String key = FindKeyParticle(trk1,trk2);
146 142
147 typedef std::map <const G4ParticleDefinition << 143 typedef std::map <G4String, G4PhysicsVector*, std::less<G4String> > StringPhysMap;
148 144
149 if (xMap.find(key)!= xMap.end()) 145 if (xMap.find(key)!= xMap.end())
150 { 146 {
151 147
152 StringPhysMap::const_iterator iter; 148 StringPhysMap::const_iterator iter;
153 for (iter = xMap.begin(); iter != xMap.e 149 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
154 { 150 {
155 const G4ParticleDefinition * str = (*iter) << 151 G4String str = (*iter).first;
156 if (str == key) 152 if (str == key)
157 { 153 {
158 G4PhysicsVector* physVector = (*iter). 154 G4PhysicsVector* physVector = (*iter).second;
159 // G4PhysicsVector* physVector = x 155 // G4PhysicsVector* physVector = xMap[key];
160 if (sqrtS >= _eMin && sqrtS <= _eMax) 156 if (sqrtS >= _eMin && sqrtS <= _eMax)
161 { 157 {
162 sigma = physVector->GetValue(sqrtS,dummy 158 sigma = physVector->GetValue(sqrtS,dummy);
163 } else if ( sqrtS < _eMin ) 159 } else if ( sqrtS < _eMin )
164 { 160 {
165 sigma = physVector->GetValue 161 sigma = physVector->GetValue(_eMin,dummy);
166 } 162 }
167 //G4cout << " sqrtS / sigma " << sqrtS/GeV <<
168 // sigma/millibarn << G4endl; <<
169 } 163 }
170 } 164 }
171 } 165 }
172 return sigma; 166 return sigma;
173 } 167 }
174 168
175 169
176 void G4XNNElasticLowE::Print() const 170 void G4XNNElasticLowE::Print() const
177 { 171 {
178 // Dump the pp cross-section table 172 // Dump the pp cross-section table
179 173
180 G4cout << Name() << ", pp cross-section: " < 174 G4cout << Name() << ", pp cross-section: " << G4endl;
181 175
182 G4bool dummy = false; 176 G4bool dummy = false;
183 G4int i; 177 G4int i;
184 const G4ParticleDefinition * key = G4Proton: << 178 G4String key = G4Proton::ProtonDefinition()->GetParticleName();
185 G4PhysicsVector* pp = 0; 179 G4PhysicsVector* pp = 0;
186 180
187 typedef std::map <const G4ParticleDefinition << 181 typedef std::map <G4String, G4PhysicsVector*, std::less<G4String> > StringPhysMap;
188 StringPhysMap::const_iterator iter; 182 StringPhysMap::const_iterator iter;
189 183
190 for (iter = xMap.begin(); iter != xMap.end() 184 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
191 { 185 {
192 const G4ParticleDefinition * str = (*ite << 186 G4String str = (*iter).first;
193 if (str == key) 187 if (str == key)
194 { 188 {
195 pp = (*iter).second; 189 pp = (*iter).second;
196 } 190 }
197 } 191 }
198 192
199 if (pp != 0) 193 if (pp != 0)
200 { 194 {
201 for (i=0; i<tableSize; i++) 195 for (i=0; i<tableSize; i++)
202 { 196 {
203 G4double e = pp->GetLowEdgeEnergy(i); 197 G4double e = pp->GetLowEdgeEnergy(i);
204 G4double sigma = pp->GetValue(e,dummy) / m 198 G4double sigma = pp->GetValue(e,dummy) / millibarn;
205 G4cout << i << ") e = " << e / GeV << " Ge 199 G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
206 } 200 }
207 } 201 }
208 202
209 // Dump the np cross-section table 203 // Dump the np cross-section table
210 204
211 G4cout << Name() << ", np cross-section: " < 205 G4cout << Name() << ", np cross-section: " << G4endl;
212 206
213 key = G4Neutron::NeutronDefinition(); << 207 key = G4Neutron::NeutronDefinition()->GetParticleName();
214 G4PhysicsVector* np = 0; 208 G4PhysicsVector* np = 0;
215 for (iter = xMap.begin(); iter != xMap.end() 209 for (iter = xMap.begin(); iter != xMap.end(); ++iter)
216 { 210 {
217 const G4ParticleDefinition * str = (*ite << 211 G4String str = (*iter).first;
218 if (str == key) 212 if (str == key)
219 { 213 {
220 np = (*iter).second; 214 np = (*iter).second;
221 } 215 }
222 } 216 }
223 217
224 // G4PhysicsVector* np = xMap[G4Neutron::Ne 218 // G4PhysicsVector* np = xMap[G4Neutron::NeutronDefinition()->GetParticleName()];
225 219
226 if (np != 0) 220 if (np != 0)
227 { 221 {
228 for (i=0; i<tableSize; i++) 222 for (i=0; i<tableSize; i++)
229 { 223 {
230 G4double e = np->GetLowEdgeEnergy(i); 224 G4double e = np->GetLowEdgeEnergy(i);
231 G4double sigma = np->GetValue(e,dummy) / m 225 G4double sigma = np->GetValue(e,dummy) / millibarn;
232 G4cout << i << ") e = " << e / GeV << " Ge 226 G4cout << i << ") e = " << e / GeV << " GeV ---- Cross section = " << sigma << " mb " << G4endl;
233 } 227 }
234 } 228 }
235 G4VCrossSectionSource::Print(); 229 G4VCrossSectionSource::Print();
236 } 230 }
237 231
238 232
239 G4String G4XNNElasticLowE::Name() const 233 G4String G4XNNElasticLowE::Name() const
240 { 234 {
241 G4String name("NNElasticLowE"); 235 G4String name("NNElasticLowE");
242 return name; 236 return name;
243 } 237 }
244 238
245 239
246 240
247 G4bool G4XNNElasticLowE::IsValid(G4double e) c 241 G4bool G4XNNElasticLowE::IsValid(G4double e) const
248 { 242 {
249 G4bool answer = InLimits(e,_lowLimit,_highLi 243 G4bool answer = InLimits(e,_lowLimit,_highLimit);
250 244
251 return answer; 245 return answer;
252 } 246 }
253 247
254 248
255 249