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