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 // ------------ G4GammaConversionToMuo 27 // ------------ G4GammaConversionToMuons physics process ------ 28 // by H.Burkhardt, S. Kelner and R. Ko 28 // by H.Burkhardt, S. Kelner and R. Kokoulin, April 2002 29 // 29 // 30 // 30 // 31 // 07-08-02: missprint in OR condition in DoIt 31 // 07-08-02: missprint in OR condition in DoIt : f1<0 || f1>f1_max ..etc ... 32 // 25-10-04: migrade to new interfaces of Part 32 // 25-10-04: migrade to new interfaces of ParticleChange (vi) 33 // ------------------------------------------- 33 // --------------------------------------------------------------------------- 34 34 35 #include "G4GammaConversionToMuons.hh" 35 #include "G4GammaConversionToMuons.hh" 36 << 36 #include "G4PhysicalConstants.hh" 37 #include "G4BetheHeitler5DModel.hh" << 37 #include "G4SystemOfUnits.hh" 38 #include "G4Electron.hh" << 38 #include "G4UnitsTable.hh" 39 #include "G4EmParameters.hh" << 39 #include "G4MuonPlus.hh" >> 40 #include "G4MuonMinus.hh" 40 #include "G4EmProcessSubType.hh" 41 #include "G4EmProcessSubType.hh" 41 #include "G4Exp.hh" << 42 #include "G4EmParameters.hh" 42 #include "G4Gamma.hh" << 43 #include "G4Log.hh" << 44 #include "G4LossTableManager.hh" 43 #include "G4LossTableManager.hh" 45 #include "G4MuonMinus.hh" << 44 #include "G4BetheHeitler5DModel.hh" 46 #include "G4MuonPlus.hh" << 45 #include "G4Gamma.hh" 47 #include "G4NistManager.hh" << 46 #include "G4Electron.hh" 48 #include "G4PhysicalConstants.hh" << 49 #include "G4Positron.hh" 47 #include "G4Positron.hh" >> 48 #include "G4NistManager.hh" >> 49 #include "G4Log.hh" >> 50 #include "G4Exp.hh" 50 #include "G4ProductionCutsTable.hh" 51 #include "G4ProductionCutsTable.hh" 51 #include "G4SystemOfUnits.hh" << 52 #include "G4UnitsTable.hh" << 53 52 54 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 55 54 56 static const G4double sqrte = std::sqrt(std::e << 55 using namespace std; 57 static const G4double PowSat = -0.88; << 56 >> 57 static const G4double sqrte=sqrt(exp(1.)); >> 58 static const G4double PowSat=-0.88; 58 59 59 G4GammaConversionToMuons::G4GammaConversionToM 60 G4GammaConversionToMuons::G4GammaConversionToMuons(const G4String& processName, 60 G4ProcessType type) 61 G4ProcessType type) 61 : G4VDiscreteProcess (processName, type), 62 : G4VDiscreteProcess (processName, type), 62 Mmuon(G4MuonPlus::MuonPlus()->GetPDGMass() 63 Mmuon(G4MuonPlus::MuonPlus()->GetPDGMass()), 63 Rc(CLHEP::elm_coupling / Mmuon), << 64 Rc(elm_coupling/Mmuon), 64 LimitEnergy(5. * Mmuon), << 65 LimitEnergy (5.*Mmuon), 65 LowestEnergyLimit(2. * Mmuon), << 66 LowestEnergyLimit (2.*Mmuon), 66 HighestEnergyLimit(1e12 * CLHEP::GeV), // << 67 HighestEnergyLimit(1e12*GeV), // ok to 1e12GeV, then LPM suppression >> 68 Energy5DLimit(0.0), >> 69 CrossSecFactor(1.), >> 70 f5Dmodel(nullptr), 67 theGamma(G4Gamma::Gamma()), 71 theGamma(G4Gamma::Gamma()), 68 theMuonPlus(G4MuonPlus::MuonPlus()), 72 theMuonPlus(G4MuonPlus::MuonPlus()), 69 theMuonMinus(G4MuonMinus::MuonMinus()) 73 theMuonMinus(G4MuonMinus::MuonMinus()) 70 { << 74 { 71 SetProcessSubType(fGammaConversionToMuMu); 75 SetProcessSubType(fGammaConversionToMuMu); >> 76 MeanFreePath = DBL_MAX; 72 fManager = G4LossTableManager::Instance(); 77 fManager = G4LossTableManager::Instance(); 73 fManager->Register(this); 78 fManager->Register(this); 74 } 79 } 75 80 76 //....oooOO0OOooo........oooOO0OOooo........oo 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 77 82 78 G4GammaConversionToMuons::~G4GammaConversionTo << 83 G4GammaConversionToMuons::~G4GammaConversionToMuons() 79 { 84 { 80 fManager->DeRegister(this); 85 fManager->DeRegister(this); 81 } 86 } 82 87 83 //....oooOO0OOooo........oooOO0OOooo........oo 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 84 89 85 G4bool G4GammaConversionToMuons::IsApplicable( 90 G4bool G4GammaConversionToMuons::IsApplicable(const G4ParticleDefinition& part) 86 { 91 { 87 return (&part == theGamma); 92 return (&part == theGamma); 88 } 93 } 89 94 90 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 91 96 92 void G4GammaConversionToMuons::BuildPhysicsTab 97 void G4GammaConversionToMuons::BuildPhysicsTable(const G4ParticleDefinition& p) 93 { << 98 // Build cross section and mean free path tables >> 99 { //here no tables, just calling PrintInfoDefinition 94 Energy5DLimit = G4EmParameters::Instance()-> 100 Energy5DLimit = G4EmParameters::Instance()->MaxEnergyFor5DMuPair(); 95 << 101 if(Energy5DLimit > 0.0 && !f5Dmodel) { 96 auto table = G4Material::GetMaterialTable(); << 97 std::size_t nelm = 0; << 98 for (auto const& mat : *table) { << 99 std::size_t n = mat->GetNumberOfElements() << 100 nelm = std::max(nelm, n); << 101 } << 102 temp.resize(nelm, 0); << 103 << 104 if (Energy5DLimit > 0.0 && nullptr != f5Dmod << 105 f5Dmodel = new G4BetheHeitler5DModel(); 102 f5Dmodel = new G4BetheHeitler5DModel(); 106 f5Dmodel->SetLeptonPair(theMuonPlus, theMu 103 f5Dmodel->SetLeptonPair(theMuonPlus, theMuonMinus); 107 const std::size_t numElems = G4ProductionC << 104 const size_t numElems = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize(); 108 const G4DataVector cuts(numElems); 105 const G4DataVector cuts(numElems); 109 f5Dmodel->Initialise(&p, cuts); 106 f5Dmodel->Initialise(&p, cuts); 110 } 107 } 111 PrintInfoDefinition(); 108 PrintInfoDefinition(); 112 } 109 } 113 110 114 //....oooOO0OOooo........oooOO0OOooo........oo 111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 115 112 116 G4double G4GammaConversionToMuons::GetMeanFree << 113 G4double G4GammaConversionToMuons::GetMeanFreePath(const G4Track& aTrack, 117 << 114 G4double, G4ForceCondition*) >> 115 118 // returns the photon mean free path in GEANT4 116 // returns the photon mean free path in GEANT4 internal units >> 117 // (MeanFreePath is a private member of the class) >> 118 119 { 119 { 120 const G4DynamicParticle* aDynamicGamma = aTr << 120 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle(); 121 G4double GammaEnergy = aDynamicGamma->GetKin << 121 G4double GammaEnergy = aDynamicGamma->GetKineticEnergy(); 122 const G4Material* aMaterial = aTrack.GetMate << 122 const G4Material* aMaterial = aTrack.GetMaterial(); 123 return ComputeMeanFreePath(GammaEnergy, aMat << 123 >> 124 MeanFreePath = (GammaEnergy <= LowestEnergyLimit) >> 125 ? DBL_MAX : ComputeMeanFreePath(GammaEnergy,aMaterial); >> 126 >> 127 return MeanFreePath; 124 } 128 } 125 129 126 //....oooOO0OOooo........oooOO0OOooo........oo 130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 127 131 128 G4double 132 G4double 129 G4GammaConversionToMuons::ComputeMeanFreePath( 133 G4GammaConversionToMuons::ComputeMeanFreePath(G4double GammaEnergy, 130 134 const G4Material* aMaterial) 131 135 132 // computes and returns the photon mean free p 136 // computes and returns the photon mean free path in GEANT4 internal units 133 { 137 { 134 if(GammaEnergy <= LowestEnergyLimit) { retur 138 if(GammaEnergy <= LowestEnergyLimit) { return DBL_MAX; } 135 const G4ElementVector* theElementVector = aM 139 const G4ElementVector* theElementVector = aMaterial->GetElementVector(); 136 const G4double* NbOfAtomsPerVolume = aMateri 140 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume(); 137 141 138 G4double SIGMA = 0.0; 142 G4double SIGMA = 0.0; 139 G4double fact = 1.0; 143 G4double fact = 1.0; 140 G4double e = GammaEnergy; 144 G4double e = GammaEnergy; 141 // low energy approximation as in Bethe-Heit 145 // low energy approximation as in Bethe-Heitler model 142 if(e < LimitEnergy) { 146 if(e < LimitEnergy) { 143 G4double y = (e - LowestEnergyLimit)/(Limi 147 G4double y = (e - LowestEnergyLimit)/(LimitEnergy - LowestEnergyLimit); 144 fact = y*y; 148 fact = y*y; 145 e = LimitEnergy; 149 e = LimitEnergy; 146 } 150 } 147 151 148 for ( std::size_t i=0 ; i < aMaterial->GetNu << 152 for ( size_t i=0 ; i < aMaterial->GetNumberOfElements(); ++i) 149 { 153 { 150 SIGMA += NbOfAtomsPerVolume[i] * fact * 154 SIGMA += NbOfAtomsPerVolume[i] * fact * 151 ComputeCrossSectionPerAtom(e, (*theEleme 155 ComputeCrossSectionPerAtom(e, (*theElementVector)[i]->GetZasInt()); 152 } 156 } 153 return (SIGMA > 0.0) ? 1./SIGMA : DBL_MAX; 157 return (SIGMA > 0.0) ? 1./SIGMA : DBL_MAX; 154 } 158 } 155 159 156 //....oooOO0OOooo........oooOO0OOooo........oo 160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 157 161 158 G4double G4GammaConversionToMuons::GetCrossSec 162 G4double G4GammaConversionToMuons::GetCrossSectionPerAtom( 159 const G4Dyn 163 const G4DynamicParticle* aDynamicGamma, 160 const G4Ele 164 const G4Element* anElement) 161 165 162 // gives the total cross section per atom in G 166 // gives the total cross section per atom in GEANT4 internal units 163 { 167 { 164 return ComputeCrossSectionPerAtom(aDynamicG 168 return ComputeCrossSectionPerAtom(aDynamicGamma->GetKineticEnergy(), 165 anElement 169 anElement->GetZasInt()); 166 } 170 } 167 171 168 //....oooOO0OOooo........oooOO0OOooo........oo 172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 169 173 170 G4double G4GammaConversionToMuons::ComputeCros 174 G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom( 171 G4double Egam, G4int 175 G4double Egam, G4int Z) 172 176 173 // Calculates the microscopic cross section in 177 // Calculates the microscopic cross section in GEANT4 internal units. 174 // Total cross section parametrisation from H. 178 // Total cross section parametrisation from H.Burkhardt 175 // It gives a good description at any energy ( 179 // It gives a good description at any energy (from 0 to 10**21 eV) 176 { 180 { 177 if(Egam <= LowestEnergyLimit) { return 0.0; 181 if(Egam <= LowestEnergyLimit) { return 0.0; } 178 182 179 G4NistManager* nist = G4NistManager::Instanc 183 G4NistManager* nist = G4NistManager::Instance(); 180 184 181 G4double PowThres, Ecor, B, Dn, Zthird, Winf << 185 G4double PowThres,Ecor,B,Dn,Zthird,Winfty,WMedAppr, 182 << 186 Wsatur,sigfac; 183 if (Z == 1) { // special case of Hydrogen << 187 184 B = 202.4; << 188 if(Z==1) // special case of Hydrogen 185 Dn = 1.49; << 189 { B=202.4; 186 } << 190 Dn=1.49; 187 else { << 191 } 188 B = 183.; << 192 else 189 Dn = 1.54 * nist->GetA27(Z); << 193 { B=183.; 190 } << 194 Dn=1.54*nist->GetA27(Z); 191 Zthird = 1. / nist->GetZ13(Z); // Z**(-1/3) << 195 } 192 Winfty = B * Zthird * Mmuon / (Dn * electron << 196 Zthird=1./nist->GetZ13(Z); // Z**(-1/3) 193 WMedAppr = 1. / (4. * Dn * sqrte * Mmuon); << 197 Winfty=B*Zthird*Mmuon/(Dn*electron_mass_c2); 194 Wsatur = Winfty / WMedAppr; << 198 WMedAppr=1./(4.*Dn*sqrte*Mmuon); 195 sigfac = 4. * fine_structure_const * Z * Z * << 199 Wsatur=Winfty/WMedAppr; 196 PowThres = 1.479 + 0.00799 * Dn; << 200 sigfac=4.*fine_structure_const*Z*Z*Rc*Rc; 197 Ecor = -18. + 4347. / (B * Zthird); << 201 PowThres=1.479+0.00799*Dn; 198 << 202 Ecor=-18.+4347./(B*Zthird); 199 G4double CorFuc = 1. + .04 * G4Log(1. + Ecor << 203 200 G4double Eg = << 204 G4double CorFuc=1.+.04*G4Log(1.+Ecor/Egam); 201 G4Exp(G4Log(1. - 4. * Mmuon / Egam) * PowT << 205 //G4double Eg=pow(1.-4.*Mmuon/Egam,PowThres)*pow( pow(Wsatur,PowSat)+ 202 * G4Exp(G4Log(G4Exp(G4Log(Wsatur) * PowSat << 206 // pow(Egam,PowSat),1./PowSat); // threshold and saturation 203 << 207 G4double Eg=G4Exp(G4Log(1.-4.*Mmuon/Egam)*PowThres)* 204 G4double CrossSection = 7. / 9. * sigfac * G << 208 G4Exp(G4Log( G4Exp(G4Log(Wsatur)*PowSat)+G4Exp(G4Log(Egam)*PowSat))/PowSat); 205 CrossSection *= CrossSecFactor; // increase << 209 G4double CrossSection=7./9.*sigfac*G4Log(1.+WMedAppr*CorFuc*Eg); >> 210 CrossSection*=CrossSecFactor; // increase the CrossSection by (by default 1) 206 return CrossSection; 211 return CrossSection; 207 } 212 } 208 213 209 //....oooOO0OOooo........oooOO0OOooo........oo 214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 210 215 211 void G4GammaConversionToMuons::SetCrossSecFact 216 void G4GammaConversionToMuons::SetCrossSecFactor(G4double fac) 212 // Set the factor to artificially increase the 217 // Set the factor to artificially increase the cross section 213 { << 218 { 214 if (fac < 0.0) return; << 219 CrossSecFactor=fac; 215 CrossSecFactor = fac; << 220 G4cout << "The cross section for GammaConversionToMuons is artificially " 216 if (verboseLevel > 1) { << 221 << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl; 217 G4cout << "The cross section for GammaConv << 218 << "increased by the CrossSecFactor << 219 } << 220 } 222 } 221 223 222 //....oooOO0OOooo........oooOO0OOooo........oo 224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 223 225 224 G4VParticleChange* G4GammaConversionToMuons::P 226 G4VParticleChange* G4GammaConversionToMuons::PostStepDoIt( 225 227 const G4Track& aTrack, 226 228 const G4Step& aStep) 227 // 229 // 228 // generation of gamma->mu+mu- 230 // generation of gamma->mu+mu- 229 // 231 // 230 { 232 { 231 aParticleChange.Initialize(aTrack); 233 aParticleChange.Initialize(aTrack); 232 const G4Material* aMaterial = aTrack.GetMate 234 const G4Material* aMaterial = aTrack.GetMaterial(); 233 235 234 // current Gamma energy and direction, retur 236 // current Gamma energy and direction, return if energy too low 235 const G4DynamicParticle* aDynamicGamma = aTr << 237 const G4DynamicParticle *aDynamicGamma = aTrack.GetDynamicParticle(); 236 G4double Egam = aDynamicGamma->GetKineticEne 238 G4double Egam = aDynamicGamma->GetKineticEnergy(); 237 if (Egam <= LowestEnergyLimit) { 239 if (Egam <= LowestEnergyLimit) { 238 return G4VDiscreteProcess::PostStepDoIt(aT 240 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep); 239 } 241 } 240 // 242 // 241 // Kill the incident photon 243 // Kill the incident photon 242 // 244 // 243 aParticleChange.ProposeMomentumDirection( 0. 245 aParticleChange.ProposeMomentumDirection( 0., 0., 0. ) ; 244 aParticleChange.ProposeEnergy( 0. ) ; 246 aParticleChange.ProposeEnergy( 0. ) ; 245 aParticleChange.ProposeTrackStatus( fStopAnd 247 aParticleChange.ProposeTrackStatus( fStopAndKill ) ; 246 248 247 if (Egam <= Energy5DLimit) { 249 if (Egam <= Energy5DLimit) { 248 std::vector<G4DynamicParticle*> fvect; 250 std::vector<G4DynamicParticle*> fvect; 249 f5Dmodel->SampleSecondaries(&fvect, aTrack 251 f5Dmodel->SampleSecondaries(&fvect, aTrack.GetMaterialCutsCouple(), 250 aTrack.GetDynamicParticle(), 0.0, DBL_ 252 aTrack.GetDynamicParticle(), 0.0, DBL_MAX); >> 253 aParticleChange.SetNumberOfSecondaries(fvect.size()); 251 for(auto dp : fvect) { aParticleChange.Add 254 for(auto dp : fvect) { aParticleChange.AddSecondary(dp); } 252 return G4VDiscreteProcess::PostStepDoIt(aT 255 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep); 253 } 256 } 254 257 255 G4ParticleMomentum GammaDirection = aDynamic 258 G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection(); 256 259 257 // select randomly one element constituting 260 // select randomly one element constituting the material 258 const G4Element* anElement = SelectRandomAto 261 const G4Element* anElement = SelectRandomAtom(aDynamicGamma, aMaterial); 259 G4int Z = anElement->GetZasInt(); 262 G4int Z = anElement->GetZasInt(); 260 G4NistManager* nist = G4NistManager::Instanc 263 G4NistManager* nist = G4NistManager::Instance(); 261 264 262 G4double B, Dn; << 265 G4double B,Dn; 263 G4double A027 = nist->GetA27(Z); 266 G4double A027 = nist->GetA27(Z); 264 267 265 if (Z == 1) { // special case of Hydrogen << 268 if(Z==1) // special case of Hydrogen 266 B = 202.4; << 269 { B=202.4; 267 Dn = 1.49; << 270 Dn=1.49; 268 } << 271 } 269 else { << 272 else 270 B = 183.; << 273 { B=183.; 271 Dn = 1.54 * A027; << 274 Dn=1.54*A027; 272 } << 275 } 273 G4double Zthird = 1. / nist->GetZ13(Z); // << 276 G4double Zthird=1./nist->GetZ13(Z); // Z**(-1/3) 274 G4double Winfty = B * Zthird * Mmuon / (Dn * << 277 G4double Winfty=B*Zthird*Mmuon/(Dn*electron_mass_c2); 275 278 276 G4double C1Num = 0.138 * A027; << 279 G4double C1Num=0.138*A027; 277 G4double C1Num2 = C1Num * C1Num; << 280 G4double C1Num2=C1Num*C1Num; 278 G4double C2Term2 = electron_mass_c2 / (183. << 281 G4double C2Term2=electron_mass_c2/(183.*Zthird*Mmuon); 279 282 280 G4double GammaMuonInv = Mmuon / Egam; << 283 G4double GammaMuonInv=Mmuon/Egam; 281 284 282 // generate xPlus according to the different 285 // generate xPlus according to the differential cross section by rejection 283 G4double xmin = (Egam <= LimitEnergy) ? 0.5 << 286 G4double xmin=(Egam < LimitEnergy) ? GammaMuonInv : .5-sqrt(.25-GammaMuonInv); 284 G4double xmax = 1. - xmin; << 287 G4double xmax=1.-xmin; 285 << 286 G4double Ds2 = (Dn * sqrte - 2.); << 287 G4double sBZ = sqrte * B * Zthird / electron << 288 G4double LogWmaxInv = << 289 1. / G4Log(Winfty * (1. + 2. * Ds2 * Gamma << 290 G4double xPlus = 0.5; << 291 G4double xMinus = 0.5; << 292 G4double xPM = 0.25; << 293 288 >> 289 G4double Ds2=(Dn*sqrte-2.); >> 290 G4double sBZ=sqrte*B*Zthird/electron_mass_c2; >> 291 G4double LogWmaxInv=1./G4Log(Winfty*(1.+2.*Ds2*GammaMuonInv) >> 292 /(1.+2.*sBZ*Mmuon*GammaMuonInv)); >> 293 G4double xPlus,xMinus,xPM,result,W; 294 G4int nn = 0; 294 G4int nn = 0; 295 const G4int nmax = 1000; 295 const G4int nmax = 1000; 296 << 296 do { 297 // sampling for Egam > LimitEnergy << 297 xPlus=xmin+G4UniformRand()*(xmax-xmin); 298 if (xmin < 0.5) { << 298 xMinus=1.-xPlus; 299 G4double result, W; << 299 xPM=xPlus*xMinus; 300 do { << 300 G4double del=Mmuon*Mmuon/(2.*Egam*xPM); 301 xPlus = xmin + G4UniformRand() * (xmax - << 301 W=Winfty*(1.+Ds2*del/Mmuon)/(1.+sBZ*del); 302 xMinus = 1. - xPlus; << 302 G4double xxp=1.-4./3.*xPM; // the main xPlus dependence 303 xPM = xPlus * xMinus; << 303 result=(xxp > 0.) ? xxp*G4Log(W)*LogWmaxInv : 0.0; 304 G4double del = Mmuon * Mmuon / (2. * Ega << 304 if(result>1.) { 305 W = Winfty * (1. + Ds2 * del / Mmuon) / << 305 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:" 306 G4double xxp = 1. - 4. / 3. * xPM; // t << 306 << " in dSigxPlusGen, result=" << result << " > 1" << G4endl; 307 result = (xxp > 0.) ? xxp * G4Log(W) * L << 308 if (result > 1.) { << 309 G4cout << "G4GammaConversionToMuons::P << 310 << " in dSigxPlusGen, result=" << 311 } << 312 ++nn; << 313 if(nn >= nmax) { break; } << 314 } 307 } 315 // Loop checking, 07-Aug-2015, Vladimir Iv << 308 ++nn; 316 while (G4UniformRand() > result); << 309 if(nn >= nmax) { break; } 317 } 310 } >> 311 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko >> 312 while (G4UniformRand() > result); 318 313 319 // now generate the angular variables via th 314 // now generate the angular variables via the auxilary variables t,psi,rho 320 G4double t; 315 G4double t; 321 G4double psi; 316 G4double psi; 322 G4double rho; 317 G4double rho; 323 318 324 G4double a3 = (GammaMuonInv / (2. * xPM)); << 319 G4double a3 = (GammaMuonInv/(2.*xPM)); 325 G4double a33 = a3 * a3; << 320 G4double a33 = a3*a3; 326 G4double f1; 321 G4double f1; 327 G4double b1 = 1./(4.*C1Num2); 322 G4double b1 = 1./(4.*C1Num2); 328 G4double b3 = b1*b1*b1; 323 G4double b3 = b1*b1*b1; 329 G4double a21 = a33 + b1; 324 G4double a21 = a33 + b1; 330 325 331 G4double f1_max=-(1.-xPM)*(2.*b1+(a21+a33)*G 326 G4double f1_max=-(1.-xPM)*(2.*b1+(a21+a33)*G4Log(a33/a21))/(2*b3); 332 327 333 G4double thetaPlus,thetaMinus,phiHalf; // fi 328 G4double thetaPlus,thetaMinus,phiHalf; // final angular variables 334 nn = 0; 329 nn = 0; 335 // t, psi, rho generation start (while angl 330 // t, psi, rho generation start (while angle < pi) 336 do { 331 do { 337 //generate t by the rejection method 332 //generate t by the rejection method 338 do { 333 do { 339 ++nn; 334 ++nn; 340 t=G4UniformRand(); 335 t=G4UniformRand(); 341 G4double a34=a33/(t*t); 336 G4double a34=a33/(t*t); 342 G4double a22 = a34 + b1; 337 G4double a22 = a34 + b1; 343 if(std::abs(b1)<0.0001*a34) { << 338 if(std::abs(b1)<0.0001*a34) 344 // special case of a34=a22 because of << 339 // special case of a34=a22 because of logarithm accuracy 345 f1=(1.-2.*xPM+4.*xPM*t*(1.-t))/(12.*a3 << 340 { 346 } << 341 f1=(1.-2.*xPM+4.*xPM*t*(1.-t))/(12.*a34*a34*a34*a34); 347 else { << 342 } 348 f1=-(1.-2.*xPM+4.*xPM*t*(1.-t))*(2.*b1 << 343 else 349 } << 344 { 350 if (f1 < 0.0 || f1 > f1_max) { // shoul << 345 f1=-(1.-2.*xPM+4.*xPM*t*(1.-t))*(2.*b1+(a22+a34)*G4Log(a34/a22))/(2*b3); 351 G4cout << "G4GammaConversionToMuons::P << 346 } 352 << "outside allowed range f1=" << f1 << 347 if(f1<0.0 || f1> f1_max) // should never happend 353 << " is set to zero, a34 = "<< a34 << << 348 { 354 << G4endl; << 349 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:" 355 f1 = 0.0; << 350 << "outside allowed range f1=" << f1 356 } << 351 << " is set to zero, a34 = "<< a34 << " a22 = "<<a22<<"." >> 352 << G4endl; >> 353 f1 = 0.0; >> 354 } 357 if(nn > nmax) { break; } 355 if(nn > nmax) { break; } 358 // Loop checking, 07-Aug-2015, Vladimir 356 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko 359 } while ( G4UniformRand()*f1_max > f1); 357 } while ( G4UniformRand()*f1_max > f1); 360 // generate psi by the rejection method 358 // generate psi by the rejection method 361 G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t)) 359 G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t)); 362 // long version 360 // long version 363 G4double f2; 361 G4double f2; 364 do { 362 do { 365 ++nn; 363 ++nn; 366 psi=twopi*G4UniformRand(); 364 psi=twopi*G4UniformRand(); 367 f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+std::co << 365 f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+cos(2.*psi)); 368 if(f2<0 || f2> f2_max) { // should never << 366 if(f2<0 || f2> f2_max) // should never happend 369 G4cout << "G4GammaConversionToMuons::P << 367 { 370 << "outside allowed range f2=" << 368 G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:" 371 f2 = 0.0; << 369 << "outside allowed range f2=" << f2 << " is set to zero" 372 } << 370 << G4endl; >> 371 f2 = 0.0; >> 372 } 373 if(nn >= nmax) { break; } 373 if(nn >= nmax) { break; } 374 // Loop checking, 07-Aug-2015, Vladimir 374 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko 375 } while ( G4UniformRand()*f2_max > f2); 375 } while ( G4UniformRand()*f2_max > f2); 376 376 377 // generate rho by direct transformation 377 // generate rho by direct transformation 378 G4double C2Term1=GammaMuonInv/(2.*xPM*t); 378 G4double C2Term1=GammaMuonInv/(2.*xPM*t); 379 G4double C22 = C2Term1*C2Term1+C2Term2*C2T 379 G4double C22 = C2Term1*C2Term1+C2Term2*C2Term2; 380 G4double C2=4.*C22*C22/std::sqrt(xPM); << 380 G4double C2=4.*C22*C22/sqrt(xPM); 381 G4double rhomax=(1./t-1.)*1.9/A027; 381 G4double rhomax=(1./t-1.)*1.9/A027; 382 G4double beta=G4Log( (C2+rhomax*rhomax*rho 382 G4double beta=G4Log( (C2+rhomax*rhomax*rhomax*rhomax)/C2 ); 383 rho=G4Exp(G4Log(C2 *( G4Exp(beta*G4Uniform 383 rho=G4Exp(G4Log(C2 *( G4Exp(beta*G4UniformRand())-1. ))*0.25); 384 384 385 //now get from t and psi the kinematical v 385 //now get from t and psi the kinematical variables 386 G4double u=std::sqrt(1./t-1.); << 386 G4double u=sqrt(1./t-1.); 387 G4double xiHalf=0.5*rho*std::cos(psi); << 387 G4double xiHalf=0.5*rho*cos(psi); 388 phiHalf=0.5*rho/u*std::sin(psi); << 388 phiHalf=0.5*rho/u*sin(psi); 389 389 390 thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus; 390 thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus; 391 thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus; 391 thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus; 392 392 393 // protection against infinite loop 393 // protection against infinite loop 394 if(nn > nmax) { 394 if(nn > nmax) { 395 if(std::abs(thetaPlus)>pi) { thetaPlus = 395 if(std::abs(thetaPlus)>pi) { thetaPlus = 0.0; } 396 if(std::abs(thetaMinus)>pi) { thetaMinus 396 if(std::abs(thetaMinus)>pi) { thetaMinus = 0.0; } 397 } 397 } 398 398 399 // Loop checking, 07-Aug-2015, Vladimir Iv 399 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko 400 } while ( std::abs(thetaPlus)>pi || std::abs 400 } while ( std::abs(thetaPlus)>pi || std::abs(thetaMinus) >pi); 401 401 402 // now construct the vectors 402 // now construct the vectors 403 // azimuthal symmetry, take phi0 at random b 403 // azimuthal symmetry, take phi0 at random between 0 and 2 pi 404 G4double phi0=twopi*G4UniformRand(); 404 G4double phi0=twopi*G4UniformRand(); 405 G4double EPlus=xPlus*Egam; 405 G4double EPlus=xPlus*Egam; 406 G4double EMinus=xMinus*Egam; 406 G4double EMinus=xMinus*Egam; 407 407 408 // mu+ mu- directions for gamma in z-directi 408 // mu+ mu- directions for gamma in z-direction 409 G4ThreeVector MuPlusDirection ( std::sin(th << 409 G4ThreeVector MuPlusDirection ( sin(thetaPlus) *cos(phi0+phiHalf), 410 std::sin(thetaPlus) *std:: << 410 sin(thetaPlus) *sin(phi0+phiHalf), cos(thetaPlus) ); 411 G4ThreeVector MuMinusDirection (-std::sin(th << 411 G4ThreeVector MuMinusDirection (-sin(thetaMinus)*cos(phi0-phiHalf), 412 -std::sin(thetaMinus) *std:: << 412 -sin(thetaMinus) *sin(phi0-phiHalf), cos(thetaMinus) ); 413 // rotate to actual gamma direction 413 // rotate to actual gamma direction 414 MuPlusDirection.rotateUz(GammaDirection); 414 MuPlusDirection.rotateUz(GammaDirection); 415 MuMinusDirection.rotateUz(GammaDirection); 415 MuMinusDirection.rotateUz(GammaDirection); 416 << 416 aParticleChange.SetNumberOfSecondaries(2); 417 // create G4DynamicParticle object for the p 417 // create G4DynamicParticle object for the particle1 418 auto aParticle1 = new G4DynamicParticle(theM << 418 G4DynamicParticle* aParticle1 = >> 419 new G4DynamicParticle(theMuonPlus,MuPlusDirection,EPlus-Mmuon); 419 aParticleChange.AddSecondary(aParticle1); 420 aParticleChange.AddSecondary(aParticle1); 420 // create G4DynamicParticle object for the p 421 // create G4DynamicParticle object for the particle2 421 auto aParticle2 = new G4DynamicParticle(theM << 422 G4DynamicParticle* aParticle2 = >> 423 new G4DynamicParticle(theMuonMinus,MuMinusDirection,EMinus-Mmuon); 422 aParticleChange.AddSecondary(aParticle2); 424 aParticleChange.AddSecondary(aParticle2); 423 // Reset NbOfInteractionLengthLeft and retu 425 // Reset NbOfInteractionLengthLeft and return aParticleChange 424 return G4VDiscreteProcess::PostStepDoIt( aTr 426 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep ); 425 } 427 } 426 428 427 //....oooOO0OOooo........oooOO0OOooo........oo 429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 428 430 429 const G4Element* G4GammaConversionToMuons::Sel 431 const G4Element* G4GammaConversionToMuons::SelectRandomAtom( 430 const G4DynamicParticle* aDynamicGamma, 432 const G4DynamicParticle* aDynamicGamma, 431 const G4Material* aMaterial) 433 const G4Material* aMaterial) 432 { 434 { 433 // select randomly 1 element within the mate 435 // select randomly 1 element within the material, invoked by PostStepDoIt 434 436 435 const std::size_t NumberOfElements = aM << 437 const G4int NumberOfElements = aMaterial->GetNumberOfElements(); 436 const G4ElementVector* theElementVector = aM 438 const G4ElementVector* theElementVector = aMaterial->GetElementVector(); 437 const G4Element* elm = (*theElementVector)[0 439 const G4Element* elm = (*theElementVector)[0]; 438 440 439 if (NumberOfElements > 1) { << 441 if (NumberOfElements > 1) { 440 G4double e = std::max(aDynamicGamma->GetKi << 442 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume(); 441 const G4double* natom = aMaterial->GetVecN << 443 >> 444 G4double PartialSumSigma = 0.; >> 445 G4double rval = G4UniformRand()/MeanFreePath; 442 446 443 G4double sum = 0.; << 447 for (G4int i=0; i<NumberOfElements; ++i) 444 for (std::size_t i=0; i<NumberOfElements; << 448 { 445 elm = (*theElementVector)[i]; 449 elm = (*theElementVector)[i]; 446 sum += natom[i]*ComputeCrossSectionPerAt << 450 PartialSumSigma += NbOfAtomsPerVolume[i] 447 temp[i] = sum; << 451 *GetCrossSectionPerAtom(aDynamicGamma, elm); 448 } << 452 if (rval <= PartialSumSigma) { break; } 449 sum *= G4UniformRand(); << 450 for (std::size_t i=0; i<NumberOfElements; << 451 if(sum <= temp[i]) { << 452 elm = (*theElementVector)[i]; << 453 break; << 454 } << 455 } 453 } 456 } 454 } 457 return elm; 455 return elm; 458 } 456 } 459 457 460 //....oooOO0OOooo........oooOO0OOooo........oo 458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 461 459 462 void G4GammaConversionToMuons::PrintInfoDefini 460 void G4GammaConversionToMuons::PrintInfoDefinition() 463 { 461 { 464 G4String comments = "gamma->mu+mu- Bethe Hei << 462 G4String comments ="gamma->mu+mu- Bethe Heitler process, SubType= "; 465 G4cout << G4endl << GetProcessName() << ": << 463 G4cout << G4endl << GetProcessName() << ": " << comments >> 464 << GetProcessSubType() << G4endl; 466 G4cout << " good cross section parame 465 G4cout << " good cross section parametrization from " 467 << G4BestUnit(LowestEnergyLimit, "Ene << 466 << G4BestUnit(LowestEnergyLimit,"Energy") 468 << " GeV for all Z." << G4endl; << 467 << " to " << HighestEnergyLimit/GeV << " GeV for all Z." << G4endl; 469 G4cout << " cross section factor: " < << 470 } 468 } 471 469 472 //....oooOO0OOooo........oooOO0OOooo........oo 470 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 473 471