| 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 // 26 //
27 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class file 29 // GEANT4 Class file
30 // 30 //
31 // 31 //
32 // File name: G4MuPairProductionModel 32 // File name: G4MuPairProductionModel
33 // 33 //
34 // Author: Vladimir Ivanchenko on base 34 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
35 // 35 //
36 // Creation date: 24.06.2002 36 // Creation date: 24.06.2002
37 // 37 //
38 // Modifications: 38 // Modifications:
39 // 39 //
40 // 04-12-02 Change G4DynamicParticle construct 40 // 04-12-02 Change G4DynamicParticle constructor in PostStep (V.Ivanchenko)
41 // 23-12-02 Change interface in order to move 41 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
42 // 24-01-03 Fix for compounds (V.Ivanchenko) 42 // 24-01-03 Fix for compounds (V.Ivanchenko)
43 // 27-01-03 Make models region aware (V.Ivanch 43 // 27-01-03 Make models region aware (V.Ivanchenko)
44 // 13-02-03 Add model (V.Ivanchenko) 44 // 13-02-03 Add model (V.Ivanchenko)
45 // 06-06-03 Fix in cross section calculation f 45 // 06-06-03 Fix in cross section calculation for high energy (V.Ivanchenko)
46 // 20-10-03 2*xi in ComputeDDMicroscopicCrossS 46 // 20-10-03 2*xi in ComputeDDMicroscopicCrossSection (R.Kokoulin)
47 // 8 integration points in ComputeDMi 47 // 8 integration points in ComputeDMicroscopicCrossSection
48 // 12-01-04 Take min cut of e- and e+ not its 48 // 12-01-04 Take min cut of e- and e+ not its sum (V.Ivanchenko)
49 // 10-02-04 Update parameterisation using R.Ko 49 // 10-02-04 Update parameterisation using R.Kokoulin model (V.Ivanchenko)
50 // 28-04-04 For complex materials repeat calcu 50 // 28-04-04 For complex materials repeat calculation of max energy for each
51 // material (V.Ivanchenko) 51 // material (V.Ivanchenko)
52 // 01-11-04 Fix bug inside ComputeDMicroscopic 52 // 01-11-04 Fix bug inside ComputeDMicroscopicCrossSection (R.Kokoulin)
53 // 08-04-05 Major optimisation of internal int 53 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
54 // 03-08-05 Add SetParticle method (V.Ivantche 54 // 03-08-05 Add SetParticle method (V.Ivantchenko)
55 // 23-10-05 Add protection in sampling of e+e- 55 // 23-10-05 Add protection in sampling of e+e- pair energy needed for
56 // low cuts (V.Ivantchenko) 56 // low cuts (V.Ivantchenko)
57 // 13-02-06 Add ComputeCrossSectionPerAtom (mm 57 // 13-02-06 Add ComputeCrossSectionPerAtom (mma)
58 // 24-04-07 Add protection in SelectRandomAtom 58 // 24-04-07 Add protection in SelectRandomAtom method (V.Ivantchenko)
59 // 12-05-06 Updated sampling (use cut) in Sele 59 // 12-05-06 Updated sampling (use cut) in SelectRandomAtom (A.Bogdanov)
60 // 11-10-07 Add ignoreCut flag (V.Ivanchenko) 60 // 11-10-07 Add ignoreCut flag (V.Ivanchenko)
61 61
62 // 62 //
63 // Class Description: 63 // Class Description:
64 // 64 //
65 // 65 //
66 // ------------------------------------------- 66 // -------------------------------------------------------------------
67 // 67 //
68 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 70
71 #include "G4MuPairProductionModel.hh" 71 #include "G4MuPairProductionModel.hh"
72 #include "G4PhysicalConstants.hh" 72 #include "G4PhysicalConstants.hh"
73 #include "G4SystemOfUnits.hh" 73 #include "G4SystemOfUnits.hh"
74 #include "G4EmParameters.hh" 74 #include "G4EmParameters.hh"
75 #include "G4Electron.hh" 75 #include "G4Electron.hh"
76 #include "G4Positron.hh" 76 #include "G4Positron.hh"
77 #include "G4MuonMinus.hh" 77 #include "G4MuonMinus.hh"
78 #include "G4MuonPlus.hh" 78 #include "G4MuonPlus.hh"
79 #include "Randomize.hh" 79 #include "Randomize.hh"
80 #include "G4Material.hh" 80 #include "G4Material.hh"
81 #include "G4Element.hh" 81 #include "G4Element.hh"
82 #include "G4ElementVector.hh" 82 #include "G4ElementVector.hh"
83 #include "G4ElementDataRegistry.hh" <<
84 #include "G4ProductionCutsTable.hh" 83 #include "G4ProductionCutsTable.hh"
85 #include "G4ParticleChangeForLoss.hh" 84 #include "G4ParticleChangeForLoss.hh"
86 #include "G4ModifiedMephi.hh" 85 #include "G4ModifiedMephi.hh"
87 #include "G4Log.hh" 86 #include "G4Log.hh"
88 #include "G4Exp.hh" 87 #include "G4Exp.hh"
89 #include "G4AutoLock.hh" <<
90 <<
91 #include <iostream> 88 #include <iostream>
92 #include <fstream> 89 #include <fstream>
93 90
94 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 92
96 const G4int G4MuPairProductionModel::ZDATPAIR[ << 93 // static members
97 << 94 //
98 const G4double G4MuPairProductionModel::xgi[] << 95 static const G4double ak1 = 6.9;
99 0.0198550717512320, 0.1016667612931865, 0. << 96 static const G4double ak2 = 1.0;
100 0.5917173212478250, 0.7627662049581645, 0. << 97 static const G4int nzdat = 5;
101 }; << 98 static const G4int zdat[5] = {1, 4, 13, 29, 92};
102 99
103 const G4double G4MuPairProductionModel::wgi[] << 100 static const G4double xgi[] =
104 0.0506142681451880, 0.1111905172266870, 0. << 101 { 0.0198550717512320, 0.1016667612931865, 0.2372337950418355, 0.4082826787521750,
105 0.1813418916891810, 0.1568533229389435, 0. << 102 0.5917173212478250, 0.7627662049581645, 0.8983332387068135, 0.9801449282487680 };
106 }; <<
107 <<
108 namespace <<
109 { <<
110 G4Mutex theMuPairMutex = G4MUTEX_INITIALIZER <<
111 103
112 const G4double ak1 = 6.9; << 104 static const G4double wgi[] =
113 const G4double ak2 = 1.0; << 105 { 0.0506142681451880, 0.1111905172266870, 0.1568533229389435, 0.1813418916891810,
114 } << 106 0.1813418916891810, 0.1568533229389435, 0.1111905172266870, 0.0506142681451880 };
115 107
116 //....oooOO0OOooo........oooOO0OOooo........oo 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 109
>> 110 using namespace std;
>> 111
118 G4MuPairProductionModel::G4MuPairProductionMod 112 G4MuPairProductionModel::G4MuPairProductionModel(const G4ParticleDefinition* p,
119 113 const G4String& nam)
120 : G4VEmModel(nam), 114 : G4VEmModel(nam),
121 factorForCross(CLHEP::fine_structure_const*C << 115 factorForCross(CLHEP::fine_structure_const*CLHEP::fine_structure_const*
122 CLHEP::classic_electr_radius*CLHEP::class << 116 CLHEP::classic_electr_radius*CLHEP::classic_electr_radius*
123 4./(3.*CLHEP::pi)), << 117 4./(3.*CLHEP::pi)),
124 sqrte(std::sqrt(G4Exp(1.))), << 118 sqrte(sqrt(G4Exp(1.))),
125 minPairEnergy(4.*CLHEP::electron_mass_c2), << 119 minPairEnergy(4.*CLHEP::electron_mass_c2),
126 lowestKinEnergy(0.85*CLHEP::GeV) << 120 lowestKinEnergy(0.85*CLHEP::GeV)
127 { 121 {
128 nist = G4NistManager::Instance(); 122 nist = G4NistManager::Instance();
129 123
130 theElectron = G4Electron::Electron(); 124 theElectron = G4Electron::Electron();
131 thePositron = G4Positron::Positron(); 125 thePositron = G4Positron::Positron();
132 126
133 if(nullptr != p) { 127 if(nullptr != p) {
134 SetParticle(p); 128 SetParticle(p);
135 lowestKinEnergy = std::max(lowestKinEnergy 129 lowestKinEnergy = std::max(lowestKinEnergy, p->GetPDGMass()*8.0);
136 } 130 }
137 emin = lowestKinEnergy; 131 emin = lowestKinEnergy;
138 emax = emin*10000.; 132 emax = emin*10000.;
139 SetAngularDistribution(new G4ModifiedMephi() 133 SetAngularDistribution(new G4ModifiedMephi());
140 } 134 }
141 135
142 //....oooOO0OOooo........oooOO0OOooo........oo 136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
143 137
144 G4double G4MuPairProductionModel::MinPrimaryEn 138 G4double G4MuPairProductionModel::MinPrimaryEnergy(const G4Material*,
145 139 const G4ParticleDefinition*,
146 140 G4double cut)
147 { 141 {
148 return std::max(lowestKinEnergy, cut); 142 return std::max(lowestKinEnergy, cut);
149 } 143 }
150 144
151 //....oooOO0OOooo........oooOO0OOooo........oo 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152 146
153 void G4MuPairProductionModel::Initialise(const 147 void G4MuPairProductionModel::Initialise(const G4ParticleDefinition* p,
154 const 148 const G4DataVector& cuts)
155 { 149 {
156 SetParticle(p); 150 SetParticle(p);
157 151
158 if (nullptr == fParticleChange) { << 152 if(nullptr == fParticleChange) {
159 fParticleChange = GetParticleChangeForLoss 153 fParticleChange = GetParticleChangeForLoss();
160 <<
161 // define scale of internal table for each <<
162 if (0 == nbine) { <<
163 emin = std::max(lowestKinEnergy, LowEner <<
164 emax = std::max(HighEnergyLimit(), emin* <<
165 nbine = std::size_t(nYBinPerDecade*std:: <<
166 if(nbine < 3) { nbine = 3; } <<
167 <<
168 ymin = G4Log(minPairEnergy/emin); <<
169 dy = -ymin/G4double(nbiny); <<
170 } <<
171 if (p == particle) { <<
172 G4int pdg = std::abs(p->GetPDGEncoding() <<
173 if (pdg == 2212) { <<
174 dataName = "pEEPairProd"; <<
175 } else if (pdg == 321) { <<
176 dataName = "kaonEEPairProd"; <<
177 } else if (pdg == 211) { <<
178 dataName = "pionEEPairProd"; <<
179 } else if (pdg == 11) { <<
180 dataName = "eEEPairProd"; <<
181 } else if (pdg == 13) { <<
182 if (GetName() == "muToMuonPairProd") { <<
183 dataName = "muMuMuPairProd"; <<
184 } else { <<
185 dataName = "muEEPairProd"; <<
186 } <<
187 } <<
188 } <<
189 } 154 }
190 155
191 // for low-energy application this process s 156 // for low-energy application this process should not work
192 if(lowestKinEnergy >= HighEnergyLimit()) { r 157 if(lowestKinEnergy >= HighEnergyLimit()) { return; }
193 158
194 if (p == particle) { << 159 // define scale of internal table for each thread only once
195 fElementData = << 160 if(0 == nbine) {
196 G4ElementDataRegistry::Instance()->GetEl << 161 emin = std::max(lowestKinEnergy, LowEnergyLimit());
197 if (nullptr == fElementData) { << 162 emax = std::max(HighEnergyLimit(), emin*2);
198 G4AutoLock l(&theMuPairMutex); << 163 nbine = size_t(nYBinPerDecade*std::log10(emax/emin));
199 fElementData = << 164 if(nbine < 3) { nbine = 3; }
200 G4ElementDataRegistry::Instance()->GetElemen << 165
201 if (nullptr == fElementData) { << 166 ymin = G4Log(minPairEnergy/emin);
202 fElementData = new G4ElementData(NZDATPAIR); << 167 dy = -ymin/G4double(nbiny);
203 fElementData->SetName(dataName); << 168 }
204 } << 169
205 G4bool useDataFile = G4EmParameters::Ins << 170 if(IsMaster() && p == particle) {
206 if (useDataFile) { useDataFile = Retrie << 171 if(nullptr == fElementData) {
207 if (!useDataFile) { MakeSamplingTables() << 172 fElementData = new G4ElementData();
208 if (fTableToFile) { StoreTables(); } << 173 G4bool dataFile = G4EmParameters::Instance()->RetrieveMuDataFromFile();
209 l.unlock(); << 174 if(dataFile) { dataFile = RetrieveTables(); }
210 } << 175 if(!dataFile) { MakeSamplingTables(); }
211 if (IsMaster()) { << 176 if(fTableToFile) { StoreTables(); }
212 InitialiseElementSelectors(p, cuts); << 177 }
213 } << 178 InitialiseElementSelectors(p, cuts);
214 } 179 }
215 } 180 }
216 181
217 //....oooOO0OOooo........oooOO0OOooo........oo 182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
218 183
219 void G4MuPairProductionModel::InitialiseLocal( 184 void G4MuPairProductionModel::InitialiseLocal(const G4ParticleDefinition* p,
220 185 G4VEmModel* masterModel)
221 { 186 {
222 if(p == particle && lowestKinEnergy < HighEn 187 if(p == particle && lowestKinEnergy < HighEnergyLimit()) {
223 SetElementSelectors(masterModel->GetElemen 188 SetElementSelectors(masterModel->GetElementSelectors());
>> 189 fElementData = masterModel->GetElementData();
224 } 190 }
225 } 191 }
226 192
227 //....oooOO0OOooo........oooOO0OOooo........oo 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
228 194
229 G4double G4MuPairProductionModel::ComputeDEDXP 195 G4double G4MuPairProductionModel::ComputeDEDXPerVolume(
230 196 const G4Material* material,
231 197 const G4ParticleDefinition*,
232 198 G4double kineticEnergy,
233 199 G4double cutEnergy)
234 { 200 {
235 G4double dedx = 0.0; 201 G4double dedx = 0.0;
236 if (cutEnergy <= minPairEnergy || kineticEne 202 if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
237 { return dedx; } 203 { return dedx; }
238 204
239 const G4ElementVector* theElementVector = ma 205 const G4ElementVector* theElementVector = material->GetElementVector();
240 const G4double* theAtomicNumDensityVector = 206 const G4double* theAtomicNumDensityVector =
241 material->G 207 material->GetAtomicNumDensityVector();
242 208
243 // loop for elements in the material 209 // loop for elements in the material
244 for (std::size_t i=0; i<material->GetNumberO << 210 for (size_t i=0; i<material->GetNumberOfElements(); ++i) {
245 G4double Z = (*theElementVector)[i]->GetZ 211 G4double Z = (*theElementVector)[i]->GetZ();
246 G4double tmax = MaxSecondaryEnergyForElem 212 G4double tmax = MaxSecondaryEnergyForElement(kineticEnergy, Z);
247 G4double loss = ComputMuPairLoss(Z, kinet 213 G4double loss = ComputMuPairLoss(Z, kineticEnergy, cutEnergy, tmax);
248 dedx += loss*theAtomicNumDensityVector[i] 214 dedx += loss*theAtomicNumDensityVector[i];
249 } 215 }
250 dedx = std::max(dedx, 0.0); 216 dedx = std::max(dedx, 0.0);
251 return dedx; 217 return dedx;
252 } 218 }
253 219
254 //....oooOO0OOooo........oooOO0OOooo........oo 220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255 221
256 G4double G4MuPairProductionModel::ComputMuPair 222 G4double G4MuPairProductionModel::ComputMuPairLoss(G4double Z,
257 223 G4double tkin,
258 224 G4double cutEnergy,
259 225 G4double tmax)
260 { 226 {
261 G4double loss = 0.0; 227 G4double loss = 0.0;
262 228
263 G4double cut = std::min(cutEnergy, tmax); 229 G4double cut = std::min(cutEnergy, tmax);
264 if(cut <= minPairEnergy) { return loss; } 230 if(cut <= minPairEnergy) { return loss; }
265 231
266 // calculate the rectricted loss 232 // calculate the rectricted loss
267 // numerical integration in log(PairEnergy) 233 // numerical integration in log(PairEnergy)
268 G4double aaa = G4Log(minPairEnergy); 234 G4double aaa = G4Log(minPairEnergy);
269 G4double bbb = G4Log(cut); 235 G4double bbb = G4Log(cut);
270 236
271 G4int kkk = std::min(std::max(G4lrint((bbb-a << 237 G4int kkk = G4lrint((bbb-aaa)/ak1+ak2);
>> 238 if(kkk > 8) { kkk = 8; }
>> 239 else if (kkk < 1) { kkk = 1; }
272 G4double hhh = (bbb-aaa)/kkk; 240 G4double hhh = (bbb-aaa)/kkk;
273 G4double x = aaa; 241 G4double x = aaa;
274 242
275 for (G4int l=0 ; l<kkk; ++l) { 243 for (G4int l=0 ; l<kkk; ++l) {
276 for (G4int ll=0; ll<NINTPAIR; ++ll) { << 244 for (G4int ll=0; ll<8; ++ll) {
277 G4double ep = G4Exp(x+xgi[ll]*hhh); 245 G4double ep = G4Exp(x+xgi[ll]*hhh);
278 loss += wgi[ll]*ep*ep*ComputeDMicroscopi 246 loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
279 } 247 }
280 x += hhh; 248 x += hhh;
281 } 249 }
282 loss *= hhh; 250 loss *= hhh;
283 loss = std::max(loss, 0.0); 251 loss = std::max(loss, 0.0);
284 return loss; 252 return loss;
285 } 253 }
286 254
287 //....oooOO0OOooo........oooOO0OOooo........oo 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 256
289 G4double G4MuPairProductionModel::ComputeMicro 257 G4double G4MuPairProductionModel::ComputeMicroscopicCrossSection(
290 G4d 258 G4double tkin,
291 G4d 259 G4double Z,
292 G4d 260 G4double cutEnergy)
293 { 261 {
294 G4double cross = 0.; 262 G4double cross = 0.;
295 G4double tmax = MaxSecondaryEnergyForElement 263 G4double tmax = MaxSecondaryEnergyForElement(tkin, Z);
296 G4double cut = std::max(cutEnergy, minPairE 264 G4double cut = std::max(cutEnergy, minPairEnergy);
297 if (tmax <= cut) { return cross; } 265 if (tmax <= cut) { return cross; }
298 266
299 G4double aaa = G4Log(cut); 267 G4double aaa = G4Log(cut);
300 G4double bbb = G4Log(tmax); 268 G4double bbb = G4Log(tmax);
301 G4int kkk = std::min(std::max(G4lrint((bbb-a << 269 G4int kkk = G4lrint((bbb-aaa)/ak1 + ak2);
>> 270 if(kkk > 8) { kkk = 8; }
>> 271 else if (kkk < 1) { kkk = 1; }
302 272
303 G4double hhh = (bbb-aaa)/(kkk); 273 G4double hhh = (bbb-aaa)/(kkk);
304 G4double x = aaa; 274 G4double x = aaa;
305 275
306 for (G4int l=0; l<kkk; ++l) { << 276 for(G4int l=0; l<kkk; ++l) {
307 for (G4int i=0; i<NINTPAIR; ++i) { << 277 for(G4int i=0; i<8; ++i) {
308 G4double ep = G4Exp(x + xgi[i]*hhh); 278 G4double ep = G4Exp(x + xgi[i]*hhh);
309 cross += ep*wgi[i]*ComputeDMicroscopicCr 279 cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
310 } 280 }
311 x += hhh; 281 x += hhh;
312 } 282 }
313 283
314 cross *= hhh; 284 cross *= hhh;
315 cross = std::max(cross, 0.0); 285 cross = std::max(cross, 0.0);
316 return cross; 286 return cross;
317 } 287 }
318 288
319 //....oooOO0OOooo........oooOO0OOooo........oo 289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
320 290
321 G4double G4MuPairProductionModel::ComputeDMicr 291 G4double G4MuPairProductionModel::ComputeDMicroscopicCrossSection(
322 G4d 292 G4double tkin,
323 G4d 293 G4double Z,
324 G4d 294 G4double pairEnergy)
325 // Calculates the differential (D) microscopi 295 // Calculates the differential (D) microscopic cross section
326 // using the cross section formula of R.P. Kok 296 // using the cross section formula of R.P. Kokoulin (18/01/98)
327 // Code modified by R.P. Kokoulin, V.N. Ivanch 297 // Code modified by R.P. Kokoulin, V.N. Ivanchenko (27/01/04)
328 { 298 {
329 static const G4double bbbtf= 183. ; 299 static const G4double bbbtf= 183. ;
330 static const G4double bbbh = 202.4 ; 300 static const G4double bbbh = 202.4 ;
331 static const G4double g1tf = 1.95e-5 ; 301 static const G4double g1tf = 1.95e-5 ;
332 static const G4double g2tf = 5.3e-5 ; 302 static const G4double g2tf = 5.3e-5 ;
333 static const G4double g1h = 4.4e-5 ; 303 static const G4double g1h = 4.4e-5 ;
334 static const G4double g2h = 4.8e-5 ; 304 static const G4double g2h = 4.8e-5 ;
335 305
336 if (pairEnergy <= minPairEnergy) 306 if (pairEnergy <= minPairEnergy)
337 return 0.0; 307 return 0.0;
338 308
339 G4double totalEnergy = tkin + particleMass; 309 G4double totalEnergy = tkin + particleMass;
340 G4double residEnergy = totalEnergy - pairEn 310 G4double residEnergy = totalEnergy - pairEnergy;
341 311
342 if (residEnergy <= 0.75*sqrte*z13*particleMa 312 if (residEnergy <= 0.75*sqrte*z13*particleMass)
343 return 0.0; 313 return 0.0;
344 314
345 G4double a0 = 1.0 / (totalEnergy * residEner 315 G4double a0 = 1.0 / (totalEnergy * residEnergy);
346 G4double alf = 4.0 * electron_mass_c2 / pair 316 G4double alf = 4.0 * electron_mass_c2 / pairEnergy;
347 G4double rt = std::sqrt(1.0 - alf); << 317 G4double rt = sqrt(1.0 - alf);
348 G4double delta = 6.0 * particleMass * partic 318 G4double delta = 6.0 * particleMass * particleMass * a0;
349 G4double tmnexp = alf/(1.0 + rt) + delta*rt; 319 G4double tmnexp = alf/(1.0 + rt) + delta*rt;
350 320
351 if(tmnexp >= 1.0) { return 0.0; } 321 if(tmnexp >= 1.0) { return 0.0; }
352 322
353 G4double tmn = G4Log(tmnexp); 323 G4double tmn = G4Log(tmnexp);
354 324
355 G4double massratio = particleMass/CLHEP::ele << 325 G4double massratio = particleMass/electron_mass_c2;
356 G4double massratio2 = massratio*massratio; << 326 G4double massratio2 = massratio*massratio;
357 G4double inv_massratio2 = 1.0 / massratio2; 327 G4double inv_massratio2 = 1.0 / massratio2;
358 328
359 // zeta calculation 329 // zeta calculation
360 G4double bbb,g1,g2; 330 G4double bbb,g1,g2;
361 if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = 331 if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
362 else { bbb = bbbtf; g1 = g1tf; g2 = 332 else { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
363 333
364 G4double zeta = 0.0; 334 G4double zeta = 0.0;
365 G4double z1exp = totalEnergy / (particleMass 335 G4double z1exp = totalEnergy / (particleMass + g1*z23*totalEnergy);
366 336
367 // 35.221047195922 is the root of zeta1(x) = 337 // 35.221047195922 is the root of zeta1(x) = 0.073 * log(x) - 0.26, so the
368 // condition below is the same as zeta1 > 0. 338 // condition below is the same as zeta1 > 0.0, but without calling log(x)
369 if (z1exp > 35.221047195922) 339 if (z1exp > 35.221047195922)
370 { 340 {
371 G4double z2exp = totalEnergy / (particleMa 341 G4double z2exp = totalEnergy / (particleMass + g2*z13*totalEnergy);
372 zeta = (0.073 * G4Log(z1exp) - 0.26) / (0. 342 zeta = (0.073 * G4Log(z1exp) - 0.26) / (0.058 * G4Log(z2exp) - 0.14);
373 } 343 }
374 344
375 G4double z2 = Z*(Z+zeta); 345 G4double z2 = Z*(Z+zeta);
376 G4double screen0 = 2.*electron_mass_c2*sqrte 346 G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
377 G4double beta = 0.5*pairEnergy*pairEnergy*a0 347 G4double beta = 0.5*pairEnergy*pairEnergy*a0;
378 G4double xi0 = 0.5*massratio2*beta; 348 G4double xi0 = 0.5*massratio2*beta;
379 349
380 // Gaussian integration in ln(1-ro) ( with 8 350 // Gaussian integration in ln(1-ro) ( with 8 points)
381 G4double rho[NINTPAIR]; << 351 G4double rho[8];
382 G4double rho2[NINTPAIR]; << 352 G4double rho2[8];
383 G4double xi[NINTPAIR]; << 353 G4double xi[8];
384 G4double xi1[NINTPAIR]; << 354 G4double xi1[8];
385 G4double xii[NINTPAIR]; << 355 G4double xii[8];
386 356
387 for (G4int i = 0; i < NINTPAIR; ++i) << 357 for (G4int i = 0; i < 8; ++i)
388 { 358 {
389 rho[i] = G4Exp(tmn*xgi[i]) - 1.0; // rho = 359 rho[i] = G4Exp(tmn*xgi[i]) - 1.0; // rho = -asymmetry
390 rho2[i] = rho[i] * rho[i]; 360 rho2[i] = rho[i] * rho[i];
391 xi[i] = xi0*(1.0-rho2[i]); 361 xi[i] = xi0*(1.0-rho2[i]);
392 xi1[i] = 1.0 + xi[i]; 362 xi1[i] = 1.0 + xi[i];
393 xii[i] = 1.0 / xi[i]; 363 xii[i] = 1.0 / xi[i];
394 } 364 }
395 365
396 G4double ye1[NINTPAIR]; << 366 G4double ye1[8];
397 G4double ym1[NINTPAIR]; << 367 G4double ym1[8];
398 368
399 G4double b40 = 4.0 * beta; 369 G4double b40 = 4.0 * beta;
400 G4double b62 = 6.0 * beta + 2.0; 370 G4double b62 = 6.0 * beta + 2.0;
401 371
402 for (G4int i = 0; i < NINTPAIR; ++i) << 372 for (G4int i = 0; i < 8; ++i)
403 { 373 {
404 G4double yeu = (b40 + 5.0) + (b40 - 1.0) * 374 G4double yeu = (b40 + 5.0) + (b40 - 1.0) * rho2[i];
405 G4double yed = b62*G4Log(3.0 + xii[i]) + ( 375 G4double yed = b62*G4Log(3.0 + xii[i]) + (2.0 * beta - 1.0)*rho2[i] - b40;
406 376
407 G4double ymu = b62 * (1.0 + rho2[i]) + 6.0 377 G4double ymu = b62 * (1.0 + rho2[i]) + 6.0;
408 G4double ymd = (b40 + 3.0)*(1.0 + rho2[i]) 378 G4double ymd = (b40 + 3.0)*(1.0 + rho2[i])*G4Log(3.0 + xi[i])
409 + 2.0 - 3.0 * rho2[i]; 379 + 2.0 - 3.0 * rho2[i];
410 380
411 ye1[i] = 1.0 + yeu / yed; 381 ye1[i] = 1.0 + yeu / yed;
412 ym1[i] = 1.0 + ymu / ymd; 382 ym1[i] = 1.0 + ymu / ymd;
413 } 383 }
414 384
415 G4double be[NINTPAIR]; << 385 G4double be[8];
416 G4double bm[NINTPAIR]; << 386 G4double bm[8];
417 387
418 for(G4int i = 0; i < NINTPAIR; ++i) { << 388 for(G4int i = 0; i < 8; ++i) {
419 if(xi[i] <= 1000.0) { 389 if(xi[i] <= 1000.0) {
420 be[i] = ((2.0 + rho2[i])*(1.0 + beta) + 390 be[i] = ((2.0 + rho2[i])*(1.0 + beta) +
421 xi[i]*(3.0 + rho2[i]))*G4Log(1.0 + xi 391 xi[i]*(3.0 + rho2[i]))*G4Log(1.0 + xii[i]) +
422 (1.0 - rho2[i] - beta)/xi1[i] - (3.0 + rho2[ 392 (1.0 - rho2[i] - beta)/xi1[i] - (3.0 + rho2[i]);
423 } else { 393 } else {
424 be[i] = 0.5*(3.0 - rho2[i] + 2.0*beta*(1 394 be[i] = 0.5*(3.0 - rho2[i] + 2.0*beta*(1.0 + rho2[i]))*xii[i];
425 } 395 }
426 396
427 if(xi[i] >= 0.001) { 397 if(xi[i] >= 0.001) {
428 G4double a10 = (1.0 + 2.0 * beta) * (1.0 398 G4double a10 = (1.0 + 2.0 * beta) * (1.0 - rho2[i]);
429 bm[i] = ((1.0 + rho2[i])*(1.0 + 1.5 * be 399 bm[i] = ((1.0 + rho2[i])*(1.0 + 1.5 * beta) - a10*xii[i])*G4Log(xi1[i]) +
430 xi[i] * (1.0 - rho2[i] - beta) 400 xi[i] * (1.0 - rho2[i] - beta)/xi1[i] + a10;
431 } else { 401 } else {
432 bm[i] = 0.5*(5.0 - rho2[i] + beta * (3.0 402 bm[i] = 0.5*(5.0 - rho2[i] + beta * (3.0 + rho2[i]))*xi[i];
433 } 403 }
434 } 404 }
435 405
436 G4double sum = 0.0; 406 G4double sum = 0.0;
437 407
438 for (G4int i = 0; i < NINTPAIR; ++i) { << 408 for (G4int i = 0; i < 8; ++i) {
439 G4double screen = screen0*xi1[i]/(1.0 - rh 409 G4double screen = screen0*xi1[i]/(1.0 - rho2[i]);
440 G4double ale = G4Log(bbb/z13*std::sqrt(xi1 << 410 G4double ale = G4Log(bbb/z13*sqrt(xi1[i]*ye1[i])/(1. + screen*ye1[i]));
441 G4double cre = 0.5*G4Log(1. + 2.25*z23*xi1 411 G4double cre = 0.5*G4Log(1. + 2.25*z23*xi1[i]*ye1[i]*inv_massratio2);
442 412
443 G4double fe = (ale-cre)*be[i]; 413 G4double fe = (ale-cre)*be[i];
444 fe = std::max(fe, 0.0); 414 fe = std::max(fe, 0.0);
445 415
446 G4double alm_crm = G4Log(bbb*massratio/(1. 416 G4double alm_crm = G4Log(bbb*massratio/(1.5*z23*(1. + screen*ym1[i])));
447 G4double fm = std::max(alm_crm*bm[i], 0.0) 417 G4double fm = std::max(alm_crm*bm[i], 0.0)*inv_massratio2;
448 418
449 sum += wgi[i]*(1.0 + rho[i])*(fe + fm); 419 sum += wgi[i]*(1.0 + rho[i])*(fe + fm);
450 } 420 }
451 421
452 return -tmn*sum*factorForCross*z2*residEnerg 422 return -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
453 } 423 }
454 424
455 //....oooOO0OOooo........oooOO0OOooo........oo 425 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
456 426
457 G4double G4MuPairProductionModel::ComputeCross 427 G4double G4MuPairProductionModel::ComputeCrossSectionPerAtom(
458 con 428 const G4ParticleDefinition*,
459 429 G4double kineticEnergy,
460 430 G4double Z, G4double,
461 431 G4double cutEnergy,
462 432 G4double maxEnergy)
463 { 433 {
464 G4double cross = 0.0; 434 G4double cross = 0.0;
465 if (kineticEnergy <= lowestKinEnergy) { retu 435 if (kineticEnergy <= lowestKinEnergy) { return cross; }
466 436
467 G4double maxPairEnergy = MaxSecondaryEnergyF 437 G4double maxPairEnergy = MaxSecondaryEnergyForElement(kineticEnergy, Z);
468 G4double tmax = std::min(maxEnergy, maxPairE 438 G4double tmax = std::min(maxEnergy, maxPairEnergy);
469 G4double cut = std::max(cutEnergy, minPairE 439 G4double cut = std::max(cutEnergy, minPairEnergy);
470 if (cut >= tmax) { return cross; } 440 if (cut >= tmax) { return cross; }
471 441
472 cross = ComputeMicroscopicCrossSection(kinet 442 cross = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut);
473 if(tmax < kineticEnergy) { 443 if(tmax < kineticEnergy) {
474 cross -= ComputeMicroscopicCrossSection(ki 444 cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
475 } 445 }
476 return cross; 446 return cross;
477 } 447 }
478 448
479 //....oooOO0OOooo........oooOO0OOooo........oo 449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
480 450
481 void G4MuPairProductionModel::MakeSamplingTabl 451 void G4MuPairProductionModel::MakeSamplingTables()
482 { 452 {
483 G4double factore = G4Exp(G4Log(emax/emin)/G4 453 G4double factore = G4Exp(G4Log(emax/emin)/G4double(nbine));
484 454
485 for (G4int iz=0; iz<NZDATPAIR; ++iz) { << 455 for (G4int iz=0; iz<nzdat; ++iz) {
486 456
487 G4double Z = ZDATPAIR[iz]; << 457 G4double Z = zdat[iz];
488 G4Physics2DVector* pv = new G4Physics2DVec 458 G4Physics2DVector* pv = new G4Physics2DVector(nbiny+1,nbine+1);
489 G4double kinEnergy = emin; 459 G4double kinEnergy = emin;
490 460
491 for (std::size_t it=0; it<=nbine; ++it) { << 461 for (size_t it=0; it<=nbine; ++it) {
492 462
493 pv->PutY(it, G4Log(kinEnergy/CLHEP::MeV) 463 pv->PutY(it, G4Log(kinEnergy/CLHEP::MeV));
494 G4double maxPairEnergy = MaxSecondaryEne 464 G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy, Z);
495 /* 465 /*
496 G4cout << "it= " << it << " E= " << kinE 466 G4cout << "it= " << it << " E= " << kinEnergy
497 << " " << particle->GetParticleN 467 << " " << particle->GetParticleName()
498 << " maxE= " << maxPairEnergy << 468 << " maxE= " << maxPairEnergy << " minE= " << minPairEnergy
499 << " ymin= " << ymin << G4endl; 469 << " ymin= " << ymin << G4endl;
500 */ 470 */
501 G4double coef = G4Log(minPairEnergy/kinE 471 G4double coef = G4Log(minPairEnergy/kinEnergy)/ymin;
502 G4double ymax = G4Log(maxPairEnergy/kinE 472 G4double ymax = G4Log(maxPairEnergy/kinEnergy)/coef;
503 G4double fac = (ymax - ymin)/dy; 473 G4double fac = (ymax - ymin)/dy;
504 std::size_t imax = (std::size_t)fac; << 474 size_t imax = (size_t)fac;
505 fac -= (G4double)imax; 475 fac -= (G4double)imax;
506 476
507 G4double xSec = 0.0; 477 G4double xSec = 0.0;
508 G4double x = ymin; 478 G4double x = ymin;
509 /* 479 /*
510 G4cout << "Z= " << currentZ << " z13= " 480 G4cout << "Z= " << currentZ << " z13= " << z13
511 << " mE= " << maxPairEnergy << " 481 << " mE= " << maxPairEnergy << " ymin= " << ymin
512 << " dy= " << dy << " c= " << co 482 << " dy= " << dy << " c= " << coef << G4endl;
513 */ 483 */
514 // start from zero 484 // start from zero
515 pv->PutValue(0, it, 0.0); 485 pv->PutValue(0, it, 0.0);
516 if(0 == it) { pv->PutX(nbiny, 0.0); } 486 if(0 == it) { pv->PutX(nbiny, 0.0); }
517 487
518 for (std::size_t i=0; i<nbiny; ++i) { << 488 for (size_t i=0; i<nbiny; ++i) {
519 489
520 if(0 == it) { pv->PutX(i, x); } 490 if(0 == it) { pv->PutX(i, x); }
521 491
522 if(i < imax) { 492 if(i < imax) {
523 G4double ep = kinEnergy*G4Exp(coef*( 493 G4double ep = kinEnergy*G4Exp(coef*(x + dy*0.5));
524 494
525 // not multiplied by interval, becau 495 // not multiplied by interval, because table
526 // will be used only for sampling 496 // will be used only for sampling
527 //G4cout << "i= " << i << " x= " << 497 //G4cout << "i= " << i << " x= " << x << "E= " << kinEnergy
528 // << " Egamma= " << ep << G 498 // << " Egamma= " << ep << G4endl;
529 xSec += ep*ComputeDMicroscopicCrossS 499 xSec += ep*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
530 500
531 // last bin before the kinematic lim 501 // last bin before the kinematic limit
532 } else if(i == imax) { 502 } else if(i == imax) {
533 G4double ep = kinEnergy*G4Exp(coef*( 503 G4double ep = kinEnergy*G4Exp(coef*(x + fac*dy*0.5));
534 xSec += ep*fac*ComputeDMicroscopicCr 504 xSec += ep*fac*ComputeDMicroscopicCrossSection(kinEnergy, Z, ep);
535 } 505 }
536 pv->PutValue(i + 1, it, xSec); 506 pv->PutValue(i + 1, it, xSec);
537 x += dy; 507 x += dy;
538 } 508 }
539 kinEnergy *= factore; 509 kinEnergy *= factore;
540 510
541 // to avoid precision lost 511 // to avoid precision lost
542 if(it+1 == nbine) { kinEnergy = emax; } 512 if(it+1 == nbine) { kinEnergy = emax; }
543 } 513 }
544 fElementData->InitialiseForElement(iz, pv) << 514 fElementData->InitialiseForElement(zdat[iz], pv);
545 } 515 }
546 } 516 }
547 517
548 //....oooOO0OOooo........oooOO0OOooo........oo 518 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
549 519
550 void G4MuPairProductionModel::SampleSecondarie 520 void G4MuPairProductionModel::SampleSecondaries(
551 std::vector<G4Dy 521 std::vector<G4DynamicParticle*>* vdp,
552 const G4Material 522 const G4MaterialCutsCouple* couple,
553 const G4DynamicP 523 const G4DynamicParticle* aDynamicParticle,
554 G4double tmin, 524 G4double tmin,
555 G4double tmax) 525 G4double tmax)
556 { 526 {
557 G4double kinEnergy = aDynamicParticle->GetKi 527 G4double kinEnergy = aDynamicParticle->GetKineticEnergy();
558 //G4cout << "------- G4MuPairProductionModel 528 //G4cout << "------- G4MuPairProductionModel::SampleSecondaries E(MeV)= "
559 // << kinEnergy << " " 529 // << kinEnergy << " "
560 // << aDynamicParticle->GetDefinitio 530 // << aDynamicParticle->GetDefinition()->GetParticleName() << G4endl;
561 G4double totalEnergy = kinEnergy + particl 531 G4double totalEnergy = kinEnergy + particleMass;
562 G4double totalMomentum = 532 G4double totalMomentum =
563 std::sqrt(kinEnergy*(kinEnergy + 2.0*parti << 533 sqrt(kinEnergy*(kinEnergy + 2.0*particleMass));
564 534
565 G4ThreeVector partDirection = aDynamicPartic 535 G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
566 536
567 // select randomly one element constituing t 537 // select randomly one element constituing the material
568 const G4Element* anElement = SelectRandomAto 538 const G4Element* anElement = SelectRandomAtom(couple,particle,kinEnergy);
569 539
570 // define interval of energy transfer 540 // define interval of energy transfer
571 G4double maxPairEnergy = MaxSecondaryEnergyF 541 G4double maxPairEnergy = MaxSecondaryEnergyForElement(kinEnergy,
572 542 anElement->GetZ());
573 G4double maxEnergy = std::min(tmax, maxPairE 543 G4double maxEnergy = std::min(tmax, maxPairEnergy);
574 G4double minEnergy = std::max(tmin, minPairE 544 G4double minEnergy = std::max(tmin, minPairEnergy);
575 545
576 if (minEnergy >= maxEnergy) { return; } << 546 if(minEnergy >= maxEnergy) { return; }
577 //G4cout << "emin= " << minEnergy << " emax= 547 //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy
578 // << " minPair= " << minPairEnergy << " max 548 // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy
579 // << " ymin= " << ymin << " dy= " << dy 549 // << " ymin= " << ymin << " dy= " << dy << G4endl;
580 550
581 G4double coeff = G4Log(minPairEnergy/kinEner 551 G4double coeff = G4Log(minPairEnergy/kinEnergy)/ymin;
582 552
583 // compute limits 553 // compute limits
584 G4double yymin = G4Log(minEnergy/kinEnergy)/ 554 G4double yymin = G4Log(minEnergy/kinEnergy)/coeff;
585 G4double yymax = G4Log(maxEnergy/kinEnergy)/ 555 G4double yymax = G4Log(maxEnergy/kinEnergy)/coeff;
586 556
587 //G4cout << "yymin= " << yymin << " yymax= 557 //G4cout << "yymin= " << yymin << " yymax= " << yymax << G4endl;
588 558
589 // units should not be used, bacause table w 559 // units should not be used, bacause table was built without
590 G4double logTkin = G4Log(kinEnergy/CLHEP::Me 560 G4double logTkin = G4Log(kinEnergy/CLHEP::MeV);
591 561
592 // sample e-e+ energy, pair energy first 562 // sample e-e+ energy, pair energy first
593 563
594 // select sample table via Z 564 // select sample table via Z
595 G4int iz1(0), iz2(0); 565 G4int iz1(0), iz2(0);
596 for (G4int iz=0; iz<NZDATPAIR; ++iz) { << 566 for(G4int iz=0; iz<nzdat; ++iz) {
597 if(currentZ == ZDATPAIR[iz]) { << 567 if(currentZ == zdat[iz]) {
598 iz1 = iz2 = iz; << 568 iz1 = iz2 = currentZ;
599 break; 569 break;
600 } else if(currentZ < ZDATPAIR[iz]) { << 570 } else if(currentZ < zdat[iz]) {
601 iz2 = iz; << 571 iz2 = zdat[iz];
602 if(iz > 0) { iz1 = iz-1; } << 572 if(iz > 0) { iz1 = zdat[iz-1]; }
603 else { iz1 = iz2; } 573 else { iz1 = iz2; }
604 break; 574 break;
605 } 575 }
606 } 576 }
607 if (0 == iz1) { iz1 = iz2 = NZDATPAIR-1; } << 577 if(0 == iz1) { iz1 = iz2 = zdat[nzdat-1]; }
608 578
609 G4double pairEnergy = 0.0; 579 G4double pairEnergy = 0.0;
610 G4int count = 0; 580 G4int count = 0;
611 //G4cout << "start loop Z1= " << iz1 << " Z2 581 //G4cout << "start loop Z1= " << iz1 << " Z2= " << iz2 << G4endl;
612 do { 582 do {
613 ++count; 583 ++count;
614 // sampling using only one random number 584 // sampling using only one random number
615 G4double rand = G4UniformRand(); 585 G4double rand = G4UniformRand();
616 586
617 G4double x = FindScaledEnergy(iz1, rand, l 587 G4double x = FindScaledEnergy(iz1, rand, logTkin, yymin, yymax);
618 if(iz1 != iz2) { 588 if(iz1 != iz2) {
619 G4double x2 = FindScaledEnergy(iz2, rand 589 G4double x2 = FindScaledEnergy(iz2, rand, logTkin, yymin, yymax);
620 G4double lz1= nist->GetLOGZ(ZDATPAIR[iz1 << 590 G4double lz1= nist->GetLOGZ(iz1);
621 G4double lz2= nist->GetLOGZ(ZDATPAIR[iz2 << 591 G4double lz2= nist->GetLOGZ(iz2);
622 //G4cout << count << ". x= " << x << " 592 //G4cout << count << ". x= " << x << " x2= " << x2
623 // << " Z1= " << iz1 << " Z2 593 // << " Z1= " << iz1 << " Z2= " << iz2 << G4endl;
624 x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1); 594 x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);
625 } 595 }
626 //G4cout << "x= " << x << " coeff= " << c 596 //G4cout << "x= " << x << " coeff= " << coeff << G4endl;
627 pairEnergy = kinEnergy*G4Exp(x*coeff); 597 pairEnergy = kinEnergy*G4Exp(x*coeff);
628 598
629 // Loop checking, 03-Aug-2015, Vladimir Iv 599 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
630 } while((pairEnergy < minEnergy || pairEnerg 600 } while((pairEnergy < minEnergy || pairEnergy > maxEnergy) && 10 > count);
631 601
632 //G4cout << "## pairEnergy(GeV)= " << pairEn 602 //G4cout << "## pairEnergy(GeV)= " << pairEnergy/GeV
633 // << " Etot(GeV)= " << totalEnergy/ 603 // << " Etot(GeV)= " << totalEnergy/GeV << G4endl;
634 604
635 // sample r=(E+-E-)/pairEnergy ( uniformly 605 // sample r=(E+-E-)/pairEnergy ( uniformly .....)
636 G4double rmax = 606 G4double rmax =
637 (1.-6.*particleMass*particleMass/(totalEne 607 (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-pairEnergy)))
638 *std::sqrt(1.-minPairEnergy/pairEnergy); << 608 *sqrt(1.-minPairEnergy/pairEnergy);
639 G4double r = rmax * (-1.+2.*G4UniformRand()) 609 G4double r = rmax * (-1.+2.*G4UniformRand()) ;
640 610
641 // compute energies from pairEnergy,r 611 // compute energies from pairEnergy,r
642 G4double eEnergy = (1.-r)*pairEnergy*0.5; 612 G4double eEnergy = (1.-r)*pairEnergy*0.5;
643 G4double pEnergy = pairEnergy - eEnergy; 613 G4double pEnergy = pairEnergy - eEnergy;
644 614
645 // Sample angles 615 // Sample angles
646 G4ThreeVector eDirection, pDirection; 616 G4ThreeVector eDirection, pDirection;
647 // 617 //
648 GetAngularDistribution()->SamplePairDirectio 618 GetAngularDistribution()->SamplePairDirections(aDynamicParticle,
649 619 eEnergy, pEnergy,
650 620 eDirection, pDirection);
651 // create G4DynamicParticle object for e+e- 621 // create G4DynamicParticle object for e+e-
652 eEnergy = std::max(eEnergy - CLHEP::electron 622 eEnergy = std::max(eEnergy - CLHEP::electron_mass_c2, 0.0);
653 pEnergy = std::max(pEnergy - CLHEP::electron 623 pEnergy = std::max(pEnergy - CLHEP::electron_mass_c2, 0.0);
654 auto aParticle1 = new G4DynamicParticle(theE << 624 G4DynamicParticle* aParticle1 =
655 auto aParticle2 = new G4DynamicParticle(theP << 625 new G4DynamicParticle(theElectron,eDirection,eEnergy);
>> 626 G4DynamicParticle* aParticle2 =
>> 627 new G4DynamicParticle(thePositron,pDirection,pEnergy);
656 // Fill output vector 628 // Fill output vector
657 vdp->push_back(aParticle1); 629 vdp->push_back(aParticle1);
658 vdp->push_back(aParticle2); 630 vdp->push_back(aParticle2);
659 631
660 // primary change 632 // primary change
661 kinEnergy -= pairEnergy; 633 kinEnergy -= pairEnergy;
662 partDirection *= totalMomentum; 634 partDirection *= totalMomentum;
663 partDirection -= (aParticle1->GetMomentum() 635 partDirection -= (aParticle1->GetMomentum() + aParticle2->GetMomentum());
664 partDirection = partDirection.unit(); 636 partDirection = partDirection.unit();
665 637
666 // if energy transfer is higher than thresho 638 // if energy transfer is higher than threshold (very high by default)
667 // then stop tracking the primary particle a 639 // then stop tracking the primary particle and create a new secondary
668 if (pairEnergy > SecondaryThreshold()) { 640 if (pairEnergy > SecondaryThreshold()) {
669 fParticleChange->ProposeTrackStatus(fStopA 641 fParticleChange->ProposeTrackStatus(fStopAndKill);
670 fParticleChange->SetProposedKineticEnergy( 642 fParticleChange->SetProposedKineticEnergy(0.0);
671 auto newdp = new G4DynamicParticle(particl << 643 G4DynamicParticle* newdp =
>> 644 new G4DynamicParticle(particle, partDirection, kinEnergy);
672 vdp->push_back(newdp); 645 vdp->push_back(newdp);
673 } else { // continue tracking the primary e- 646 } else { // continue tracking the primary e-/e+ otherwise
674 fParticleChange->SetProposedMomentumDirect 647 fParticleChange->SetProposedMomentumDirection(partDirection);
675 fParticleChange->SetProposedKineticEnergy( 648 fParticleChange->SetProposedKineticEnergy(kinEnergy);
676 } 649 }
677 //G4cout << "-- G4MuPairProductionModel::Sam 650 //G4cout << "-- G4MuPairProductionModel::SampleSecondaries done" << G4endl;
678 } 651 }
679 652
680 //....oooOO0OOooo........oooOO0OOooo........oo 653 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
681 654
682 G4double 655 G4double
683 G4MuPairProductionModel::FindScaledEnergy(G4in << 656 G4MuPairProductionModel::FindScaledEnergy(G4int Z, G4double rand,
684 G4double logTkin, 657 G4double logTkin,
685 G4double yymin, G4double yymax) 658 G4double yymin, G4double yymax)
686 { 659 {
687 G4double res = yymin; 660 G4double res = yymin;
688 G4Physics2DVector* pv = fElementData->GetEle << 661 G4Physics2DVector* pv = fElementData->GetElement2DData(Z);
689 if (nullptr != pv) { << 662 if(nullptr != pv) {
690 G4double pmin = pv->Value(yymin, logTkin); 663 G4double pmin = pv->Value(yymin, logTkin);
691 G4double pmax = pv->Value(yymax, logTkin); 664 G4double pmax = pv->Value(yymax, logTkin);
692 G4double p0 = pv->Value(0.0, logTkin); 665 G4double p0 = pv->Value(0.0, logTkin);
693 if(p0 <= 0.0) { DataCorrupted(ZDATPAIR[iz] << 666 if(p0 <= 0.0) { DataCorrupted(Z, logTkin); }
694 else { res = pv->FindLinearX((pmin + rand* 667 else { res = pv->FindLinearX((pmin + rand*(pmax - pmin))/p0, logTkin); }
695 } else { 668 } else {
696 DataCorrupted(ZDATPAIR[iz], logTkin); << 669 DataCorrupted(Z, logTkin);
697 } 670 }
698 return res; 671 return res;
699 } 672 }
700 673
701 //....oooOO0OOooo........oooOO0OOooo........oo 674 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
702 675
703 void G4MuPairProductionModel::DataCorrupted(G4 676 void G4MuPairProductionModel::DataCorrupted(G4int Z, G4double logTkin) const
704 { 677 {
705 G4ExceptionDescription ed; 678 G4ExceptionDescription ed;
706 ed << "G4ElementData is not properly initial 679 ed << "G4ElementData is not properly initialized Z= " << Z
707 << " Ekin(MeV)= " << G4Exp(logTkin) 680 << " Ekin(MeV)= " << G4Exp(logTkin)
708 << " IsMasterThread= " << IsMaster() 681 << " IsMasterThread= " << IsMaster()
709 << " Model " << GetName(); 682 << " Model " << GetName();
710 G4Exception("G4MuPairProductionModel::()", " 683 G4Exception("G4MuPairProductionModel::()", "em0033", FatalException, ed, "");
711 } 684 }
712 685
713 //....oooOO0OOooo........oooOO0OOooo........oo 686 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
714 687
715 void G4MuPairProductionModel::StoreTables() co 688 void G4MuPairProductionModel::StoreTables() const
716 { 689 {
717 for (G4int iz=0; iz<NZDATPAIR; ++iz) { << 690 for (G4int iz=0; iz<nzdat; ++iz) {
718 G4int Z = ZDATPAIR[iz]; << 691 G4int Z = zdat[iz];
719 G4Physics2DVector* pv = fElementData->GetE 692 G4Physics2DVector* pv = fElementData->GetElement2DData(Z);
720 if(nullptr == pv) { 693 if(nullptr == pv) {
721 DataCorrupted(Z, 1.0); 694 DataCorrupted(Z, 1.0);
722 return; 695 return;
723 } 696 }
724 std::ostringstream ss; 697 std::ostringstream ss;
725 ss << "mupair/" << particle->GetParticleNa 698 ss << "mupair/" << particle->GetParticleName() << Z << ".dat";
726 std::ofstream outfile(ss.str()); 699 std::ofstream outfile(ss.str());
727 pv->Store(outfile); 700 pv->Store(outfile);
728 } 701 }
729 } 702 }
730 703
731 //....oooOO0OOooo........oooOO0OOooo........oo 704 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
732 705
733 G4bool G4MuPairProductionModel::RetrieveTables 706 G4bool G4MuPairProductionModel::RetrieveTables()
734 { 707 {
735 for (G4int iz=0; iz<NZDATPAIR; ++iz) { << 708 const char* path = G4FindDataDir("G4LEDATA");
736 G4double Z = ZDATPAIR[iz]; << 709 G4String dir("");
>> 710 if (path) {
>> 711 std::ostringstream ost;
>> 712 ost << path << "/mupair/";
>> 713 dir = ost.str();
>> 714 } else {
>> 715 dir = "./mupair/";
>> 716 }
>> 717
>> 718 for (G4int iz=0; iz<nzdat; ++iz) {
>> 719 G4double Z = zdat[iz];
737 G4Physics2DVector* pv = new G4Physics2DVec 720 G4Physics2DVector* pv = new G4Physics2DVector(nbiny+1,nbine+1);
738 std::ostringstream ss; 721 std::ostringstream ss;
739 ss << G4EmParameters::Instance()->GetDirLE << 722 ss << dir << particle->GetParticleName() << Z << ".dat";
740 << particle->GetParticleName() << Z << <<
741 std::ifstream infile(ss.str(), std::ios::i 723 std::ifstream infile(ss.str(), std::ios::in);
742 if(!pv->Retrieve(infile)) { 724 if(!pv->Retrieve(infile)) {
743 delete pv; 725 delete pv;
744 return false; 726 return false;
745 } 727 }
746 fElementData->InitialiseForElement(iz, pv) << 728 fElementData->InitialiseForElement(Z, pv);
747 } 729 }
748 return true; 730 return true;
749 } 731 }
750 732
751 //....oooOO0OOooo........oooOO0OOooo........oo 733 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
752 734