| 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 // Author: Sebastien Incerti 26 // Author: Sebastien Incerti
27 // 22 January 2012 27 // 22 January 2012
28 // on base of G4LivermoreNuclearGammaC 28 // on base of G4LivermoreNuclearGammaConversionModel (original version)
29 // and G4LivermoreRayleighModel (MT ve 29 // and G4LivermoreRayleighModel (MT version)
30 30
31 #include "G4LivermoreNuclearGammaConversionMod 31 #include "G4LivermoreNuclearGammaConversionModel.hh"
32 #include "G4PhysicalConstants.hh" 32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh"
34 #include "G4Log.hh" 34 #include "G4Log.hh"
35 #include "G4Exp.hh" 35 #include "G4Exp.hh"
36 #include "G4AutoLock.hh" <<
37 36
38 //....oooOO0OOooo........oooOO0OOooo........oo 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 38
40 using namespace std; 39 using namespace std;
41 namespace { G4Mutex LivermoreNuclearGammaConve <<
42 40
43 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 42
45 G4PhysicsFreeVector* G4LivermoreNuclearGammaCo << 43 G4int G4LivermoreNuclearGammaConversionModel::maxZ = 100;
>> 44 G4LPhysicsFreeVector* G4LivermoreNuclearGammaConversionModel::data[] = {0};
46 45
47 G4LivermoreNuclearGammaConversionModel::G4Live 46 G4LivermoreNuclearGammaConversionModel::G4LivermoreNuclearGammaConversionModel
48 (const G4ParticleDefinition*, const G4String& 47 (const G4ParticleDefinition*, const G4String& nam)
49 :G4VEmModel(nam),smallEnergy(2.*MeV), << 48 :G4VEmModel(nam),isInitialised(false),smallEnergy(2.*MeV)
50 isInitialised(false) <<
51 { 49 {
52 fParticleChange = nullptr; << 50 fParticleChange = 0;
53 51
54 lowEnergyLimit = 2.0*electron_mass_c2; 52 lowEnergyLimit = 2.0*electron_mass_c2;
55 53
56 verboseLevel= 0; 54 verboseLevel= 0;
57 // Verbosity scale for debugging purposes: 55 // Verbosity scale for debugging purposes:
58 // 0 = nothing 56 // 0 = nothing
59 // 1 = calculation of cross sections, file o 57 // 1 = calculation of cross sections, file openings...
60 // 2 = entering in methods 58 // 2 = entering in methods
61 59
62 if(verboseLevel > 0) 60 if(verboseLevel > 0)
63 { 61 {
64 G4cout << "G4LivermoreNuclearGammaConversi 62 G4cout << "G4LivermoreNuclearGammaConversionModel is constructed " << G4endl;
65 } 63 }
66 } 64 }
67 65
68 //....oooOO0OOooo........oooOO0OOooo........oo 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 67
70 G4LivermoreNuclearGammaConversionModel::~G4Liv 68 G4LivermoreNuclearGammaConversionModel::~G4LivermoreNuclearGammaConversionModel()
71 { 69 {
72 if(IsMaster()) { 70 if(IsMaster()) {
73 for(G4int i=0; i<maxZ; ++i) { 71 for(G4int i=0; i<maxZ; ++i) {
74 if(data[i]) { 72 if(data[i]) {
75 delete data[i]; 73 delete data[i];
76 data[i] = 0; 74 data[i] = 0;
77 } 75 }
78 } 76 }
79 } 77 }
80 } 78 }
81 79
82 //....oooOO0OOooo........oooOO0OOooo........oo 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 81
84 void G4LivermoreNuclearGammaConversionModel::I 82 void G4LivermoreNuclearGammaConversionModel::Initialise(
85 const G4Partic 83 const G4ParticleDefinition* particle,
86 const G4DataVector& cuts) 84 const G4DataVector& cuts)
87 { 85 {
>> 86
88 if (verboseLevel > 1) 87 if (verboseLevel > 1)
89 { << 88 {
90 G4cout << "Calling Initialise() of G4Liver 89 G4cout << "Calling Initialise() of G4LivermoreNuclearGammaConversionModel."
91 << G4endl 90 << G4endl
92 << "Energy range: " 91 << "Energy range: "
93 << LowEnergyLimit() / MeV << " MeV - " 92 << LowEnergyLimit() / MeV << " MeV - "
94 << HighEnergyLimit() / GeV << " GeV" 93 << HighEnergyLimit() / GeV << " GeV"
95 << G4endl; 94 << G4endl;
96 } << 95 }
97 96
98 if(IsMaster()) 97 if(IsMaster())
99 { 98 {
100 99
101 // Initialise element selector 100 // Initialise element selector
>> 101
102 InitialiseElementSelectors(particle, cuts) 102 InitialiseElementSelectors(particle, cuts);
103 103
104 // Access to elements << 104 // Access to elements
105 const char* path = G4FindDataDir("G4LEDATA << 105
>> 106 char* path = std::getenv("G4LEDATA");
106 107
107 G4ProductionCutsTable* theCoupleTable = 108 G4ProductionCutsTable* theCoupleTable =
108 G4ProductionCutsTable::GetProductionCuts 109 G4ProductionCutsTable::GetProductionCutsTable();
109 110
110 G4int numOfCouples = (G4int)theCoupleTable << 111 G4int numOfCouples = theCoupleTable->GetTableSize();
111 112
112 for(G4int i=0; i<numOfCouples; ++i) 113 for(G4int i=0; i<numOfCouples; ++i)
113 { 114 {
114 const G4Material* material = 115 const G4Material* material =
115 theCoupleTable->GetMaterialCutsCouple( 116 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
116 const G4ElementVector* theElementVector 117 const G4ElementVector* theElementVector = material->GetElementVector();
117 std::size_t nelm = material->GetNumberOf << 118 G4int nelm = material->GetNumberOfElements();
118 119
119 for (std::size_t j=0; j<nelm; ++j) << 120 for (G4int j=0; j<nelm; ++j)
120 { << 121 {
121 G4int Z = (G4int)(*theElementVector)[j << 122 G4int Z = (G4int)(*theElementVector)[j]->GetZ();
122 if(Z < 1) { Z = 1; } << 123 if(Z < 1) { Z = 1; }
123 else if(Z > maxZ) { Z = maxZ; } << 124 else if(Z > maxZ) { Z = maxZ; }
124 if(!data[Z]) { ReadData(Z, path); } << 125 if(!data[Z]) { ReadData(Z, path); }
125 } << 126 }
126 } 127 }
127 } 128 }
128 if(isInitialised) { return; } 129 if(isInitialised) { return; }
129 fParticleChange = GetParticleChangeForGamma( 130 fParticleChange = GetParticleChangeForGamma();
130 isInitialised = true; 131 isInitialised = true;
131 } 132 }
132 133
133 //....oooOO0OOooo........oooOO0OOooo........oo 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
134 135
135 void G4LivermoreNuclearGammaConversionModel::I 136 void G4LivermoreNuclearGammaConversionModel::InitialiseLocal(
136 const G4ParticleDefinition*, G4VEmModel* 137 const G4ParticleDefinition*, G4VEmModel* masterModel)
137 { 138 {
138 SetElementSelectors(masterModel->GetElementS 139 SetElementSelectors(masterModel->GetElementSelectors());
139 } 140 }
140 141
141 //....oooOO0OOooo........oooOO0OOooo........oo 142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142 143
143 G4double 144 G4double
144 G4LivermoreNuclearGammaConversionModel::MinPri 145 G4LivermoreNuclearGammaConversionModel::MinPrimaryEnergy(const G4Material*,
145 const G4ParticleDefinition*, 146 const G4ParticleDefinition*,
146 G4double) 147 G4double)
147 { 148 {
148 return lowEnergyLimit; 149 return lowEnergyLimit;
149 } 150 }
150 151
151 //....oooOO0OOooo........oooOO0OOooo........oo 152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152 153
153 void G4LivermoreNuclearGammaConversionModel::R 154 void G4LivermoreNuclearGammaConversionModel::ReadData(size_t Z, const char* path)
154 { 155 {
155 if (verboseLevel > 1) 156 if (verboseLevel > 1)
156 { 157 {
157 G4cout << "Calling ReadData() of G4Livermo 158 G4cout << "Calling ReadData() of G4LivermoreNuclearGammaConversionModel"
158 << G4endl; 159 << G4endl;
159 } 160 }
160 161
161 162
162 if(data[Z]) { return; } 163 if(data[Z]) { return; }
163 164
164 const char* datadir = path; 165 const char* datadir = path;
165 166
166 if(!datadir) 167 if(!datadir)
167 { 168 {
168 datadir = G4FindDataDir("G4LEDATA"); << 169 datadir = std::getenv("G4LEDATA");
169 if(!datadir) 170 if(!datadir)
170 { 171 {
171 G4Exception("G4LivermoreNuclearGammaConv 172 G4Exception("G4LivermoreNuclearGammaConversionModel::ReadData()",
172 "em0006",FatalException, 173 "em0006",FatalException,
173 "Environment variable G4LEDATA not defin 174 "Environment variable G4LEDATA not defined");
174 return; 175 return;
175 } 176 }
176 } 177 }
177 178
178 data[Z] = new G4PhysicsFreeVector(0,/*spline << 179 //
>> 180
>> 181 data[Z] = new G4LPhysicsFreeVector();
>> 182
>> 183 //
179 184
180 std::ostringstream ost; 185 std::ostringstream ost;
181 ost << datadir << "/livermore/pairdata/pp-pa 186 ost << datadir << "/livermore/pairdata/pp-pair-cs-" << Z <<".dat";
182 std::ifstream fin(ost.str().c_str()); 187 std::ifstream fin(ost.str().c_str());
183 188
184 if( !fin.is_open()) 189 if( !fin.is_open())
185 { 190 {
186 G4ExceptionDescription ed; 191 G4ExceptionDescription ed;
187 ed << "G4LivermoreNuclearGammaConversionMo 192 ed << "G4LivermoreNuclearGammaConversionModel data file <" << ost.str().c_str()
188 << "> is not opened!" << G4endl; 193 << "> is not opened!" << G4endl;
189 G4Exception("G4LivermoreNuclearGammaConver 194 G4Exception("G4LivermoreNuclearGammaConversionModel::ReadData()",
190 "em0003",FatalException, 195 "em0003",FatalException,
191 ed,"G4LEDATA version should be G4EMLOW8.0 << 196 ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
192 return; 197 return;
193 } 198 }
194 else <<
195 { <<
196 <<
197 if(verboseLevel > 3) { G4cout << "File " <<
198 << " is opened by G4LivermoreNucle <<
199 <<
200 data[Z]->Retrieve(fin, true); <<
201 } <<
202 199
>> 200 else
>> 201 {
>> 202
>> 203 if(verboseLevel > 3) { G4cout << "File " << ost.str()
>> 204 << " is opened by G4LivermoreNuclearGammaConversionModel" << G4endl;}
>> 205
>> 206 data[Z]->Retrieve(fin, true);
>> 207 }
>> 208
203 // Activation of spline interpolation 209 // Activation of spline interpolation
204 data[Z] ->FillSecondDerivatives(); << 210 data[Z] ->SetSpline(true);
>> 211
205 } 212 }
206 213
207 //....oooOO0OOooo........oooOO0OOooo........oo 214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
208 215
209 G4double 216 G4double
210 G4LivermoreNuclearGammaConversionModel::Comput 217 G4LivermoreNuclearGammaConversionModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
211 G4double GammaEnergy, 218 G4double GammaEnergy,
212 G4double Z, G4double, 219 G4double Z, G4double,
213 G4double, G4double) 220 G4double, G4double)
214 { 221 {
215 if (verboseLevel > 1) 222 if (verboseLevel > 1)
216 { 223 {
217 G4cout << "Calling ComputeCrossSectionPerA 224 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreNuclearGammaConversionModel"
218 << G4endl; 225 << G4endl;
219 } 226 }
220 227
221 if (GammaEnergy < lowEnergyLimit) { return 0 228 if (GammaEnergy < lowEnergyLimit) { return 0.0; }
222 229
223 G4double xs = 0.0; 230 G4double xs = 0.0;
224 231
225 G4int intZ=G4int(Z); 232 G4int intZ=G4int(Z);
226 233
227 if(intZ < 1 || intZ > maxZ) { return xs; } 234 if(intZ < 1 || intZ > maxZ) { return xs; }
228 235
229 G4PhysicsFreeVector* pv = data[intZ]; << 236 G4LPhysicsFreeVector* pv = data[intZ];
230 237
231 // if element was not initialised 238 // if element was not initialised
232 // do initialisation safely for MT mode 239 // do initialisation safely for MT mode
233 if(!pv) 240 if(!pv)
234 { 241 {
235 InitialiseForElement(0, intZ); 242 InitialiseForElement(0, intZ);
236 pv = data[intZ]; 243 pv = data[intZ];
237 if(!pv) { return xs; } 244 if(!pv) { return xs; }
238 } 245 }
239 // x-section is taken from the table 246 // x-section is taken from the table
240 xs = pv->Value(GammaEnergy); 247 xs = pv->Value(GammaEnergy);
241 248
242 if(verboseLevel > 0) 249 if(verboseLevel > 0)
243 { << 250 {
244 std::size_t n = pv->GetVectorLength() - 1; << 251 G4int n = pv->GetVectorLength() - 1;
245 G4cout << "****** DEBUG: tcs value for Z 252 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
246 << GammaEnergy/MeV << G4endl; 253 << GammaEnergy/MeV << G4endl;
247 G4cout << " cs (Geant4 internal unit)=" 254 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
248 G4cout << " -> first cs value in EADL 255 G4cout << " -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl;
249 G4cout << " -> last cs value in EADL 256 G4cout << " -> last cs value in EADL data file (iu) =" << (*pv)[n] << G4endl;
250 G4cout << "***************************** 257 G4cout << "*********************************************************" << G4endl;
251 } << 258 }
>> 259
252 return xs; 260 return xs;
>> 261
253 } 262 }
254 263
255 //....oooOO0OOooo........oooOO0OOooo........oo 264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
256 265
257 void G4LivermoreNuclearGammaConversionModel::S 266 void G4LivermoreNuclearGammaConversionModel::SampleSecondaries(
258 std::vector<G 267 std::vector<G4DynamicParticle*>* fvect,
259 const G4MaterialCutsCouple* couple, 268 const G4MaterialCutsCouple* couple,
260 const G4DynamicParticle* aDynamicGamm 269 const G4DynamicParticle* aDynamicGamma,
261 G4double, G4double) 270 G4double, G4double)
262 { 271 {
263 // The energies of the e+ e- secondaries are << 272
264 // cross sections with Coulomb correction. A << 273 // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
265 // number techniques of Butcher & Messel is << 274 // cross sections with Coulomb correction. A modified version of the random
266 << 275 // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
267 // Note 1 : Effects due to the breakdown of << 276
268 // energy are ignored. << 277 // Note 1 : Effects due to the breakdown of the Born approximation at low
269 // Note 2 : The differential cross section i << 278 // energy are ignored.
270 // pair creation in both nuclear and atomic << 279 // Note 2 : The differential cross section implicitly takes account of
271 // prodution is not generated. << 280 // pair creation in both nuclear and atomic electron fields. However triplet
272 << 281 // prodution is not generated.
>> 282
273 if (verboseLevel > 1) { 283 if (verboseLevel > 1) {
274 G4cout << "Calling SampleSecondaries() of 284 G4cout << "Calling SampleSecondaries() of G4LivermoreNuclearGammaConversionModel"
275 << G4endl; 285 << G4endl;
276 } 286 }
277 287
278 G4double photonEnergy = aDynamicGamma->GetKi 288 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
279 G4ParticleMomentum photonDirection = aDynami 289 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
280 290
281 G4double epsilon ; 291 G4double epsilon ;
282 G4double epsilon0Local = electron_mass_c2 / 292 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
283 293
284 // Do it fast if photon energy < 2. MeV 294 // Do it fast if photon energy < 2. MeV
285 if (photonEnergy < smallEnergy ) 295 if (photonEnergy < smallEnergy )
286 { 296 {
287 epsilon = epsilon0Local + (0.5 - epsilon0L 297 epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
288 } 298 }
289 else 299 else
290 { 300 {
291 // Select randomly one element in the curr 301 // Select randomly one element in the current material
>> 302
292 const G4ParticleDefinition* particle = aD 303 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
293 const G4Element* element = SelectRandomAto 304 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
294 305
295 if (element == nullptr) << 306 if (element == 0)
296 { 307 {
297 G4cout << "G4LivermoreNuclearGammaConversion 308 G4cout << "G4LivermoreNuclearGammaConversionModel::SampleSecondaries - element = 0"
298 << G4endl; 309 << G4endl;
299 return; 310 return;
300 } 311 }
301 G4IonisParamElm* ionisation = element->Get 312 G4IonisParamElm* ionisation = element->GetIonisation();
302 if (ionisation == nullptr) << 313 if (ionisation == 0)
303 { 314 {
304 G4cout << "G4LivermoreNuclearGammaConversion 315 G4cout << "G4LivermoreNuclearGammaConversionModel::SampleSecondaries - ionisation = 0"
305 << G4endl; 316 << G4endl;
306 return; 317 return;
307 } 318 }
308 319
309 // Extract Coulomb factor for this Element 320 // Extract Coulomb factor for this Elements
310 G4double fZ = 8. * (ionisation->GetlogZ3() 321 G4double fZ = 8. * (ionisation->GetlogZ3());
311 if (photonEnergy > 50. * MeV) fZ += 8. * ( 322 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
312 323
313 // Limits of the screening variable 324 // Limits of the screening variable
314 G4double screenFactor = 136. * epsilon0Loc 325 G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
315 G4double screenMax = G4Exp ((42.24 - fZ)/8 326 G4double screenMax = G4Exp ((42.24 - fZ)/8.368) - 0.952 ;
316 G4double screenMin = std::min(4.*screenFac 327 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
317 328
318 // Limits of the energy sampling 329 // Limits of the energy sampling
319 G4double epsilon1 = 0.5 - 0.5 * std::sqrt( 330 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
320 G4double epsilonMin = std::max(epsilon0Loc 331 G4double epsilonMin = std::max(epsilon0Local,epsilon1);
321 G4double epsilonRange = 0.5 - epsilonMin ; 332 G4double epsilonRange = 0.5 - epsilonMin ;
322 333
323 // Sample the energy rate of the created e 334 // Sample the energy rate of the created electron (or positron)
324 G4double screen; 335 G4double screen;
325 G4double gReject ; 336 G4double gReject ;
326 337
327 G4double f10 = ScreenFunction1(screenMin) 338 G4double f10 = ScreenFunction1(screenMin) - fZ;
328 G4double f20 = ScreenFunction2(screenMin) 339 G4double f20 = ScreenFunction2(screenMin) - fZ;
329 G4double normF1 = std::max(f10 * epsilonRa 340 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
330 G4double normF2 = std::max(1.5 * f20,0.); 341 G4double normF2 = std::max(1.5 * f20,0.);
331 342
332 do 343 do
333 { 344 {
334 if (normF1 / (normF1 + normF2) > G4UniformRa 345 if (normF1 / (normF1 + normF2) > G4UniformRand() )
335 { 346 {
336 epsilon = 0.5 - epsilonRange * std::pow( 347 epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.333333) ;
337 screen = screenFactor / (epsilon * (1. - 348 screen = screenFactor / (epsilon * (1. - epsilon));
338 gReject = (ScreenFunction1(screen) - fZ) 349 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
339 } 350 }
340 else 351 else
341 { 352 {
342 epsilon = epsilonMin + epsilonRange * G4 353 epsilon = epsilonMin + epsilonRange * G4UniformRand();
343 screen = screenFactor / (epsilon * (1 - 354 screen = screenFactor / (epsilon * (1 - epsilon));
344 gReject = (ScreenFunction2(screen) - fZ) 355 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
345 } 356 }
346 } while ( gReject < G4UniformRand() ); << 357 } while ( gReject < G4UniformRand() );
>> 358
347 } // End of epsilon sampling 359 } // End of epsilon sampling
348 360
349 // Fix charges randomly 361 // Fix charges randomly
>> 362
350 G4double electronTotEnergy; 363 G4double electronTotEnergy;
351 G4double positronTotEnergy; 364 G4double positronTotEnergy;
352 365
353 if (G4UniformRand() > 0.5) 366 if (G4UniformRand() > 0.5)
354 { 367 {
355 electronTotEnergy = (1. - epsilon) * pho 368 electronTotEnergy = (1. - epsilon) * photonEnergy;
356 positronTotEnergy = epsilon * photonEner 369 positronTotEnergy = epsilon * photonEnergy;
357 } 370 }
358 else 371 else
359 { 372 {
360 positronTotEnergy = (1. - epsilon) * pho 373 positronTotEnergy = (1. - epsilon) * photonEnergy;
361 electronTotEnergy = epsilon * photonEner 374 electronTotEnergy = epsilon * photonEnergy;
362 } 375 }
363 376
364 // Scattered electron (positron) angles. ( Z 377 // Scattered electron (positron) angles. ( Z - axis along the parent photon)
365 // Universal distribution suggested by L. Ur 378 // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
366 // derived from Tsai distribution (Rev. Mod. 379 // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
367 << 380
368 G4double u; 381 G4double u;
369 const G4double a1 = 0.625; 382 const G4double a1 = 0.625;
370 G4double a2 = 3. * a1; 383 G4double a2 = 3. * a1;
>> 384 // G4double d = 27. ;
371 385
>> 386 // if (9. / (9. + d) > G4UniformRand())
372 if (0.25 > G4UniformRand()) 387 if (0.25 > G4UniformRand())
373 { 388 {
374 u = - G4Log(G4UniformRand() * G4UniformR 389 u = - G4Log(G4UniformRand() * G4UniformRand()) / a1 ;
375 } 390 }
376 else 391 else
377 { 392 {
378 u = - G4Log(G4UniformRand() * G4UniformR 393 u = - G4Log(G4UniformRand() * G4UniformRand()) / a2 ;
379 } 394 }
380 395
381 G4double thetaEle = u*electron_mass_c2/elect 396 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
382 G4double thetaPos = u*electron_mass_c2/posit 397 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
383 G4double phi = twopi * G4UniformRand(); 398 G4double phi = twopi * G4UniformRand();
384 399
385 G4double dxEle= std::sin(thetaEle)*std::cos( 400 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
386 G4double dxPos=-std::sin(thetaPos)*std::cos( 401 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
387 << 402
>> 403
388 // Kinematics of the created pair: 404 // Kinematics of the created pair:
389 // the electron and positron are assumed to 405 // the electron and positron are assumed to have a symetric angular
390 // distribution with respect to the Z axis a 406 // distribution with respect to the Z axis along the parent photon
391 407
392 G4double electronKineEnergy = std::max(0.,el 408 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
393 409
394 G4ThreeVector electronDirection (dxEle, dyEl 410 G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
395 electronDirection.rotateUz(photonDirection); 411 electronDirection.rotateUz(photonDirection);
396 412
397 G4DynamicParticle* particle1 = new G4Dynamic 413 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
398 electronDirection, 414 electronDirection,
399 electronKineEnergy); 415 electronKineEnergy);
400 416
401 // The e+ is always created 417 // The e+ is always created
402 G4double positronKineEnergy = std::max(0.,po 418 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
403 419
404 G4ThreeVector positronDirection (dxPos, dyPo 420 G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
405 positronDirection.rotateUz(photonDirection); 421 positronDirection.rotateUz(photonDirection);
406 422
407 // Create G4DynamicParticle object for the p 423 // Create G4DynamicParticle object for the particle2
408 G4DynamicParticle* particle2 = new G4Dynamic 424 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
409 positronDirection, 425 positronDirection,
410 positronKineEnergy); 426 positronKineEnergy);
411 // Fill output vector 427 // Fill output vector
412 fvect->push_back(particle1); 428 fvect->push_back(particle1);
413 fvect->push_back(particle2); 429 fvect->push_back(particle2);
414 430
415 // kill incident photon 431 // kill incident photon
416 fParticleChange->SetProposedKineticEnergy(0. 432 fParticleChange->SetProposedKineticEnergy(0.);
417 fParticleChange->ProposeTrackStatus(fStopAnd 433 fParticleChange->ProposeTrackStatus(fStopAndKill);
418 434
419 } 435 }
420 436
421 //....oooOO0OOooo........oooOO0OOooo........oo 437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
422 438
423 G4double 439 G4double
424 G4LivermoreNuclearGammaConversionModel::Screen 440 G4LivermoreNuclearGammaConversionModel::ScreenFunction1(G4double screenVariable)
425 { 441 {
426 // Compute the value of the screening functi 442 // Compute the value of the screening function 3*phi1 - phi2
427 443
428 G4double value; 444 G4double value;
429 445
430 if (screenVariable > 1.) 446 if (screenVariable > 1.)
431 value = 42.24 - 8.368 * G4Log(screenVariab 447 value = 42.24 - 8.368 * G4Log(screenVariable + 0.952);
432 else 448 else
433 value = 42.392 - screenVariable * (7.796 - 449 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
434 450
435 return value; 451 return value;
436 } 452 }
437 453
438 //....oooOO0OOooo........oooOO0OOooo........oo 454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
439 455
440 G4double 456 G4double
441 G4LivermoreNuclearGammaConversionModel::Screen 457 G4LivermoreNuclearGammaConversionModel::ScreenFunction2(G4double screenVariable)
442 { 458 {
443 // Compute the value of the screening functi 459 // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
>> 460
444 G4double value; 461 G4double value;
445 462
446 if (screenVariable > 1.) 463 if (screenVariable > 1.)
447 value = 42.24 - 8.368 * G4Log(screenVariab 464 value = 42.24 - 8.368 * G4Log(screenVariable + 0.952);
448 else 465 else
449 value = 41.405 - screenVariable * (5.828 - 466 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
450 467
451 return value; 468 return value;
452 } 469 }
453 470
454 //....oooOO0OOooo........oooOO0OOooo........oo 471 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
455 472
>> 473 #include "G4AutoLock.hh"
>> 474 namespace { G4Mutex LivermoreNuclearGammaConversionModelMutex = G4MUTEX_INITIALIZER; }
>> 475
456 void G4LivermoreNuclearGammaConversionModel::I 476 void G4LivermoreNuclearGammaConversionModel::InitialiseForElement(
457 const G4ParticleDefinition 477 const G4ParticleDefinition*,
458 G4int Z) 478 G4int Z)
459 { 479 {
460 G4AutoLock l(&LivermoreNuclearGammaConversio << 480 G4AutoLock l(&LivermoreNuclearGammaConversionModelMutex);
>> 481 // G4cout << "G4LivermoreNuclearGammaConversionModel::InitialiseForElement Z= "
>> 482 // << Z << G4endl;
461 if(!data[Z]) { ReadData(Z); } 483 if(!data[Z]) { ReadData(Z); }
462 l.unlock(); 484 l.unlock();
463 } 485 }
464 486
465 //....oooOO0OOooo........oooOO0OOooo........oo 487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
466 488