| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // Author: Luciano Pandola
28 //
29 // History:
30 // --------
31 // 08 Mar 2012 L Pandola First complete implementation
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::G4PenelopeIonisationXSHandler(size_t nb)
47 :fXSTableElectron(nullptr),fXSTablePositron(nullptr),
48 fDeltaTable(nullptr),fEnergyGrid(nullptr)
49 {
50 fNBins = nb;
51 G4double LowEnergyLimit = 100.0*eV;
52 G4double HighEnergyLimit = 100.0*GeV;
53 fOscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
54 fXSTableElectron = new
55 std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
56 fXSTablePositron = new
57 std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
58
59 fDeltaTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
60 fEnergyGrid = new G4PhysicsLogVector(LowEnergyLimit,
61 HighEnergyLimit,
62 fNBins-1); //one hidden bin is added
63 fVerboseLevel = 0;
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67
68 G4PenelopeIonisationXSHandler::~G4PenelopeIonisationXSHandler()
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. Tables have been cleared"
103 << G4endl;
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107
108 const G4PenelopeCrossSection*
109 G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple(const G4ParticleDefinition* part,
110 const G4Material* mat,
111 const G4double cut) const
112 {
113 if (part != G4Electron::Electron() && part != G4Positron::Positron())
114 {
115 G4ExceptionDescription ed;
116 ed << "Invalid particle: " << part->GetParticleName() << G4endl;
117 G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
118 "em0001",FatalException,ed);
119 return nullptr;
120 }
121
122 if (part == G4Electron::Electron())
123 {
124 if (!fXSTableElectron)
125 {
126 G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
127 "em0028",FatalException,
128 "The Cross Section Table for e- was not initialized correctly!");
129 return nullptr;
130 }
131 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
132 if (fXSTableElectron->count(theKey)) //table already built
133 return fXSTableElectron->find(theKey)->second;
134 else
135 return nullptr;
136 }
137
138 if (part == G4Positron::Positron())
139 {
140 if (!fXSTablePositron)
141 {
142 G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
143 "em0028",FatalException,
144 "The Cross Section Table for e+ was not initialized correctly!");
145 return nullptr;
146 }
147 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
148 if (fXSTablePositron->count(theKey)) //table already built
149 return fXSTablePositron->find(theKey)->second;
150 else
151 return nullptr;
152 }
153 return nullptr;
154 }
155
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157
158 void G4PenelopeIonisationXSHandler::BuildXSTable(const G4Material* mat,G4double cut,
159 const G4ParticleDefinition* part,
160 G4bool isMaster)
161 {
162 //Just to check
163 if (!isMaster)
164 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
165 "em0100",FatalException,"Worker thread in this method");
166
167 //
168 //This method fills the G4PenelopeCrossSection containers for electrons or positrons
169 //and for the given material/cut couple. The calculation is done as sum over the
170 //individual shells.
171 //Equivalent of subroutines EINaT and PINaT of Penelope
172 //
173 if (fVerboseLevel > 2)
174 {
175 G4cout << "G4PenelopeIonisationXSHandler: going to build cross section table " << G4endl;
176 G4cout << "for " << part->GetParticleName() << " in " << mat->GetName() << G4endl;
177 G4cout << "Cut= " << cut/keV << " keV" << G4endl;
178 }
179
180 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
181 //Check if the table already exists
182 if (part == G4Electron::Electron())
183 {
184 if (fXSTableElectron->count(theKey)) //table already built
185 return;
186 }
187 if (part == G4Positron::Positron())
188 {
189 if (fXSTablePositron->count(theKey)) //table already built
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 by GetCrossSectionTableForCouple)
199 G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableIonisation(mat);
200 size_t numberOfOscillators = theTable->size();
201
202 if (fEnergyGrid->GetVectorLength() != fNBins)
203 {
204 G4ExceptionDescription ed;
205 ed << "Energy Grid looks not initialized" << G4endl;
206 ed << fNBins << " " << fEnergyGrid->GetVectorLength() << G4endl;
207 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
208 "em2030",FatalException,ed);
209 }
210
211 G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(fNBins,numberOfOscillators);
212
213 //loop on the energy grid
214 for (size_t bin=0;bin<fNBins;bin++)
215 {
216 G4double energy = fEnergyGrid->GetLowEdgeEnergy(bin);
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<numberOfOscillators;iosc++)
222 {
223 G4DataVector* tempStorage = nullptr;
224
225 G4PenelopeOscillator* theOsc = (*theTable)[iosc];
226 G4double delta = GetDensityCorrection(mat,energy);
227 if (part == G4Electron::Electron())
228 tempStorage = ComputeShellCrossSectionsElectron(theOsc,energy,cut,delta);
229 else if (part == G4Positron::Positron())
230 tempStorage = ComputeShellCrossSectionsPositron(theOsc,energy,cut,delta);
231 //check results are all right
232 if (!tempStorage)
233 {
234 G4ExceptionDescription ed;
235 ed << "Problem in calculating the shell XS for shell # "
236 << iosc << G4endl;
237 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
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 shell XS " << G4endl;
246 ed << "Result has dimension " << tempStorage->size() << " instead of 6" << G4endl;
247 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
248 "em2031",FatalException,ed);
249 }
250 G4double stre = theOsc->GetOscillatorStrength();
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,iosc,energy,stre*(*tempStorage)[0]);
259 if (tempStorage)
260 {
261 delete tempStorage;
262 tempStorage = nullptr;
263 }
264 }
265 XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2);
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(theKey,XSEntry));
273 else if (part == G4Positron::Positron())
274 fXSTablePositron->insert(std::make_pair(theKey,XSEntry));
275 else
276 delete XSEntry;
277
278 return;
279 }
280
281
282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
283
284 G4double G4PenelopeIonisationXSHandler::GetDensityCorrection(const G4Material* mat,
285 const G4double energy) const
286 {
287 G4double result = 0;
288 if (!fDeltaTable)
289 {
290 G4Exception("G4PenelopeIonisationXSHandler::GetDensityCorrection()",
291 "em2032",FatalException,
292 "Delta Table not initialized. Was Initialise() run?");
293 return 0;
294 }
295 if (energy <= 0*eV)
296 {
297 G4cout << "G4PenelopeIonisationXSHandler::GetDensityCorrection()" << G4endl;
298 G4cout << "Invalid energy " << energy/eV << " eV " << G4endl;
299 return 0;
300 }
301 G4double logene = G4Log(energy);
302
303 if (fDeltaTable->count(mat))
304 {
305 const G4PhysicsFreeVector* vec = fDeltaTable->find(mat)->second;
306 result = vec->Value(logene); //the table has delta vs. ln(E)
307 }
308 else
309 {
310 G4ExceptionDescription ed;
311 ed << "Unable to build table for " << mat->GetName() << G4endl;
312 G4Exception("G4PenelopeIonisationXSHandler::GetDensityCorrection()",
313 "em2033",FatalException,ed);
314 }
315
316 return result;
317 }
318
319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
320
321 void G4PenelopeIonisationXSHandler::BuildDeltaTable(const G4Material* mat)
322 {
323 G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableIonisation(mat);
324 G4double plasmaSq = fOscManager->GetPlasmaEnergySquared(mat);
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 not initialized" << G4endl;
332 ed << fNBins << " " << fEnergyGrid->GetVectorLength() << G4endl;
333 G4Exception("G4PenelopeIonisationXSHandler::BuildDeltaTable()",
334 "em2030",FatalException,ed);
335 }
336
337 G4PhysicsFreeVector* theVector = new G4PhysicsFreeVector(fNBins);
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->GetLowEdgeEnergy(bin);
344
345 //Here calculate delta
346 G4double gam = 1.0+(energy/electron_mass_c2);
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)[i];
357 G4double wri = theOsc->GetResonanceEnergy();
358 fdel += theOsc->GetOscillatorStrength()/(wri*wri+wl2);
359 }
360 if (fdel >= TST) //if fdel < TST, delta = 0
361 {
362 //get last oscillator
363 G4PenelopeOscillator* theOsc = (*theTable)[numberOfOscillators-1];
364 wl2 = theOsc->GetResonanceEnergy()*theOsc->GetResonanceEnergy();
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;i++)
374 {
375 G4PenelopeOscillator* theOscLocal1 = (*theTable)[i];
376 G4double wri = theOscLocal1->GetResonanceEnergy();
377 fdel += theOscLocal1->GetOscillatorStrength()/(wri*wri+wl2);
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;i++)
392 {
393 G4PenelopeOscillator* theOscLocal2 = (*theTable)[i];
394 G4double wri = theOscLocal2->GetResonanceEnergy();
395 fdel += theOscLocal2->GetOscillatorStrength()/(wri*wri+wl2);
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 = (*theTable)[i];
410 G4double wri = theOscLocal3->GetResonanceEnergy();
411 delta += theOscLocal3->GetOscillatorStrength()*
412 G4Log(1.0+(wl2/(wri*wri)));
413 }
414 delta = (delta/totalZ)-wl2/(gamSq*plasmaSq);
415 }
416 energy = std::max(1e-9*eV,energy); //prevents log(0)
417 theVector->PutValue(bin,G4Log(energy),delta);
418 }
419 fDeltaTable->insert(std::make_pair(mat,theVector));
420 return;
421 }
422
423 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
424 G4DataVector* G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsElectron(G4PenelopeOscillator* theOsc,
425 G4double energy,
426 G4double cut,
427 G4double delta)
428 {
429 //
430 //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for
431 //the given oscillator/cut and at the given energy.
432 //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2)
433 //Equivalent of subroutines EINaT1 of Penelope
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->GetIonisationEnergy();
441
442 //return a set of zero's if the energy it too low to excite the current oscillator
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 calculation
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*classic_electr_radius; //pi*re^2
454 G4double constant = pielr2*2.0*electron_mass_c2/beta;
455 G4double XHDT0 = G4Log(gammaSq)-beta;
456
457 G4double cpSq = energy*(energy+2.0*electron_mass_c2);
458 G4double cp = std::sqrt(cpSq);
459 G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));
460
461 //
462 // Distant interactions
463 //
464 G4double resEne = theOsc->GetResonanceEnergy();
465 G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy();
466 if (energy > resEne)
467 {
468 G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
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_mass_c2*electron_mass_c2)-electron_mass_c2;
475 else
476 {
477 QM = resEne*resEne/(beta*2.0*electron_mass_c2);
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_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
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.0/wu) + (1.0/wl) +
512 (1.0-amol)*G4Log(((ee-wu)*wl)/((ee-wl)*wu))/ee +
513 amol*(wu-wl)/(ee*ee);
514 H1 += G4Log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) +
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)/(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) + (1.0/wl) +
537 (1.0-amol)*G4Log(((ee-wu)*wl)/((ee-wl)*wu))/ee +
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-wu)) -
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........oooOO0OOooo........oooOO0OOooo....
556 G4DataVector* G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsPositron(G4PenelopeOscillator* theOsc,
557 G4double energy,
558 G4double cut,
559 G4double delta)
560 {
561 //
562 //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for
563 //the given oscillator/cut and at the given energy.
564 //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2)
565 //Equivalent of subroutines PINaT1 of Penelope
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->GetIonisationEnergy();
573
574 //return a set of zero's if the energy it too low to excite the current oscillator
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 calculation
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*classic_electr_radius; //pi*re^2
586 G4double constant = pielr2*2.0*electron_mass_c2/beta;
587 G4double XHDT0 = G4Log(gammaSq)-beta;
588
589 G4double cpSq = energy*(energy+2.0*electron_mass_c2);
590 G4double cp = std::sqrt(cpSq);
591 G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2));
592 G4double g12 = (gamma+1.0)*(gamma+1.0);
593 //Bhabha coefficients
594 G4double bha1 = amol*(2.0*g12-1.0)/(gammaSq-1.0);
595 G4double bha2 = amol*(3.0+1.0/g12);
596 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/g12;
597 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/g12;
598
599 //
600 // Distant interactions
601 //
602 G4double resEne = theOsc->GetResonanceEnergy();
603 G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy();
604 if (energy > resEne)
605 {
606 G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2);
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_mass_c2*electron_mass_c2)-electron_mass_c2;
613 else
614 {
615 QM = resEne*resEne/(beta*2.0*electron_mass_c2);
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_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)));
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/wl)/energy
652 + bha2*(wu-wl)/energySq
653 - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
654 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);
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*energy)
658 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq);
659 H2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy)
660 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq)
661 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy)
662 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);
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)/energy
682 + bha2*(wu-wl)/energySq
683 - bha3*(wuSq-wlSq)/(2.0*energySq*energy)
684 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq);
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*energy)
689 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*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*energy)
694 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq);
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