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