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