Geant4  10.02
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 83584 2014-09-02 08:45:37Z 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  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 
101  delete theDeltaTable;
102  theDeltaTable = 0;
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 const G4PenelopeCrossSection*
115 G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple(const G4ParticleDefinition* part,
116  const G4Material* mat,
117  const G4double cut) const
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  return NULL;
142  }
143 
144  if (part == G4Positron::Positron())
145  {
146  if (!XSTablePositron)
147  {
148  G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()",
149  "em0028",FatalException,
150  "The Cross Section Table for e+ was not initialized correctly!");
151  return NULL;
152  }
153  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
154  if (XSTablePositron->count(theKey)) //table already built
155  return XSTablePositron->find(theKey)->second;
156  else
157  return NULL;
158  }
159  return NULL;
160 }
161 
162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
163 
164 void G4PenelopeIonisationXSHandler::BuildXSTable(const G4Material* mat,G4double cut,
165  const G4ParticleDefinition* part,
166  G4bool isMaster)
167 {
168  //Just to check
169  if (!isMaster)
170  G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
171  "em0100",FatalException,"Worker thread in this method");
172 
173  //
174  //This method fills the G4PenelopeCrossSection containers for electrons or positrons
175  //and for the given material/cut couple. The calculation is done as sum over the
176  //individual shells.
177  //Equivalent of subroutines EINaT and PINaT of Penelope
178  //
179  if (verboseLevel > 2)
180  {
181  G4cout << "G4PenelopeIonisationXSHandler: going to build cross section table " << G4endl;
182  G4cout << "for " << part->GetParticleName() << " in " << mat->GetName() << G4endl;
183  G4cout << "Cut= " << cut/keV << " keV" << G4endl;
184  }
185 
186  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
187  //Check if the table already exists
188  if (part == G4Electron::Electron())
189  {
190  if (XSTableElectron->count(theKey)) //table already built
191  return;
192  }
193  if (part == G4Positron::Positron())
194  {
195  if (XSTablePositron->count(theKey)) //table already built
196  return;
197  }
198 
199  //check if the material has been built
200  if (!(theDeltaTable->count(mat)))
201  BuildDeltaTable(mat);
202 
203 
204  //Tables have been already created (checked by GetCrossSectionTableForCouple)
205  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
206  size_t numberOfOscillators = theTable->size();
207 
208  if (energyGrid->GetVectorLength() != nBins)
209  {
210  G4ExceptionDescription ed;
211  ed << "Energy Grid looks not initialized" << G4endl;
212  ed << nBins << " " << energyGrid->GetVectorLength() << G4endl;
213  G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
214  "em2030",FatalException,ed);
215  }
216 
217  G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins,numberOfOscillators);
218 
219  //loop on the energy grid
220  for (size_t bin=0;bin<nBins;bin++)
221  {
222  G4double energy = energyGrid->GetLowEdgeEnergy(bin);
223  G4double XH0=0, XH1=0, XH2=0;
224  G4double XS0=0, XS1=0, XS2=0;
225 
226  //oscillator loop
227  for (size_t iosc=0;iosc<numberOfOscillators;iosc++)
228  {
229  G4DataVector* tempStorage = 0;
230 
231  G4PenelopeOscillator* theOsc = (*theTable)[iosc];
232  G4double delta = GetDensityCorrection(mat,energy);
233  if (part == G4Electron::Electron())
234  tempStorage = ComputeShellCrossSectionsElectron(theOsc,energy,cut,delta);
235  else if (part == G4Positron::Positron())
236  tempStorage = ComputeShellCrossSectionsPositron(theOsc,energy,cut,delta);
237  //check results are all right
238  if (!tempStorage)
239  {
240  G4ExceptionDescription ed;
241  ed << "Problem in calculating the shell XS for shell # "
242  << iosc << G4endl;
243  G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
244  "em2031",FatalException,ed);
245  delete XSEntry;
246  return;
247  }
248  if (tempStorage->size() != 6)
249  {
250  G4ExceptionDescription ed;
251  ed << "Problem in calculating the shell XS " << G4endl;
252  ed << "Result has dimension " << tempStorage->size() << " instead of 6" << G4endl;
253  G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()",
254  "em2031",FatalException,ed);
255  }
256  G4double stre = theOsc->GetOscillatorStrength();
257 
258  XH0 += stre*(*tempStorage)[0];
259  XH1 += stre*(*tempStorage)[1];
260  XH2 += stre*(*tempStorage)[2];
261  XS0 += stre*(*tempStorage)[3];
262  XS1 += stre*(*tempStorage)[4];
263  XS2 += stre*(*tempStorage)[5];
264  XSEntry->AddShellCrossSectionPoint(bin,iosc,energy,stre*(*tempStorage)[0]);
265  if (tempStorage)
266  {
267  delete tempStorage;
268  tempStorage = 0;
269  }
270  }
271  XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2);
272  }
273  //Do (only once) the final normalization
274  XSEntry->NormalizeShellCrossSections();
275 
276  //Insert in the appropriate table
277  if (part == G4Electron::Electron())
278  XSTableElectron->insert(std::make_pair(theKey,XSEntry));
279  else if (part == G4Positron::Positron())
280  XSTablePositron->insert(std::make_pair(theKey,XSEntry));
281  else
282  delete XSEntry;
283 
284  return;
285 }
286 
287 
288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
289 
290 G4double G4PenelopeIonisationXSHandler::GetDensityCorrection(const G4Material* mat,
291  const G4double energy) const
292 {
293  G4double result = 0;
294  if (!theDeltaTable)
295  {
296  G4Exception("G4PenelopeIonisationXSHandler::GetDensityCorrection()",
297  "em2032",FatalException,
298  "Delta Table not initialized. Was Initialise() run?");
299  return 0;
300  }
301  if (energy <= 0*eV)
302  {
303  G4cout << "G4PenelopeIonisationXSHandler::GetDensityCorrection()" << G4endl;
304  G4cout << "Invalid energy " << energy/eV << " eV " << G4endl;
305  return 0;
306  }
307  G4double logene = std::log(energy);
308 
309  if (theDeltaTable->count(mat))
310  {
311  const 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 }
G4PenelopeIonisationXSHandler::XSTableElectron
std::map< std::pair< const G4Material *, G4double >, G4PenelopeCrossSection * > * XSTableElectron
Definition: G4PenelopeIonisationXSHandler.hh:103
G4PenelopeOscillatorManager::GetOscillatorTableIonisation
G4PenelopeOscillatorTable * GetOscillatorTableIonisation(const G4Material *)
Definition: G4PenelopeOscillatorManager.cc:316
G4PhysicsFreeVector::PutValue
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
Definition: G4PhysicsFreeVector.cc:94
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4PenelopeOscillatorManager::GetPlasmaEnergySquared
G4double GetPlasmaEnergySquared(const G4Material *)
Returns the squared plasma energy.
Definition: G4PenelopeOscillatorManager.cc:1170
G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsElectron
G4DataVector * ComputeShellCrossSectionsElectron(G4PenelopeOscillator *, G4double energy, G4double cut, G4double delta)
Definition: G4PenelopeIonisationXSHandler.cc:430
G4PenelopeIonisationXSHandler::XSTablePositron
std::map< std::pair< const G4Material *, G4double >, G4PenelopeCrossSection * > * XSTablePositron
Definition: G4PenelopeIonisationXSHandler.hh:104
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:178
G4PenelopeCrossSection::AddShellCrossSectionPoint
void AddShellCrossSectionPoint(size_t binNumber, size_t shellID, G4double energy, G4double xs)
Definition: G4PenelopeCrossSection.cc:183
G4PhysicsLogVector.hh
G4PenelopeOscillatorManager.hh
G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple
const G4PenelopeCrossSection * GetCrossSectionTableForCouple(const G4ParticleDefinition *, const G4Material *, const G4double cut) const
Returns the table of cross sections for the given particle, given material and given cut as a G4Penel...
Definition: G4PenelopeIonisationXSHandler.cc:115
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:153
G4PenelopeCrossSection.hh
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:159
G4PenelopeCrossSection::AddCrossSectionPoint
void AddCrossSectionPoint(size_t binNumber, G4double energy, G4double XH0, G4double XH1, G4double XH2, G4double XS0, G4double XS1, G4double XS2)
Public interface for the master thread.
Definition: G4PenelopeCrossSection.cc:123
G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsPositron
G4DataVector * ComputeShellCrossSectionsPositron(G4PenelopeOscillator *, G4double energy, G4double cut, G4double delta)
Definition: G4PenelopeIonisationXSHandler.cc:562
G4PenelopeIonisationXSHandler::theDeltaTable
std::map< const G4Material *, G4PhysicsFreeVector * > * theDeltaTable
Definition: G4PenelopeIonisationXSHandler.hh:107
G4PenelopeOscillator.hh
G4cout
G4GLOB_DLL std::ostream G4cout
G4PenelopeIonisationXSHandler::G4PenelopeIonisationXSHandler
G4PenelopeIonisationXSHandler(size_t nBins=200)
Constructor. nBins is the number of intervals in the energy grid. By default the energy grid goes fro...
Definition: G4PenelopeIonisationXSHandler.cc:47
G4bool
bool G4bool
Definition: G4Types.hh:79
tab
std::ostream & tab(std::ostream &)
Definition: CCalutils.cc:89
G4Positron.hh
GeV
static const double GeV
Definition: G4SIunits.hh:214
G4ParticleDefinition.hh
G4PenelopeOscillatorManager::GetOscillatorManager
static G4PenelopeOscillatorManager * GetOscillatorManager()
Definition: G4PenelopeOscillatorManager.cc:89
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
G4PenelopeIonisationXSHandler::oscManager
G4PenelopeOscillatorManager * oscManager
Definition: G4PenelopeIonisationXSHandler.hh:100
eV
static const double eV
Definition: G4SIunits.hh:212
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4Electron.hh
pi
static const double pi
Definition: G4SIunits.hh:74
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:290
FatalException
Definition: G4ExceptionSeverity.hh:60
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
G4PhysicsFreeVector.hh
keV
static const double keV
Definition: G4SIunits.hh:213
G4PenelopeIonisationXSHandler::BuildXSTable
void BuildXSTable(const G4Material *, G4double cut, const G4ParticleDefinition *, G4bool isMaster=true)
This can be inkoved only by the master.
Definition: G4PenelopeIonisationXSHandler.cc:164
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4PenelopeOscillatorManager::GetTotalZ
G4double GetTotalZ(const G4Material *)
These are cumulative for the molecule Returns the total Z for the molecule.
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