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 if (theMicroRoughnessTable) delete theMicroRoughnessTable; 81 delete maxMicroRoughnessTable; << 82 if (maxMicroRoughnessTable) delete maxMicroRoughnessTable; 82 delete theMicroRoughnessTransTable; << 83 if (theMicroRoughnessTransTable) delete theMicroRoughnessTransTable; 83 delete maxMicroRoughnessTransTable; << 84 if (maxMicroRoughnessTransTable) 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")) 111 Nthetadim = G4int(GetConstProperty("MR_NBT << 117 Nthetadim = G4int(GetConstProperty("MR_NBTHETA")+0.1); 112 } << 113 118 114 // Checks if the number of energies is avail 119 // Checks if the number of energies is available and stores it 115 120 116 if (ConstPropertyExists("MR_NBE")) { << 121 if (ConstPropertyExists("MR_NBE")) 117 NEdim = G4int(GetConstProperty("MR_NBE") + << 122 NEdim = G4int(GetConstProperty("MR_NBE")+0.1); 118 } << 123 >> 124 //G4cout << "thetadim: " << Nthetadim << " , Edim: " << NEdim << G4endl; 119 125 120 // If both dimensions of the lookup-table ar 126 // If both dimensions of the lookup-table are non-trivial: 121 // delete old tables if existing and allocat 127 // delete old tables if existing and allocate memory for new tables 122 128 123 if (Nthetadim * NEdim > 0) { << 129 if (Nthetadim*NEdim > 0) { 124 delete theMicroRoughnessTable; << 130 if (theMicroRoughnessTable) delete theMicroRoughnessTable; 125 theMicroRoughnessTable = new G4double[Nthe << 131 theMicroRoughnessTable = new G4double[Nthetadim*NEdim]; 126 delete maxMicroRoughnessTable; << 132 if (maxMicroRoughnessTable) delete maxMicroRoughnessTable; 127 maxMicroRoughnessTable = new G4double[Nthe << 133 maxMicroRoughnessTable = new G4double[Nthetadim*NEdim]; 128 delete theMicroRoughnessTransTable; << 134 if (theMicroRoughnessTransTable) delete theMicroRoughnessTransTable; 129 theMicroRoughnessTransTable = new G4double << 135 theMicroRoughnessTransTable = new G4double[Nthetadim*NEdim]; 130 delete maxMicroRoughnessTransTable; << 136 if (maxMicroRoughnessTransTable) delete maxMicroRoughnessTransTable; 131 maxMicroRoughnessTransTable = new G4double << 137 maxMicroRoughnessTransTable = new G4double[Nthetadim*NEdim]; 132 } 138 } 133 } 139 } 134 140 135 void G4UCNMaterialPropertiesTable::ComputeMicr 141 void G4UCNMaterialPropertiesTable::ComputeMicroRoughnessTables() 136 { 142 { 137 // Reads the parameters for the mr-probabili << 143 // Reads the parameters for the mr-probability computation from the 138 // corresponding material properties and sto << 144 // corresponding material properties and stores it in the appropriate 139 // variables << 145 // variables 140 146 141 b = GetConstProperty("MR_RRMS"); 147 b = GetConstProperty("MR_RRMS"); 142 G4double b2 = b * b; << 148 G4double b2 = b*b; 143 w = GetConstProperty("MR_CORRLEN"); 149 w = GetConstProperty("MR_CORRLEN"); 144 G4double w2 = w * w; << 150 G4double w2 = w*w; 145 151 146 no_theta_i = G4int(GetConstProperty("MR_NBTH << 152 no_theta_i = G4int(GetConstProperty("MR_NBTHETA")+0.1); 147 noE = G4int(GetConstProperty("MR_NBE") + 0.1 << 153 //G4cout << "no_theta: " << no_theta_i << G4endl; >> 154 noE = G4int(GetConstProperty("MR_NBE")+0.1); >> 155 //G4cout << "noE: " << noE << G4endl; 148 156 149 theta_i_min = GetConstProperty("MR_THETAMIN" 157 theta_i_min = GetConstProperty("MR_THETAMIN"); 150 theta_i_max = GetConstProperty("MR_THETAMAX" 158 theta_i_max = GetConstProperty("MR_THETAMAX"); 151 Emin = GetConstProperty("MR_EMIN"); 159 Emin = GetConstProperty("MR_EMIN"); 152 Emax = GetConstProperty("MR_EMAX"); 160 Emax = GetConstProperty("MR_EMAX"); 153 auto AngNoTheta = G4int(GetConstProperty("MR << 161 G4int AngNoTheta = G4int(GetConstProperty("MR_ANGNOTHETA")+0.1); 154 auto AngNoPhi = G4int(GetConstProperty("MR_A << 162 G4int AngNoPhi = G4int(GetConstProperty("MR_ANGNOPHI")+0.1); 155 AngCut = GetConstProperty("MR_ANGCUT"); 163 AngCut = GetConstProperty("MR_ANGCUT"); 156 164 157 // The Fermi potential was saved in neV by P 165 // The Fermi potential was saved in neV by P.F. 158 166 159 G4double fermipot = GetConstProperty("FERMIP << 167 G4double fermipot = GetConstProperty("FERMIPOT")*(1.e-9*eV); >> 168 >> 169 //G4cout << "Fermipot: " << fermipot/(1.e-9*eV) << "neV" << G4endl; 160 170 161 G4double theta_i, E; 171 G4double theta_i, E; 162 172 163 // Calculates the increment in theta_i in th 173 // Calculates the increment in theta_i in the lookup-table 164 theta_i_step = (theta_i_max - theta_i_min) / << 174 theta_i_step = (theta_i_max-theta_i_min)/(no_theta_i-1); >> 175 >> 176 //G4cout << "theta_i_step: " << theta_i_step << G4endl; 165 177 166 // Calculates the increment in energy in the 178 // Calculates the increment in energy in the lookup-table 167 E_step = (Emax - Emin) / (noE - 1); << 179 E_step = (Emax-Emin)/(noE-1); 168 180 169 // Runs the lookup-table memory allocation 181 // Runs the lookup-table memory allocation 170 InitMicroRoughnessTables(); 182 InitMicroRoughnessTables(); 171 183 172 G4int counter = 0; 184 G4int counter = 0; 173 185 >> 186 //G4cout << "Calculating Microroughnesstables..."; >> 187 174 // Writes the mr-lookup-tables to files for 188 // 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 189 178 // G4cout << theMicroRoughnessTable << G4end << 190 std::ofstream dateir("MRrefl.dat",std::ios::out); >> 191 std::ofstream dateit("MRtrans.dat",std::ios::out); >> 192 >> 193 //G4cout << theMicroRoughnessTable << G4endl; 179 194 180 for (theta_i = theta_i_min; theta_i <= theta << 195 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 << 196 // Calculation for each cell in the lookup-table 182 for (E = Emin; E <= Emax; E += E_step) { << 197 for (E=Emin; E<=Emax; E+=E_step) { 183 *(theMicroRoughnessTable + counter) = G4 << 198 *(theMicroRoughnessTable+counter) = 184 fermipot, theta_i, AngNoTheta, AngNoPh << 199 G4UCNMicroRoughnessHelper::GetInstance() -> 185 << 200 IntIplus(E, fermipot, theta_i, AngNoTheta, AngNoPhi, 186 *(theMicroRoughnessTransTable + counter) << 201 b2, w2, maxMicroRoughnessTable+counter, AngCut); 187 G4UCNMicroRoughnessHelper::GetInstance << 202 188 AngNoPhi, b2, w2, maxMicroRoughnessT << 203 *(theMicroRoughnessTransTable+counter) = >> 204 G4UCNMicroRoughnessHelper::GetInstance() -> >> 205 IntIminus(E, fermipot, theta_i, AngNoTheta, AngNoPhi, >> 206 b2, w2, maxMicroRoughnessTransTable+counter, >> 207 AngCut); 189 208 190 dateir << *(theMicroRoughnessTable + cou << 209 dateir << *(theMicroRoughnessTable+counter) << G4endl; 191 dateit << *(theMicroRoughnessTransTable << 210 dateit << *(theMicroRoughnessTransTable+counter) << G4endl; 192 211 193 counter++; << 212 counter++; 194 } << 213 >> 214 //G4cout << counter << G4endl; >> 215 } 195 } 216 } 196 217 197 dateir.close(); 218 dateir.close(); 198 dateit.close(); 219 dateit.close(); 199 220 200 // Writes the mr-lookup-tables to files for 221 // Writes the mr-lookup-tables to files for immediate control 201 222 202 std::ofstream dateic("MRcheck.dat", std::ios << 223 std::ofstream dateic("MRcheck.dat",std::ios::out); 203 std::ofstream dateimr("MRmaxrefl.dat", std:: << 224 std::ofstream dateimr("MRmaxrefl.dat",std::ios::out); 204 std::ofstream dateimt("MRmaxtrans.dat", std: << 225 std::ofstream dateimt("MRmaxtrans.dat",std::ios::out); 205 << 226 206 for (theta_i = theta_i_min; theta_i <= theta << 227 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) { << 228 for (E=Emin; E<=Emax; E+=E_step) { 208 // tests the GetXXProbability functions << 229 209 // of the lookup tables to files << 230 // tests the GetXXProbability functions by writing the entries 210 << 231 // of the lookup tables to files 211 dateic << GetMRIntProbability(theta_i, E << 232 212 dateimr << GetMRMaxProbability(theta_i, << 233 dateic << GetMRIntProbability(theta_i,E) << G4endl; 213 dateimt << GetMRMaxTransProbability(thet << 234 dateimr << GetMRMaxProbability(theta_i,E) << G4endl; 214 } << 235 dateimt << GetMRMaxTransProbability(theta_i,E) << G4endl; >> 236 } 215 } 237 } 216 238 217 dateic.close(); 239 dateic.close(); 218 dateimr.close(); 240 dateimr.close(); 219 dateimt.close(); 241 dateimt.close(); 220 } 242 } 221 243 222 G4double G4UCNMaterialPropertiesTable::GetMRIn << 244 G4double G4UCNMaterialPropertiesTable:: >> 245 GetMRIntProbability(G4double theta_i, G4double Energy) 223 { 246 { 224 if (theMicroRoughnessTable == nullptr) { << 247 if (!theMicroRoughnessTable) { 225 G4cout << "Do not have theMicroRoughnessTa << 248 G4cout << "Dont have theMicroRoughnessTable" << G4endl; 226 return 0.; << 249 return 0.; 227 } 250 } 228 251 229 // if theta_i or energy outside the range fo 252 // if theta_i or energy outside the range for which the lookup table is 230 // calculated, the probability is set to zer 253 // calculated, the probability is set to zero 231 if (theta_i < theta_i_min || theta_i > theta << 254 232 return 0.; << 255 //G4cout << "theta_i: " << theta_i/degree << "degree" 233 } << 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.; 234 264 235 // Determines the nearest cell in the lookup 265 // Determines the nearest cell in the lookup table which contains 236 // the probability 266 // the probability 237 267 238 auto theta_i_pos = G4int((theta_i - theta_i_ << 268 G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5); 239 auto E_pos = G4int((Energy - Emin) / E_step << 269 G4int E_pos = G4int((Energy-Emin)/E_step+0.5); 240 270 241 // lookup table is onedimensional (1 row), e 271 // lookup table is onedimensional (1 row), energy is in rows, 242 // theta_i in columns 272 // theta_i in columns 243 return *(theMicroRoughnessTable + E_pos + th << 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)); 244 } 278 } 245 279 246 G4double G4UCNMaterialPropertiesTable::GetMRIn << 280 G4double G4UCNMaterialPropertiesTable:: >> 281 GetMRIntTransProbability(G4double theta_i, G4double Energy) 247 { 282 { 248 if (theMicroRoughnessTransTable == nullptr) << 283 if (!theMicroRoughnessTransTable) return 0.; 249 return 0.; << 250 } << 251 284 252 // if theta_i or energy outside the range fo 285 // if theta_i or energy outside the range for which the lookup table 253 // is calculated, the probability is set to 286 // is calculated, the probability is set to zero 254 287 255 if (theta_i < theta_i_min || theta_i > theta << 288 if (theta_i<theta_i_min || theta_i>theta_i_max || Energy<Emin || Energy>Emax) 256 return 0.; << 289 return 0.; 257 } << 258 290 259 // Determines the nearest cell in the lookup 291 // Determines the nearest cell in the lookup table which contains 260 // the probability 292 // the probability 261 293 262 auto theta_i_pos = G4int((theta_i - theta_i_ << 294 G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5); 263 auto E_pos = G4int((Energy - Emin) / E_step << 295 G4int E_pos = G4int((Energy-Emin)/E_step+0.5); 264 296 265 // lookup table is onedimensional (1 row), e 297 // lookup table is onedimensional (1 row), energy is in rows, 266 // theta_i in columns 298 // theta_i in columns 267 299 268 return *(theMicroRoughnessTransTable + E_pos << 300 return *(theMicroRoughnessTransTable+E_pos+theta_i_pos*(noE - 1)); 269 } 301 } 270 302 271 G4double G4UCNMaterialPropertiesTable::GetMRMa << 303 G4double G4UCNMaterialPropertiesTable:: >> 304 GetMRMaxProbability(G4double theta_i, G4double Energy) 272 { 305 { 273 if (maxMicroRoughnessTable == nullptr) { << 306 if (!maxMicroRoughnessTable) return 0.; 274 return 0.; << 275 } << 276 307 277 // if theta_i or energy outside the range fo 308 // if theta_i or energy outside the range for which the lookup table 278 // is calculated, the probability is set to 309 // is calculated, the probability is set to zero 279 310 280 if (theta_i < theta_i_min || theta_i > theta << 311 if (theta_i<theta_i_min || theta_i>theta_i_max || Energy<Emin || Energy>Emax) 281 return 0.; << 312 return 0.; 282 } << 283 313 284 // Determines the nearest cell in the lookup 314 // Determines the nearest cell in the lookup table which contains 285 // the probability 315 // the probability 286 316 287 auto theta_i_pos = G4int((theta_i - theta_i_ << 317 G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5); 288 auto E_pos = G4int((Energy - Emin) / E_step << 318 G4int E_pos = G4int((Energy-Emin)/E_step+0.5); 289 319 290 // lookup table is onedimensional (1 row), e 320 // lookup table is onedimensional (1 row), energy is in rows, 291 // theta_i in columns 321 // theta_i in columns 292 322 293 return *(maxMicroRoughnessTable + E_pos + th << 323 return *(maxMicroRoughnessTable+E_pos+theta_i_pos*noE); 294 } 324 } 295 325 296 void G4UCNMaterialPropertiesTable::SetMRMaxPro << 326 void G4UCNMaterialPropertiesTable:: 297 G4double theta_i, G4double Energy, G4double << 327 SetMRMaxProbability(G4double theta_i, G4double Energy, G4double value) 298 { 328 { 299 if (maxMicroRoughnessTable != nullptr) { << 329 if (maxMicroRoughnessTable) { 300 if (theta_i < theta_i_min || theta_i > the << 301 } << 302 else { << 303 // Determines the nearest cell in the lo << 304 // the probability << 305 330 306 auto theta_i_pos = G4int((theta_i - thet << 331 if (theta_i<theta_i_min || theta_i>theta_i_max || 307 auto E_pos = G4int((Energy - Emin) / E_s << 332 Energy<Emin || Energy>Emax) { >> 333 } else { 308 334 309 // lookup table is onedimensional (1 row << 335 // Determines the nearest cell in the lookup table which contains 310 // theta_i in columns << 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); 311 340 312 *(maxMicroRoughnessTable + E_pos + theta << 341 // lookup table is onedimensional (1 row), energy is in rows, 313 } << 342 // theta_i in columns >> 343 >> 344 *(maxMicroRoughnessTable+E_pos+theta_i_pos*noE) = value; >> 345 } 314 } 346 } 315 } 347 } 316 348 317 G4double G4UCNMaterialPropertiesTable::GetMRMa << 349 G4double G4UCNMaterialPropertiesTable:: >> 350 GetMRMaxTransProbability(G4double theta_i, G4double Energy) 318 { 351 { 319 if (maxMicroRoughnessTransTable == nullptr) << 352 if (!maxMicroRoughnessTransTable) return 0.; 320 return 0.; << 321 } << 322 353 323 // if theta_i or energy outside the range fo 354 // if theta_i or energy outside the range for which the lookup table 324 // is calculated, the probability is set to 355 // is calculated, the probability is set to zero 325 356 326 if (theta_i < theta_i_min || theta_i > theta << 357 if (theta_i<theta_i_min || theta_i>theta_i_max || Energy<Emin || Energy>Emax) 327 return 0.; << 358 return 0.; 328 } << 329 359 330 // Determines the nearest cell in the lookup 360 // Determines the nearest cell in the lookup table which contains 331 // the probability 361 // the probability 332 362 333 auto theta_i_pos = G4int((theta_i - theta_i_ << 363 G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5); 334 auto E_pos = G4int((Energy - Emin) / E_step << 364 G4int E_pos = G4int((Energy-Emin)/E_step+0.5); 335 365 336 // lookup table is onedimensional (1 row), e 366 // lookup table is onedimensional (1 row), energy is in rows, 337 // theta_i in columns 367 // theta_i in columns 338 368 339 return *(maxMicroRoughnessTransTable + E_pos << 369 return *(maxMicroRoughnessTransTable+E_pos+theta_i_pos*noE); 340 } 370 } 341 371 342 void G4UCNMaterialPropertiesTable::SetMRMaxTra << 372 void G4UCNMaterialPropertiesTable:: 343 G4double theta_i, G4double Energy, G4double << 373 SetMRMaxTransProbability(G4double theta_i, G4double Energy, G4double value) 344 { 374 { 345 if (maxMicroRoughnessTransTable != nullptr) << 375 if (maxMicroRoughnessTransTable) { 346 if (theta_i < theta_i_min || theta_i > the << 376 347 } << 377 if (theta_i<theta_i_min || theta_i>theta_i_max || 348 else { << 378 Energy<Emin || Energy>Emax) { 349 // Determines the nearest cell in the lo << 379 } else { 350 // the probability << 351 380 352 auto theta_i_pos = G4int((theta_i - thet << 381 // Determines the nearest cell in the lookup table which contains 353 auto E_pos = G4int((Energy - Emin) / E_s << 382 // the probability 354 383 355 // lookup table is onedimensional (1 row << 384 G4int theta_i_pos = G4int((theta_i-theta_i_min)/theta_i_step+0.5); 356 // theta_i in columns << 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 357 389 358 *(maxMicroRoughnessTransTable + E_pos + << 390 *(maxMicroRoughnessTransTable+E_pos+theta_i_pos*noE) = value; 359 } << 391 } 360 } 392 } 361 } 393 } 362 394 363 G4double G4UCNMaterialPropertiesTable::GetMRPr << 395 G4double G4UCNMaterialPropertiesTable:: 364 G4double theta_i, G4double Energy, G4double << 396 GetMRProbability(G4double theta_i, G4double Energy, >> 397 G4double fermipot, >> 398 G4double theta_o, G4double phi_o) 365 { 399 { 366 return G4UCNMicroRoughnessHelper::GetInstanc << 400 return G4UCNMicroRoughnessHelper::GetInstance()-> 367 Energy, fermipot, theta_i, theta_o, phi_o, << 401 ProbIplus(Energy, fermipot, theta_i, theta_o, phi_o, b, w, AngCut); 368 } 402 } 369 403 370 G4double G4UCNMaterialPropertiesTable::GetMRTr << 404 G4double G4UCNMaterialPropertiesTable:: 371 G4double theta_i, G4double Energy, G4double << 405 GetMRTransProbability(G4double theta_i, G4double Energy, >> 406 G4double fermipot, >> 407 G4double theta_o, G4double phi_o) 372 { 408 { 373 return G4UCNMicroRoughnessHelper::GetInstanc << 409 return G4UCNMicroRoughnessHelper::GetInstance()-> 374 Energy, fermipot, theta_i, theta_o, phi_o, << 410 ProbIminus(Energy, fermipot,theta_i, theta_o, phi_o, b, w, AngCut); 375 } 411 } 376 412 377 G4bool G4UCNMaterialPropertiesTable::Condition << 413 G4bool G4UCNMaterialPropertiesTable::ConditionsValid(G4double E, >> 414 G4double VFermi, >> 415 G4double theta_i) 378 { 416 { 379 G4double k = std::sqrt(2 * neutron_mass_c2 * << 417 G4double k = std::sqrt(2*neutron_mass_c2*E / hbarc_squared); 380 G4double k_l = std::sqrt(2 * neutron_mass_c2 << 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; 381 426 382 // see eq. 17 of the Steyerl paper 427 // see eq. 17 of the Steyerl paper 383 return 2 * b * k * std::cos(theta_i) < 1 && << 428 >> 429 if (2*b*k*std::cos(theta_i) < 1 && 2*b*k_l < 1) return true; >> 430 else return false; 384 } 431 } 385 432 386 G4bool G4UCNMaterialPropertiesTable::TransCond << 433 G4bool G4UCNMaterialPropertiesTable::TransConditionsValid(G4double E, 387 G4double E, G4double VFermi, G4double theta_ << 434 G4double VFermi, >> 435 G4double theta_i) 388 { 436 { 389 G4double k2 = 2 * neutron_mass_c2 * E / hbar << 437 G4double k2 = 2*neutron_mass_c2*E / hbarc_squared; 390 G4double k_l2 = 2 * neutron_mass_c2 * VFermi << 438 G4double k_l2 = 2*neutron_mass_c2*VFermi / hbarc_squared; 391 439 392 if (E * (std::cos(theta_i) * std::cos(theta_ << 440 if (E*(std::cos(theta_i)*std::cos(theta_i)) < VFermi) return false; 393 return false; << 394 } << 395 441 396 G4double kS2 = k_l2 - k2; 442 G4double kS2 = k_l2 - k2; 397 443 >> 444 //G4cout << "Conditions; 2bk' cos(theta): " << 2*b*sqrt(kS2)*cos(theta_i) << >> 445 // ", 2bk_l: " << 2*b*sqrt(k_l2) << G4endl; >> 446 398 // see eq. 18 of the Steyerl paper 447 // see eq. 18 of the Steyerl paper 399 return 2 * b * std::sqrt(kS2) * std::cos(the << 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; 400 } 451 } 401 452 402 void G4UCNMaterialPropertiesTable::SetMicroRou << 453 void G4UCNMaterialPropertiesTable:: 403 G4int no_theta, G4int no_E, G4double theta_m << 454 SetMicroRoughnessParameters(G4double ww, G4double bb, 404 G4double E_max, G4int AngNoTheta, G4int AngN << 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) 405 { 460 { >> 461 //G4cout << "Setting Microroughness Parameters..."; >> 462 406 // Removes an existing RMS roughness 463 // Removes an existing RMS roughness 407 if (ConstPropertyExists("MR_RRMS")) { << 464 if (ConstPropertyExists("MR_RRMS")) RemoveConstProperty("MR_RRMS"); 408 RemoveConstProperty("MR_RRMS"); << 409 } << 410 465 411 // Adds a new RMS roughness 466 // Adds a new RMS roughness 412 AddConstProperty("MR_RRMS", bb); 467 AddConstProperty("MR_RRMS", bb); 413 468 >> 469 //G4cout << "b: " << bb << G4endl; >> 470 414 // Removes an existing correlation length 471 // Removes an existing correlation length 415 if (ConstPropertyExists("MR_CORRLEN")) { << 472 if (ConstPropertyExists("MR_CORRLEN")) RemoveConstProperty("MR_CORRLEN"); 416 RemoveConstProperty("MR_CORRLEN"); << 417 } << 418 473 419 // Adds a new correlation length 474 // Adds a new correlation length 420 AddConstProperty("MR_CORRLEN", ww); 475 AddConstProperty("MR_CORRLEN", ww); 421 476 >> 477 //G4cout << "w: " << ww << G4endl; >> 478 422 // Removes an existing number of thetas 479 // Removes an existing number of thetas 423 if (ConstPropertyExists("MR_NBTHETA")) { << 480 if (ConstPropertyExists("MR_NBTHETA")) RemoveConstProperty("MR_NBTHETA"); 424 RemoveConstProperty("MR_NBTHETA"); << 425 } << 426 481 427 // Adds a new number of thetas 482 // Adds a new number of thetas 428 AddConstProperty("MR_NBTHETA", (G4double)no_ 483 AddConstProperty("MR_NBTHETA", (G4double)no_theta); 429 484 >> 485 //G4cout << "no_theta: " << no_theta << G4endl; >> 486 430 // Removes an existing number of Energies 487 // Removes an existing number of Energies 431 if (ConstPropertyExists("MR_NBE")) { << 488 if (ConstPropertyExists("MR_NBE")) RemoveConstProperty("MR_NBE"); 432 RemoveConstProperty("MR_NBE"); << 433 } << 434 489 435 // Adds a new number of Energies 490 // Adds a new number of Energies 436 AddConstProperty("MR_NBE", (G4double)no_E); 491 AddConstProperty("MR_NBE", (G4double)no_E); 437 492 >> 493 //G4cout << "no_E: " << no_E << G4endl; >> 494 438 // Removes an existing minimum theta 495 // Removes an existing minimum theta 439 if (ConstPropertyExists("MR_THETAMIN")) { << 496 if (ConstPropertyExists("MR_THETAMIN")) RemoveConstProperty("MR_THETAMIN"); 440 RemoveConstProperty("MR_THETAMIN"); << 441 } << 442 497 443 // Adds a new minimum theta 498 // Adds a new minimum theta 444 AddConstProperty("MR_THETAMIN", theta_min); 499 AddConstProperty("MR_THETAMIN", theta_min); 445 500 >> 501 //G4cout << "theta_min: " << theta_min << G4endl; >> 502 446 // Removes an existing maximum theta 503 // Removes an existing maximum theta 447 if (ConstPropertyExists("MR_THETAMAX")) { << 504 if (ConstPropertyExists("MR_THETAMAX")) RemoveConstProperty("MR_THETAMAX"); 448 RemoveConstProperty("MR_THETAMAX"); << 449 } << 450 505 451 // Adds a new maximum theta 506 // Adds a new maximum theta 452 AddConstProperty("MR_THETAMAX", theta_max); 507 AddConstProperty("MR_THETAMAX", theta_max); 453 508 >> 509 //G4cout << "theta_max: " << theta_max << G4endl; >> 510 454 // Removes an existing minimum energy 511 // Removes an existing minimum energy 455 if (ConstPropertyExists("MR_EMIN")) { << 512 if (ConstPropertyExists("MR_EMIN")) RemoveConstProperty("MR_EMIN"); 456 RemoveConstProperty("MR_EMIN"); << 457 } << 458 513 459 // Adds a new minimum energy 514 // Adds a new minimum energy 460 AddConstProperty("MR_EMIN", E_min); 515 AddConstProperty("MR_EMIN", E_min); 461 516 >> 517 //G4cout << "Emin: " << E_min << G4endl; >> 518 462 // Removes an existing maximum energy 519 // Removes an existing maximum energy 463 if (ConstPropertyExists("MR_EMAX")) { << 520 if (ConstPropertyExists("MR_EMAX")) RemoveConstProperty("MR_EMAX"); 464 RemoveConstProperty("MR_EMAX"); << 465 } << 466 521 467 // Adds a new maximum energy 522 // Adds a new maximum energy 468 AddConstProperty("MR_EMAX", E_max); 523 AddConstProperty("MR_EMAX", E_max); 469 524 >> 525 //G4cout << "Emax: " << E_max << G4endl; >> 526 470 // Removes an existing Theta angle number 527 // Removes an existing Theta angle number 471 if (ConstPropertyExists("MR_ANGNOTHETA")) { << 528 if(ConstPropertyExists("MR_ANGNOTHETA"))RemoveConstProperty("MR_ANGNOTHETA"); 472 RemoveConstProperty("MR_ANGNOTHETA"); << 473 } << 474 529 475 // Adds a new Theta angle number 530 // Adds a new Theta angle number 476 AddConstProperty("MR_ANGNOTHETA", (G4double) 531 AddConstProperty("MR_ANGNOTHETA", (G4double)AngNoTheta); 477 532 >> 533 //G4cout << "AngNoTheta: " << AngNoTheta << G4endl; >> 534 478 // Removes an existing Phi angle number 535 // Removes an existing Phi angle number 479 if (ConstPropertyExists("MR_ANGNOPHI")) { << 536 if (ConstPropertyExists("MR_ANGNOPHI")) RemoveConstProperty("MR_ANGNOPHI"); 480 RemoveConstProperty("MR_ANGNOPHI"); << 481 } << 482 537 483 // Adds a new Phi angle number 538 // Adds a new Phi angle number 484 AddConstProperty("MR_ANGNOPHI", (G4double)An 539 AddConstProperty("MR_ANGNOPHI", (G4double)AngNoPhi); 485 540 >> 541 //G4cout << "AngNoPhi: " << AngNoPhi << G4endl; >> 542 486 // Removes an existing angular cut 543 // Removes an existing angular cut 487 if (ConstPropertyExists("MR_ANGCUT")) { << 544 if (ConstPropertyExists("MR_ANGCUT")) RemoveConstProperty("MR_ANGCUT"); 488 RemoveConstProperty("MR_ANGCUT"); << 489 } << 490 545 491 // Adds a new angle number 546 // Adds a new angle number 492 AddConstProperty("MR_ANGCUT", AngularCut); 547 AddConstProperty("MR_ANGCUT", AngularCut); 493 548 >> 549 //G4cout << "AngularCut: " << AngularCut/degree << "degree" << G4endl; >> 550 494 // Starts the lookup table calculation 551 // Starts the lookup table calculation >> 552 495 ComputeMicroRoughnessTables(); 553 ComputeMicroRoughnessTables(); >> 554 >> 555 //G4cout << " *** DONE! ***" << G4endl; 496 } 556 } 497 557