| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // 27 //
28 // 28 //
29 //....oooOO0OOooo........oooOO0OOooo........oo 29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30 // 30 //
31 // G4UCNMaterialPropertiesTable 31 // G4UCNMaterialPropertiesTable
32 // ---------------------------- 32 // ----------------------------
33 // 33 //
34 // Derives from G4MaterialPropertiesTable and 34 // Derives from G4MaterialPropertiesTable and adds a double pointer to the
35 // MicroRoughnessTable 35 // MicroRoughnessTable
36 // 36 //
37 // This file includes the initialization and c 37 // This file includes the initialization and calculation of the mr-lookup
38 // tables. It also provides the functions to r 38 // tables. It also provides the functions to read from these tables and to
39 // calculate the probability for certain angle 39 // calculate the probability for certain angles.
40 // 40 //
41 // For a closer description see the header fil 41 // For a closer description see the header file
42 // 42 //
43 //....oooOO0OOooo........oooOO0OOooo........oo 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 44
45 #include "G4UCNMaterialPropertiesTable.hh" << 45 #include <fstream>
46 46
47 #include "G4PhysicalConstants.hh" << 47 #include "G4UCNMaterialPropertiesTable.hh"
48 #include "G4SystemOfUnits.hh" <<
49 #include "G4UCNMicroRoughnessHelper.hh" 48 #include "G4UCNMicroRoughnessHelper.hh"
50 49
51 #include <fstream> << 50 #include "G4SystemOfUnits.hh"
>> 51 #include "G4PhysicalConstants.hh"
52 52
53 G4UCNMaterialPropertiesTable::G4UCNMaterialPro 53 G4UCNMaterialPropertiesTable::G4UCNMaterialPropertiesTable()
>> 54 : G4MaterialPropertiesTable()
54 { 55 {
55 theMicroRoughnessTable = nullptr; 56 theMicroRoughnessTable = nullptr;
56 maxMicroRoughnessTable = nullptr; 57 maxMicroRoughnessTable = nullptr;
57 theMicroRoughnessTransTable = nullptr; 58 theMicroRoughnessTransTable = nullptr;
58 maxMicroRoughnessTransTable = nullptr; 59 maxMicroRoughnessTransTable = nullptr;
59 60
60 theta_i_min = 0. * degree; << 61 theta_i_min = 0.*degree;
61 theta_i_max = 90. * degree; << 62 theta_i_max = 90.*degree;
62 63
63 Emin = 0.e-9 * eV; << 64 Emin = 0.e-9*eV;
64 Emax = 1000.e-9 * eV; << 65 Emax = 1000.e-9*eV;
65 66
66 no_theta_i = 90; 67 no_theta_i = 90;
67 noE = 100; 68 noE = 100;
68 69
69 theta_i_step = (theta_i_max - theta_i_min) / << 70 theta_i_step = (theta_i_max-theta_i_min)/(no_theta_i-1);
70 E_step = (Emax - Emin) / (noE - 1); << 71 E_step = (Emax-Emin)/(noE-1);
71 72
72 b = 1 * nm; << 73 b = 1*nm;
73 w = 30 * nm; << 74 w = 30*nm;
74 75
75 AngCut = 0.01 * degree; << 76 AngCut = 0.01*degree;
76 } 77 }
77 78
78 G4UCNMaterialPropertiesTable::~G4UCNMaterialPr 79 G4UCNMaterialPropertiesTable::~G4UCNMaterialPropertiesTable()
79 { 80 {
80 delete theMicroRoughnessTable; 81 delete theMicroRoughnessTable;
81 delete maxMicroRoughnessTable; 82 delete maxMicroRoughnessTable;
82 delete theMicroRoughnessTransTable; 83 delete theMicroRoughnessTransTable;
83 delete maxMicroRoughnessTransTable; 84 delete maxMicroRoughnessTransTable;
84 } 85 }
85 86
86 G4double* G4UCNMaterialPropertiesTable::GetMic << 87 G4double* G4UCNMaterialPropertiesTable::GetMicroRoughnessTable ()
>> 88 {
>> 89 return theMicroRoughnessTable;
>> 90 }
87 91
88 G4double* G4UCNMaterialPropertiesTable::GetMic << 92 G4double* G4UCNMaterialPropertiesTable::GetMicroRoughnessTransTable ()
89 { 93 {
90 return theMicroRoughnessTransTable; 94 return theMicroRoughnessTransTable;
91 } 95 }
92 96
93 void G4UCNMaterialPropertiesTable::LoadMicroRo << 97 void G4UCNMaterialPropertiesTable::
94 G4double* pmaxMicroRoughnessTable, G4double* << 98 LoadMicroRoughnessTables(G4double* pMicroRoughnessTable,
95 G4double* pmaxMicroRoughnessTransTable) << 99 G4double* pmaxMicroRoughnessTable,
>> 100 G4double* pMicroRoughnessTransTable,
>> 101 G4double* pmaxMicroRoughnessTransTable)
96 { 102 {
97 theMicroRoughnessTable = pMicroRoughnessTabl << 103 theMicroRoughnessTable = pMicroRoughnessTable;
98 maxMicroRoughnessTable = pmaxMicroRoughnessT << 104 maxMicroRoughnessTable = pmaxMicroRoughnessTable;
99 theMicroRoughnessTransTable = pMicroRoughnes 105 theMicroRoughnessTransTable = pMicroRoughnessTransTable;
100 maxMicroRoughnessTransTable = pmaxMicroRough 106 maxMicroRoughnessTransTable = pmaxMicroRoughnessTransTable;
101 } 107 }
102 108
103 void G4UCNMaterialPropertiesTable::InitMicroRo 109 void G4UCNMaterialPropertiesTable::InitMicroRoughnessTables()
104 { 110 {
105 G4int NEdim = 0; << 111 G4int NEdim = 0;
106 G4int Nthetadim = 0; 112 G4int Nthetadim = 0;
107 113
108 // Checks if the number of angles is availab 114 // Checks if the number of angles is available and stores it
109 115
110 if (ConstPropertyExists("MR_NBTHETA")) { << 116 if(ConstPropertyExists("MR_NBTHETA"))
>> 117 {
111 Nthetadim = G4int(GetConstProperty("MR_NBT 118 Nthetadim = G4int(GetConstProperty("MR_NBTHETA") + 0.1);
112 } 119 }
113 120
114 // Checks if the number of energies is avail 121 // Checks if the number of energies is available and stores it
115 122
116 if (ConstPropertyExists("MR_NBE")) { << 123 if(ConstPropertyExists("MR_NBE"))
>> 124 {
117 NEdim = G4int(GetConstProperty("MR_NBE") + 125 NEdim = G4int(GetConstProperty("MR_NBE") + 0.1);
118 } 126 }
119 127
>> 128 //G4cout << "thetadim: " << Nthetadim << " , Edim: " << NEdim << G4endl;
>> 129
120 // If both dimensions of the lookup-table ar 130 // If both dimensions of the lookup-table are non-trivial:
121 // delete old tables if existing and allocat 131 // delete old tables if existing and allocate memory for new tables
122 132
123 if (Nthetadim * NEdim > 0) { << 133 if (Nthetadim*NEdim > 0) {
124 delete theMicroRoughnessTable; 134 delete theMicroRoughnessTable;
125 theMicroRoughnessTable = new G4double[Nthe 135 theMicroRoughnessTable = new G4double[Nthetadim * NEdim];
126 delete maxMicroRoughnessTable; 136 delete maxMicroRoughnessTable;
127 maxMicroRoughnessTable = new G4double[Nthe 137 maxMicroRoughnessTable = new G4double[Nthetadim * NEdim];
128 delete theMicroRoughnessTransTable; 138 delete theMicroRoughnessTransTable;
129 theMicroRoughnessTransTable = new G4double 139 theMicroRoughnessTransTable = new G4double[Nthetadim * NEdim];
130 delete maxMicroRoughnessTransTable; 140 delete maxMicroRoughnessTransTable;
131 maxMicroRoughnessTransTable = new G4double 141 maxMicroRoughnessTransTable = new G4double[Nthetadim * NEdim];
132 } 142 }
133 } 143 }
134 144
135 void G4UCNMaterialPropertiesTable::ComputeMicr 145 void G4UCNMaterialPropertiesTable::ComputeMicroRoughnessTables()
136 { 146 {
137 // Reads the parameters for the mr-probabili << 147 // Reads the parameters for the mr-probability computation from the
138 // corresponding material properties and sto << 148 // corresponding material properties and stores it in the appropriate
139 // variables << 149 // variables
140 150
141 b = GetConstProperty("MR_RRMS"); 151 b = GetConstProperty("MR_RRMS");
142 G4double b2 = b * b; << 152 G4double b2 = b*b;
143 w = GetConstProperty("MR_CORRLEN"); 153 w = GetConstProperty("MR_CORRLEN");
144 G4double w2 = w * w; << 154 G4double w2 = w*w;
145 155
146 no_theta_i = G4int(GetConstProperty("MR_NBTH << 156 no_theta_i = G4int(GetConstProperty("MR_NBTHETA")+0.1);
147 noE = G4int(GetConstProperty("MR_NBE") + 0.1 << 157 //G4cout << "no_theta: " << no_theta_i << G4endl;
>> 158 noE = G4int(GetConstProperty("MR_NBE")+0.1);
>> 159 //G4cout << "noE: " << noE << G4endl;
148 160
149 theta_i_min = GetConstProperty("MR_THETAMIN" 161 theta_i_min = GetConstProperty("MR_THETAMIN");
150 theta_i_max = GetConstProperty("MR_THETAMAX" 162 theta_i_max = GetConstProperty("MR_THETAMAX");
151 Emin = GetConstProperty("MR_EMIN"); 163 Emin = GetConstProperty("MR_EMIN");
152 Emax = GetConstProperty("MR_EMAX"); 164 Emax = GetConstProperty("MR_EMAX");
153 auto AngNoTheta = G4int(GetConstProperty("MR << 165 G4int AngNoTheta = G4int(GetConstProperty("MR_ANGNOTHETA")+0.1);
154 auto AngNoPhi = G4int(GetConstProperty("MR_A << 166 G4int AngNoPhi = G4int(GetConstProperty("MR_ANGNOPHI")+0.1);
155 AngCut = GetConstProperty("MR_ANGCUT"); 167 AngCut = GetConstProperty("MR_ANGCUT");
156 168
157 // The Fermi potential was saved in neV by P 169 // The Fermi potential was saved in neV by P.F.
158 170
159 G4double fermipot = GetConstProperty("FERMIP << 171 G4double fermipot = GetConstProperty("FERMIPOT")*(1.e-9*eV);
>> 172
>> 173 //G4cout << "Fermipot: " << fermipot/(1.e-9*eV) << "neV" << G4endl;
160 174
161 G4double theta_i, E; 175 G4double theta_i, E;
162 176
163 // Calculates the increment in theta_i in th 177 // Calculates the increment in theta_i in the lookup-table
164 theta_i_step = (theta_i_max - theta_i_min) / << 178 theta_i_step = (theta_i_max-theta_i_min)/(no_theta_i-1);
>> 179
>> 180 //G4cout << "theta_i_step: " << theta_i_step << G4endl;
165 181
166 // Calculates the increment in energy in the 182 // Calculates the increment in energy in the lookup-table
167 E_step = (Emax - Emin) / (noE - 1); << 183 E_step = (Emax-Emin)/(noE-1);
168 184
169 // Runs the lookup-table memory allocation 185 // Runs the lookup-table memory allocation
170 InitMicroRoughnessTables(); 186 InitMicroRoughnessTables();
171 187
172 G4int counter = 0; 188 G4int counter = 0;
173 189
>> 190 //G4cout << "Calculating Microroughnesstables...";
>> 191
174 // Writes the mr-lookup-tables to files for 192 // Writes the mr-lookup-tables to files for immediate control
175 std::ofstream dateir("MRrefl.dat", std::ios: <<
176 std::ofstream dateit("MRtrans.dat", std::ios <<
177 193
178 // G4cout << theMicroRoughnessTable << G4end << 194 std::ofstream dateir("MRrefl.dat",std::ios::out);
>> 195 std::ofstream dateit("MRtrans.dat",std::ios::out);
179 196
180 for (theta_i = theta_i_min; theta_i <= theta << 197 //G4cout << theMicroRoughnessTable << G4endl;
181 // Calculation for each cell in the lookup <<
182 for (E = Emin; E <= Emax; E += E_step) { <<
183 *(theMicroRoughnessTable + counter) = G4 <<
184 fermipot, theta_i, AngNoTheta, AngNoPh <<
185 198
186 *(theMicroRoughnessTransTable + counter) << 199 for (theta_i=theta_i_min; theta_i<=theta_i_max+1e-6; theta_i+=theta_i_step) {
187 G4UCNMicroRoughnessHelper::GetInstance << 200 // Calculation for each cell in the lookup-table
188 AngNoPhi, b2, w2, maxMicroRoughnessT << 201 for (E=Emin; E<=Emax; E+=E_step) {
>> 202 *(theMicroRoughnessTable+counter) =
>> 203 G4UCNMicroRoughnessHelper::GetInstance() ->
>> 204 IntIplus(E, fermipot, theta_i, AngNoTheta, AngNoPhi,
>> 205 b2, w2, maxMicroRoughnessTable+counter, AngCut);
>> 206
>> 207 *(theMicroRoughnessTransTable+counter) =
>> 208 G4UCNMicroRoughnessHelper::GetInstance() ->
>> 209 IntIminus(E, fermipot, theta_i, AngNoTheta, AngNoPhi,
>> 210 b2, w2, maxMicroRoughnessTransTable+counter,
>> 211 AngCut);
189 212
190 dateir << *(theMicroRoughnessTable + cou << 213 dateir << *(theMicroRoughnessTable+counter) << G4endl;
191 dateit << *(theMicroRoughnessTransTable << 214 dateit << *(theMicroRoughnessTransTable+counter) << G4endl;
192 215
193 counter++; << 216 counter++;
194 } << 217
>> 218 //G4cout << counter << G4endl;
>> 219 }
195 } 220 }
196 221
197 dateir.close(); 222 dateir.close();
198 dateit.close(); 223 dateit.close();
199 224
200 // Writes the mr-lookup-tables to files for 225 // Writes the mr-lookup-tables to files for immediate control
201 226
202 std::ofstream dateic("MRcheck.dat", std::ios << 227 std::ofstream dateic("MRcheck.dat",std::ios::out);
203 std::ofstream dateimr("MRmaxrefl.dat", std:: << 228 std::ofstream dateimr("MRmaxrefl.dat",std::ios::out);
204 std::ofstream dateimt("MRmaxtrans.dat", std: << 229 std::ofstream dateimt("MRmaxtrans.dat",std::ios::out);
205 << 230
206 for (theta_i = theta_i_min; theta_i <= theta << 231 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) { << 232 for (E=Emin; E<=Emax; E+=E_step) {
208 // tests the GetXXProbability functions << 233
209 // of the lookup tables to files << 234 // tests the GetXXProbability functions by writing the entries
210 << 235 // of the lookup tables to files
211 dateic << GetMRIntProbability(theta_i, E << 236
212 dateimr << GetMRMaxProbability(theta_i, << 237 dateic << GetMRIntProbability(theta_i,E) << G4endl;
213 dateimt << GetMRMaxTransProbability(thet << 238 dateimr << GetMRMaxProbability(theta_i,E) << G4endl;
214 } << 239 dateimt << GetMRMaxTransProbability(theta_i,E) << G4endl;
>> 240 }
215 } 241 }
216 242
217 dateic.close(); 243 dateic.close();
218 dateimr.close(); 244 dateimr.close();
219 dateimt.close(); 245 dateimt.close();
220 } 246 }
221 247
222 G4double G4UCNMaterialPropertiesTable::GetMRIn << 248 G4double G4UCNMaterialPropertiesTable::
>> 249 GetMRIntProbability(G4double theta_i, G4double Energy)
223 { 250 {
224 if (theMicroRoughnessTable == nullptr) { << 251 if(theMicroRoughnessTable == nullptr)
>> 252 {
225 G4cout << "Do not have theMicroRoughnessTa 253 G4cout << "Do not have theMicroRoughnessTable" << G4endl;
226 return 0.; 254 return 0.;
227 } 255 }
228 256
229 // if theta_i or energy outside the range fo 257 // if theta_i or energy outside the range for which the lookup table is
230 // calculated, the probability is set to zer 258 // calculated, the probability is set to zero
231 if (theta_i < theta_i_min || theta_i > theta << 259
>> 260 //G4cout << "theta_i: " << theta_i/degree << "degree"
>> 261 // << " theta_i_min: " << theta_i_min/degree << "degree"
>> 262 // << " theta_i_max: " << theta_i_max/degree << "degree"
>> 263 // << " Energy: " << Energy/(1.e-9*eV) << "neV"
>> 264 // << " Emin: " << Emin/(1.e-9*eV) << "neV"
>> 265 // << " Emax: " << Emax/(1.e-9*eV) << "neV" << G4endl;
>> 266
>> 267 if(theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin ||
>> 268 Energy > Emax)
>> 269 {
232 return 0.; 270 return 0.;
233 } 271 }
234 272
235 // Determines the nearest cell in the lookup 273 // Determines the nearest cell in the lookup table which contains
236 // the probability 274 // the probability
237 275
238 auto theta_i_pos = G4int((theta_i - theta_i_ << 276 G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5);
239 auto E_pos = G4int((Energy - Emin) / E_step << 277 G4int E_pos = G4int((Energy-Emin)/E_step+0.5);
240 278
241 // lookup table is onedimensional (1 row), e 279 // lookup table is onedimensional (1 row), energy is in rows,
242 // theta_i in columns 280 // theta_i in columns
243 return *(theMicroRoughnessTable + E_pos + th << 281
>> 282 //G4cout << "E_pos: " << E_pos << " theta_i_pos: " << theta_i_pos << G4endl;
>> 283 //G4cout << "Probability: " << *(theMicroRoughnessTable+E_pos+theta_i_pos*(noE-1)) << G4endl;
>> 284
>> 285 return *(theMicroRoughnessTable+E_pos+theta_i_pos*(noE - 1));
244 } 286 }
245 287
246 G4double G4UCNMaterialPropertiesTable::GetMRIn << 288 G4double G4UCNMaterialPropertiesTable::
>> 289 GetMRIntTransProbability(G4double theta_i, G4double Energy)
247 { 290 {
248 if (theMicroRoughnessTransTable == nullptr) << 291 if(theMicroRoughnessTransTable == nullptr)
>> 292 {
249 return 0.; 293 return 0.;
250 } 294 }
251 295
252 // if theta_i or energy outside the range fo 296 // if theta_i or energy outside the range for which the lookup table
253 // is calculated, the probability is set to 297 // is calculated, the probability is set to zero
254 298
255 if (theta_i < theta_i_min || theta_i > theta << 299 if(theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin ||
>> 300 Energy > Emax)
>> 301 {
256 return 0.; 302 return 0.;
257 } 303 }
258 304
259 // Determines the nearest cell in the lookup 305 // Determines the nearest cell in the lookup table which contains
260 // the probability 306 // the probability
261 307
262 auto theta_i_pos = G4int((theta_i - theta_i_ << 308 G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5);
263 auto E_pos = G4int((Energy - Emin) / E_step << 309 G4int E_pos = G4int((Energy-Emin)/E_step+0.5);
264 310
265 // lookup table is onedimensional (1 row), e 311 // lookup table is onedimensional (1 row), energy is in rows,
266 // theta_i in columns 312 // theta_i in columns
267 313
268 return *(theMicroRoughnessTransTable + E_pos << 314 return *(theMicroRoughnessTransTable+E_pos+theta_i_pos*(noE - 1));
269 } 315 }
270 316
271 G4double G4UCNMaterialPropertiesTable::GetMRMa << 317 G4double G4UCNMaterialPropertiesTable::
>> 318 GetMRMaxProbability(G4double theta_i, G4double Energy)
272 { 319 {
273 if (maxMicroRoughnessTable == nullptr) { << 320 if(maxMicroRoughnessTable == nullptr)
>> 321 {
274 return 0.; 322 return 0.;
275 } 323 }
276 324
277 // if theta_i or energy outside the range fo 325 // if theta_i or energy outside the range for which the lookup table
278 // is calculated, the probability is set to 326 // is calculated, the probability is set to zero
279 327
280 if (theta_i < theta_i_min || theta_i > theta << 328 if(theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin ||
>> 329 Energy > Emax)
>> 330 {
281 return 0.; 331 return 0.;
282 } 332 }
283 333
284 // Determines the nearest cell in the lookup 334 // Determines the nearest cell in the lookup table which contains
285 // the probability 335 // the probability
286 336
287 auto theta_i_pos = G4int((theta_i - theta_i_ << 337 G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5);
288 auto E_pos = G4int((Energy - Emin) / E_step << 338 G4int E_pos = G4int((Energy-Emin)/E_step+0.5);
289 339
290 // lookup table is onedimensional (1 row), e 340 // lookup table is onedimensional (1 row), energy is in rows,
291 // theta_i in columns 341 // theta_i in columns
292 342
293 return *(maxMicroRoughnessTable + E_pos + th << 343 return *(maxMicroRoughnessTable+E_pos+theta_i_pos*noE);
294 } 344 }
295 345
296 void G4UCNMaterialPropertiesTable::SetMRMaxPro << 346 void G4UCNMaterialPropertiesTable::
297 G4double theta_i, G4double Energy, G4double << 347 SetMRMaxProbability(G4double theta_i, G4double Energy, G4double value)
298 { 348 {
299 if (maxMicroRoughnessTable != nullptr) { << 349 if(maxMicroRoughnessTable != nullptr)
300 if (theta_i < theta_i_min || theta_i > the << 350 {
301 } << 351 if(theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin ||
302 else { << 352 Energy > Emax)
>> 353 {}
>> 354 else
>> 355 {
303 // Determines the nearest cell in the lo 356 // Determines the nearest cell in the lookup table which contains
304 // the probability 357 // the probability
305 358
306 auto theta_i_pos = G4int((theta_i - thet << 359 G4int theta_i_pos = G4int((theta_i - theta_i_min) / theta_i_step + 0.5);
307 auto E_pos = G4int((Energy - Emin) / E_s << 360 G4int E_pos = G4int((Energy - Emin) / E_step + 0.5);
308 361
309 // lookup table is onedimensional (1 row 362 // lookup table is onedimensional (1 row), energy is in rows,
310 // theta_i in columns 363 // theta_i in columns
311 364
312 *(maxMicroRoughnessTable + E_pos + theta 365 *(maxMicroRoughnessTable + E_pos + theta_i_pos * noE) = value;
313 } 366 }
314 } 367 }
315 } 368 }
316 369
317 G4double G4UCNMaterialPropertiesTable::GetMRMa << 370 G4double G4UCNMaterialPropertiesTable::
>> 371 GetMRMaxTransProbability(G4double theta_i, G4double Energy)
318 { 372 {
319 if (maxMicroRoughnessTransTable == nullptr) << 373 if(maxMicroRoughnessTransTable == nullptr)
>> 374 {
320 return 0.; 375 return 0.;
321 } 376 }
322 377
323 // if theta_i or energy outside the range fo 378 // if theta_i or energy outside the range for which the lookup table
324 // is calculated, the probability is set to 379 // is calculated, the probability is set to zero
325 380
326 if (theta_i < theta_i_min || theta_i > theta << 381 if(theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin ||
>> 382 Energy > Emax)
>> 383 {
327 return 0.; 384 return 0.;
328 } 385 }
329 386
330 // Determines the nearest cell in the lookup 387 // Determines the nearest cell in the lookup table which contains
331 // the probability 388 // the probability
332 389
333 auto theta_i_pos = G4int((theta_i - theta_i_ << 390 G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5);
334 auto E_pos = G4int((Energy - Emin) / E_step << 391 G4int E_pos = G4int((Energy-Emin)/E_step+0.5);
335 392
336 // lookup table is onedimensional (1 row), e 393 // lookup table is onedimensional (1 row), energy is in rows,
337 // theta_i in columns 394 // theta_i in columns
338 395
339 return *(maxMicroRoughnessTransTable + E_pos << 396 return *(maxMicroRoughnessTransTable+E_pos+theta_i_pos*noE);
340 } 397 }
341 398
342 void G4UCNMaterialPropertiesTable::SetMRMaxTra << 399 void G4UCNMaterialPropertiesTable::
343 G4double theta_i, G4double Energy, G4double << 400 SetMRMaxTransProbability(G4double theta_i, G4double Energy, G4double value)
344 { 401 {
345 if (maxMicroRoughnessTransTable != nullptr) << 402 if(maxMicroRoughnessTransTable != nullptr)
346 if (theta_i < theta_i_min || theta_i > the << 403 {
347 } << 404 if(theta_i < theta_i_min || theta_i > theta_i_max || Energy < Emin ||
348 else { << 405 Energy > Emax)
>> 406 {}
>> 407 else
>> 408 {
349 // Determines the nearest cell in the lo 409 // Determines the nearest cell in the lookup table which contains
350 // the probability 410 // the probability
351 411
352 auto theta_i_pos = G4int((theta_i - thet << 412 G4int theta_i_pos = G4int((theta_i - theta_i_min) / theta_i_step + 0.5);
353 auto E_pos = G4int((Energy - Emin) / E_s << 413 G4int E_pos = G4int((Energy - Emin) / E_step + 0.5);
354 414
355 // lookup table is onedimensional (1 row 415 // lookup table is onedimensional (1 row), energy is in rows,
356 // theta_i in columns 416 // theta_i in columns
357 417
358 *(maxMicroRoughnessTransTable + E_pos + 418 *(maxMicroRoughnessTransTable + E_pos + theta_i_pos * noE) = value;
359 } 419 }
360 } 420 }
361 } 421 }
362 422
363 G4double G4UCNMaterialPropertiesTable::GetMRPr << 423 G4double G4UCNMaterialPropertiesTable::
364 G4double theta_i, G4double Energy, G4double << 424 GetMRProbability(G4double theta_i, G4double Energy,
>> 425 G4double fermipot,
>> 426 G4double theta_o, G4double phi_o)
365 { 427 {
366 return G4UCNMicroRoughnessHelper::GetInstanc << 428 return G4UCNMicroRoughnessHelper::GetInstance()->
367 Energy, fermipot, theta_i, theta_o, phi_o, << 429 ProbIplus(Energy, fermipot, theta_i, theta_o, phi_o, b, w, AngCut);
368 } 430 }
369 431
370 G4double G4UCNMaterialPropertiesTable::GetMRTr << 432 G4double G4UCNMaterialPropertiesTable::
371 G4double theta_i, G4double Energy, G4double << 433 GetMRTransProbability(G4double theta_i, G4double Energy,
>> 434 G4double fermipot,
>> 435 G4double theta_o, G4double phi_o)
372 { 436 {
373 return G4UCNMicroRoughnessHelper::GetInstanc << 437 return G4UCNMicroRoughnessHelper::GetInstance()->
374 Energy, fermipot, theta_i, theta_o, phi_o, << 438 ProbIminus(Energy, fermipot,theta_i, theta_o, phi_o, b, w, AngCut);
375 } 439 }
376 440
377 G4bool G4UCNMaterialPropertiesTable::Condition << 441 G4bool G4UCNMaterialPropertiesTable::ConditionsValid(G4double E,
>> 442 G4double VFermi,
>> 443 G4double theta_i)
378 { 444 {
379 G4double k = std::sqrt(2 * neutron_mass_c2 * << 445 G4double k = std::sqrt(2*neutron_mass_c2*E / hbarc_squared);
380 G4double k_l = std::sqrt(2 * neutron_mass_c2 << 446 G4double k_l = std::sqrt(2*neutron_mass_c2*VFermi / hbarc_squared);
>> 447
>> 448 //G4cout << " Energy: " << E/(1.e-9*eV) << "neV"
>> 449 // << " VFermi: " << VFermi/(1.e-9*eV) << "neV"
>> 450 // << " theta: " << theta_i/degree << "degree" << G4endl;
>> 451
>> 452 //G4cout << " Conditions - 2*b*k*cos(theta): " << 2*b*k*cos(theta_i)
>> 453 // << ", 2*b*kl: " << 2*b*k_l << G4endl;
381 454
382 // see eq. 17 of the Steyerl paper 455 // see eq. 17 of the Steyerl paper
>> 456
383 return 2 * b * k * std::cos(theta_i) < 1 && 457 return 2 * b * k * std::cos(theta_i) < 1 && 2 * b * k_l < 1;
384 } 458 }
385 459
386 G4bool G4UCNMaterialPropertiesTable::TransCond << 460 G4bool G4UCNMaterialPropertiesTable::TransConditionsValid(G4double E,
387 G4double E, G4double VFermi, G4double theta_ << 461 G4double VFermi,
>> 462 G4double theta_i)
388 { 463 {
389 G4double k2 = 2 * neutron_mass_c2 * E / hbar << 464 G4double k2 = 2*neutron_mass_c2*E / hbarc_squared;
390 G4double k_l2 = 2 * neutron_mass_c2 * VFermi << 465 G4double k_l2 = 2*neutron_mass_c2*VFermi / hbarc_squared;
391 466
392 if (E * (std::cos(theta_i) * std::cos(theta_ << 467 if(E * (std::cos(theta_i) * std::cos(theta_i)) < VFermi)
>> 468 {
393 return false; 469 return false;
394 } 470 }
395 471
396 G4double kS2 = k_l2 - k2; 472 G4double kS2 = k_l2 - k2;
397 473
>> 474 //G4cout << "Conditions; 2bk' cos(theta): " << 2*b*sqrt(kS2)*cos(theta_i) <<
>> 475 // ", 2bk_l: " << 2*b*sqrt(k_l2) << G4endl;
>> 476
398 // see eq. 18 of the Steyerl paper 477 // see eq. 18 of the Steyerl paper
399 return 2 * b * std::sqrt(kS2) * std::cos(the << 478
>> 479 return 2 * b * std::sqrt(kS2) * std::cos(theta_i) < 1 &&
>> 480 2 * b * std::sqrt(k_l2) < 1;
400 } 481 }
401 482
402 void G4UCNMaterialPropertiesTable::SetMicroRou << 483 void G4UCNMaterialPropertiesTable::
403 G4int no_theta, G4int no_E, G4double theta_m << 484 SetMicroRoughnessParameters(G4double ww, G4double bb,
404 G4double E_max, G4int AngNoTheta, G4int AngN << 485 G4int no_theta, G4int no_E,
>> 486 G4double theta_min, G4double theta_max,
>> 487 G4double E_min, G4double E_max,
>> 488 G4int AngNoTheta, G4int AngNoPhi,
>> 489 G4double AngularCut)
405 { 490 {
>> 491 //G4cout << "Setting Microroughness Parameters...";
>> 492
406 // Removes an existing RMS roughness 493 // Removes an existing RMS roughness
407 if (ConstPropertyExists("MR_RRMS")) { << 494 if(ConstPropertyExists("MR_RRMS"))
>> 495 {
408 RemoveConstProperty("MR_RRMS"); 496 RemoveConstProperty("MR_RRMS");
409 } 497 }
410 498
411 // Adds a new RMS roughness 499 // Adds a new RMS roughness
412 AddConstProperty("MR_RRMS", bb); 500 AddConstProperty("MR_RRMS", bb);
413 501
>> 502 //G4cout << "b: " << bb << G4endl;
>> 503
414 // Removes an existing correlation length 504 // Removes an existing correlation length
415 if (ConstPropertyExists("MR_CORRLEN")) { << 505 if(ConstPropertyExists("MR_CORRLEN"))
>> 506 {
416 RemoveConstProperty("MR_CORRLEN"); 507 RemoveConstProperty("MR_CORRLEN");
417 } 508 }
418 509
419 // Adds a new correlation length 510 // Adds a new correlation length
420 AddConstProperty("MR_CORRLEN", ww); 511 AddConstProperty("MR_CORRLEN", ww);
421 512
>> 513 //G4cout << "w: " << ww << G4endl;
>> 514
422 // Removes an existing number of thetas 515 // Removes an existing number of thetas
423 if (ConstPropertyExists("MR_NBTHETA")) { << 516 if(ConstPropertyExists("MR_NBTHETA"))
>> 517 {
424 RemoveConstProperty("MR_NBTHETA"); 518 RemoveConstProperty("MR_NBTHETA");
425 } 519 }
426 520
427 // Adds a new number of thetas 521 // Adds a new number of thetas
428 AddConstProperty("MR_NBTHETA", (G4double)no_ 522 AddConstProperty("MR_NBTHETA", (G4double)no_theta);
429 523
>> 524 //G4cout << "no_theta: " << no_theta << G4endl;
>> 525
430 // Removes an existing number of Energies 526 // Removes an existing number of Energies
431 if (ConstPropertyExists("MR_NBE")) { << 527 if(ConstPropertyExists("MR_NBE"))
>> 528 {
432 RemoveConstProperty("MR_NBE"); 529 RemoveConstProperty("MR_NBE");
433 } 530 }
434 531
435 // Adds a new number of Energies 532 // Adds a new number of Energies
436 AddConstProperty("MR_NBE", (G4double)no_E); 533 AddConstProperty("MR_NBE", (G4double)no_E);
437 534
>> 535 //G4cout << "no_E: " << no_E << G4endl;
>> 536
438 // Removes an existing minimum theta 537 // Removes an existing minimum theta
439 if (ConstPropertyExists("MR_THETAMIN")) { << 538 if(ConstPropertyExists("MR_THETAMIN"))
>> 539 {
440 RemoveConstProperty("MR_THETAMIN"); 540 RemoveConstProperty("MR_THETAMIN");
441 } 541 }
442 542
443 // Adds a new minimum theta 543 // Adds a new minimum theta
444 AddConstProperty("MR_THETAMIN", theta_min); 544 AddConstProperty("MR_THETAMIN", theta_min);
445 545
>> 546 //G4cout << "theta_min: " << theta_min << G4endl;
>> 547
446 // Removes an existing maximum theta 548 // Removes an existing maximum theta
447 if (ConstPropertyExists("MR_THETAMAX")) { << 549 if(ConstPropertyExists("MR_THETAMAX"))
>> 550 {
448 RemoveConstProperty("MR_THETAMAX"); 551 RemoveConstProperty("MR_THETAMAX");
449 } 552 }
450 553
451 // Adds a new maximum theta 554 // Adds a new maximum theta
452 AddConstProperty("MR_THETAMAX", theta_max); 555 AddConstProperty("MR_THETAMAX", theta_max);
453 556
>> 557 //G4cout << "theta_max: " << theta_max << G4endl;
>> 558
454 // Removes an existing minimum energy 559 // Removes an existing minimum energy
455 if (ConstPropertyExists("MR_EMIN")) { << 560 if(ConstPropertyExists("MR_EMIN"))
>> 561 {
456 RemoveConstProperty("MR_EMIN"); 562 RemoveConstProperty("MR_EMIN");
457 } 563 }
458 564
459 // Adds a new minimum energy 565 // Adds a new minimum energy
460 AddConstProperty("MR_EMIN", E_min); 566 AddConstProperty("MR_EMIN", E_min);
461 567
>> 568 //G4cout << "Emin: " << E_min << G4endl;
>> 569
462 // Removes an existing maximum energy 570 // Removes an existing maximum energy
463 if (ConstPropertyExists("MR_EMAX")) { << 571 if(ConstPropertyExists("MR_EMAX"))
>> 572 {
464 RemoveConstProperty("MR_EMAX"); 573 RemoveConstProperty("MR_EMAX");
465 } 574 }
466 575
467 // Adds a new maximum energy 576 // Adds a new maximum energy
468 AddConstProperty("MR_EMAX", E_max); 577 AddConstProperty("MR_EMAX", E_max);
469 578
>> 579 //G4cout << "Emax: " << E_max << G4endl;
>> 580
470 // Removes an existing Theta angle number 581 // Removes an existing Theta angle number
471 if (ConstPropertyExists("MR_ANGNOTHETA")) { << 582 if(ConstPropertyExists("MR_ANGNOTHETA"))
>> 583 {
472 RemoveConstProperty("MR_ANGNOTHETA"); 584 RemoveConstProperty("MR_ANGNOTHETA");
473 } 585 }
474 586
475 // Adds a new Theta angle number 587 // Adds a new Theta angle number
476 AddConstProperty("MR_ANGNOTHETA", (G4double) 588 AddConstProperty("MR_ANGNOTHETA", (G4double)AngNoTheta);
477 589
>> 590 //G4cout << "AngNoTheta: " << AngNoTheta << G4endl;
>> 591
478 // Removes an existing Phi angle number 592 // Removes an existing Phi angle number
479 if (ConstPropertyExists("MR_ANGNOPHI")) { << 593 if(ConstPropertyExists("MR_ANGNOPHI"))
>> 594 {
480 RemoveConstProperty("MR_ANGNOPHI"); 595 RemoveConstProperty("MR_ANGNOPHI");
481 } 596 }
482 597
483 // Adds a new Phi angle number 598 // Adds a new Phi angle number
484 AddConstProperty("MR_ANGNOPHI", (G4double)An 599 AddConstProperty("MR_ANGNOPHI", (G4double)AngNoPhi);
485 600
>> 601 //G4cout << "AngNoPhi: " << AngNoPhi << G4endl;
>> 602
486 // Removes an existing angular cut 603 // Removes an existing angular cut
487 if (ConstPropertyExists("MR_ANGCUT")) { << 604 if(ConstPropertyExists("MR_ANGCUT"))
>> 605 {
488 RemoveConstProperty("MR_ANGCUT"); 606 RemoveConstProperty("MR_ANGCUT");
489 } 607 }
490 608
491 // Adds a new angle number 609 // Adds a new angle number
492 AddConstProperty("MR_ANGCUT", AngularCut); 610 AddConstProperty("MR_ANGCUT", AngularCut);
493 611
>> 612 //G4cout << "AngularCut: " << AngularCut/degree << "degree" << G4endl;
>> 613
494 // Starts the lookup table calculation 614 // Starts the lookup table calculation
>> 615
495 ComputeMicroRoughnessTables(); 616 ComputeMicroRoughnessTables();
>> 617
>> 618 //G4cout << " *** DONE! ***" << G4endl;
496 } 619 }
497 620