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