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