| 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 //
27 //
28 //
29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30 //
31 // G4UCNMaterialPropertiesTable
32 // ----------------------------
33 //
34 // Derives from G4MaterialPropertiesTable and adds a double pointer to the
35 // MicroRoughnessTable
36 //
37 // This file includes the initialization and calculation of the mr-lookup
38 // tables. It also provides the functions to read from these tables and to
39 // calculate the probability for certain angles.
40 //
41 // For a closer description see the header file
42 //
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44
45 #include "G4UCNMaterialPropertiesTable.hh"
46
47 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4UCNMicroRoughnessHelper.hh"
50
51 #include <fstream>
52
53 G4UCNMaterialPropertiesTable::G4UCNMaterialPropertiesTable()
54 {
55 theMicroRoughnessTable = nullptr;
56 maxMicroRoughnessTable = nullptr;
57 theMicroRoughnessTransTable = nullptr;
58 maxMicroRoughnessTransTable = nullptr;
59
60 theta_i_min = 0. * degree;
61 theta_i_max = 90. * degree;
62
63 Emin = 0.e-9 * eV;
64 Emax = 1000.e-9 * eV;
65
66 no_theta_i = 90;
67 noE = 100;
68
69 theta_i_step = (theta_i_max - theta_i_min) / (no_theta_i - 1);
70 E_step = (Emax - Emin) / (noE - 1);
71
72 b = 1 * nm;
73 w = 30 * nm;
74
75 AngCut = 0.01 * degree;
76 }
77
78 G4UCNMaterialPropertiesTable::~G4UCNMaterialPropertiesTable()
79 {
80 delete theMicroRoughnessTable;
81 delete maxMicroRoughnessTable;
82 delete theMicroRoughnessTransTable;
83 delete maxMicroRoughnessTransTable;
84 }
85
86 G4double* G4UCNMaterialPropertiesTable::GetMicroRoughnessTable() { return theMicroRoughnessTable; }
87
88 G4double* G4UCNMaterialPropertiesTable::GetMicroRoughnessTransTable()
89 {
90 return theMicroRoughnessTransTable;
91 }
92
93 void G4UCNMaterialPropertiesTable::LoadMicroRoughnessTables(G4double* pMicroRoughnessTable,
94 G4double* pmaxMicroRoughnessTable, G4double* pMicroRoughnessTransTable,
95 G4double* pmaxMicroRoughnessTransTable)
96 {
97 theMicroRoughnessTable = pMicroRoughnessTable;
98 maxMicroRoughnessTable = pmaxMicroRoughnessTable;
99 theMicroRoughnessTransTable = pMicroRoughnessTransTable;
100 maxMicroRoughnessTransTable = pmaxMicroRoughnessTransTable;
101 }
102
103 void G4UCNMaterialPropertiesTable::InitMicroRoughnessTables()
104 {
105 G4int NEdim = 0;
106 G4int Nthetadim = 0;
107
108 // Checks if the number of angles is available and stores it
109
110 if (ConstPropertyExists("MR_NBTHETA")) {
111 Nthetadim = G4int(GetConstProperty("MR_NBTHETA") + 0.1);
112 }
113
114 // Checks if the number of energies is available and stores it
115
116 if (ConstPropertyExists("MR_NBE")) {
117 NEdim = G4int(GetConstProperty("MR_NBE") + 0.1);
118 }
119
120 // If both dimensions of the lookup-table are non-trivial:
121 // delete old tables if existing and allocate memory for new tables
122
123 if (Nthetadim * NEdim > 0) {
124 delete theMicroRoughnessTable;
125 theMicroRoughnessTable = new G4double[Nthetadim * NEdim];
126 delete maxMicroRoughnessTable;
127 maxMicroRoughnessTable = new G4double[Nthetadim * NEdim];
128 delete theMicroRoughnessTransTable;
129 theMicroRoughnessTransTable = new G4double[Nthetadim * NEdim];
130 delete maxMicroRoughnessTransTable;
131 maxMicroRoughnessTransTable = new G4double[Nthetadim * NEdim];
132 }
133 }
134
135 void G4UCNMaterialPropertiesTable::ComputeMicroRoughnessTables()
136 {
137 // Reads the parameters for the mr-probability computation from the
138 // corresponding material properties and stores it in the appropriate
139 // variables
140
141 b = GetConstProperty("MR_RRMS");
142 G4double b2 = b * b;
143 w = GetConstProperty("MR_CORRLEN");
144 G4double w2 = w * w;
145
146 no_theta_i = G4int(GetConstProperty("MR_NBTHETA") + 0.1);
147 noE = G4int(GetConstProperty("MR_NBE") + 0.1);
148
149 theta_i_min = GetConstProperty("MR_THETAMIN");
150 theta_i_max = GetConstProperty("MR_THETAMAX");
151 Emin = GetConstProperty("MR_EMIN");
152 Emax = GetConstProperty("MR_EMAX");
153 auto AngNoTheta = G4int(GetConstProperty("MR_ANGNOTHETA") + 0.1);
154 auto AngNoPhi = G4int(GetConstProperty("MR_ANGNOPHI") + 0.1);
155 AngCut = GetConstProperty("MR_ANGCUT");
156
157 // The Fermi potential was saved in neV by P.F.
158
159 G4double fermipot = GetConstProperty("FERMIPOT") * (1.e-9 * eV);
160
161 G4double theta_i, E;
162
163 // Calculates the increment in theta_i in the lookup-table
164 theta_i_step = (theta_i_max - theta_i_min) / (no_theta_i - 1);
165
166 // Calculates the increment in energy in the lookup-table
167 E_step = (Emax - Emin) / (noE - 1);
168
169 // Runs the lookup-table memory allocation
170 InitMicroRoughnessTables();
171
172 G4int counter = 0;
173
174 // Writes the mr-lookup-tables to files for immediate control
175 std::ofstream dateir("MRrefl.dat", std::ios::out);
176 std::ofstream dateit("MRtrans.dat", std::ios::out);
177
178 // G4cout << theMicroRoughnessTable << G4endl;
179
180 for (theta_i = theta_i_min; theta_i <= theta_i_max + 1e-6; theta_i += theta_i_step) {
181 // Calculation for each cell in the lookup-table
182 for (E = Emin; E <= Emax; E += E_step) {
183 *(theMicroRoughnessTable + counter) = G4UCNMicroRoughnessHelper::GetInstance()->IntIplus(E,
184 fermipot, theta_i, AngNoTheta, AngNoPhi, b2, w2, maxMicroRoughnessTable + counter, AngCut);
185
186 *(theMicroRoughnessTransTable + counter) =
187 G4UCNMicroRoughnessHelper::GetInstance()->IntIminus(E, fermipot, theta_i, AngNoTheta,
188 AngNoPhi, b2, w2, maxMicroRoughnessTransTable + counter, AngCut);
189
190 dateir << *(theMicroRoughnessTable + counter) << G4endl;
191 dateit << *(theMicroRoughnessTransTable + counter) << G4endl;
192
193 counter++;
194 }
195 }
196
197 dateir.close();
198 dateit.close();
199
200 // Writes the mr-lookup-tables to files for immediate control
201
202 std::ofstream dateic("MRcheck.dat", std::ios::out);
203 std::ofstream dateimr("MRmaxrefl.dat", std::ios::out);
204 std::ofstream dateimt("MRmaxtrans.dat", std::ios::out);
205
206 for (theta_i = theta_i_min; theta_i <= theta_i_max + 1e-6; theta_i += theta_i_step) {
207 for (E = Emin; E <= Emax; E += E_step) {
208 // tests the GetXXProbability functions by writing the entries
209 // of the lookup tables to files
210
211 dateic << GetMRIntProbability(theta_i, E) << G4endl;
212 dateimr << GetMRMaxProbability(theta_i, E) << G4endl;
213 dateimt << GetMRMaxTransProbability(theta_i, E) << G4endl;
214 }
215 }
216
217 dateic.close();
218 dateimr.close();
219 dateimt.close();
220 }
221
222 G4double G4UCNMaterialPropertiesTable::GetMRIntProbability(G4double theta_i, G4double Energy)
223 {
224 if (theMicroRoughnessTable == nullptr) {
225 G4cout << "Do not have theMicroRoughnessTable" << G4endl;
226 return 0.;
227 }
228
229 // if theta_i or energy outside the range for which the lookup table is
230 // calculated, the probability is set to zero
231 if (theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin || Energy > Emax) {
232 return 0.;
233 }
234
235 // Determines the nearest cell in the lookup table which contains
236 // the probability
237
238 auto theta_i_pos = G4int((theta_i - theta_i_min) / theta_i_step + 0.5);
239 auto E_pos = G4int((Energy - Emin) / E_step + 0.5);
240
241 // lookup table is onedimensional (1 row), energy is in rows,
242 // theta_i in columns
243 return *(theMicroRoughnessTable + E_pos + theta_i_pos * (noE - 1));
244 }
245
246 G4double G4UCNMaterialPropertiesTable::GetMRIntTransProbability(G4double theta_i, G4double Energy)
247 {
248 if (theMicroRoughnessTransTable == nullptr) {
249 return 0.;
250 }
251
252 // if theta_i or energy outside the range for which the lookup table
253 // is calculated, the probability is set to zero
254
255 if (theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin || Energy > Emax) {
256 return 0.;
257 }
258
259 // Determines the nearest cell in the lookup table which contains
260 // the probability
261
262 auto theta_i_pos = G4int((theta_i - theta_i_min) / theta_i_step + 0.5);
263 auto E_pos = G4int((Energy - Emin) / E_step + 0.5);
264
265 // lookup table is onedimensional (1 row), energy is in rows,
266 // theta_i in columns
267
268 return *(theMicroRoughnessTransTable + E_pos + theta_i_pos * (noE - 1));
269 }
270
271 G4double G4UCNMaterialPropertiesTable::GetMRMaxProbability(G4double theta_i, G4double Energy)
272 {
273 if (maxMicroRoughnessTable == nullptr) {
274 return 0.;
275 }
276
277 // if theta_i or energy outside the range for which the lookup table
278 // is calculated, the probability is set to zero
279
280 if (theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin || Energy > Emax) {
281 return 0.;
282 }
283
284 // Determines the nearest cell in the lookup table which contains
285 // the probability
286
287 auto theta_i_pos = G4int((theta_i - theta_i_min) / theta_i_step + 0.5);
288 auto E_pos = G4int((Energy - Emin) / E_step + 0.5);
289
290 // lookup table is onedimensional (1 row), energy is in rows,
291 // theta_i in columns
292
293 return *(maxMicroRoughnessTable + E_pos + theta_i_pos * noE);
294 }
295
296 void G4UCNMaterialPropertiesTable::SetMRMaxProbability(
297 G4double theta_i, G4double Energy, G4double value)
298 {
299 if (maxMicroRoughnessTable != nullptr) {
300 if (theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin || Energy > Emax) {
301 }
302 else {
303 // Determines the nearest cell in the lookup table which contains
304 // the probability
305
306 auto theta_i_pos = G4int((theta_i - theta_i_min) / theta_i_step + 0.5);
307 auto E_pos = G4int((Energy - Emin) / E_step + 0.5);
308
309 // lookup table is onedimensional (1 row), energy is in rows,
310 // theta_i in columns
311
312 *(maxMicroRoughnessTable + E_pos + theta_i_pos * noE) = value;
313 }
314 }
315 }
316
317 G4double G4UCNMaterialPropertiesTable::GetMRMaxTransProbability(G4double theta_i, G4double Energy)
318 {
319 if (maxMicroRoughnessTransTable == nullptr) {
320 return 0.;
321 }
322
323 // if theta_i or energy outside the range for which the lookup table
324 // is calculated, the probability is set to zero
325
326 if (theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin || Energy > Emax) {
327 return 0.;
328 }
329
330 // Determines the nearest cell in the lookup table which contains
331 // the probability
332
333 auto theta_i_pos = G4int((theta_i - theta_i_min) / theta_i_step + 0.5);
334 auto E_pos = G4int((Energy - Emin) / E_step + 0.5);
335
336 // lookup table is onedimensional (1 row), energy is in rows,
337 // theta_i in columns
338
339 return *(maxMicroRoughnessTransTable + E_pos + theta_i_pos * noE);
340 }
341
342 void G4UCNMaterialPropertiesTable::SetMRMaxTransProbability(
343 G4double theta_i, G4double Energy, G4double value)
344 {
345 if (maxMicroRoughnessTransTable != nullptr) {
346 if (theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin || Energy > Emax) {
347 }
348 else {
349 // Determines the nearest cell in the lookup table which contains
350 // the probability
351
352 auto theta_i_pos = G4int((theta_i - theta_i_min) / theta_i_step + 0.5);
353 auto E_pos = G4int((Energy - Emin) / E_step + 0.5);
354
355 // lookup table is onedimensional (1 row), energy is in rows,
356 // theta_i in columns
357
358 *(maxMicroRoughnessTransTable + E_pos + theta_i_pos * noE) = value;
359 }
360 }
361 }
362
363 G4double G4UCNMaterialPropertiesTable::GetMRProbability(
364 G4double theta_i, G4double Energy, G4double fermipot, G4double theta_o, G4double phi_o)
365 {
366 return G4UCNMicroRoughnessHelper::GetInstance()->ProbIplus(
367 Energy, fermipot, theta_i, theta_o, phi_o, b, w, AngCut);
368 }
369
370 G4double G4UCNMaterialPropertiesTable::GetMRTransProbability(
371 G4double theta_i, G4double Energy, G4double fermipot, G4double theta_o, G4double phi_o)
372 {
373 return G4UCNMicroRoughnessHelper::GetInstance()->ProbIminus(
374 Energy, fermipot, theta_i, theta_o, phi_o, b, w, AngCut);
375 }
376
377 G4bool G4UCNMaterialPropertiesTable::ConditionsValid(G4double E, G4double VFermi, G4double theta_i)
378 {
379 G4double k = std::sqrt(2 * neutron_mass_c2 * E / hbarc_squared);
380 G4double k_l = std::sqrt(2 * neutron_mass_c2 * VFermi / hbarc_squared);
381
382 // see eq. 17 of the Steyerl paper
383 return 2 * b * k * std::cos(theta_i) < 1 && 2 * b * k_l < 1;
384 }
385
386 G4bool G4UCNMaterialPropertiesTable::TransConditionsValid(
387 G4double E, G4double VFermi, G4double theta_i)
388 {
389 G4double k2 = 2 * neutron_mass_c2 * E / hbarc_squared;
390 G4double k_l2 = 2 * neutron_mass_c2 * VFermi / hbarc_squared;
391
392 if (E * (std::cos(theta_i) * std::cos(theta_i)) < VFermi) {
393 return false;
394 }
395
396 G4double kS2 = k_l2 - k2;
397
398 // see eq. 18 of the Steyerl paper
399 return 2 * b * std::sqrt(kS2) * std::cos(theta_i) < 1 && 2 * b * std::sqrt(k_l2) < 1;
400 }
401
402 void G4UCNMaterialPropertiesTable::SetMicroRoughnessParameters(G4double ww, G4double bb,
403 G4int no_theta, G4int no_E, G4double theta_min, G4double theta_max, G4double E_min,
404 G4double E_max, G4int AngNoTheta, G4int AngNoPhi, G4double AngularCut)
405 {
406 // Removes an existing RMS roughness
407 if (ConstPropertyExists("MR_RRMS")) {
408 RemoveConstProperty("MR_RRMS");
409 }
410
411 // Adds a new RMS roughness
412 AddConstProperty("MR_RRMS", bb);
413
414 // Removes an existing correlation length
415 if (ConstPropertyExists("MR_CORRLEN")) {
416 RemoveConstProperty("MR_CORRLEN");
417 }
418
419 // Adds a new correlation length
420 AddConstProperty("MR_CORRLEN", ww);
421
422 // Removes an existing number of thetas
423 if (ConstPropertyExists("MR_NBTHETA")) {
424 RemoveConstProperty("MR_NBTHETA");
425 }
426
427 // Adds a new number of thetas
428 AddConstProperty("MR_NBTHETA", (G4double)no_theta);
429
430 // Removes an existing number of Energies
431 if (ConstPropertyExists("MR_NBE")) {
432 RemoveConstProperty("MR_NBE");
433 }
434
435 // Adds a new number of Energies
436 AddConstProperty("MR_NBE", (G4double)no_E);
437
438 // Removes an existing minimum theta
439 if (ConstPropertyExists("MR_THETAMIN")) {
440 RemoveConstProperty("MR_THETAMIN");
441 }
442
443 // Adds a new minimum theta
444 AddConstProperty("MR_THETAMIN", theta_min);
445
446 // Removes an existing maximum theta
447 if (ConstPropertyExists("MR_THETAMAX")) {
448 RemoveConstProperty("MR_THETAMAX");
449 }
450
451 // Adds a new maximum theta
452 AddConstProperty("MR_THETAMAX", theta_max);
453
454 // Removes an existing minimum energy
455 if (ConstPropertyExists("MR_EMIN")) {
456 RemoveConstProperty("MR_EMIN");
457 }
458
459 // Adds a new minimum energy
460 AddConstProperty("MR_EMIN", E_min);
461
462 // Removes an existing maximum energy
463 if (ConstPropertyExists("MR_EMAX")) {
464 RemoveConstProperty("MR_EMAX");
465 }
466
467 // Adds a new maximum energy
468 AddConstProperty("MR_EMAX", E_max);
469
470 // Removes an existing Theta angle number
471 if (ConstPropertyExists("MR_ANGNOTHETA")) {
472 RemoveConstProperty("MR_ANGNOTHETA");
473 }
474
475 // Adds a new Theta angle number
476 AddConstProperty("MR_ANGNOTHETA", (G4double)AngNoTheta);
477
478 // Removes an existing Phi angle number
479 if (ConstPropertyExists("MR_ANGNOPHI")) {
480 RemoveConstProperty("MR_ANGNOPHI");
481 }
482
483 // Adds a new Phi angle number
484 AddConstProperty("MR_ANGNOPHI", (G4double)AngNoPhi);
485
486 // Removes an existing angular cut
487 if (ConstPropertyExists("MR_ANGCUT")) {
488 RemoveConstProperty("MR_ANGCUT");
489 }
490
491 // Adds a new angle number
492 AddConstProperty("MR_ANGCUT", AngularCut);
493
494 // Starts the lookup table calculation
495 ComputeMicroRoughnessTables();
496 }
497