ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeIonisationXSHandler.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PenelopeIonisationXSHandler.cc
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  :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 
69 G4PenelopeIonisationXSHandler::~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 
109 const G4PenelopeCrossSection*
110 G4PenelopeIonisationXSHandler::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 
159 void 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 
285 G4double 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 = std::log(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 
322 void 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  std::log(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,std::log(energy),delta);
419  }
420  theDeltaTable->insert(std::make_pair(mat,theVector));
421  return;
422 }
423 
424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
425 G4DataVector* 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 = std::log(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 = std::log(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)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee +
514  amol*(wu-wl)/(ee*ee);
515  H1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) +
516  (2.0-amol)*std::log((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*std::log((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)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee +
539  amol*(wu-wl)/(ee*ee);
540  S1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) +
541  (2.0-amol)*std::log((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*std::log((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....
557 G4DataVector* 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 = std::log(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 = std::log(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*std::log(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 += std::log(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*std::log(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 += std::log(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 }