| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file processes/phonon/src/G4PhononDowncon
27 /// \brief Implementation of the G4PhononDownc
28 //
29 //
30 // 20131111 Add verbose output for MFP calcul
31 // 20131115 Initialize data buffers in ctor
32
33 #include "G4PhononDownconversion.hh"
34 #include "G4LatticePhysical.hh"
35 #include "G4PhononLong.hh"
36 #include "G4PhononPolarization.hh"
37 #include "G4PhononTrackMap.hh"
38 #include "G4PhononTransFast.hh"
39 #include "G4PhononTransSlow.hh"
40 #include "G4PhysicalConstants.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
50 G4PhononDownconversion::G4PhononDownconversion
51 : G4VPhononProcess(aName), fBeta(0.), fGamma
52
53 G4PhononDownconversion::~G4PhononDownconversio
54
55 //....oooOO0OOooo........oooOO0OOooo........oo
56
57 G4double G4PhononDownconversion::GetMeanFreePa
58 G4double /*previousStepSize*/,
59 G4ForceCondition* condition) {
60 //Determines mean free path for longitudinal
61 G4double A = theLattice->GetAnhDecConstant()
62 G4double Eoverh = aTrack.GetKineticEnergy()/
63
64 //Calculate mean free path for anh. decay
65 G4double mfp = aTrack.GetVelocity()/(Eoverh*
66
67 if (verboseLevel > 1)
68 G4cout << "G4PhononDownconversion::GetMean
69
70 *condition = NotForced;
71 return mfp;
72 }
73
74 //....oooOO0OOooo........oooOO0OOooo........oo
75
76
77 G4VParticleChange* G4PhononDownconversion::Pos
78 const G4Step&) {
79 aParticleChange.Initialize(aTrack);
80
81 //Obtain dynamical constants from this volum
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 d
88 //74% chance that daughter phonons are both
89 //26% Transverse and Longitudinal
90 if (G4UniformRand()>0.740) MakeLTSecondaries
91 else MakeTTSecondaries(aTrack);
92
93 aParticleChange.ProposeEnergy(0.);
94 aParticleChange.ProposeTrackStatus(fStopAndK
95
96 return &aParticleChange;
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oo
100
101 G4bool G4PhononDownconversion::IsApplicable(co
102 //Only L-phonons decay
103 return (&aPD==G4PhononLong::PhononDefinition
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oo
107
108 //probability density of energy distribution o
109
110 G4double G4PhononDownconversion::GetLTDecayPro
111 //d=delta= ratio of group velocities vl/vt a
112 return (1/(x*x))*(1-x*x)*(1-x*x)*((1+x)*(1+x
113 }
114
115 //....oooOO0OOooo........oooOO0OOooo........oo
116
117 //probability density of energy distribution o
118 G4double G4PhononDownconversion::GetTTDecayPro
119 //dynamic constants from Tamura, PRL31, 1985
120 G4double A = 0.5*(1-d*d)*(fBeta+fLambda+(1+d
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+fLamb
124
125 return (A+B*d*x-B*x*x)*(A+B*d*x-B*x*x)+(C*x*
126 }
127
128 //....oooOO0OOooo........oooOO0OOooo........oo
129
130
131 G4double G4PhononDownconversion::MakeLDeviatio
132 //change in L'-phonon propagation direction
133
134 return std::acos((1+(x*x)-((d*d)*(1-x)*(1-x)
135 }
136
137 //....oooOO0OOooo........oooOO0OOooo........oo
138
139
140 G4double G4PhononDownconversion::MakeTDeviatio
141 //change in T-phonon propagation direction a
142
143 return std::acos((1-x*x+d*d*(1-x)*(1-x))/(2*
144 }
145
146 //....oooOO0OOooo........oooOO0OOooo........oo
147
148
149 G4double G4PhononDownconversion::MakeTTDeviati
150 //change in T-phonon propagation direction a
151
152 return std::acos((1-d*d*(1-x)*(1-x)+d*d*x*x)
153 }
154
155 //....oooOO0OOooo........oooOO0OOooo........oo
156
157
158 //Generate daughter phonons from L->T+T proces
159
160 void G4PhononDownconversion::MakeTTSecondaries
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 distr
167 //if a random point on the energy-probabilit
168 //smaller that the curve of the probability
169 //then accept that point.
170 //x=fraction of parent phonon energy in firs
171 G4double x = G4UniformRand()*(upperBound-low
172 G4double p = 1.5*G4UniformRand();
173 while(p >= GetTTDecayProb(d, x*d)) {
174 x = G4UniformRand()*(upperBound-lowerBound
175 p = 1.5*G4UniformRand();
176 }
177
178 //using energy fraction x to calculate daugh
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
185 G4ThreeVector ran = G4RandomDirection(); //
186
187 G4double ph=G4UniformRand()*twopi;
188 dir1 = dir1.rotate(dir1.orthogonal(),theta1)
189 dir2 = dir2.rotate(dir2.orthogonal(),-theta2
190
191 G4double E=aTrack.GetKineticEnergy();
192 G4double Esec1 = x*E, Esec2 = E-Esec1;
193
194 // Make FT or ST phonon (0. means no longitu
195 G4int polarization1 = ChoosePolarization(0.,
196 theLattice->GetFTDOS());
197
198 // Make FT or ST phonon (0. means no longitu
199 G4int polarization2 = ChoosePolarization(0.,
200 theLattice->GetFTDOS());
201
202 // Construct the secondaries and set their w
203 G4Track* sec1 = CreateSecondary(polarization
204 G4Track* sec2 = CreateSecondary(polarization
205
206 aParticleChange.SetNumberOfSecondaries(2);
207 aParticleChange.AddSecondary(sec1);
208 aParticleChange.AddSecondary(sec2);
209 }
210
211 //....oooOO0OOooo........oooOO0OOooo........oo
212
213
214 //Generate daughter phonons from L->L'+T proce
215
216 void G4PhononDownconversion::MakeLTSecondaries
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 distr
223 //if a random point on the energy-probabilit
224 //smaller that the curve of the probability
225 //then accept that point.
226 //x=fraction of parent phonon energy in L ph
227 G4double x = G4UniformRand()*(upperBound-low
228 G4double p = 4.0*G4UniformRand();
229 while(p >= GetLTDecayProb(d, x)) {
230 x = G4UniformRand()*(upperBound-lowerBound
231 p = 4.0*G4UniformRand(); //4.0 is
232 }
233
234 //using energy fraction x to calculate daugh
235 G4double thetaL=MakeLDeviation(d, x);
236 G4double thetaT=MakeTDeviation(d, x); // F
237 G4ThreeVector dir1=trackKmap->GetK(aTrack);
238 G4ThreeVector dir2=dir1;
239
240 G4double ph=G4UniformRand()*twopi;
241 dir1 = dir1.rotate(dir1.orthogonal(),thetaL)
242 dir2 = dir2.rotate(dir2.orthogonal(),-thetaT
243
244 G4double E=aTrack.GetKineticEnergy();
245 G4double Esec1 = x*E, Esec2 = E-Esec1;
246
247 // First secondary is longitudnal
248 G4int polarization1 = G4PhononPolarization::
249
250 // Make FT or ST phonon (0. means no longitu
251 G4int polarization2 = ChoosePolarization(0.,
252 theLattice->GetFTDOS());
253
254 // Construct the secondaries and set their w
255 G4Track* sec1 = CreateSecondary(polarization
256 G4Track* sec2 = CreateSecondary(polarization
257
258 aParticleChange.SetNumberOfSecondaries(2);
259 aParticleChange.AddSecondary(sec1);
260 aParticleChange.AddSecondary(sec2);
261 }
262
263 //....oooOO0OOooo........oooOO0OOooo........oo
264
265