| 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 // Author: Luciano Pandola
28 //
29 // History:
30 // --------
31 // 08 Mar 2012 L Pandola First complete i
32 //
33
34 #include "G4PenelopeIonisationXSHandler.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4ParticleDefinition.hh"
38 #include "G4Electron.hh"
39 #include "G4Positron.hh"
40 #include "G4PenelopeOscillatorManager.hh"
41 #include "G4PenelopeOscillator.hh"
42 #include "G4PenelopeCrossSection.hh"
43 #include "G4PhysicsFreeVector.hh"
44 #include "G4PhysicsLogVector.hh"
45
46 G4PenelopeIonisationXSHandler::G4PenelopeIonis
47 :fXSTableElectron(nullptr),fXSTablePositron(
48 fDeltaTable(nullptr),fEnergyGrid(nullptr)
49 {
50 fNBins = nb;
51 G4double LowEnergyLimit = 100.0*eV;
52 G4double HighEnergyLimit = 100.0*GeV;
53 fOscManager = G4PenelopeOscillatorManager::G
54 fXSTableElectron = new
55 std::map< std::pair<const G4Material*,G4do
56 fXSTablePositron = new
57 std::map< std::pair<const G4Material*,G4do
58
59 fDeltaTable = new std::map<const G4Material*
60 fEnergyGrid = new G4PhysicsLogVector(LowEner
61 HighEnergyLimit,
62 fNBins-1); //one hidden bin is a
63 fVerboseLevel = 0;
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oo
67
68 G4PenelopeIonisationXSHandler::~G4PenelopeIoni
69 {
70 if (fXSTableElectron)
71 {
72 for (auto& item : (*fXSTableElectron))
73 {
74 //G4PenelopeCrossSection* tab = i->second;
75 delete item.second;
76 }
77 delete fXSTableElectron;
78 fXSTableElectron = nullptr;
79 }
80
81 if (fXSTablePositron)
82 {
83 for (auto& item : (*fXSTablePositron))
84 {
85 //G4PenelopeCrossSection* tab = i->second;
86 delete item.second;
87 }
88 delete fXSTablePositron;
89 fXSTablePositron = nullptr;
90 }
91 if (fDeltaTable)
92 {
93 for (auto& item : (*fDeltaTable))
94 delete item.second;
95 delete fDeltaTable;
96 fDeltaTable = nullptr;
97 }
98 if (fEnergyGrid)
99 delete fEnergyGrid;
100
101 if (fVerboseLevel > 2)
102 G4cout << "G4PenelopeIonisationXSHandler.
103 << G4endl;
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oo
107
108 const G4PenelopeCrossSection*
109 G4PenelopeIonisationXSHandler::GetCrossSection
110 const G4Material* mat,
111 const G4double cut) const
112 {
113 if (part != G4Electron::Electron() && part !
114 {
115 G4ExceptionDescription ed;
116 ed << "Invalid particle: " << part->GetP
117 G4Exception("G4PenelopeIonisationXSHandl
118 "em0001",FatalException,ed);
119 return nullptr;
120 }
121
122 if (part == G4Electron::Electron())
123 {
124 if (!fXSTableElectron)
125 {
126 G4Exception("G4PenelopeIonisationXSHandler
127 "em0028",FatalException,
128 "The Cross Section Table for e- was
129 return nullptr;
130 }
131 std::pair<const G4Material*,G4double> th
132 if (fXSTableElectron->count(theKey)) //t
133 return fXSTableElectron->find(theKey)->secon
134 else
135 return nullptr;
136 }
137
138 if (part == G4Positron::Positron())
139 {
140 if (!fXSTablePositron)
141 {
142 G4Exception("G4PenelopeIonisationXSHandler
143 "em0028",FatalException,
144 "The Cross Section Table for e+ was
145 return nullptr;
146 }
147 std::pair<const G4Material*,G4double> th
148 if (fXSTablePositron->count(theKey)) //t
149 return fXSTablePositron->find(theKey)->secon
150 else
151 return nullptr;
152 }
153 return nullptr;
154 }
155
156 //....oooOO0OOooo........oooOO0OOooo........oo
157
158 void G4PenelopeIonisationXSHandler::BuildXSTab
159 const G4ParticleDefinition* part,
160 G4bool isMaster)
161 {
162 //Just to check
163 if (!isMaster)
164 G4Exception("G4PenelopeIonisationXSHandler
165 "em0100",FatalException,"Worke
166
167 //
168 //This method fills the G4PenelopeCrossSecti
169 //and for the given material/cut couple. The
170 //individual shells.
171 //Equivalent of subroutines EINaT and PINaT
172 //
173 if (fVerboseLevel > 2)
174 {
175 G4cout << "G4PenelopeIonisationXSHandler
176 G4cout << "for " << part->GetParticleNam
177 G4cout << "Cut= " << cut/keV << " keV" <
178 }
179
180 std::pair<const G4Material*,G4double> theKey
181 //Check if the table already exists
182 if (part == G4Electron::Electron())
183 {
184 if (fXSTableElectron->count(theKey)) //t
185 return;
186 }
187 if (part == G4Positron::Positron())
188 {
189 if (fXSTablePositron->count(theKey)) //t
190 return;
191 }
192
193 //check if the material has been built
194 if (!(fDeltaTable->count(mat)))
195 BuildDeltaTable(mat);
196
197
198 //Tables have been already created (checked
199 G4PenelopeOscillatorTable* theTable = fOscMa
200 size_t numberOfOscillators = theTable->size(
201
202 if (fEnergyGrid->GetVectorLength() != fNBins
203 {
204 G4ExceptionDescription ed;
205 ed << "Energy Grid looks not initialized
206 ed << fNBins << " " << fEnergyGrid->GetV
207 G4Exception("G4PenelopeIonisationXSHandl
208 "em2030",FatalException,ed);
209 }
210
211 G4PenelopeCrossSection* XSEntry = new G4Pene
212
213 //loop on the energy grid
214 for (size_t bin=0;bin<fNBins;bin++)
215 {
216 G4double energy = fEnergyGrid->GetLowEd
217 G4double XH0=0, XH1=0, XH2=0;
218 G4double XS0=0, XS1=0, XS2=0;
219
220 //oscillator loop
221 for (size_t iosc=0;iosc<numberOfOscilla
222 {
223 G4DataVector* tempStorage = nullptr;
224
225 G4PenelopeOscillator* theOsc = (*theTable
226 G4double delta = GetDensityCorrection(mat
227 if (part == G4Electron::Electron())
228 tempStorage = ComputeShellCrossSections
229 else if (part == G4Positron::Positron())
230 tempStorage = ComputeShellCrossSections
231 //check results are all right
232 if (!tempStorage)
233 {
234 G4ExceptionDescription ed;
235 ed << "Problem in calculating the she
236 << iosc << G4endl;
237 G4Exception("G4PenelopeIonisationXSHa
238 "em2031",FatalException,ed);
239 delete XSEntry;
240 return;
241 }
242 if (tempStorage->size() != 6)
243 {
244 G4ExceptionDescription ed;
245 ed << "Problem in calculating the she
246 ed << "Result has dimension " << temp
247 G4Exception("G4PenelopeIonisationXSHa
248 "em2031",FatalException,ed);
249 }
250 G4double stre = theOsc->GetOscillatorStre
251
252 XH0 += stre*(*tempStorage)[0];
253 XH1 += stre*(*tempStorage)[1];
254 XH2 += stre*(*tempStorage)[2];
255 XS0 += stre*(*tempStorage)[3];
256 XS1 += stre*(*tempStorage)[4];
257 XS2 += stre*(*tempStorage)[5];
258 XSEntry->AddShellCrossSectionPoint(bin,io
259 if (tempStorage)
260 {
261 delete tempStorage;
262 tempStorage = nullptr;
263 }
264 }
265 XSEntry->AddCrossSectionPoint(bin,energ
266 }
267 //Do (only once) the final normalization
268 XSEntry->NormalizeShellCrossSections();
269
270 //Insert in the appropriate table
271 if (part == G4Electron::Electron())
272 fXSTableElectron->insert(std::make_pair(th
273 else if (part == G4Positron::Positron())
274 fXSTablePositron->insert(std::make_pair(th
275 else
276 delete XSEntry;
277
278 return;
279 }
280
281
282 //....oooOO0OOooo........oooOO0OOooo........oo
283
284 G4double G4PenelopeIonisationXSHandler::GetDen
285 const G4double energy) cons
286 {
287 G4double result = 0;
288 if (!fDeltaTable)
289 {
290 G4Exception("G4PenelopeIonisationXSHandl
291 "em2032",FatalException,
292 "Delta Table not initialized. Was Initia
293 return 0;
294 }
295 if (energy <= 0*eV)
296 {
297 G4cout << "G4PenelopeIonisationXSHandler
298 G4cout << "Invalid energy " << energy/eV
299 return 0;
300 }
301 G4double logene = G4Log(energy);
302
303 if (fDeltaTable->count(mat))
304 {
305 const G4PhysicsFreeVector* vec = fDeltaT
306 result = vec->Value(logene); //the table
307 }
308 else
309 {
310 G4ExceptionDescription ed;
311 ed << "Unable to build table for " << ma
312 G4Exception("G4PenelopeIonisationXSHandl
313 "em2033",FatalException,ed);
314 }
315
316 return result;
317 }
318
319 //....oooOO0OOooo........oooOO0OOooo........oo
320
321 void G4PenelopeIonisationXSHandler::BuildDelta
322 {
323 G4PenelopeOscillatorTable* theTable = fOscMa
324 G4double plasmaSq = fOscManager->GetPlasmaEn
325 G4double totalZ = fOscManager->GetTotalZ(mat
326 size_t numberOfOscillators = theTable->size(
327
328 if (fEnergyGrid->GetVectorLength() != fNBins
329 {
330 G4ExceptionDescription ed;
331 ed << "Energy Grid for Delta table looks
332 ed << fNBins << " " << fEnergyGrid->GetV
333 G4Exception("G4PenelopeIonisationXSHandl
334 "em2030",FatalException,ed);
335 }
336
337 G4PhysicsFreeVector* theVector = new G4Physi
338
339 //loop on the energy grid
340 for (size_t bin=0;bin<fNBins;bin++)
341 {
342 G4double delta = 0.;
343 G4double energy = fEnergyGrid->GetLowEdg
344
345 //Here calculate delta
346 G4double gam = 1.0+(energy/electron_mass
347 G4double gamSq = gam*gam;
348
349 G4double TST = totalZ/(gamSq*plasmaSq);
350 G4double wl2 = 0;
351 G4double fdel = 0;
352
353 //loop on oscillators
354 for (size_t i=0;i<numberOfOscillators;i+
355 {
356 G4PenelopeOscillator* theOsc = (*theTable)
357 G4double wri = theOsc->GetResonanceEnergy(
358 fdel += theOsc->GetOscillatorStrength()/(w
359 }
360 if (fdel >= TST) //if fdel < TST, delta
361 {
362 //get last oscillator
363 G4PenelopeOscillator* theOsc = (*theTable)
364 wl2 = theOsc->GetResonanceEnergy()*theOsc-
365
366 //First iteration
367 G4bool loopAgain = false;
368 do
369 {
370 loopAgain = false;
371 wl2 += wl2;
372 fdel = 0.;
373 for (size_t i=0;i<numberOfOscillators;
374 {
375 G4PenelopeOscillator* theOscLocal1 = (*t
376 G4double wri = theOscLocal1->GetResonanc
377 fdel += theOscLocal1->GetOscillatorStren
378 }
379 if (fdel > TST)
380 loopAgain = true;
381 }while(loopAgain);
382
383 G4double wl2l = 0;
384 G4double wl2u = wl2;
385 //second iteration
386 do
387 {
388 loopAgain = false;
389 wl2 = 0.5*(wl2l+wl2u);
390 fdel = 0;
391 for (size_t i=0;i<numberOfOscillators;
392 {
393 G4PenelopeOscillator* theOscLocal2 = (*t
394 G4double wri = theOscLocal2->GetResonanc
395 fdel += theOscLocal2->GetOscillatorStren
396 }
397 if (fdel > TST)
398 wl2l = wl2;
399 else
400 wl2u = wl2;
401 if ((wl2u-wl2l)>1e-12*wl2)
402 loopAgain = true;
403 }while(loopAgain);
404
405 //Eventually get density correction
406 delta = 0.;
407 for (size_t i=0;i<numberOfOscillators;i++)
408 {
409 G4PenelopeOscillator* theOscLocal3 = (
410 G4double wri = theOscLocal3->GetResona
411 delta += theOscLocal3->GetOscillatorSt
412 G4Log(1.0+(wl2/(wri*wri)));
413 }
414 delta = (delta/totalZ)-wl2/(gamSq*plasmaSq
415 }
416 energy = std::max(1e-9*eV,energy); //pre
417 theVector->PutValue(bin,G4Log(energy),de
418 }
419 fDeltaTable->insert(std::make_pair(mat,theVe
420 return;
421 }
422
423 //....oooOO0OOooo........oooOO0OOooo........oo
424 G4DataVector* G4PenelopeIonisationXSHandler::C
425 G4double energy,
426 G4double cut,
427 G4double delta)
428 {
429 //
430 //This method calculates the hard and soft c
431 //the given oscillator/cut and at the given
432 //It returns a G4DataVector* with 6 entries
433 //Equivalent of subroutines EINaT1 of Penelo
434 //
435 // Results are _per target electron_
436 //
437 G4DataVector* result = new G4DataVector();
438 for (size_t i=0;i<6;i++)
439 result->push_back(0.);
440 G4double ionEnergy = theOsc->GetIonisationEn
441
442 //return a set of zero's if the energy it to
443 if (energy < ionEnergy)
444 return result;
445
446 G4double H0=0.,H1=0.,H2=0.;
447 G4double S0=0.,S1=0.,S2=0.;
448
449 //Define useful constants to be used in the
450 G4double gamma = 1.0+energy/electron_mass_c2
451 G4double gammaSq = gamma*gamma;
452 G4double beta = (gammaSq-1.0)/gammaSq;
453 G4double pielr2 = pi*classic_electr_radius*c
454 G4double constant = pielr2*2.0*electron_mass
455 G4double XHDT0 = G4Log(gammaSq)-beta;
456
457 G4double cpSq = energy*(energy+2.0*electron_
458 G4double cp = std::sqrt(cpSq);
459 G4double amol = (energy/(energy+electron_mas
460
461 //
462 // Distant interactions
463 //
464 G4double resEne = theOsc->GetResonanceEnergy
465 G4double cutoffEne = theOsc->GetCutoffRecoil
466 if (energy > resEne)
467 {
468 G4double cp1Sq = (energy-resEne)*(energy
469 G4double cp1 = std::sqrt(cp1Sq);
470
471 //Distant longitudinal interactions
472 G4double QM = 0;
473 if (resEne > 1e-6*energy)
474 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_ma
475 else
476 {
477 QM = resEne*resEne/(beta*2.0*electron_mass
478 QM = QM*(1.0-0.5*QM/electron_mass_c2);
479 }
480 G4double SDL1 = 0;
481 if (QM < cutoffEne)
482 SDL1 = G4Log(cutoffEne*(QM+2.0*electron_mass
483
484 //Distant transverse interactions
485 if (SDL1)
486 {
487 G4double SDT1 = std::max(XHDT0-delta,0.0);
488 G4double SD1 = SDL1+SDT1;
489 if (cut > resEne)
490 {
491 S1 = SD1; //XS1
492 S0 = SD1/resEne; //XS0
493 S2 = SD1*resEne; //XS2
494 }
495 else
496 {
497 H1 = SD1; //XH1
498 H0 = SD1/resEne; //XH0
499 H2 = SD1*resEne; //XH2
500 }
501 }
502 }
503 //
504 // Close collisions (Moller's cross section)
505 //
506 G4double wl = std::max(cut,cutoffEne);
507 G4double ee = energy + ionEnergy;
508 G4double wu = 0.5*ee;
509 if (wl < wu-(1e-5*eV))
510 {
511 H0 += (1.0/(ee-wu)) - (1.0/(ee-wl)) - (1
512 (1.0-amol)*G4Log(((ee-wu)*wl)/((ee-wl)*wu))/
513 amol*(wu-wl)/(ee*ee);
514 H1 += G4Log(wu/wl)+(ee/(ee-wu))-(ee/(ee-
515 (2.0-amol)*G4Log((ee-wu)/(ee-wl)) +
516 amol*(wu*wu-wl*wl)/(2.0*ee*ee);
517 H2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)
518 (wl*(2.0*ee-wl)/(ee-wl)) +
519 (3.0-amol)*ee*G4Log((ee-wu)/(ee-wl)) +
520 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);
521 wu = wl;
522 }
523 wl = cutoffEne;
524
525 if (wl > wu-(1e-5*eV))
526 {
527 (*result)[0] = constant*H0;
528 (*result)[1] = constant*H1;
529 (*result)[2] = constant*H2;
530 (*result)[3] = constant*S0;
531 (*result)[4] = constant*S1;
532 (*result)[5] = constant*S2;
533 return result;
534 }
535
536 S0 += (1.0/(ee-wu))-(1.0/(ee-wl)) - (1.0/wu)
537 (1.0-amol)*G4Log(((ee-wu)*wl)/((ee-wl)*wu)
538 amol*(wu-wl)/(ee*ee);
539 S1 += G4Log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl))
540 (2.0-amol)*G4Log((ee-wu)/(ee-wl)) +
541 amol*(wu*wu-wl*wl)/(2.0*ee*ee);
542 S2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee
543 (wl*(2.0*ee-wl)/(ee-wl)) +
544 (3.0-amol)*ee*G4Log((ee-wu)/(ee-wl)) +
545 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee);
546
547 (*result)[0] = constant*H0;
548 (*result)[1] = constant*H1;
549 (*result)[2] = constant*H2;
550 (*result)[3] = constant*S0;
551 (*result)[4] = constant*S1;
552 (*result)[5] = constant*S2;
553 return result;
554 }
555 //....oooOO0OOooo........oooOO0OOooo........oo
556 G4DataVector* G4PenelopeIonisationXSHandler::C
557 G4double energy,
558 G4double cut,
559 G4double delta)
560 {
561 //
562 //This method calculates the hard and soft c
563 //the given oscillator/cut and at the given
564 //It returns a G4DataVector* with 6 entries
565 //Equivalent of subroutines PINaT1 of Penelo
566 //
567 // Results are _per target electron_
568 //
569 G4DataVector* result = new G4DataVector();
570 for (size_t i=0;i<6;i++)
571 result->push_back(0.);
572 G4double ionEnergy = theOsc->GetIonisationEn
573
574 //return a set of zero's if the energy it to
575 if (energy < ionEnergy)
576 return result;
577
578 G4double H0=0.,H1=0.,H2=0.;
579 G4double S0=0.,S1=0.,S2=0.;
580
581 //Define useful constants to be used in the
582 G4double gamma = 1.0+energy/electron_mass_c2
583 G4double gammaSq = gamma*gamma;
584 G4double beta = (gammaSq-1.0)/gammaSq;
585 G4double pielr2 = pi*classic_electr_radius*c
586 G4double constant = pielr2*2.0*electron_mass
587 G4double XHDT0 = G4Log(gammaSq)-beta;
588
589 G4double cpSq = energy*(energy+2.0*electron_
590 G4double cp = std::sqrt(cpSq);
591 G4double amol = (energy/(energy+electron_mas
592 G4double g12 = (gamma+1.0)*(gamma+1.0);
593 //Bhabha coefficients
594 G4double bha1 = amol*(2.0*g12-1.0)/(gammaSq-
595 G4double bha2 = amol*(3.0+1.0/g12);
596 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/g
597 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)
598
599 //
600 // Distant interactions
601 //
602 G4double resEne = theOsc->GetResonanceEnergy
603 G4double cutoffEne = theOsc->GetCutoffRecoil
604 if (energy > resEne)
605 {
606 G4double cp1Sq = (energy-resEne)*(energy
607 G4double cp1 = std::sqrt(cp1Sq);
608
609 //Distant longitudinal interactions
610 G4double QM = 0;
611 if (resEne > 1e-6*energy)
612 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_ma
613 else
614 {
615 QM = resEne*resEne/(beta*2.0*electron_mass
616 QM = QM*(1.0-0.5*QM/electron_mass_c2);
617 }
618 G4double SDL1 = 0;
619 if (QM < cutoffEne)
620 SDL1 = G4Log(cutoffEne*(QM+2.0*electron_mass
621
622 //Distant transverse interactions
623 if (SDL1)
624 {
625 G4double SDT1 = std::max(XHDT0-delta,0.0);
626 G4double SD1 = SDL1+SDT1;
627 if (cut > resEne)
628 {
629 S1 = SD1; //XS1
630 S0 = SD1/resEne; //XS0
631 S2 = SD1*resEne; //XS2
632 }
633 else
634 {
635 H1 = SD1; //XH1
636 H0 = SD1/resEne; //XH0
637 H2 = SD1*resEne; //XH2
638 }
639 }
640 }
641 //
642 // Close collisions (Bhabha's cross section)
643 //
644 G4double wl = std::max(cut,cutoffEne);
645 G4double wu = energy;
646 G4double energySq = energy*energy;
647 if (wl < wu-(1e-5*eV))
648 {
649 G4double wlSq = wl*wl;
650 G4double wuSq = wu*wu;
651 H0 += (1.0/wl) - (1.0/wu)- bha1*G4Log(wu
652 + bha2*(wu-wl)/energySq
653 - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
654 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energ
655 H1 += G4Log(wu/wl) - bha1*(wu-wl)/energy
656 + bha2*(wuSq-wlSq)/(2.0*energySq)
657 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energ
658 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*e
659 H2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*en
660 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
661 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*e
662 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*ener
663 wu = wl;
664 }
665 wl = cutoffEne;
666
667 if (wl > wu-(1e-5*eV))
668 {
669 (*result)[0] = constant*H0;
670 (*result)[1] = constant*H1;
671 (*result)[2] = constant*H2;
672 (*result)[3] = constant*S0;
673 (*result)[4] = constant*S1;
674 (*result)[5] = constant*S2;
675 return result;
676 }
677
678 G4double wlSq = wl*wl;
679 G4double wuSq = wu*wu;
680
681 S0 += (1.0/wl) - (1.0/wu) - bha1*G4Log(wu/wl
682 + bha2*(wu-wl)/energySq
683 - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
684 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*ene
685
686 S1 += G4Log(wu/wl) - bha1*(wu-wl)/energy
687 + bha2*(wuSq-wlSq)/(2.0*energySq)
688 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*ene
689 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq
690
691 S2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy
692 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
693 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq
694 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*en
695
696 (*result)[0] = constant*H0;
697 (*result)[1] = constant*H1;
698 (*result)[2] = constant*H2;
699 (*result)[3] = constant*S0;
700 (*result)[4] = constant*S1;
701 (*result)[5] = constant*S2;
702
703 return result;
704 }
705