Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // 28 // ------------ G4AnnihiToMuPair physics process ------ 29 // by H.Burkhardt, S. Kelner and R. Kokoulin, November 2002 30 // ----------------------------------------------------------------------------- 31 // 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......// 33 // 34 // 27.01.03 : first implementation (hbu) 35 // 04.02.03 : cosmetic simplifications (mma) 36 // 25.10.04 : migrade to new interfaces of ParticleChange (vi) 37 // 28.02.18 : cross section now including SSS threshold factor 38 // 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 40 41 #include "G4AnnihiToMuPair.hh" 42 43 #include "G4Exp.hh" 44 #include "G4LossTableManager.hh" 45 #include "G4Material.hh" 46 #include "G4MuonMinus.hh" 47 #include "G4MuonPlus.hh" 48 #include "G4PhysicalConstants.hh" 49 #include "G4Positron.hh" 50 #include "G4Step.hh" 51 #include "G4SystemOfUnits.hh" 52 #include "G4TauMinus.hh" 53 #include "G4TauPlus.hh" 54 #include "G4ios.hh" 55 #include "Randomize.hh" 56 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 58 59 G4AnnihiToMuPair::G4AnnihiToMuPair(const G4String& processName, 60 G4ProcessType type):G4VDiscreteProcess (processName, type) 61 { 62 //e+ Energy threshold 63 if(processName == "AnnihiToTauPair") { 64 SetProcessSubType(fAnnihilationToTauTau); 65 part1 = G4TauPlus::TauPlus(); 66 part2 = G4TauMinus::TauMinus(); 67 fInfo = "e+e->tau+tau-"; 68 } else { 69 SetProcessSubType(fAnnihilationToMuMu); 70 part1 = G4MuonPlus::MuonPlus(); 71 part2 = G4MuonMinus::MuonMinus(); 72 } 73 fMass = part1->GetPDGMass(); 74 fLowEnergyLimit = 2. * fMass * fMass / CLHEP::electron_mass_c2 - CLHEP::electron_mass_c2; 75 76 // model is ok up to 1000 TeV due to neglected Z-interference 77 fHighEnergyLimit = 1000. * TeV; 78 79 fCurrentSigma = 0.0; 80 fCrossSecFactor = 1.; 81 fManager = G4LossTableManager::Instance(); 82 fManager->Register(this); 83 } 84 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 86 87 G4AnnihiToMuPair::~G4AnnihiToMuPair() // (empty) destructor 88 { 89 fManager->DeRegister(this); 90 } 91 92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 93 94 G4bool G4AnnihiToMuPair::IsApplicable(const G4ParticleDefinition& particle) 95 { 96 return ( &particle == G4Positron::Positron() ); 97 } 98 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 100 101 void G4AnnihiToMuPair::BuildPhysicsTable(const G4ParticleDefinition&) 102 { 103 PrintInfoDefinition(); 104 } 105 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 107 108 void G4AnnihiToMuPair::SetCrossSecFactor(G4double fac) 109 // Set the factor to artificially increase the cross section 110 { 111 fCrossSecFactor = fac; 112 //G4cout << "The cross section for AnnihiToMuPair is artificially " 113 // << "increased by the CrossSecFactor=" << fCrossSecFactor << G4endl; 114 } 115 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 117 118 G4double G4AnnihiToMuPair::ComputeCrossSectionPerElectron(const G4double e) 119 // Calculates the microscopic cross section in GEANT4 internal units. 120 // It gives a good description from threshold to 1000 GeV 121 { 122 G4double rmuon = CLHEP::elm_coupling/fMass; //classical particle radius 123 G4double sig0 = CLHEP::pi*rmuon*rmuon/3.; //constant in crossSection 124 const G4double pial = CLHEP::pi*CLHEP::fine_structure_const; // pi * alphaQED 125 126 if (e <= fLowEnergyLimit) return 0.0; 127 128 const G4double xi = fLowEnergyLimit/e; 129 const G4double piaxi = pial * std::sqrt(xi); 130 G4double sigma = sig0 * xi * (1. + xi*0.5); 131 //G4cout << "### xi= " << xi << " piaxi=" << piaxi << G4endl; 132 133 // argument of the exponent below 0.1 or above 10 134 // Sigma per electron * number of electrons per atom 135 if(xi <= 1.0 - 100*piaxi*piaxi) { 136 sigma *= std::sqrt(1.0 - xi); 137 } 138 else if (xi >= 1.0 - 0.01 * piaxi * piaxi) { 139 sigma *= piaxi; 140 } 141 else { 142 sigma *= piaxi / (1. - G4Exp(-piaxi / std::sqrt(1 - xi))); 143 } 144 // G4cout << "### sigma= " << sigma << G4endl; 145 return sigma; 146 } 147 148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 149 150 G4double G4AnnihiToMuPair::ComputeCrossSectionPerAtom(const G4double energy, 151 const G4double Z) 152 { 153 return ComputeCrossSectionPerElectron(energy)*Z; 154 } 155 156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 157 158 G4double G4AnnihiToMuPair::CrossSectionPerVolume(G4double energy, 159 const G4Material* aMaterial) 160 { 161 return ComputeCrossSectionPerElectron(energy)*aMaterial->GetTotNbOfElectPerVolume(); 162 } 163 164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 165 166 G4double G4AnnihiToMuPair::GetMeanFreePath(const G4Track& aTrack, 167 G4double, G4ForceCondition*) 168 // returns the positron mean free path in GEANT4 internal units 169 { 170 const G4DynamicParticle* aDynamicPositron = aTrack.GetDynamicParticle(); 171 G4double energy = aDynamicPositron->GetTotalEnergy(); 172 const G4Material* aMaterial = aTrack.GetMaterial(); 173 174 // cross section before step 175 fCurrentSigma = CrossSectionPerVolume(energy, aMaterial); 176 177 // increase the CrossSection by CrossSecFactor (default 1) 178 return (fCurrentSigma > 0.0) ? 1.0/(fCurrentSigma*fCrossSecFactor) : DBL_MAX; 179 } 180 181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 182 183 G4VParticleChange* G4AnnihiToMuPair::PostStepDoIt(const G4Track& aTrack, 184 const G4Step& aStep) 185 // 186 // generation of e+e- -> mu+mu- 187 // 188 { 189 aParticleChange.Initialize(aTrack); 190 191 // current Positron energy and direction, return if energy too low 192 const G4DynamicParticle *aDynamicPositron = aTrack.GetDynamicParticle(); 193 const G4double Mele = CLHEP::electron_mass_c2; 194 G4double Epos = aDynamicPositron->GetTotalEnergy(); 195 G4double xs = CrossSectionPerVolume(Epos, aTrack.GetMaterial()); 196 197 // test of cross section 198 if(xs > 0.0 && fCurrentSigma*G4UniformRand() > xs) { 199 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep); 200 } 201 202 const G4ThreeVector PosiDirection = aDynamicPositron->GetMomentumDirection(); 203 G4double xi = fLowEnergyLimit/Epos; // xi is always less than 1, 204 // goes to 0 at high Epos 205 206 // generate cost; probability function 1+cost**2 at high Epos 207 // 208 G4double cost; 209 do { cost = 2.*G4UniformRand()-1.; } 210 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko 211 while (2.*G4UniformRand() > 1.+xi+cost*cost*(1.-xi) ); 212 G4double sint = std::sqrt(1.-cost*cost); 213 214 // generate phi 215 // 216 G4double phi = 2.*CLHEP::pi*G4UniformRand(); 217 218 G4double Ecm = std::sqrt(0.5*Mele*(Epos+Mele)); 219 G4double Pcm = std::sqrt(Ecm*Ecm - fMass*fMass); 220 G4double beta = std::sqrt((Epos-Mele)/(Epos+Mele)); 221 G4double gamma = Ecm/Mele; 222 G4double Pt = Pcm*sint; 223 224 // energy and momentum of the muons in the Lab 225 // 226 G4double EmuPlus = gamma*(Ecm + cost*beta*Pcm); 227 G4double EmuMinus = gamma*(Ecm - cost*beta*Pcm); 228 G4double PmuPlusZ = gamma*(beta*Ecm + cost*Pcm); 229 G4double PmuMinusZ = gamma*(beta*Ecm - cost*Pcm); 230 G4double PmuPlusX = Pt*std::cos(phi); 231 G4double PmuPlusY = Pt*std::sin(phi); 232 G4double PmuMinusX =-PmuPlusX; 233 G4double PmuMinusY =-PmuPlusY; 234 // absolute momenta 235 G4double PmuPlus = std::sqrt(Pt*Pt+PmuPlusZ *PmuPlusZ ); 236 G4double PmuMinus = std::sqrt(Pt*Pt+PmuMinusZ*PmuMinusZ); 237 238 // mu+ mu- directions for Positron in z-direction 239 // 240 G4ThreeVector MuPlusDirection(PmuPlusX / PmuPlus, PmuPlusY / PmuPlus, PmuPlusZ / PmuPlus); 241 G4ThreeVector MuMinusDirection(PmuMinusX / PmuMinus, PmuMinusY / PmuMinus, PmuMinusZ / PmuMinus); 242 243 // rotate to actual Positron direction 244 // 245 MuPlusDirection.rotateUz(PosiDirection); 246 MuMinusDirection.rotateUz(PosiDirection); 247 248 aParticleChange.SetNumberOfSecondaries(2); 249 250 // create G4DynamicParticle object for the particle1 251 auto aParticle1 = new G4DynamicParticle(part1, MuPlusDirection, EmuPlus - fMass); 252 aParticleChange.AddSecondary(aParticle1); 253 // create G4DynamicParticle object for the particle2 254 auto aParticle2 = new G4DynamicParticle(part2, MuMinusDirection, EmuMinus - fMass); 255 aParticleChange.AddSecondary(aParticle2); 256 257 // Kill the incident positron 258 // 259 aParticleChange.ProposeEnergy(0.); 260 aParticleChange.ProposeTrackStatus(fStopAndKill); 261 262 return &aParticleChange; 263 } 264 265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 266 267 void G4AnnihiToMuPair::PrintInfoDefinition() 268 { 269 G4String comments = fInfo + " annihilation, atomic e- at rest, SubType="; 270 G4cout << G4endl << GetProcessName() << ": " << comments << GetProcessSubType() << G4endl; 271 G4cout << " threshold at " << fLowEnergyLimit / CLHEP::GeV << " GeV" 272 << " good description up to " << fHighEnergyLimit / CLHEP::TeV << " TeV for all Z." 273 << G4endl; 274 } 275 276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 277