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