| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //------------------ G4XrayReflection physics
27 //
28 // History:
29 // 14-09-23 H. Burkhardt, initial implementati
30 //
31 // -------------------------------------------
32
33 #include "G4XrayReflection.hh"
34
35 #include "G4EmProcessSubType.hh"
36 #include "G4Exp.hh"
37 #include "G4Gamma.hh"
38 #include "G4Step.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4ThreeVector.hh"
41 #include "G4Track.hh"
42 #include "G4TransportationManager.hh"
43 #include "G4VDiscreteProcess.hh"
44
45 G4double G4XrayReflection::fSurfaceRoughness =
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 G4XrayReflection::G4XrayReflection(const G4Str
50 G4ProcessTy
51 : G4VDiscreteProcess(processName, type)
52 {
53 SetProcessSubType(fGammaReflection);
54 SaveHenkeDataAsMaterialProperty();
55 }
56
57 //....oooOO0OOooo........oooOO0OOooo........oo
58
59 G4bool G4XrayReflection::IsApplicable(const G4
60 {
61 return (&particle == G4Gamma::Gamma()); //
62 }
63
64 //....oooOO0OOooo........oooOO0OOooo........oo
65
66 G4double G4XrayReflection::Reflectivity(const
67 const
68 {
69 G4double theReflectivity = 0;
70 const G4MaterialPropertiesTable* theMatProp
71 if (SinIncidentAngle < 0.9 && theMatProp !=
72 { // avoid perpendicular refl. at straight
73 // data available
74 G4MaterialPropertyVector* RealIndex = theM
75 G4MaterialPropertyVector* ImagIndex = theM
76 if (nullptr == RealIndex || nullptr == Ima
77 const G4double delta = RealIndex->Value(Ga
78 const G4double beta = ImagIndex->Value(Gam
79 const G4double sin2 = std::pow(SinIncident
80 const G4double rho2 =
81 0.5 * (sin2 - 2 * delta + std::sqrt(std:
82 const G4double rho = std::sqrt(rho2);
83 const G4double Refl_sigma = (rho2 * std::p
84 / (rho2 * std:
85 const G4double coscot = std::sqrt(1 - sin2
86 const G4double pi_over_sigma = (rho2 * std
87 / (rho2 * s
88 const G4double Refl_pi = Refl_sigma * pi_o
89 theReflectivity = 0.5 * (Refl_sigma + Refl
90 G4double RoughAtten = 1;
91 if (fSurfaceRoughness > 0) {
92 G4double kiz = SinIncidentAngle * GamEne
93 G4double kjz = SinIncidentAngle * (1 - d
94 RoughAtten = G4Exp(-2 * kiz * kjz * fSur
95 theReflectivity *= RoughAtten;
96 }
97 if (GetVerboseLevel() > 1)
98 G4cout << std::left << std::setw(12) <<
99 << std::right << std::setw(4) <<
100 << " fSurfaceRoughness=" << G4Bes
101 << " RoughAtten=" << RoughAtten <
102 << " delta=" << delta << " beta="
103 << " Refl_pi=" << Refl_pi << " th
104 }
105 return theReflectivity;
106 }
107
108 //....oooOO0OOooo........oooOO0OOooo........oo
109
110 G4double G4XrayReflection::GetMeanFreePath(con
111 G4F
112 {
113 *condition = NotForced;
114 G4double GamEner = aTrack.GetDynamicParticle
115 if (GamEner < 30. * eV || GamEner > 30. * ke
116 return DBL_MAX; // do nothing below and a
117
118 if (GetVerboseLevel() > 2)
119 G4cout << std::left << std::setw(12) << __
120 << std::right << std::setw(4) << __
121 << " keV previousStepSize=" << prev
122 << " TrackLength=" << aTrack.GetTra
123 << G4endl;
124
125 G4double MeanFreePath = DBL_MAX; // by defa
126 G4VPhysicalVolume* Volume = aTrack.GetVolume
127 if (fLastVolume && Volume != fLastVolume &&
128 const G4Material* theLastMat = fLastVolume
129 const G4Material* theMat = Volume->GetLogi
130
131 G4double last_density = theLastMat->GetDen
132 G4double density = theMat->GetDensity();
133 if (density > last_density) { // density
134 G4Navigator* theNavigator =
135 G4TransportationManager::GetTransporta
136 G4bool valid = false;
137 G4ThreeVector theSurfaceNormal =
138 theNavigator->GetGlobalExitNormal(aTra
139 if (valid) fSurfaceNormal = theSurfaceNo
140 G4double SinIncidentAngle =
141 aTrack.GetDynamicParticle()->GetMoment
142 if (G4UniformRand() < Reflectivity(GamEn
143 MeanFreePath = 0;
144 }
145 G4ThreeVector Position = aTrack.GetPosit
146 const G4VSolid* LastSolid_Volume = fLast
147 if (GetVerboseLevel() > 1 && MeanFreePat
148 G4cout << std::left << std::setw(12) <
149 << std::right << std::setw(4) <
150 << " trigger reflection SinInci
151 << " at z=" << Position.getZ()
152 else if (GetVerboseLevel() > 2)
153 G4cout << std::left << std::setw(12) <
154 << std::right << std::setw(4) <
155 << " last logical volume name =
156 << " last logical volume materi
157 << " last density=" << last_den
158 << " logical volume name =" <<
159 << " logical volume material na
160 << " part/cm3 ? "
161 << " LastSolid_Volume->Inside(P
162 << " sin(IncidentAngle)=" << Si
163 << G4endl;
164 }
165 }
166 fLastVolume = Volume;
167 return MeanFreePath;
168 }
169
170 //....oooOO0OOooo........oooOO0OOooo........oo
171
172 G4VParticleChange* G4XrayReflection::PostStepD
173 {
174 aParticleChange.Initialize(aTrack); // copy
175 G4ThreeVector PhotDir = aTrack.GetDynamicPar
176 G4ThreeVector para_part = (PhotDir * fSurfac
177 G4ThreeVector photon_reflected = PhotDir - 2
178 if (GetVerboseLevel() > 1)
179 G4cout << std::left << std::setw(12) << __
180 << std::right << std::setw(4) << __
181 << " StepLength=" << aStep.GetStepL
182 << " photon_reflected=" << photon_r
183 << " aParticleChange.GetTrackStatus
184 << G4endl;
185
186 aParticleChange.ProposeTrackStatus(
187 fStopAndKill); // needed when working wit
188 // primary
189 auto ReflectedPhoton = new G4DynamicParticle
190
191 aParticleChange.AddSecondary(ReflectedPhoton
192 return &aParticleChange;
193 }
194
195 //....oooOO0OOooo........oooOO0OOooo........oo
196
197 void G4XrayReflection::BuildPhysicsTable(const
198 {
199 ProcessDescription(G4cout);
200
201 if (GetVerboseLevel() > 2)
202 G4cout << std::left << std::setw(12) << __
203 << std::right << std::setw(4) << __
204 << " is gamma=" << (&part == G4Gamm
205 << " fSurfaceRoughness=" << G4BestU
206 }
207
208 //....oooOO0OOooo........oooOO0OOooo........oo
209
210 void G4XrayReflection::ProcessDescription(std:
211 {
212 if (G4Threading::IsMasterThread())
213 out << '\n' << GetProcessName() << ": Gamm
214 }
215
216 //....oooOO0OOooo........oooOO0OOooo........oo
217
218 G4int G4XrayReflection::ReadHenkeXrayData(std:
219 std:
220 {
221 std::transform(ElName.begin(), ElName.end(),
222 ::tolower); // henke_physica
223 const G4String DataDir = G4EmParameters::Ins
224 const G4String InpFname = DataDir + ElName +
225 std::ifstream infile(InpFname);
226 if (!infile.is_open()) {
227 G4cout << "ReadHenkeXrayReflData " << InpF
228 return 1; // failure
229 }
230 std::vector<std::string> VarName(3);
231 infile >> VarName[0] >> VarName[1] >> VarNam
232 if (GetVerboseLevel())
233 G4cout << "ReadHenkeXrayData variable name
234 << VarName[2] << G4endl;
235 G4double E_eV_i, f1_i, f2_i;
236 Ephot.resize(0);
237 f1.resize(0);
238 f2.resize(0);
239 for (;;) {
240 infile >> E_eV_i >> f1_i >> f2_i;
241 if (infile.eof()) break;
242 Ephot.push_back(E_eV_i * eV);
243 f1.push_back(f1_i);
244 f2.push_back(f2_i);
245 }
246 return 0; // success
247 }
248
249 //....oooOO0OOooo........oooOO0OOooo........oo
250
251 void G4XrayReflection::SaveHenkeDataAsMaterial
252 {
253 // loop through the material table and load
254 // with Henke data used to calculate the ref
255 auto materialTable = G4Material::GetMaterial
256 for (auto a_material : *materialTable) {
257 auto N = a_material->GetTotNbOfAtomsPerVol
258 if (GetVerboseLevel() > 2)
259 if (GetVerboseLevel() > 2)
260 G4cout << std::left << std::setw(12) <
261 << std::right << std::setw(4) <
262 << " NbOfAtomsPerVolume()=" <<
263 << " NumberOfElements()=" << a_
264 // calculate the reflectivity from input d
265 // materials of a single element
266 if (a_material->GetNumberOfElements() == 1
267 G4double factor = N * CLHEP::classic_ele
268 std::vector<G4double> Ephot, f1, f2;
269 const G4Element* theElement = a_material
270 G4int iret = ReadHenkeXrayData(theElemen
271 if (iret) {
272 if (GetVerboseLevel() > 2)
273 G4cout << std::left << std::setw(12)
274 << std::right << std::setw(4)
275 << a_material->GetName() << "
276 }
277 else {
278 std::vector<G4double> RealIndex(Ephot.
279 for (std::size_t i = 0; i < Ephot.size
280 G4double lambda = CLHEP::twopi * CLH
281 G4double lambda_sqr = lambda * lambd
282 RealIndex[i] = fmax(0, factor * lamb
283 ImagIndex[i] = factor * lambda_sqr *
284 if (GetVerboseLevel() > 2)
285 G4cout << "Ephot=" << std::setw(10
286 << RealIndex[i] << " beta="
287 } // photon energy
288 G4MaterialPropertiesTable* proptab = a
289 if(proptab == nullptr) {
290 proptab = new G4MaterialPropertiesTa
291 a_material->SetMaterialPropertiesTab
292 }
293 proptab->AddProperty("REALRINDEX", Eph
294 proptab->AddProperty("IMAGINARYRINDEX"
295 if (GetVerboseLevel() > 2)
296 G4cout << std::left << std::setw(12)
297 << std::right << std::setw(4)
298 << " " << theElement->GetName
299 << " reflection data saved in
300 } // data found
301 }
302 }
303 }
304
305 //....oooOO0OOooo........oooOO0OOooo........oo
306
307 void G4XrayReflection::SetSurfaceRoughness(con
308 {
309 fSurfaceRoughness = value;
310 }
311
312 //....oooOO0OOooo........oooOO0OOooo........oo
313