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