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