| 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: G4PenelopeGammaConversionModel.cc 79186 2014-02-20 09:20:02Z gcosmo $
26 // 27 //
27 // Author: Luciano Pandola 28 // Author: Luciano Pandola
28 // 29 //
29 // History: 30 // History:
30 // -------- 31 // --------
31 // 13 Jan 2010 L Pandola First implementa 32 // 13 Jan 2010 L Pandola First implementation (updated to Penelope08)
32 // 24 May 2011 L Pandola Renamed (make v2 33 // 24 May 2011 L Pandola Renamed (make v2008 as default Penelope)
33 // 18 Sep 2013 L Pandola Migration to MT << 34 // 18 Sep 2013 L Pandola Migration to MT paradigm. Only master model deals with
34 // data and create 35 // data and creates shared tables
35 // 36 //
36 37
37 #include "G4PenelopeGammaConversionModel.hh" 38 #include "G4PenelopeGammaConversionModel.hh"
38 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
40 #include "G4ParticleDefinition.hh" 41 #include "G4ParticleDefinition.hh"
41 #include "G4MaterialCutsCouple.hh" 42 #include "G4MaterialCutsCouple.hh"
42 #include "G4ProductionCutsTable.hh" 43 #include "G4ProductionCutsTable.hh"
43 #include "G4DynamicParticle.hh" 44 #include "G4DynamicParticle.hh"
44 #include "G4Element.hh" 45 #include "G4Element.hh"
45 #include "G4Gamma.hh" 46 #include "G4Gamma.hh"
46 #include "G4Electron.hh" 47 #include "G4Electron.hh"
47 #include "G4Positron.hh" 48 #include "G4Positron.hh"
48 #include "G4PhysicsFreeVector.hh" 49 #include "G4PhysicsFreeVector.hh"
49 #include "G4MaterialCutsCouple.hh" 50 #include "G4MaterialCutsCouple.hh"
50 #include "G4AutoLock.hh" 51 #include "G4AutoLock.hh"
51 #include "G4Exp.hh" <<
52 <<
53 //....oooOO0OOooo........oooOO0OOooo........oo <<
54 const G4int G4PenelopeGammaConversionModel::fM <<
55 G4PhysicsFreeVector* G4PenelopeGammaConversion <<
56 G4double G4PenelopeGammaConversionModel::fAtom <<
57 1.2281e+02,7.3167e+01,6.9 <<
58 6.4696e+01,6.1228e+01,5.7 <<
59 5.0787e+01,4.7851e+01,4.6 <<
60 4.4503e+01,4.3815e+01,4.3 <<
61 4.1586e+01,4.0953e+01,4.0 <<
62 3.9756e+01,3.9144e+01,3.8 <<
63 3.7174e+01,3.6663e+01,3.5 <<
64 3.4688e+01,3.4197e+01,3.3 <<
65 3.3068e+01,3.2740e+01,3.2 <<
66 3.1884e+01,3.1622e+01,3.1 <<
67 3.0950e+01,3.0758e+01,3.0 <<
68 3.0097e+01,2.9832e+01,2.9 <<
69 2.9247e+01,2.9085e+01,2.8 <<
70 2.8580e+01,2.8442e+01,2.8 <<
71 2.7973e+01,2.7819e+01,2.7 <<
72 2.7285e+01,2.7093e+01,2.6 <<
73 2.6516e+01,2.6304e+01,2.6 <<
74 2.5730e+01,2.5577e+01,2.5 <<
75 2.5100e+01,2.4941e+01,2.4 <<
76 2.4506e+01,2.4391e+01,2.4 <<
77 2.4039e+01,2.3922e+01,2.3 <<
78 2.3621e+01,2.3523e+01,2.3 <<
79 2.3238e+01,2.3139e+01,2.3 <<
80 2.2833e+01,2.2694e+01,2.2 <<
81 2.2446e+01,2.2358e+01,2.2 <<
82 52
83 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 54
85 G4PenelopeGammaConversionModel::G4PenelopeGamm 55 G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel(const G4ParticleDefinition* part,
86 const G4String& nam) 56 const G4String& nam)
87 :G4VEmModel(nam),fParticleChange(nullptr),fP << 57 :G4VEmModel(nam),fParticleChange(0),fParticle(0),
88 fEffectiveCharge(nullptr),fMaterialInvScree << 58 logAtomicCrossSection(0),
89 fScreeningFunction(nullptr),fIsInitialised( << 59 fEffectiveCharge(0),fMaterialInvScreeningRadius(0),
90 { << 60 fScreeningFunction(0),isInitialised(false),fLocalTable(false)
>> 61
>> 62 {
91 fIntrinsicLowEnergyLimit = 2.0*electron_mass 63 fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2;
92 fIntrinsicHighEnergyLimit = 100.0*GeV; 64 fIntrinsicHighEnergyLimit = 100.0*GeV;
93 fSmallEnergy = 1.1*MeV; 65 fSmallEnergy = 1.1*MeV;
>> 66
>> 67 InitializeScreeningRadii();
94 68
95 if (part) 69 if (part)
96 SetParticle(part); 70 SetParticle(part);
97 71
98 // SetLowEnergyLimit(fIntrinsicLowEnergyLim 72 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
99 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 73 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
100 // 74 //
101 fVerboseLevel= 0; << 75 verboseLevel= 0;
102 // Verbosity scale: 76 // Verbosity scale:
103 // 0 = nothing << 77 // 0 = nothing
104 // 1 = warning for energy non-conservation << 78 // 1 = warning for energy non-conservation
105 // 2 = details of energy budget 79 // 2 = details of energy budget
106 // 3 = calculation of cross sections, file o 80 // 3 = calculation of cross sections, file openings, sampling of atoms
107 // 4 = entering in methods 81 // 4 = entering in methods
108 } 82 }
109 83
110 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 85
112 G4PenelopeGammaConversionModel::~G4PenelopeGam 86 G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel()
113 { 87 {
114 //Delete shared tables, they exist only in t 88 //Delete shared tables, they exist only in the master model
115 if (IsMaster() || fLocalTable) 89 if (IsMaster() || fLocalTable)
116 { << 90 {
117 for(G4int i=0; i<=fMaxZ; ++i) << 91 if (logAtomicCrossSection)
118 { 92 {
119 if(fLogAtomicCrossSection[i]) { << 93 /*
120 delete fLogAtomicCrossSection[i]; << 94 std::map <G4int,G4PhysicsFreeVector*>::iterator i;
121 fLogAtomicCrossSection[i] = nullptr; << 95 for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++)
122 } << 96 if (i->second) delete i->second;
>> 97 */
>> 98 delete logAtomicCrossSection;
123 } 99 }
124 if (fEffectiveCharge) 100 if (fEffectiveCharge)
125 delete fEffectiveCharge; 101 delete fEffectiveCharge;
126 if (fMaterialInvScreeningRadius) 102 if (fMaterialInvScreeningRadius)
127 delete fMaterialInvScreeningRadius; 103 delete fMaterialInvScreeningRadius;
128 if (fScreeningFunction) 104 if (fScreeningFunction)
129 delete fScreeningFunction; 105 delete fScreeningFunction;
130 } 106 }
131 } 107 }
132 108
>> 109
133 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
134 111
135 void G4PenelopeGammaConversionModel::Initialis 112 void G4PenelopeGammaConversionModel::Initialise(const G4ParticleDefinition* part,
136 const G4DataVector&) 113 const G4DataVector&)
137 { 114 {
138 if (fVerboseLevel > 3) << 115 if (verboseLevel > 3)
139 G4cout << "Calling G4PenelopeGammaConvers 116 G4cout << "Calling G4PenelopeGammaConversionModel::Initialise()" << G4endl;
140 117
141 SetParticle(part); 118 SetParticle(part);
142 119
143 //Only the master model creates/fills/destro 120 //Only the master model creates/fills/destroys the tables
144 if (IsMaster() && part == fParticle) 121 if (IsMaster() && part == fParticle)
145 { 122 {
>> 123 // logAtomicCrossSection is created only once, since it is never cleared
>> 124 if (!logAtomicCrossSection)
>> 125 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
>> 126
146 //delete old material data... 127 //delete old material data...
147 if (fEffectiveCharge) 128 if (fEffectiveCharge)
148 { 129 {
149 delete fEffectiveCharge; 130 delete fEffectiveCharge;
150 fEffectiveCharge = nullptr; << 131 fEffectiveCharge = 0;
151 } 132 }
152 if (fMaterialInvScreeningRadius) 133 if (fMaterialInvScreeningRadius)
153 { 134 {
154 delete fMaterialInvScreeningRadius; 135 delete fMaterialInvScreeningRadius;
155 fMaterialInvScreeningRadius = nullptr; << 136 fMaterialInvScreeningRadius = 0;
156 } 137 }
157 if (fScreeningFunction) 138 if (fScreeningFunction)
158 { 139 {
159 delete fScreeningFunction; 140 delete fScreeningFunction;
160 fScreeningFunction = nullptr; << 141 fScreeningFunction = 0;
161 } << 142 }
162 //and create new ones 143 //and create new ones
163 fEffectiveCharge = new std::map<const G4 144 fEffectiveCharge = new std::map<const G4Material*,G4double>;
164 fMaterialInvScreeningRadius = new std::m 145 fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
165 fScreeningFunction = new std::map<const 146 fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
166 << 147
167 G4ProductionCutsTable* theCoupleTable = << 148 G4ProductionCutsTable* theCoupleTable =
168 G4ProductionCutsTable::GetProductionCutsTabl 149 G4ProductionCutsTable::GetProductionCutsTable();
169 150
170 for (G4int i=0;i<(G4int)theCoupleTable-> << 151 for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
171 { 152 {
172 const G4Material* material = << 153 const G4Material* material =
173 theCoupleTable->GetMaterialCutsCouple(i) 154 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
174 const G4ElementVector* theElementVector = 155 const G4ElementVector* theElementVector = material->GetElementVector();
175 << 156
176 for (std::size_t j=0;j<material->GetNumber << 157 for (size_t j=0;j<material->GetNumberOfElements();j++)
177 { 158 {
178 G4int iZ = theElementVector->at(j)->Ge << 159 G4int iZ = (G4int) theElementVector->at(j)->GetZ();
179 //read data files only in the master 160 //read data files only in the master
180 if (iZ <= fMaxZ && !fLogAtomicCrossSe << 161 if (!logAtomicCrossSection->count(iZ))
181 ReadDataFile(iZ); 162 ReadDataFile(iZ);
182 } 163 }
183 164
184 //check if material data are available 165 //check if material data are available
185 if (!fEffectiveCharge->count(material)) 166 if (!fEffectiveCharge->count(material))
186 InitializeScreeningFunctions(material); << 167 InitializeScreeningFunctions(material);
187 } 168 }
188 if (fVerboseLevel > 0) { << 169
>> 170
>> 171 if (verboseLevel > 0) {
189 G4cout << "Penelope Gamma Conversion model v 172 G4cout << "Penelope Gamma Conversion model v2008 is initialized " << G4endl
190 << "Energy range: " 173 << "Energy range: "
191 << LowEnergyLimit() / MeV << " MeV - 174 << LowEnergyLimit() / MeV << " MeV - "
192 << HighEnergyLimit() / GeV << " GeV" 175 << HighEnergyLimit() / GeV << " GeV"
193 << G4endl; 176 << G4endl;
194 } 177 }
>> 178
195 } 179 }
196 if(fIsInitialised) return; << 180
>> 181
>> 182 if(isInitialised) return;
197 fParticleChange = GetParticleChangeForGamma( 183 fParticleChange = GetParticleChangeForGamma();
198 fIsInitialised = true; << 184 isInitialised = true;
199 } 185 }
200 186
201 //....oooOO0OOooo........oooOO0OOooo........oo 187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
202 188
203 void G4PenelopeGammaConversionModel::Initialis 189 void G4PenelopeGammaConversionModel::InitialiseLocal(const G4ParticleDefinition* part,
204 G4VEmModel *masterModel) 190 G4VEmModel *masterModel)
205 { 191 {
206 if (fVerboseLevel > 3) << 192 if (verboseLevel > 3)
207 G4cout << "Calling G4PenelopeGammaConvers 193 G4cout << "Calling G4PenelopeGammaConversionModel::InitialiseLocal()" << G4endl;
>> 194
208 // 195 //
209 //Check that particle matches: one might hav << 196 //Check that particle matches: one might have multiple master models (e.g.
210 //for e+ and e-). 197 //for e+ and e-).
211 // 198 //
212 if (part == fParticle) 199 if (part == fParticle)
213 { 200 {
214 //Get the const table pointers from the 201 //Get the const table pointers from the master to the workers
215 const G4PenelopeGammaConversionModel* th << 202 const G4PenelopeGammaConversionModel* theModel =
216 static_cast<G4PenelopeGammaConversionModel*> 203 static_cast<G4PenelopeGammaConversionModel*> (masterModel);
217 << 204
218 //Copy pointers to the data tables 205 //Copy pointers to the data tables
219 fEffectiveCharge = theModel->fEffectiveC 206 fEffectiveCharge = theModel->fEffectiveCharge;
220 fMaterialInvScreeningRadius = theModel-> 207 fMaterialInvScreeningRadius = theModel->fMaterialInvScreeningRadius;
221 fScreeningFunction = theModel->fScreenin << 208 fScreeningFunction = theModel->fScreeningFunction;
222 for(G4int i=0; i<=fMaxZ; ++i) << 209 logAtomicCrossSection = theModel->logAtomicCrossSection;
223 fLogAtomicCrossSection[i] = theModel->fLogAt << 210
224 <<
225 //Same verbosity for all workers, as the 211 //Same verbosity for all workers, as the master
226 fVerboseLevel = theModel->fVerboseLevel; << 212 verboseLevel = theModel->verboseLevel;
227 } << 213 }
228 214
229 return; 215 return;
230 } 216 }
231 217
232 //....oooOO0OOooo........oooOO0OOooo........oo 218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233 namespace { G4Mutex PenelopeGammaConversionMo 219 namespace { G4Mutex PenelopeGammaConversionModelMutex = G4MUTEX_INITIALIZER; }
234 220
235 G4double G4PenelopeGammaConversionModel::Compu 221 G4double G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom(
236 const G4ParticleDefinition 222 const G4ParticleDefinition*,
237 G4double energy, 223 G4double energy,
238 G4double Z, G4double, 224 G4double Z, G4double,
239 G4double, G4double) 225 G4double, G4double)
240 { 226 {
241 // 227 //
242 // Penelope model v2008. 228 // Penelope model v2008.
243 // Cross section (including triplet producti << 229 // Cross section (including triplet production) read from database and managed
244 // through the G4CrossSectionHandler utility 230 // through the G4CrossSectionHandler utility. Cross section data are from
245 // M.J. Berger and J.H. Hubbel (XCOM), Repor 231 // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598
246 // 232 //
247 << 233
248 if (energy < fIntrinsicLowEnergyLimit) 234 if (energy < fIntrinsicLowEnergyLimit)
249 return 0; 235 return 0;
250 236
251 G4int iZ = G4int(Z); << 237 G4int iZ = (G4int) Z;
252 << 238
253 if (!fLogAtomicCrossSection[iZ]) << 239 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
>> 240 //not invoked
>> 241 if (!logAtomicCrossSection)
>> 242 {
>> 243 //create a **thread-local** version of the table. Used only for G4EmCalculator and
>> 244 //Unit Tests
>> 245 fLocalTable = true;
>> 246 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
>> 247 }
>> 248 //now it should be ok
>> 249 if (!logAtomicCrossSection->count(iZ))
254 { 250 {
255 //If we are here, it means that Initial << 251 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
256 //not filled up. This can happen in a U 252 //not filled up. This can happen in a UnitTest or via G4EmCalculator
257 if (fVerboseLevel > 0) << 253 if (verboseLevel > 0)
258 { 254 {
259 //Issue a G4Exception (warning) only in v 255 //Issue a G4Exception (warning) only in verbose mode
260 G4ExceptionDescription ed; 256 G4ExceptionDescription ed;
261 ed << "Unable to retrieve the cross secti 257 ed << "Unable to retrieve the cross section table for Z=" << iZ << G4endl;
262 ed << "This can happen only in Unit Tests << 258 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
263 G4Exception("G4PenelopeGammaConversionMod 259 G4Exception("G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom()",
264 "em2018",JustWarning,ed); 260 "em2018",JustWarning,ed);
265 } 261 }
266 //protect file reading via autolock 262 //protect file reading via autolock
267 G4AutoLock lock(&PenelopeGammaConversio 263 G4AutoLock lock(&PenelopeGammaConversionModelMutex);
268 ReadDataFile(iZ); << 264 ReadDataFile(iZ);
269 lock.unlock(); 265 lock.unlock();
270 fLocalTable = true; <<
271 } 266 }
>> 267
272 G4double cs = 0; 268 G4double cs = 0;
273 G4double logene = G4Log(energy); << 269 G4double logene = std::log(energy);
274 G4PhysicsFreeVector* theVec = fLogAtomicCros << 270
>> 271 G4PhysicsFreeVector* theVec = logAtomicCrossSection->find(iZ)->second;
>> 272
275 G4double logXS = theVec->Value(logene); 273 G4double logXS = theVec->Value(logene);
276 cs = G4Exp(logXS); << 274 cs = std::exp(logXS);
277 275
278 if (fVerboseLevel > 2) << 276 if (verboseLevel > 2)
279 G4cout << "Gamma conversion cross section << 277 G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z <<
280 " = " << cs/barn << " barn" << G4endl; 278 " = " << cs/barn << " barn" << G4endl;
281 return cs; 279 return cs;
282 } 280 }
283 281
284 //....oooOO0OOooo........oooOO0OOooo........oo 282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
285 283
286 void << 284 void
287 G4PenelopeGammaConversionModel::SampleSecondar 285 G4PenelopeGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
288 const G4MaterialCutsCouple* coup 286 const G4MaterialCutsCouple* couple,
289 const G4DynamicParticle* aDynami 287 const G4DynamicParticle* aDynamicGamma,
290 G4double, 288 G4double,
291 G4double) 289 G4double)
292 { 290 {
293 // 291 //
294 // Penelope model v2008. 292 // Penelope model v2008.
295 // Final state is sampled according to the B 293 // Final state is sampled according to the Bethe-Heitler model with Coulomb
296 // corrections, according to the semi-empiri 294 // corrections, according to the semi-empirical model of
297 // J. Baro' et al., Radiat. Phys. Chem. 44 295 // J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531.
298 // 296 //
299 // The model uses the high energy Coulomb co << 297 // The model uses the high energy Coulomb correction from
300 // H. Davies et al., Phys. Rev. 93 (1954) 7 298 // H. Davies et al., Phys. Rev. 93 (1954) 788
301 // and atomic screening radii tabulated from << 299 // and atomic screening radii tabulated from
302 // J.H. Hubbel et al., J. Phys. Chem. Ref. 300 // J.H. Hubbel et al., J. Phys. Chem. Ref. Data 9 (1980) 1023
303 // for Z= 1 to 92. << 301 // for Z= 1 to 92.
304 // 302 //
305 if (fVerboseLevel > 3) << 303 if (verboseLevel > 3)
306 G4cout << "Calling SamplingSecondaries() o 304 G4cout << "Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" << G4endl;
307 305
308 G4double photonEnergy = aDynamicGamma->GetKi 306 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
309 307
310 // Always kill primary 308 // Always kill primary
311 fParticleChange->ProposeTrackStatus(fStopAnd 309 fParticleChange->ProposeTrackStatus(fStopAndKill);
312 fParticleChange->SetProposedKineticEnergy(0. 310 fParticleChange->SetProposedKineticEnergy(0.);
313 311
314 if (photonEnergy <= fIntrinsicLowEnergyLimit 312 if (photonEnergy <= fIntrinsicLowEnergyLimit)
315 { 313 {
316 fParticleChange->ProposeLocalEnergyDepos 314 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
317 return ; 315 return ;
318 } 316 }
319 317
320 G4ParticleMomentum photonDirection = aDynami 318 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
321 const G4Material* mat = couple->GetMaterial( 319 const G4Material* mat = couple->GetMaterial();
322 320
323 //Either Initialize() was not called, or we << 321 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
324 //not invoked 322 //not invoked
325 if (!fEffectiveCharge) 323 if (!fEffectiveCharge)
326 { 324 {
327 //create a **thread-local** version of t << 325 //create a **thread-local** version of the table. Used only for G4EmCalculator and
328 //Unit Tests 326 //Unit Tests
329 fLocalTable = true; << 327 fLocalTable = true;
330 fEffectiveCharge = new std::map<const G4 328 fEffectiveCharge = new std::map<const G4Material*,G4double>;
331 fMaterialInvScreeningRadius = new std::m 329 fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
332 fScreeningFunction = new std::map<const 330 fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
333 } 331 }
334 332
335 if (!fEffectiveCharge->count(mat)) 333 if (!fEffectiveCharge->count(mat))
336 { 334 {
337 //If we are here, it means that Initiali << 335 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
338 //not filled up. This can happen in a Un 336 //not filled up. This can happen in a UnitTest or via G4EmCalculator
339 if (fVerboseLevel > 0) << 337 if (verboseLevel > 0)
340 { 338 {
341 //Issue a G4Exception (warning) only in ve 339 //Issue a G4Exception (warning) only in verbose mode
342 G4ExceptionDescription ed; 340 G4ExceptionDescription ed;
343 ed << "Unable to allocate the EffectiveCha << 341 ed << "Unable to allocate the EffectiveCharge data for " <<
344 mat->GetName() << G4endl; 342 mat->GetName() << G4endl;
345 ed << "This can happen only in Unit Tests" << 343 ed << "This can happen only in Unit Tests" << G4endl;
346 G4Exception("G4PenelopeGammaConversionMode 344 G4Exception("G4PenelopeGammaConversionModel::SampleSecondaries()",
347 "em2019",JustWarning,ed); 345 "em2019",JustWarning,ed);
348 } 346 }
349 //protect file reading via autolock 347 //protect file reading via autolock
350 G4AutoLock lock(&PenelopeGammaConversion 348 G4AutoLock lock(&PenelopeGammaConversionModelMutex);
351 InitializeScreeningFunctions(mat); 349 InitializeScreeningFunctions(mat);
352 lock.unlock(); 350 lock.unlock();
353 } 351 }
354 352
355 // eps is the fraction of the photon energy 353 // eps is the fraction of the photon energy assigned to e- (including rest mass)
356 G4double eps = 0; 354 G4double eps = 0;
357 G4double eki = electron_mass_c2/photonEnergy 355 G4double eki = electron_mass_c2/photonEnergy;
358 356
359 //Do it fast for photon energy < 1.1 MeV (cl 357 //Do it fast for photon energy < 1.1 MeV (close to threshold)
360 if (photonEnergy < fSmallEnergy) 358 if (photonEnergy < fSmallEnergy)
361 eps = eki + (1.0-2.0*eki)*G4UniformRand(); 359 eps = eki + (1.0-2.0*eki)*G4UniformRand();
362 else 360 else
363 { 361 {
364 //Complete calculation 362 //Complete calculation
365 G4double effC = fEffectiveCharge->find(m 363 G4double effC = fEffectiveCharge->find(mat)->second;
366 G4double alz = effC*fine_structure_const 364 G4double alz = effC*fine_structure_const;
367 G4double T = std::sqrt(2.0*eki); 365 G4double T = std::sqrt(2.0*eki);
368 G4double F00=(-1.774-1.210e1*alz+1.118e1 366 G4double F00=(-1.774-1.210e1*alz+1.118e1*alz*alz)*T
369 +(8.523+7.326e1*alz-4.441e1*alz*alz)* 367 +(8.523+7.326e1*alz-4.441e1*alz*alz)*T*T
370 -(1.352e1+1.211e2*alz-9.641e1*alz*alz 368 -(1.352e1+1.211e2*alz-9.641e1*alz*alz)*T*T*T
371 +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T 369 +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T;
372 << 370
373 G4double F0b = fScreeningFunction->find( 371 G4double F0b = fScreeningFunction->find(mat)->second.second;
374 G4double g0 = F0b + F00; 372 G4double g0 = F0b + F00;
375 G4double invRad = fMaterialInvScreeningR 373 G4double invRad = fMaterialInvScreeningRadius->find(mat)->second;
376 G4double bmin = 4.0*eki/invRad; 374 G4double bmin = 4.0*eki/invRad;
377 std::pair<G4double,G4double> scree = Ge 375 std::pair<G4double,G4double> scree = GetScreeningFunctions(bmin);
378 G4double g1 = scree.first; 376 G4double g1 = scree.first;
379 G4double g2 = scree.second; 377 G4double g2 = scree.second;
380 G4double g1min = g1+g0; 378 G4double g1min = g1+g0;
381 G4double g2min = g2+g0; 379 G4double g2min = g2+g0;
382 G4double xr = 0.5-eki; 380 G4double xr = 0.5-eki;
383 G4double a1 = 2.*g1min*xr*xr/3.; << 381 G4double a1 = 2.*g1min*xr*xr/3.;
384 G4double p1 = a1/(a1+g2min); 382 G4double p1 = a1/(a1+g2min);
385 383
386 G4bool loopAgain = false; 384 G4bool loopAgain = false;
387 //Random sampling of eps 385 //Random sampling of eps
388 do{ 386 do{
389 loopAgain = false; 387 loopAgain = false;
390 if (G4UniformRand() <= p1) 388 if (G4UniformRand() <= p1)
391 { 389 {
392 G4double ru2m1 = 2.0*G4UniformRand()-1. 390 G4double ru2m1 = 2.0*G4UniformRand()-1.0;
393 if (ru2m1 < 0) 391 if (ru2m1 < 0)
394 eps = 0.5-xr*std::pow(std::abs(ru2m1), 392 eps = 0.5-xr*std::pow(std::abs(ru2m1),1./3.);
395 else 393 else
396 eps = 0.5+xr*std::pow(ru2m1,1./3.); 394 eps = 0.5+xr*std::pow(ru2m1,1./3.);
397 G4double B = eki/(invRad*eps*(1.0-eps)); 395 G4double B = eki/(invRad*eps*(1.0-eps));
398 scree = GetScreeningFunctions(B); 396 scree = GetScreeningFunctions(B);
399 g1 = scree.first; 397 g1 = scree.first;
400 g1 = std::max(g1+g0,0.); 398 g1 = std::max(g1+g0,0.);
401 if (G4UniformRand()*g1min > g1) << 399 if (G4UniformRand()*g1min > g1)
402 loopAgain = true; 400 loopAgain = true;
403 } 401 }
404 else 402 else
405 { 403 {
406 eps = eki+2.0*xr*G4UniformRand(); 404 eps = eki+2.0*xr*G4UniformRand();
407 G4double B = eki/(invRad*eps*(1.0-eps)); 405 G4double B = eki/(invRad*eps*(1.0-eps));
408 scree = GetScreeningFunctions(B); 406 scree = GetScreeningFunctions(B);
409 g2 = scree.second; 407 g2 = scree.second;
410 g2 = std::max(g2+g0,0.); 408 g2 = std::max(g2+g0,0.);
411 if (G4UniformRand()*g2min > g2) 409 if (G4UniformRand()*g2min > g2)
412 loopAgain = true; 410 loopAgain = true;
413 } << 411 }
414 }while(loopAgain); 412 }while(loopAgain);
>> 413
415 } 414 }
416 if (fVerboseLevel > 4) << 415 if (verboseLevel > 4)
417 G4cout << "Sampled eps = " << eps << G4end 416 G4cout << "Sampled eps = " << eps << G4endl;
418 417
419 G4double electronTotEnergy = eps*photonEnerg 418 G4double electronTotEnergy = eps*photonEnergy;
420 G4double positronTotEnergy = (1.0-eps)*photo 419 G4double positronTotEnergy = (1.0-eps)*photonEnergy;
421 << 420
422 // Scattered electron (positron) angles. ( Z 421 // Scattered electron (positron) angles. ( Z - axis along the parent photon)
423 422
424 //electron kinematics 423 //electron kinematics
425 G4double electronKineEnergy = std::max(0.,el << 424 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
426 G4double costheta_el = G4UniformRand()*2.0-1 425 G4double costheta_el = G4UniformRand()*2.0-1.0;
427 G4double kk = std::sqrt(electronKineEnergy*( 426 G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
428 costheta_el = (costheta_el*electronTotEnergy 427 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
429 G4double phi_el = twopi * G4UniformRand() ; 428 G4double phi_el = twopi * G4UniformRand() ;
430 G4double dirX_el = std::sqrt(1.-costheta_el* 429 G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
431 G4double dirY_el = std::sqrt(1.-costheta_el* 430 G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
432 G4double dirZ_el = costheta_el; 431 G4double dirZ_el = costheta_el;
433 432
434 //positron kinematics 433 //positron kinematics
435 G4double positronKineEnergy = std::max(0.,po 434 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
436 G4double costheta_po = G4UniformRand()*2.0-1 435 G4double costheta_po = G4UniformRand()*2.0-1.0;
437 kk = std::sqrt(positronKineEnergy*(positronK 436 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
438 costheta_po = (costheta_po*positronTotEnergy 437 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
439 G4double phi_po = twopi * G4UniformRand() ; 438 G4double phi_po = twopi * G4UniformRand() ;
440 G4double dirX_po = std::sqrt(1.-costheta_po* 439 G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
441 G4double dirY_po = std::sqrt(1.-costheta_po* 440 G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
442 G4double dirZ_po = costheta_po; 441 G4double dirZ_po = costheta_po;
443 442
444 // Kinematics of the created pair: 443 // Kinematics of the created pair:
445 // the electron and positron are assumed to 444 // the electron and positron are assumed to have a symetric angular
446 // distribution with respect to the Z axis a 445 // distribution with respect to the Z axis along the parent photon
447 G4double localEnergyDeposit = 0. ; 446 G4double localEnergyDeposit = 0. ;
448 << 447
449 if (electronKineEnergy > 0.0) 448 if (electronKineEnergy > 0.0)
450 { 449 {
451 G4ThreeVector electronDirection ( dirX_e 450 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
452 electronDirection.rotateUz(photonDirecti 451 electronDirection.rotateUz(photonDirection);
453 G4DynamicParticle* electron = new G4Dyna 452 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
454 electronDirection, 453 electronDirection,
455 electronKineEnergy); 454 electronKineEnergy);
456 fvect->push_back(electron); 455 fvect->push_back(electron);
457 } 456 }
458 else 457 else
459 { 458 {
460 localEnergyDeposit += electronKineEnergy 459 localEnergyDeposit += electronKineEnergy;
461 electronKineEnergy = 0; 460 electronKineEnergy = 0;
462 } 461 }
463 462
464 //Generate the positron. Real particle in an 463 //Generate the positron. Real particle in any case, because it will annihilate. If below
465 //threshold, produce it at rest 464 //threshold, produce it at rest
466 if (positronKineEnergy < 0.0) 465 if (positronKineEnergy < 0.0)
467 { 466 {
468 localEnergyDeposit += positronKineEnergy 467 localEnergyDeposit += positronKineEnergy;
469 positronKineEnergy = 0; //produce it at 468 positronKineEnergy = 0; //produce it at rest
470 } 469 }
471 G4ThreeVector positronDirection(dirX_po,dirY 470 G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po);
472 positronDirection.rotateUz(photonDirection); 471 positronDirection.rotateUz(photonDirection);
473 G4DynamicParticle* positron = new G4DynamicP 472 G4DynamicParticle* positron = new G4DynamicParticle(G4Positron::Positron(),
474 positronDirection, positronK 473 positronDirection, positronKineEnergy);
475 fvect->push_back(positron); 474 fvect->push_back(positron);
476 475
477 //Add rest of energy to the local energy dep 476 //Add rest of energy to the local energy deposit
478 fParticleChange->ProposeLocalEnergyDeposit(l 477 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
479 << 478
480 if (fVerboseLevel > 1) << 479 if (verboseLevel > 1)
481 { 480 {
482 G4cout << "----------------------------- 481 G4cout << "-----------------------------------------------------------" << G4endl;
483 G4cout << "Energy balance from G4Penelop 482 G4cout << "Energy balance from G4PenelopeGammaConversion" << G4endl;
484 G4cout << "Incoming photon energy: " << 483 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
485 G4cout << "----------------------------- 484 G4cout << "-----------------------------------------------------------" << G4endl;
486 if (electronKineEnergy) 485 if (electronKineEnergy)
487 G4cout << "Electron (explicitly produced) " << 486 G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV"
488 << G4endl; 487 << G4endl;
489 if (positronKineEnergy) 488 if (positronKineEnergy)
490 G4cout << "Positron (not at rest) " << posit 489 G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl;
491 G4cout << "Rest masses of e+/- " << 2.0* 490 G4cout << "Rest masses of e+/- " << 2.0*electron_mass_c2/keV << " keV" << G4endl;
492 if (localEnergyDeposit) 491 if (localEnergyDeposit)
493 G4cout << "Local energy deposit " << localEn 492 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
494 G4cout << "Total final state: " << (elec 493 G4cout << "Total final state: " << (electronKineEnergy+positronKineEnergy+
495 localEnergyDeposit+2.0*electron_ma 494 localEnergyDeposit+2.0*electron_mass_c2)/keV <<
496 " keV" << G4endl; 495 " keV" << G4endl;
497 G4cout << "----------------------------- 496 G4cout << "-----------------------------------------------------------" << G4endl;
498 } 497 }
499 if (fVerboseLevel > 0) << 498 if (verboseLevel > 0)
500 { 499 {
501 G4double energyDiff = std::fabs(electron 500 G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
502 localEnergyDeposit+2.0*electron_ 501 localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
503 if (energyDiff > 0.05*keV) 502 if (energyDiff > 0.05*keV)
504 G4cout << "Warning from G4PenelopeGammaConve << 503 G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: "
505 << (electronKineEnergy+positronKineEn 504 << (electronKineEnergy+positronKineEnergy+
506 localEnergyDeposit+2.0*electron_mass_c2 << 505 localEnergyDeposit+2.0*electron_mass_c2)/keV
507 << " keV (final) vs. " << photonEnerg 506 << " keV (final) vs. " << photonEnergy/keV << " keV (initial)" << G4endl;
508 } << 507 }
509 } 508 }
510 509
511 //....oooOO0OOooo........oooOO0OOooo........oo 510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
512 511
513 void G4PenelopeGammaConversionModel::ReadDataF 512 void G4PenelopeGammaConversionModel::ReadDataFile(const G4int Z)
514 { 513 {
515 if (!IsMaster()) << 514 if (!IsMaster())
516 //Should not be here! 515 //Should not be here!
517 G4Exception("G4PenelopeGammaConversionMode 516 G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
518 "em0100",FatalException,"Worker thread in << 517 "em0100",FatalException,"Worker thread in this method");
519 518
520 if (fVerboseLevel > 2) << 519 if (verboseLevel > 2)
521 { 520 {
522 G4cout << "G4PenelopeGammaConversionMode 521 G4cout << "G4PenelopeGammaConversionModel::ReadDataFile()" << G4endl;
523 G4cout << "Going to read Gamma Conversio 522 G4cout << "Going to read Gamma Conversion data files for Z=" << Z << G4endl;
524 } 523 }
525 << 524
526 const char* path = G4FindDataDir("G4LEDATA << 525 char* path = getenv("G4LEDATA");
527 if(!path) << 526 if (!path)
528 { 527 {
529 G4String excep = << 528 G4String excep =
530 "G4PenelopeGammaConversionModel - G4LEDATA e 529 "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!";
531 G4Exception("G4PenelopeGammaConversionMo 530 G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
532 "em0006",FatalException,excep); 531 "em0006",FatalException,excep);
533 return; 532 return;
534 } 533 }
535 << 534
536 /* 535 /*
537 Read the cross section file 536 Read the cross section file
538 */ 537 */
539 std::ostringstream ost; 538 std::ostringstream ost;
540 if (Z>9) 539 if (Z>9)
541 ost << path << "/penelope/pairproduction/p 540 ost << path << "/penelope/pairproduction/pdgpp" << Z << ".p08";
542 else 541 else
543 ost << path << "/penelope/pairproduction/p 542 ost << path << "/penelope/pairproduction/pdgpp0" << Z << ".p08";
544 std::ifstream file(ost.str().c_str()); 543 std::ifstream file(ost.str().c_str());
545 if (!file.is_open()) 544 if (!file.is_open())
546 { 545 {
547 G4String excep = "G4PenelopeGammaConvers << 546 G4String excep = "G4PenelopeGammaConversionModel - data file " +
548 G4String(ost.str()) + " not found!"; 547 G4String(ost.str()) + " not found!";
549 G4Exception("G4PenelopeGammaConversionMo 548 G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
550 "em0003",FatalException,excep); 549 "em0003",FatalException,excep);
551 } 550 }
552 551
553 //I have to know in advance how many points 552 //I have to know in advance how many points are in the data list
554 //to initialize the G4PhysicsFreeVector() 553 //to initialize the G4PhysicsFreeVector()
555 std::size_t ndata=0; << 554 size_t ndata=0;
556 G4String line; 555 G4String line;
557 while( getline(file, line) ) 556 while( getline(file, line) )
558 ndata++; 557 ndata++;
559 ndata -= 1; //remove one header line 558 ndata -= 1; //remove one header line
>> 559 //G4cout << "Found: " << ndata << " lines" << G4endl;
560 560
561 file.clear(); 561 file.clear();
562 file.close(); 562 file.close();
563 file.open(ost.str().c_str()); 563 file.open(ost.str().c_str());
564 G4int readZ =0; 564 G4int readZ =0;
565 file >> readZ; << 565 file >> readZ;
566 566
567 if (fVerboseLevel > 3) << 567 if (verboseLevel > 3)
568 G4cout << "Element Z=" << Z << G4endl; 568 G4cout << "Element Z=" << Z << G4endl;
569 569
570 //check the right file is opened. 570 //check the right file is opened.
571 if (readZ != Z) 571 if (readZ != Z)
572 { 572 {
573 G4ExceptionDescription ed; 573 G4ExceptionDescription ed;
574 ed << "Corrupted data file for Z=" << Z << 574 ed << "Corrupted data file for Z=" << Z << G4endl;
575 G4Exception("G4PenelopeGammaConversionMo 575 G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
576 "em0005",FatalException,ed); 576 "em0005",FatalException,ed);
577 } 577 }
578 578
579 fLogAtomicCrossSection[Z] = new G4PhysicsFre << 579 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(ndata);
580 G4double ene=0,xs=0; 580 G4double ene=0,xs=0;
581 for (std::size_t i=0;i<ndata;++i) << 581 for (size_t i=0;i<ndata;i++)
582 { 582 {
583 file >> ene >> xs; 583 file >> ene >> xs;
584 //dimensional quantities 584 //dimensional quantities
585 ene *= eV; 585 ene *= eV;
586 xs *= barn; 586 xs *= barn;
587 if (xs < 1e-40*cm2) //protection against 587 if (xs < 1e-40*cm2) //protection against log(0)
588 xs = 1e-40*cm2; 588 xs = 1e-40*cm2;
589 fLogAtomicCrossSection[Z]->PutValue(i,G4 << 589 theVec->PutValue(i,std::log(ene),std::log(xs));
590 } 590 }
591 file.close(); 591 file.close();
592 592
>> 593 if (!logAtomicCrossSection)
>> 594 {
>> 595 G4ExceptionDescription ed;
>> 596 ed << "Problem with allocation of logAtomicCrossSection data table " << G4endl;
>> 597 G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
>> 598 "em2020",FatalException,ed);
>> 599 delete theVec;
>> 600 return;
>> 601 }
>> 602 logAtomicCrossSection->insert(std::make_pair(Z,theVec));
>> 603
593 return; 604 return;
>> 605
>> 606 }
>> 607
>> 608 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 609
>> 610 void G4PenelopeGammaConversionModel::InitializeScreeningRadii()
>> 611 {
>> 612 G4double temp[99] = {1.2281e+02,7.3167e+01,6.9228e+01,6.7301e+01,6.4696e+01,
>> 613 6.1228e+01,5.7524e+01,5.4033e+01,5.0787e+01,4.7851e+01,4.6373e+01,
>> 614 4.5401e+01,4.4503e+01,4.3815e+01,4.3074e+01,4.2321e+01,4.1586e+01,
>> 615 4.0953e+01,4.0524e+01,4.0256e+01,3.9756e+01,3.9144e+01,3.8462e+01,
>> 616 3.7778e+01,3.7174e+01,3.6663e+01,3.5986e+01,3.5317e+01,3.4688e+01,
>> 617 3.4197e+01,3.3786e+01,3.3422e+01,3.3068e+01,3.2740e+01,3.2438e+01,
>> 618 3.2143e+01,3.1884e+01,3.1622e+01,3.1438e+01,3.1142e+01,3.0950e+01,
>> 619 3.0758e+01,3.0561e+01,3.0285e+01,3.0097e+01,2.9832e+01,2.9581e+01,
>> 620 2.9411e+01,2.9247e+01,2.9085e+01,2.8930e+01,2.8721e+01,2.8580e+01,
>> 621 2.8442e+01,2.8312e+01,2.8139e+01,2.7973e+01,2.7819e+01,2.7675e+01,
>> 622 2.7496e+01,2.7285e+01,2.7093e+01,2.6911e+01,2.6705e+01,2.6516e+01,
>> 623 2.6304e+01,2.6108e+01,2.5929e+01,2.5730e+01,2.5577e+01,2.5403e+01,
>> 624 2.5245e+01,2.5100e+01,2.4941e+01,2.4790e+01,2.4655e+01,2.4506e+01,
>> 625 2.4391e+01,2.4262e+01,2.4145e+01,2.4039e+01,2.3922e+01,2.3813e+01,
>> 626 2.3712e+01,2.3621e+01,2.3523e+01,2.3430e+01,2.3331e+01,2.3238e+01,
>> 627 2.3139e+01,2.3048e+01,2.2967e+01,2.2833e+01,2.2694e+01,2.2624e+01,
>> 628 2.2545e+01,2.2446e+01,2.2358e+01,2.2264e+01};
>> 629
>> 630 //copy temporary vector in class data member
>> 631 for (G4int i=0;i<99;i++)
>> 632 fAtomicScreeningRadius[i] = temp[i];
594 } 633 }
595 634
596 //....oooOO0OOooo........oooOO0OOooo........oo 635 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
597 636
598 void G4PenelopeGammaConversionModel::Initializ 637 void G4PenelopeGammaConversionModel::InitializeScreeningFunctions(const G4Material* material)
599 { 638 {
>> 639 /*
>> 640 if (!IsMaster())
>> 641 //Should not be here!
>> 642 G4Exception("G4PenelopeGammaConversionModel::InitializeScreeningFunctions()",
>> 643 "em01001",FatalException,"Worker thread in this method");
>> 644 */
>> 645
600 // This is subroutine GPPa0 of Penelope 646 // This is subroutine GPPa0 of Penelope
601 // 647 //
602 // 1) calculate the effective Z for the purp 648 // 1) calculate the effective Z for the purpose
603 // 649 //
604 G4double zeff = 0; 650 G4double zeff = 0;
605 G4int intZ = 0; 651 G4int intZ = 0;
606 G4int nElements = (G4int)material->GetNumber << 652 G4int nElements = material->GetNumberOfElements();
607 const G4ElementVector* elementVector = mater 653 const G4ElementVector* elementVector = material->GetElementVector();
608 654
609 //avoid calculations if only one building el 655 //avoid calculations if only one building element!
610 if (nElements == 1) 656 if (nElements == 1)
611 { 657 {
612 zeff = (*elementVector)[0]->GetZ(); 658 zeff = (*elementVector)[0]->GetZ();
613 intZ = (G4int) zeff; 659 intZ = (G4int) zeff;
614 } 660 }
615 else // many elements...let's do the calcula 661 else // many elements...let's do the calculation
616 { 662 {
617 const G4double* fractionVector = materia 663 const G4double* fractionVector = material->GetVecNbOfAtomsPerVolume();
618 << 664
619 G4double atot = 0; 665 G4double atot = 0;
620 for (G4int i=0;i<nElements;i++) 666 for (G4int i=0;i<nElements;i++)
621 { 667 {
622 G4double Zelement = (*elementVector)[i]->G 668 G4double Zelement = (*elementVector)[i]->GetZ();
623 G4double Aelement = (*elementVector)[i]->G 669 G4double Aelement = (*elementVector)[i]->GetAtomicMassAmu();
624 atot += Aelement*fractionVector[i]; 670 atot += Aelement*fractionVector[i];
625 zeff += Zelement*Aelement*fractionVector[i 671 zeff += Zelement*Aelement*fractionVector[i]; //average with the number of nuclei
626 } 672 }
627 atot /= material->GetTotNbOfAtomsPerVolu 673 atot /= material->GetTotNbOfAtomsPerVolume();
628 zeff /= (material->GetTotNbOfAtomsPerVol 674 zeff /= (material->GetTotNbOfAtomsPerVolume()*atot);
629 << 675
630 intZ = (G4int) (zeff+0.25); 676 intZ = (G4int) (zeff+0.25);
631 if (intZ <= 0) 677 if (intZ <= 0)
632 intZ = 1; 678 intZ = 1;
633 if (intZ > fMaxZ) << 679 if (intZ > 99)
634 intZ = fMaxZ; << 680 intZ = 99;
635 } 681 }
636 682
637 if (fEffectiveCharge) 683 if (fEffectiveCharge)
638 fEffectiveCharge->insert(std::make_pair(ma 684 fEffectiveCharge->insert(std::make_pair(material,zeff));
639 685
640 // 686 //
641 // 2) Calculate Coulomb Correction 687 // 2) Calculate Coulomb Correction
642 // 688 //
643 G4double alz = fine_structure_const*zeff; 689 G4double alz = fine_structure_const*zeff;
644 G4double alzSquared = alz*alz; 690 G4double alzSquared = alz*alz;
645 G4double fc = alzSquared*(0.202059-alzSquar 691 G4double fc = alzSquared*(0.202059-alzSquared*
646 (0.03693-alzSquared* 692 (0.03693-alzSquared*
647 (0.00835-alzSquared*(0.00201-alzSq 693 (0.00835-alzSquared*(0.00201-alzSquared*
648 (0.00049-alzSquared* 694 (0.00049-alzSquared*
649 (0.00012-alzSquared*0.00003))) 695 (0.00012-alzSquared*0.00003)))))
650 +1.0/(alzSquared+1.0)); 696 +1.0/(alzSquared+1.0));
651 // 697 //
652 // 3) Screening functions and low-energy cor 698 // 3) Screening functions and low-energy corrections
653 // 699 //
654 G4double matRadius = 2.0/ fAtomicScreeningRa << 700 G4double matRadius = 2.0/ fAtomicScreeningRadius[intZ-1];
655 if (fMaterialInvScreeningRadius) 701 if (fMaterialInvScreeningRadius)
656 fMaterialInvScreeningRadius->insert(std::m 702 fMaterialInvScreeningRadius->insert(std::make_pair(material,matRadius));
657 703
658 std::pair<G4double,G4double> myPair(0,0); 704 std::pair<G4double,G4double> myPair(0,0);
659 G4double f0a = 4.0*G4Log(fAtomicScreeningRad << 705 G4double f0a = 4.0*std::log(fAtomicScreeningRadius[intZ-1]);
660 G4double f0b = f0a - 4.0*fc; 706 G4double f0b = f0a - 4.0*fc;
661 myPair.first = f0a; 707 myPair.first = f0a;
662 myPair.second = f0b; 708 myPair.second = f0b;
663 709
664 if (fScreeningFunction) 710 if (fScreeningFunction)
665 fScreeningFunction->insert(std::make_pair( 711 fScreeningFunction->insert(std::make_pair(material,myPair));
666 712
667 if (fVerboseLevel > 2) << 713 if (verboseLevel > 2)
668 { 714 {
669 G4cout << "Average Z for material " << m << 715 G4cout << "Average Z for material " << material->GetName() << " = " <<
670 zeff << G4endl; 716 zeff << G4endl;
671 G4cout << "Effective radius for material << 717 G4cout << "Effective radius for material " << material->GetName() << " = " <<
672 fAtomicScreeningRadius[intZ] << " m_e*c/hbar << 718 fAtomicScreeningRadius[intZ-1] << " m_e*c/hbar --> BCB = " <<
673 matRadius << G4endl; 719 matRadius << G4endl;
674 G4cout << "Screening parameters F0 for m << 720 G4cout << "Screening parameters F0 for material " << material->GetName() << " = " <<
675 f0a << "," << f0b << G4endl; 721 f0a << "," << f0b << G4endl;
676 } 722 }
677 return; 723 return;
678 } 724 }
679 725
680 //....oooOO0OOooo........oooOO0OOooo........oo 726 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
681 727
682 std::pair<G4double,G4double> << 728 std::pair<G4double,G4double>
683 G4PenelopeGammaConversionModel::GetScreeningFu 729 G4PenelopeGammaConversionModel::GetScreeningFunctions(G4double B)
684 { 730 {
685 // This is subroutine SCHIFF of Penelope 731 // This is subroutine SCHIFF of Penelope
686 // 732 //
687 // Screening Functions F1(B) and F2(B) in th << 733 // Screening Functions F1(B) and F2(B) in the Bethe-Heitler differential cross
688 // section for pair production 734 // section for pair production
689 // 735 //
690 std::pair<G4double,G4double> result(0.,0.); 736 std::pair<G4double,G4double> result(0.,0.);
691 G4double BSquared = B*B; 737 G4double BSquared = B*B;
692 G4double f1 = 2.0-2.0*G4Log(1.0+BSquared); << 738 G4double f1 = 2.0-2.0*std::log(1.0+BSquared);
693 G4double f2 = f1 - 6.66666666e-1; // (-2/3) 739 G4double f2 = f1 - 6.66666666e-1; // (-2/3)
694 if (B < 1.0e-10) 740 if (B < 1.0e-10)
695 f1 = f1-twopi*B; 741 f1 = f1-twopi*B;
696 else 742 else
697 { 743 {
698 G4double a0 = 4.0*B*std::atan(1./B); 744 G4double a0 = 4.0*B*std::atan(1./B);
699 f1 = f1 - a0; 745 f1 = f1 - a0;
700 f2 += 2.0*BSquared*(4.0-a0-3.0*G4Log((1. << 746 f2 += 2.0*BSquared*(4.0-a0-3.0*std::log((1.0+BSquared)/BSquared));
701 } 747 }
702 G4double g1 = 0.5*(3.0*f1-f2); 748 G4double g1 = 0.5*(3.0*f1-f2);
703 G4double g2 = 0.25*(3.0*f1+f2); 749 G4double g2 = 0.25*(3.0*f1+f2);
704 750
705 result.first = g1; 751 result.first = g1;
706 result.second = g2; 752 result.second = g2;
707 753
708 return result; 754 return result;
709 } 755 }
710 756
711 //....oooOO0OOooo........oooOO0OOooo........oo 757 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
712 758
713 void G4PenelopeGammaConversionModel::SetPartic 759 void G4PenelopeGammaConversionModel::SetParticle(const G4ParticleDefinition* p)
714 { 760 {
715 if(!fParticle) { 761 if(!fParticle) {
716 fParticle = p; << 762 fParticle = p;
717 } 763 }
718 } 764 }
>> 765
719 766