ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AnnihiToMuPair.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
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 "G4ios.hh"
44 #include "Randomize.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 
48 #include "G4Positron.hh"
49 #include "G4MuonPlus.hh"
50 #include "G4MuonMinus.hh"
51 #include "G4Material.hh"
52 #include "G4Step.hh"
53 #include "G4LossTableManager.hh"
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
57 using namespace std;
58 
59 G4AnnihiToMuPair::G4AnnihiToMuPair(const G4String& processName,
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.;
71  SetProcessSubType(6);
72  G4LossTableManager::Instance()->Register(this);
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
76 
77 G4AnnihiToMuPair::~G4AnnihiToMuPair() // (empty) destructor
78 {
79  G4LossTableManager::Instance()->DeRegister(this);
80 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83 
84 G4bool G4AnnihiToMuPair::IsApplicable(const G4ParticleDefinition& particle)
85 {
86  return ( &particle == G4Positron::Positron() );
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 void G4AnnihiToMuPair::BuildPhysicsTable(const G4ParticleDefinition&)
92 // Build cross section and mean free path tables
93 //here no tables, just calling PrintInfoDefinition
94 {
95  CurrentSigma = 0.0;
96  PrintInfoDefinition();
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 
101 void G4AnnihiToMuPair::SetCrossSecFactor(G4double fac)
102 // Set the factor to artificially increase the cross section
103 {
104  CrossSecFactor = fac;
105  G4cout << "The cross section for AnnihiToMuPair is artificially "
106  << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 
111 G4double G4AnnihiToMuPair::ComputeCrossSectionPerAtom(G4double Epos, G4double Z)
112 // Calculates the microscopic cross section in GEANT4 internal units.
113 // It gives a good description from threshold to 1000 GeV
114 {
115  static const G4double Mmuon = G4MuonPlus::MuonPlus()->GetPDGMass();
116  static const G4double Rmuon = CLHEP::elm_coupling/Mmuon; //classical particle radius
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
119 
120  G4double CrossSection = 0.;
121  if (Epos < LowestEnergyLimit) return CrossSection;
122 
123  G4double xi = LowestEnergyLimit/Epos;
124  G4double piaxi = pia * sqrt(xi);
125  G4double SigmaEl = Sig0 * xi * (1.+xi/2.) * piaxi;
126  if( Epos>LowestEnergyLimit+1.e-5 ) SigmaEl /= (1.-std::exp( -piaxi/std::sqrt(1-xi) ));
127  CrossSection = SigmaEl*Z; // SigmaEl per electron * number of electrons per atom
128  return CrossSection;
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 
133 G4double G4AnnihiToMuPair::CrossSectionPerVolume(G4double PositronEnergy,
134  const G4Material* aMaterial)
135 {
136  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
137  const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
138 
139  G4double SIGMA = 0.0;
140 
141  for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; ++i )
142  {
143  G4double AtomicZ = (*theElementVector)[i]->GetZ();
144  SIGMA += NbOfAtomsPerVolume[i] *
145  ComputeCrossSectionPerAtom(PositronEnergy,AtomicZ);
146  }
147  return SIGMA;
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
151 
152 G4double G4AnnihiToMuPair::GetMeanFreePath(const G4Track& aTrack,
153  G4double, G4ForceCondition*)
154 
155 // returns the positron mean free path in GEANT4 internal units
156 
157 {
158  const G4DynamicParticle* aDynamicPositron = aTrack.GetDynamicParticle();
159  G4double PositronEnergy = aDynamicPositron->GetKineticEnergy()
160  +electron_mass_c2;
161  G4Material* aMaterial = aTrack.GetMaterial();
162  CurrentSigma = CrossSectionPerVolume(PositronEnergy, aMaterial);
163 
164  // increase the CrossSection by CrossSecFactor (default 1)
165  G4double mfp = DBL_MAX;
166  if(CurrentSigma > DBL_MIN) mfp = 1.0/(CurrentSigma*CrossSecFactor);
167 
168  return mfp;
169 }
170 
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172 
173 G4VParticleChange* G4AnnihiToMuPair::PostStepDoIt(const G4Track& aTrack,
174  const G4Step& aStep)
175 //
176 // generation of e+e- -> mu+mu-
177 //
178 {
179 
180  aParticleChange.Initialize(aTrack);
181  static const G4double Mele=electron_mass_c2;
182  static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass();
183 
184  // current Positron energy and direction, return if energy too low
185  const G4DynamicParticle *aDynamicPositron = aTrack.GetDynamicParticle();
186  G4double Epos = aDynamicPositron->GetKineticEnergy() + Mele;
187 
188  // test of cross section
189  if(CurrentSigma*G4UniformRand() >
190  CrossSectionPerVolume(Epos, aTrack.GetMaterial()))
191  {
192  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
193  }
194 
195  if (Epos < LowestEnergyLimit) {
196  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
197  }
198 
199  G4ParticleMomentum PositronDirection =
200  aDynamicPositron->GetMomentumDirection();
201  G4double xi = LowestEnergyLimit/Epos; // xi is always less than 1,
202  // goes to 0 at high Epos
203 
204  // generate cost
205  //
206  G4double cost;
207  do { cost = 2.*G4UniformRand()-1.; }
208  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
209  while (2.*G4UniformRand() > 1.+xi+cost*cost*(1.-xi) );
210  //1+cost**2 at high Epos
211  G4double sint = sqrt(1.-cost*cost);
212 
213  // generate phi
214  //
215  G4double phi=2.*pi*G4UniformRand();
216 
217  G4double Ecm = sqrt(0.5*Mele*(Epos+Mele));
218  G4double Pcm = sqrt(Ecm*Ecm-Mmuon*Mmuon);
219  G4double beta = sqrt((Epos-Mele)/(Epos+Mele));
220  G4double gamma = Ecm/Mele; // =sqrt((Epos+Mele)/(2.*Mele));
221  G4double Pt = Pcm*sint;
222 
223  // energy and momentum of the muons in the Lab
224  //
225  G4double EmuPlus = gamma*( Ecm+cost*beta*Pcm);
226  G4double EmuMinus = gamma*( Ecm-cost*beta*Pcm);
227  G4double PmuPlusZ = gamma*(beta*Ecm+cost* Pcm);
228  G4double PmuMinusZ = gamma*(beta*Ecm-cost* Pcm);
229  G4double PmuPlusX = Pt*cos(phi);
230  G4double PmuPlusY = Pt*sin(phi);
231  G4double PmuMinusX =-Pt*cos(phi);
232  G4double PmuMinusY =-Pt*sin(phi);
233  // absolute momenta
234  G4double PmuPlus = sqrt(Pt*Pt+PmuPlusZ *PmuPlusZ );
235  G4double PmuMinus = sqrt(Pt*Pt+PmuMinusZ*PmuMinusZ);
236 
237  // mu+ mu- directions for Positron in z-direction
238  //
239  G4ThreeVector
240  MuPlusDirection ( PmuPlusX/PmuPlus, PmuPlusY/PmuPlus, PmuPlusZ/PmuPlus );
241  G4ThreeVector
242  MuMinusDirection(PmuMinusX/PmuMinus,PmuMinusY/PmuMinus,PmuMinusZ/PmuMinus);
243 
244  // rotate to actual Positron direction
245  //
246  MuPlusDirection.rotateUz(PositronDirection);
247  MuMinusDirection.rotateUz(PositronDirection);
248 
249  aParticleChange.SetNumberOfSecondaries(2);
250  // create G4DynamicParticle object for the particle1
251  G4DynamicParticle* aParticle1= new G4DynamicParticle(
252  G4MuonPlus::MuonPlus(),MuPlusDirection,EmuPlus-Mmuon);
253  aParticleChange.AddSecondary(aParticle1);
254  // create G4DynamicParticle object for the particle2
255  G4DynamicParticle* aParticle2= new G4DynamicParticle(
256  G4MuonMinus::MuonMinus(),MuMinusDirection,EmuMinus-Mmuon);
257  aParticleChange.AddSecondary(aParticle2);
258 
259  // Kill the incident positron
260  //
261  aParticleChange.ProposeEnergy(0.);
262  aParticleChange.ProposeTrackStatus(fStopAndKill);
263 
264  return &aParticleChange;
265 }
266 
267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
268 
269 void G4AnnihiToMuPair::PrintInfoDefinition()
270 {
271  G4String comments ="e+e->mu+mu- annihilation, atomic e- at rest, SubType=.";
272  G4cout << G4endl << GetProcessName() << ": " << comments
273  << GetProcessSubType() << G4endl;
274  G4cout << " threshold at " << LowestEnergyLimit/GeV << " GeV"
275  << " good description up to "
276  << HighestEnergyLimit/TeV << " TeV for all Z." << G4endl;
277 }
278 
279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......