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