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