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