| Geant4 Cross Reference |
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 //------------------ G4XrayReflection physics process -----------------------
27 //
28 // History:
29 // 14-09-23 H. Burkhardt, initial implementation
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 = 0;
46
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48
49 G4XrayReflection::G4XrayReflection(const G4String& processName,
50 G4ProcessType type) // Constructor
51 : G4VDiscreteProcess(processName, type)
52 {
53 SetProcessSubType(fGammaReflection);
54 SaveHenkeDataAsMaterialProperty();
55 }
56
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59 G4bool G4XrayReflection::IsApplicable(const G4ParticleDefinition& particle)
60 {
61 return (&particle == G4Gamma::Gamma()); // apply only to gamma
62 }
63
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65
66 G4double G4XrayReflection::Reflectivity(const G4double GamEner, const G4double SinIncidentAngle,
67 const G4Material* theMat) const
68 {
69 G4double theReflectivity = 0;
70 const G4MaterialPropertiesTable* theMatProp = theMat->GetMaterialPropertiesTable();
71 if (SinIncidentAngle < 0.9 && theMatProp != nullptr)
72 { // avoid perpendicular refl. at straight entry and require
73 // data available
74 G4MaterialPropertyVector* RealIndex = theMatProp->GetProperty(kREALRINDEX);
75 G4MaterialPropertyVector* ImagIndex = theMatProp->GetProperty(kIMAGINARYRINDEX);
76 if (nullptr == RealIndex || nullptr == ImagIndex) { return theReflectivity; }
77 const G4double delta = RealIndex->Value(GamEner);
78 const G4double beta = ImagIndex->Value(GamEner);
79 const G4double sin2 = std::pow(SinIncidentAngle, 2);
80 const G4double rho2 =
81 0.5 * (sin2 - 2 * delta + std::sqrt(std::pow(sin2 - 2 * delta, 2) + 4 * beta * beta));
82 const G4double rho = std::sqrt(rho2);
83 const G4double Refl_sigma = (rho2 * std::pow(SinIncidentAngle - rho, 2) + std::pow(beta, 2))
84 / (rho2 * std::pow(SinIncidentAngle + rho, 2) + std::pow(beta, 2));
85 const G4double coscot = std::sqrt(1 - sin2) / SinIncidentAngle;
86 const G4double pi_over_sigma = (rho2 * std::pow(rho - coscot, 2) + std::pow(beta, 2))
87 / (rho2 * std::pow(rho + coscot, 2) + std::pow(beta, 2));
88 const G4double Refl_pi = Refl_sigma * pi_over_sigma;
89 theReflectivity = 0.5 * (Refl_sigma + Refl_pi); // unpolarized
90 G4double RoughAtten = 1;
91 if (fSurfaceRoughness > 0) {
92 G4double kiz = SinIncidentAngle * GamEner / CLHEP::hbarc;
93 G4double kjz = SinIncidentAngle * (1 - delta) * GamEner / CLHEP::hbarc;
94 RoughAtten = G4Exp(-2 * kiz * kjz * fSurfaceRoughness * fSurfaceRoughness); // Nevot–Croce
95 theReflectivity *= RoughAtten;
96 }
97 if (GetVerboseLevel() > 1)
98 G4cout << std::left << std::setw(12) << __FILE__ << " " << __FUNCTION__ << " line "
99 << std::right << std::setw(4) << __LINE__ << " GamEner=" << GamEner
100 << " fSurfaceRoughness=" << G4BestUnit(fSurfaceRoughness, "Length")
101 << " RoughAtten=" << RoughAtten << " SinIncidentAngle=" << SinIncidentAngle
102 << " delta=" << delta << " beta=" << beta << " Refl_sigma=" << Refl_sigma
103 << " Refl_pi=" << Refl_pi << " theReflectivity=" << theReflectivity << G4endl;
104 }
105 return theReflectivity;
106 }
107
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109
110 G4double G4XrayReflection::GetMeanFreePath(const G4Track& aTrack, G4double previousStepSize,
111 G4ForceCondition* condition)
112 {
113 *condition = NotForced;
114 G4double GamEner = aTrack.GetDynamicParticle()->GetTotalEnergy();
115 if (GamEner < 30. * eV || GamEner > 30. * keV)
116 return DBL_MAX; // do nothing below and above the limits
117
118 if (GetVerboseLevel() > 2)
119 G4cout << std::left << std::setw(12) << __FILE__ << " " << __FUNCTION__ << " line "
120 << std::right << std::setw(4) << __LINE__ << " GamEner=" << GamEner / keV
121 << " keV previousStepSize=" << previousStepSize
122 << " TrackLength=" << aTrack.GetTrackLength() << " StepLength=" << aTrack.GetStepLength()
123 << G4endl;
124
125 G4double MeanFreePath = DBL_MAX; // by default no reflection
126 G4VPhysicalVolume* Volume = aTrack.GetVolume();
127 if (fLastVolume && Volume != fLastVolume && aTrack.GetTrackLength() > 0) { // at a boundary
128 const G4Material* theLastMat = fLastVolume->GetLogicalVolume()->GetMaterial();
129 const G4Material* theMat = Volume->GetLogicalVolume()->GetMaterial();
130
131 G4double last_density = theLastMat->GetDensity();
132 G4double density = theMat->GetDensity();
133 if (density > last_density) { // density has increased
134 G4Navigator* theNavigator =
135 G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
136 G4bool valid = false;
137 G4ThreeVector theSurfaceNormal =
138 theNavigator->GetGlobalExitNormal(aTrack.GetPosition(), &valid);
139 if (valid) fSurfaceNormal = theSurfaceNormal;
140 G4double SinIncidentAngle =
141 aTrack.GetDynamicParticle()->GetMomentumDirection() * fSurfaceNormal;
142 if (G4UniformRand() < Reflectivity(GamEner, SinIncidentAngle, theMat)) {
143 MeanFreePath = 0;
144 }
145 G4ThreeVector Position = aTrack.GetPosition(); // only for info
146 const G4VSolid* LastSolid_Volume = fLastVolume->GetLogicalVolume()->GetSolid(); // for info
147 if (GetVerboseLevel() > 1 && MeanFreePath == 0)
148 G4cout << std::left << std::setw(12) << __FILE__ << " " << __FUNCTION__ << " line "
149 << std::right << std::setw(4) << __LINE__
150 << " trigger reflection SinIncidentAngle=" << SinIncidentAngle
151 << " at z=" << Position.getZ() / meter << " m" << G4endl;
152 else if (GetVerboseLevel() > 2)
153 G4cout << std::left << std::setw(12) << __FILE__ << " " << __FUNCTION__ << " line "
154 << std::right << std::setw(4) << __LINE__ << " volume has changed "
155 << " last logical volume name =" << fLastVolume->GetLogicalVolume()->GetName()
156 << " last logical volume material name =" << theLastMat->GetName()
157 << " last density=" << last_density << " part/cm3 ? "
158 << " logical volume name =" << Volume->GetLogicalVolume()->GetName()
159 << " logical volume material name =" << theMat->GetName() << " density=" << density
160 << " part/cm3 ? "
161 << " LastSolid_Volume->Inside(Position)=" << LastSolid_Volume->Inside(Position)
162 << " sin(IncidentAngle)=" << SinIncidentAngle << " MeanFreePath=" << MeanFreePath
163 << G4endl;
164 }
165 }
166 fLastVolume = Volume;
167 return MeanFreePath;
168 }
169
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171
172 G4VParticleChange* G4XrayReflection::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
173 {
174 aParticleChange.Initialize(aTrack); // copy the current position to the changed particle
175 G4ThreeVector PhotDir = aTrack.GetDynamicParticle()->GetMomentumDirection();
176 G4ThreeVector para_part = (PhotDir * fSurfaceNormal) * fSurfaceNormal;
177 G4ThreeVector photon_reflected = PhotDir - 2 * para_part; // invert the parallel component
178 if (GetVerboseLevel() > 1)
179 G4cout << std::left << std::setw(12) << __FILE__ << " " << __FUNCTION__ << " line "
180 << std::right << std::setw(4) << __LINE__ << " fSurfaceNormal=" << fSurfaceNormal
181 << " StepLength=" << aStep.GetStepLength() << " PhotDir=" << PhotDir
182 << " photon_reflected=" << photon_reflected << " para_part=" << para_part
183 << " aParticleChange.GetTrackStatus()=" << aParticleChange.GetTrackStatus()
184 << G4endl;
185
186 aParticleChange.ProposeTrackStatus(
187 fStopAndKill); // needed when working with primary gamma to get rid of
188 // primary
189 auto ReflectedPhoton = new G4DynamicParticle(G4Gamma::Gamma(), photon_reflected,
190 aTrack.GetDynamicParticle()->GetTotalEnergy());
191 aParticleChange.AddSecondary(ReflectedPhoton);
192 return &aParticleChange;
193 }
194
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
196
197 void G4XrayReflection::BuildPhysicsTable(const G4ParticleDefinition& part)
198 {
199 ProcessDescription(G4cout);
200
201 if (GetVerboseLevel() > 2)
202 G4cout << std::left << std::setw(12) << __FILE__ << " " << __FUNCTION__ << " line "
203 << std::right << std::setw(4) << __LINE__
204 << " is gamma=" << (&part == G4Gamma::Definition())
205 << " fSurfaceRoughness=" << G4BestUnit(fSurfaceRoughness, "Length") << G4endl;
206 }
207
208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
209
210 void G4XrayReflection::ProcessDescription(std::ostream& out) const
211 {
212 if (G4Threading::IsMasterThread())
213 out << '\n' << GetProcessName() << ": Gamma specular reflection for energies > 30 eV.\n";
214 }
215
216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
217
218 G4int G4XrayReflection::ReadHenkeXrayData(std::string ElName, std::vector<G4double>& Ephot,
219 std::vector<G4double>& f1, std::vector<G4double>& f2)
220 {
221 std::transform(ElName.begin(), ElName.end(), ElName.begin(),
222 ::tolower); // henke_physical_reference uses lower case filanames
223 const G4String DataDir = G4EmParameters::Instance()->GetDirLEDATA() + "/XRayReflection_data/";
224 const G4String InpFname = DataDir + ElName + ".nff";
225 std::ifstream infile(InpFname);
226 if (!infile.is_open()) {
227 G4cout << "ReadHenkeXrayReflData " << InpFname << " not found" << G4endl;
228 return 1; // failure
229 }
230 std::vector<std::string> VarName(3);
231 infile >> VarName[0] >> VarName[1] >> VarName[2];
232 if (GetVerboseLevel())
233 G4cout << "ReadHenkeXrayData variable names " << VarName[0] << " " << VarName[1] << " "
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........oooOO0OOooo........oooOO0OOooo......
250
251 void G4XrayReflection::SaveHenkeDataAsMaterialProperty()
252 {
253 // loop through the material table and load set up MaterialPropertiesTable
254 // with Henke data used to calculate the reflection
255 auto materialTable = G4Material::GetMaterialTable();
256 for (auto a_material : *materialTable) {
257 auto N = a_material->GetTotNbOfAtomsPerVolume();
258 if (GetVerboseLevel() > 2)
259 if (GetVerboseLevel() > 2)
260 G4cout << std::left << std::setw(12) << __FILE__ << " " << __FUNCTION__ << " line "
261 << std::right << std::setw(4) << __LINE__ << " " << a_material->GetName()
262 << " NbOfAtomsPerVolume()=" << N
263 << " NumberOfElements()=" << a_material->GetNumberOfElements() << G4endl;
264 // calculate the reflectivity from input data. Implemented for dense
265 // materials of a single element
266 if (a_material->GetNumberOfElements() == 1 && a_material->GetDensity() > 1) {
267 G4double factor = N * CLHEP::classic_electr_radius / CLHEP::twopi;
268 std::vector<G4double> Ephot, f1, f2;
269 const G4Element* theElement = a_material->GetElement(0);
270 G4int iret = ReadHenkeXrayData(theElement->GetName(), Ephot, f1, f2);
271 if (iret) {
272 if (GetVerboseLevel() > 2)
273 G4cout << std::left << std::setw(12) << __FILE__ << " " << __FUNCTION__ << " line "
274 << std::right << std::setw(4) << __LINE__ << " no Henke data found for "
275 << a_material->GetName() << " " << theElement->GetName() << G4endl;
276 }
277 else {
278 std::vector<G4double> RealIndex(Ephot.size()), ImagIndex(Ephot.size());
279 for (std::size_t i = 0; i < Ephot.size(); ++i) {
280 G4double lambda = CLHEP::twopi * CLHEP::hbarc / Ephot[i];
281 G4double lambda_sqr = lambda * lambda;
282 RealIndex[i] = fmax(0, factor * lambda_sqr * f1[i]); // delta or 1-RealIndex
283 ImagIndex[i] = factor * lambda_sqr * f2[i]; // beta or -ImagIndex
284 if (GetVerboseLevel() > 2)
285 G4cout << "Ephot=" << std::setw(10) << Ephot[i] / eV << " eV delta=" << std::setw(10)
286 << RealIndex[i] << " beta=" << std::setw(10) << ImagIndex[i] << G4endl;
287 } // photon energy
288 G4MaterialPropertiesTable* proptab = a_material->GetMaterialPropertiesTable();
289 if(proptab == nullptr) {
290 proptab = new G4MaterialPropertiesTable();
291 a_material->SetMaterialPropertiesTable(proptab);
292 }
293 proptab->AddProperty("REALRINDEX", Ephot, RealIndex); // 1-RealIndex
294 proptab->AddProperty("IMAGINARYRINDEX", Ephot, ImagIndex);
295 if (GetVerboseLevel() > 2)
296 G4cout << std::left << std::setw(12) << __FILE__ << " " << __FUNCTION__ << " line "
297 << std::right << std::setw(4) << __LINE__ << " " << a_material->GetName()
298 << " " << theElement->GetName()
299 << " reflection data saved in PropertiesTable" << G4endl;
300 } // data found
301 }
302 }
303 }
304
305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
306
307 void G4XrayReflection::SetSurfaceRoughness(const G4double value)
308 {
309 fSurfaceRoughness = value;
310 }
311
312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
313