Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4AnnihiToMuPair.cc
이 파일의 문서화 페이지로 가기
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 // $Id: G4AnnihiToMuPair.cc 108750 2018-03-02 15:26:50Z gcosmo $
28 //
29 // ------------ G4AnnihiToMuPair physics process ------
30 // by H.Burkhardt, S. Kelner and R. Kokoulin, November 2002
31 // -----------------------------------------------------------------------------
32 //
33 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......//
34 //
35 // 27.01.03 : first implementation (hbu)
36 // 04.02.03 : cosmetic simplifications (mma)
37 // 25.10.04 : migrade to new interfaces of ParticleChange (vi)
38 // 28.02.18 : cross section now including SSS threshold factor
39 //
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
41 
42 #include "G4AnnihiToMuPair.hh"
43 
44 #include "G4ios.hh"
45 #include "Randomize.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 
49 #include "G4Positron.hh"
50 #include "G4MuonPlus.hh"
51 #include "G4MuonMinus.hh"
52 #include "G4Material.hh"
53 #include "G4Step.hh"
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
57 using namespace std;
58 
60  G4ProcessType type):G4VDiscreteProcess (processName, type)
61 {
62  //e+ Energy threshold
63  const G4double Mu_massc2 = G4MuonPlus::MuonPlus()->GetPDGMass();
64  LowestEnergyLimit = 2.*Mu_massc2*Mu_massc2/electron_mass_c2 - electron_mass_c2;
65 
66  //modele ok up to 1000 TeV due to neglected Z-interference
67  HighestEnergyLimit = 1000.*TeV;
68 
69  CurrentSigma = 0.0;
70  CrossSecFactor = 1.;
72 
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
76 
77 G4AnnihiToMuPair::~G4AnnihiToMuPair() // (empty) destructor
78 { }
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81 
83 {
84  return ( &particle == G4Positron::Positron() );
85 }
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 
90 // Build cross section and mean free path tables
91 //here no tables, just calling PrintInfoDefinition
92 {
93  CurrentSigma = 0.0;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 
100 // Set the factor to artificially increase the cross section
101 {
103  G4cout << "The cross section for AnnihiToMuPair is artificially "
104  << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108 
110 // Calculates the microscopic cross section in GEANT4 internal units.
111 // It gives a good description from threshold to 1000 GeV
112 {
113  static const G4double Mmuon = G4MuonPlus::MuonPlus()->GetPDGMass();
114  static const G4double Rmuon = CLHEP::elm_coupling/Mmuon; //classical particle radius
115  static const G4double Sig0 = CLHEP::pi*Rmuon*Rmuon/3.; //constant in crossSection
116  static const G4double pia = CLHEP::pi * CLHEP::fine_structure_const; // pi * alphaQED
117 
118  G4double CrossSection = 0.;
119  if (Epos < LowestEnergyLimit) return CrossSection;
120 
121  G4double xi = LowestEnergyLimit/Epos;
122  G4double piaxi = pia * sqrt(xi);
123  G4double SigmaEl = Sig0 * xi * (1.+xi/2.) * piaxi;
124  if( Epos>LowestEnergyLimit+1.e-5 ) SigmaEl /= (1.-std::exp( -piaxi/std::sqrt(1-xi) ));
125  CrossSection = SigmaEl*Z; // SigmaEl per electron * number of electrons per atom
126  return CrossSection;
127 }
128 
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130 
132  const G4Material* aMaterial)
133 {
134  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
135  const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
136 
137  G4double SIGMA = 0.0;
138 
139  for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; ++i )
140  {
141  G4double AtomicZ = (*theElementVector)[i]->GetZ();
142  SIGMA += NbOfAtomsPerVolume[i] *
143  ComputeCrossSectionPerAtom(PositronEnergy,AtomicZ);
144  }
145  return SIGMA;
146 }
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149 
152 
153 // returns the positron mean free path in GEANT4 internal units
154 
155 {
156  const G4DynamicParticle* aDynamicPositron = aTrack.GetDynamicParticle();
157  G4double PositronEnergy = aDynamicPositron->GetKineticEnergy()
159  G4Material* aMaterial = aTrack.GetMaterial();
160  CurrentSigma = CrossSectionPerVolume(PositronEnergy, aMaterial);
161 
162  // increase the CrossSection by CrossSecFactor (default 1)
163  G4double mfp = DBL_MAX;
165 
166  return mfp;
167 }
168 
169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
170 
172  const G4Step& aStep)
173 //
174 // generation of e+e- -> mu+mu-
175 //
176 {
177 
178  aParticleChange.Initialize(aTrack);
179  static const G4double Mele=electron_mass_c2;
180  static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass();
181 
182  // current Positron energy and direction, return if energy too low
183  const G4DynamicParticle *aDynamicPositron = aTrack.GetDynamicParticle();
184  G4double Epos = aDynamicPositron->GetKineticEnergy() + Mele;
185 
186  // test of cross section
188  CrossSectionPerVolume(Epos, aTrack.GetMaterial()))
189  {
190  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
191  }
192 
193  if (Epos < LowestEnergyLimit) {
194  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
195  }
196 
197  G4ParticleMomentum PositronDirection =
198  aDynamicPositron->GetMomentumDirection();
199  G4double xi = LowestEnergyLimit/Epos; // xi is always less than 1,
200  // goes to 0 at high Epos
201 
202  // generate cost
203  //
204  G4double cost;
205  do { cost = 2.*G4UniformRand()-1.; }
206  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
207  while (2.*G4UniformRand() > 1.+xi+cost*cost*(1.-xi) );
208  //1+cost**2 at high Epos
209  G4double sint = sqrt(1.-cost*cost);
210 
211  // generate phi
212  //
213  G4double phi=2.*pi*G4UniformRand();
214 
215  G4double Ecm = sqrt(0.5*Mele*(Epos+Mele));
216  G4double Pcm = sqrt(Ecm*Ecm-Mmuon*Mmuon);
217  G4double beta = sqrt((Epos-Mele)/(Epos+Mele));
218  G4double gamma = Ecm/Mele; // =sqrt((Epos+Mele)/(2.*Mele));
219  G4double Pt = Pcm*sint;
220 
221  // energy and momentum of the muons in the Lab
222  //
223  G4double EmuPlus = gamma*( Ecm+cost*beta*Pcm);
224  G4double EmuMinus = gamma*( Ecm-cost*beta*Pcm);
225  G4double PmuPlusZ = gamma*(beta*Ecm+cost* Pcm);
226  G4double PmuMinusZ = gamma*(beta*Ecm-cost* Pcm);
227  G4double PmuPlusX = Pt*cos(phi);
228  G4double PmuPlusY = Pt*sin(phi);
229  G4double PmuMinusX =-Pt*cos(phi);
230  G4double PmuMinusY =-Pt*sin(phi);
231  // absolute momenta
232  G4double PmuPlus = sqrt(Pt*Pt+PmuPlusZ *PmuPlusZ );
233  G4double PmuMinus = sqrt(Pt*Pt+PmuMinusZ*PmuMinusZ);
234 
235  // mu+ mu- directions for Positron in z-direction
236  //
238  MuPlusDirection ( PmuPlusX/PmuPlus, PmuPlusY/PmuPlus, PmuPlusZ/PmuPlus );
240  MuMinusDirection(PmuMinusX/PmuMinus,PmuMinusY/PmuMinus,PmuMinusZ/PmuMinus);
241 
242  // rotate to actual Positron direction
243  //
244  MuPlusDirection.rotateUz(PositronDirection);
245  MuMinusDirection.rotateUz(PositronDirection);
246 
248  // create G4DynamicParticle object for the particle1
249  G4DynamicParticle* aParticle1= new G4DynamicParticle(
250  G4MuonPlus::MuonPlus(),MuPlusDirection,EmuPlus-Mmuon);
251  aParticleChange.AddSecondary(aParticle1);
252  // create G4DynamicParticle object for the particle2
253  G4DynamicParticle* aParticle2= new G4DynamicParticle(
254  G4MuonMinus::MuonMinus(),MuMinusDirection,EmuMinus-Mmuon);
255  aParticleChange.AddSecondary(aParticle2);
256 
257  // Kill the incident positron
258  //
261 
262  return &aParticleChange;
263 }
264 
265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
266 
268 {
269  G4String comments ="e+e->mu+mu- annihilation, atomic e- at rest, SubType=.";
270  G4cout << G4endl << GetProcessName() << ": " << comments
271  << GetProcessSubType() << G4endl;
272  G4cout << " threshold at " << LowestEnergyLimit/GeV << " GeV"
273  << " good description up to "
274  << HighestEnergyLimit/TeV << " TeV for all Z." << G4endl;
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *) override
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & GetMomentumDirection() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:191
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
static G4MuonMinus * MuonMinus()
Definition: G4MuonMinus.cc:100
Double_t beta
G4double GetPDGMass() const
void AddSecondary(G4Track *aSecondary)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
Float_t Z
static constexpr double TeV
Definition: G4SIunits.hh:218
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
void BuildPhysicsTable(const G4ParticleDefinition &) override
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep) override
virtual void Initialize(const G4Track &)
G4ProcessType
static constexpr double electron_mass_c2
static constexpr double elm_coupling
G4double CrossSectionPerVolume(G4double PositronEnergy, const G4Material *)
#define G4UniformRand()
Definition: Randomize.hh:53
const G4String & GetProcessName() const
Definition: G4VProcess.hh:411
static G4Positron * Positron()
Definition: G4Positron.cc:94
Definition: G4Step.hh:76
void ProposeEnergy(G4double finalEnergy)
G4Material * GetMaterial() const
G4double ComputeCrossSectionPerAtom(G4double PositronEnergy, G4double AtomicZ)
G4bool IsApplicable(const G4ParticleDefinition &) override
void SetCrossSecFactor(G4double fac)
std::vector< G4Element * > G4ElementVector
static const G4double fac
#define DBL_MIN
Definition: templates.hh:75
static G4MuonPlus * MuonPlus()
Definition: G4MuonPlus.cc:99
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:435
G4ForceCondition
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
static constexpr double pi
Definition: G4SIunits.hh:75
static constexpr double fine_structure_const
#define DBL_MAX
Definition: templates.hh:83
G4AnnihiToMuPair(const G4String &processName="AnnihiToMuPair", G4ProcessType type=fElectromagnetic)
void ProposeTrackStatus(G4TrackStatus status)
static constexpr double GeV
Definition: G4SIunits.hh:217
const G4DynamicParticle * GetDynamicParticle() const
static constexpr double pi
Definition: SystemOfUnits.h:54
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
G4int GetProcessSubType() const
Definition: G4VProcess.hh:429
void SetNumberOfSecondaries(G4int totSecondaries)