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