Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/solidstate/phonon/src/G4PhononDownconversion.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 
 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(const G4String& aName)
 51   : G4VPhononProcess(aName), fBeta(0.), fGamma(0.), fLambda(0.), fMu(0.) {;}
 52 
 53 G4PhononDownconversion::~G4PhononDownconversion() {;}
 54 
 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 56 
 57 G4double G4PhononDownconversion::GetMeanFreePath(const G4Track& aTrack,
 58              G4double /*previousStepSize*/,
 59              G4ForceCondition* condition) {
 60   //Determines mean free path for longitudinal phonons to split
 61   G4double A = theLattice->GetAnhDecConstant();
 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   
 70   *condition = NotForced;
 71   return mfp;
 72 }
 73 
 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75 
 76 
 77 G4VParticleChange* G4PhononDownconversion::PostStepDoIt( const G4Track& aTrack,
 78                const G4Step&) {
 79   aParticleChange.Initialize(aTrack);
 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 
 93   aParticleChange.ProposeEnergy(0.);
 94   aParticleChange.ProposeTrackStatus(fStopAndKill);    
 95        
 96   return &aParticleChange;
 97 }
 98 
 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
101 G4bool G4PhononDownconversion::IsApplicable(const G4ParticleDefinition& aPD) {
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 
110 G4double 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
118 G4double 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 
131 G4double 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 
140 G4double 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 
149 G4double 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    
160 void 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(),
196              theLattice->GetFTDOS());
197 
198   // Make FT or ST phonon (0. means no longitudinal)
199   G4int polarization2 = ChoosePolarization(0., theLattice->GetSTDOS(),
200              theLattice->GetFTDOS());
201 
202   // Construct the secondaries and set their wavevectors
203   G4Track* sec1 = CreateSecondary(polarization1, dir1, Esec1);
204   G4Track* sec2 = CreateSecondary(polarization2, dir2, Esec2);
205 
206   aParticleChange.SetNumberOfSecondaries(2);
207   aParticleChange.AddSecondary(sec1);
208   aParticleChange.AddSecondary(sec2);
209 }
210 
211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
212 
213 
214 //Generate daughter phonons from L->L'+T process
215    
216 void 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(),
252              theLattice->GetFTDOS());
253 
254   // Construct the secondaries and set their wavevectors
255   G4Track* sec1 = CreateSecondary(polarization1, dir1, Esec1);
256   G4Track* sec2 = CreateSecondary(polarization2, dir2, Esec2);
257 
258   aParticleChange.SetNumberOfSecondaries(2);
259   aParticleChange.AddSecondary(sec1);
260   aParticleChange.AddSecondary(sec2); 
261 }
262 
263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
264 
265