| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4ICRU73QOModel
33 //
34 // Author: Alexander Bagulya
35 //
36 // Creation date: 21.05.2010
37 //
38 // Modifications:
39 //
40 //
41 // -------------------------------------------
42 //
43 //....oooOO0OOooo........oooOO0OOooo........oo
44 //....oooOO0OOooo........oooOO0OOooo........oo
45
46 #include "G4ICRU73QOModel.hh"
47 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "Randomize.hh"
50 #include "G4Electron.hh"
51 #include "G4ParticleChangeForLoss.hh"
52 #include "G4LossTableManager.hh"
53 #include "G4AntiProton.hh"
54 #include "G4DeltaAngle.hh"
55 #include "G4Log.hh"
56 #include "G4Exp.hh"
57
58 //....oooOO0OOooo........oooOO0OOooo........oo
59
60 using namespace std;
61
62 //....oooOO0OOooo........oooOO0OOooo........oo
63
64 G4ICRU73QOModel::G4ICRU73QOModel(const G4Parti
65 : G4VEmModel(nam),
66 particle(nullptr),
67 isInitialised(false)
68 {
69 mass = charge = chargeSquare = massRate = ra
70 if(p) { SetParticle(p); }
71 SetHighEnergyLimit(10.0*MeV);
72
73 lowestKinEnergy = 5.0*keV;
74
75 sizeL0 = 67;
76 sizeL1 = 22;
77 sizeL2 = 14;
78
79 theElectron = G4Electron::Electron();
80
81 for (G4int i = 0; i < 100; ++i)
82 {
83 indexZ[i] = -1;
84 }
85 for(G4int i = 0; i < NQOELEM; ++i)
86 {
87 if(ZElementAvailable[i] > 0) {
88 indexZ[ZElementAvailable[i]] = i;
89 }
90 }
91 fParticleChange = nullptr;
92 denEffData = nullptr;
93 }
94
95 //....oooOO0OOooo........oooOO0OOooo........oo
96
97 void G4ICRU73QOModel::Initialise(const G4Parti
98 const G4DataV
99 {
100 if(p != particle) SetParticle(p);
101
102 // always false before the run
103 SetDeexcitationFlag(false);
104
105 if(!isInitialised) {
106 isInitialised = true;
107
108 if(UseAngularGeneratorFlag() && !GetAngula
109 SetAngularDistribution(new G4DeltaAngle(
110 }
111
112 G4String pname = particle->GetParticleName
113 fParticleChange = GetParticleChangeForLoss
114 const G4MaterialTable* mtab = G4Material::
115 denEffData = (*mtab)[0]->GetIonisation()->
116 }
117 }
118
119 //....oooOO0OOooo........oooOO0OOooo........oo
120
121 G4double G4ICRU73QOModel::ComputeCrossSectionP
122 con
123
124
125
126 {
127 G4double cross = 0.0;
128 const G4double tmax = MaxSecondaryEnerg
129 const G4double maxEnergy = std::min(tmax, ma
130 const G4double cutEnergy = std::max(cut, low
131 if(cutEnergy < maxEnergy) {
132
133 const G4double energy = kineticEnergy + m
134 const G4double energy2 = energy*energy;
135 const G4double beta2 = kineticEnergy*(ki
136 cross = (maxEnergy - cutEnergy)/(cutEnergy
137 - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
138
139 cross *= CLHEP::twopi_mc2_rcl2*chargeSquar
140 }
141
142 return cross;
143 }
144
145 //....oooOO0OOooo........oooOO0OOooo........oo
146
147 G4double G4ICRU73QOModel::ComputeCrossSectionP
148 con
149
150
151
152
153 {
154 G4double cross = Z*ComputeCrossSectionPerEle
155 (p,ki
156 return cross;
157 }
158
159 //....oooOO0OOooo........oooOO0OOooo........oo
160
161 G4double G4ICRU73QOModel::CrossSectionPerVolum
162 con
163 con
164
165
166
167 {
168 G4double eDensity = material->GetElectronDen
169 G4double cross = eDensity*ComputeCrossSectio
170 (p,ki
171 return cross;
172 }
173
174 //....oooOO0OOooo........oooOO0OOooo........oo
175
176 G4double G4ICRU73QOModel::ComputeDEDXPerVolume
177
178
179
180 {
181 SetParticle(p);
182 const G4double tmax = MaxSecondaryEnergy(p,
183 const G4double tkin = kineticEnergy/massRat
184 const G4double cutEnergy = std::max(cut, low
185 G4double dedx = 0.0;
186 if(tkin > lowestKinEnergy) { dedx = DEDX(mat
187 else { dedx = DEDX(material, lowestKinEnergy
188
189 if (cutEnergy < tmax) {
190
191 const G4double tau = kineticEnergy/mass;
192 const G4double x = cutEnergy/tmax;
193 dedx += (G4Log(x)*(tau + 1.)*(tau + 1.)/(t
194 CLHEP::twopi_mc2_rcl2 *chargeSquare * ma
195 }
196 dedx = std::max(dedx, 0.0);
197 return dedx;
198 }
199
200 //....oooOO0OOooo........oooOO0OOooo........oo
201
202 G4double G4ICRU73QOModel::DEDX(const G4Materia
203 G4double kineti
204 {
205 G4double eloss = 0.0;
206 const std::size_t numberOfElements = materia
207 const G4double* theAtomicNumDensityVector =
208 material->Get
209
210 // Bragg's rule calculation
211 const G4ElementVector* theElementVector =
212 material->GetElemen
213
214 // loop for the elements in the material
215 for (std::size_t i=0; i<numberOfElements; ++
216 {
217 const G4Element* element = (*theElementV
218 eloss += DEDXPerElement(element->GetZasI
219 * theAtomicNumDensityVector[i] * eleme
220 }
221 return eloss;
222 }
223
224 //....oooOO0OOooo........oooOO0OOooo........oo
225
226 G4double G4ICRU73QOModel::DEDXPerElement(G4int
227 G4dou
228 {
229 G4int Z = std::min(AtomicNumber, 97);
230 G4int nbOfShells = std::max(GetNumberOfShell
231
232 G4double v = CLHEP::c_light * std::sqrt( 2.0
233
234 G4double fBetheVelocity = CLHEP::fine_struct
235
236 G4double tau = kineticEnergy/proton_mass_c
237 G4double gam = tau + 1.0;
238 G4double bg2 = tau * (tau+2.0);
239 G4double beta2 = bg2/(gam*gam);
240
241 G4double l0Term = 0, l1Term = 0, l2Term = 0;
242
243 for (G4int nos = 0; nos < nbOfShells; ++nos)
244
245 G4double NormalizedEnergy = (2.0*CLHEP::el
246 GetShellEnergy(Z,nos);
247
248 G4double shStrength = GetShellStrength(Z,n
249
250 G4double l0 = GetL0(NormalizedEnergy);
251 l0Term += shStrength * l0;
252
253 G4double l1 = GetL1(NormalizedEnergy);
254 l1Term += shStrength * l1;
255
256 G4double l2 = GetL2(NormalizedEnergy);
257 l2Term += shStrength * l2;
258
259 }
260 G4double dedx = 2*CLHEP::twopi_mc2_rcl2*cha
261 (l0Term + charge*fBetheVelocity*l1Term
262 + chargeSquare*fBetheVelocity*fBetheVeloc
263 return dedx;
264 }
265
266
267 //....oooOO0OOooo........oooOO0OOooo........oo
268
269 G4double G4ICRU73QOModel::GetOscillatorEnergy(
270
271 {
272 G4int idx = denEffData->GetElementIndex(Z, k
273 if(idx == -1) { idx = denEffData->GetElement
274 G4double PlasmaEnergy = denEffData->GetPlasm
275
276 G4double PlasmaEnergy2 = PlasmaEnergy * Plas
277
278 G4double plasmonTerm = 0.66667
279 * G4AtomicShells::GetNumberOfElectrons(Z,n
280 * PlasmaEnergy2 / (Z*Z) ;
281
282 static const G4double exphalf = G4Exp(0.5);
283 G4double ionTerm = exphalf *
284 (G4AtomicShells::GetBindingEnergy(Z,nbOfTh
285 G4double ionTerm2 = ionTerm*ionTerm ;
286
287 G4double oscShellEnergy = std::sqrt( ionTerm
288
289 return oscShellEnergy;
290 }
291
292 //....oooOO0OOooo........oooOO0OOooo........oo
293
294 G4int G4ICRU73QOModel::GetNumberOfShells(G4int
295 {
296 G4int nShell = 0;
297
298 if(indexZ[Z] >= 0) {
299 nShell = nbofShellsForElement[indexZ[Z]];
300 } else {
301 nShell = G4AtomicShells::GetNumberOfShells
302 }
303 return nShell;
304 }
305
306 //....oooOO0OOooo........oooOO0OOooo........oo
307
308 G4double G4ICRU73QOModel::GetShellEnergy(G4int
309 {
310 G4double shellEnergy = 0.;
311
312 G4int idx = indexZ[Z];
313
314 if(idx >= 0) {
315 shellEnergy = ShellEnergy[startElemIndex[i
316 } else {
317 shellEnergy = GetOscillatorEnergy(Z, nbOfT
318 }
319
320 return shellEnergy;
321 }
322
323 //....oooOO0OOooo........oooOO0OOooo........oo
324
325 G4double G4ICRU73QOModel::GetShellStrength(G4i
326 {
327 G4double shellStrength = 0.;
328
329 G4int idx = indexZ[Z];
330
331 if(idx >= 0) {
332 shellStrength = SubShellOccupation[startEl
333 } else {
334 shellStrength = G4double(G4AtomicShells::G
335 }
336
337 return shellStrength;
338 }
339
340 //....oooOO0OOooo........oooOO0OOooo........oo
341
342 G4double G4ICRU73QOModel::GetL0(G4double normE
343 {
344 G4int n;
345
346 for(n = 0; n < sizeL0; n++) {
347 if( normEnergy < L0[n][0] ) break;
348 }
349 if(0 == n) { n = 1; }
350 if(n >= sizeL0) { n = sizeL0 - 1; }
351
352 G4double l0 = L0[n][1];
353 G4double l0p = L0[n-1][1];
354 G4double bethe = l0p + (l0 - l0p) * ( normEn
355 (L0[n][0] - L0[n-1][0]);
356
357 return bethe ;
358 }
359
360 //....oooOO0OOooo........oooOO0OOooo........oo
361
362 G4double G4ICRU73QOModel::GetL1(G4double normE
363 {
364 G4int n;
365
366 for(n = 0; n < sizeL1; n++) {
367 if( normEnergy < L1[n][0] ) break;
368 }
369 if(0 == n) n = 1 ;
370 if(n >= sizeL1) n = sizeL1 - 1 ;
371
372 G4double l1 = L1[n][1];
373 G4double l1p = L1[n-1][1];
374 G4double barkas= l1p + (l1 - l1p) * ( normEn
375 (L1[n][0] - L1[n-1][0]);
376
377 return barkas;
378 }
379
380 //....oooOO0OOooo........oooOO0OOooo........oo
381
382 G4double G4ICRU73QOModel::GetL2(G4double normE
383 {
384 G4int n;
385 for(n = 0; n < sizeL2; n++) {
386 if( normEnergy < L2[n][0] ) break;
387 }
388 if(0 == n) n = 1 ;
389 if(n >= sizeL2) n = sizeL2 - 1 ;
390
391 G4double l2 = L2[n][1];
392 G4double l2p = L2[n-1][1];
393 G4double bloch = l2p + (l2 - l2p) * ( normEn
394 (L2[n][0] - L2[n-1][0]);
395
396 return bloch;
397 }
398
399 //....oooOO0OOooo........oooOO0OOooo........oo
400
401 void G4ICRU73QOModel::CorrectionsAlongStep(con
402 con
403 con
404 G4d
405 {}
406
407 //....oooOO0OOooo........oooOO0OOooo........oo
408
409 void G4ICRU73QOModel::SampleSecondaries(std::v
410 const
411 const
412 G4doub
413 G4doub
414 {
415 const G4double tmax = MaxSecondaryKinEnergy(
416 const G4double xmax = std::min(tmax, maxEner
417 const G4double xmin = std::max(minEnergy, lo
418 if(xmin >= xmax) { return; }
419
420 G4double kineticEnergy = dp->GetKineticEnerg
421 const G4double energy = kineticEnergy + mas
422 const G4double energy2 = energy*energy;
423 const G4double beta2 = kineticEnergy*(kine
424 G4double grej = 1.0;
425 G4double deltaKinEnergy, f;
426
427 G4ThreeVector direction = dp->GetMomentumDir
428
429 // sampling follows ...
430 do {
431 G4double x = G4UniformRand();
432 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - x)
433
434 f = 1.0 - beta2*deltaKinEnergy/tmax;
435
436 if(f > grej) {
437 G4cout << "G4ICRU73QOModel::SampleSeco
438 << "Majorant " << grej << " < "
439 << f << " for e= " << deltaKinE
440 << G4endl;
441 }
442
443 // Loop checking, 03-Aug-2015, Vladimir Iv
444 } while( grej*G4UniformRand() >= f );
445
446 G4ThreeVector deltaDirection;
447
448 if(UseAngularGeneratorFlag()) {
449 const G4Material* mat = couple->GetMateri
450 G4int Z = SelectRandomAtomNumber(mat);
451
452 deltaDirection =
453 GetAngularDistribution()->SampleDirectio
454
455 } else {
456
457 G4double deltaMomentum =
458 sqrt(deltaKinEnergy * (deltaKinEner
459 G4double totMomentum = energy*sqrt(beta2);
460 G4double cost = deltaKinEnergy * (energy +
461 (deltaMomen
462 if(cost > 1.0) { cost = 1.0; }
463 G4double sint = sqrt((1.0 - cost)*(1.0 + c
464
465 G4double phi = twopi * G4UniformRand() ;
466
467 deltaDirection.set(sint*cos(phi),sint*sin(
468 deltaDirection.rotateUz(direction);
469 }
470 // create G4DynamicParticle object for delta
471 auto delta = new G4DynamicParticle(theElectr
472
473 // Change kinematics of primary particle
474 kineticEnergy -= deltaKinEnergy;
475 G4ThreeVector finalP = dp->GetMomentum() - d
476 finalP = finalP.unit();
477
478 fParticleChange->SetProposedKineticEnergy(ki
479 fParticleChange->SetProposedMomentumDirectio
480
481 vdp->push_back(delta);
482 }
483
484 //....oooOO0OOooo........oooOO0OOooo........oo
485
486 G4double G4ICRU73QOModel::MaxSecondaryEnergy(c
487 G
488 {
489 if(pd != particle) { SetParticle(pd); }
490 G4double tau = kinEnergy/mass;
491 G4double tmax = 2.0*electron_mass_c2*tau*(ta
492 (1. + 2.0*(tau + 1.)*ratio +
493 return tmax;
494 }
495
496 //....oooOO0OOooo........oooOO0OOooo........oo
497
498 const G4int G4ICRU73QOModel::ZElementAvailable
499 {1,2,4,6,7,8,10,13,14,-18,
500 22,26,28,29,32,36,42,47,
501 50,54,73,74,78,79,82,92};
502
503 const G4int G4ICRU73QOModel::nbofShellsForElem
504 {1,1,2,3,3,3,3,4,5,4,
505 5,5,5,5,6,4,6,6,
506 7,6,6,8,7,7,9,9};
507
508 const G4int G4ICRU73QOModel::startElemIndex[NQ
509 {0,1,2,4,7,10,13,16,20,25,
510 29,34,39,44,49,55,59,65,
511 71,78,84,90,98,105,112,121};
512
513 //....oooOO0OOooo........oooOO0OOooo........oo
514
515 // SubShellOccupation = Z * ShellStrength
516 const G4double G4ICRU73QOModel::SubShellOccupa
517 {
518 1.000, // H 0
519 2.000, // He 1
520 1.930, 2.070, // Be 2-3
521 1.992, 1.841, 2.167, // C 4-6
522 1.741, 1.680, 3.579, // N 7-9
523 1.802, 1.849, 4.349, // O 10-12
524 1.788, 2.028, 6.184, // Ne 13-15
525 1.623, 2.147, 6.259, 2.971, // Al 16-19
526 1.631, 2.094, 6.588, 2.041, 1.646, // Si 2
527 1.535, 8.655, 1.706, 6.104, // Ar 25-28
528 1.581, 8.358, 8.183, 2.000, 1.878, // Ti 2
529 1.516, 8.325, 8.461, 6.579, 1.119, // Fe 3
530 1.422, 7.81, 8.385, 8.216, 2.167, // Ni 39
531 1.458, 8.049, 8.79, 9.695, 1.008, // Cu 44
532 1.442, 7.791, 7.837, 10.122, 2.463, 2.345,
533 1.645, 7.765, 19.192, 7.398, // Kr 55-58
534 1.313, 6.409, 19.229, 8.633, 5.036, 1.380,
535 1.295, 6.219, 18.751, 8.748, 10.184, 1.803
536 1.277, 6.099, 20.386, 8.011, 10.007, 2.272
537 1.563, 6.312, 21.868, 5.762, 11.245, 7.250
538 0.9198, 6.5408, 18.9727, 24.9149, 15.0161,
539 1.202, 5.582, 19.527, 18.741, 8.411, 14.38
540 1.159, 5.467, 18.802, 33.905, 8.300, 9.342
541 1.124, 5.331, 18.078, 34.604, 8.127, 10.41
542 2.000, 8.000, 18.000, 18.000, 14.000, 8.00
543 2.000, 8.000, 18.000, 32.000, 18.000, 8.00
544 };
545
546 //....oooOO0OOooo........oooOO0OOooo........oo
547
548 // ShellEnergy in eV
549 const G4double G4ICRU73QOModel::ShellEnergy[NQ
550 {
551 19.2, // H
552 41.8, // He
553 209.11, 21.68, // Be
554 486.2, 60.95, 23.43, // C
555 732.61, 100.646, 23.550, // N
556 965.1, 129.85, 31.60, // O
557 1525.9, 234.9, 56.18, // Ne
558 2701, 476.5, 150.42, 16.89, // Al
559 3206.1, 586.4, 186.8, 23.52, 14.91, // Si
560 5551.6, 472.43, 124.85, 22.332, // Ar
561 8554.6, 850.58, 93.47, 39.19, 19.46, // Ti
562 12254.7, 1279.29, 200.35, 49.19, 17.66, //
563 14346.9, 1532.28, 262.71, 74.37, 23.03, //
564 15438.5, 1667.96, 294.1, 70.69, 16.447, //
565 19022.1, 2150.79, 455.79, 179.87, 57.89, 2
566 24643, 2906.4, 366.85, 22.24, // Kr
567 34394, 4365.3, 589.36, 129.42, 35.59, 18.4
568 43664.3, 5824.91, 909.79, 175.47, 54.89, 1
569 49948, 6818.2, 1036.1, 172.65, 70.89, 33.8
570 58987, 8159, 1296.6, 356.75, 101.03, 16.52
571 88926, 18012, 3210, 575, 108.7, 30.8, // T
572 115025.9, 17827.44, 3214.36, 750.41, 305.2
573 128342, 20254, 3601.8, 608.1, 115.0, 42.75
574 131872, 20903, 3757.4, 682.1, 105.2, 44.89
575 154449, 25067, 5105.0, 987.44, 247.59, 188
576 167282, 27868, 6022.7, 1020.4, 244.81, 51.
577 };
578
579 //....oooOO0OOooo........oooOO0OOooo........oo
580
581 // Data for L0 from: Sigmund P., Haagerup U. P
582 const G4double G4ICRU73QOModel::L0[67][2] =
583 {
584 {0.00, 0.000001},
585 {0.10, 0.000001},
586 {0.12, 0.00001},
587 {0.14, 0.00005},
588 {0.16, 0.00014},
589 {0.18, 0.00030},
590 {0.20, 0.00057},
591 {0.25, 0.00189},
592 {0.30, 0.00429},
593 {0.35, 0.00784},
594 {0.40, 0.01248},
595 {0.45, 0.01811},
596 {0.50, 0.02462},
597 {0.60, 0.03980},
598 {0.70, 0.05731},
599 {0.80, 0.07662},
600 {0.90, 0.09733},
601 {1.00, 0.11916},
602 {1.20, 0.16532},
603 {1.40, 0.21376},
604 {1.60, 0.26362},
605 {1.80, 0.31428},
606 {2.00, 0.36532},
607 {2.50, 0.49272},
608 {3.00, 0.61765},
609 {3.50, 0.73863},
610 {4.00, 0.85496},
611 {4.50, 0.96634},
612 {5.00, 1.07272},
613 {6.00, 1.27086},
614 {7.00, 1.45075},
615 {8.00, 1.61412},
616 {9.00, 1.76277},
617 {10.00, 1.89836},
618 {12.00, 2.13625},
619 {14.00, 2.33787},
620 {16.00, 2.51093},
621 {18.00, 2.66134},
622 {20.00, 2.79358},
623 {25.00, 3.06539},
624 {30.00, 3.27902},
625 {35.00, 3.45430},
626 {40.00, 3.60281},
627 {45.00, 3.73167},
628 {50.00, 3.84555},
629 {60.00, 4.04011},
630 {70.00, 4.20264},
631 {80.00, 4.34229},
632 {90.00, 4.46474},
633 {100.00, 4.57378},
634 {120.00, 4.76155},
635 {140.00, 4.91953},
636 {160.00, 5.05590},
637 {180.00, 5.17588},
638 {200.00, 5.28299},
639 {250.00, 5.50925},
640 {300.00, 5.69364},
641 {350.00, 5.84926},
642 {400.00, 5.98388},
643 {450.00, 6.10252},
644 {500.00, 6.20856},
645 {600.00, 6.39189},
646 {700.00, 6.54677},
647 {800.00, 6.68084},
648 {900.00, 6.79905},
649 {1000.00, 6.90474}
650 };
651
652 //....oooOO0OOooo........oooOO0OOooo........oo
653
654 // Data for L1 from: Mikkelsen H.H., Sigmund P
655 const G4double G4ICRU73QOModel::L1[22][2] =
656 {
657 {0.00, -0.000001},
658 {0.10, -0.00001},
659 {0.20, -0.00049},
660 {0.30, -0.00084},
661 {0.40, 0.00085},
662 {0.50, 0.00519},
663 {0.60, 0.01198},
664 {0.70, 0.02074},
665 {0.80, 0.03133},
666 {0.90, 0.04369},
667 {1.00, 0.06035},
668 {2.00, 0.24023},
669 {3.00, 0.44284},
670 {4.00, 0.62012},
671 {5.00, 0.77031},
672 {6.00, 0.90390},
673 {7.00, 1.02705},
674 {8.00, 1.10867},
675 {9.00, 1.17546},
676 {10.00, 1.21599},
677 {15.00, 1.24349},
678 {20.00, 1.16752}
679 };
680
681 //....oooOO0OOooo........oooOO0OOooo........oo
682
683 // Data for L2 from: Mikkelsen H.H. Nucl. Inst
684 const G4double G4ICRU73QOModel::L2[14][2] =
685 {
686 {0.00, 0.000001},
687 {0.10, 0.00001},
688 {0.20, 0.00000},
689 {0.40, -0.00120},
690 {0.60, -0.00036},
691 {0.80, 0.00372},
692 {1.00, 0.01298},
693 {2.00, 0.08296},
694 {4.00, 0.21953},
695 {6.00, 0.23903},
696 {8.00, 0.20893},
697 {10.00, 0.10879},
698 {20.00, -0.88409},
699 {40.00, -1.13902}
700 };
701
702 //....oooOO0OOooo........oooOO0OOooo........oo
703
704 // Correction obtained by V.Ivanchenko using G
705 const G4double G4ICRU73QOModel::factorBethe[99
706 0.9637, 0.9872, 0.9469, 0.9875, 0.91, 0.989, 0
707 0.8357, 0.868, 0.9417, 0.9466, 0.8911, 0.905,
708 0.9098, 0.976, 0.8425, 0.8099, 0.7858, 0.947,
709 0.7223, 0.9784, 0.774, 0.7953, 0.829, 0.9405,
710 0.8207, 0.9033, 0.8063, 0.7837, 0.7818, 0.744,
711 0.8231, 0.8462, 0.853, 0.8736, 0.856, 0.8762,
712 0.7533, 0.7273, 0.7093, 0.7157, 0.6823, 0.6612
713 0.6497, 0.6746, 0.8568, 0.8541, 0.6958, 0.6962
714 0.7454, 0.78, 0.7783, 0.7996, 0.8216, 0.8632,
715 0.8321, 0.8272, 0.7999, 0.7934, 0.7787, 0.7851
716