| 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 // $Id: G4DNASancheExcitationModel.cc,v 1.4 2010-11-11 22:32:22 sincerti Exp $
>> 27 // GEANT4 tag $Name: not supported by cvs2svn $
26 // 28 //
27 29
28 // Created by Z. Francis 30 // Created by Z. Francis
29 31
30 #include "G4DNASancheExcitationModel.hh" 32 #include "G4DNASancheExcitationModel.hh"
31 #include "G4SystemOfUnits.hh" <<
32 #include "G4DNAMolecularMaterial.hh" <<
33 33
34 //....oooOO0OOooo........oooOO0OOooo........oo << 34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35 35
36 using namespace std; 36 using namespace std;
37 37
38 //#define SANCHE_VERBOSE << 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 <<
40 //....oooOO0OOooo........oooOO0OOooo........oo <<
41 39
42 G4DNASancheExcitationModel::G4DNASancheExcitat 40 G4DNASancheExcitationModel::G4DNASancheExcitationModel(const G4ParticleDefinition*,
43 << 41 const G4String& nam)
44 G4VEmModel(nam) << 42 :G4VEmModel(nam),isInitialised(false)
45 { 43 {
46 fpWaterDensity = nullptr; << 44 nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
47 45
48 SetLowEnergyLimit(2.*eV); << 46 lowEnergyLimit = 2 * eV;
49 SetHighEnergyLimit(100*eV); << 47 highEnergyLimit = 100 * eV;
>> 48 SetLowEnergyLimit(lowEnergyLimit);
>> 49 SetHighEnergyLimit(highEnergyLimit);
50 nLevels = 9; 50 nLevels = 9;
51 51
52 verboseLevel = 0; << 52 verboseLevel= 0;
53 // Verbosity scale: 53 // Verbosity scale:
54 // 0 = nothing << 54 // 0 = nothing
55 // 1 = warning for energy non-conservation << 55 // 1 = warning for energy non-conservation
56 // 2 = details of energy budget 56 // 2 = details of energy budget
57 // 3 = calculation of cross sections, file o 57 // 3 = calculation of cross sections, file openings, sampling of atoms
58 // 4 = entering in methods 58 // 4 = entering in methods
59 << 59
60 #ifdef SANCHE_VERBOSE <<
61 if (verboseLevel > 0) 60 if (verboseLevel > 0)
62 { 61 {
63 G4cout << "Sanche Excitation model is cons << 62 G4cout << "Sanche Excitation model is constructed " << G4endl
64 << G4endl <<
65 << "Energy range: " 63 << "Energy range: "
66 << LowEnergyLimit() / eV << " eV - << 64 << lowEnergyLimit / eV << " eV - "
67 << HighEnergyLimit() / eV << " eV" << 65 << highEnergyLimit / eV << " eV"
68 << G4endl; 66 << G4endl;
69 } 67 }
70 #endif << 68 fParticleChangeForGamma = 0;
71 <<
72 fParticleChangeForGamma = nullptr; <<
73 fpWaterDensity = nullptr; <<
74 <<
75 // Selection of stationary mode <<
76 <<
77 statCode = false; <<
78 } 69 }
79 70
80 //....oooOO0OOooo........oooOO0OOooo........oo << 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 72
82 G4DNASancheExcitationModel::~G4DNASancheExcita 73 G4DNASancheExcitationModel::~G4DNASancheExcitationModel()
83 = default; << 74 {}
84 75
85 //....oooOO0OOooo........oooOO0OOooo........oo << 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86 77
87 void << 78 void G4DNASancheExcitationModel::Initialise(const G4ParticleDefinition* /*particle*/,
88 G4DNASancheExcitationModel:: << 79 const G4DataVector& /*cuts*/)
89 Initialise(const G4ParticleDefinition* /*parti <<
90 const G4DataVector& /*cuts*/) <<
91 { 80 {
92 81
93 #ifdef SANCHE_VERBOSE << 82
94 if (verboseLevel > 3) 83 if (verboseLevel > 3)
95 { << 84 G4cout << "Calling G4DNASancheExcitationModel::Initialise()" << G4endl;
96 G4cout << "Calling G4DNASancheExcitationMo <<
97 << G4endl; <<
98 } <<
99 #endif <<
100 <<
101 // Energy limits <<
102 85
103 if (LowEnergyLimit() < 2.*eV) << 86 // Energy limits
>> 87
>> 88 if (LowEnergyLimit() < lowEnergyLimit)
104 { 89 {
105 G4Exception("*** WARNING : the G4DNASanche << 90 G4cout << "G4DNASancheExcitationModel: low energy limit increased from " <<
106 "validated below 2 eV !", << 91 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
107 "", JustWarning, ""); << 92 SetLowEnergyLimit(lowEnergyLimit);
108 } << 93 }
109 94
110 if (HighEnergyLimit() > 100.*eV) << 95 if (HighEnergyLimit() > highEnergyLimit)
111 { 96 {
112 G4cout << "G4DNASancheExcitationModel: hig << 97 G4cout << "G4DNASancheExcitationModel: high energy limit decreased from " <<
113 HighEnergyLimit()/eV << " eV to " << 100. << 98 HighEnergyLimit()/eV << " eV to " << highEnergyLimit/eV << " eV" << G4endl;
114 SetHighEnergyLimit(100.*eV); << 99 SetHighEnergyLimit(highEnergyLimit);
115 } 100 }
116 101
117 // 102 //
118 #ifdef SANCHE_VERBOSE <<
119 if (verboseLevel > 0) <<
120 { <<
121 G4cout << "Sanche Excitation model is init <<
122 << "Energy range: " <<
123 << LowEnergyLimit() / eV << " eV - " <<
124 << HighEnergyLimit() / eV << " eV" <<
125 << G4endl; <<
126 } <<
127 #endif <<
128 103
129 // Initialize water density pointer << 104 if (verboseLevel > 0)
130 fpWaterDensity = G4DNAMolecularMaterial::Ins <<
131 GetNumMolPerVolTableFor(G4Material::GetM <<
132 105
133 if (isInitialised) {return;} << 106 G4cout << "Sanche Excitation model is initialized " << G4endl
>> 107 << "Energy range: "
>> 108 << LowEnergyLimit() / eV << " eV - "
>> 109 << HighEnergyLimit() / eV << " eV"
>> 110 << G4endl;
134 111
>> 112 if (isInitialised) { return; }
135 fParticleChangeForGamma = GetParticleChangeF 113 fParticleChangeForGamma = GetParticleChangeForGamma();
136 isInitialised = true; 114 isInitialised = true;
137 115
138 const char *path = G4FindDataDir("G4LEDATA") << 116 char *path = getenv("G4LEDATA");
139 std::ostringstream eFullFileName; 117 std::ostringstream eFullFileName;
140 eFullFileName << path << "/dna/sigma_excitat 118 eFullFileName << path << "/dna/sigma_excitationvib_e_sanche.dat";
141 std::ifstream input(eFullFileName.str().c_st 119 std::ifstream input(eFullFileName.str().c_str());
142 120
143 if (!input) 121 if (!input)
>> 122 {
>> 123 G4Exception("G4DNASancheExcitationModel::Initialise","em0003",
>> 124 FatalException,"Missing data file:/dna/sigma_excitationvib_e_sanche.dat");
>> 125 }
>> 126
>> 127 while(!input.eof())
144 { 128 {
145 G4Exception("G4DNASancheExcitationModel::I << 129 double t;
146 FatalException,"Missing data file:/dna << 130 input>>t;
>> 131 tdummyVec.push_back(t);
>> 132 input>>map1[t][0]>>map1[t][1]>>map1[t][2]>>map1[t][3]>>map1[t][4]>>map1[t][5]>>map1[t][6]>>map1[t][7]>>map1[t][8];
>> 133 //G4cout<<t<<" "<<map1[t][0]<<map1[t][1]<<map1[t][2]<<map1[t][3]<<map1[t][4]<<map1[t][5]<<map1[t][6]<<map1[t][7]<<map1[t][8]<<G4endl;
147 } 134 }
148 135
149 // March 25th, 2014 - Vaclav Stepan, Sebasti << 136 }
150 // Added clear for MT <<
151 tdummyVec.clear(); <<
152 // <<
153 137
154 G4double t; << 138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155 G4double xs; <<
156 139
157 while(!input.eof()) << 140 G4double G4DNASancheExcitationModel::CrossSectionPerVolume(const G4Material* material,
158 { << 141 const G4ParticleDefinition* particleDefinition,
159 input>>t; << 142 G4double ekin,
160 tdummyVec.push_back(t); << 143 G4double,
>> 144 G4double)
>> 145 {
>> 146 if (verboseLevel > 3)
>> 147 G4cout << "Calling CrossSectionPerVolume() of G4DNASancheExcitationModel" << G4endl;
161 148
162 fEnergyLevelXS.emplace_back(); << 149 // Calculate total cross section for model
163 fEnergyTotalXS.push_back(0); <<
164 std::vector<G4double>& levelXS = fEnergyLe <<
165 levelXS.reserve(9); <<
166 150
167 // G4cout<<t; << 151 G4double sigma=0;
>> 152
>> 153 if (material == nistwater || material->GetBaseMaterial() == nistwater)
>> 154 {
168 155
169 for(size_t i = 0 ; i < 9 ;++i) << 156 if (particleDefinition == G4Electron::ElectronDefinition())
>> 157 {
>> 158 if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
170 { 159 {
171 input>>xs; << 160 sigma = Sum(ekin);
172 levelXS.push_back(xs); <<
173 fEnergyTotalXS.back() += xs; <<
174 // G4cout <<" " << levelXS[i]; <<
175 } 161 }
176 <<
177 // G4cout << G4endl; <<
178 } 162 }
179 } << 163
180 <<
181 //....oooOO0OOooo........oooOO0OOooo........oo <<
182 <<
183 G4double G4DNASancheExcitationModel::CrossSect <<
184 <<
185 #ifdef SANCHE_VERBOSE <<
186 <<
187 #endif <<
188 <<
189 <<
190 <<
191 <<
192 { <<
193 #ifdef SANCHE_VERBOSE <<
194 if (verboseLevel > 3) 164 if (verboseLevel > 3)
195 { 165 {
196 G4cout << "Calling CrossSectionPerVolume() << 166 G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
197 << G4endl; << 167 G4cout << " - Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
198 } << 168 G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
199 #endif << 169 }
200 <<
201 // Calculate total cross section for model <<
202 170
203 G4double sigma = 0.; << 171 } // if water
204 172
205 G4double waterDensity = (*fpWaterDensity)[ma <<
206 173
207 if (ekin >= LowEnergyLimit() && ekin <= High << 174 return sigma*2*material->GetAtomicNumDensityVector()[1];
208 sigma = TotalCrossSection(ekin); << 175 // see papers for factor 2 description
209 176
210 #ifdef SANCHE_VERBOSE <<
211 if (verboseLevel > 2) <<
212 { <<
213 G4cout << "_______________________________ <<
214 G4cout << "=== G4DNASancheExcitationModel <<
215 G4cout << "=== Kinetic energy(eV)=" << eki <<
216 G4cout << "=== Cross section per water mol <<
217 G4cout << "=== Cross section per water mol <<
218 G4cout << "=== G4DNASancheExcitationModel <<
219 } <<
220 #endif <<
221 <<
222 return sigma*2.*waterDensity; <<
223 // see papers for factor 2 description (liqu <<
224 } 177 }
225 178
226 //....oooOO0OOooo........oooOO0OOooo........oo << 179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227 180
228 void G4DNASancheExcitationModel::SampleSeconda << 181 void G4DNASancheExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
229 << 182 const G4MaterialCutsCouple*,
230 << 183 const G4DynamicParticle* aDynamicElectron,
231 << 184 G4double,
232 << 185 G4double)
233 <<
234 { 186 {
235 #ifdef SANCHE_VERBOSE << 187
236 if (verboseLevel > 3) 188 if (verboseLevel > 3)
237 { << 189 G4cout << "Calling SampleSecondaries() of G4DNASancheExcitationModel" << G4endl;
238 G4cout << "Calling SampleSecondaries() of <<
239 << G4endl; <<
240 } <<
241 #endif <<
242 190
243 G4double electronEnergy0 = aDynamicElectron- 191 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
244 G4int level = RandomSelect(electronEnergy0); 192 G4int level = RandomSelect(electronEnergy0);
245 G4double excitationEnergy = VibrationEnergy( << 193 G4double excitationEnergy = VibrationEnergy(level); // levels go from 0 to 8
246 G4double newEnergy = electronEnergy0 - excit 194 G4double newEnergy = electronEnergy0 - excitationEnergy;
247 << 195
248 /* << 196 /*
249 if (electronEnergy0 < highEnergyLimit) << 197 if (electronEnergy0 < highEnergyLimit)
250 { <<
251 if (newEnergy >= lowEnergyLimit) <<
252 { <<
253 fParticleChangeForGamma->ProposeMomentu <<
254 fParticleChangeForGamma->SetProposedKin <<
255 fParticleChangeForGamma->ProposeLocalEn <<
256 } <<
257 <<
258 else <<
259 { <<
260 fParticleChangeForGamma->ProposeTrackSt <<
261 fParticleChangeForGamma->ProposeLocalEn <<
262 } <<
263 } <<
264 */ <<
265 <<
266 if (electronEnergy0 <= HighEnergyLimit() && <<
267 { 198 {
268 << 199 if (newEnergy >= lowEnergyLimit)
269 if (!statCode) << 200 {
270 { <<
271 fParticleChangeForGamma->ProposeMomentum 201 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
272 fParticleChangeForGamma->SetProposedKine 202 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
273 fParticleChangeForGamma->ProposeLocalEne 203 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
274 } 204 }
275 << 205
276 else << 206 else
277 { 207 {
278 fParticleChangeForGamma->ProposeMomentum << 208 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
279 fParticleChangeForGamma->SetProposedKine << 209 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
280 fParticleChangeForGamma->ProposeLocalEne <<
281 } 210 }
282 <<
283 } 211 }
>> 212 */
284 213
285 // << 214 if (electronEnergy0 < highEnergyLimit && newEnergy>0.)
286 } << 215 {
287 << 216 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
288 //....oooOO0OOooo........oooOO0OOooo........oo << 217 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
289 << 218 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
290 G4double G4DNASancheExcitationModel::PartialCr << 219 }
291 <<
292 { <<
293 // Protection against out of boundary access <<
294 if (t/eV==tdummyVec.back()) t=t*(1.-1e-12); <<
295 // <<
296 <<
297 auto t2 = std::upper_bound(tdummyVec.begin() <<
298 <<
299 auto t1 = t2 - 1; <<
300 <<
301 size_t i1 = t1 - tdummyVec.begin(); <<
302 size_t i2 = t2 - tdummyVec.begin(); <<
303 <<
304 G4double sigma = LinInterpolate((*t1), (*t2) <<
305 t / eV, <<
306 fEnergyLevelXS <<
307 fEnergyLevelXS <<
308 <<
309 static const G4double conv_factor = 1e-16 * <<
310 220
311 sigma *= conv_factor; << 221 //
312 if (sigma == 0.) sigma = 1e-30; <<
313 return (sigma); <<
314 } 222 }
315 223
316 //....oooOO0OOooo........oooOO0OOooo........oo 224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
317 225
318 G4double G4DNASancheExcitationModel::TotalCros << 226 G4double G4DNASancheExcitationModel::PartialCrossSection(G4double t, G4int level)
319 { 227 {
320 // Protection against out of boundary access << 228 std::vector<double>::iterator t2 = std::upper_bound(tdummyVec.begin(),tdummyVec.end(), t/eV);
321 if (t/eV==tdummyVec.back()) t=t*(1.-1e-12); << 229 std::vector<double>::iterator t1 = t2-1;
322 // <<
323 <<
324 auto t2 = std::upper_bound(tdummyVec.begin() <<
325 <<
326 auto t1 = t2 - 1; <<
327 <<
328 size_t i1 = t1 - tdummyVec.begin(); <<
329 size_t i2 = t2 - tdummyVec.begin(); <<
330 <<
331 G4double sigma = LinInterpolate((*t1), (*t2) <<
332 t / eV, <<
333 fEnergyTotalXS <<
334 fEnergyTotalXS <<
335 230
336 static const G4double conv_factor = 1e-16 * << 231 double sigma = LinInterpolate((*t1), (*t2), t/eV, map1[*t1][level], map1[*t2][level]);
337 << 232 sigma*=1e-16*cm*cm;
338 sigma *= conv_factor; << 233 if(sigma==0.)sigma=1e-30;
339 if (sigma == 0.) sigma = 1e-30; <<
340 return (sigma); 234 return (sigma);
341 } 235 }
342 236
343 //....oooOO0OOooo........oooOO0OOooo........oo 237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
344 238
345 G4double G4DNASancheExcitationModel::Vibration 239 G4double G4DNASancheExcitationModel::VibrationEnergy(G4int level)
346 { 240 {
347 static G4double energies[9] = { 0.01, 0.024, << 241 G4double energies[9] = {0.01, 0.024, 0.061, 0.092, 0.204, 0.417, 0.460, 0.500, 0.835};
348 0.500, 0.835 }; << 242 return(energies[level]*eV);
349 return (energies[level] * eV); <<
350 } 243 }
351 244
352 //....oooOO0OOooo........oooOO0OOooo........oo 245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
353 246
354 G4int G4DNASancheExcitationModel::RandomSelect 247 G4int G4DNASancheExcitationModel::RandomSelect(G4double k)
355 { 248 {
356 249
357 // Level Selection Counting can be done here 250 // Level Selection Counting can be done here !
358 251
359 G4int i = nLevels; 252 G4int i = nLevels;
360 G4double value = 0.; 253 G4double value = 0.;
361 std::deque<G4double> values; << 254 std::deque<double> values;
362 << 255
363 while (i > 0) 256 while (i > 0)
364 { 257 {
365 i--; 258 i--;
366 G4double partial = PartialCrossSection(k, << 259 G4double partial = PartialCrossSection(k,i);
367 values.push_front(partial); 260 values.push_front(partial);
368 value += partial; 261 value += partial;
369 } 262 }
370 263
371 value *= G4UniformRand(); 264 value *= G4UniformRand();
372 << 265
373 i = nLevels; 266 i = nLevels;
374 267
375 while (i > 0) 268 while (i > 0)
376 { 269 {
377 i--; 270 i--;
378 if (values[i] > value) 271 if (values[i] > value)
379 { 272 {
380 //outcount<<i<<" "<<VibrationEnergy(i)< << 273 //outcount<<i<<" "<<VibrationEnergy(i)<<G4endl;
381 return i; 274 return i;
382 } 275 }
383 value -= values[i]; 276 value -= values[i];
384 } 277 }
385 << 278
386 //outcount<<0<<" "<<VibrationEnergy(0)<<G4e 279 //outcount<<0<<" "<<VibrationEnergy(0)<<G4endl;
387 << 280
388 return 0; 281 return 0;
389 } 282 }
390 283
391 //....oooOO0OOooo........oooOO0OOooo........oo 284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
392 285
393 G4double G4DNASancheExcitationModel::Sum(G4dou 286 G4double G4DNASancheExcitationModel::Sum(G4double k)
394 { 287 {
395 G4double totalCrossSection = 0.; 288 G4double totalCrossSection = 0.;
396 289
397 for (G4int i = 0; i < nLevels; i++) << 290 for (G4int i=0; i<nLevels; i++)
398 { 291 {
399 totalCrossSection += PartialCrossSection(k << 292 totalCrossSection += PartialCrossSection(k,i);
400 } 293 }
401 <<
402 return totalCrossSection; 294 return totalCrossSection;
403 } 295 }
404 296
405 //....oooOO0OOooo........oooOO0OOooo........oo 297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
406 298
407 G4double G4DNASancheExcitationModel::LinInterp << 299 G4double G4DNASancheExcitationModel::LinInterpolate(G4double e1,
408 << 300 G4double e2,
409 << 301 G4double e,
410 << 302 G4double xs1,
411 << 303 G4double xs2)
412 { 304 {
413 G4double a = (xs2 - xs1) / (e2 - e1); 305 G4double a = (xs2 - xs1) / (e2 - e1);
414 G4double b = xs2 - a * e2; << 306 G4double b = xs2 - a*e2;
415 G4double value = a * e + b; << 307 G4double value = a*e + b;
416 // G4cout<<"interP >> "<<e1<<" "<<e2<<" " << 308 // G4cout<<"interP >> "<<e1<<" "<<e2<<" "<<e<<" "<<xs1<<" "<<xs2<<" "<<a<<" "<<b<<" "<<value<<G4endl;
417 // <<xs1<<" "<<xs2<<" "<<a<<" "<<b<<" "< << 309
418 <<
419 return value; 310 return value;
420 } 311 }
421 312
422 313