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