| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 << 25 //
26 /// \file materials/src/G4LatticeLogical.cc 26 /// \file materials/src/G4LatticeLogical.cc
27 /// \brief Implementation of the G4LatticeLogi 27 /// \brief Implementation of the G4LatticeLogical class
28 << 28 //
>> 29 //
29 #include "G4LatticeLogical.hh" 30 #include "G4LatticeLogical.hh"
30 <<
31 #include "G4PhysicalConstants.hh" <<
32 #include "G4SystemOfUnits.hh" 31 #include "G4SystemOfUnits.hh"
33 << 32 #include "G4PhysicalConstants.hh"
34 #include <cmath> 33 #include <cmath>
35 #include <fstream> 34 #include <fstream>
36 35
37 36
>> 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 38
38 G4LatticeLogical::G4LatticeLogical() 39 G4LatticeLogical::G4LatticeLogical()
39 { << 40 : verboseLevel(0), fVresTheta(0), fVresPhi(0), fDresTheta(0), fDresPhi(0),
40 for (G4int i = 0; i < 3; i++) { << 41 fA(0), fB(0), fLDOS(0), fSTDOS(0), fFTDOS(0),
41 for (G4int j = 0; j < MAXRES; j++) { << 42 fBeta(0), fGamma(0), fLambda(0), fMu(0) {
42 for (G4int k = 0; k < MAXRES; k++) { << 43 for (G4int i=0; i<3; i++) {
43 fMap[i][j][k] = 0.; << 44 for (G4int j=0; j<MAXRES; j++) {
44 fN_map[i][j][k].set(0., 0., 0.); << 45 for (G4int k=0; k<MAXRES; k++) {
>> 46 fMap[i][j][k] = 0.;
>> 47 fN_map[i][j][k].set(0.,0.,0.);
45 } 48 }
46 } 49 }
47 } 50 }
48 } 51 }
49 52
>> 53 G4LatticeLogical::~G4LatticeLogical() {;}
>> 54
>> 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 56
>> 57
50 /////////////////////////////////////////// 58 ///////////////////////////////////////////
51 // Load map of group velocity scalars (m/s) << 59 //Load map of group velocity scalars (m/s)
52 //////////////////////////////////////////// 60 ////////////////////////////////////////////
53 G4bool G4LatticeLogical::LoadMap(G4int tRes, G << 61 G4bool G4LatticeLogical::LoadMap(G4int tRes, G4int pRes,
54 { << 62 G4int polarizationState, G4String map) {
55 if (tRes > MAXRES || pRes > MAXRES) { << 63 if (tRes>MAXRES || pRes>MAXRES) {
56 G4cerr << "G4LatticeLogical::LoadMap excee << 64 G4cerr << "G4LatticeLogical::LoadMap exceeds maximum resolution of "
57 << MAXRES << ". terminating." << G4 << 65 << MAXRES << " by " << MAXRES << ". terminating." << G4endl;
58 return false; // terminate if resolution << 66 return false; //terminate if resolution out of bounds.
59 } 67 }
60 68
61 std::ifstream fMapFile(map.data()); 69 std::ifstream fMapFile(map.data());
62 if (! fMapFile.is_open()) { << 70 if(!fMapFile.is_open())
>> 71 {
63 return false; 72 return false;
64 } 73 }
65 74
66 G4double vgrp = 0.; 75 G4double vgrp = 0.;
67 for (G4int theta = 0; theta < tRes; theta++) << 76 for (G4int theta = 0; theta<tRes; theta++) {
68 for (G4int phi = 0; phi < pRes; phi++) { << 77 for (G4int phi = 0; phi<pRes; phi++) {
69 fMapFile >> vgrp; 78 fMapFile >> vgrp;
70 fMap[polarizationState][theta][phi] = vg << 79 fMap[polarizationState][theta][phi] = vgrp * (m/s);
71 } 80 }
72 } 81 }
73 82
74 if (verboseLevel != 0) { << 83 if(verboseLevel != 0)
>> 84 {
75 G4cout << "\nG4LatticeLogical::LoadMap(" < 85 G4cout << "\nG4LatticeLogical::LoadMap(" << map << ") successful"
76 << " (Vg scalars " << tRes << " x " << 86 << " (Vg scalars " << tRes << " x " << pRes << " for polarization "
77 << ")." << G4endl; << 87 << polarizationState << ")." << G4endl;
78 } 88 }
79 89
80 fVresTheta = tRes; // store map dimensions << 90 fVresTheta=tRes; //store map dimensions
81 fVresPhi = pRes; << 91 fVresPhi=pRes;
82 return true; 92 return true;
83 } 93 }
84 94
>> 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 96
85 97
86 //////////////////////////////////// 98 ////////////////////////////////////
87 // Load map of group velocity unit vectors << 99 //Load map of group velocity unit vectors
88 /////////////////////////////////// 100 ///////////////////////////////////
89 G4bool G4LatticeLogical::Load_NMap(G4int tRes, << 101 G4bool G4LatticeLogical::Load_NMap(G4int tRes, G4int pRes,
90 { << 102 G4int polarizationState, G4String map) {
91 if (tRes > MAXRES || pRes > MAXRES) { << 103 if (tRes>MAXRES || pRes>MAXRES) {
92 G4cerr << "G4LatticeLogical::LoadMap excee << 104 G4cerr << "G4LatticeLogical::LoadMap exceeds maximum resolution of "
93 << MAXRES << ". terminating." << G4 << 105 << MAXRES << " by " << MAXRES << ". terminating." << G4endl;
94 return false; // terminate if resolution << 106 return false; //terminate if resolution out of bounds.
95 } 107 }
96 108
97 std::ifstream fMapFile(map.data()); 109 std::ifstream fMapFile(map.data());
98 if (! fMapFile.is_open()) { << 110 if(!fMapFile.is_open())
>> 111 {
99 return false; 112 return false;
100 } 113 }
101 114
102 G4double x, y, z; // Buffers to read coordi << 115 G4double x,y,z; // Buffers to read coordinates from file
103 G4ThreeVector dir; 116 G4ThreeVector dir;
104 for (G4int theta = 0; theta < tRes; theta++) << 117 for (G4int theta = 0; theta<tRes; theta++) {
105 for (G4int phi = 0; phi < pRes; phi++) { << 118 for (G4int phi = 0; phi<pRes; phi++) {
106 fMapFile >> x >> y >> z; 119 fMapFile >> x >> y >> z;
107 dir.set(x, y, z); << 120 dir.set(x,y,z);
108 fN_map[polarizationState][theta][phi] = << 121 fN_map[polarizationState][theta][phi] = dir.unit(); // Enforce unity
109 } 122 }
110 } 123 }
111 124
112 if (verboseLevel != 0) { << 125 if(verboseLevel != 0)
>> 126 {
113 G4cout << "\nG4LatticeLogical::Load_NMap(" 127 G4cout << "\nG4LatticeLogical::Load_NMap(" << map << ") successful"
114 << " (Vdir " << tRes << " x " << pR << 128 << " (Vdir " << tRes << " x " << pRes << " for polarization "
115 << ")." << G4endl; << 129 << polarizationState << ")." << G4endl;
116 } 130 }
117 131
118 fDresTheta = tRes; // store map dimensions << 132 fDresTheta=tRes; //store map dimensions
119 fDresPhi = pRes; << 133 fDresPhi=pRes;
120 return true; 134 return true;
121 } 135 }
122 136
>> 137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 138
124 // Given the phonon wave vector k, phonon phys << 139 //Given the phonon wave vector k, phonon physical volume Vol
125 // and polarizationState(0=LON, 1=FT, 2=ST), << 140 //and polarizationState(0=LON, 1=FT, 2=ST),
126 // returns phonon velocity in m/s << 141 //returns phonon velocity in m/s
127 142
128 G4double G4LatticeLogical::MapKtoV(G4int polar << 143 G4double G4LatticeLogical::MapKtoV(G4int polarizationState,
129 { << 144 const G4ThreeVector& k) const {
130 G4double theta, phi, tRes, pRes; 145 G4double theta, phi, tRes, pRes;
131 146
132 tRes = pi / fVresTheta; << 147 tRes=pi/fVresTheta;
133 pRes = twopi / fVresPhi; << 148 pRes=twopi/fVresPhi;
134 << 149
135 theta = k.getTheta(); << 150 theta=k.getTheta();
136 phi = k.getPhi(); << 151 phi=k.getPhi();
137 152
138 if (phi < 0) { << 153 if(phi < 0)
>> 154 {
139 phi = phi + twopi; 155 phi = phi + twopi;
140 } 156 }
141 if (theta > pi) { << 157 if(theta > pi)
>> 158 {
142 theta = theta - pi; 159 theta = theta - pi;
143 } 160 }
144 161
145 G4double Vg = fMap[polarizationState][int(th << 162 G4double Vg = fMap[polarizationState][int(theta/tRes)][int(phi/pRes)];
146 163
147 if (Vg == 0) { << 164 if(Vg == 0){
148 G4cout << "\nFound v=0 for polarization " << 165 G4cout<<"\nFound v=0 for polarization "<<polarizationState
149 << phi << " translating to map coor << 166 <<" theta "<<theta<<" phi "<<phi<< " translating to map coords "
150 << "theta " << int(theta / tRes) << << 167 <<"theta "<< int(theta/tRes) << " phi " << int(phi/pRes)<<G4endl;
151 } 168 }
152 169
153 if (verboseLevel > 1) { << 170 if (verboseLevel>1) {
154 G4cout << "G4LatticeLogical::MapKtoV theta << 171 G4cout << "G4LatticeLogical::MapKtoV theta,phi=" << theta << " " << phi
155 << int(theta / tRes) << " " << int( << 172 << " : ith,iph " << int(theta/tRes) << " " << int(phi/pRes)
>> 173 << " : V " << Vg << G4endl;
156 } 174 }
157 175
158 return Vg; 176 return Vg;
159 } 177 }
160 178
>> 179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161 180
162 // Given the phonon wave vector k, phonon phys << 181 //Given the phonon wave vector k, phonon physical volume Vol
163 // and polarizationState(0=LON, 1=FT, 2=ST), << 182 //and polarizationState(0=LON, 1=FT, 2=ST),
164 // returns phonon propagation direction as dim << 183 //returns phonon propagation direction as dimensionless unit vector
165 184
166 G4ThreeVector G4LatticeLogical::MapKtoVDir(G4i << 185 G4ThreeVector G4LatticeLogical::MapKtoVDir(G4int polarizationState,
167 { << 186 const G4ThreeVector& k) const {
168 G4double theta, phi, tRes, pRes; 187 G4double theta, phi, tRes, pRes;
169 188
170 tRes = pi / (fDresTheta - 1); // The summan << 189 tRes=pi/(fDresTheta-1);//The summant "-1" is required:index=[0:array length-1]
171 pRes = 2 * pi / (fDresPhi - 1); << 190 pRes=2*pi/(fDresPhi-1);
172 191
173 theta = k.getTheta(); << 192 theta=k.getTheta();
174 phi = k.getPhi(); << 193 phi=k.getPhi();
175 194
176 if (theta > pi) { << 195 if(theta > pi)
>> 196 {
177 theta = theta - pi; 197 theta = theta - pi;
178 } 198 }
179 // phi=[0 to 2 pi] in accordance with DMC // << 199 //phi=[0 to 2 pi] in accordance with DMC //if(phi>pi/2) phi=phi-pi/2;
180 if (phi < 0) { << 200 if(phi < 0)
>> 201 {
181 phi = phi + 2 * pi; 202 phi = phi + 2 * pi;
182 } 203 }
183 204
184 auto iTheta = int(theta / tRes + 0.5); << 205 G4int iTheta = int(theta/tRes+0.5);
185 auto iPhi = int(phi / pRes + 0.5); << 206 G4int iPhi = int(phi/pRes+0.5);
186 207
187 if (verboseLevel > 1) { << 208 if (verboseLevel>1) {
188 G4cout << "G4LatticeLogical::MapKtoVDir th << 209 G4cout << "G4LatticeLogical::MapKtoVDir theta,phi=" << theta << " " << phi
189 << iTheta << " " << iPhi << " : dir << 210 << " : ith,iph " << iTheta << " " << iPhi
190 << G4endl; << 211 << " : dir " << fN_map[polarizationState][iTheta][iPhi] << G4endl;
191 } 212 }
192 213
193 return fN_map[polarizationState][iTheta][iPh 214 return fN_map[polarizationState][iTheta][iPhi];
194 } 215 }
195 216
>> 217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196 218
197 // Dump structure in format compatible with re 219 // Dump structure in format compatible with reading back
198 220
199 void G4LatticeLogical::Dump(std::ostream& os) << 221 void G4LatticeLogical::Dump(std::ostream& os) const {
200 { << 222 os << "dyn " << fBeta << " " << fGamma << " " << fLambda << " " << fMu
201 os << "dyn " << fBeta << " " << fGamma << " << 223 << "\nscat " << fB << " decay " << fA
202 << " decay " << fA << "\nLDOS " << fLDOS << 224 << "\nLDOS " << fLDOS << " STDOS " << fSTDOS << " FTDOS " << fFTDOS
203 << std::endl; 225 << std::endl;
204 226
205 Dump_NMap(os, 0, "LVec.ssv"); 227 Dump_NMap(os, 0, "LVec.ssv");
206 Dump_NMap(os, 1, "FTVec.ssv"); 228 Dump_NMap(os, 1, "FTVec.ssv");
207 Dump_NMap(os, 2, "STVec.ssv"); 229 Dump_NMap(os, 2, "STVec.ssv");
208 230
209 DumpMap(os, 0, "L.ssv"); 231 DumpMap(os, 0, "L.ssv");
210 DumpMap(os, 1, "FT.ssv"); 232 DumpMap(os, 1, "FT.ssv");
211 DumpMap(os, 2, "ST.ssv"); 233 DumpMap(os, 2, "ST.ssv");
212 } 234 }
213 235
214 void G4LatticeLogical::DumpMap(std::ostream& o << 236 void G4LatticeLogical::DumpMap(std::ostream& os, G4int pol,
215 { << 237 const G4String& name) const {
216 os << "VG " << name << " " << 238 os << "VG " << name << " " << (pol==0?"L":pol==1?"FT":pol==2?"ST":"??")
217 << (pol == 0 ? "L" <<
218 : pol == 1 ? "FT" <<
219 : pol == 2 ? "ST" <<
220 : "??") <<
221 << " " << fVresTheta << " " << fVresPhi < 239 << " " << fVresTheta << " " << fVresPhi << std::endl;
222 240
223 for (G4int iTheta = 0; iTheta < fVresTheta; << 241 for (G4int iTheta=0; iTheta<fVresTheta; iTheta++) {
224 for (G4int iPhi = 0; iPhi < fVresPhi; iPhi << 242 for (G4int iPhi=0; iPhi<fVresPhi; iPhi++) {
225 os << fMap[pol][iTheta][iPhi] << std::en 243 os << fMap[pol][iTheta][iPhi] << std::endl;
226 } 244 }
227 } 245 }
228 } 246 }
229 247
230 void G4LatticeLogical::Dump_NMap(std::ostream& << 248 void G4LatticeLogical::Dump_NMap(std::ostream& os, G4int pol,
231 { << 249 const G4String& name) const {
232 os << "VDir " << name << " " << 250 os << "VDir " << name << " " << (pol==0?"L":pol==1?"FT":pol==2?"ST":"??")
233 << (pol == 0 ? "L" <<
234 : pol == 1 ? "FT" <<
235 : pol == 2 ? "ST" <<
236 : "??") <<
237 << " " << fDresTheta << " " << fDresPhi < 251 << " " << fDresTheta << " " << fDresPhi << std::endl;
238 252
239 for (G4int iTheta = 0; iTheta < fDresTheta; << 253 for (G4int iTheta=0; iTheta<fDresTheta; iTheta++) {
240 for (G4int iPhi = 0; iPhi < fDresPhi; iPhi << 254 for (G4int iPhi=0; iPhi<fDresPhi; iPhi++) {
241 os << fN_map[pol][iTheta][iPhi].x() << " << 255 os << fN_map[pol][iTheta][iPhi].x()
242 << fN_map[pol][iTheta][iPhi].z() << s << 256 << " " << fN_map[pol][iTheta][iPhi].y()
>> 257 << " " << fN_map[pol][iTheta][iPhi].z()
>> 258 << std::endl;
243 } 259 }
244 } 260 }
245 } 261 }
>> 262
246 263