Geant4  9.6.p02
 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$
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  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 
114 G4PenelopeCrossSection*
115 G4PenelopeIonisationXSHandler::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 
186 void 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 
286 G4double 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 
327 void 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....
430 G4DataVector* 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....
562 G4DataVector* 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 }
Generated on Sat May 25 2013 14:33:34 for Geant4 by   doxygen 1.8.4