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 // ------------ G4AnnihiToMuPair physi 28 // ------------ G4AnnihiToMuPair physics process ------ 29 // by H.Burkhardt, S. Kelner and R. Ko 29 // by H.Burkhardt, S. Kelner and R. Kokoulin, November 2002 30 // ------------------------------------------- 30 // ----------------------------------------------------------------------------- 31 // 31 // 32 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......// 33 // 33 // 34 // 27.01.03 : first implementation (hbu) 34 // 27.01.03 : first implementation (hbu) 35 // 04.02.03 : cosmetic simplifications (mma) 35 // 04.02.03 : cosmetic simplifications (mma) 36 // 25.10.04 : migrade to new interfaces of Par 36 // 25.10.04 : migrade to new interfaces of ParticleChange (vi) 37 // 28.02.18 : cross section now including SSS 37 // 28.02.18 : cross section now including SSS threshold factor 38 // 38 // 39 //....oooOO0OOooo........oooOO0OOooo........oo 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 40 40 41 #include "G4AnnihiToMuPair.hh" 41 #include "G4AnnihiToMuPair.hh" 42 42 43 #include "G4Exp.hh" << 43 #include "G4ios.hh" 44 #include "G4LossTableManager.hh" << 44 #include "Randomize.hh" 45 #include "G4Material.hh" << 46 #include "G4MuonMinus.hh" << 47 #include "G4MuonPlus.hh" << 48 #include "G4PhysicalConstants.hh" 45 #include "G4PhysicalConstants.hh" >> 46 #include "G4SystemOfUnits.hh" >> 47 49 #include "G4Positron.hh" 48 #include "G4Positron.hh" >> 49 #include "G4MuonPlus.hh" >> 50 #include "G4MuonMinus.hh" >> 51 #include "G4Material.hh" 50 #include "G4Step.hh" 52 #include "G4Step.hh" 51 #include "G4SystemOfUnits.hh" << 53 #include "G4LossTableManager.hh" 52 #include "G4TauMinus.hh" << 53 #include "G4TauPlus.hh" << 54 #include "G4ios.hh" << 55 #include "Randomize.hh" << 56 54 57 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 58 56 >> 57 using namespace std; >> 58 59 G4AnnihiToMuPair::G4AnnihiToMuPair(const G4Str 59 G4AnnihiToMuPair::G4AnnihiToMuPair(const G4String& processName, 60 G4ProcessType type):G4VDiscreteProcess (pr 60 G4ProcessType type):G4VDiscreteProcess (processName, type) 61 { 61 { 62 //e+ Energy threshold 62 //e+ Energy threshold 63 if(processName == "AnnihiToTauPair") { << 63 const G4double Mu_massc2 = G4MuonPlus::MuonPlus()->GetPDGMass(); 64 SetProcessSubType(fAnnihilationToTauTau); << 64 LowestEnergyLimit = 2.*Mu_massc2*Mu_massc2/electron_mass_c2 - electron_mass_c2; 65 part1 = G4TauPlus::TauPlus(); << 65 66 part2 = G4TauMinus::TauMinus(); << 66 //modele ok up to 1000 TeV due to neglected Z-interference 67 fInfo = "e+e->tau+tau-"; << 67 HighestEnergyLimit = 1000.*TeV; 68 } else { << 68 69 SetProcessSubType(fAnnihilationToMuMu); << 69 CurrentSigma = 0.0; 70 part1 = G4MuonPlus::MuonPlus(); << 70 CrossSecFactor = 1.; 71 part2 = G4MuonMinus::MuonMinus(); << 71 SetProcessSubType(6); 72 } << 72 G4LossTableManager::Instance()->Register(this); 73 fMass = part1->GetPDGMass(); << 74 fLowEnergyLimit = 2. * fMass * fMass / CLHEP << 75 << 76 // model is ok up to 1000 TeV due to neglect << 77 fHighEnergyLimit = 1000. * TeV; << 78 << 79 fCurrentSigma = 0.0; << 80 fCrossSecFactor = 1.; << 81 fManager = G4LossTableManager::Instance(); << 82 fManager->Register(this); << 83 } 73 } 84 74 85 //....oooOO0OOooo........oooOO0OOooo........oo 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 86 76 87 G4AnnihiToMuPair::~G4AnnihiToMuPair() // (empt 77 G4AnnihiToMuPair::~G4AnnihiToMuPair() // (empty) destructor 88 { 78 { 89 fManager->DeRegister(this); << 79 G4LossTableManager::Instance()->DeRegister(this); 90 } 80 } 91 81 92 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 93 83 94 G4bool G4AnnihiToMuPair::IsApplicable(const G4 84 G4bool G4AnnihiToMuPair::IsApplicable(const G4ParticleDefinition& particle) 95 { 85 { 96 return ( &particle == G4Positron::Positron() 86 return ( &particle == G4Positron::Positron() ); 97 } 87 } 98 88 99 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 100 90 101 void G4AnnihiToMuPair::BuildPhysicsTable(const 91 void G4AnnihiToMuPair::BuildPhysicsTable(const G4ParticleDefinition&) >> 92 // Build cross section and mean free path tables >> 93 //here no tables, just calling PrintInfoDefinition 102 { 94 { >> 95 CurrentSigma = 0.0; 103 PrintInfoDefinition(); 96 PrintInfoDefinition(); 104 } 97 } 105 98 106 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 107 100 108 void G4AnnihiToMuPair::SetCrossSecFactor(G4dou 101 void G4AnnihiToMuPair::SetCrossSecFactor(G4double fac) 109 // Set the factor to artificially increase the 102 // Set the factor to artificially increase the cross section 110 { 103 { 111 fCrossSecFactor = fac; << 104 CrossSecFactor = fac; 112 //G4cout << "The cross section for AnnihiToM << 105 G4cout << "The cross section for AnnihiToMuPair is artificially " 113 // << "increased by the CrossSecFactor << 106 << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl; 114 } 107 } 115 108 116 //....oooOO0OOooo........oooOO0OOooo........oo 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 117 110 118 G4double G4AnnihiToMuPair::ComputeCrossSection << 111 G4double G4AnnihiToMuPair::ComputeCrossSectionPerAtom(G4double Epos, G4double Z) 119 // Calculates the microscopic cross section in 112 // Calculates the microscopic cross section in GEANT4 internal units. 120 // It gives a good description from threshold 113 // It gives a good description from threshold to 1000 GeV 121 { 114 { 122 G4double rmuon = CLHEP::elm_coupling/fMass; << 115 static const G4double Mmuon = G4MuonPlus::MuonPlus()->GetPDGMass(); 123 G4double sig0 = CLHEP::pi*rmuon*rmuon/3.; << 116 static const G4double Rmuon = CLHEP::elm_coupling/Mmuon; //classical particle radius 124 const G4double pial = CLHEP::pi*CLHEP::fine_ << 117 static const G4double Sig0 = CLHEP::pi*Rmuon*Rmuon/3.; //constant in crossSection >> 118 static const G4double pia = CLHEP::pi * CLHEP::fine_structure_const; // pi * alphaQED 125 119 126 if (e <= fLowEnergyLimit) return 0.0; << 120 G4double CrossSection = 0.; >> 121 if (Epos < LowestEnergyLimit) return CrossSection; 127 122 128 const G4double xi = fLowEnergyLimit/e; << 123 G4double xi = LowestEnergyLimit/Epos; 129 const G4double piaxi = pial * std::sqrt(xi); << 124 G4double piaxi = pia * sqrt(xi); 130 G4double sigma = sig0 * xi * (1. + xi*0.5); << 125 G4double SigmaEl = Sig0 * xi * (1.+xi/2.) * piaxi; 131 //G4cout << "### xi= " << xi << " piaxi=" << << 126 if( Epos>LowestEnergyLimit+1.e-5 ) SigmaEl /= (1.-std::exp( -piaxi/std::sqrt(1-xi) )); 132 << 127 CrossSection = SigmaEl*Z; // SigmaEl per electron * number of electrons per atom 133 // argument of the exponent below 0.1 or abo << 128 return CrossSection; 134 // Sigma per electron * number of electrons << 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: << 143 } << 144 // G4cout << "### sigma= " << sigma << G4end << 145 return sigma; << 146 } 129 } 147 130 148 //....oooOO0OOooo........oooOO0OOooo........oo 131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 149 132 150 G4double G4AnnihiToMuPair::ComputeCrossSection << 133 G4double G4AnnihiToMuPair::CrossSectionPerVolume(G4double PositronEnergy, 151 << 134 const G4Material* aMaterial) 152 { 135 { 153 return ComputeCrossSectionPerElectron(energy << 136 const G4ElementVector* theElementVector = aMaterial->GetElementVector(); 154 } << 137 const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume(); 155 138 156 //....oooOO0OOooo........oooOO0OOooo........oo << 139 G4double SIGMA = 0.0; 157 140 158 G4double G4AnnihiToMuPair::CrossSectionPerVolu << 141 for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; ++i ) 159 const G4Material* aMaterial) << 142 { 160 { << 143 G4double AtomicZ = (*theElementVector)[i]->GetZ(); 161 return ComputeCrossSectionPerElectron(energy << 144 SIGMA += NbOfAtomsPerVolume[i] * >> 145 ComputeCrossSectionPerAtom(PositronEnergy,AtomicZ); >> 146 } >> 147 return SIGMA; 162 } 148 } 163 149 164 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 165 151 166 G4double G4AnnihiToMuPair::GetMeanFreePath(con 152 G4double G4AnnihiToMuPair::GetMeanFreePath(const G4Track& aTrack, 167 G4d 153 G4double, G4ForceCondition*) >> 154 168 // returns the positron mean free path in GEAN 155 // returns the positron mean free path in GEANT4 internal units >> 156 169 { 157 { 170 const G4DynamicParticle* aDynamicPositron = 158 const G4DynamicParticle* aDynamicPositron = aTrack.GetDynamicParticle(); 171 G4double energy = aDynamicPositron->GetTotal << 159 G4double PositronEnergy = aDynamicPositron->GetKineticEnergy() 172 const G4Material* aMaterial = aTrack.GetMate << 160 +electron_mass_c2; 173 << 161 G4Material* aMaterial = aTrack.GetMaterial(); 174 // cross section before step << 162 CurrentSigma = CrossSectionPerVolume(PositronEnergy, aMaterial); 175 fCurrentSigma = CrossSectionPerVolume(energy << 176 163 177 // increase the CrossSection by CrossSecFact 164 // increase the CrossSection by CrossSecFactor (default 1) 178 return (fCurrentSigma > 0.0) ? 1.0/(fCurrent << 165 G4double mfp = DBL_MAX; >> 166 if(CurrentSigma > DBL_MIN) mfp = 1.0/(CurrentSigma*CrossSecFactor); >> 167 >> 168 return mfp; 179 } 169 } 180 170 181 //....oooOO0OOooo........oooOO0OOooo........oo 171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 182 172 183 G4VParticleChange* G4AnnihiToMuPair::PostStepD 173 G4VParticleChange* G4AnnihiToMuPair::PostStepDoIt(const G4Track& aTrack, 184 174 const G4Step& aStep) 185 // 175 // 186 // generation of e+e- -> mu+mu- 176 // generation of e+e- -> mu+mu- 187 // 177 // 188 { 178 { >> 179 189 aParticleChange.Initialize(aTrack); 180 aParticleChange.Initialize(aTrack); >> 181 static const G4double Mele=electron_mass_c2; >> 182 static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass(); 190 183 191 // current Positron energy and direction, re 184 // current Positron energy and direction, return if energy too low 192 const G4DynamicParticle *aDynamicPositron = 185 const G4DynamicParticle *aDynamicPositron = aTrack.GetDynamicParticle(); 193 const G4double Mele = CLHEP::electron_mass_c << 186 G4double Epos = aDynamicPositron->GetKineticEnergy() + Mele; 194 G4double Epos = aDynamicPositron->GetTotalEn << 195 G4double xs = CrossSectionPerVolume(Epos, aT << 196 187 197 // test of cross section 188 // test of cross section 198 if(xs > 0.0 && fCurrentSigma*G4UniformRand() << 189 if(CurrentSigma*G4UniformRand() > 199 return G4VDiscreteProcess::PostStepDoIt(aT << 190 CrossSectionPerVolume(Epos, aTrack.GetMaterial())) >> 191 { >> 192 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep); >> 193 } >> 194 >> 195 if (Epos < LowestEnergyLimit) { >> 196 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep); 200 } 197 } 201 198 202 const G4ThreeVector PosiDirection = aDynamic << 199 G4ParticleMomentum PositronDirection = 203 G4double xi = fLowEnergyLimit/Epos; // xi is << 200 aDynamicPositron->GetMomentumDirection(); 204 // goes << 201 G4double xi = LowestEnergyLimit/Epos; // xi is always less than 1, >> 202 // goes to 0 at high Epos 205 203 206 // generate cost; probability function 1+cos << 204 // generate cost 207 // 205 // 208 G4double cost; 206 G4double cost; 209 do { cost = 2.*G4UniformRand()-1.; } 207 do { cost = 2.*G4UniformRand()-1.; } 210 // Loop checking, 07-Aug-2015, Vladimir Ivan 208 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko 211 while (2.*G4UniformRand() > 1.+xi+cost*cost* 209 while (2.*G4UniformRand() > 1.+xi+cost*cost*(1.-xi) ); 212 G4double sint = std::sqrt(1.-cost*cost); << 210 //1+cost**2 at high Epos >> 211 G4double sint = sqrt(1.-cost*cost); 213 212 214 // generate phi 213 // generate phi 215 // 214 // 216 G4double phi = 2.*CLHEP::pi*G4UniformRand(); << 215 G4double phi=2.*pi*G4UniformRand(); 217 216 218 G4double Ecm = std::sqrt(0.5*Mele*(Epos+Me << 217 G4double Ecm = sqrt(0.5*Mele*(Epos+Mele)); 219 G4double Pcm = std::sqrt(Ecm*Ecm - fMass*f << 218 G4double Pcm = sqrt(Ecm*Ecm-Mmuon*Mmuon); 220 G4double beta = std::sqrt((Epos-Mele)/(Epos << 219 G4double beta = sqrt((Epos-Mele)/(Epos+Mele)); 221 G4double gamma = Ecm/Mele; << 220 G4double gamma = Ecm/Mele; // =sqrt((Epos+Mele)/(2.*Mele)); 222 G4double Pt = Pcm*sint; 221 G4double Pt = Pcm*sint; 223 222 224 // energy and momentum of the muons in the L 223 // energy and momentum of the muons in the Lab 225 // 224 // 226 G4double EmuPlus = gamma*(Ecm + cost*beta* << 225 G4double EmuPlus = gamma*( Ecm+cost*beta*Pcm); 227 G4double EmuMinus = gamma*(Ecm - cost*beta* << 226 G4double EmuMinus = gamma*( Ecm-cost*beta*Pcm); 228 G4double PmuPlusZ = gamma*(beta*Ecm + cost* << 227 G4double PmuPlusZ = gamma*(beta*Ecm+cost* Pcm); 229 G4double PmuMinusZ = gamma*(beta*Ecm - cost* << 228 G4double PmuMinusZ = gamma*(beta*Ecm-cost* Pcm); 230 G4double PmuPlusX = Pt*std::cos(phi); << 229 G4double PmuPlusX = Pt*cos(phi); 231 G4double PmuPlusY = Pt*std::sin(phi); << 230 G4double PmuPlusY = Pt*sin(phi); 232 G4double PmuMinusX =-PmuPlusX; << 231 G4double PmuMinusX =-Pt*cos(phi); 233 G4double PmuMinusY =-PmuPlusY; << 232 G4double PmuMinusY =-Pt*sin(phi); 234 // absolute momenta 233 // absolute momenta 235 G4double PmuPlus = std::sqrt(Pt*Pt+PmuPlusZ << 234 G4double PmuPlus = sqrt(Pt*Pt+PmuPlusZ *PmuPlusZ ); 236 G4double PmuMinus = std::sqrt(Pt*Pt+PmuMinus << 235 G4double PmuMinus = sqrt(Pt*Pt+PmuMinusZ*PmuMinusZ); 237 236 238 // mu+ mu- directions for Positron in z-dire 237 // mu+ mu- directions for Positron in z-direction 239 // 238 // 240 G4ThreeVector MuPlusDirection(PmuPlusX / Pmu << 239 G4ThreeVector 241 G4ThreeVector MuMinusDirection(PmuMinusX / P << 240 MuPlusDirection ( PmuPlusX/PmuPlus, PmuPlusY/PmuPlus, PmuPlusZ/PmuPlus ); >> 241 G4ThreeVector >> 242 MuMinusDirection(PmuMinusX/PmuMinus,PmuMinusY/PmuMinus,PmuMinusZ/PmuMinus); 242 243 243 // rotate to actual Positron direction 244 // rotate to actual Positron direction 244 // 245 // 245 MuPlusDirection.rotateUz(PosiDirection); << 246 MuPlusDirection.rotateUz(PositronDirection); 246 MuMinusDirection.rotateUz(PosiDirection); << 247 MuMinusDirection.rotateUz(PositronDirection); 247 248 248 aParticleChange.SetNumberOfSecondaries(2); 249 aParticleChange.SetNumberOfSecondaries(2); 249 << 250 // create G4DynamicParticle object for the p 250 // create G4DynamicParticle object for the particle1 251 auto aParticle1 = new G4DynamicParticle(part << 251 G4DynamicParticle* aParticle1= new G4DynamicParticle( >> 252 G4MuonPlus::MuonPlus(),MuPlusDirection,EmuPlus-Mmuon); 252 aParticleChange.AddSecondary(aParticle1); 253 aParticleChange.AddSecondary(aParticle1); 253 // create G4DynamicParticle object for the p 254 // create G4DynamicParticle object for the particle2 254 auto aParticle2 = new G4DynamicParticle(part << 255 G4DynamicParticle* aParticle2= new G4DynamicParticle( >> 256 G4MuonMinus::MuonMinus(),MuMinusDirection,EmuMinus-Mmuon); 255 aParticleChange.AddSecondary(aParticle2); 257 aParticleChange.AddSecondary(aParticle2); 256 258 257 // Kill the incident positron 259 // Kill the incident positron 258 // 260 // 259 aParticleChange.ProposeEnergy(0.); 261 aParticleChange.ProposeEnergy(0.); 260 aParticleChange.ProposeTrackStatus(fStopAndK 262 aParticleChange.ProposeTrackStatus(fStopAndKill); 261 263 262 return &aParticleChange; 264 return &aParticleChange; 263 } 265 } 264 266 265 //....oooOO0OOooo........oooOO0OOooo........oo 267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 266 268 267 void G4AnnihiToMuPair::PrintInfoDefinition() 269 void G4AnnihiToMuPair::PrintInfoDefinition() 268 { 270 { 269 G4String comments = fInfo + " annihilation, << 271 G4String comments ="e+e->mu+mu- annihilation, atomic e- at rest, SubType=."; 270 G4cout << G4endl << GetProcessName() << ": << 272 G4cout << G4endl << GetProcessName() << ": " << comments 271 G4cout << " threshold at " << fLowEne << 273 << GetProcessSubType() << G4endl; 272 << " good description up to " << fHig << 274 G4cout << " threshold at " << LowestEnergyLimit/GeV << " GeV" 273 << G4endl; << 275 << " good description up to " >> 276 << HighestEnergyLimit/TeV << " TeV for all Z." << G4endl; 274 } 277 } 275 278 276 //....oooOO0OOooo........oooOO0OOooo........oo 279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 277 280