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