| 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()) return false;
63 return false; <<
64 } <<
65 71
66 G4double vgrp = 0.; 72 G4double vgrp = 0.;
67 for (G4int theta = 0; theta < tRes; theta++) << 73 for (G4int theta = 0; theta<tRes; theta++) {
68 for (G4int phi = 0; phi < pRes; phi++) { << 74 for (G4int phi = 0; phi<pRes; phi++) {
69 fMapFile >> vgrp; 75 fMapFile >> vgrp;
70 fMap[polarizationState][theta][phi] = vg << 76 fMap[polarizationState][theta][phi] = vgrp * (m/s);
71 } 77 }
72 } 78 }
73 79
74 if (verboseLevel != 0) { << 80 if (verboseLevel) {
75 G4cout << "\nG4LatticeLogical::LoadMap(" < 81 G4cout << "\nG4LatticeLogical::LoadMap(" << map << ") successful"
76 << " (Vg scalars " << tRes << " x " << 82 << " (Vg scalars " << tRes << " x " << pRes << " for polarization "
77 << ")." << G4endl; << 83 << polarizationState << ")." << G4endl;
78 } 84 }
79 85
80 fVresTheta = tRes; // store map dimensions << 86 fVresTheta=tRes; //store map dimensions
81 fVresPhi = pRes; << 87 fVresPhi=pRes;
82 return true; 88 return true;
83 } 89 }
84 90
>> 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 92
85 93
86 //////////////////////////////////// 94 ////////////////////////////////////
87 // Load map of group velocity unit vectors << 95 //Load map of group velocity unit vectors
88 /////////////////////////////////// 96 ///////////////////////////////////
89 G4bool G4LatticeLogical::Load_NMap(G4int tRes, << 97 G4bool G4LatticeLogical::Load_NMap(G4int tRes, G4int pRes,
90 { << 98 G4int polarizationState, G4String map) {
91 if (tRes > MAXRES || pRes > MAXRES) { << 99 if (tRes>MAXRES || pRes>MAXRES) {
92 G4cerr << "G4LatticeLogical::LoadMap excee << 100 G4cerr << "G4LatticeLogical::LoadMap exceeds maximum resolution of "
93 << MAXRES << ". terminating." << G4 << 101 << MAXRES << " by " << MAXRES << ". terminating." << G4endl;
94 return false; // terminate if resolution << 102 return false; //terminate if resolution out of bounds.
95 } 103 }
96 104
97 std::ifstream fMapFile(map.data()); 105 std::ifstream fMapFile(map.data());
98 if (! fMapFile.is_open()) { << 106 if(!fMapFile.is_open()) return false;
99 return false; <<
100 } <<
101 107
102 G4double x, y, z; // Buffers to read coordi << 108 G4double x,y,z; // Buffers to read coordinates from file
103 G4ThreeVector dir; 109 G4ThreeVector dir;
104 for (G4int theta = 0; theta < tRes; theta++) << 110 for (G4int theta = 0; theta<tRes; theta++) {
105 for (G4int phi = 0; phi < pRes; phi++) { << 111 for (G4int phi = 0; phi<pRes; phi++) {
106 fMapFile >> x >> y >> z; 112 fMapFile >> x >> y >> z;
107 dir.set(x, y, z); << 113 dir.set(x,y,z);
108 fN_map[polarizationState][theta][phi] = << 114 fN_map[polarizationState][theta][phi] = dir.unit(); // Enforce unity
109 } 115 }
110 } 116 }
111 117
112 if (verboseLevel != 0) { << 118 if (verboseLevel) {
113 G4cout << "\nG4LatticeLogical::Load_NMap(" 119 G4cout << "\nG4LatticeLogical::Load_NMap(" << map << ") successful"
114 << " (Vdir " << tRes << " x " << pR << 120 << " (Vdir " << tRes << " x " << pRes << " for polarization "
115 << ")." << G4endl; << 121 << polarizationState << ")." << G4endl;
116 } 122 }
117 123
118 fDresTheta = tRes; // store map dimensions << 124 fDresTheta=tRes; //store map dimensions
119 fDresPhi = pRes; << 125 fDresPhi=pRes;
120 return true; 126 return true;
121 } 127 }
122 128
>> 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 130
124 // Given the phonon wave vector k, phonon phys << 131 //Given the phonon wave vector k, phonon physical volume Vol
125 // and polarizationState(0=LON, 1=FT, 2=ST), << 132 //and polarizationState(0=LON, 1=FT, 2=ST),
126 // returns phonon velocity in m/s << 133 //returns phonon velocity in m/s
127 134
128 G4double G4LatticeLogical::MapKtoV(G4int polar << 135 G4double G4LatticeLogical::MapKtoV(G4int polarizationState,
129 { << 136 const G4ThreeVector& k) const {
130 G4double theta, phi, tRes, pRes; 137 G4double theta, phi, tRes, pRes;
131 138
132 tRes = pi / fVresTheta; << 139 tRes=pi/fVresTheta;
133 pRes = twopi / fVresPhi; << 140 pRes=twopi/fVresPhi;
134 << 141
135 theta = k.getTheta(); << 142 theta=k.getTheta();
136 phi = k.getPhi(); << 143 phi=k.getPhi();
137 << 144
138 if (phi < 0) { << 145 if(phi<0) phi = phi + twopi;
139 phi = phi + twopi; << 146 if(theta>pi) theta=theta-pi;
140 } << 147
141 if (theta > pi) { << 148 G4double Vg = fMap[polarizationState][int(theta/tRes)][int(phi/pRes)];
142 theta = theta - pi; << 149
143 } << 150 if(Vg == 0){
144 << 151 G4cout<<"\nFound v=0 for polarization "<<polarizationState
145 G4double Vg = fMap[polarizationState][int(th << 152 <<" theta "<<theta<<" phi "<<phi<< " translating to map coords "
146 << 153 <<"theta "<< int(theta/tRes) << " phi " << int(phi/pRes)<<G4endl;
147 if (Vg == 0) { << 154 }
148 G4cout << "\nFound v=0 for polarization " << 155
149 << phi << " translating to map coor << 156 if (verboseLevel>1) {
150 << "theta " << int(theta / tRes) << << 157 G4cout << "G4LatticeLogical::MapKtoV theta,phi=" << theta << " " << phi
151 } << 158 << " : ith,iph " << int(theta/tRes) << " " << int(phi/pRes)
152 << 159 << " : V " << Vg << G4endl;
153 if (verboseLevel > 1) { <<
154 G4cout << "G4LatticeLogical::MapKtoV theta <<
155 << int(theta / tRes) << " " << int( <<
156 } 160 }
157 161
158 return Vg; 162 return Vg;
159 } 163 }
160 164
>> 165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161 166
162 // Given the phonon wave vector k, phonon phys << 167 //Given the phonon wave vector k, phonon physical volume Vol
163 // and polarizationState(0=LON, 1=FT, 2=ST), << 168 //and polarizationState(0=LON, 1=FT, 2=ST),
164 // returns phonon propagation direction as dim << 169 //returns phonon propagation direction as dimensionless unit vector
165 170
166 G4ThreeVector G4LatticeLogical::MapKtoVDir(G4i << 171 G4ThreeVector G4LatticeLogical::MapKtoVDir(G4int polarizationState,
167 { << 172 const G4ThreeVector& k) const {
168 G4double theta, phi, tRes, pRes; 173 G4double theta, phi, tRes, pRes;
169 174
170 tRes = pi / (fDresTheta - 1); // The summan << 175 tRes=pi/(fDresTheta-1);//The summant "-1" is required:index=[0:array length-1]
171 pRes = 2 * pi / (fDresPhi - 1); << 176 pRes=2*pi/(fDresPhi-1);
172 <<
173 theta = k.getTheta(); <<
174 phi = k.getPhi(); <<
175 <<
176 if (theta > pi) { <<
177 theta = theta - pi; <<
178 } <<
179 // phi=[0 to 2 pi] in accordance with DMC // <<
180 if (phi < 0) { <<
181 phi = phi + 2 * pi; <<
182 } <<
183 177
184 auto iTheta = int(theta / tRes + 0.5); << 178 theta=k.getTheta();
185 auto iPhi = int(phi / pRes + 0.5); << 179 phi=k.getPhi();
186 180
187 if (verboseLevel > 1) { << 181 if(theta>pi) theta=theta-pi;
188 G4cout << "G4LatticeLogical::MapKtoVDir th << 182 //phi=[0 to 2 pi] in accordance with DMC //if(phi>pi/2) phi=phi-pi/2;
189 << iTheta << " " << iPhi << " : dir << 183 if(phi<0) phi = phi + 2*pi;
190 << G4endl; << 184
>> 185 G4int iTheta = int(theta/tRes+0.5);
>> 186 G4int iPhi = int(phi/pRes+0.5);
>> 187
>> 188 if (verboseLevel>1) {
>> 189 G4cout << "G4LatticeLogical::MapKtoVDir theta,phi=" << theta << " " << phi
>> 190 << " : ith,iph " << iTheta << " " << iPhi
>> 191 << " : dir " << fN_map[polarizationState][iTheta][iPhi] << G4endl;
191 } 192 }
192 193
193 return fN_map[polarizationState][iTheta][iPh 194 return fN_map[polarizationState][iTheta][iPhi];
194 } 195 }
195 196
>> 197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196 198
197 // Dump structure in format compatible with re 199 // Dump structure in format compatible with reading back
198 200
199 void G4LatticeLogical::Dump(std::ostream& os) << 201 void G4LatticeLogical::Dump(std::ostream& os) const {
200 { << 202 os << "dyn " << fBeta << " " << fGamma << " " << fLambda << " " << fMu
201 os << "dyn " << fBeta << " " << fGamma << " << 203 << "\nscat " << fB << " decay " << fA
202 << " decay " << fA << "\nLDOS " << fLDOS << 204 << "\nLDOS " << fLDOS << " STDOS " << fSTDOS << " FTDOS " << fFTDOS
203 << std::endl; 205 << std::endl;
204 206
205 Dump_NMap(os, 0, "LVec.ssv"); 207 Dump_NMap(os, 0, "LVec.ssv");
206 Dump_NMap(os, 1, "FTVec.ssv"); 208 Dump_NMap(os, 1, "FTVec.ssv");
207 Dump_NMap(os, 2, "STVec.ssv"); 209 Dump_NMap(os, 2, "STVec.ssv");
208 210
209 DumpMap(os, 0, "L.ssv"); 211 DumpMap(os, 0, "L.ssv");
210 DumpMap(os, 1, "FT.ssv"); 212 DumpMap(os, 1, "FT.ssv");
211 DumpMap(os, 2, "ST.ssv"); 213 DumpMap(os, 2, "ST.ssv");
212 } 214 }
213 215
214 void G4LatticeLogical::DumpMap(std::ostream& o << 216 void G4LatticeLogical::DumpMap(std::ostream& os, G4int pol,
215 { << 217 const G4String& name) const {
216 os << "VG " << name << " " << 218 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 < 219 << " " << fVresTheta << " " << fVresPhi << std::endl;
222 220
223 for (G4int iTheta = 0; iTheta < fVresTheta; << 221 for (G4int iTheta=0; iTheta<fVresTheta; iTheta++) {
224 for (G4int iPhi = 0; iPhi < fVresPhi; iPhi << 222 for (G4int iPhi=0; iPhi<fVresPhi; iPhi++) {
225 os << fMap[pol][iTheta][iPhi] << std::en 223 os << fMap[pol][iTheta][iPhi] << std::endl;
226 } 224 }
227 } 225 }
228 } 226 }
229 227
230 void G4LatticeLogical::Dump_NMap(std::ostream& << 228 void G4LatticeLogical::Dump_NMap(std::ostream& os, G4int pol,
231 { << 229 const G4String& name) const {
232 os << "VDir " << name << " " << 230 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 < 231 << " " << fDresTheta << " " << fDresPhi << std::endl;
238 232
239 for (G4int iTheta = 0; iTheta < fDresTheta; << 233 for (G4int iTheta=0; iTheta<fDresTheta; iTheta++) {
240 for (G4int iPhi = 0; iPhi < fDresPhi; iPhi << 234 for (G4int iPhi=0; iPhi<fDresPhi; iPhi++) {
241 os << fN_map[pol][iTheta][iPhi].x() << " << 235 os << fN_map[pol][iTheta][iPhi].x()
242 << fN_map[pol][iTheta][iPhi].z() << s << 236 << " " << fN_map[pol][iTheta][iPhi].y()
>> 237 << " " << fN_map[pol][iTheta][iPhi].z()
>> 238 << std::endl;
243 } 239 }
244 } 240 }
245 } 241 }
>> 242
246 243