Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4PhononDownconversion.cc
Go to the documentation of this file.
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/// \file processes/phonon/src/G4PhononDownconversion.cc
27/// \brief Implementation of the G4PhononDownconversion class
28//
29//
30// 20131111 Add verbose output for MFP calculation
31// 20131115 Initialize data buffers in ctor
32
34#include "G4LatticePhysical.hh"
35#include "G4PhononLong.hh"
37#include "G4PhononTrackMap.hh"
38#include "G4PhononTransFast.hh"
39#include "G4PhononTransSlow.hh"
41#include "G4RandomDirection.hh"
42#include "G4Step.hh"
43#include "G4SystemOfUnits.hh"
44#include "G4VParticleChange.hh"
45#include "Randomize.hh"
46#include <cmath>
47
48
49
51 : G4VPhononProcess(aName), fBeta(0.), fGamma(0.), fLambda(0.), fMu(0.) {;}
52
54
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56
58 G4double /*previousStepSize*/,
60 //Determines mean free path for longitudinal phonons to split
62 G4double Eoverh = aTrack.GetKineticEnergy()/h_Planck;
63
64 //Calculate mean free path for anh. decay
65 G4double mfp = aTrack.GetVelocity()/(Eoverh*Eoverh*Eoverh*Eoverh*Eoverh*A);
66
67 if (verboseLevel > 1)
68 G4cout << "G4PhononDownconversion::GetMeanFreePath = " << mfp << G4endl;
69
71 return mfp;
72}
73
74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
76
78 const G4Step&) {
80
81 //Obtain dynamical constants from this volume's lattice
82 fBeta=theLattice->GetBeta();
83 fGamma=theLattice->GetGamma();
84 fLambda=theLattice->GetLambda();
85 fMu=theLattice->GetMu();
86
87 //Destroy the parent phonon and create the daughter phonons.
88 //74% chance that daughter phonons are both transverse
89 //26% Transverse and Longitudinal
90 if (G4UniformRand()>0.740) MakeLTSecondaries(aTrack);
91 else MakeTTSecondaries(aTrack);
92
95
96 return &aParticleChange;
97}
98
99//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100
102 //Only L-phonons decay
103 return (&aPD==G4PhononLong::PhononDefinition());
104}
105
106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107
108//probability density of energy distribution of L'-phonon in L->L'+T process
109
110G4double G4PhononDownconversion::GetLTDecayProb(G4double d, G4double x) const {
111 //d=delta= ratio of group velocities vl/vt and x is the fraction of energy in the longitudinal mode, i.e. x=EL'/EL
112 return (1/(x*x))*(1-x*x)*(1-x*x)*((1+x)*(1+x)-d*d*((1-x)*(1-x)))*(1+x*x-d*d*(1-x)*(1-x))*(1+x*x-d*d*(1-x)*(1-x));
113}
114
115//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116
117//probability density of energy distribution of T-phonon in L->T+T process
118G4double G4PhononDownconversion::GetTTDecayProb(G4double d, G4double x) const {
119 //dynamic constants from Tamura, PRL31, 1985
120 G4double A = 0.5*(1-d*d)*(fBeta+fLambda+(1+d*d)*(fGamma+fMu));
121 G4double B = fBeta+fLambda+2*d*d*(fGamma+fMu);
122 G4double C = fBeta + fLambda + 2*(fGamma+fMu);
123 G4double D = (1-d*d)*(2*fBeta+4*fGamma+fLambda+3*fMu);
124
125 return (A+B*d*x-B*x*x)*(A+B*d*x-B*x*x)+(C*x*(d-x)-D/(d-x)*(x-d-(1-d*d)/(4*x)))*(C*x*(d-x)-D/(d-x)*(x-d-(1-d*d)/(4*x)));
126}
127
128//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129
130
131G4double G4PhononDownconversion::MakeLDeviation(G4double d, G4double x) const {
132 //change in L'-phonon propagation direction after decay
133
134 return std::acos((1+(x*x)-((d*d)*(1-x)*(1-x)))/(2*x));
135}
136
137//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138
139
140G4double G4PhononDownconversion::MakeTDeviation(G4double d, G4double x) const {
141 //change in T-phonon propagation direction after decay (L->L+T process)
142
143 return std::acos((1-x*x+d*d*(1-x)*(1-x))/(2*d*(1-x)));
144}
145
146//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147
148
149G4double G4PhononDownconversion::MakeTTDeviation(G4double d, G4double x) const {
150 //change in T-phonon propagation direction after decay (L->T+T process)
151
152 return std::acos((1-d*d*(1-x)*(1-x)+d*d*x*x)/(2*d*x));
153}
154
155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156
157
158//Generate daughter phonons from L->T+T process
159
160void G4PhononDownconversion::MakeTTSecondaries(const G4Track& aTrack) {
161 //d is the velocity ratio vL/vT
162 G4double d=1.6338;
163 G4double upperBound=(1+(1/d))/2;
164 G4double lowerBound=(1-(1/d))/2;
165
166 //Use MC method to generate point from distribution:
167 //if a random point on the energy-probability plane is
168 //smaller that the curve of the probability density,
169 //then accept that point.
170 //x=fraction of parent phonon energy in first T phonon
171 G4double x = G4UniformRand()*(upperBound-lowerBound) + lowerBound;
172 G4double p = 1.5*G4UniformRand();
173 while(p >= GetTTDecayProb(d, x*d)) {
174 x = G4UniformRand()*(upperBound-lowerBound) + lowerBound;
175 p = 1.5*G4UniformRand();
176 }
177
178 //using energy fraction x to calculate daughter phonon directions
179 G4double theta1=MakeTTDeviation(d, x);
180 G4double theta2=MakeTTDeviation(d, 1-x);
181 G4ThreeVector dir1=trackKmap->GetK(aTrack);
182 G4ThreeVector dir2=dir1;
183
184 // FIXME: These extra randoms change timing and causting outputs of example!
185 G4ThreeVector ran = G4RandomDirection(); // FIXME: Drop this line
186
187 G4double ph=G4UniformRand()*twopi;
188 dir1 = dir1.rotate(dir1.orthogonal(),theta1).rotate(dir1, ph);
189 dir2 = dir2.rotate(dir2.orthogonal(),-theta2).rotate(dir2,ph);
190
191 G4double E=aTrack.GetKineticEnergy();
192 G4double Esec1 = x*E, Esec2 = E-Esec1;
193
194 // Make FT or ST phonon (0. means no longitudinal)
195 G4int polarization1 = ChoosePolarization(0., theLattice->GetSTDOS(),
197
198 // Make FT or ST phonon (0. means no longitudinal)
199 G4int polarization2 = ChoosePolarization(0., theLattice->GetSTDOS(),
201
202 // Construct the secondaries and set their wavevectors
203 G4Track* sec1 = CreateSecondary(polarization1, dir1, Esec1);
204 G4Track* sec2 = CreateSecondary(polarization2, dir2, Esec2);
205
209}
210
211//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
212
213
214//Generate daughter phonons from L->L'+T process
215
216void G4PhononDownconversion::MakeLTSecondaries(const G4Track& aTrack) {
217 //d is the velocity ratio vL/v
218 G4double d=1.6338;
219 G4double upperBound=1;
220 G4double lowerBound=(d-1)/(d+1);
221
222 //Use MC method to generate point from distribution:
223 //if a random point on the energy-probability plane is
224 //smaller that the curve of the probability density,
225 //then accept that point.
226 //x=fraction of parent phonon energy in L phonon
227 G4double x = G4UniformRand()*(upperBound-lowerBound) + lowerBound;
228 G4double p = 4.0*G4UniformRand();
229 while(p >= GetLTDecayProb(d, x)) {
230 x = G4UniformRand()*(upperBound-lowerBound) + lowerBound;
231 p = 4.0*G4UniformRand(); //4.0 is about the max in the PDF
232 }
233
234 //using energy fraction x to calculate daughter phonon directions
235 G4double thetaL=MakeLDeviation(d, x);
236 G4double thetaT=MakeTDeviation(d, x); // FIXME: Should be 1-x?
237 G4ThreeVector dir1=trackKmap->GetK(aTrack);
238 G4ThreeVector dir2=dir1;
239
240 G4double ph=G4UniformRand()*twopi;
241 dir1 = dir1.rotate(dir1.orthogonal(),thetaL).rotate(dir1, ph);
242 dir2 = dir2.rotate(dir2.orthogonal(),-thetaT).rotate(dir2,ph);
243
244 G4double E=aTrack.GetKineticEnergy();
245 G4double Esec1 = x*E, Esec2 = E-Esec1;
246
247 // First secondary is longitudnal
248 G4int polarization1 = G4PhononPolarization::Long;
249
250 // Make FT or ST phonon (0. means no longitudinal)
251 G4int polarization2 = ChoosePolarization(0., theLattice->GetSTDOS(),
253
254 // Construct the secondaries and set their wavevectors
255 G4Track* sec1 = CreateSecondary(polarization1, dir1, Esec1);
256 G4Track* sec2 = CreateSecondary(polarization2, dir2, Esec2);
257
261}
262
263//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
264
double B(double temperature)
double C(double temp)
double A(double temperature)
double D(double temp)
G4double condition(const G4ErrorSymMatrix &m)
G4ForceCondition
@ NotForced
Definition of the G4LatticePhysical class.
G4ThreeVector G4RandomDirection()
@ fStopAndKill
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
Hep3Vector orthogonal() const
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:24
G4double GetFTDOS() const
G4double GetAnhDecConstant() const
G4double GetGamma() const
G4double GetMu() const
G4double GetSTDOS() const
G4double GetBeta() const
G4double GetLambda() const
void AddSecondary(G4Track *aSecondary)
void ProposeEnergy(G4double finalEnergy)
virtual void Initialize(const G4Track &)
virtual G4double GetMeanFreePath(const G4Track &, G4double, G4ForceCondition *)
virtual G4bool IsApplicable(const G4ParticleDefinition &)
G4PhononDownconversion(const G4String &processName="phononDownconversion")
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
static G4PhononLong * PhononDefinition()
Definition: G4PhononLong.cc:73
const G4ThreeVector & GetK(const G4Track *track) const
Definition: G4Step.hh:62
G4double GetVelocity() const
G4double GetKineticEnergy() const
void ProposeTrackStatus(G4TrackStatus status)
void SetNumberOfSecondaries(G4int totSecondaries)
virtual G4int ChoosePolarization(G4double Ldos, G4double STdos, G4double FTdos) const
virtual G4Track * CreateSecondary(G4int polarization, const G4ThreeVector &K, G4double energy) const
const G4LatticePhysical * theLattice
G4PhononTrackMap * trackKmap
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:327
G4int verboseLevel
Definition: G4VProcess.hh:356