| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4PenelopeBremsstrahlungModel.cc,v 1.8 2010-11-25 09:44:05 pandola Exp $
>> 27 // GEANT4 tag $Name: not supported by cvs2svn $
26 // 28 //
27 // Author: Luciano Pandola 29 // Author: Luciano Pandola
28 // <<
29 // History: <<
30 // -------- 30 // --------
31 // 23 Nov 2010 L Pandola First complete i << 31 // 05 Dec 2008 L Pandola Migration from process to model
32 // 24 May 2011 L. Pandola Renamed (make de << 32 // 25 Mar 2008 L Pandola Fixed .unit() call
33 // 13 Mar 2012 L. Pandola Updated the inte << 33 // 16 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
34 // 18 Jul 2012 L. Pandola Migrate to the n << 34 // - apply internal high-energy limit only in constructor
35 // now provides the << 35 // - do not apply low-energy limit (default is 0)
36 // 02 Oct 2013 L. Pandola Migrated to MT << 36 // - added MinEnergyCut method
37 // 17 Oct 2013 L. Pandola Partially revert << 37 // - do not change track status
38 // kept as thread- << 38 // 14 May 2009 L Pandola Explicitely set to zero pointers deleted in
>> 39 // Initialise(), since they are checked later on
39 // 40 //
40 <<
41 #include "G4PenelopeBremsstrahlungModel.hh" 41 #include "G4PenelopeBremsstrahlungModel.hh"
42 #include "G4PhysicalConstants.hh" << 42 #include "G4PenelopeBremsstrahlungContinuous.hh"
43 #include "G4SystemOfUnits.hh" <<
44 #include "G4PenelopeBremsstrahlungFS.hh" <<
45 #include "G4PenelopeBremsstrahlungAngular.hh" 43 #include "G4PenelopeBremsstrahlungAngular.hh"
46 #include "G4ParticleDefinition.hh" << 44 #include "G4eBremsstrahlungSpectrum.hh"
47 #include "G4MaterialCutsCouple.hh" << 45 #include "G4CrossSectionHandler.hh"
48 #include "G4ProductionCutsTable.hh" << 46 #include "G4VEMDataSet.hh"
49 #include "G4DynamicParticle.hh" << 47 #include "G4DataVector.hh"
50 #include "G4Gamma.hh" <<
51 #include "G4Electron.hh" <<
52 #include "G4Positron.hh" 48 #include "G4Positron.hh"
53 #include "G4PenelopeOscillatorManager.hh" << 49 #include "G4Electron.hh"
54 #include "G4PenelopeCrossSection.hh" << 50 #include "G4Gamma.hh"
55 #include "G4PhysicsFreeVector.hh" << 51 #include "G4MaterialCutsCouple.hh"
56 #include "G4PhysicsLogVector.hh" << 52 #include "G4LogLogInterpolation.hh"
57 #include "G4PhysicsTable.hh" <<
58 #include "G4Exp.hh" <<
59 53
60 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 namespace { G4Mutex PenelopeBremsstrahlungMod <<
62 55
63 G4PenelopeBremsstrahlungModel::G4PenelopeBrems << 56 G4PenelopeBremsstrahlungModel::G4PenelopeBremsstrahlungModel(const G4ParticleDefinition*,
64 const G4String& nam) 57 const G4String& nam)
65 :G4VEmModel(nam),fParticleChange(nullptr),fP << 58 :G4VEmModel(nam),isInitialised(false),energySpectrum(0),
66 fPenelopeFSHelper(nullptr),fPenelopeAngular << 59 angularData(0),stoppingPowerData(0),crossSectionHandler(0)
67 fXSTableElectron(nullptr),fXSTablePositron( <<
68 fIsInitialised(false),fLocalTable(false) <<
69 <<
70 { 60 {
71 fIntrinsicLowEnergyLimit = 100.0*eV; 61 fIntrinsicLowEnergyLimit = 100.0*eV;
72 fIntrinsicHighEnergyLimit = 100.0*GeV; 62 fIntrinsicHighEnergyLimit = 100.0*GeV;
73 nBins = 200; << 63 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
74 <<
75 if (part) <<
76 SetParticle(part); <<
77 <<
78 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 64 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
79 // 65 //
80 fOscManager = G4PenelopeOscillatorManager::G << 66
81 // << 67 verboseLevel= 0;
82 fVerboseLevel= 0; << 68
83 // Verbosity scale: 69 // Verbosity scale:
84 // 0 = nothing 70 // 0 = nothing
85 // 1 = warning for energy non-conservation 71 // 1 = warning for energy non-conservation
86 // 2 = details of energy budget 72 // 2 = details of energy budget
87 // 3 = calculation of cross sections, file o 73 // 3 = calculation of cross sections, file openings, sampling of atoms
88 // 4 = entering in methods 74 // 4 = entering in methods
89 << 75
90 // Atomic deexcitation model activated by de << 76 //These vectors do not change when materials or cut change.
91 SetDeexcitationFlag(true); << 77 //Therefore I can read it at the constructor
92 << 78 angularData = new std::map<G4int,G4PenelopeBremsstrahlungAngular*>;
>> 79
>> 80 //These data do not depend on materials and cuts.
>> 81 G4DataVector eBins;
>> 82
>> 83 eBins.push_back(1.0e-12);
>> 84 eBins.push_back(0.05);
>> 85 eBins.push_back(0.075);
>> 86 eBins.push_back(0.1);
>> 87 eBins.push_back(0.125);
>> 88 eBins.push_back(0.15);
>> 89 eBins.push_back(0.2);
>> 90 eBins.push_back(0.25);
>> 91 eBins.push_back(0.3);
>> 92 eBins.push_back(0.35);
>> 93 eBins.push_back(0.40);
>> 94 eBins.push_back(0.45);
>> 95 eBins.push_back(0.50);
>> 96 eBins.push_back(0.55);
>> 97 eBins.push_back(0.60);
>> 98 eBins.push_back(0.65);
>> 99 eBins.push_back(0.70);
>> 100 eBins.push_back(0.75);
>> 101 eBins.push_back(0.80);
>> 102 eBins.push_back(0.85);
>> 103 eBins.push_back(0.90);
>> 104 eBins.push_back(0.925);
>> 105 eBins.push_back(0.95);
>> 106 eBins.push_back(0.97);
>> 107 eBins.push_back(0.99);
>> 108 eBins.push_back(0.995);
>> 109 eBins.push_back(0.999);
>> 110 eBins.push_back(0.9995);
>> 111 eBins.push_back(0.9999);
>> 112 eBins.push_back(0.99995);
>> 113 eBins.push_back(0.99999);
>> 114 eBins.push_back(1.0);
>> 115
>> 116 const G4String dataName("/penelope/br-sp-pen.dat");
>> 117 energySpectrum = new G4eBremsstrahlungSpectrum(eBins,dataName);
93 } 118 }
94 119
95 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 121
97 G4PenelopeBremsstrahlungModel::~G4PenelopeBrem 122 G4PenelopeBremsstrahlungModel::~G4PenelopeBremsstrahlungModel()
98 { 123 {
99 if (IsMaster() || fLocalTable) << 124 if (crossSectionHandler)
>> 125 delete crossSectionHandler;
>> 126
>> 127 if (energySpectrum)
>> 128 delete energySpectrum;
>> 129
>> 130 if (angularData)
>> 131 {
>> 132 std::map <G4int,G4PenelopeBremsstrahlungAngular*>::iterator i;
>> 133 for (i=angularData->begin();i != angularData->end();i++)
>> 134 if (i->second) delete i->second;
>> 135 delete angularData;
>> 136 }
>> 137
>> 138 if (stoppingPowerData)
100 { 139 {
101 ClearTables(); << 140 std::map <std::pair<G4int,G4double>,G4PenelopeBremsstrahlungContinuous*>::iterator j;
102 if (fPenelopeFSHelper) << 141 for (j=stoppingPowerData->begin();j != stoppingPowerData->end();j++)
103 delete fPenelopeFSHelper; << 142 if (j->second) delete j->second;
104 } << 143 delete stoppingPowerData;
105 // This is thread-local at the moment << 144 }
106 if (fPenelopeAngular) <<
107 delete fPenelopeAngular; <<
108 } 145 }
109 146
110 //....oooOO0OOooo........oooOO0OOooo........oo 147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 148
112 void G4PenelopeBremsstrahlungModel::Initialise << 149 void G4PenelopeBremsstrahlungModel::Initialise(const G4ParticleDefinition* particle,
113 c << 150 const G4DataVector& cuts)
114 { 151 {
115 if (fVerboseLevel > 3) << 152 if (verboseLevel > 3)
116 G4cout << "Calling G4PenelopeBremsstrahlun 153 G4cout << "Calling G4PenelopeBremsstrahlungModel::Initialise()" << G4endl;
>> 154
>> 155 // Delete everything, but angular data (do not depend on cuts)
>> 156 if (crossSectionHandler)
>> 157 {
>> 158 crossSectionHandler->Clear();
>> 159 delete crossSectionHandler;
>> 160 crossSectionHandler = 0;
>> 161 }
>> 162
>> 163 if (stoppingPowerData)
>> 164 {
>> 165 std::map <std::pair<G4int,G4double>,G4PenelopeBremsstrahlungContinuous*>::iterator j;
>> 166 for (j=stoppingPowerData->begin();j != stoppingPowerData->end();j++)
>> 167 if (j->second)
>> 168 {
>> 169 delete j->second;
>> 170 j->second = 0;
>> 171 }
>> 172 delete stoppingPowerData;
>> 173 stoppingPowerData = 0;
>> 174 }
>> 175
>> 176 crossSectionHandler = new G4CrossSectionHandler();
>> 177 crossSectionHandler->Clear();
>> 178 //
>> 179 if (particle==G4Electron::Electron())
>> 180 crossSectionHandler->LoadData("brem/br-cs-");
>> 181 else
>> 182 crossSectionHandler->LoadData("penelope/br-cs-pos-"); //cross section for positrons
>> 183
>> 184 //This is used to retrieve cross section values later on
>> 185 G4VEMDataSet* emdata =
>> 186 crossSectionHandler->BuildMeanFreePathForMaterials();
>> 187 //The method BuildMeanFreePathForMaterials() is required here only to force
>> 188 //the building of an internal table: the output pointer can be deleted
>> 189 delete emdata;
>> 190
>> 191 if (verboseLevel > 2)
>> 192 G4cout << "Loaded cross section files for PenelopeBremsstrahlungModel" << G4endl;
>> 193
>> 194 if (verboseLevel > 0) {
>> 195 G4cout << "Penelope Bremsstrahlung model is initialized " << G4endl
>> 196 << "Energy range: "
>> 197 << LowEnergyLimit() / keV << " keV - "
>> 198 << HighEnergyLimit() / GeV << " GeV"
>> 199 << G4endl;
>> 200 }
>> 201
>> 202 //This has to be invoked AFTER the crossSectionHandler has been created,
>> 203 //because it makes use of ComputeCrossSectionPerAtom()
>> 204 InitialiseElementSelectors(particle,cuts);
117 205
118 SetParticle(part); << 206 if(isInitialised) return;
119 <<
120 if (IsMaster() && part == fParticle) <<
121 { <<
122 if (!fPenelopeFSHelper) <<
123 fPenelopeFSHelper = new G4PenelopeBremsstrah <<
124 if (!fPenelopeAngular) <<
125 fPenelopeAngular = new G4PenelopeBremsstrahl <<
126 //Clear and re-build the tables <<
127 ClearTables(); <<
128 <<
129 //forces the cleaning of tables, in this <<
130 if (fPenelopeAngular) <<
131 fPenelopeAngular->Initialize(); <<
132 <<
133 //Set the number of bins for the tables. <<
134 nBins = (std::size_t) (20*std::log10(Hig <<
135 nBins = std::max(nBins,(std::size_t)100) <<
136 fEnergyGrid = new G4PhysicsLogVector(Low <<
137 HighEner <<
138 nBins-1) <<
139 <<
140 fXSTableElectron = new <<
141 std::map< std::pair<const G4Material*,G4doub <<
142 fXSTablePositron = new <<
143 std::map< std::pair<const G4Material*,G4doub <<
144 <<
145 G4ProductionCutsTable* theCoupleTable = <<
146 G4ProductionCutsTable::GetProductionCutsTabl <<
147 <<
148 //Build tables for all materials <<
149 for (G4int i=0;i<(G4int)theCoupleTable-> <<
150 { <<
151 const G4Material* theMat = <<
152 theCoupleTable->GetMaterialCutsCouple(i) <<
153 //Forces the building of the helper tables <<
154 fPenelopeFSHelper->BuildScaledXSTable(theM <<
155 fPenelopeAngular->PrepareTables(theMat,IsM <<
156 BuildXSTable(theMat,theCuts.at(i)); <<
157 <<
158 } <<
159 <<
160 if (fVerboseLevel > 2) { <<
161 G4cout << "Penelope Bremsstrahlung model v20 <<
162 << "Energy range: " <<
163 << LowEnergyLimit() / keV << " keV - <<
164 << HighEnergyLimit() / GeV << " GeV." <<
165 << G4endl; <<
166 } <<
167 } <<
168 <<
169 if(fIsInitialised) return; <<
170 fParticleChange = GetParticleChangeForLoss() 207 fParticleChange = GetParticleChangeForLoss();
171 fIsInitialised = true; << 208 isInitialised = true;
172 } 209 }
173 210
174 //....oooOO0OOooo........oooOO0OOooo........oo 211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
175 212
176 void G4PenelopeBremsstrahlungModel::Initialise << 213 G4double G4PenelopeBremsstrahlungModel::MinEnergyCut(const G4ParticleDefinition*,
177 G4VEmModel *masterModel) << 214 const G4MaterialCutsCouple*)
178 { 215 {
179 if (fVerboseLevel > 3) << 216 return 250.*eV;
180 G4cout << "Calling G4PenelopeBremsstrahlu <<
181 // <<
182 //Check that particle matches: one might hav <<
183 //for e+ and e-). <<
184 // <<
185 if (part == fParticle) <<
186 { <<
187 //Get the const table pointers from the <<
188 const G4PenelopeBremsstrahlungModel* the <<
189 static_cast<G4PenelopeBremsstrahlungModel*> <<
190 <<
191 //Copy pointers to the data tables <<
192 fEnergyGrid = theModel->fEnergyGrid; <<
193 fXSTableElectron = theModel->fXSTableEle <<
194 fXSTablePositron = theModel->fXSTablePos <<
195 fPenelopeFSHelper = theModel->fPenelopeF <<
196 <<
197 //created in each thread and initialized <<
198 if (!fPenelopeAngular) <<
199 fPenelopeAngular = new G4PenelopeBremsstrahl <<
200 //forces the cleaning of tables, in this <<
201 if (fPenelopeAngular) <<
202 fPenelopeAngular->Initialize(); <<
203 <<
204 G4ProductionCutsTable* theCoupleTable = <<
205 G4ProductionCutsTable::GetProductionCutsTabl <<
206 //Build tables for all materials <<
207 for (G4int i=0;i<(G4int)theCoupleTable-> <<
208 { <<
209 const G4Material* theMat = <<
210 theCoupleTable->GetMaterialCutsCouple(i) <<
211 fPenelopeAngular->PrepareTables(theMat,IsM <<
212 } <<
213 <<
214 //copy the data <<
215 nBins = theModel->nBins; <<
216 <<
217 //Same verbosity for all workers, as the <<
218 fVerboseLevel = theModel->fVerboseLevel; <<
219 } <<
220 return; <<
221 } 217 }
222 218
223 //....oooOO0OOooo........oooOO0OOooo........oo 219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
224 220
225 G4double G4PenelopeBremsstrahlungModel::CrossS << 221 G4double G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
226 con << 222 G4double kinEnergy,
227 G4d << 223 G4double Z,
228 G4d << 224 G4double,
229 G4d << 225 G4double cutEnergy,
>> 226 G4double)
230 { 227 {
>> 228 // Penelope model to calculate cross section for hard bremsstrahlung emission
>> 229 // (gamma energy > cutEnergy).
>> 230 //
>> 231 // The total bremsstrahlung cross section is read from database, following data
>> 232 // reported in the EEDL library
>> 233 // D.E.Cullen et al., Report UCRL-50400 (Lawrence Livermore National Laboratory) (1989)
>> 234 // The probability to have photon emission above a given threshold is calculated
>> 235 // analytically using the differential cross section model dSigma/dW = F(x)/x, where
>> 236 // W is the outgoing photon energy and x = W/E is the ratio of the photon energy to the
>> 237 // incident energy. The function F(x) is tabulated (for all elements) using 32 points in x
>> 238 // ranging from 1e-12 to 1. Data are derived from
>> 239 // S.M.Seltzer and M.J.Berger, At.Data Nucl.Data Tables 35,345 (1986)
>> 240 // Differential cross sections for electrons and positrons dSigma/dW are assumed to scale
>> 241 // with a function S(Z,E) which does not depend on W; therefore, only overall cross section
>> 242 // changes but not the shape of the photon energy spectrum.
231 // 243 //
232 if (fVerboseLevel > 3) <<
233 G4cout << "Calling CrossSectionPerVolume() <<
234 <<
235 SetupForMaterial(theParticle, material, ener <<
236 G4double crossPerMolecule = 0.; <<
237 244
238 G4PenelopeCrossSection* theXS = GetCrossSect << 245 if (verboseLevel > 3)
239 << 246 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeBremsstrahlungModel" << G4endl;
240 if (theXS) << 247
241 crossPerMolecule = theXS->GetHardCrossSect << 248 G4int iZ = (G4int) Z;
242 <<
243 G4double atomDensity = material->GetTotNbOfA <<
244 G4double atPerMol = fOscManager->GetAtomsPe <<
245 <<
246 if (fVerboseLevel > 3) <<
247 G4cout << "Material " << material->GetName <<
248 "atoms per molecule" << G4endl; <<
249 <<
250 G4double moleculeDensity = 0.; <<
251 if (atPerMol) <<
252 moleculeDensity = atomDensity/atPerMol; <<
253 <<
254 G4double crossPerVolume = crossPerMolecule*m <<
255 <<
256 if (fVerboseLevel > 2) <<
257 { <<
258 G4cout << "G4PenelopeBremsstrahlungModel " <<
259 G4cout << "Mean free path for gamma emissi <<
260 energy/keV << " keV = " << (1./crossPerV <<
261 } <<
262 249
263 return crossPerVolume; << 250 // VI - not needed in run time
264 } << 251 // if (!crossSectionHandler)
>> 252 // {
>> 253 // G4cout << "G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom" << G4endl;
>> 254 // G4cout << "The cross section handler is not correctly initialized" << G4endl;
>> 255 // G4Exception();
>> 256 // }
>> 257 G4double totalCs = crossSectionHandler->FindValue(iZ,kinEnergy);
>> 258 G4double cs = totalCs * energySpectrum->Probability(iZ,cutEnergy,kinEnergy,kinEnergy);
265 259
266 //....oooOO0OOooo........oooOO0OOooo........oo << 260 if (verboseLevel > 2)
>> 261 {
>> 262 G4cout << "Bremsstrahlung cross section at " << kinEnergy/MeV << " MeV for Z=" << Z <<
>> 263 " and energy > " << cutEnergy/keV << " keV --> " << cs/barn << " barn" << G4endl;
>> 264 G4cout << "Total bremsstrahlung cross section at " << kinEnergy/MeV << " MeV for Z=" <<
>> 265 Z << " --> " << totalCs/barn << " barn" << G4endl;
>> 266 }
>> 267 return cs;
267 268
268 //This is a dummy method. Never inkoved by the <<
269 //a warning if one tries to get Cross Sections <<
270 //G4EmCalculator. <<
271 G4double G4PenelopeBremsstrahlungModel::Comput <<
272 G4double, <<
273 G4double, <<
274 G4double, <<
275 G4double, <<
276 G4double) <<
277 { <<
278 G4cout << "*** G4PenelopeBremsstrahlungModel <<
279 G4cout << "Penelope Bremsstrahlung model v20 <<
280 G4cout << "so the result is always zero. For <<
281 G4cout << "GetCrossSectionPerVolume() or Get <<
282 return 0; <<
283 } 269 }
284 270
285 //....oooOO0OOooo........oooOO0OOooo........oo 271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286 << 272 G4double
287 G4double G4PenelopeBremsstrahlungModel::Comput << 273 G4PenelopeBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* theMaterial,
288 const G4ParticleDefinition* << 274 const G4ParticleDefinition* theParticle,
289 G4double kineticEnergy, << 275 G4double kineticEnergy,
290 G4double cutEnergy) << 276 G4double cutEnergy)
291 { 277 {
292 if (fVerboseLevel > 3) << 278 // Penelope model to calculate the stopping power (in [Energy]/[Length]) for soft
293 G4cout << "Calling ComputeDEDX() of G4Pene << 279 // bremsstrahlung emission (gamma energy < cutEnergy).
294 << 280 //
295 G4PenelopeCrossSection* theXS = GetCrossSect << 281 // The actual calculation is performed by the helper class
296 << 282 // G4PenelopeBremsstrahlungContinuous and its method CalculateStopping(). Notice:
297 G4double sPowerPerMolecule = 0.0; << 283 // CalculateStopping() gives the stopping cross section, namely the first momentum of
298 if (theXS) << 284 // dSigma/dW, restricted to W < cut (W = gamma energy) This is dimensionally:
299 sPowerPerMolecule = theXS->GetSoftStopping << 285 // [Energy]*[Surface]
300 << 286 // The calculation is performed by interpolation (in E = incident energy and
301 G4double atomDensity = material->GetTotNbOfA << 287 // x=W/E) from the tabulated data derived from
302 G4double atPerMol = fOscManager->GetAtomsPe << 288 // M.J.Berger and S.M.Seltzer, Report NBSIR 82-2550 (National Bureau of Standards) (1982);
303 << 289 // for electrons.
304 G4double moleculeDensity = 0.; << 290 // For positrons, dSigma/dW are assumed to scale with a function S(Z,E) with respect to electrons.
305 if (atPerMol) << 291 // An analytical approximation for the scaling function S(Z,E) is given in
306 moleculeDensity = atomDensity/atPerMol; << 292 // L.Kim et al., Phys.Rev.A 33,3002 (1986)
>> 293 //
>> 294 if (!stoppingPowerData)
>> 295 stoppingPowerData = new std::map<std::pair<G4int,G4double>,
>> 296 G4PenelopeBremsstrahlungContinuous*>;
>> 297
>> 298 const G4ElementVector* theElementVector = theMaterial->GetElementVector();
>> 299 const G4double* theAtomicNumDensityVector = theMaterial->GetAtomicNumDensityVector();
307 300
308 G4double sPowerPerVolume = sPowerPerMolecule << 301 G4double sPower = 0.0;
309 302
310 if (fVerboseLevel > 2) << 303 //Loop on the elements of the material
>> 304 for (size_t iel=0;iel<theMaterial->GetNumberOfElements();iel++)
311 { 305 {
312 G4cout << "G4PenelopeBremsstrahlungModel << 306 G4int iZ = (G4int) ((*theElementVector)[iel]->GetZ());
313 G4cout << "Stopping power < " << cutEner << 307 G4PenelopeBremsstrahlungContinuous* theContinuousCalculator =
314 kineticEnergy/keV << " keV = " << << 308 GetStoppingPowerData(iZ,cutEnergy,theParticle);
315 sPowerPerVolume/(keV/mm) << " keV/mm" << 309 sPower += theContinuousCalculator->CalculateStopping(kineticEnergy)*
>> 310 theAtomicNumDensityVector[iel];
316 } 311 }
317 return sPowerPerVolume; << 312
>> 313 if (verboseLevel > 2)
>> 314 {
>> 315 G4cout << "Bremsstrahlung stopping power at " << kineticEnergy/MeV
>> 316 << " MeV for material " << theMaterial->GetName()
>> 317 << " and energy < " << cutEnergy/keV << " keV --> "
>> 318 << sPower/(keV/mm) << " keV/mm" << G4endl;
>> 319 }
>> 320
>> 321 return sPower;
318 } 322 }
319 323
320 //....oooOO0OOooo........oooOO0OOooo........oo 324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
321 325
322 void G4PenelopeBremsstrahlungModel::SampleSeco << 326 void G4PenelopeBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
323 const G4MaterialCutsCouple* 327 const G4MaterialCutsCouple* couple,
324 const G4DynamicParticle* aDy 328 const G4DynamicParticle* aDynamicParticle,
325 G4double cutG, << 329 G4double cutG,G4double)
326 G4double) <<
327 { 330 {
328 if (fVerboseLevel > 3) << 331 // Penelope model to sample the final state for hard bremsstrahlung emission
>> 332 // (gamma energy < cutEnergy).
>> 333 // The energy distributionof the emitted photons is sampled according to the F(x)
>> 334 // function tabulated in the database from
>> 335 // S.M.Seltzer and M.J.Berger, At.Data Nucl.Data Tables 35,345 (1986)
>> 336 // The database contains the function F(x) (32 points) for 57 energies of the
>> 337 // incident electron between 1 keV and 100 GeV. For other primary energies,
>> 338 // logarithmic interpolation is used to obtain the values of the function F(x).
>> 339 // The double differential cross section dSigma/(dW dOmega), with W=photon energy,
>> 340 // is described as
>> 341 // dSigma/(dW dOmega) = dSigma/dW * p(Z,E,x,cosTheta)
>> 342 // where the shape function p depends on atomic number, incident energy and x=W/E.
>> 343 // Numerical values of the shape function, calculated by partial-waves methods, have been
>> 344 // reported in
>> 345 // L.Kissel et al., At.Data Nucl.Data.Tab. 28,381 (1983);
>> 346 // for Z=2,8,13,47,79 and 92; E=1,5,10,50,100 and 500 keV; x=0,0.6,0.8 and 0.95. The
>> 347 // function p(Z,E,x,cosTheta) is approximated by a Lorentz-dipole function as reported in
>> 348 // Penelope - A Code System for Monte Carlo Simulation of Electron and Photon Transport,
>> 349 // Workshop Proceedings Issy-les-Moulineaux, France, 5-7 November 2001, AEN-NEA;
>> 350 // The analytical function contains two adjustable parameters that are obtained by fitting
>> 351 // the complete solution from
>> 352 // L.Kissel et al., At.Data Nucl.Data.Tab. 28,381 (1983);
>> 353 // This allows the evaluation of p(Z,E,x,cosTheta) for any choice of Z, E and x. The actual
>> 354 // sampling of cos(theta) is performed in the helper class
>> 355 // G4PenelopeBremsstrahlungAngular, method ExtractCosTheta()
>> 356 // Energy and direction of the primary particle are updated according to energy-momentum
>> 357 // conservation. For positrons, it is sampled the same final state as for electrons.
>> 358 //
>> 359 if (verboseLevel > 3)
329 G4cout << "Calling SampleSecondaries() of 360 G4cout << "Calling SampleSecondaries() of G4PenelopeBremsstrahlungModel" << G4endl;
330 361
331 G4double kineticEnergy = aDynamicParticle->G 362 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
332 const G4Material* material = couple->GetMate << 363 const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
333 364
334 if (kineticEnergy <= fIntrinsicLowEnergyLimi 365 if (kineticEnergy <= fIntrinsicLowEnergyLimit)
335 { 366 {
336 fParticleChange->SetProposedKineticEnergy 367 fParticleChange->SetProposedKineticEnergy(0.);
337 fParticleChange->ProposeLocalEnergyDeposi 368 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
338 return ; 369 return ;
339 } 370 }
340 << 371
341 G4ParticleMomentum particleDirection0 = aDyn 372 G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
342 //This is the momentum 373 //This is the momentum
343 G4ThreeVector initialMomentum = aDynamicPar 374 G4ThreeVector initialMomentum = aDynamicParticle->GetMomentum();
>> 375
>> 376 //One can use Vladimir's selector!
>> 377 if (verboseLevel > 2)
>> 378 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
>> 379 // atom can be selected effitiantly if element selectors are initialised
>> 380 const G4Element* anElement = SelectRandomAtom(couple,theParticle,kineticEnergy);
>> 381 G4int iZ = (G4int) anElement->GetZ();
>> 382 if (verboseLevel > 2)
>> 383 G4cout << "Selected " << anElement->GetName() << G4endl;
>> 384 //
344 385
345 //Not enough energy to produce a secondary! << 386 //Sample gamma's energy according to the spectrum
346 if (kineticEnergy < cutG) << 387 G4double gammaEnergy = energySpectrum->SampleEnergy(iZ,cutG,kineticEnergy,kineticEnergy);
347 return; << 388
348 << 389 //Now sample cosTheta for the Gamma
349 if (fVerboseLevel > 3) << 390 G4double cosThetaPrimary = GetAngularDataForZ(iZ)->ExtractCosTheta(kineticEnergy,gammaEnergy);
350 G4cout << "Going to sample gamma energy fo << 391
351 "energy = " << kineticEnergy/keV << ", c <<
352 <<
353 //Sample gamma's energy according to the sp <<
354 G4double gammaEnergy = <<
355 fPenelopeFSHelper->SampleGammaEnergy(kinet <<
356 <<
357 if (fVerboseLevel > 3) <<
358 G4cout << "Sampled gamma energy: " << gamm <<
359 <<
360 //Now sample the direction for the Gamma. No <<
361 //Z is unused here, I plug 0. The informatio <<
362 G4ThreeVector gammaDirection1 = <<
363 fPenelopeAngular->SampleDirection(aDynami <<
364 <<
365 if (fVerboseLevel > 3) <<
366 G4cout << "Sampled cosTheta for e-: " << g <<
367 <<
368 G4double residualPrimaryEnergy = kineticEner 392 G4double residualPrimaryEnergy = kineticEnergy-gammaEnergy;
369 if (residualPrimaryEnergy < 0) 393 if (residualPrimaryEnergy < 0)
370 { 394 {
371 //Ok we have a problem, all energy goes 395 //Ok we have a problem, all energy goes with the gamma
372 gammaEnergy += residualPrimaryEnergy; 396 gammaEnergy += residualPrimaryEnergy;
373 residualPrimaryEnergy = 0.0; 397 residualPrimaryEnergy = 0.0;
374 } 398 }
375 399
>> 400 //Get primary kinematics
>> 401 G4double sinTheta = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
>> 402 G4double phi = twopi * G4UniformRand();
>> 403 G4ThreeVector gammaDirection1(sinTheta* std::cos(phi),
>> 404 sinTheta* std::sin(phi),
>> 405 cosThetaPrimary);
>> 406
>> 407 gammaDirection1.rotateUz(particleDirection0);
>> 408
376 //Produce final state according to momentum 409 //Produce final state according to momentum conservation
377 G4ThreeVector particleDirection1 = initialMo 410 G4ThreeVector particleDirection1 = initialMomentum - gammaEnergy*gammaDirection1;
378 particleDirection1 = particleDirection1.unit << 411 particleDirection1 = particleDirection1.unit(); //normalize
379 412
380 //Update the primary particle 413 //Update the primary particle
381 if (residualPrimaryEnergy > 0.) 414 if (residualPrimaryEnergy > 0.)
382 { 415 {
383 fParticleChange->ProposeMomentumDirectio 416 fParticleChange->ProposeMomentumDirection(particleDirection1);
384 fParticleChange->SetProposedKineticEnerg 417 fParticleChange->SetProposedKineticEnergy(residualPrimaryEnergy);
385 } 418 }
386 else 419 else
387 { 420 {
388 fParticleChange->SetProposedKineticEnerg 421 fParticleChange->SetProposedKineticEnergy(0.);
389 } 422 }
390 423
391 //Now produce the photon 424 //Now produce the photon
392 G4DynamicParticle* theGamma = new G4DynamicP 425 G4DynamicParticle* theGamma = new G4DynamicParticle(G4Gamma::Gamma(),
393 << 426 gammaDirection1,
394 << 427 gammaEnergy);
395 fvect->push_back(theGamma); 428 fvect->push_back(theGamma);
396 429
397 if (fVerboseLevel > 1) << 430 if (verboseLevel > 1)
398 { 431 {
399 G4cout << "----------------------------- 432 G4cout << "-----------------------------------------------------------" << G4endl;
400 G4cout << "Energy balance from G4Penelop 433 G4cout << "Energy balance from G4PenelopeBremsstrahlung" << G4endl;
401 G4cout << "Incoming primary energy: " << 434 G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
402 G4cout << "----------------------------- 435 G4cout << "-----------------------------------------------------------" << G4endl;
403 G4cout << "Outgoing primary energy: " << 436 G4cout << "Outgoing primary energy: " << residualPrimaryEnergy/keV << " keV" << G4endl;
404 G4cout << "Bremsstrahlung photon " << ga 437 G4cout << "Bremsstrahlung photon " << gammaEnergy/keV << " keV" << G4endl;
405 G4cout << "Total final state: " << (resi << 438 G4cout << "Total final state: " << (residualPrimaryEnergy+gammaEnergy)/keV
406 << " keV" << G4endl; << 439 << " keV" << G4endl;
407 G4cout << "----------------------------- 440 G4cout << "-----------------------------------------------------------" << G4endl;
408 } 441 }
409 << 442 if (verboseLevel > 0)
410 if (fVerboseLevel > 0) <<
411 { 443 {
412 G4double energyDiff = std::fabs(residual 444 G4double energyDiff = std::fabs(residualPrimaryEnergy+gammaEnergy-kineticEnergy);
413 if (energyDiff > 0.05*keV) 445 if (energyDiff > 0.05*keV)
414 G4cout << "Warning from G4PenelopeBrem << 446 G4cout << "Warning from G4PenelopeBremsstrahlung: problem with energy conservation: " <<
415 << <<
416 (residualPrimaryEnergy+gammaEnergy)/ 447 (residualPrimaryEnergy+gammaEnergy)/keV <<
417 " keV (final) vs. " << 448 " keV (final) vs. " <<
418 kineticEnergy/keV << " keV (initial) 449 kineticEnergy/keV << " keV (initial)" << G4endl;
419 } 450 }
420 return; <<
421 } 451 }
422 452
423 //....oooOO0OOooo........oooOO0OOooo........oo 453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
424 454
425 void G4PenelopeBremsstrahlungModel::ClearTable << 455 G4PenelopeBremsstrahlungAngular* G4PenelopeBremsstrahlungModel::GetAngularDataForZ(G4int iZ)
426 { << 456 {
427 if (!IsMaster() && !fLocalTable) << 457 if (!angularData)
428 //Should not be here! << 458 angularData = new std::map<G4int,G4PenelopeBremsstrahlungAngular*>;
429 G4Exception("G4PenelopeBremsstrahlungModel << 459
430 "em0100",FatalException,"Worker thread in << 460 if (angularData->count(iZ)) //the material already exists
431 << 461 return angularData->find(iZ)->second;
432 if (fXSTableElectron) << 462
433 { << 463 //Otherwise create a new object, store it and return it
434 for (auto& item : (*fXSTableElectron)) << 464 G4PenelopeBremsstrahlungAngular* theAngular = new G4PenelopeBremsstrahlungAngular(iZ);
435 delete item.second; << 465 angularData->insert(std::make_pair(iZ,theAngular));
436 delete fXSTableElectron; <<
437 fXSTableElectron = nullptr; <<
438 } <<
439 if (fXSTablePositron) <<
440 { <<
441 for (auto& item : (*fXSTablePositron)) <<
442 delete item.second; <<
443 delete fXSTablePositron; <<
444 fXSTablePositron = nullptr; <<
445 } <<
446 /* <<
447 if (fEnergyGrid) <<
448 delete fEnergyGrid; <<
449 */ <<
450 if (fPenelopeFSHelper) <<
451 fPenelopeFSHelper->ClearTables(IsMaster()) <<
452 <<
453 if (fVerboseLevel > 2) <<
454 G4cout << "G4PenelopeBremsstrahlungModel: <<
455 <<
456 return; <<
457 } <<
458 <<
459 //....oooOO0OOooo........oooOO0OOooo........oo <<
460 466
461 G4double G4PenelopeBremsstrahlungModel::MinEne << 467 if (angularData->count(iZ)) //the material should exist now
462 const G4MaterialCutsCouple* << 468 return angularData->find(iZ)->second;
463 { << 469 else
464 return fIntrinsicLowEnergyLimit; << 470 {
>> 471 G4Exception("Problem in G4PenelopeBremsstrahlungModel::GetAngularDataForZ()");
>> 472 return 0;
>> 473 }
465 } 474 }
466 475
467 //....oooOO0OOooo........oooOO0OOooo........oo 476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
468 477
469 void G4PenelopeBremsstrahlungModel::BuildXSTab << 478 G4PenelopeBremsstrahlungContinuous*
470 { << 479 G4PenelopeBremsstrahlungModel::GetStoppingPowerData(G4int iZ,G4double energyCut,
471 if (!IsMaster() && !fLocalTable) << 480 const G4ParticleDefinition*
472 //Should not be here! << 481 theParticle)
473 G4Exception("G4PenelopeBremsstrahlungModel << 482 {
474 "em0100",FatalException,"Worker thread in << 483 if (!stoppingPowerData)
475 << 484 stoppingPowerData = new std::map<std::pair<G4int,G4double>,G4PenelopeBremsstrahlungContinuous*>;
476 //The key of the map <<
477 std::pair<const G4Material*,G4double> theKey <<
478 <<
479 //The table already exists <<
480 if (fXSTableElectron->count(theKey) && fXSTa <<
481 return; <<
482 <<
483 // <<
484 //This method fills the G4PenelopeCrossSecti <<
485 //and for the given material/cut couple. <<
486 //Equivalent of subroutines EBRaT and PINaT <<
487 // <<
488 if (fVerboseLevel > 2) <<
489 { <<
490 G4cout << "G4PenelopeBremsstrahlungModel <<
491 G4cout << "for e+/e- in " << mat->GetNam <<
492 cut/keV << " keV " << G4endl; <<
493 } <<
494 <<
495 //Tables have been already created (checked <<
496 if (fEnergyGrid->GetVectorLength() != nBins) <<
497 { <<
498 G4ExceptionDescription ed; <<
499 ed << "Energy Grid looks not initialized <<
500 ed << nBins << " " << fEnergyGrid->GetVe <<
501 G4Exception("G4PenelopeBremsstrahlungMod <<
502 "em2016",FatalException,ed); <<
503 } <<
504 <<
505 G4PenelopeCrossSection* XSEntry = new G4Pene <<
506 G4PenelopeCrossSection* XSEntry2 = new G4Pen <<
507 const G4PhysicsTable* table = fPenelopeFSHel <<
508 <<
509 //loop on the energy grid <<
510 for (std::size_t bin=0;bin<nBins;++bin) <<
511 { <<
512 G4double energy = fEnergyGrid->GetLowEd <<
513 G4double XH0=0, XH1=0, XH2=0; <<
514 G4double XS0=0, XS1=0, XS2=0; <<
515 <<
516 //Global xs factor <<
517 G4double fact = fPenelopeFSHelper->GetE <<
518 ((energy+electron_mass_c2)*(energy+electron <<
519 (energy*(energy+2.0*electron_mass_c2))); <<
520 <<
521 G4double restrictedCut = cut/energy; <<
522 <<
523 //Now I need the dSigma/dX profile - in <<
524 //the 32-point x grid. Interpolation is <<
525 std::size_t nBinsX = fPenelopeFSHelper- <<
526 G4double* tempData = new G4double[nBins <<
527 G4double logene = G4Log(energy); <<
528 for (std::size_t ix=0;ix<nBinsX;++ix) <<
529 { <<
530 //find dSigma/dx for the given E. X belon <<
531 G4double val = (*table)[ix]->Value(logene <<
532 tempData[ix] = G4Exp(val); //back to the <<
533 } <<
534 <<
535 G4double XH0A = 0.; <<
536 if (restrictedCut <= 1) //calculate onl <<
537 XH0A = fPenelopeFSHelper->GetMomentumIntegr <<
538 fPenelopeFSHelper->GetMomentumIntegral(te <<
539 G4double XS1A = fPenelopeFSHelper->GetM <<
540 restrictedCut,0); <<
541 G4double XS2A = fPenelopeFSHelper->GetM <<
542 restrictedCut,1); <<
543 G4double XH1A=0, XH2A=0; <<
544 if (restrictedCut <=1) <<
545 { <<
546 XH1A = fPenelopeFSHelper->GetMomentumInte <<
547 XS1A; <<
548 XH2A = fPenelopeFSHelper->GetMomentumInte <<
549 XS2A; <<
550 } <<
551 delete[] tempData; <<
552 <<
553 XH0 = XH0A*fact; <<
554 XS1 = XS1A*fact*energy; <<
555 XH1 = XH1A*fact*energy; <<
556 XS2 = XS2A*fact*energy*energy; <<
557 XH2 = XH2A*fact*energy*energy; <<
558 <<
559 XSEntry->AddCrossSectionPoint(bin,energ <<
560 <<
561 //take care of positrons <<
562 G4double posCorrection = GetPositronXSC <<
563 XSEntry2->AddCrossSectionPoint(bin,ener <<
564 XH1*posCorrection, <<
565 XH2*posCorrection, <<
566 XS0, <<
567 XS1*posCorrection, <<
568 XS2*posCorrection); <<
569 } <<
570 <<
571 //Insert in the appropriate table <<
572 fXSTableElectron->insert(std::make_pair(theK <<
573 fXSTablePositron->insert(std::make_pair(theK <<
574 485
575 return; << 486 std::pair<G4int,G4double> theKey = std::make_pair(iZ,energyCut);
576 } <<
577 487
>> 488 if (stoppingPowerData->count(theKey)) //the material already exists
>> 489 return stoppingPowerData->find(theKey)->second;
578 490
579 //....oooOO0OOooo........oooOO0OOooo........oo << 491 //Otherwise create a new object, store it and return it
>> 492 G4String theParticleName = theParticle->GetParticleName();
>> 493 G4PenelopeBremsstrahlungContinuous* theContinuous = new
>> 494 G4PenelopeBremsstrahlungContinuous(iZ,energyCut,LowEnergyLimit(),HighEnergyLimit(),theParticleName);
>> 495 stoppingPowerData->insert(std::make_pair(theKey,theContinuous));
580 496
581 G4PenelopeCrossSection* << 497 if (stoppingPowerData->count(theKey)) //the material should exist now
582 G4PenelopeBremsstrahlungModel::GetCrossSection << 498 return stoppingPowerData->find(theKey)->second;
583 const G4Material* mat, << 499 else
584 G4double cut) <<
585 { <<
586 if (part != G4Electron::Electron() && part ! <<
587 { 500 {
588 G4ExceptionDescription ed; << 501 G4Exception("Problem in G4PenelopeBremsstrahlungModel::GetStoppingPowerData()");
589 ed << "Invalid particle: " << part->GetP << 502 return 0;
590 G4Exception("G4PenelopeBremsstrahlungMod <<
591 "em0001",FatalException,ed); <<
592 return nullptr; <<
593 } <<
594 <<
595 if (part == G4Electron::Electron()) <<
596 { <<
597 //Either Initialize() was not called, or <<
598 //not invoked <<
599 if (!fXSTableElectron) <<
600 { <<
601 //create a **thread-local** version of the <<
602 //Unit Tests <<
603 G4String excep = "The Cross Section <<
604 G4Exception("G4PenelopeBremsstrahlun <<
605 "em2013",JustWarning,excep); <<
606 fLocalTable = true; <<
607 fXSTableElectron = new <<
608 std::map< std::pair<const G4Material*,G4 <<
609 if (!fEnergyGrid) <<
610 fEnergyGrid = new G4PhysicsLogVector(Low <<
611 HighEnergyLimit(), <<
612 nBins-1); //one hidden bin is adde <<
613 if (!fPenelopeFSHelper) <<
614 fPenelopeFSHelper = new G4PenelopeBremss <<
615 } <<
616 std::pair<const G4Material*,G4double> th <<
617 if (fXSTableElectron->count(theKey)) //t <<
618 return fXSTableElectron->find(theKey)- <<
619 else <<
620 { <<
621 //If we are here, it means that Initialize <<
622 //not filled up. This can happen in a Unit <<
623 if (fVerboseLevel > 0) <<
624 { <<
625 //G4Exception (warning) is issued only <<
626 G4ExceptionDescription ed; <<
627 ed << "Unable to find e- table for " < <<
628 << cut/keV << " keV " << G4endl; <<
629 ed << "This can happen only in Unit Te <<
630 G4Exception("G4PenelopeBremsstrahlungM <<
631 "em2009",JustWarning,ed); <<
632 } <<
633 //protect file reading via autolock <<
634 G4AutoLock lock(&PenelopeBremsstrahlungMod <<
635 fPenelopeFSHelper->BuildScaledXSTabl <<
636 BuildXSTable(mat,cut); <<
637 lock.unlock(); <<
638 //now it should be ok <<
639 return fXSTableElectron->find(theKey)->sec <<
640 } <<
641 } <<
642 if (part == G4Positron::Positron()) <<
643 { <<
644 //Either Initialize() was not called, or <<
645 //not invoked <<
646 if (!fXSTablePositron) <<
647 { <<
648 G4String excep = "The Cross Section Table <<
649 G4Exception("G4PenelopeBremsstrahlun <<
650 "em2013",JustWarning,excep); <<
651 fLocalTable = true; <<
652 fXSTablePositron = new <<
653 std::map< std::pair<const G4Material*,G4 <<
654 if (!fEnergyGrid) <<
655 fEnergyGrid = new G4PhysicsLogVector(Low <<
656 HighEnergyLimit(), <<
657 nBins-1); //one hidden bin is adde <<
658 if (!fPenelopeFSHelper) <<
659 fPenelopeFSHelper = new G4Penelope <<
660 } <<
661 std::pair<const G4Material*,G4double> th <<
662 if (fXSTablePositron->count(theKey)) //t <<
663 return fXSTablePositron->find(theKey)- <<
664 else <<
665 { <<
666 //If we are here, it means that Initialize <<
667 //not filled up. This can happen in a Unit <<
668 if (fVerboseLevel > 0) <<
669 { <<
670 //Issue a G4Exception (warning) only i <<
671 G4ExceptionDescription ed; <<
672 ed << "Unable to find e+ table for " < <<
673 << cut/keV << " keV " << G4endl; <<
674 ed << "This can happen only in Unit Te <<
675 G4Exception("G4PenelopeBremsstrahlungM <<
676 "em2009",JustWarning,ed); <<
677 } <<
678 //protect file reading via autolock <<
679 G4AutoLock lock(&PenelopeBremsstrahlungMod <<
680 fPenelopeFSHelper->BuildScaledXSTabl <<
681 BuildXSTable(mat,cut); <<
682 lock.unlock(); <<
683 //now it should be ok <<
684 return fXSTablePositron->find(theKey)->sec <<
685 } <<
686 } 503 }
687 return nullptr; <<
688 } <<
689 <<
690 //....oooOO0OOooo........oooOO0OOooo........oo <<
691 <<
692 G4double G4PenelopeBremsstrahlungModel::GetPos <<
693 G4double energy) <<
694 { <<
695 //The electron-to-positron correction factor <<
696 //radiative stopping powers for positrons an <<
697 //by Kim et al. (1986) (cf. Berger and Seltz <<
698 //analytical approximation which reproduces <<
699 //accuracy <<
700 G4double t=G4Log(1.0+1e6*energy/ <<
701 (electron_mass_c2*fPenelopeFSHelper- <<
702 G4double corr = 1.0-G4Exp(-t*(1.2359e-1-t*(6 <<
703 (3.1516e-2-t*(7.7446e-3-t*(1.0595 <<
704 (7.0568e-5-t* <<
705 1.8080e-6))))))); <<
706 return corr; <<
707 } <<
708 <<
709 //....oooOO0OOooo........oooOO0OOooo........oo <<
710 <<
711 void G4PenelopeBremsstrahlungModel::SetParticl <<
712 { <<
713 if(!fParticle) { <<
714 fParticle = p; <<
715 } <<
716 } 504 }
717 505