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