ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeBremsstrahlungModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PenelopeBremsstrahlungModel.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 // 23 Nov 2010 L Pandola First complete implementation, Penelope v2008
32 // 24 May 2011 L. Pandola Renamed (make default Penelope)
33 // 13 Mar 2012 L. Pandola Updated the interface for the angular generator
34 // 18 Jul 2012 L. Pandola Migrate to the new interface of the angular generator, which
35 // now provides the G4ThreeVector and takes care of rotation
36 // 02 Oct 2013 L. Pandola Migrated to MT
37 // 17 Oct 2013 L. Pandola Partially revert the MT migration: the angular generator is
38 // kept as thread-local, and created/managed by the workers.
39 //
40 
41 #include "G4PenelopeBremsstrahlungModel.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4PenelopeBremsstrahlungFS.hh"
45 #include "G4PenelopeBremsstrahlungAngular.hh"
46 #include "G4ParticleDefinition.hh"
47 #include "G4MaterialCutsCouple.hh"
48 #include "G4ProductionCutsTable.hh"
49 #include "G4DynamicParticle.hh"
50 #include "G4Gamma.hh"
51 #include "G4Electron.hh"
52 #include "G4Positron.hh"
53 #include "G4PenelopeOscillatorManager.hh"
54 #include "G4PenelopeCrossSection.hh"
55 #include "G4PhysicsFreeVector.hh"
56 #include "G4PhysicsLogVector.hh"
57 #include "G4PhysicsTable.hh"
58 #include "G4Exp.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 namespace { G4Mutex PenelopeBremsstrahlungModelMutex = G4MUTEX_INITIALIZER; }
62 
63 G4PenelopeBremsstrahlungModel::G4PenelopeBremsstrahlungModel(const G4ParticleDefinition* part,
64  const G4String& nam)
65  :G4VEmModel(nam),fParticleChange(0),fParticle(0),
66  isInitialised(false),energyGrid(0),
67  XSTableElectron(0),XSTablePositron(0),fPenelopeFSHelper(0),
68  fPenelopeAngular(0),fLocalTable(false)
69 
70 {
71  fIntrinsicLowEnergyLimit = 100.0*eV;
72  fIntrinsicHighEnergyLimit = 100.0*GeV;
73  nBins = 200;
74 
75  if (part)
76  SetParticle(part);
77 
78  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
79  //
80  oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
81  //
82  verboseLevel= 0;
83 
84  // Verbosity scale:
85  // 0 = nothing
86  // 1 = warning for energy non-conservation
87  // 2 = details of energy budget
88  // 3 = calculation of cross sections, file openings, sampling of atoms
89  // 4 = entering in methods
90 
91  // Atomic deexcitation model activated by default
92  SetDeexcitationFlag(true);
93 
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 
98 G4PenelopeBremsstrahlungModel::~G4PenelopeBremsstrahlungModel()
99 {
100  if (IsMaster() || fLocalTable)
101  {
102  ClearTables();
103  if (fPenelopeFSHelper)
104  delete fPenelopeFSHelper;
105  }
106  // This is thread-local at the moment
107  if (fPenelopeAngular)
108  delete fPenelopeAngular;
109 }
110 
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112 
113 void G4PenelopeBremsstrahlungModel::Initialise(const G4ParticleDefinition* part,
114  const G4DataVector& theCuts)
115 {
116  if (verboseLevel > 3)
117  G4cout << "Calling G4PenelopeBremsstrahlungModel::Initialise()" << G4endl;
118 
119  SetParticle(part);
120 
121  if (IsMaster() && part == fParticle)
122  {
123 
124  if (!fPenelopeFSHelper)
125  fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS(verboseLevel);
126  if (!fPenelopeAngular)
127  fPenelopeAngular = new G4PenelopeBremsstrahlungAngular();
128  //Clear and re-build the tables
129  ClearTables();
130 
131  //forces the cleaning of tables, in this specific case
132  if (fPenelopeAngular)
133  fPenelopeAngular->Initialize();
134 
135  //Set the number of bins for the tables. 20 points per decade
136  nBins = (size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit()));
137  nBins = std::max(nBins,(size_t)100);
138  energyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
139  HighEnergyLimit(),
140  nBins-1); //one hidden bin is added
141 
142 
143  XSTableElectron = new
144  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
145  XSTablePositron = new
146  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
147 
148  G4ProductionCutsTable* theCoupleTable =
149  G4ProductionCutsTable::GetProductionCutsTable();
150 
151  //Build tables for all materials
152  for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
153  {
154  const G4Material* theMat =
155  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
156  //Forces the building of the helper tables
157  fPenelopeFSHelper->BuildScaledXSTable(theMat,theCuts.at(i),IsMaster());
158  fPenelopeAngular->PrepareTables(theMat,IsMaster());
159  BuildXSTable(theMat,theCuts.at(i));
160 
161  }
162 
163 
164  if (verboseLevel > 2) {
165  G4cout << "Penelope Bremsstrahlung model v2008 is initialized " << G4endl
166  << "Energy range: "
167  << LowEnergyLimit() / keV << " keV - "
168  << HighEnergyLimit() / GeV << " GeV."
169  << G4endl;
170  }
171  }
172 
173  if(isInitialised) return;
174  fParticleChange = GetParticleChangeForLoss();
175  isInitialised = true;
176 }
177 
178 
179 void G4PenelopeBremsstrahlungModel::InitialiseLocal(const G4ParticleDefinition* part,
180  G4VEmModel *masterModel)
181 {
182  if (verboseLevel > 3)
183  G4cout << "Calling G4PenelopeBremsstrahlungModel::InitialiseLocal()" << G4endl;
184 
185  //
186  //Check that particle matches: one might have multiple master models (e.g.
187  //for e+ and e-).
188  //
189  if (part == fParticle)
190  {
191  //Get the const table pointers from the master to the workers
192  const G4PenelopeBremsstrahlungModel* theModel =
193  static_cast<G4PenelopeBremsstrahlungModel*> (masterModel);
194 
195  //Copy pointers to the data tables
196  energyGrid = theModel->energyGrid;
197  XSTableElectron = theModel->XSTableElectron;
198  XSTablePositron = theModel->XSTablePositron;
199  fPenelopeFSHelper = theModel->fPenelopeFSHelper;
200 
201  //fPenelopeAngular = theModel->fPenelopeAngular;
202 
203  //created in each thread and initialized.
204  if (!fPenelopeAngular)
205  fPenelopeAngular = new G4PenelopeBremsstrahlungAngular();
206  //forces the cleaning of tables, in this specific case
207  if (fPenelopeAngular)
208  fPenelopeAngular->Initialize();
209 
210  G4ProductionCutsTable* theCoupleTable =
211  G4ProductionCutsTable::GetProductionCutsTable();
212  //Build tables for all materials
213  for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
214  {
215  const G4Material* theMat =
216  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
217  fPenelopeAngular->PrepareTables(theMat,IsMaster());
218  }
219 
220 
221  //copy the data
222  nBins = theModel->nBins;
223 
224  //Same verbosity for all workers, as the master
225  verboseLevel = theModel->verboseLevel;
226  }
227 
228  return;
229 }
230 
231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
232 
233 G4double G4PenelopeBremsstrahlungModel::CrossSectionPerVolume(const G4Material* material,
234  const G4ParticleDefinition* theParticle,
235  G4double energy,
236  G4double cutEnergy,
237  G4double)
238 {
239  //
240  if (verboseLevel > 3)
241  G4cout << "Calling CrossSectionPerVolume() of G4PenelopeBremsstrahlungModel" << G4endl;
242 
243  SetupForMaterial(theParticle, material, energy);
244 
245  G4double crossPerMolecule = 0.;
246 
247  G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
248  cutEnergy);
249 
250  if (theXS)
251  crossPerMolecule = theXS->GetHardCrossSection(energy);
252 
253  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
254  G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
255 
256  if (verboseLevel > 3)
257  G4cout << "Material " << material->GetName() << " has " << atPerMol <<
258  "atoms per molecule" << G4endl;
259 
260  G4double moleculeDensity = 0.;
261  if (atPerMol)
262  moleculeDensity = atomDensity/atPerMol;
263 
264  G4double crossPerVolume = crossPerMolecule*moleculeDensity;
265 
266  if (verboseLevel > 2)
267  {
268  G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
269  G4cout << "Mean free path for gamma emission > " << cutEnergy/keV << " keV at " <<
270  energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl;
271  }
272 
273  return crossPerVolume;
274 }
275 
276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
277 
278 //This is a dummy method. Never inkoved by the tracking, it just issues
279 //a warning if one tries to get Cross Sections per Atom via the
280 //G4EmCalculator.
281 G4double G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
282  G4double,
283  G4double,
284  G4double,
285  G4double,
286  G4double)
287 {
288  G4cout << "*** G4PenelopeBremsstrahlungModel -- WARNING ***" << G4endl;
289  G4cout << "Penelope Bremsstrahlung model v2008 does not calculate cross section _per atom_ " << G4endl;
290  G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
291  G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
292  return 0;
293 }
294 
295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296 
297 G4double G4PenelopeBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* material,
298  const G4ParticleDefinition* theParticle,
299  G4double kineticEnergy,
300  G4double cutEnergy)
301 {
302  if (verboseLevel > 3)
303  G4cout << "Calling ComputeDEDX() of G4PenelopeBremsstrahlungModel" << G4endl;
304 
305  G4PenelopeCrossSection* theXS = GetCrossSectionTableForCouple(theParticle,material,
306  cutEnergy);
307  G4double sPowerPerMolecule = 0.0;
308  if (theXS)
309  sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy);
310 
311 
312  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
313  G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
314 
315  G4double moleculeDensity = 0.;
316  if (atPerMol)
317  moleculeDensity = atomDensity/atPerMol;
318 
319  G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
320 
321  if (verboseLevel > 2)
322  {
323  G4cout << "G4PenelopeBremsstrahlungModel " << G4endl;
324  G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
325  kineticEnergy/keV << " keV = " <<
326  sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl;
327  }
328  return sPowerPerVolume;
329 }
330 
331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
332 
333 void G4PenelopeBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>*fvect,
334  const G4MaterialCutsCouple* couple,
335  const G4DynamicParticle* aDynamicParticle,
336  G4double cutG,
337  G4double)
338 {
339  if (verboseLevel > 3)
340  G4cout << "Calling SampleSecondaries() of G4PenelopeBremsstrahlungModel" << G4endl;
341 
342  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
343  const G4Material* material = couple->GetMaterial();
344 
345  if (kineticEnergy <= fIntrinsicLowEnergyLimit)
346  {
347  fParticleChange->SetProposedKineticEnergy(0.);
348  fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
349  return ;
350  }
351 
352  G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
353  //This is the momentum
354  G4ThreeVector initialMomentum = aDynamicParticle->GetMomentum();
355 
356  //Not enough energy to produce a secondary! Return with nothing happened
357  if (kineticEnergy < cutG)
358  return;
359 
360  if (verboseLevel > 3)
361  G4cout << "Going to sample gamma energy for: " <<material->GetName() << " " <<
362  "energy = " << kineticEnergy/keV << ", cut = " << cutG/keV << G4endl;
363 
364  //Sample gamma's energy according to the spectrum
365  G4double gammaEnergy =
366  fPenelopeFSHelper->SampleGammaEnergy(kineticEnergy,material,cutG);
367 
368  if (verboseLevel > 3)
369  G4cout << "Sampled gamma energy: " << gammaEnergy/keV << " keV" << G4endl;
370 
371  //Now sample the direction for the Gamma. Notice that the rotation is already done
372  //Z is unused here, I plug 0. The information is in the material pointer
373  G4ThreeVector gammaDirection1 =
374  fPenelopeAngular->SampleDirection(aDynamicParticle,gammaEnergy,0,material);
375 
376  if (verboseLevel > 3)
377  G4cout << "Sampled cosTheta for e-: " << gammaDirection1.cosTheta() << G4endl;
378 
379  G4double residualPrimaryEnergy = kineticEnergy-gammaEnergy;
380  if (residualPrimaryEnergy < 0)
381  {
382  //Ok we have a problem, all energy goes with the gamma
383  gammaEnergy += residualPrimaryEnergy;
384  residualPrimaryEnergy = 0.0;
385  }
386 
387  //Produce final state according to momentum conservation
388  G4ThreeVector particleDirection1 = initialMomentum - gammaEnergy*gammaDirection1;
389  particleDirection1 = particleDirection1.unit(); //normalize
390 
391  //Update the primary particle
392  if (residualPrimaryEnergy > 0.)
393  {
394  fParticleChange->ProposeMomentumDirection(particleDirection1);
395  fParticleChange->SetProposedKineticEnergy(residualPrimaryEnergy);
396  }
397  else
398  {
399  fParticleChange->SetProposedKineticEnergy(0.);
400  }
401 
402  //Now produce the photon
403  G4DynamicParticle* theGamma = new G4DynamicParticle(G4Gamma::Gamma(),
404  gammaDirection1,
405  gammaEnergy);
406  fvect->push_back(theGamma);
407 
408  if (verboseLevel > 1)
409  {
410  G4cout << "-----------------------------------------------------------" << G4endl;
411  G4cout << "Energy balance from G4PenelopeBremsstrahlung" << G4endl;
412  G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
413  G4cout << "-----------------------------------------------------------" << G4endl;
414  G4cout << "Outgoing primary energy: " << residualPrimaryEnergy/keV << " keV" << G4endl;
415  G4cout << "Bremsstrahlung photon " << gammaEnergy/keV << " keV" << G4endl;
416  G4cout << "Total final state: " << (residualPrimaryEnergy+gammaEnergy)/keV
417  << " keV" << G4endl;
418  G4cout << "-----------------------------------------------------------" << G4endl;
419  }
420 
421  if (verboseLevel > 0)
422  {
423  G4double energyDiff = std::fabs(residualPrimaryEnergy+gammaEnergy-kineticEnergy);
424  if (energyDiff > 0.05*keV)
425  G4cout << "Warning from G4PenelopeBremsstrahlung: problem with energy conservation: "
426  <<
427  (residualPrimaryEnergy+gammaEnergy)/keV <<
428  " keV (final) vs. " <<
429  kineticEnergy/keV << " keV (initial)" << G4endl;
430  }
431  return;
432 }
433 
434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
435 
436 void G4PenelopeBremsstrahlungModel::ClearTables()
437 {
438  if (!IsMaster() && !fLocalTable)
439  //Should not be here!
440  G4Exception("G4PenelopeBremsstrahlungModel::ClearTables()",
441  "em0100",FatalException,"Worker thread in this method");
442 
443  if (XSTableElectron)
444  {
445  for (auto& item : (*XSTableElectron))
446  delete item.second;
447  delete XSTableElectron;
448  XSTableElectron = nullptr;
449  }
450  if (XSTablePositron)
451  {
452  for (auto& item : (*XSTablePositron))
453  delete item.second;
454  delete XSTablePositron;
455  XSTablePositron = nullptr;
456  }
457  /*
458  if (energyGrid)
459  delete energyGrid;
460  */
461  if (fPenelopeFSHelper)
462  fPenelopeFSHelper->ClearTables(IsMaster());
463 
464  if (verboseLevel > 2)
465  G4cout << "G4PenelopeBremsstrahlungModel: cleared tables" << G4endl;
466 
467  return;
468 }
469 
470 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
471 
472 G4double G4PenelopeBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
473  const G4MaterialCutsCouple*)
474 {
475  return fIntrinsicLowEnergyLimit;
476 }
477 
478 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
479 
480 void G4PenelopeBremsstrahlungModel::BuildXSTable(const G4Material* mat,G4double cut)
481 {
482  if (!IsMaster() && !fLocalTable)
483  //Should not be here!
484  G4Exception("G4PenelopeBremsstrahlungModel::BuildXSTable()",
485  "em0100",FatalException,"Worker thread in this method");
486 
487  //The key of the map
488  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
489 
490  //The table already exists
491  if (XSTableElectron->count(theKey) && XSTablePositron->count(theKey))
492  return;
493 
494  //
495  //This method fills the G4PenelopeCrossSection containers for electrons or positrons
496  //and for the given material/cut couple.
497  //Equivalent of subroutines EBRaT and PINaT of Penelope
498  //
499  if (verboseLevel > 2)
500  {
501  G4cout << "G4PenelopeBremsstrahlungModel: going to build cross section table " << G4endl;
502  G4cout << "for e+/e- in " << mat->GetName() << " for Ecut(gamma)= " <<
503  cut/keV << " keV " << G4endl;
504  }
505 
506  //Tables have been already created (checked by GetCrossSectionTableForCouple)
507  if (energyGrid->GetVectorLength() != nBins)
508  {
509  G4ExceptionDescription ed;
510  ed << "Energy Grid looks not initialized" << G4endl;
511  ed << nBins << " " << energyGrid->GetVectorLength() << G4endl;
512  G4Exception("G4PenelopeBremsstrahlungModel::BuildXSTable()",
513  "em2016",FatalException,ed);
514  }
515 
516  G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins);
517  G4PenelopeCrossSection* XSEntry2 = new G4PenelopeCrossSection(nBins);
518 
519  const G4PhysicsTable* table = fPenelopeFSHelper->GetScaledXSTable(mat,cut);
520 
521 
522  //loop on the energy grid
523  for (size_t bin=0;bin<nBins;bin++)
524  {
525  G4double energy = energyGrid->GetLowEdgeEnergy(bin);
526  G4double XH0=0, XH1=0, XH2=0;
527  G4double XS0=0, XS1=0, XS2=0;
528 
529  //Global xs factor
530  G4double fact = fPenelopeFSHelper->GetEffectiveZSquared(mat)*
531  ((energy+electron_mass_c2)*(energy+electron_mass_c2)/
532  (energy*(energy+2.0*electron_mass_c2)));
533 
534  G4double restrictedCut = cut/energy;
535 
536  //Now I need the dSigma/dX profile - interpolated on energy - for
537  //the 32-point x grid. Interpolation is log-log
538  size_t nBinsX = fPenelopeFSHelper->GetNBinsX();
539  G4double* tempData = new G4double[nBinsX];
540  G4double logene = std::log(energy);
541  for (size_t ix=0;ix<nBinsX;ix++)
542  {
543  //find dSigma/dx for the given E. X belongs to the 32-point grid.
544  G4double val = (*table)[ix]->Value(logene);
545  tempData[ix] = G4Exp(val); //back to the real value!
546  }
547 
548  G4double XH0A = 0.;
549  if (restrictedCut <= 1) //calculate only if we are above threshold!
550  XH0A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,-1) -
551  fPenelopeFSHelper->GetMomentumIntegral(tempData,restrictedCut,-1);
552  G4double XS1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
553  restrictedCut,0);
554  G4double XS2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,
555  restrictedCut,1);
556  G4double XH1A=0, XH2A=0;
557  if (restrictedCut <=1)
558  {
559  XH1A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,0) -
560  XS1A;
561  XH2A = fPenelopeFSHelper->GetMomentumIntegral(tempData,1.0,1) -
562  XS2A;
563  }
564  delete[] tempData;
565 
566  XH0 = XH0A*fact;
567  XS1 = XS1A*fact*energy;
568  XH1 = XH1A*fact*energy;
569  XS2 = XS2A*fact*energy*energy;
570  XH2 = XH2A*fact*energy*energy;
571 
572  XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2);
573 
574  //take care of positrons
575  G4double posCorrection = GetPositronXSCorrection(mat,energy);
576  XSEntry2->AddCrossSectionPoint(bin,energy,XH0*posCorrection,
577  XH1*posCorrection,
578  XH2*posCorrection,
579  XS0,
580  XS1*posCorrection,
581  XS2*posCorrection);
582  }
583 
584  //Insert in the appropriate table
585  XSTableElectron->insert(std::make_pair(theKey,XSEntry));
586  XSTablePositron->insert(std::make_pair(theKey,XSEntry2));
587 
588  return;
589 }
590 
591 
592 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
593 
594 G4PenelopeCrossSection*
595 G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple(const G4ParticleDefinition* part,
596  const G4Material* mat,
597  G4double cut)
598 {
599  if (part != G4Electron::Electron() && part != G4Positron::Positron())
600  {
601  G4ExceptionDescription ed;
602  ed << "Invalid particle: " << part->GetParticleName() << G4endl;
603  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
604  "em0001",FatalException,ed);
605  return nullptr;
606  }
607 
608  if (part == G4Electron::Electron())
609  {
610  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
611  //not invoked
612  if (!XSTableElectron)
613  {
614  //create a **thread-local** version of the table. Used only for G4EmCalculator and
615  //Unit Tests
616  G4String excep = "The Cross Section Table for e- was not initialized correctly!";
617  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
618  "em2013",JustWarning,excep);
619  fLocalTable = true;
620  XSTableElectron = new
621  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
622  if (!energyGrid)
623  energyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
624  HighEnergyLimit(),
625  nBins-1); //one hidden bin is added
626  if (!fPenelopeFSHelper)
627  fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS(verboseLevel);
628  }
629  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
630  if (XSTableElectron->count(theKey)) //table already built
631  return XSTableElectron->find(theKey)->second;
632  else
633  {
634  //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
635  //not filled up. This can happen in a UnitTest or via G4EmCalculator
636  if (verboseLevel > 0)
637  {
638  //G4Exception (warning) is issued only in verbose mode
639  G4ExceptionDescription ed;
640  ed << "Unable to find e- table for " << mat->GetName() << " at Ecut(gamma)= "
641  << cut/keV << " keV " << G4endl;
642  ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
643  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
644  "em2009",JustWarning,ed);
645  }
646  //protect file reading via autolock
647  G4AutoLock lock(&PenelopeBremsstrahlungModelMutex);
648  fPenelopeFSHelper->BuildScaledXSTable(mat,cut,true); //pretend to be a master
649  BuildXSTable(mat,cut);
650  lock.unlock();
651  //now it should be ok
652  return XSTableElectron->find(theKey)->second;
653  }
654 
655  }
656  if (part == G4Positron::Positron())
657  {
658  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
659  //not invoked
660  if (!XSTablePositron)
661  {
662  G4String excep = "The Cross Section Table for e+ was not initialized correctly!";
663  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
664  "em2013",JustWarning,excep);
665  fLocalTable = true;
666  XSTablePositron = new
667  std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>;
668  if (!energyGrid)
669  energyGrid = new G4PhysicsLogVector(LowEnergyLimit(),
670  HighEnergyLimit(),
671  nBins-1); //one hidden bin is added
672  if (!fPenelopeFSHelper)
673  fPenelopeFSHelper = new G4PenelopeBremsstrahlungFS(verboseLevel);
674  }
675  std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut);
676  if (XSTablePositron->count(theKey)) //table already built
677  return XSTablePositron->find(theKey)->second;
678  else
679  {
680  //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
681  //not filled up. This can happen in a UnitTest or via G4EmCalculator
682  if (verboseLevel > 0)
683  {
684  //Issue a G4Exception (warning) only in verbose mode
685  G4ExceptionDescription ed;
686  ed << "Unable to find e+ table for " << mat->GetName() << " at Ecut(gamma)= "
687  << cut/keV << " keV " << G4endl;
688  ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
689  G4Exception("G4PenelopeBremsstrahlungModel::GetCrossSectionTableForCouple()",
690  "em2009",JustWarning,ed);
691  }
692  //protect file reading via autolock
693  G4AutoLock lock(&PenelopeBremsstrahlungModelMutex);
694  fPenelopeFSHelper->BuildScaledXSTable(mat,cut,true); //pretend to be a master
695  BuildXSTable(mat,cut);
696  lock.unlock();
697  //now it should be ok
698  return XSTablePositron->find(theKey)->second;
699  }
700  }
701  return nullptr;
702 }
703 
704 
705 
706 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
707 
708 G4double G4PenelopeBremsstrahlungModel::GetPositronXSCorrection(const G4Material* mat,
709  G4double energy)
710 {
711  //The electron-to-positron correction factor is set equal to the ratio of the
712  //radiative stopping powers for positrons and electrons, which has been calculated
713  //by Kim et al. (1986) (cf. Berger and Seltzer, 1982). Here, it is used an
714  //analytical approximation which reproduces the tabulated values with 0.5%
715  //accuracy
716  G4double t=std::log(1.0+1e6*energy/
717  (electron_mass_c2*fPenelopeFSHelper->GetEffectiveZSquared(mat)));
718  G4double corr = 1.0-G4Exp(-t*(1.2359e-1-t*(6.1274e-2-t*
719  (3.1516e-2-t*(7.7446e-3-t*(1.0595e-3-t*
720  (7.0568e-5-t*
721  1.8080e-6)))))));
722  return corr;
723 }
724 
725 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
726 
727 void G4PenelopeBremsstrahlungModel::SetParticle(const G4ParticleDefinition* p)
728 {
729  if(!fParticle) {
730  fParticle = p;
731  }
732 }