ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PhononDownconversion.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PhononDownconversion.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 //
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"
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 
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 
70  *condition = NotForced;
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
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 
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
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 
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 
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 
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 
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 
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 
209 }
210 
211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
212 
213 
214 //Generate daughter phonons from L->L'+T process
215 
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 
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 
261 }
262 
263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
264