| 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 // ******************************************* 26 // *********************************************************************
27 // | 27 // | |
28 // | G4LowEPComptonModel-- Geant4 28 // | G4LowEPComptonModel-- Geant4 Monash University |
29 // | low energy Compton scat 29 // | low energy Compton scattering model. |
30 // | J. M. C. Brown, Monash Univer 30 // | J. M. C. Brown, Monash University, Australia |
31 // | ## Unpolarised photons 31 // | ## Unpolarised photons only ## |
32 // | 32 // | |
33 // | 33 // | |
34 // ******************************************* 34 // *********************************************************************
35 // | 35 // | |
36 // | The following is a Geant4 class to simula 36 // | The following is a Geant4 class to simulate the process of |
37 // | bound electron Compton scattering. Genera 37 // | bound electron Compton scattering. General code structure is |
38 // | based on G4LowEnergyCompton.cc and G4Live 38 // | based on G4LowEnergyCompton.cc and G4LivermoreComptonModel.cc. |
39 // | Algorithms for photon energy, and ejected 39 // | Algorithms for photon energy, and ejected Compton electron |
40 // | direction taken from: 40 // | direction taken from: |
41 // | 41 // | |
42 // | J. M. C. Brown, M. R. Dimmock, J. E. Gill 42 // | J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin, |
43 // | "A low energy bound atomic electron Compt 43 // | "A low energy bound atomic electron Compton scattering model |
44 // | for Geant4", NIMB, Vol. 338, 77-88, 2014 << 44 // | for Geant4", NIMA, submitted 2013. |
45 // | 45 // | |
46 // | The author acknowledges the work of the G 46 // | The author acknowledges the work of the Geant4 collaboration |
47 // | in developing the following algorithms th 47 // | in developing the following algorithms that have been employed |
48 // | or adapeted for the present software: 48 // | or adapeted for the present software: |
49 // | 49 // | |
50 // | # sampling of photon scattering angle, 50 // | # sampling of photon scattering angle, |
51 // | # target element selection in composite 51 // | # target element selection in composite materials, |
52 // | # target shell selection in element, 52 // | # target shell selection in element, |
53 // | # and sampling of bound electron momentu 53 // | # and sampling of bound electron momentum from Compton profiles. |
54 // | 54 // | |
55 // ******************************************* 55 // *********************************************************************
56 // | 56 // | |
57 // | History: 57 // | History: |
58 // | -------- 58 // | -------- |
59 // | 59 // | |
60 // | Nov. 2011 JMCB - First version 60 // | Nov. 2011 JMCB - First version |
61 // | Feb. 2012 JMCB - Migration to Geant 61 // | Feb. 2012 JMCB - Migration to Geant4 9.5 |
62 // | Sep. 2012 JMCB - Final fixes for Ge 62 // | Sep. 2012 JMCB - Final fixes for Geant4 9.6 |
63 // | Feb. 2013 JMCB - Geant4 9.6 FPE fix 63 // | Feb. 2013 JMCB - Geant4 9.6 FPE fix for bug 1426 |
64 // | Jan. 2015 JMCB - Migration to MT (B <<
65 // | implementation) <<
66 // | Feb. 2016 JMCB - Geant4 10.2 FPE fi <<
67 // | 64 // | |
68 // ******************************************* 65 // *********************************************************************
69 66
70 #include <limits> 67 #include <limits>
71 #include "G4LowEPComptonModel.hh" 68 #include "G4LowEPComptonModel.hh"
72 #include "G4PhysicalConstants.hh" 69 #include "G4PhysicalConstants.hh"
73 #include "G4SystemOfUnits.hh" 70 #include "G4SystemOfUnits.hh"
74 #include "G4Electron.hh" 71 #include "G4Electron.hh"
75 #include "G4ParticleChangeForGamma.hh" 72 #include "G4ParticleChangeForGamma.hh"
76 #include "G4LossTableManager.hh" 73 #include "G4LossTableManager.hh"
77 #include "G4VAtomDeexcitation.hh" 74 #include "G4VAtomDeexcitation.hh"
78 #include "G4AtomicShell.hh" 75 #include "G4AtomicShell.hh"
79 #include "G4AutoLock.hh" << 76 #include "G4CrossSectionHandler.hh"
>> 77 #include "G4CompositeEMDataSet.hh"
>> 78 #include "G4LogLogInterpolation.hh"
80 #include "G4Gamma.hh" 79 #include "G4Gamma.hh"
81 #include "G4ShellData.hh" <<
82 #include "G4DopplerProfile.hh" <<
83 #include "G4Log.hh" <<
84 #include "G4Exp.hh" <<
85 80
86 //******************************************** 81 //****************************************************************************
87 82
88 using namespace std; 83 using namespace std;
89 namespace { G4Mutex LowEPComptonModelMutex = G <<
90 84
91 G4PhysicsFreeVector* G4LowEPComptonModel::data << 85 //****************************************************************************
92 G4ShellData* G4LowEPComptonModel::shellD <<
93 G4DopplerProfile* G4LowEPComptonModel::profil <<
94 <<
95 static const G4double ln10 = G4Log(10.); <<
96 86
97 G4LowEPComptonModel::G4LowEPComptonModel(const 87 G4LowEPComptonModel::G4LowEPComptonModel(const G4ParticleDefinition*,
98 << 88 const G4String& nam)
99 : G4VEmModel(nam),isInitialised(false) << 89 :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
>> 90 scatterFunctionData(0),crossSectionHandler(0),fAtomDeexcitation(0)
100 { 91 {
101 verboseLevel=1 ; << 92 lowEnergyLimit = 250 * eV;
>> 93 highEnergyLimit = 100 * GeV;
>> 94
>> 95 verboseLevel=0 ;
102 // Verbosity scale: 96 // Verbosity scale:
103 // 0 = nothing 97 // 0 = nothing
104 // 1 = warning for energy non-conservation 98 // 1 = warning for energy non-conservation
105 // 2 = details of energy budget 99 // 2 = details of energy budget
106 // 3 = calculation of cross sections, file o 100 // 3 = calculation of cross sections, file openings, sampling of atoms
107 // 4 = entering in methods 101 // 4 = entering in methods
108 102
109 if( verboseLevel>1 ) { << 103 if( verboseLevel>0 ) {
110 G4cout << "Low energy photon Compton model << 104 G4cout << "Low energy photon Compton model is constructed " << G4endl
>> 105 << "Energy range: "
>> 106 << lowEnergyLimit / eV << " eV - "
>> 107 << highEnergyLimit / GeV << " GeV"
>> 108 << G4endl;
111 } 109 }
112 110
113 //Mark this model as "applicable" for atomic 111 //Mark this model as "applicable" for atomic deexcitation
114 SetDeexcitationFlag(true); 112 SetDeexcitationFlag(true);
115 113
116 fParticleChange = nullptr; <<
117 fAtomDeexcitation = nullptr; <<
118 } 114 }
119 115
120 //******************************************** 116 //****************************************************************************
121 117
122 G4LowEPComptonModel::~G4LowEPComptonModel() 118 G4LowEPComptonModel::~G4LowEPComptonModel()
123 { << 119 {
124 if(IsMaster()) { << 120 delete crossSectionHandler;
125 delete shellData; << 121 delete scatterFunctionData;
126 shellData = nullptr; <<
127 delete profileData; <<
128 profileData = nullptr; <<
129 } <<
130 } 122 }
131 123
132 //******************************************** 124 //****************************************************************************
133 125
134 void G4LowEPComptonModel::Initialise(const G4P 126 void G4LowEPComptonModel::Initialise(const G4ParticleDefinition* particle,
135 const G4D << 127 const G4DataVector& cuts)
136 { 128 {
137 if (verboseLevel > 1) { << 129 if (verboseLevel > 2) {
138 G4cout << "Calling G4LowEPComptonModel::In 130 G4cout << "Calling G4LowEPComptonModel::Initialise()" << G4endl;
139 } 131 }
140 132
141 // Initialise element selector << 133 if (crossSectionHandler)
142 if(IsMaster()) { << 134 {
143 << 135 crossSectionHandler->Clear();
144 // Access to elements << 136 delete crossSectionHandler;
145 const char* path = G4FindDataDir("G4LEDATA << 137 }
146 << 138 delete scatterFunctionData;
147 G4ProductionCutsTable* theCoupleTable = <<
148 G4ProductionCutsTable::GetProductionCuts <<
149 G4int numOfCouples = (G4int)theCoupleTable <<
150 <<
151 for(G4int i=0; i<numOfCouples; ++i) { <<
152 const G4Material* material = <<
153 theCoupleTable->GetMaterialCutsCouple( <<
154 const G4ElementVector* theElementVector <<
155 std::size_t nelm = material->GetNumberOf <<
156 <<
157 for (std::size_t j=0; j<nelm; ++j) { <<
158 G4int Z = G4lrint((*theElementVector)[ <<
159 if(Z < 1) { Z = 1; } <<
160 else if(Z > maxZ){ Z = maxZ; } <<
161 <<
162 if( (!data[Z]) ) { ReadData(Z, path); <<
163 } <<
164 } <<
165 139
166 // For Doppler broadening << 140 // Reading of data files - all materials are read
167 if(!shellData) { << 141 crossSectionHandler = new G4CrossSectionHandler;
168 shellData = new G4ShellData(); << 142 G4String crossSectionFile = "comp/ce-cs-";
169 shellData->SetOccupancyData(); << 143 crossSectionHandler->LoadData(crossSectionFile);
170 G4String file = "/doppler/shell-doppler" << 144
171 shellData->LoadData(file); << 145 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
172 } << 146 G4String scatterFile = "comp/ce-sf-";
173 if(!profileData) { profileData = new G4Dop << 147 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
>> 148 scatterFunctionData->LoadData(scatterFile);
>> 149
>> 150 // For Doppler broadening
>> 151 shellData.SetOccupancyData();
>> 152 G4String file = "/doppler/shell-doppler";
>> 153 shellData.LoadData(file);
174 154
175 InitialiseElementSelectors(particle, cuts) << 155 InitialiseElementSelectors(particle,cuts);
176 } <<
177 156
178 if (verboseLevel > 2) { 157 if (verboseLevel > 2) {
179 G4cout << "Loaded cross section files" << << 158 G4cout << "Loaded cross section files for low energy photon Compton model" << G4endl;
180 } <<
181 <<
182 if( verboseLevel>1 ) { <<
183 G4cout << "G4LowEPComptonModel is initiali <<
184 << "Energy range: " <<
185 << LowEnergyLimit() / eV << " eV - <<
186 << HighEnergyLimit() / GeV << " GeV <<
187 << G4endl; <<
188 } 159 }
189 160
190 if(isInitialised) { return; } 161 if(isInitialised) { return; }
191 <<
192 fParticleChange = GetParticleChangeForGamma( <<
193 fAtomDeexcitation = G4LossTableManager::Ins <<
194 isInitialised = true; 162 isInitialised = true;
195 } <<
196 <<
197 //******************************************** <<
198 <<
199 void G4LowEPComptonModel::InitialiseLocal(cons <<
200 <<
201 { <<
202 SetElementSelectors(masterModel->GetElementS <<
203 } <<
204 <<
205 //******************************************** <<
206 <<
207 void G4LowEPComptonModel::ReadData(std::size_t <<
208 { <<
209 if (verboseLevel > 1) <<
210 { <<
211 G4cout << "G4LowEPComptonModel::ReadData() <<
212 << G4endl; <<
213 } <<
214 if(data[Z]) { return; } <<
215 const char* datadir = path; <<
216 if(!datadir) <<
217 { <<
218 datadir = G4FindDataDir("G4LEDATA"); <<
219 if(!datadir) <<
220 { <<
221 G4Exception("G4LowEPComptonModel::ReadDa <<
222 "em0006",FatalException, <<
223 "Environment variable G4LEDA <<
224 return; <<
225 } <<
226 } <<
227 163
228 data[Z] = new G4PhysicsFreeVector(); << 164 fParticleChange = GetParticleChangeForGamma();
229 165
230 std::ostringstream ost; << 166 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
231 ost << datadir << "/livermore/comp/ce-cs-" < <<
232 std::ifstream fin(ost.str().c_str()); <<
233 167
234 if( !fin.is_open()) << 168 if( verboseLevel>0 ) {
235 { << 169 G4cout << "Low energy photon Compton model is initialized " << G4endl
236 G4ExceptionDescription ed; << 170 << "Energy range: "
237 ed << "G4LowEPComptonModel data file <" << 171 << LowEnergyLimit() / eV << " eV - "
238 << "> is not opened!" << G4endl; << 172 << HighEnergyLimit() / GeV << " GeV"
239 G4Exception("G4LowEPComptonModel::ReadData << 173 << G4endl;
240 "em0003",FatalException, << 174 }
241 ed,"G4LEDATA version should be <<
242 return; <<
243 } else { <<
244 if(verboseLevel > 3) { <<
245 G4cout << "File " << ost.str() <<
246 << " is opened by G4LowEPCompto <<
247 } <<
248 data[Z]->Retrieve(fin, true); <<
249 data[Z]->ScaleVector(MeV, MeV*barn); <<
250 } <<
251 fin.close(); <<
252 } 175 }
253 176
254 //******************************************** 177 //****************************************************************************
255 178
256 << 179 G4double G4LowEPComptonModel::ComputeCrossSectionPerAtom(
257 G4double << 180 const G4ParticleDefinition*,
258 G4LowEPComptonModel::ComputeCrossSectionPerAto << 181 G4double GammaEnergy,
259 << 182 G4double Z, G4double,
260 << 183 G4double, G4double)
261 <<
262 { 184 {
263 if (verboseLevel > 3) { 185 if (verboseLevel > 3) {
264 G4cout << "G4LowEPComptonModel::ComputeCro << 186 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LowEPComptonModel" << G4endl;
265 << G4endl; <<
266 } 187 }
267 G4double cs = 0.0; << 188 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) { return 0.0; }
268 << 189
269 if (GammaEnergy < LowEnergyLimit()) { return << 190 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
270 << 191 return cs;
271 G4int intZ = G4lrint(Z); << 192 }
272 if(intZ < 1 || intZ > maxZ) { return cs; } <<
273 193
274 G4PhysicsFreeVector* pv = data[intZ]; <<
275 194
276 // if element was not initialised <<
277 // do initialisation safely for MT mode <<
278 if(!pv) <<
279 { <<
280 InitialiseForElement(0, intZ); <<
281 pv = data[intZ]; <<
282 if(!pv) { return cs; } <<
283 } <<
284 195
285 G4int n = G4int(pv->GetVectorLength() - 1); <<
286 G4double e1 = pv->Energy(0); <<
287 G4double e2 = pv->Energy(n); <<
288 <<
289 if(GammaEnergy <= e1) { cs = GammaEnerg <<
290 else if(GammaEnergy <= e2) { cs = pv->Value( <<
291 else if(GammaEnergy > e2) { cs = pv->Value( <<
292 196
293 return cs; <<
294 } <<
295 197
296 //******************************************** 198 //****************************************************************************
297 199
>> 200
298 void G4LowEPComptonModel::SampleSecondaries(st 201 void G4LowEPComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
299 << 202 const G4MaterialCutsCouple* couple,
300 << 203 const G4DynamicParticle* aDynamicGamma,
301 << 204 G4double, G4double)
302 { 205 {
303 206
304 // The scattered gamma energy is sampled acc 207 // The scattered gamma energy is sampled according to Klein - Nishina formula.
305 // then accepted or rejected depending on th 208 // then accepted or rejected depending on the Scattering Function multiplied
306 // by factor from Klein - Nishina formula. 209 // by factor from Klein - Nishina formula.
307 // Expression of the angular distribution as 210 // Expression of the angular distribution as Klein Nishina
308 // angular and energy distribution and Scatt 211 // angular and energy distribution and Scattering fuctions is taken from
309 // D. E. Cullen "A simple model of photon tr 212 // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth.
310 // Phys. Res. B 101 (1995). Method of sampli 213 // Phys. Res. B 101 (1995). Method of sampling with form factors is different
311 // data are interpolated while in the articl 214 // data are interpolated while in the article they are fitted.
312 // Reference to the article is from J. Stepa 215 // Reference to the article is from J. Stepanek New Photon, Positron
313 // and Electron Interaction Data for GEANT i 216 // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10
314 // TeV (draft). 217 // TeV (draft).
315 // The random number techniques of Butcher & 218 // The random number techniques of Butcher & Messel are used
316 // (Nucl Phys 20(1960),15). 219 // (Nucl Phys 20(1960),15).
317 220
318 G4double photonEnergy0 = aDynamicGamma->GetK 221 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy()/MeV;
319 222
320 if (verboseLevel > 3) { 223 if (verboseLevel > 3) {
321 G4cout << "G4LowEPComptonModel::SampleSeco << 224 G4cout << "G4LowEPComptonModel::SampleSecondaries() E(MeV)= "
322 << photonEnergy0/MeV << " in " << c << 225 << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName()
323 << G4endl; << 226 << G4endl;
324 } 227 }
325 // do nothing below the threshold << 228
326 // should never get here because the XS is z << 229 // low-energy gamma is absorpted by this process
327 if (photonEnergy0 < LowEnergyLimit()) << 230 if (photonEnergy0 <= lowEnergyLimit)
328 return ; << 231 {
>> 232 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 233 fParticleChange->SetProposedKineticEnergy(0.);
>> 234 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
>> 235 return ;
>> 236 }
329 237
330 G4double e0m = photonEnergy0 / electron_mass 238 G4double e0m = photonEnergy0 / electron_mass_c2 ;
331 G4ParticleMomentum photonDirection0 = aDynam 239 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
332 240
333 // Select randomly one element in the curren 241 // Select randomly one element in the current material
334 const G4ParticleDefinition* particle = aDyn 242 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
335 const G4Element* elm = SelectRandomAtom(coup 243 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
336 G4int Z = (G4int)elm->GetZ(); 244 G4int Z = (G4int)elm->GetZ();
337 245
338 G4double LowEPCepsilon0 = 1. / (1. + 2. * e0 246 G4double LowEPCepsilon0 = 1. / (1. + 2. * e0m);
339 G4double LowEPCepsilon0Sq = LowEPCepsilon0 * 247 G4double LowEPCepsilon0Sq = LowEPCepsilon0 * LowEPCepsilon0;
340 G4double alpha1 = -std::log(LowEPCepsilon0); 248 G4double alpha1 = -std::log(LowEPCepsilon0);
341 G4double alpha2 = 0.5 * (1. - LowEPCepsilon0 249 G4double alpha2 = 0.5 * (1. - LowEPCepsilon0Sq);
342 250
343 G4double wlPhoton = h_Planck*c_light/photonE 251 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
344 252
345 // Sample the energy of the scattered photon 253 // Sample the energy of the scattered photon
346 G4double LowEPCepsilon; 254 G4double LowEPCepsilon;
347 G4double LowEPCepsilonSq; 255 G4double LowEPCepsilonSq;
348 G4double oneCosT; 256 G4double oneCosT;
349 G4double sinT2; 257 G4double sinT2;
350 G4double gReject; 258 G4double gReject;
351 << 259
352 if (verboseLevel > 3) { <<
353 G4cout << "Started loop to sample gamma en <<
354 } <<
355 <<
356 do 260 do
357 { 261 {
358 if ( alpha1/(alpha1+alpha2) > G4UniformR 262 if ( alpha1/(alpha1+alpha2) > G4UniformRand())
359 { << 263 {
360 LowEPCepsilon = G4Exp(-alpha1 * G4Un << 264 LowEPCepsilon = std::exp(-alpha1 * G4UniformRand());
361 LowEPCepsilonSq = LowEPCepsilon * Lo << 265 LowEPCepsilonSq = LowEPCepsilon * LowEPCepsilon;
362 } << 266 }
363 else 267 else
364 { << 268 {
365 LowEPCepsilonSq = LowEPCepsilon0Sq + << 269 LowEPCepsilonSq = LowEPCepsilon0Sq + (1. - LowEPCepsilon0Sq) * G4UniformRand();
366 LowEPCepsilon = std::sqrt(LowEPCepsi << 270 LowEPCepsilon = std::sqrt(LowEPCepsilonSq);
367 } << 271 }
368 272
369 oneCosT = (1. - LowEPCepsilon) / ( LowEP 273 oneCosT = (1. - LowEPCepsilon) / ( LowEPCepsilon * e0m);
370 sinT2 = oneCosT * (2. - oneCosT); 274 sinT2 = oneCosT * (2. - oneCosT);
371 G4double x = std::sqrt(oneCosT/2.) / (wl 275 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
372 G4double scatteringFunction = ComputeSca << 276 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
373 gReject = (1. - LowEPCepsilon * sinT2 / 277 gReject = (1. - LowEPCepsilon * sinT2 / (1. + LowEPCepsilonSq)) * scatteringFunction;
374 278
375 } while(gReject < G4UniformRand()*Z); << 279 } while(gReject < G4UniformRand()*Z);
376 280
377 G4double cosTheta = 1. - oneCosT; 281 G4double cosTheta = 1. - oneCosT;
378 G4double sinTheta = std::sqrt(sinT2); 282 G4double sinTheta = std::sqrt(sinT2);
379 G4double phi = twopi * G4UniformRand(); 283 G4double phi = twopi * G4UniformRand();
380 G4double dirx = sinTheta * std::cos(phi); 284 G4double dirx = sinTheta * std::cos(phi);
381 G4double diry = sinTheta * std::sin(phi); 285 G4double diry = sinTheta * std::sin(phi);
382 G4double dirz = cosTheta ; 286 G4double dirz = cosTheta ;
383 287
>> 288
384 // Scatter photon energy and Compton electro 289 // Scatter photon energy and Compton electron direction - Method based on:
385 // J. M. C. Brown, M. R. Dimmock, J. E. Gill 290 // J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin'
386 // "A low energy bound atomic electron Compt 291 // "A low energy bound atomic electron Compton scattering model for Geant4"
387 // NIMB, Vol. 338, 77-88, 2014. << 292 // NIMA ISSUE, PG, 2013
388 << 293
389 // Set constants and initialize scattering p 294 // Set constants and initialize scattering parameters
>> 295
390 const G4double vel_c = c_light / (m/s); 296 const G4double vel_c = c_light / (m/s);
391 const G4double momentum_au_to_nat = halfpi* 297 const G4double momentum_au_to_nat = halfpi* hbar_Planck / Bohr_radius / (kg*m/s);
392 const G4double e_mass_kg = electron_mass_c2 298 const G4double e_mass_kg = electron_mass_c2 / c_squared / kg ;
393 << 299
394 const G4int maxDopplerIterations = 1000; << 300 const G4int maxDopplerIterations = 1000;
395 G4double bindingE = 0.; << 301 G4double bindingE = 0.;
396 G4double pEIncident = photonEnergy0 ; 302 G4double pEIncident = photonEnergy0 ;
397 G4double pERecoil = -1.; 303 G4double pERecoil = -1.;
398 G4double eERecoil = -1.; 304 G4double eERecoil = -1.;
399 G4double e_alpha =0.; 305 G4double e_alpha =0.;
400 G4double e_beta = 0.; 306 G4double e_beta = 0.;
401 << 307
402 G4double CE_emission_flag = 0.; 308 G4double CE_emission_flag = 0.;
403 G4double ePAU = -1; 309 G4double ePAU = -1;
404 G4int shellIdx = 0; 310 G4int shellIdx = 0;
405 G4double u_temp = 0; 311 G4double u_temp = 0;
406 G4double cosPhiE =0; 312 G4double cosPhiE =0;
407 G4double sinThetaE =0; 313 G4double sinThetaE =0;
408 G4double cosThetaE =0; 314 G4double cosThetaE =0;
409 G4int iteration = 0; << 315 G4int iteration = 0;
410 <<
411 if (verboseLevel > 3) { <<
412 G4cout << "Started loop to sample photon e <<
413 } <<
414 <<
415 do{ 316 do{
>> 317
>> 318
416 // ************************************* 319 // ******************************************
417 // | Determine scatter photon energy 320 // | Determine scatter photon energy |
418 // ************************************* << 321 // ******************************************
419 do << 322
420 { << 323 do
421 iteration++; << 324 {
422 // ***************************************** << 325 iteration++;
423 // | Sample bound electron information << 326
424 // ***************************************** << 327
425 << 328 // ********************************************
426 // Select shell based on shell occupancy << 329 // | Sample bound electron information |
427 shellIdx = shellData->SelectRandomShell(Z); << 330 // ********************************************
428 bindingE = shellData->BindingEnergy(Z,shellI << 331
429 << 332 // Select shell based on shell occupancy
430 // Randomly sample bound electron momentum ( << 333
431 ePAU = profileData->RandomSelectMomentum(Z,s << 334 shellIdx = shellData.SelectRandomShell(Z);
432 << 335 bindingE = shellData.BindingEnergy(Z,shellIdx)/MeV;
433 // Convert to SI units << 336
434 G4double ePSI = ePAU * momentum_au_to_nat; << 337
435 << 338 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
436 //Calculate bound electron velocity and norm << 339 ePAU = profileData.RandomSelectMomentum(Z,shellIdx);
437 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / << 340
438 << 341 // Convert to SI units
439 // Sample incident electron direction, amorp << 342 G4double ePSI = ePAU * momentum_au_to_nat;
440 e_alpha = pi*G4UniformRand(); << 343
441 e_beta = twopi*G4UniformRand(); << 344 //Calculate bound electron velocity and normalise to natural units
442 << 345 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)) )/vel_c;
443 // Total energy of system << 346
444 << 347 // Sample incident electron direction, amorphous material, to scattering photon scattering plane
445 G4double eEIncident = electron_mass_c2 / sqr << 348
446 G4double systemE = eEIncident + pEIncident; << 349 e_alpha = pi*G4UniformRand();
447 G4double gamma_temp = 1.0 / sqrt( 1 - (u_tem << 350 e_beta = twopi*G4UniformRand();
448 G4double numerator = gamma_temp*electron_mas << 351
449 G4double subdenom1 = u_temp*cosTheta*std::c << 352 // Total energy of system
450 G4double subdenom2 = u_temp*sinTheta*std::si << 353
451 G4double denominator = (1.0 - cosTheta) + ( << 354 G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
452 pERecoil = (numerator/denominator); << 355 G4double systemE = eEIncident + pEIncident;
453 eERecoil = systemE - pERecoil; << 356
454 CE_emission_flag = pEIncident - pERecoil; << 357
455 } while ( (iteration <= maxDopplerIterat << 358 G4double gamma_temp = 1.0 / sqrt( 1 - (u_temp*u_temp));
456 << 359 G4double numerator = gamma_temp*electron_mass_c2*(1 - u_temp * std::cos(e_alpha));
457 // End of recalculation of photon energy w << 360 G4double subdenom1 = u_temp*cosTheta*std::cos(e_alpha);
458 << 361 G4double subdenom2 = u_temp*sinTheta*std::sin(e_alpha)*std::cos(e_beta);
459 // *************************************** << 362 G4double denominator = (1.0 - cosTheta) + (gamma_temp*electron_mass_c2*(1 - subdenom1 - subdenom2) / pEIncident);
460 // | Determine ejected Compton electro << 363 pERecoil = (numerator/denominator);
461 // *************************************** << 364 eERecoil = systemE - pERecoil;
462 << 365 CE_emission_flag = pEIncident - pERecoil;
463 // Calculate velocity of ejected Compton e << 366 } while ( (iteration <= maxDopplerIterations) && (CE_emission_flag < bindingE));
464 << 367
465 G4double a_temp = eERecoil / electron_mass << 368
466 G4double u_p_temp = sqrt(1 - (1 / (a_temp* << 369
467 << 370 // End of recalculation of photon energy with Doppler broadening
468 // Coefficients and terms from simulatenou << 371
469 G4double sinAlpha = std::sin(e_alpha); << 372
470 G4double cosAlpha = std::cos(e_alpha); << 373
471 G4double sinBeta = std::sin(e_beta); << 374 // *******************************************************
472 G4double cosBeta = std::cos(e_beta); << 375 // | Determine ejected Compton electron direction |
473 << 376 // *******************************************************
474 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_ << 377
475 G4double gamma_p = 1.0 / sqrt(1 - (u_p_tem << 378 // Calculate velocity of ejected Compton electron
476 << 379
477 G4double var_A = pERecoil*u_p_temp*sinThet << 380 G4double a_temp = eERecoil / electron_mass_c2;
478 G4double var_B = u_p_temp* (pERecoil*cosTh << 381 G4double u_p_temp = sqrt(1 - (1 / (a_temp*a_temp)));
479 G4double var_C = (pERecoil-pEIncident) - ( << 382
480 << 383 // Coefficients and terms from simulatenous equations
481 G4double var_D1 = gamma*electron_mass_c2*p << 384
482 G4double var_D2 = (1 - (u_temp*cosTheta*co << 385 G4double sinAlpha = std::sin(e_alpha);
483 G4double var_D3 = ((electron_mass_c2*elect << 386 G4double cosAlpha = std::cos(e_alpha);
484 G4double var_D = var_D1*var_D2 + var_D3; << 387 G4double sinBeta = std::sin(e_beta);
485 << 388 G4double cosBeta = std::cos(e_beta);
486 G4double var_E1 = ((gamma*gamma_p)*(electr << 389
487 G4double var_E2 = gamma_p*electron_mass_c2 << 390 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_temp));
488 G4double var_E = var_E1 - var_E2; << 391 G4double gamma_p = 1.0 / sqrt(1 - (u_p_temp*u_p_temp));
489 << 392
490 G4double var_F1 = ((gamma*gamma_p)*(electr << 393 G4double var_A = pERecoil*u_p_temp*sinTheta;
491 G4double var_F2 = (gamma_p*electron_mass_c << 394 G4double var_B = u_p_temp* (pERecoil*cosTheta-pEIncident);
492 G4double var_F = var_F1 - var_F2; << 395 G4double var_C = (pERecoil-pEIncident) - ( (pERecoil*pEIncident) / (gamma_p*electron_mass_c2))*(1 - cosTheta);
493 << 396
494 G4double var_G = (gamma*gamma_p)*(electron << 397 G4double var_D1 = gamma*electron_mass_c2*pERecoil;
495 << 398 G4double var_D2 = (1 - (u_temp*cosTheta*cosAlpha) - (u_temp*sinTheta*cosBeta*sinAlpha));
496 // Two equations form a quadratic form of << 399 G4double var_D3 = ((electron_mass_c2*electron_mass_c2)*(gamma*gamma_p - 1)) - (gamma_p*electron_mass_c2*pERecoil);
497 // Coefficents and solution to quadratic << 400 G4double var_D = var_D1*var_D2 + var_D3;
498 G4double var_W1 = (var_F*var_B - var_E*var << 401
499 G4double var_W2 = (var_G*var_G)*(var_A*var << 402 G4double var_E1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosAlpha);
500 G4double var_W = var_W1 + var_W2; << 403 G4double var_E2 = gamma_p*electron_mass_c2*pERecoil*u_p_temp*cosTheta;
501 << 404 G4double var_E = var_E1 - var_E2;
502 G4double var_Y = 2.0*(((var_A*var_D-var_F* << 405
503 << 406 G4double var_F1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosBeta*sinAlpha);
504 G4double var_Z1 = (var_A*var_D - var_F*var << 407 G4double var_F2 = (gamma_p*electron_mass_c2*pERecoil*u_p_temp*sinTheta);
505 G4double var_Z2 = (var_G*var_G)*(var_C*var << 408 G4double var_F = var_F1 - var_F2;
506 G4double var_Z = var_Z1 + var_Z2; << 409
507 G4double diff1 = var_Y*var_Y; << 410 G4double var_G = (gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*sinBeta*sinAlpha;
508 G4double diff2 = 4*var_W*var_Z; << 411
509 G4double diff = diff1 - diff2; << 412 // Two equations form a quadratic form of Wx^2 + Yx + Z = 0
510 << 413 // Coefficents and solution to quadratic
511 // Check if diff is less than zero, if so << 414
512 //Determine number of digits (in decimal b << 415 G4double var_W1 = (var_F*var_B - var_E*var_A)*(var_F*var_B - var_E*var_A);
513 G4double g4d_order = G4double(numeric_limi << 416 G4double var_W2 = (var_G*var_G)*(var_A*var_A) + (var_G*var_G)*(var_B*var_B);
514 G4double g4d_limit = std::pow(10.,-g4d_ord << 417 G4double var_W = var_W1 + var_W2;
515 //Confirm that diff less than zero is due << 418
516 //than 10^(-g4d_order), then set diff to z << 419 G4double var_Y = 2.0*(((var_A*var_D-var_F*var_C)*(var_F*var_B-var_E*var_A)) - ((var_G*var_G)*var_B*var_C));
517 << 420
518 if ((diff < 0.0) && (abs(diff / diff1) < g << 421 G4double var_Z1 = (var_A*var_D - var_F*var_C)*(var_A*var_D - var_F*var_C);
519 { << 422 G4double var_Z2 = (var_G*var_G)*(var_C*var_C) - (var_G*var_G)*(var_A*var_A);
520 diff = 0.0; << 423 G4double var_Z = var_Z1 + var_Z2;
521 } << 424 G4double diff1 = var_Y*var_Y;
522 << 425 G4double diff2 = 4*var_W*var_Z;
523 << 426 G4double diff = diff1 - diff2;
524 // Plus and minus of quadratic << 427
525 G4double X_p = (-var_Y + sqrt (diff))/(2*v << 428
526 G4double X_m = (-var_Y - sqrt (diff))/(2*v << 429 // Check if diff is less than zero, if so ensure it is due to FPE
527 << 430
528 // Floating point precision protection << 431 //Determine number of digits (in decimal base) that G4double can accurately represent
529 // Check if X_p and X_m are greater than o << 432 G4double g4d_order = G4double(numeric_limits<G4double>::digits10);
530 // Issue due to propagation of FPE and onl << 433 G4double g4d_limit = std::pow(10.,-g4d_order);
531 if(X_p >1){X_p=1;} if(X_p<-1){X_p=-1;} << 434 //Confirm that diff less than zero is due FPE, i.e if abs of diff / diff1 and diff/ diff2 is less
532 if(X_m >1){X_m=1;} if(X_m<-1){X_m=-1;} << 435 //than 10^(-g4d_order), then set diff to zero
533 // End of FP protection << 436
534 << 437 if ((diff < 0.0) && (abs(diff / diff1) < g4d_limit) && (abs(diff / diff2) < g4d_limit) )
535 G4double ThetaE = 0.; << 438 {
536 // Randomly sample one of the two possible << 439 diff = 0.0;
537 G4double sol_select = G4UniformRand(); << 440 }
538 << 441
539 if (sol_select < 0.5) << 442
>> 443 // Plus and minus of quadratic
>> 444 G4double X_p = (-var_Y + sqrt (diff))/(2*var_W);
>> 445 G4double X_m = (-var_Y - sqrt (diff))/(2*var_W);
>> 446
>> 447
>> 448 // Randomly sample one of the two possible solutions and determin theta angle of ejected Compton electron
>> 449 G4double ThetaE = 0.;
>> 450 G4double sol_select = G4UniformRand();
>> 451
>> 452 if (sol_select < 0.5)
540 { 453 {
541 ThetaE = std::acos(X_p); << 454 ThetaE = std::acos(X_p);
542 } 455 }
543 if (sol_select > 0.5) << 456 if (sol_select > 0.5)
544 { 457 {
545 ThetaE = std::acos(X_m); << 458 ThetaE = std::acos(X_m);
546 } 459 }
547 cosThetaE = std::cos(ThetaE); << 460
548 sinThetaE = std::sin(ThetaE); << 461 cosThetaE = std::cos(ThetaE);
549 G4double Theta = std::acos(cosTheta); << 462 sinThetaE = std::sin(ThetaE);
550 << 463 G4double Theta = std::acos(cosTheta);
551 //Calculate electron Phi << 464
552 G4double iSinThetaE = std::sqrt(1+std::tan << 465 //Calculate electron Phi
553 G4double iSinTheta = std::sqrt(1+std::tan( << 466 G4double iSinThetaE = std::sqrt(1+std::tan((pi/2.0)-ThetaE)*std::tan((pi/2.0)-ThetaE));
554 G4double ivar_A = iSinTheta/ (pERecoil*u_p << 467 G4double iSinTheta = std::sqrt(1+std::tan((pi/2.0)-Theta)*std::tan((pi/2.0)-Theta));
555 // Trigs << 468 G4double ivar_A = iSinTheta/ (pERecoil*u_p_temp);
556 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_ << 469 // Trigs
557 << 470 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_A*iSinThetaE);
558 // End of calculation of ejection Compton << 471
559 //Fix for floating point errors << 472 // End of calculation of ejection Compton electron direction
560 << 473
561 } while ( (iteration <= maxDopplerIterations << 474 //Fix for floating point errors
>> 475
>> 476 } while ( (iteration <= maxDopplerIterations) && (abs(cosPhiE) > 1));
>> 477
>> 478 // Revert to original if maximum number of iterations threshold has been reached
>> 479
562 480
563 // Revert to original if maximum number of i <<
564 if (iteration >= maxDopplerIterations) 481 if (iteration >= maxDopplerIterations)
565 { 482 {
566 pERecoil = photonEnergy0 ; 483 pERecoil = photonEnergy0 ;
567 bindingE = 0.; 484 bindingE = 0.;
568 dirx=0.0; 485 dirx=0.0;
569 diry=0.0; 486 diry=0.0;
570 dirz=1.0; 487 dirz=1.0;
571 } 488 }
572 << 489
573 // Set "scattered" photon direction and ener 490 // Set "scattered" photon direction and energy
>> 491
574 G4ThreeVector photonDirection1(dirx,diry,dir 492 G4ThreeVector photonDirection1(dirx,diry,dirz);
575 photonDirection1.rotateUz(photonDirection0); 493 photonDirection1.rotateUz(photonDirection0);
576 fParticleChange->ProposeMomentumDirection(ph 494 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
577 495
578 if (pERecoil > 0.) 496 if (pERecoil > 0.)
579 { 497 {
580 fParticleChange->SetProposedKineticEnergy 498 fParticleChange->SetProposedKineticEnergy(pERecoil) ;
581 499
582 // Set ejected Compton electron direction 500 // Set ejected Compton electron direction and energy
583 G4double PhiE = std::acos(cosPhiE); 501 G4double PhiE = std::acos(cosPhiE);
584 G4double eDirX = sinThetaE * std::cos(phi 502 G4double eDirX = sinThetaE * std::cos(phi+PhiE);
585 G4double eDirY = sinThetaE * std::sin(phi 503 G4double eDirY = sinThetaE * std::sin(phi+PhiE);
586 G4double eDirZ = cosThetaE; 504 G4double eDirZ = cosThetaE;
587 << 505
588 G4double eKineticEnergy = pEIncident - pE << 506 G4double eKineticEnergy = pEIncident - pERecoil - bindingE;
589 << 507
590 G4ThreeVector eDirection(eDirX,eDirY,eDir 508 G4ThreeVector eDirection(eDirX,eDirY,eDirZ);
591 eDirection.rotateUz(photonDirection0); 509 eDirection.rotateUz(photonDirection0);
592 G4DynamicParticle* dp = new G4DynamicPart 510 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
593 << 511 eDirection,eKineticEnergy) ;
594 fvect->push_back(dp); 512 fvect->push_back(dp);
595 513
596 } 514 }
597 else 515 else
598 { 516 {
599 fParticleChange->SetProposedKineticEnerg 517 fParticleChange->SetProposedKineticEnergy(0.);
600 fParticleChange->ProposeTrackStatus(fSto << 518 fParticleChange->ProposeTrackStatus(fStopAndKill);
601 } 519 }
>> 520
>> 521
602 522
603 // sample deexcitation <<
604 // <<
605 if (verboseLevel > 3) { <<
606 G4cout << "Started atomic de-excitation " <<
607 } <<
608 523
>> 524 // sample deexcitation
>> 525
609 if(fAtomDeexcitation && iteration < maxDoppl 526 if(fAtomDeexcitation && iteration < maxDopplerIterations) {
610 G4int index = couple->GetIndex(); 527 G4int index = couple->GetIndex();
611 if(fAtomDeexcitation->CheckDeexcitationAct 528 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
612 std::size_t nbefore = fvect->size(); << 529 size_t nbefore = fvect->size();
613 G4AtomicShellEnumerator as = G4AtomicShe 530 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
614 const G4AtomicShell* shell = fAtomDeexci 531 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
615 fAtomDeexcitation->GenerateParticles(fve 532 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
616 std::size_t nafter = fvect->size(); << 533 size_t nafter = fvect->size();
617 if(nafter > nbefore) { 534 if(nafter > nbefore) {
618 for (std::size_t i=nbefore; i<nafter; << 535 for (size_t i=nbefore; i<nafter; ++i) {
619 //Check if there is enough residual << 536 //Check if there is enough residual energy
620 if (bindingE >= ((*fvect)[i])->GetKi 537 if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
621 { 538 {
622 //Ok, this is a valid secondary: 539 //Ok, this is a valid secondary: keep it
623 bindingE -= ((*fvect)[i])->GetKin 540 bindingE -= ((*fvect)[i])->GetKineticEnergy();
624 } 541 }
625 else 542 else
626 { 543 {
627 //Invalid secondary: not enough e 544 //Invalid secondary: not enough energy to create it!
628 //Keep its energy in the local de 545 //Keep its energy in the local deposit
629 delete (*fvect)[i]; 546 delete (*fvect)[i];
630 (*fvect)[i]=nullptr; << 547 (*fvect)[i]=0;
631 } 548 }
632 } << 549 }
633 } 550 }
634 } 551 }
635 } 552 }
636 553
637 //This should never happen 554 //This should never happen
638 if(bindingE < 0.0) 555 if(bindingE < 0.0)
639 G4Exception("G4LowEPComptonModel::SampleS 556 G4Exception("G4LowEPComptonModel::SampleSecondaries()",
640 "em2051",FatalException,"Nega 557 "em2051",FatalException,"Negative local energy deposit");
641 558
642 fParticleChange->ProposeLocalEnergyDeposit(b 559 fParticleChange->ProposeLocalEnergyDeposit(bindingE);
>> 560
643 } 561 }
644 <<
645 //******************************************** <<
646 <<
647 G4double <<
648 G4LowEPComptonModel::ComputeScatteringFunction <<
649 { <<
650 G4double value = Z; <<
651 if (x <= ScatFuncFitParam[Z][2]) { <<
652 <<
653 G4double lgq = G4Log(x)/ln10; <<
654 <<
655 if (lgq < ScatFuncFitParam[Z][1]) { <<
656 value = ScatFuncFitParam[Z][3] + lgq*Sca <<
657 } else { <<
658 value = ScatFuncFitParam[Z][5] + lgq*Sca <<
659 lgq*lgq*ScatFuncFitParam[Z][7] + lgq*l <<
660 } <<
661 value = G4Exp(value*ln10); <<
662 } <<
663 return value; <<
664 } <<
665 <<
666 <<
667 //******************************************** <<
668 <<
669 void <<
670 G4LowEPComptonModel::InitialiseForElement(cons <<
671 <<
672 { <<
673 G4AutoLock l(&LowEPComptonModelMutex); <<
674 if(!data[Z]) { ReadData(Z); } <<
675 l.unlock(); <<
676 } <<
677 <<
678 //******************************************** <<
679 <<
680 //Fitting data to compute scattering function <<
681 <<
682 const G4double G4LowEPComptonModel::ScatFuncFi <<
683 { 0, 0., 0., 0., 0., <<
684 { 1, 6.673, 1.49968E+08, -14.352, 1.999, -143 <<
685 { 2, 6.500, 2.50035E+08, -14.215, 1.970, -53. <<
686 { 3, 6.551, 3.99945E+08, -13.555, 1.993, -62. <<
687 { 4, 6.500, 5.00035E+08, -13.746, 1.998, -127 <<
688 { 5, 6.500, 5.99791E+08, -13.800, 1.998, -131 <<
689 { 6, 6.708, 6.99842E+08, -13.885, 1.999, -128 <<
690 { 7, 6.685, 7.99834E+08, -13.885, 2.000, -131 <<
691 { 8, 6.669, 7.99834E+08, -13.962, 2.001, -128 <<
692 { 9, 6.711, 7.99834E+08, -13.999, 2.000, -122 <<
693 { 10, 6.702, 7.99834E+08, -14.044, 1.999, -110 <<
694 { 11, 6.425, 1.00000E+09, -13.423, 1.993, -41. <<
695 { 12, 6.542, 1.00000E+09, -13.389, 1.997, -53. <<
696 { 13, 6.570, 1.49968E+09, -13.401, 1.997, -66. <<
697 { 14, 6.364, 1.49968E+09, -13.452, 1.999, -78. <<
698 { 15, 6.500, 1.49968E+09, -13.488, 1.998, -85. <<
699 { 16, 6.500, 1.49968E+09, -13.532, 1.998, -93. <<
700 { 17, 6.500, 1.49968E+09, -13.584, 2.000, -98. <<
701 { 18, 6.500, 1.49968E+09, -13.618, 1.999, -100 <<
702 { 19, 6.500, 1.99986E+09, -13.185, 1.992, -53. <<
703 { 20, 6.490, 1.99986E+09, -13.123, 1.993, -52. <<
704 { 21, 6.498, 1.99986E+09, -13.157, 1.994, -55. <<
705 { 22, 6.495, 1.99986E+09, -13.183, 1.994, -57. <<
706 { 23, 6.487, 1.99986E+09, -13.216, 1.995, -58. <<
707 { 24, 6.500, 1.99986E+09, -13.330, 1.997, -62. <<
708 { 25, 6.488, 1.99986E+09, -13.277, 1.997, -58. <<
709 { 26, 6.500, 5.00035E+09, -13.292, 1.997, -61. <<
710 { 27, 6.500, 5.00035E+09, -13.321, 1.998, -61. <<
711 { 28, 6.500, 5.00035E+09, -13.340, 1.998, -62. <<
712 { 29, 6.500, 5.00035E+09, -13.439, 1.998, -67. <<
713 { 30, 6.500, 5.00035E+09, -13.383, 1.999, -62. <<
714 { 31, 6.500, 5.00035E+09, -13.349, 1.997, -61. <<
715 { 32, 6.500, 5.00035E+09, -13.373, 1.999, -63. <<
716 { 33, 6.500, 5.00035E+09, -13.395, 1.999, -65. <<
717 { 34, 6.500, 5.00035E+09, -13.417, 1.999, -69. <<
718 { 35, 6.500, 5.00035E+09, -13.442, 2.000, -72. <<
719 { 36, 6.500, 5.00035E+09, -13.451, 1.998, -74. <<
720 { 37, 6.500, 5.00035E+09, -13.082, 1.991, -46. <<
721 { 38, 6.465, 5.00035E+09, -13.022, 1.993, -41. <<
722 { 39, 6.492, 5.00035E+09, -13.043, 1.994, -44. <<
723 { 40, 6.499, 5.00035E+09, -13.064, 1.994, -47. <<
724 { 41, 6.384, 5.00035E+09, -13.156, 1.996, -53. <<
725 { 42, 6.500, 5.00035E+09, -13.176, 1.996, -54. <<
726 { 43, 6.500, 5.00035E+09, -13.133, 1.997, -51. <<
727 { 44, 6.500, 5.00035E+09, -13.220, 1.996, -58. <<
728 { 45, 6.500, 5.00035E+09, -13.246, 1.998, -59. <<
729 { 46, 6.500, 5.00035E+09, -13.407, 1.999, -72. <<
730 { 47, 6.500, 5.00035E+09, -13.277, 1.998, -60. <<
731 { 48, 6.500, 5.00035E+09, -13.222, 1.998, -56. <<
732 { 49, 6.500, 5.00035E+09, -13.199, 1.997, -56. <<
733 { 50, 6.500, 5.00035E+09, -13.215, 1.998, -58. <<
734 { 51, 6.500, 5.00035E+09, -13.230, 1.998, -60. <<
735 { 52, 6.500, 7.99834E+09, -13.246, 1.998, -63. <<
736 { 53, 6.500, 5.00035E+09, -13.262, 1.998, -66. <<
737 { 54, 6.500, 7.99834E+09, -13.279, 1.998, -67. <<
738 { 55, 6.500, 5.00035E+09, -12.951, 1.990, -45. <<
739 { 56, 6.425, 5.00035E+09, -12.882, 1.992, -39. <<
740 { 57, 6.466, 2.82488E+09, -12.903, 1.992, -38. <<
741 { 58, 6.451, 5.00035E+09, -12.915, 1.993, -41. <<
742 { 59, 6.434, 5.00035E+09, -12.914, 1.993, -40. <<
743 { 60, 6.444, 5.00035E+09, -12.922, 1.992, -39. <<
744 { 61, 6.414, 7.99834E+09, -12.930, 1.993, -42. <<
745 { 62, 6.420, 7.99834E+09, -12.938, 1.992, -42. <<
746 { 63, 6.416, 7.99834E+09, -12.946, 1.993, -42. <<
747 { 64, 6.443, 7.99834E+09, -12.963, 1.993, -43. <<
748 { 65, 6.449, 7.99834E+09, -12.973, 1.993, -43. <<
749 { 66, 6.419, 7.99834E+09, -12.966, 1.993, -42. <<
750 { 67, 6.406, 7.99834E+09, -12.976, 1.993, -42. <<
751 { 68, 6.424, 7.99834E+09, -12.986, 1.993, -41. <<
752 { 69, 6.417, 7.99834E+09, -12.989, 1.993, -42. <<
753 { 70, 6.405, 7.99834E+09, -13.000, 1.994, -42. <<
754 { 71, 6.449, 7.99834E+09, -13.015, 1.994, -42. <<
755 { 72, 6.465, 7.99834E+09, -13.030, 1.994, -43. <<
756 { 73, 6.447, 7.99834E+09, -13.048, 1.996, -44. <<
757 { 74, 6.452, 7.99834E+09, -13.073, 1.997, -45. <<
758 { 75, 6.432, 7.99834E+09, -13.082, 1.997, -46. <<
759 { 76, 6.439, 7.99834E+09, -13.100, 1.997, -47. <<
760 { 77, 6.432, 7.99834E+09, -13.110, 1.997, -48. <<
761 { 78, 6.500, 7.99834E+09, -13.185, 1.997, -53. <<
762 { 79, 6.500, 7.99834E+09, -13.200, 1.997, -53. <<
763 { 80, 6.500, 7.99834E+09, -13.156, 1.998, -49. <<
764 { 81, 6.500, 7.99834E+09, -13.128, 1.997, -49. <<
765 { 82, 6.500, 7.99834E+09, -13.134, 1.997, -51. <<
766 { 83, 6.500, 7.99834E+09, -13.148, 1.998, -52. <<
767 { 84, 6.500, 7.99834E+09, -13.161, 1.998, -54. <<
768 { 85, 6.500, 7.99834E+09, -13.175, 1.998, -56. <<
769 { 86, 6.500, 7.99834E+09, -13.189, 1.998, -57. <<
770 { 87, 6.500, 7.99834E+09, -12.885, 1.990, -39. <<
771 { 88, 6.417, 7.99834E+09, -12.816, 1.991, -34. <<
772 { 89, 6.442, 7.99834E+09, -12.831, 1.992, -36. <<
773 { 90, 6.463, 7.99834E+09, -12.850, 1.993, -37. <<
774 { 91, 6.447, 7.99834E+09, -12.852, 1.993, -37. <<
775 { 92, 6.439, 7.99834E+09, -12.858, 1.993, -37. <<
776 { 93, 6.437, 1.00000E+10, -12.866, 1.993, -39. <<
777 { 94, 6.432, 7.99834E+09, -12.862, 1.993, -37. <<
778 { 95, 6.435, 7.99834E+09, -12.869, 1.993, -37. <<
779 { 96, 6.449, 1.00000E+10, -12.886, 1.993, -39. <<
780 { 97, 6.446, 1.00000E+10, -12.892, 1.993, -40. <<
781 { 98, 6.421, 1.00000E+10, -12.887, 1.993, -39. <<
782 { 99, 6.414, 1.00000E+10, -12.894, 1.993, -39. <<
783 {100, 6.412, 1.00000E+10, -12.900, 1.993, -39. <<
784 }; <<
785 <<
786 <<
787 562
788 563