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 // >> 27 // 26 ////////////////////////////////////////////// 28 //////////////////////////////////////////////////////////////////////// 27 // Optical Surface Class Implementation 29 // Optical Surface Class Implementation 28 ////////////////////////////////////////////// 30 //////////////////////////////////////////////////////////////////////// 29 // 31 // 30 // File: G4OpticalSurface.cc 32 // File: G4OpticalSurface.cc 31 // Description: An optical surface class for u 33 // Description: An optical surface class for use in G4OpBoundaryProcess 32 // Version: 2.0 34 // Version: 2.0 33 // Created: 1997-06-26 35 // Created: 1997-06-26 34 // Author: Peter Gumplinger 36 // Author: Peter Gumplinger 35 // updated: 2017-02-24 Mariele Stockhoff a << 37 // updated: 2017-02-24 Mariele Stockhoff add DAVIS model >> 38 // mail: gum@triumf.ca 36 ////////////////////////////////////////////// 39 //////////////////////////////////////////////////////////////////////// 37 40 38 #include "G4OpticalSurface.hh" << 41 #include <iostream> >> 42 #include <fstream> 39 43 >> 44 //#include "G4ios.hh" 40 #include "globals.hh" 45 #include "globals.hh" >> 46 #include "G4OpticalSurface.hh" 41 47 42 #include <zlib.h> << 48 ///////////////////////// 43 << 49 // Class Implementation 44 #include <fstream> << 50 ///////////////////////// 45 #include <iostream> << 51 >> 52 ////////////// >> 53 // Operators >> 54 ////////////// 46 55 47 G4OpticalSurface& G4OpticalSurface::operator=( 56 G4OpticalSurface& G4OpticalSurface::operator=(const G4OpticalSurface& right) 48 { 57 { 49 if (this != &right) { << 58 if (this != &right) 50 theName = right.theName; << 59 { 51 theType = right.theType; << 60 theName = right.theName; 52 theModel = right.theModel; << 61 theType = right.theType; 53 theFinish = right.theFinish; << 62 theModel = right.theModel; 54 sigma_alpha = right.sigma_alpha; << 63 theFinish = right.theFinish; 55 polish = right.polish; << 64 sigma_alpha = right.sigma_alpha; 56 theMaterialPropertiesTable = right.theMate << 65 polish = right.polish; 57 << 66 theMaterialPropertiesTable = right.theMaterialPropertiesTable; 58 delete[] AngularDistribution; << 67 59 AngularDistribution = new G4float[incident << 68 if (AngularDistribution) delete [] AngularDistribution; 60 *(AngularDistribution) = *(right.AngularDi << 69 AngularDistribution = 61 << 70 new G4float[incidentIndexMax*thetaIndexMax*phiIndexMax]; 62 delete[] AngularDistributionLUT; << 71 *(AngularDistribution) = *(right.AngularDistribution); 63 AngularDistributionLUT = new G4float[index << 72 64 *(AngularDistributionLUT) = *(right.Angula << 73 if (AngularDistributionLUT) delete [] AngularDistributionLUT; 65 << 74 AngularDistributionLUT = new G4float[indexmax]; 66 delete[] Reflectivity; << 75 *(AngularDistributionLUT) = *(right.AngularDistributionLUT); 67 Reflectivity = new G4float[RefMax]; << 76 68 *(Reflectivity) = *(right.Reflectivity); << 77 if (Reflectivity) delete [] Reflectivity; 69 << 78 Reflectivity = new G4float[RefMax]; 70 delete DichroicVector; << 79 *(Reflectivity) = *(right.Reflectivity); 71 DichroicVector = new G4Physics2DVector(); << 80 72 *DichroicVector = *(right.DichroicVector); << 81 if (DichroicVector) delete DichroicVector; 73 } << 82 DichroicVector = new G4Physics2DVector(); >> 83 *DichroicVector = *(right.DichroicVector); >> 84 } 74 return *this; 85 return *this; 75 } 86 } 76 87 77 G4OpticalSurface::G4OpticalSurface(const G4Str << 88 ///////////////// 78 G4OpticalSurfaceFinish finish, G4SurfaceType << 89 // Constructors 79 : G4SurfaceProperty(name, type), theModel(mo << 90 ///////////////// 80 { << 91 81 AngularDistribution = nullptr; << 92 G4OpticalSurface::G4OpticalSurface(const G4String& name, 82 << 93 G4OpticalSurfaceModel model, 83 AngularDistributionLUT = nullptr; << 94 G4OpticalSurfaceFinish finish, 84 Reflectivity = nullptr; << 95 G4SurfaceType type, 85 << 96 G4double value) 86 DichroicVector = nullptr; << 97 : G4SurfaceProperty(name,type), 87 << 98 theModel(model), 88 switch (theModel) { << 99 theFinish(finish), 89 case glisur: << 100 theMaterialPropertiesTable(0) 90 polish = value; << 101 { 91 sigma_alpha = 0.0; << 102 if ( model == glisur ){ 92 break; << 103 polish = value; 93 case LUT: << 104 sigma_alpha = 0.0; 94 case dichroic: << 105 } 95 case DAVIS: << 106 else if ( model == unified ) { 96 ReadDataFile(); << 107 sigma_alpha = value; 97 // fall through << 108 polish = 0.0; 98 case unified: << 109 } 99 sigma_alpha = value; << 110 else if ( model == LUT ) { 100 polish = 0.0; << 111 sigma_alpha = value; 101 break; << 112 polish = 0.0; 102 default: << 113 } 103 G4Exception("G4OpticalSurface::G4Optical << 114 else if ( model == dichroic ) { 104 "Constructor called with INVALID model << 115 sigma_alpha = value; 105 } << 116 polish = 0.0; >> 117 } >> 118 else if ( model == DAVIS ) { >> 119 sigma_alpha = value; >> 120 polish = 0.0; >> 121 } >> 122 else { >> 123 G4Exception("G4OpticalSurface::G4OpticalSurface()", "mat309", >> 124 FatalException, >> 125 "Constructor called with INVALID model."); >> 126 } >> 127 >> 128 AngularDistribution = NULL; >> 129 >> 130 AngularDistributionLUT = NULL; >> 131 Reflectivity = NULL; >> 132 >> 133 DichroicVector = NULL; >> 134 >> 135 if (type == dielectric_LUT) { >> 136 AngularDistribution = >> 137 new G4float[incidentIndexMax*thetaIndexMax*phiIndexMax]; >> 138 ReadLUTFile(); >> 139 } >> 140 >> 141 if (type == dielectric_LUTDAVIS) { >> 142 AngularDistributionLUT = new G4float[indexmax]; >> 143 ReadLUTDAVISFile(); >> 144 >> 145 Reflectivity = new G4float[RefMax]; >> 146 ReadReflectivityLUTFile(); >> 147 } >> 148 >> 149 if (type == dielectric_dichroic) { >> 150 DichroicVector = new G4Physics2DVector(); >> 151 ReadDichroicFile(); >> 152 } 106 } 153 } 107 154 108 G4OpticalSurface::~G4OpticalSurface() 155 G4OpticalSurface::~G4OpticalSurface() 109 { 156 { 110 delete[] AngularDistribution; << 157 if (AngularDistribution) delete [] AngularDistribution; 111 158 112 delete[] AngularDistributionLUT; << 159 if (AngularDistributionLUT) delete [] AngularDistributionLUT; >> 160 if (Reflectivity) delete [] Reflectivity; 113 161 114 delete[] Reflectivity; << 162 if (DichroicVector) delete DichroicVector; 115 << 116 delete DichroicVector; << 117 } 163 } 118 164 119 G4OpticalSurface::G4OpticalSurface(const G4Opt << 165 G4OpticalSurface::G4OpticalSurface(const G4OpticalSurface &right) 120 : G4SurfaceProperty(right.theName, right.the << 166 : G4SurfaceProperty(right.theName,right.theType) 121 { 167 { 122 *this = right; << 168 *this = right; 123 this->theName = right.theName; << 169 this->theName = right.theName; 124 this->theType = right.theType; << 170 this->theType = right.theType; 125 this->theModel = right.theModel; << 171 this->theModel = right.theModel; 126 this->theFinish = right.theFinish; << 172 this->theFinish = right.theFinish; 127 this->sigma_alpha = right.sigma_alpha; << 173 this->sigma_alpha = right.sigma_alpha; 128 this->polish = right.polish; << 174 this->polish = right.polish; 129 this->theMaterialPropertiesTable = right.the << 175 this->theMaterialPropertiesTable = right.theMaterialPropertiesTable; 130 << 176 131 delete[] AngularDistribution; << 177 if (this->AngularDistribution) delete [] AngularDistribution; 132 this->AngularDistribution = new G4float[inci << 178 this->AngularDistribution = 133 *(this->AngularDistribution) = *(right.Angul << 179 new G4float[incidentIndexMax*thetaIndexMax*phiIndexMax]; 134 << 180 *(this->AngularDistribution) = *(right.AngularDistribution); 135 delete[] AngularDistributionLUT; << 181 136 this->AngularDistributionLUT = new G4float[i << 182 if (this->AngularDistributionLUT) delete [] AngularDistributionLUT; 137 *(this->AngularDistributionLUT) = *(right.An << 183 this->AngularDistributionLUT = new G4float[indexmax]; >> 184 *(this->AngularDistributionLUT) = *(right.AngularDistributionLUT); >> 185 if (this->Reflectivity) delete [] Reflectivity; >> 186 this->Reflectivity = new G4float[RefMax]; >> 187 *(this->Reflectivity) = *(right.Reflectivity); 138 188 139 delete[] Reflectivity; << 189 if (this->DichroicVector) delete DichroicVector; 140 this->Reflectivity = new G4float[RefMax]; << 190 this->DichroicVector = new G4Physics2DVector(); 141 *(this->Reflectivity) = *(right.Reflectivity << 191 *(this->DichroicVector) = *(right.DichroicVector); 142 << 143 delete DichroicVector; << 144 this->DichroicVector = new G4Physics2DVector << 145 *(this->DichroicVector) = *(right.DichroicVe << 146 } 192 } 147 193 148 G4bool G4OpticalSurface::operator==(const G4Op << 194 G4bool G4OpticalSurface::operator==(const G4OpticalSurface &right) const 149 { 195 { 150 return (this == (G4OpticalSurface*)&right); << 196 return (this == (G4OpticalSurface *) &right); 151 } 197 } 152 198 153 G4bool G4OpticalSurface::operator!=(const G4Op << 199 G4bool G4OpticalSurface::operator!=(const G4OpticalSurface &right) const 154 { 200 { 155 return (this != (G4OpticalSurface*)&right); << 201 return (this != (G4OpticalSurface *) &right); 156 } 202 } >> 203 //////////// >> 204 // Methods >> 205 //////////// 157 206 158 G4int G4OpticalSurface::GetInmax() const { ret << 207 G4int G4OpticalSurface::GetInmax(void) const 159 << 208 { 160 G4int G4OpticalSurface::GetLUTbins() const { r << 209 return indexmax; >> 210 } 161 211 162 G4int G4OpticalSurface::GetRefMax() const { re << 212 G4int G4OpticalSurface::GetLUTbins(void) const >> 213 { >> 214 return LUTbins; >> 215 } 163 216 164 G4int G4OpticalSurface::GetThetaIndexMax() con << 217 G4int G4OpticalSurface::GetRefMax(void) const >> 218 { >> 219 return RefMax; >> 220 } >> 221 >> 222 G4int G4OpticalSurface::GetThetaIndexMax(void) const >> 223 { >> 224 return thetaIndexMax; >> 225 } 165 226 166 G4int G4OpticalSurface::GetPhiIndexMax() const << 227 G4int G4OpticalSurface::GetPhiIndexMax(void) const >> 228 { >> 229 return phiIndexMax; >> 230 } 167 231 168 void G4OpticalSurface::DumpInfo() const << 232 void G4OpticalSurface::DumpInfo() const 169 { 233 { 170 // Dump info for surface << 234 // Dump info for surface 171 235 172 G4cout << " Surface type = " << G4int(the << 236 G4cout << 173 << " Surface finish = " << G4int(the << 237 " Surface type = " << G4int(theType) << G4endl << 174 << " Surface model = " << G4int(the << 238 " Surface finish = " << G4int(theFinish) << G4endl << 175 << G4endl << " ----------------- " < << 239 " Surface model = " << G4int(theModel) << G4endl; 176 240 177 if (theModel == glisur) { << 241 G4cout << G4endl; 178 G4cout << " polish: " << polish << G4endl; << 242 179 } << 243 G4cout << " Surface parameter " << G4endl; 180 else { << 244 G4cout << " ----------------- " << G4endl; 181 G4cout << " sigma_alpha: " << sigma_alpha << 245 if (theModel == glisur ){ 182 } << 246 G4cout << polish << G4endl; 183 G4cout << G4endl; << 247 } >> 248 else if (theModel == LUT ){ >> 249 G4cout << sigma_alpha << G4endl; >> 250 } >> 251 else { >> 252 G4cout << sigma_alpha << G4endl; >> 253 } >> 254 G4cout << G4endl; 184 } 255 } 185 256 186 void G4OpticalSurface::SetType(const G4Surface 257 void G4OpticalSurface::SetType(const G4SurfaceType& type) 187 { 258 { 188 theType = type; 259 theType = type; 189 ReadDataFile(); << 260 >> 261 if (type == dielectric_LUTDAVIS) { >> 262 if (!AngularDistributionLUT) AngularDistributionLUT = >> 263 new G4float[indexmax]; >> 264 ReadLUTDAVISFile(); >> 265 } >> 266 >> 267 if (type == dielectric_LUT) { >> 268 if (!AngularDistribution) AngularDistribution = >> 269 new G4float[incidentIndexMax*thetaIndexMax*phiIndexMax]; >> 270 >> 271 ReadLUTFile(); >> 272 } >> 273 if (type == dielectric_dichroic) { >> 274 if (!DichroicVector) DichroicVector = new G4Physics2DVector(); >> 275 ReadDichroicFile(); >> 276 } 190 } 277 } 191 278 192 void G4OpticalSurface::SetFinish(const G4Optic 279 void G4OpticalSurface::SetFinish(const G4OpticalSurfaceFinish finish) 193 { 280 { 194 theFinish = finish; 281 theFinish = finish; 195 ReadDataFile(); << 282 if (theType == dielectric_LUT) { 196 } << 283 if (!AngularDistribution) AngularDistribution = 197 << 284 new G4float[incidentIndexMax*thetaIndexMax*phiIndexMax]; 198 void G4OpticalSurface::ReadDataFile() << 285 ReadLUTFile(); 199 { << 286 } 200 // type and finish can be set in either orde << 287 201 // for consistency. Need to read file on set << 288 if (theType == dielectric_LUTDAVIS) { 202 switch (theType) { << 289 if (!AngularDistributionLUT) AngularDistributionLUT = 203 case dielectric_LUT: << 290 new G4float[indexmax]; 204 if (AngularDistribution == nullptr) { << 291 ReadLUTDAVISFile(); 205 AngularDistribution = new G4float[inci << 292 206 } << 293 if (!Reflectivity) Reflectivity = new G4float[RefMax]; 207 ReadLUTFile(); << 294 ReadReflectivityLUTFile(); 208 break; << 295 } 209 case dielectric_LUTDAVIS: << 296 210 if (AngularDistributionLUT == nullptr) { << 297 if (theType == dielectric_dichroic) { 211 AngularDistributionLUT = new G4float[i << 298 if (!DichroicVector) DichroicVector = new G4Physics2DVector(); 212 } << 299 ReadDichroicFile(); 213 ReadLUTDAVISFile(); << 214 << 215 if (Reflectivity == nullptr) { << 216 Reflectivity = new G4float[RefMax]; << 217 } << 218 ReadReflectivityLUTFile(); << 219 break; << 220 case dielectric_dichroic: << 221 if (DichroicVector == nullptr) { << 222 DichroicVector = new G4Physics2DVector << 223 } << 224 ReadDichroicFile(); << 225 break; << 226 default: << 227 break; << 228 } 300 } 229 } 301 } 230 302 231 void G4OpticalSurface::ReadLUTFile() 303 void G4OpticalSurface::ReadLUTFile() 232 { 304 { 233 G4String readLUTFileName; << 305 G4String readLUTFileName = " "; >> 306 >> 307 if (theFinish == polishedlumirrorglue) { >> 308 readLUTFileName = "PolishedLumirrorGlue.dat"; >> 309 } >> 310 else if (theFinish == polishedlumirrorair) { >> 311 readLUTFileName = "PolishedLumirror.dat"; >> 312 } >> 313 else if (theFinish == polishedteflonair) { >> 314 readLUTFileName = "PolishedTeflon.dat"; >> 315 } >> 316 else if (theFinish == polishedtioair) { >> 317 readLUTFileName = "PolishedTiO.dat"; >> 318 } >> 319 else if (theFinish == polishedtyvekair) { >> 320 readLUTFileName = "PolishedTyvek.dat"; >> 321 } >> 322 else if (theFinish == polishedvm2000glue) { >> 323 readLUTFileName = "PolishedVM2000Glue.dat"; >> 324 } >> 325 else if (theFinish == polishedvm2000air) { >> 326 readLUTFileName = "PolishedVM2000.dat"; >> 327 } >> 328 else if (theFinish == etchedlumirrorglue) { >> 329 readLUTFileName = "EtchedLumirrorGlue.dat"; >> 330 } >> 331 else if (theFinish == etchedlumirrorair) { >> 332 readLUTFileName = "EtchedLumirror.dat"; >> 333 } >> 334 else if (theFinish == etchedteflonair) { >> 335 readLUTFileName = "EtchedTeflon.dat"; >> 336 } >> 337 else if (theFinish == etchedtioair) { >> 338 readLUTFileName = "EtchedTiO.dat"; >> 339 } >> 340 else if (theFinish == etchedtyvekair) { >> 341 readLUTFileName = "EtchedTyvek.dat"; >> 342 } >> 343 else if (theFinish == etchedvm2000glue) { >> 344 readLUTFileName = "EtchedVM2000Glue.dat"; >> 345 } >> 346 else if (theFinish == etchedvm2000air) { >> 347 readLUTFileName = "EtchedVM2000.dat"; >> 348 } >> 349 else if (theFinish == groundlumirrorglue) { >> 350 readLUTFileName = "GroundLumirrorGlue.dat"; >> 351 } >> 352 else if (theFinish == groundlumirrorair) { >> 353 readLUTFileName = "GroundLumirror.dat"; >> 354 } >> 355 else if (theFinish == groundteflonair) { >> 356 readLUTFileName = "GroundTeflon.dat"; >> 357 } >> 358 else if (theFinish == groundtioair) { >> 359 readLUTFileName = "GroundTiO.dat"; >> 360 } >> 361 else if (theFinish == groundtyvekair) { >> 362 readLUTFileName = "GroundTyvek.dat"; >> 363 } >> 364 else if (theFinish == groundvm2000glue) { >> 365 readLUTFileName = "GroundVM2000Glue.dat"; >> 366 } >> 367 else if (theFinish == groundvm2000air) { >> 368 readLUTFileName = "GroundVM2000.dat"; >> 369 } 234 370 235 switch (theFinish) { << 371 if (readLUTFileName == " ") return; 236 case polishedlumirrorglue: << 372 237 readLUTFileName = "PolishedLumirrorGlue. << 373 char* path = std::getenv("G4REALSURFACEDATA"); 238 break; << 374 if (!path) { 239 case polishedlumirrorair: << 375 G4String excep = 240 readLUTFileName = "PolishedLumirror.z"; << 376 "G4OpBoundaryProcess - G4REALSURFACEDATA environment variable not set"; 241 break; << 377 G4Exception("G4OpticalSurface::ReadLUTFile()", "mat310", 242 case polishedteflonair: << 378 FatalException, excep); 243 readLUTFileName = "PolishedTeflon.z"; << 379 return; 244 break; << 380 } 245 case polishedtioair: << 381 G4String pathString(path); 246 readLUTFileName = "PolishedTiO.z"; << 382 247 break; << 383 readLUTFileName = pathString + "/" + readLUTFileName; 248 case polishedtyvekair: << 384 249 readLUTFileName = "PolishedTyvek.z"; << 385 std::ifstream readLUTFileHandle(readLUTFileName, std::ios::in); 250 break; << 386 251 case polishedvm2000glue: << 387 if (readLUTFileHandle) { 252 readLUTFileName = "PolishedVM2000Glue.z" << 388 G4int idxmax = incidentIndexMax*thetaIndexMax*phiIndexMax; 253 break; << 389 for (G4int i = 0; i<idxmax; i++) { 254 case polishedvm2000air: << 390 if (readLUTFileHandle.eof()) break; 255 readLUTFileName = "PolishedVM2000.z"; << 391 readLUTFileHandle >> AngularDistribution[i]; 256 break; << 392 } 257 case etchedlumirrorglue: << 393 if (!readLUTFileHandle.bad()) { 258 readLUTFileName = "EtchedLumirrorGlue.z" << 394 G4cout <<"LUT - data file: "<< readLUTFileName <<" read in! "<< G4endl; 259 break; << 395 } 260 case etchedlumirrorair: << 396 else { 261 readLUTFileName = "EtchedLumirror.z"; << 397 G4String excep="LUT - data file: "+readLUTFileName+" not read propery"; 262 break; << 398 G4Exception("G4OpticalSurface::ReadLUTFile()", "mat312", 263 case etchedteflonair: << 399 FatalException, excep); 264 readLUTFileName = "EtchedTeflon.z"; << 400 return; 265 break; << 401 } 266 case etchedtioair: << 267 readLUTFileName = "EtchedTiO.z"; << 268 break; << 269 case etchedtyvekair: << 270 readLUTFileName = "EtchedTyvek.z"; << 271 break; << 272 case etchedvm2000glue: << 273 readLUTFileName = "EtchedVM2000Glue.z"; << 274 break; << 275 case etchedvm2000air: << 276 readLUTFileName = "EtchedVM2000.z"; << 277 break; << 278 case groundlumirrorglue: << 279 readLUTFileName = "GroundLumirrorGlue.z" << 280 break; << 281 case groundlumirrorair: << 282 readLUTFileName = "GroundLumirror.z"; << 283 break; << 284 case groundteflonair: << 285 readLUTFileName = "GroundTeflon.z"; << 286 break; << 287 case groundtioair: << 288 readLUTFileName = "GroundTiO.z"; << 289 break; << 290 case groundtyvekair: << 291 readLUTFileName = "GroundTyvek.z"; << 292 break; << 293 case groundvm2000glue: << 294 readLUTFileName = "GroundVM2000Glue.z"; << 295 break; << 296 case groundvm2000air: << 297 readLUTFileName = "GroundVM2000.z"; << 298 break; << 299 default: << 300 return; << 301 } << 302 << 303 std::istringstream iss; << 304 ReadCompressedFile(readLUTFileName, iss); << 305 << 306 size_t idxmax = incidentIndexMax * thetaInde << 307 for (size_t i = 0; i < idxmax; ++i) { << 308 iss >> AngularDistribution[i]; << 309 } 402 } 310 G4cout << "LUT - data file: " << readLUTFile << 403 else { >> 404 G4String excep ="LUT - data file: "+readLUTFileName+" not found"; >> 405 G4Exception("G4OpticalSurface::ReadLUTFile()", "mat311", >> 406 FatalException, excep); >> 407 return; >> 408 } >> 409 readLUTFileHandle.close(); 311 } 410 } 312 411 313 void G4OpticalSurface::ReadLUTDAVISFile() 412 void G4OpticalSurface::ReadLUTDAVISFile() 314 { 413 { 315 G4String readLUTDAVISFileName; << 414 G4String readLUTDAVISFileName = " "; 316 415 317 switch (theFinish) { << 416 if (theFinish == Rough_LUT) { 318 case Rough_LUT: << 417 readLUTDAVISFileName = "Rough_LUT.dat"; 319 readLUTDAVISFileName = "Rough_LUT.z"; << 418 }else if (theFinish == RoughTeflon_LUT) { 320 break; << 419 readLUTDAVISFileName = "RoughTeflon_LUT.dat"; 321 case RoughTeflon_LUT: << 420 }else if (theFinish == RoughESR_LUT) { 322 readLUTDAVISFileName = "RoughTeflon_LUT. << 421 readLUTDAVISFileName = "RoughESR_LUT.dat"; 323 break; << 422 }else if (theFinish == RoughESRGrease_LUT) { 324 case RoughESR_LUT: << 423 readLUTDAVISFileName = "RoughESRGrease_LUT.dat"; 325 readLUTDAVISFileName = "RoughESR_LUT.z"; << 424 }else if (theFinish == Polished_LUT) { 326 break; << 425 readLUTDAVISFileName = "Polished_LUT.dat"; 327 case RoughESRGrease_LUT: << 426 }else if (theFinish == PolishedTeflon_LUT) { 328 readLUTDAVISFileName = "RoughESRGrease_L << 427 readLUTDAVISFileName = "PolishedTeflon_LUT.dat"; 329 break; << 428 }else if (theFinish == PolishedESR_LUT) { 330 case Polished_LUT: << 429 readLUTDAVISFileName = "PolishedESR_LUT.dat"; 331 readLUTDAVISFileName = "Polished_LUT.z"; << 430 }else if (theFinish == PolishedESRGrease_LUT) { 332 break; << 431 readLUTDAVISFileName = "PolishedESRGrease_LUT.dat"; 333 case PolishedTeflon_LUT: << 432 }else if (theFinish == Detector_LUT) { 334 readLUTDAVISFileName = "PolishedTeflon_L << 433 readLUTDAVISFileName = "Detector_LUT.dat"; 335 break; << 434 } 336 case PolishedESR_LUT: << 435 337 readLUTDAVISFileName = "PolishedESR_LUT. << 436 if (readLUTDAVISFileName == " ") return; 338 break; << 437 339 case PolishedESRGrease_LUT: << 438 char* path = std::getenv("G4REALSURFACEDATA"); 340 readLUTDAVISFileName = "PolishedESRGreas << 439 if (!path) { 341 break; << 440 G4String excep = 342 case Detector_LUT: << 441 "G4OpBoundaryProcess - G4REALSURFACEDATA environment variable not set"; 343 readLUTDAVISFileName = "Detector_LUT.z"; << 442 G4Exception("G4OpticalSurface::ReadLUTDAVISFile()", "mat310", 344 break; << 443 FatalException, excep); 345 default: << 444 return; 346 return; << 445 } >> 446 G4String pathString(path); >> 447 >> 448 readLUTDAVISFileName = pathString + "/" + readLUTDAVISFileName; >> 449 >> 450 std::ifstream readLUTDAVISFileHandle(readLUTDAVISFileName, std::ios::in); >> 451 >> 452 if (readLUTDAVISFileHandle) { >> 453 >> 454 for (G4int i = 0; i<indexmax; i++) { >> 455 if (readLUTDAVISFileHandle.eof()) break; >> 456 readLUTDAVISFileHandle >> AngularDistributionLUT[i]; >> 457 } >> 458 if (!readLUTDAVISFileHandle.bad()) { >> 459 G4cout <<"LUT DAVIS - data file: " << >> 460 readLUTDAVISFileName << " read in! " << G4endl; >> 461 } >> 462 else { >> 463 G4String excep="LUT DAVIS- data file: "+readLUTDAVISFileName+" not read propery"; >> 464 G4Exception("G4OpticalSurface::ReadLUTFile()", "mat312", >> 465 FatalException, excep); >> 466 return; >> 467 } 347 } 468 } 348 << 469 else { 349 std::istringstream iss; << 470 G4String excep ="LUT DAVIS- data file: "+readLUTDAVISFileName+" not found"; 350 ReadCompressedFile(readLUTDAVISFileName, iss << 471 G4Exception("G4OpticalSurface::ReadLUTFile()", "mat311", 351 << 472 FatalException, excep); 352 for (size_t i = 0; i < indexmax; ++i) { << 473 return; 353 iss >> AngularDistributionLUT[i]; << 354 } 474 } 355 G4cout << "LUT DAVIS - data file: " << readL << 475 readLUTDAVISFileHandle.close(); 356 } 476 } 357 477 358 void G4OpticalSurface::ReadReflectivityLUTFile 478 void G4OpticalSurface::ReadReflectivityLUTFile() 359 { 479 { 360 G4String readReflectivityLUTFileName; << 480 G4String readReflectivityLUTFileName = " "; 361 481 362 switch (theFinish) { << 482 if (theFinish == Rough_LUT) { 363 case Rough_LUT: << 483 readReflectivityLUTFileName = "Rough_LUTR.dat"; 364 readReflectivityLUTFileName = "Rough_LUT << 365 break; << 366 case RoughTeflon_LUT: << 367 readReflectivityLUTFileName = "RoughTefl << 368 break; << 369 case RoughESR_LUT: << 370 readReflectivityLUTFileName = "RoughESR_ << 371 break; << 372 case RoughESRGrease_LUT: << 373 readReflectivityLUTFileName = "RoughESRG << 374 break; << 375 case Polished_LUT: << 376 readReflectivityLUTFileName = "Polished_ << 377 break; << 378 case PolishedTeflon_LUT: << 379 readReflectivityLUTFileName = "PolishedT << 380 break; << 381 case PolishedESR_LUT: << 382 readReflectivityLUTFileName = "PolishedE << 383 break; << 384 case PolishedESRGrease_LUT: << 385 readReflectivityLUTFileName = "PolishedE << 386 break; << 387 case Detector_LUT: << 388 readReflectivityLUTFileName = "Detector_ << 389 break; << 390 default: << 391 return; << 392 } << 393 << 394 std::istringstream iss; << 395 ReadCompressedFile(readReflectivityLUTFileNa << 396 << 397 for (size_t i = 0; i < RefMax; ++i) { << 398 iss >> Reflectivity[i]; << 399 } << 400 G4cout << "LUT DAVIS - reflectivity data fil << 401 << G4endl; << 402 } << 403 << 404 // uncompress one data file into the input str << 405 void G4OpticalSurface::ReadCompressedFile(cons << 406 { << 407 G4String* dataString = nullptr; << 408 G4String path = G4FindDataDir("G4REALSURFACE << 409 G4String compfilename = path + "/" + filenam << 410 // create input stream with binary mode oper << 411 std::ifstream in(compfilename, std::ios::bin << 412 if (in.good()) { << 413 // get current position in the stream (was << 414 G4int fileSize = (G4int)in.tellg(); << 415 // set current position being the beginnin << 416 in.seekg(0, std::ios::beg); << 417 // create (zlib) byte buffer for the data << 418 auto compdata = new Bytef[fileSize]; << 419 while (in) { << 420 in.read((char*)compdata, fileSize); << 421 } << 422 // create (zlib) byte buffer for the uncom << 423 auto complen = (uLongf)(fileSize * 4); << 424 auto uncompdata = new Bytef[complen]; << 425 while (Z_OK != uncompress(uncompdata, &com << 426 // increase uncompressed byte buffer << 427 delete[] uncompdata; << 428 complen *= 2; << 429 uncompdata = new Bytef[complen]; << 430 } << 431 // delete the compressed data buffer << 432 delete[] compdata; << 433 // create a string from uncompressed data << 434 dataString = new G4String((char*)uncompdat << 435 // delete the uncompressed data buffer << 436 delete[] uncompdata; << 437 } 484 } 438 else { << 485 if (theFinish == RoughTeflon_LUT) { 439 G4ExceptionDescription ed; << 486 readReflectivityLUTFileName = "RoughTeflon_LUTR.dat"; 440 ed << "Problem while trying to read " + co << 487 } 441 G4Exception("G4OpticalSurface::ReadCompres << 488 if (theFinish== RoughESR_LUT) { 442 return; << 489 readReflectivityLUTFileName = "RoughESR_LUTR.dat"; >> 490 } >> 491 if (theFinish == RoughESRGrease_LUT) { >> 492 readReflectivityLUTFileName = "RoughESRGrease_LUTR.dat"; >> 493 } >> 494 if (theFinish== Polished_LUT) { >> 495 readReflectivityLUTFileName = "Polished_LUTR.dat"; >> 496 } >> 497 if (theFinish == PolishedTeflon_LUT) { >> 498 readReflectivityLUTFileName = "PolishedTeflon_LUTR.dat"; >> 499 } >> 500 if (theFinish== PolishedESR_LUT) { >> 501 readReflectivityLUTFileName = "PolishedESR_LUTR.dat"; >> 502 } >> 503 if (theFinish== PolishedESRGrease_LUT) { >> 504 readReflectivityLUTFileName = "PolishedESRGrease_LUTR.dat"; >> 505 } >> 506 if (theFinish== Detector_LUT) { >> 507 readReflectivityLUTFileName = "Detector_LUTR.dat"; >> 508 } >> 509 >> 510 if (readReflectivityLUTFileName == " ") return; >> 511 >> 512 char* path = std::getenv("G4REALSURFACEDATA"); >> 513 if (!path) { >> 514 G4String excep = >> 515 "G4OpBoundaryProcess - G4REALSURFACEDATA environment variable not set"; >> 516 G4Exception("G4OpticalSurface::ReadReflectivityLUTFile()", "mat310", >> 517 FatalException, excep); >> 518 return; >> 519 } >> 520 G4String pathString(path); >> 521 >> 522 readReflectivityLUTFileName = pathString + "/" + readReflectivityLUTFileName; >> 523 >> 524 std::ifstream readReflectivityLUTFileHandle(readReflectivityLUTFileName, >> 525 std::ios::in); >> 526 >> 527 if (readReflectivityLUTFileHandle) { >> 528 for (G4int i = 0; i<RefMax; i++) { >> 529 if (readReflectivityLUTFileHandle.eof()) break; >> 530 readReflectivityLUTFileHandle >> Reflectivity[i]; >> 531 } >> 532 if (!readReflectivityLUTFileHandle.bad()) { >> 533 G4cout << "Reflectivity LUT DAVIS - data file: " << >> 534 readReflectivityLUTFileName << " read in! " << G4endl; >> 535 } >> 536 else { >> 537 G4String excep="Reflectivity LUT DAVIS- data file: "+readReflectivityLUTFileName+" not read propery"; >> 538 G4Exception("G4OpticalSurface::ReadReflectivityLUTFile()", "mat312", >> 539 FatalException, excep); >> 540 return; >> 541 } 443 } 542 } 444 // create the input string stream from the d << 543 else { 445 if (dataString != nullptr) { << 544 G4String excep ="Reflectivity LUT DAVIS - data file: "+readReflectivityLUTFileName+" not found"; 446 iss.str(*dataString); << 545 G4Exception("G4OpticalSurface::ReadReflectivityLUTFile()", "mat311", 447 in.close(); << 546 FatalException, excep); 448 delete dataString; << 547 return; 449 G4cout << "G4OpticalSurface: data file " < << 450 } 548 } >> 549 readReflectivityLUTFileHandle.close(); 451 } 550 } 452 551 453 void G4OpticalSurface::ReadDichroicFile() 552 void G4OpticalSurface::ReadDichroicFile() 454 { 553 { 455 const char* datadir = G4FindDataDir("G4DICHR << 554 const char* datadir = std::getenv("G4DICHROICDATA"); 456 555 457 if (datadir == nullptr) { << 556 if(!datadir) { 458 G4Exception("G4OpticalSurface::ReadDichroi << 557 G4Exception("G4OpticalSurface::ReadDichroicFile()","mat313", 459 "Environment variable G4DICHROICDATA not << 558 FatalException,"Environment variable G4DICHROICDATA not defined"); 460 return; 559 return; 461 } 560 } 462 561 463 std::ostringstream ost; 562 std::ostringstream ost; 464 ost << datadir; 563 ost << datadir; 465 std::ifstream fin(ost.str().c_str()); 564 std::ifstream fin(ost.str().c_str()); 466 if (! fin.is_open()) { << 565 if( !fin.is_open()) { 467 G4ExceptionDescription ed; 566 G4ExceptionDescription ed; 468 ed << "Dichroic surface data file <" << os << 567 ed << "Dichroic surface data file <" << ost.str().c_str() 469 G4Exception("G4OpticalSurface::ReadDichroi << 568 << "> is not opened!" << G4endl; >> 569 G4Exception("G4OpticalSurface::ReadDichroicFile()","mat314", >> 570 FatalException,ed," "); 470 return; 571 return; 471 } 572 } 472 573 473 if (! (DichroicVector->Retrieve(fin))) { << 574 if( !(DichroicVector->Retrieve(fin)) ) { 474 G4ExceptionDescription ed; 575 G4ExceptionDescription ed; 475 ed << "Dichroic surface data file <" << os << 576 ed << "Dichroic surface data file <" << ost.str().c_str() 476 G4Exception("G4OpticalSurface::ReadDichroi << 577 << "> is not opened!" << G4endl; >> 578 G4Exception("G4OpticalSurface::ReadDichroicFile()","mat315", >> 579 FatalException,ed," "); 477 return; 580 return; 478 } 581 } 479 582 480 // DichroicVector->SetBicubicInterpolation( << 583 // DichroicVector->SetBicubicInterpolation(true); 481 584 482 G4cout << " *** Dichroic surface data file * 585 G4cout << " *** Dichroic surface data file *** " << G4endl; 483 586 484 auto numberOfXNodes = (G4int)DichroicVector- << 587 G4int numberOfXNodes = DichroicVector->GetLengthX(); 485 auto numberOfYNodes = (G4int)DichroicVector- << 588 G4int numberOfYNodes = DichroicVector->GetLengthY(); 486 589 487 G4cout << "numberOfXNodes: " << numberOfXNod 590 G4cout << "numberOfXNodes: " << numberOfXNodes << G4endl; 488 G4cout << "numberOfYNodes: " << numberOfYNod 591 G4cout << "numberOfYNodes: " << numberOfYNodes << G4endl; 489 592 490 if (0 > numberOfXNodes || numberOfXNodes >= << 593 if (0 > numberOfXNodes || numberOfXNodes >= INT_MAX) numberOfXNodes = 0; 491 numberOfXNodes = 0; << 594 if (0 > numberOfYNodes || numberOfYNodes >= INT_MAX) numberOfYNodes = 0; 492 } << 493 if (0 > numberOfYNodes || numberOfYNodes >= << 494 numberOfYNodes = 0; << 495 } << 496 595 497 G4PV2DDataVector xVector; << 596 G4PV2DDataVector xVector; 498 G4PV2DDataVector yVector; << 597 G4PV2DDataVector yVector; 499 598 500 xVector.resize(numberOfXNodes, 0.); << 599 xVector.resize(numberOfXNodes,0.); 501 yVector.resize(numberOfYNodes, 0.); << 600 yVector.resize(numberOfYNodes,0.); 502 601 503 for (G4int i = 0; i < numberOfXNodes; ++i) { << 602 for(G4int i = 0; i<numberOfXNodes; ++i) { 504 G4cout << "i: " << DichroicVector->GetX(i) << 603 G4cout << "i: " << DichroicVector->GetX(i) << G4endl; 505 xVector[i] = DichroicVector->GetX(i); << 604 xVector[i] = DichroicVector->GetX(i); 506 } 605 } 507 for (G4int j = 0; j < numberOfYNodes; ++j) { << 606 for(G4int j = 0; j<numberOfYNodes; ++j) { 508 G4cout << "j: " << DichroicVector->GetY(j) << 607 G4cout << "j: " << DichroicVector->GetY(j) << G4endl; 509 yVector[j] = DichroicVector->GetY(j); << 608 yVector[j] = DichroicVector->GetY(j); 510 } 609 } 511 610 512 for (G4int j = 0; j < numberOfYNodes; ++j) { << 611 for(G4int j = 0; j<numberOfYNodes; ++j) { 513 for (G4int i = 0; i < numberOfXNodes; ++i) << 612 for(G4int i = 0; i<numberOfXNodes; ++i) { 514 G4cout << " i: " << i << " j: " << j << << 613 G4cout << " i: " << i << " j: " << j << " " 515 } << 614 << DichroicVector->GetValue(i,j) << G4endl; >> 615 } 516 } 616 } >> 617 >> 618 // G4int idx, idy; >> 619 >> 620 // for(G4int j = 0; j<numberOfYNodes-1; ++j) { >> 621 // G4double y = (yVector[j] + yVector[j+1])/2.; >> 622 // for(G4int i = 0; i<numberOfXNodes-1; ++i) { >> 623 // G4double x = (xVector[i] + xVector[i+1])/2.; >> 624 // G4cout << " x: " << x << " y: " << y << " " >> 625 // << DichroicVector->Value(x,y,idx,idy) << G4endl; >> 626 // } >> 627 // } >> 628 517 } 629 } 518 630