Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UCNMaterialPropertiesTable.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 //
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 <fstream>
46 
49 
50 #include "G4SystemOfUnits.hh"
51 #include "G4PhysicalConstants.hh"
52 
55 {
56  theMicroRoughnessTable = nullptr;
57  maxMicroRoughnessTable = nullptr;
58  theMicroRoughnessTransTable = nullptr;
59  maxMicroRoughnessTransTable = nullptr;
60 
61  theta_i_min = 0.*degree;
62  theta_i_max = 90.*degree;
63 
64  Emin = 0.e-9*eV;
65  Emax = 1000.e-9*eV;
66 
67  no_theta_i = 90;
68  noE = 100;
69 
70  theta_i_step = (theta_i_max-theta_i_min)/(no_theta_i-1);
71  E_step = (Emax-Emin)/(noE-1);
72 
73  b = 1*nm;
74  w = 30*nm;
75 
76  AngCut = 0.01*degree;
77 }
78 
80 {
81  if (theMicroRoughnessTable) delete theMicroRoughnessTable;
82  if (maxMicroRoughnessTable) delete maxMicroRoughnessTable;
83  if (theMicroRoughnessTransTable) delete theMicroRoughnessTransTable;
84  if (maxMicroRoughnessTransTable) delete maxMicroRoughnessTransTable;
85 }
86 
88 {
89  return theMicroRoughnessTable;
90 }
91 
93 {
94  return theMicroRoughnessTransTable;
95 }
96 
98  LoadMicroRoughnessTables(G4double* pMicroRoughnessTable,
99  G4double* pmaxMicroRoughnessTable,
100  G4double* pMicroRoughnessTransTable,
101  G4double* pmaxMicroRoughnessTransTable)
102 {
103  theMicroRoughnessTable = pMicroRoughnessTable;
104  maxMicroRoughnessTable = pmaxMicroRoughnessTable;
105  theMicroRoughnessTransTable = pMicroRoughnessTransTable;
106  maxMicroRoughnessTransTable = pmaxMicroRoughnessTransTable;
107 }
108 
110 {
111  G4int NEdim = 0;
112  G4int Nthetadim = 0;
113 
114  // Checks if the number of angles is available and stores it
115 
116  if (ConstPropertyExists("MR_NBTHETA"))
117  Nthetadim = G4int(GetConstProperty("MR_NBTHETA")+0.1);
118 
119  // Checks if the number of energies is available and stores it
120 
121  if (ConstPropertyExists("MR_NBE"))
122  NEdim = G4int(GetConstProperty("MR_NBE")+0.1);
123 
124  //G4cout << "thetadim: " << Nthetadim << " , Edim: " << NEdim << G4endl;
125 
126  // If both dimensions of the lookup-table are non-trivial:
127  // delete old tables if existing and allocate memory for new tables
128 
129  if (Nthetadim*NEdim > 0) {
130  if (theMicroRoughnessTable) delete theMicroRoughnessTable;
131  theMicroRoughnessTable = new G4double[Nthetadim*NEdim];
132  if (maxMicroRoughnessTable) delete maxMicroRoughnessTable;
133  maxMicroRoughnessTable = new G4double[Nthetadim*NEdim];
134  if (theMicroRoughnessTransTable) delete theMicroRoughnessTransTable;
135  theMicroRoughnessTransTable = new G4double[Nthetadim*NEdim];
136  if (maxMicroRoughnessTransTable) delete maxMicroRoughnessTransTable;
137  maxMicroRoughnessTransTable = new G4double[Nthetadim*NEdim];
138  }
139 }
140 
142 {
143 // Reads the parameters for the mr-probability computation from the
144 // corresponding material properties and stores it in the appropriate
145 // variables
146 
147  b = GetConstProperty("MR_RRMS");
148  G4double b2 = b*b;
149  w = GetConstProperty("MR_CORRLEN");
150  G4double w2 = w*w;
151 
152  no_theta_i = G4int(GetConstProperty("MR_NBTHETA")+0.1);
153  //G4cout << "no_theta: " << no_theta_i << G4endl;
154  noE = G4int(GetConstProperty("MR_NBE")+0.1);
155  //G4cout << "noE: " << noE << G4endl;
156 
157  theta_i_min = GetConstProperty("MR_THETAMIN");
158  theta_i_max = GetConstProperty("MR_THETAMAX");
159  Emin = GetConstProperty("MR_EMIN");
160  Emax = GetConstProperty("MR_EMAX");
161  G4int AngNoTheta = G4int(GetConstProperty("MR_ANGNOTHETA")+0.1);
162  G4int AngNoPhi = G4int(GetConstProperty("MR_ANGNOPHI")+0.1);
163  AngCut = GetConstProperty("MR_ANGCUT");
164 
165  // The Fermi potential was saved in neV by P.F.
166 
167  G4double fermipot = GetConstProperty("FERMIPOT")*(1.e-9*eV);
168 
169  //G4cout << "Fermipot: " << fermipot/(1.e-9*eV) << "neV" << G4endl;
170 
171  G4double theta_i, E;
172 
173  // Calculates the increment in theta_i in the lookup-table
174  theta_i_step = (theta_i_max-theta_i_min)/(no_theta_i-1);
175 
176  //G4cout << "theta_i_step: " << theta_i_step << G4endl;
177 
178  // Calculates the increment in energy in the lookup-table
179  E_step = (Emax-Emin)/(noE-1);
180 
181  // Runs the lookup-table memory allocation
183 
184  G4int counter = 0;
185 
186  //G4cout << "Calculating Microroughnesstables...";
187 
188  // Writes the mr-lookup-tables to files for immediate control
189 
190  std::ofstream dateir("MRrefl.dat",std::ios::out);
191  std::ofstream dateit("MRtrans.dat",std::ios::out);
192 
193  //G4cout << theMicroRoughnessTable << G4endl;
194 
195  for (theta_i=theta_i_min; theta_i<=theta_i_max+1e-6; theta_i+=theta_i_step) {
196  // Calculation for each cell in the lookup-table
197  for (E=Emin; E<=Emax; E+=E_step) {
198  *(theMicroRoughnessTable+counter) =
200  IntIplus(E, fermipot, theta_i, AngNoTheta, AngNoPhi,
201  b2, w2, maxMicroRoughnessTable+counter, AngCut);
202 
203  *(theMicroRoughnessTransTable+counter) =
205  IntIminus(E, fermipot, theta_i, AngNoTheta, AngNoPhi,
206  b2, w2, maxMicroRoughnessTransTable+counter,
207  AngCut);
208 
209  dateir << *(theMicroRoughnessTable+counter) << G4endl;
210  dateit << *(theMicroRoughnessTransTable+counter) << G4endl;
211 
212  counter++;
213 
214  //G4cout << counter << G4endl;
215  }
216  }
217 
218  dateir.close();
219  dateit.close();
220 
221  // Writes the mr-lookup-tables to files for immediate control
222 
223  std::ofstream dateic("MRcheck.dat",std::ios::out);
224  std::ofstream dateimr("MRmaxrefl.dat",std::ios::out);
225  std::ofstream dateimt("MRmaxtrans.dat",std::ios::out);
226 
227  for (theta_i=theta_i_min; theta_i<=theta_i_max+1e-6; theta_i+=theta_i_step) {
228  for (E=Emin; E<=Emax; E+=E_step) {
229 
230  // tests the GetXXProbability functions by writing the entries
231  // of the lookup tables to files
232 
233  dateic << GetMRIntProbability(theta_i,E) << G4endl;
234  dateimr << GetMRMaxProbability(theta_i,E) << G4endl;
235  dateimt << GetMRMaxTransProbability(theta_i,E) << G4endl;
236  }
237  }
238 
239  dateic.close();
240  dateimr.close();
241  dateimt.close();
242 }
243 
246 {
247  if (!theMicroRoughnessTable) {
248  G4cout << "Dont have theMicroRoughnessTable" << G4endl;
249  return 0.;
250  }
251 
252  // if theta_i or energy outside the range for which the lookup table is
253  // calculated, the probability is set to zero
254 
255  //G4cout << "theta_i: " << theta_i/degree << "degree"
256  // << " theta_i_min: " << theta_i_min/degree << "degree"
257  // << " theta_i_max: " << theta_i_max/degree << "degree"
258  // << " Energy: " << Energy/(1.e-9*eV) << "neV"
259  // << " Emin: " << Emin/(1.e-9*eV) << "neV"
260  // << " Emax: " << Emax/(1.e-9*eV) << "neV" << G4endl;
261 
262  if (theta_i<theta_i_min || theta_i>theta_i_max || Energy<Emin || Energy>Emax)
263  return 0.;
264 
265  // Determines the nearest cell in the lookup table which contains
266  // the probability
267 
268  G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5);
269  G4int E_pos = G4int((Energy-Emin)/E_step+0.5);
270 
271  // lookup table is onedimensional (1 row), energy is in rows,
272  // theta_i in columns
273 
274  //G4cout << "E_pos: " << E_pos << " theta_i_pos: " << theta_i_pos << G4endl;
275  //G4cout << "Probability: " << *(theMicroRoughnessTable+E_pos+theta_i_pos*(noE-1)) << G4endl;
276 
277  return *(theMicroRoughnessTable+E_pos+theta_i_pos*(noE - 1));
278 }
279 
282 {
283  if (!theMicroRoughnessTransTable) return 0.;
284 
285  // if theta_i or energy outside the range for which the lookup table
286  // is calculated, the probability is set to zero
287 
288  if (theta_i<theta_i_min || theta_i>theta_i_max || Energy<Emin || Energy>Emax)
289  return 0.;
290 
291  // Determines the nearest cell in the lookup table which contains
292  // the probability
293 
294  G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5);
295  G4int E_pos = G4int((Energy-Emin)/E_step+0.5);
296 
297  // lookup table is onedimensional (1 row), energy is in rows,
298  // theta_i in columns
299 
300  return *(theMicroRoughnessTransTable+E_pos+theta_i_pos*(noE - 1));
301 }
302 
305 {
306  if (!maxMicroRoughnessTable) return 0.;
307 
308  // if theta_i or energy outside the range for which the lookup table
309  // is calculated, the probability is set to zero
310 
311  if (theta_i<theta_i_min || theta_i>theta_i_max || Energy<Emin || Energy>Emax)
312  return 0.;
313 
314  // Determines the nearest cell in the lookup table which contains
315  // the probability
316 
317  G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5);
318  G4int E_pos = G4int((Energy-Emin)/E_step+0.5);
319 
320  // lookup table is onedimensional (1 row), energy is in rows,
321  // theta_i in columns
322 
323  return *(maxMicroRoughnessTable+E_pos+theta_i_pos*noE);
324 }
325 
328 {
329  if (maxMicroRoughnessTable) {
330 
331  if (theta_i<theta_i_min || theta_i>theta_i_max ||
332  Energy<Emin || Energy>Emax) {
333  } else {
334 
335  // Determines the nearest cell in the lookup table which contains
336  // the probability
337 
338  G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5);
339  G4int E_pos = G4int((Energy-Emin)/E_step+0.5);
340 
341  // lookup table is onedimensional (1 row), energy is in rows,
342  // theta_i in columns
343 
344  *(maxMicroRoughnessTable+E_pos+theta_i_pos*noE) = value;
345  }
346  }
347 }
348 
351 {
352  if (!maxMicroRoughnessTransTable) return 0.;
353 
354  // if theta_i or energy outside the range for which the lookup table
355  // is calculated, the probability is set to zero
356 
357  if (theta_i<theta_i_min || theta_i>theta_i_max || Energy<Emin || Energy>Emax)
358  return 0.;
359 
360  // Determines the nearest cell in the lookup table which contains
361  // the probability
362 
363  G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5);
364  G4int E_pos = G4int((Energy-Emin)/E_step+0.5);
365 
366  // lookup table is onedimensional (1 row), energy is in rows,
367  // theta_i in columns
368 
369  return *(maxMicroRoughnessTransTable+E_pos+theta_i_pos*noE);
370 }
371 
374 {
375  if (maxMicroRoughnessTransTable) {
376 
377  if (theta_i<theta_i_min || theta_i>theta_i_max ||
378  Energy<Emin || Energy>Emax) {
379  } else {
380 
381  // Determines the nearest cell in the lookup table which contains
382  // the probability
383 
384  G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5);
385  G4int E_pos = G4int((Energy-Emin)/E_step+0.5);
386 
387  // lookup table is onedimensional (1 row), energy is in rows,
388  // theta_i in columns
389 
390  *(maxMicroRoughnessTransTable+E_pos+theta_i_pos*noE) = value;
391  }
392  }
393 }
394 
397  G4double fermipot,
398  G4double theta_o, G4double phi_o)
399 {
401  ProbIplus(Energy, fermipot, theta_i, theta_o, phi_o, b, w, AngCut);
402 }
403 
406  G4double fermipot,
407  G4double theta_o, G4double phi_o)
408 {
410  ProbIminus(Energy, fermipot,theta_i, theta_o, phi_o, b, w, AngCut);
411 }
412 
414  G4double VFermi,
415  G4double theta_i)
416 {
417  G4double k = std::sqrt(2*neutron_mass_c2*E / hbarc_squared);
418  G4double k_l = std::sqrt(2*neutron_mass_c2*VFermi / hbarc_squared);
419 
420  //G4cout << " Energy: " << E/(1.e-9*eV) << "neV"
421  // << " VFermi: " << VFermi/(1.e-9*eV) << "neV"
422  // << " theta: " << theta_i/degree << "degree" << G4endl;
423 
424  //G4cout << " Conditions - 2*b*k*cos(theta): " << 2*b*k*cos(theta_i)
425  // << ", 2*b*kl: " << 2*b*k_l << G4endl;
426 
427  // see eq. 17 of the Steyerl paper
428 
429  if (2*b*k*std::cos(theta_i) < 1 && 2*b*k_l < 1) return true;
430  else return false;
431 }
432 
434  G4double VFermi,
435  G4double theta_i)
436 {
438  G4double k_l2 = 2*neutron_mass_c2*VFermi / hbarc_squared;
439 
440  if (E*(std::cos(theta_i)*std::cos(theta_i)) < VFermi) return false;
441 
442  G4double kS2 = k_l2 - k2;
443 
444  //G4cout << "Conditions; 2bk' cos(theta): " << 2*b*sqrt(kS2)*cos(theta_i) <<
445  // ", 2bk_l: " << 2*b*sqrt(k_l2) << G4endl;
446 
447  // see eq. 18 of the Steyerl paper
448 
449  if (2*b*std::sqrt(kS2)*std::cos(theta_i) < 1 && 2*b*std::sqrt(k_l2) < 1) return true;
450  else return false;
451 }
452 
455  G4int no_theta, G4int no_E,
456  G4double theta_min, G4double theta_max,
457  G4double E_min, G4double E_max,
458  G4int AngNoTheta, G4int AngNoPhi,
459  G4double AngularCut)
460 {
461  //G4cout << "Setting Microroughness Parameters...";
462 
463  // Removes an existing RMS roughness
464  if (ConstPropertyExists("MR_RRMS")) RemoveConstProperty("MR_RRMS");
465 
466  // Adds a new RMS roughness
467  AddConstProperty("MR_RRMS", bb);
468 
469  //G4cout << "b: " << bb << G4endl;
470 
471  // Removes an existing correlation length
472  if (ConstPropertyExists("MR_CORRLEN")) RemoveConstProperty("MR_CORRLEN");
473 
474  // Adds a new correlation length
475  AddConstProperty("MR_CORRLEN", ww);
476 
477  //G4cout << "w: " << ww << G4endl;
478 
479  // Removes an existing number of thetas
480  if (ConstPropertyExists("MR_NBTHETA")) RemoveConstProperty("MR_NBTHETA");
481 
482  // Adds a new number of thetas
483  AddConstProperty("MR_NBTHETA", (G4double)no_theta);
484 
485  //G4cout << "no_theta: " << no_theta << G4endl;
486 
487  // Removes an existing number of Energies
488  if (ConstPropertyExists("MR_NBE")) RemoveConstProperty("MR_NBE");
489 
490  // Adds a new number of Energies
491  AddConstProperty("MR_NBE", (G4double)no_E);
492 
493  //G4cout << "no_E: " << no_E << G4endl;
494 
495  // Removes an existing minimum theta
496  if (ConstPropertyExists("MR_THETAMIN")) RemoveConstProperty("MR_THETAMIN");
497 
498  // Adds a new minimum theta
499  AddConstProperty("MR_THETAMIN", theta_min);
500 
501  //G4cout << "theta_min: " << theta_min << G4endl;
502 
503  // Removes an existing maximum theta
504  if (ConstPropertyExists("MR_THETAMAX")) RemoveConstProperty("MR_THETAMAX");
505 
506  // Adds a new maximum theta
507  AddConstProperty("MR_THETAMAX", theta_max);
508 
509  //G4cout << "theta_max: " << theta_max << G4endl;
510 
511  // Removes an existing minimum energy
512  if (ConstPropertyExists("MR_EMIN")) RemoveConstProperty("MR_EMIN");
513 
514  // Adds a new minimum energy
515  AddConstProperty("MR_EMIN", E_min);
516 
517  //G4cout << "Emin: " << E_min << G4endl;
518 
519  // Removes an existing maximum energy
520  if (ConstPropertyExists("MR_EMAX")) RemoveConstProperty("MR_EMAX");
521 
522  // Adds a new maximum energy
523  AddConstProperty("MR_EMAX", E_max);
524 
525  //G4cout << "Emax: " << E_max << G4endl;
526 
527  // Removes an existing Theta angle number
528  if(ConstPropertyExists("MR_ANGNOTHETA"))RemoveConstProperty("MR_ANGNOTHETA");
529 
530  // Adds a new Theta angle number
531  AddConstProperty("MR_ANGNOTHETA", (G4double)AngNoTheta);
532 
533  //G4cout << "AngNoTheta: " << AngNoTheta << G4endl;
534 
535  // Removes an existing Phi angle number
536  if (ConstPropertyExists("MR_ANGNOPHI")) RemoveConstProperty("MR_ANGNOPHI");
537 
538  // Adds a new Phi angle number
539  AddConstProperty("MR_ANGNOPHI", (G4double)AngNoPhi);
540 
541  //G4cout << "AngNoPhi: " << AngNoPhi << G4endl;
542 
543  // Removes an existing angular cut
544  if (ConstPropertyExists("MR_ANGCUT")) RemoveConstProperty("MR_ANGCUT");
545 
546  // Adds a new angle number
547  AddConstProperty("MR_ANGCUT", AngularCut);
548 
549  //G4cout << "AngularCut: " << AngularCut/degree << "degree" << G4endl;
550 
551  // Starts the lookup table calculation
552 
554 
555  //G4cout << " *** DONE! ***" << G4endl;
556 }
G4double GetMRMaxTransProbability(G4double, G4double)
G4double GetMRProbability(G4double, G4double, G4double, G4double, G4double)
int G4int
Definition: G4Types.hh:78
void SetMRMaxProbability(G4double, G4double, G4double)
static ulg bb
Definition: csz_inflate.cc:344
void LoadMicroRoughnessTables(G4double *, G4double *, G4double *, G4double *)
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
void RemoveConstProperty(const char *key)
G4double GetMRTransProbability(G4double, G4double, G4double, G4double, G4double)
static constexpr double degree
Definition: G4SIunits.hh:144
bool G4bool
Definition: G4Types.hh:79
static constexpr double eV
Definition: G4SIunits.hh:215
static G4UCNMicroRoughnessHelper * GetInstance()
G4bool TransConditionsValid(G4double E, G4double VFermi, G4double theta_i)
void AddConstProperty(const char *key, G4double PropertyValue)
float neutron_mass_c2
Definition: hepunit.py:276
G4bool ConstPropertyExists(const char *key) const
G4double GetMRIntProbability(G4double, G4double)
static constexpr double nm
Definition: G4SIunits.hh:112
G4double GetMRIntTransProbability(G4double, G4double)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4bool ConditionsValid(G4double E, G4double VFermi, G4double theta_i)
void SetMicroRoughnessParameters(G4double, G4double, G4int, G4int, G4double, G4double, G4double, G4double, G4int, G4int, G4double)
G4double GetMRMaxProbability(G4double, G4double)
void SetMRMaxTransProbability(G4double, G4double, G4double)
G4double GetConstProperty(const char *key) const