Geant4  10.02.p01
XVCrystalPlanarAnalytical.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 //
26 
28 #include "G4PhysicsLinearVector.hh"
29 
31  fNumberOfPlanes = 4;
32 }
33 
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35 
37 }
38 
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 
42  fNumberOfPlanes = vNumberOfPlanes;
43 }
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 
48  return fNumberOfPlanes;
49 }
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52 
54 ComputeEC(G4ThreeVector vPositionVector,
55  XPhysicalLattice* vLattice){
56  G4VPhysicalVolume* vVolume = vLattice->GetVolume();
57 
58  G4double vInterplanarDistance =
60 
61  G4double vPosition = ComputePositionInUnitCell(vPositionVector,
62  vLattice).x();
63 
64  G4double vValue = 0.;
65  for(int i=-int(GetNumberOfPlanes()/2);i<=+int(GetNumberOfPlanes()/2);i++){
66  vValue += ComputeECForSinglePlane( ( vPosition + G4double(i) ) *
67  vInterplanarDistance ,
68  vLattice );
69  }
70 
71  return G4ThreeVector(vValue,0.,0.);
72 }
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
78  G4double vInterplanarPeriod = fPhysicalLattice->ComputeInterplanarPeriod();
79  G4double vX = vPosition.x();
80  if (vX < 0.0) {
81  vX += ((int( - vX / vInterplanarPeriod) + 1.0 )
82  * vInterplanarPeriod);
83  }
84  else if( vX > vInterplanarPeriod ){
85  vX -= ( int( vX / vInterplanarPeriod) * vInterplanarPeriod );
86  }
87  return G4ThreeVector(fVectorEC->Value(vX),0.,0.);
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
94  XPhysicalLattice* vLattice){
95  G4VPhysicalVolume* vVolume = vLattice->GetVolume();
96 
97  G4double vInterplanarPeriod =
99 
100  G4double vPositionX = vPosition.x();
101 
102  if((vPositionX >= 0.) &&
103  (vPositionX < vInterplanarPeriod)){
104  return G4ThreeVector(vPositionX/vInterplanarPeriod,0.,0.);
105  }
106  else if(vPositionX == vInterplanarPeriod){
107  return G4ThreeVector(0.,0.,0.);
108  }
109  else{
110  vPositionX -= std::fmod(vPosition.x(),vInterplanarPeriod)
111  * vInterplanarPeriod;
112  return G4ThreeVector(vPositionX/vInterplanarPeriod,0.,0.);
113  }
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117 
119  unsigned int vPrecision = 1024;
120  G4VPhysicalVolume* vVolume = vLattice->GetVolume();
122  / vPrecision;
123 
124  G4double vMaximum = -DBL_MAX;
125  G4double vValue;
126 
127  for(unsigned int i=0;i<vPrecision;i++){
128  if( (vValue = GetEC(G4ThreeVector(vStep * i,0.,0.),vLattice).x() )
129  > vMaximum) {vMaximum = vValue;}
130  }
131  return vMaximum;
132 }
133 
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135 
137  unsigned int vPrecision = 1024;
138  G4VPhysicalVolume* vVolume = vLattice->GetVolume();
140  / vPrecision;
141 
142  G4double vMinimum = +DBL_MAX;
143  G4double vValue;
144 
145  for(unsigned int i=0;i<vPrecision;i++){
146  if( (vValue = GetEC(G4ThreeVector(vStep * i,0.,0.),vLattice).x() )
147  < vMinimum) {vMinimum = vValue;}
148  }
149  return vMinimum;
150 }
151 
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
153 
155  XPhysicalLattice* vLattice,
156  G4double vUnit){
157  std::ofstream vFileOut;
158  vFileOut.open(filename);
159  vFileOut << "pos val" << std::endl;
160 
161  G4int imax = 8192;
162  G4double vXposition = 0.;
163  G4double vInterplanarPeriod = vLattice->ComputeInterplanarPeriod();
164 
165  for(G4int i = 0;i<imax;i++){
166  vXposition = double(i) / double(imax) * vInterplanarPeriod;
167  vFileOut << vXposition / CLHEP::angstrom << " ";
168  vFileOut << GetEC(G4ThreeVector(vXposition,0.,0.),vLattice).x() / vUnit;
169  vFileOut << std::endl;
170  }
171 
172  vFileOut.close();
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
176 
179  std::ifstream vFileIn;
180  vFileIn.open(filename);
181  std::string vTempString;
182  vFileIn >> vTempString >> vTempString;
183 
184  double vTempX;
185  double vTempVal;
186  std::vector<double> vTempXvector;
187  std::vector<double> vTempValvector;
188 
189  while(!vFileIn.eof()){
190  vFileIn >> vTempX;
191  vFileIn >> vTempVal;
192  vTempXvector.push_back(vTempX);
193  vTempValvector.push_back(vTempVal);
194  }
195 
196  vFileIn.close();
197 
198  G4int imax = vTempXvector.size();
199  G4double vInterplanarPeriod = vTempXvector.at(imax-1);
200 
202  vInterplanarPeriod*(imax-1)/imax,
203  imax);
204  for(G4int i = 0;i<imax;i++){
205  fVectorEC->PutValue(i,vTempValvector.at(i));
206  }
207 }
208 
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
210 
212  G4double vConversion){
213  std::ifstream vFileIn;
214  vFileIn.open(filename);
215 
216  G4int imax;
217  G4double xmax;
218 
219  vFileIn >> imax;
220  vFileIn >> xmax;
221 
222  xmax *= CLHEP::meter;
223  fMaximum = -DBL_MAX;
224  fMinimum = +DBL_MAX;
225 
226  fVectorEC = new G4PhysicsLinearVector(0,xmax,imax);
227 
228  for(G4int i=0;i<imax; i++){
229  double vTempX;
230  vFileIn >> vTempX;
231 
232  vTempX *= vConversion;
233  if(vTempX > fMaximum) {fMaximum = vTempX;}
234  if(vTempX < fMinimum) {fMinimum = vTempX;}
235  fVectorEC->PutValue(i,vTempX);
236  }
237 
238  G4cout << "XVCrystalPlanarAnalytical::ReadFromECHARM() - " <<
239  vConversion << " " << imax << " " << xmax << " " <<
240  fMinimum << " " << fMaximum << G4endl;
241 
242  vFileIn.close();
243 
244 }
245 
246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
247 
249  G4int imax = 4096;
250  G4double vXposition = 0.;
251  G4double vInterplanarPeriod = fPhysicalLattice->ComputeInterplanarPeriod();
252 
254  vInterplanarPeriod*(imax-1)/imax,
255  imax);
256  for(G4int i = 0;i<imax;i++){
257  vXposition = double(i) / double(imax) * vInterplanarPeriod;
258  fVectorEC->PutValue(i,ComputeEC(G4ThreeVector(vXposition,0.,0.),
259  fPhysicalLattice).x());
260  }
261 }
262 
263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
virtual G4double ComputeECForSinglePlane(G4double, XPhysicalLattice *)=0
CLHEP::Hep3Vector G4ThreeVector
virtual void ReadFromFile(const G4String &, XPhysicalLattice *, G4double=1)
XPhysicalLattice * GetXPhysicalLattice(G4VPhysicalVolume *)
int G4int
Definition: G4Types.hh:78
virtual G4double ComputeMaximum(XPhysicalLattice *)
G4GLOB_DLL std::ostream G4cout
static const double meter
Definition: G4SIunits.hh:81
XPhysicalLattice * fPhysicalLattice
G4ThreeVector ComputeECFromVector(G4ThreeVector)
void PutValue(size_t index, G4double theValue)
G4double Value(G4double theEnergy, size_t &lastidx) const
G4ThreeVector ComputePositionInUnitCell(G4ThreeVector, XPhysicalLattice *)
virtual void PrintOnFile(const G4String &, XPhysicalLattice *, G4double)
G4ThreeVector ComputeEC(G4ThreeVector, XPhysicalLattice *)
G4ThreeVector GetEC(G4ThreeVector, XPhysicalLattice *)
virtual G4double ComputeMinimum(XPhysicalLattice *)
const G4double x[NPOINTSGL]
#define G4endl
Definition: G4ios.hh:61
static const G4int imax
double G4double
Definition: G4Types.hh:76
virtual void ReadFromECHARM(const G4String &, G4double=1)
#define DBL_MAX
Definition: templates.hh:83
G4double ComputeInterplanarPeriod()
G4VPhysicalVolume * GetVolume()
static const double angstrom
Definition: G4SIunits.hh:101