| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25
26 ////////////////////////////////////////////////////////////////////////
27 // Optical Surface Class Implementation
28 ////////////////////////////////////////////////////////////////////////
29 //
30 // File: G4OpticalSurface.cc
31 // Description: An optical surface class for use in G4OpBoundaryProcess
32 // Version: 2.0
33 // Created: 1997-06-26
34 // Author: Peter Gumplinger
35 // updated: 2017-02-24 Mariele Stockhoff add DAVIS model
36 ////////////////////////////////////////////////////////////////////////
37
38 #include "G4OpticalSurface.hh"
39
40 #include "globals.hh"
41
42 #include <zlib.h>
43
44 #include <fstream>
45 #include <iostream>
46
47 G4OpticalSurface& G4OpticalSurface::operator=(const G4OpticalSurface& right)
48 {
49 if (this != &right) {
50 theName = right.theName;
51 theType = right.theType;
52 theModel = right.theModel;
53 theFinish = right.theFinish;
54 sigma_alpha = right.sigma_alpha;
55 polish = right.polish;
56 theMaterialPropertiesTable = right.theMaterialPropertiesTable;
57
58 delete[] AngularDistribution;
59 AngularDistribution = new G4float[incidentIndexMax * thetaIndexMax * phiIndexMax];
60 *(AngularDistribution) = *(right.AngularDistribution);
61
62 delete[] AngularDistributionLUT;
63 AngularDistributionLUT = new G4float[indexmax];
64 *(AngularDistributionLUT) = *(right.AngularDistributionLUT);
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;
75 }
76
77 G4OpticalSurface::G4OpticalSurface(const G4String& name, G4OpticalSurfaceModel model,
78 G4OpticalSurfaceFinish finish, G4SurfaceType type, G4double value)
79 : G4SurfaceProperty(name, type), theModel(model), theFinish(finish)
80 {
81 AngularDistribution = nullptr;
82
83 AngularDistributionLUT = nullptr;
84 Reflectivity = nullptr;
85
86 DichroicVector = nullptr;
87
88 switch (theModel) {
89 case glisur:
90 polish = value;
91 sigma_alpha = 0.0;
92 break;
93 case LUT:
94 case dichroic:
95 case DAVIS:
96 ReadDataFile();
97 // fall through
98 case unified:
99 sigma_alpha = value;
100 polish = 0.0;
101 break;
102 default:
103 G4Exception("G4OpticalSurface::G4OpticalSurface()", "mat309", FatalException,
104 "Constructor called with INVALID model.");
105 }
106 }
107
108 G4OpticalSurface::~G4OpticalSurface()
109 {
110 delete[] AngularDistribution;
111
112 delete[] AngularDistributionLUT;
113
114 delete[] Reflectivity;
115
116 delete DichroicVector;
117 }
118
119 G4OpticalSurface::G4OpticalSurface(const G4OpticalSurface& right)
120 : G4SurfaceProperty(right.theName, right.theType)
121 {
122 *this = right;
123 this->theName = right.theName;
124 this->theType = right.theType;
125 this->theModel = right.theModel;
126 this->theFinish = right.theFinish;
127 this->sigma_alpha = right.sigma_alpha;
128 this->polish = right.polish;
129 this->theMaterialPropertiesTable = right.theMaterialPropertiesTable;
130
131 delete[] AngularDistribution;
132 this->AngularDistribution = new G4float[incidentIndexMax * thetaIndexMax * phiIndexMax];
133 *(this->AngularDistribution) = *(right.AngularDistribution);
134
135 delete[] AngularDistributionLUT;
136 this->AngularDistributionLUT = new G4float[indexmax];
137 *(this->AngularDistributionLUT) = *(right.AngularDistributionLUT);
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.DichroicVector);
146 }
147
148 G4bool G4OpticalSurface::operator==(const G4OpticalSurface& right) const
149 {
150 return (this == (G4OpticalSurface*)&right);
151 }
152
153 G4bool G4OpticalSurface::operator!=(const G4OpticalSurface& right) const
154 {
155 return (this != (G4OpticalSurface*)&right);
156 }
157
158 G4int G4OpticalSurface::GetInmax() const { return indexmax; }
159
160 G4int G4OpticalSurface::GetLUTbins() const { return LUTbins; }
161
162 G4int G4OpticalSurface::GetRefMax() const { return RefMax; }
163
164 G4int G4OpticalSurface::GetThetaIndexMax() const { return thetaIndexMax; }
165
166 G4int G4OpticalSurface::GetPhiIndexMax() const { return phiIndexMax; }
167
168 void G4OpticalSurface::DumpInfo() const
169 {
170 // Dump info for surface
171
172 G4cout << " Surface type = " << G4int(theType) << G4endl
173 << " Surface finish = " << G4int(theFinish) << G4endl
174 << " Surface model = " << G4int(theModel) << G4endl << G4endl << " Surface parameter "
175 << G4endl << " ----------------- " << G4endl;
176
177 if (theModel == glisur) {
178 G4cout << " polish: " << polish << G4endl;
179 }
180 else {
181 G4cout << " sigma_alpha: " << sigma_alpha << G4endl;
182 }
183 G4cout << G4endl;
184 }
185
186 void G4OpticalSurface::SetType(const G4SurfaceType& type)
187 {
188 theType = type;
189 ReadDataFile();
190 }
191
192 void G4OpticalSurface::SetFinish(const G4OpticalSurfaceFinish finish)
193 {
194 theFinish = finish;
195 ReadDataFile();
196 }
197
198 void G4OpticalSurface::ReadDataFile()
199 {
200 // type and finish can be set in either order. Thus, we can't check
201 // for consistency. Need to read file on setting either type or finish.
202 switch (theType) {
203 case dielectric_LUT:
204 if (AngularDistribution == nullptr) {
205 AngularDistribution = new G4float[incidentIndexMax * thetaIndexMax * phiIndexMax];
206 }
207 ReadLUTFile();
208 break;
209 case dielectric_LUTDAVIS:
210 if (AngularDistributionLUT == nullptr) {
211 AngularDistributionLUT = new G4float[indexmax];
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 }
229 }
230
231 void G4OpticalSurface::ReadLUTFile()
232 {
233 G4String readLUTFileName;
234
235 switch (theFinish) {
236 case polishedlumirrorglue:
237 readLUTFileName = "PolishedLumirrorGlue.z";
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 * thetaIndexMax * phiIndexMax;
307 for (size_t i = 0; i < idxmax; ++i) {
308 iss >> AngularDistribution[i];
309 }
310 G4cout << "LUT - data file: " << readLUTFileName << " read in! " << G4endl;
311 }
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.z";
323 break;
324 case RoughESR_LUT:
325 readLUTDAVISFileName = "RoughESR_LUT.z";
326 break;
327 case RoughESRGrease_LUT:
328 readLUTDAVISFileName = "RoughESRGrease_LUT.z";
329 break;
330 case Polished_LUT:
331 readLUTDAVISFileName = "Polished_LUT.z";
332 break;
333 case PolishedTeflon_LUT:
334 readLUTDAVISFileName = "PolishedTeflon_LUT.z";
335 break;
336 case PolishedESR_LUT:
337 readLUTDAVISFileName = "PolishedESR_LUT.z";
338 break;
339 case PolishedESRGrease_LUT:
340 readLUTDAVISFileName = "PolishedESRGrease_LUT.z";
341 break;
342 case Detector_LUT:
343 readLUTDAVISFileName = "Detector_LUT.z";
344 break;
345 default:
346 return;
347 }
348
349 std::istringstream iss;
350 ReadCompressedFile(readLUTDAVISFileName, iss);
351
352 for (size_t i = 0; i < indexmax; ++i) {
353 iss >> AngularDistributionLUT[i];
354 }
355 G4cout << "LUT DAVIS - data file: " << readLUTDAVISFileName << " read in! " << G4endl;
356 }
357
358 void G4OpticalSurface::ReadReflectivityLUTFile()
359 {
360 G4String readReflectivityLUTFileName;
361
362 switch (theFinish) {
363 case Rough_LUT:
364 readReflectivityLUTFileName = "Rough_LUTR.z";
365 break;
366 case RoughTeflon_LUT:
367 readReflectivityLUTFileName = "RoughTeflon_LUTR.z";
368 break;
369 case RoughESR_LUT:
370 readReflectivityLUTFileName = "RoughESR_LUTR.z";
371 break;
372 case RoughESRGrease_LUT:
373 readReflectivityLUTFileName = "RoughESRGrease_LUTR.z";
374 break;
375 case Polished_LUT:
376 readReflectivityLUTFileName = "Polished_LUTR.z";
377 break;
378 case PolishedTeflon_LUT:
379 readReflectivityLUTFileName = "PolishedTeflon_LUTR.z";
380 break;
381 case PolishedESR_LUT:
382 readReflectivityLUTFileName = "PolishedESR_LUTR.z";
383 break;
384 case PolishedESRGrease_LUT:
385 readReflectivityLUTFileName = "PolishedESRGrease_LUTR.z";
386 break;
387 case Detector_LUT:
388 readReflectivityLUTFileName = "Detector_LUTR.z";
389 break;
390 default:
391 return;
392 }
393
394 std::istringstream iss;
395 ReadCompressedFile(readReflectivityLUTFileName, iss);
396
397 for (size_t i = 0; i < RefMax; ++i) {
398 iss >> Reflectivity[i];
399 }
400 G4cout << "LUT DAVIS - reflectivity data file: " << readReflectivityLUTFileName << " read in! "
401 << G4endl;
402 }
403
404 // uncompress one data file into the input string stream
405 void G4OpticalSurface::ReadCompressedFile(const G4String& filename, std::istringstream& iss)
406 {
407 G4String* dataString = nullptr;
408 G4String path = G4FindDataDir("G4REALSURFACEDATA");
409 G4String compfilename = path + "/" + filename;
410 // create input stream with binary mode operation and position at end of file
411 std::ifstream in(compfilename, std::ios::binary | std::ios::ate);
412 if (in.good()) {
413 // get current position in the stream (was set to the end)
414 G4int fileSize = (G4int)in.tellg();
415 // set current position being the beginning of the stream
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 uncompressed data
423 auto complen = (uLongf)(fileSize * 4);
424 auto uncompdata = new Bytef[complen];
425 while (Z_OK != uncompress(uncompdata, &complen, compdata, fileSize)) {
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 (will be deallocated by caller)
434 dataString = new G4String((char*)uncompdata, (long)complen);
435 // delete the uncompressed data buffer
436 delete[] uncompdata;
437 }
438 else {
439 G4ExceptionDescription ed;
440 ed << "Problem while trying to read " + compfilename + " data file.\n";
441 G4Exception("G4OpticalSurface::ReadCompressedFile", "mat316", FatalException, ed);
442 return;
443 }
444 // create the input string stream from the data string
445 if (dataString != nullptr) {
446 iss.str(*dataString);
447 in.close();
448 delete dataString;
449 G4cout << "G4OpticalSurface: data file " << compfilename << " successfully read in." << G4endl;
450 }
451 }
452
453 void G4OpticalSurface::ReadDichroicFile()
454 {
455 const char* datadir = G4FindDataDir("G4DICHROICDATA");
456
457 if (datadir == nullptr) {
458 G4Exception("G4OpticalSurface::ReadDichroicFile()", "mat313", FatalException,
459 "Environment variable G4DICHROICDATA not defined");
460 return;
461 }
462
463 std::ostringstream ost;
464 ost << datadir;
465 std::ifstream fin(ost.str().c_str());
466 if (! fin.is_open()) {
467 G4ExceptionDescription ed;
468 ed << "Dichroic surface data file <" << ost.str().c_str() << "> is not opened!" << G4endl;
469 G4Exception("G4OpticalSurface::ReadDichroicFile()", "mat314", FatalException, ed, " ");
470 return;
471 }
472
473 if (! (DichroicVector->Retrieve(fin))) {
474 G4ExceptionDescription ed;
475 ed << "Dichroic surface data file <" << ost.str().c_str() << "> is not opened!" << G4endl;
476 G4Exception("G4OpticalSurface::ReadDichroicFile()", "mat315", FatalException, ed, " ");
477 return;
478 }
479
480 // DichroicVector->SetBicubicInterpolation(true);
481
482 G4cout << " *** Dichroic surface data file *** " << G4endl;
483
484 auto numberOfXNodes = (G4int)DichroicVector->GetLengthX();
485 auto numberOfYNodes = (G4int)DichroicVector->GetLengthY();
486
487 G4cout << "numberOfXNodes: " << numberOfXNodes << G4endl;
488 G4cout << "numberOfYNodes: " << numberOfYNodes << G4endl;
489
490 if (0 > numberOfXNodes || numberOfXNodes >= INT_MAX) {
491 numberOfXNodes = 0;
492 }
493 if (0 > numberOfYNodes || numberOfYNodes >= INT_MAX) {
494 numberOfYNodes = 0;
495 }
496
497 G4PV2DDataVector xVector;
498 G4PV2DDataVector yVector;
499
500 xVector.resize(numberOfXNodes, 0.);
501 yVector.resize(numberOfYNodes, 0.);
502
503 for (G4int i = 0; i < numberOfXNodes; ++i) {
504 G4cout << "i: " << DichroicVector->GetX(i) << G4endl;
505 xVector[i] = DichroicVector->GetX(i);
506 }
507 for (G4int j = 0; j < numberOfYNodes; ++j) {
508 G4cout << "j: " << DichroicVector->GetY(j) << G4endl;
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 << " " << DichroicVector->GetValue(i, j) << G4endl;
515 }
516 }
517 }
518