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