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