| 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: G4LivermoreGammaConversionModel.cc,v 1.9 2010-12-27 17:45:12 vnivanch Exp $
>> 27 // GEANT4 tag $Name: not supported by cvs2svn $
>> 28 //
>> 29 //
26 // Author: Sebastien Incerti 30 // Author: Sebastien Incerti
27 // 22 January 2012 << 31 // 30 October 2008
28 // on base of G4LivermoreGammaConversi << 32 // on base of G4LowEnergyGammaConversion developed by A.Forti and M.G.Pia
29 // and G4LivermoreRayleighModel (MT ve <<
30 // 33 //
31 // Modifications: Zhuxin Li@CENBG << 34 // History:
32 // 11 March 2020 << 35 // --------
33 // derives from G4PairProductio << 36 // 12 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
34 // ------------------------------------------- << 37 // - apply internal high-energy limit only in constructor
>> 38 // - do not apply low-energy limit (default is 0)
>> 39 // - use CLHEP electron mass for low-enegry limit
>> 40 // - remove MeanFreePath method and table
>> 41 // 26 Dec 2010 V Ivanchenko Load data tables only once to avoid memory leak
>> 42
35 43
36 #include "G4LivermoreGammaConversionModel.hh" 44 #include "G4LivermoreGammaConversionModel.hh"
37 45
38 #include "G4AutoLock.hh" << 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 #include "G4Electron.hh" <<
40 #include "G4EmParameters.hh" <<
41 #include "G4Exp.hh" <<
42 #include "G4ParticleChangeForGamma.hh" <<
43 #include "G4PhysicalConstants.hh" <<
44 #include "G4PhysicsFreeVector.hh" <<
45 #include "G4SystemOfUnits.hh" <<
46 47
47 namespace << 48 using namespace std;
48 { <<
49 G4Mutex LivermoreGammaConversionModelMutex = G <<
50 } <<
51 49
52 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53 //....oooOO0OOooo........oooOO0OOooo........oo <<
54 <<
55 G4PhysicsFreeVector* G4LivermoreGammaConversio <<
56 G4String G4LivermoreGammaConversionModel::gDat <<
57 51
58 G4LivermoreGammaConversionModel::G4LivermoreGa << 52 G4LivermoreGammaConversionModel::G4LivermoreGammaConversionModel(const G4ParticleDefinition*,
59 << 53 const G4String& nam)
60 : G4PairProductionRelModel(p, nam) << 54 :G4VEmModel(nam),smallEnergy(2.*MeV),isInitialised(false),
61 { << 55 crossSectionHandler(0),meanFreePathTable(0)
62 fParticleChange = nullptr; << 56 {
63 lowEnergyLimit = 2. * CLHEP::electron_mass_c << 57 lowEnergyLimit = 2.0*electron_mass_c2;
64 verboseLevel = 0; << 58 highEnergyLimit = 100 * GeV;
65 // Verbosity scale for debugging purposes: << 59 SetHighEnergyLimit(highEnergyLimit);
66 // 0 = nothing << 60
67 // 1 = calculation of cross sections, file o << 61 verboseLevel= 0;
68 // 2 = entering in methods << 62 // Verbosity scale:
69 if (verboseLevel > 0) { << 63 // 0 = nothing
70 G4cout << "G4LivermoreGammaConversionModel << 64 // 1 = warning for energy non-conservation
>> 65 // 2 = details of energy budget
>> 66 // 3 = calculation of cross sections, file openings, sampling of atoms
>> 67 // 4 = entering in methods
>> 68
>> 69 if(verboseLevel > 0) {
>> 70 G4cout << "Livermore Gamma conversion is constructed " << G4endl
>> 71 << "Energy range: "
>> 72 << lowEnergyLimit / MeV << " MeV - "
>> 73 << highEnergyLimit / GeV << " GeV"
>> 74 << G4endl;
71 } 75 }
72 } 76 }
73 77
74 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 79
76 G4LivermoreGammaConversionModel::~G4LivermoreG 80 G4LivermoreGammaConversionModel::~G4LivermoreGammaConversionModel()
77 { << 81 {
78 if (IsMaster()) { << 82 if (crossSectionHandler) { delete crossSectionHandler; }
79 for (G4int i = 0; i <= maxZ; ++i) { <<
80 if (data[i]) { <<
81 delete data[i]; <<
82 data[i] = nullptr; <<
83 } <<
84 } <<
85 } <<
86 } 83 }
87 84
88 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 86
90 void G4LivermoreGammaConversionModel::Initiali << 87 void
91 << 88 G4LivermoreGammaConversionModel::Initialise(const G4ParticleDefinition*,
>> 89 const G4DataVector&)
92 { 90 {
93 G4PairProductionRelModel::Initialise(particl << 91 if (verboseLevel > 3) {
94 if (verboseLevel > 1) { << 92 G4cout << "Calling G4LivermoreGammaConversionModel::Initialise()" << G4endl;
95 G4cout << "Calling Initialise() of G4Liver <<
96 << "Energy range: " << LowEnergyLim <<
97 << " GeV isMater: " << IsMaster() < <<
98 } <<
99 <<
100 if (IsMaster()) { <<
101 // Initialise element selector <<
102 InitialiseElementSelectors(particle, cuts) <<
103 <<
104 // Access to elements <<
105 const G4ElementTable* elemTable = G4Elemen <<
106 std::size_t numElems = (*elemTable).size() <<
107 for (std::size_t ie = 0; ie < numElems; ++ <<
108 const G4Element* elem = (*elemTable)[ie] <<
109 const G4int Z = std::min(maxZ, elem->Get <<
110 if (data[Z] == nullptr) { <<
111 ReadData(Z); <<
112 } <<
113 } <<
114 } 93 }
115 if (isInitialised) { << 94
116 return; << 95 if (crossSectionHandler)
>> 96 {
>> 97 crossSectionHandler->Clear();
>> 98 delete crossSectionHandler;
>> 99 }
>> 100
>> 101 // Read data tables for all materials
>> 102
>> 103 crossSectionHandler = new G4CrossSectionHandler();
>> 104 crossSectionHandler->Initialise(0,lowEnergyLimit,100.*GeV,400);
>> 105 G4String crossSectionFile = "pair/pp-cs-";
>> 106 crossSectionHandler->LoadData(crossSectionFile);
>> 107 //
>> 108
>> 109 if (verboseLevel > 2) {
>> 110 G4cout << "Loaded cross section files for Livermore Gamma Conversion model" << G4endl;
>> 111 }
>> 112 if (verboseLevel > 0) {
>> 113 G4cout << "Livermore Gamma Conversion model is initialized " << G4endl
>> 114 << "Energy range: "
>> 115 << LowEnergyLimit() / MeV << " MeV - "
>> 116 << HighEnergyLimit() / GeV << " GeV"
>> 117 << G4endl;
117 } 118 }
>> 119
>> 120 if(isInitialised) { return; }
118 fParticleChange = GetParticleChangeForGamma( 121 fParticleChange = GetParticleChangeForGamma();
119 isInitialised = true; 122 isInitialised = true;
120 } 123 }
121 124
122 //....oooOO0OOooo........oooOO0OOooo........oo 125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 126
124 const G4String& G4LivermoreGammaConversionMode << 127 G4double
>> 128 G4LivermoreGammaConversionModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
>> 129 G4double GammaEnergy,
>> 130 G4double Z, G4double,
>> 131 G4double, G4double)
125 { 132 {
126 // no check in this method - environment var << 133 if (verboseLevel > 3) {
127 if (gDataDirectory.empty()) { << 134 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModel"
128 auto param = G4EmParameters::Instance(); << 135 << G4endl;
129 std::ostringstream ost; <<
130 if (param->LivermoreDataDir() == "livermor <<
131 ost << param->GetDirLEDATA() << "/liverm <<
132 useSpline = true; <<
133 } <<
134 else { <<
135 ost << param->GetDirLEDATA() << "/epics2 <<
136 } <<
137 gDataDirectory = ost.str(); <<
138 } 136 }
139 return gDataDirectory; << 137 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) { return 0.0; }
>> 138
>> 139 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
>> 140 return cs;
140 } 141 }
141 142
142 //....oooOO0OOooo........oooOO0OOooo........oo 143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
143 144
144 void G4LivermoreGammaConversionModel::ReadData << 145 void G4LivermoreGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
>> 146 const G4MaterialCutsCouple* couple,
>> 147 const G4DynamicParticle* aDynamicGamma,
>> 148 G4double,
>> 149 G4double)
145 { 150 {
146 if (verboseLevel > 1) { <<
147 G4cout << "Calling ReadData() of G4Livermo <<
148 } <<
149 <<
150 if (data[Z] != nullptr) { <<
151 return; <<
152 } <<
153 <<
154 std::ostringstream ost; <<
155 ost << FindDirectoryPath() << "pp-cs-" << Z <<
156 151
157 data[Z] = new G4PhysicsFreeVector(useSpline) << 152 // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
158 << 153 // cross sections with Coulomb correction. A modified version of the random
159 std::ifstream fin(ost.str().c_str()); << 154 // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
>> 155
>> 156 // Note 1 : Effects due to the breakdown of the Born approximation at low
>> 157 // energy are ignored.
>> 158 // Note 2 : The differential cross section implicitly takes account of
>> 159 // pair creation in both nuclear and atomic electron fields. However triplet
>> 160 // prodution is not generated.
>> 161
>> 162 if (verboseLevel > 3)
>> 163 G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModel" << G4endl;
>> 164
>> 165 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
>> 166 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
>> 167
>> 168 G4double epsilon ;
>> 169 G4double epsilon0 = electron_mass_c2 / photonEnergy ;
>> 170
>> 171 // Do it fast if photon energy < 2. MeV
>> 172 if (photonEnergy < smallEnergy )
>> 173 {
>> 174 epsilon = epsilon0 + (0.5 - epsilon0) * G4UniformRand();
>> 175 }
>> 176 else
>> 177 {
>> 178 // Select randomly one element in the current material
>> 179 //const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
>> 180 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
>> 181 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
>> 182
>> 183 if (element == 0)
>> 184 {
>> 185 G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - element = 0"
>> 186 << G4endl;
>> 187 return;
>> 188 }
>> 189 G4IonisParamElm* ionisation = element->GetIonisation();
>> 190 if (ionisation == 0)
>> 191 {
>> 192 G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0"
>> 193 << G4endl;
>> 194 return;
>> 195 }
>> 196
>> 197 // Extract Coulomb factor for this Element
>> 198 G4double fZ = 8. * (ionisation->GetlogZ3());
>> 199 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
>> 200
>> 201 // Limits of the screening variable
>> 202 G4double screenFactor = 136. * epsilon0 / (element->GetIonisation()->GetZ3()) ;
>> 203 G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
>> 204 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
>> 205
>> 206 // Limits of the energy sampling
>> 207 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
>> 208 G4double epsilonMin = std::max(epsilon0,epsilon1);
>> 209 G4double epsilonRange = 0.5 - epsilonMin ;
>> 210
>> 211 // Sample the energy rate of the created electron (or positron)
>> 212 G4double screen;
>> 213 G4double gReject ;
>> 214
>> 215 G4double f10 = ScreenFunction1(screenMin) - fZ;
>> 216 G4double f20 = ScreenFunction2(screenMin) - fZ;
>> 217 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
>> 218 G4double normF2 = std::max(1.5 * f20,0.);
>> 219
>> 220 do {
>> 221 if (normF1 / (normF1 + normF2) > G4UniformRand() )
>> 222 {
>> 223 epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.3333) ;
>> 224 screen = screenFactor / (epsilon * (1. - epsilon));
>> 225 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
>> 226 }
>> 227 else
>> 228 {
>> 229 epsilon = epsilonMin + epsilonRange * G4UniformRand();
>> 230 screen = screenFactor / (epsilon * (1 - epsilon));
>> 231 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
>> 232 }
>> 233 } while ( gReject < G4UniformRand() );
>> 234
>> 235 } // End of epsilon sampling
>> 236
>> 237 // Fix charges randomly
>> 238
>> 239 G4double electronTotEnergy;
>> 240 G4double positronTotEnergy;
>> 241
>> 242 if (G4int(2*G4UniformRand()))
>> 243 {
>> 244 electronTotEnergy = (1. - epsilon) * photonEnergy;
>> 245 positronTotEnergy = epsilon * photonEnergy;
>> 246 }
>> 247 else
>> 248 {
>> 249 positronTotEnergy = (1. - epsilon) * photonEnergy;
>> 250 electronTotEnergy = epsilon * photonEnergy;
>> 251 }
160 252
161 if (!fin.is_open()) { << 253 // Scattered electron (positron) angles. ( Z - axis along the parent photon)
162 G4ExceptionDescription ed; << 254 // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
163 ed << "G4LivermoreGammaConversionModel dat << 255 // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
164 << G4endl; << 256
165 G4Exception("G4LivermoreGammaConversionMod << 257 G4double u;
166 "G4LEDATA version should be G4 << 258 const G4double a1 = 0.625;
167 return; << 259 G4double a2 = 3. * a1;
168 } << 260 // G4double d = 27. ;
169 else { << 261
170 if (verboseLevel > 1) { << 262 // if (9. / (9. + d) > G4UniformRand())
171 G4cout << "File " << ost.str() << " is o << 263 if (0.25 > G4UniformRand())
>> 264 {
>> 265 u = - std::log(G4UniformRand() * G4UniformRand()) / a1 ;
>> 266 }
>> 267 else
>> 268 {
>> 269 u = - std::log(G4UniformRand() * G4UniformRand()) / a2 ;
172 } 270 }
173 271
174 data[Z]->Retrieve(fin, true); << 272 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
175 } << 273 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
176 // Activation of spline interpolation << 274 G4double phi = twopi * G4UniformRand();
177 if (useSpline) data[Z]->FillSecondDerivative << 275
>> 276 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
>> 277 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
>> 278
>> 279
>> 280 // Kinematics of the created pair:
>> 281 // the electron and positron are assumed to have a symetric angular
>> 282 // distribution with respect to the Z axis along the parent photon
>> 283
>> 284 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
>> 285
>> 286 // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
>> 287
>> 288 G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
>> 289 electronDirection.rotateUz(photonDirection);
>> 290
>> 291 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
>> 292 electronDirection,
>> 293 electronKineEnergy);
>> 294
>> 295 // The e+ is always created (even with kinetic energy = 0) for further annihilation
>> 296 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
>> 297
>> 298 // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
>> 299
>> 300 G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
>> 301 positronDirection.rotateUz(photonDirection);
>> 302
>> 303 // Create G4DynamicParticle object for the particle2
>> 304 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
>> 305 positronDirection, positronKineEnergy);
>> 306 // Fill output vector
>> 307 fvect->push_back(particle1);
>> 308 fvect->push_back(particle2);
>> 309
>> 310 // kill incident photon
>> 311 fParticleChange->SetProposedKineticEnergy(0.);
>> 312 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 313
178 } 314 }
179 315
180 //....oooOO0OOooo........oooOO0OOooo........oo 316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
181 317
182 G4double << 318 G4double G4LivermoreGammaConversionModel::ScreenFunction1(G4double screenVariable)
183 G4LivermoreGammaConversionModel::ComputeCrossS <<
184 <<
185 <<
186 { 319 {
187 if (verboseLevel > 1) { << 320 // Compute the value of the screening function 3*phi1 - phi2
188 G4cout << "G4LivermoreGammaConversionModel <<
189 } <<
190 <<
191 if (GammaEnergy < lowEnergyLimit) { <<
192 return 0.0; <<
193 } <<
194 321
195 G4double xs = 0.0; << 322 G4double value;
>> 323
>> 324 if (screenVariable > 1.)
>> 325 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
>> 326 else
>> 327 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
>> 328
>> 329 return value;
>> 330 }
196 331
197 G4int intZ = std::max(1, std::min(G4lrint(Z) << 332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198 <<
199 G4PhysicsFreeVector* pv = data[intZ]; <<
200 <<
201 // if element was not initialised <<
202 // do initialisation safely for MT mode <<
203 if (pv == nullptr) { <<
204 InitialiseForElement(particle, intZ); <<
205 pv = data[intZ]; <<
206 if (pv == nullptr) { <<
207 return xs; <<
208 } <<
209 } <<
210 // x-section is taken from the table <<
211 xs = pv->Value(GammaEnergy); <<
212 <<
213 if (verboseLevel > 0) { <<
214 G4cout << "*** Gamma conversion xs for Z=" <<
215 << " cs=" << xs / millibarn << " m <<
216 } <<
217 return xs; <<
218 } <<
219 <<
220 //....oooOO0OOooo........oooOO0OOooo........oo <<
221 333
222 void G4LivermoreGammaConversionModel::Initiali << 334 G4double G4LivermoreGammaConversionModel::ScreenFunction2(G4double screenVariable)
223 { 335 {
224 G4AutoLock l(&LivermoreGammaConversionModelM << 336 // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
225 if (data[Z] == nullptr) { << 337
226 ReadData(Z); << 338 G4double value;
227 } << 339
228 l.unlock(); << 340 if (screenVariable > 1.)
229 } << 341 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
>> 342 else
>> 343 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
>> 344
>> 345 return value;
>> 346 }
230 347
231 //....oooOO0OOooo........oooOO0OOooo........oo <<
232 348