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