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