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