| 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 // G4MaterialPropertiesTable Implementation
28 ////////////////////////////////////////////////////////////////////////
29 //
30 // File: G4MaterialPropertiesTable.cc
31 // Version: 1.0
32 // Created: 1996-02-08
33 // Author: Juliet Armstrong
34 // Updated: 2005-05-12 add SetGROUPVEL(), courtesy of
35 // Horton-Smith (bug report #741), by P. Gumplinger
36 // 2002-11-05 add named material constants by P. Gumplinger
37 // 1999-11-05 Migration from G4RWTPtrHashDictionary to STL
38 // by John Allison
39 // 1997-03-26 by Peter Gumplinger
40 // > cosmetics (only)
41 //
42 ////////////////////////////////////////////////////////////////////////
43
44 #include "G4MaterialPropertiesTable.hh"
45
46 #include "G4Log.hh"
47 #include "G4OpticalMaterialProperties.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "globals.hh"
50
51 #include <algorithm>
52 #include <cassert>
53
54 #ifdef G4MULTITHREADED
55 # include "G4AutoLock.hh"
56 namespace
57 {
58 G4Mutex materialPropertyTableMutex = G4MUTEX_INITIALIZER;
59 }
60 #endif
61
62 G4MaterialPropertiesTable::G4MaterialPropertiesTable()
63 {
64 // elements of these 2 vectors must be in same order as
65 // the corresponding enums in G4MaterialPropertiesIndex.hh
66 fMatPropNames.assign(kNumberOfPropertyIndex, "");
67 fMatPropNames[kRINDEX] = "RINDEX";
68 fMatPropNames[kREFLECTIVITY] = "REFLECTIVITY";
69 fMatPropNames[kREALRINDEX] = "REALRINDEX";
70 fMatPropNames[kIMAGINARYRINDEX] = "IMAGINARYRINDEX";
71 fMatPropNames[kEFFICIENCY] = "EFFICIENCY";
72 fMatPropNames[kTRANSMITTANCE] = "TRANSMITTANCE";
73 fMatPropNames[kSPECULARLOBECONSTANT] = "SPECULARLOBECONSTANT";
74 fMatPropNames[kSPECULARSPIKECONSTANT] = "SPECULARSPIKECONSTANT";
75 fMatPropNames[kBACKSCATTERCONSTANT] = "BACKSCATTERCONSTANT";
76 fMatPropNames[kGROUPVEL] = "GROUPVEL";
77 fMatPropNames[kMIEHG] = "MIEHG";
78 fMatPropNames[kRAYLEIGH] = "RAYLEIGH";
79 fMatPropNames[kWLSCOMPONENT] = "WLSCOMPONENT";
80 fMatPropNames[kWLSABSLENGTH] = "WLSABSLENGTH";
81 fMatPropNames[kWLSCOMPONENT2] = "WLSCOMPONENT2";
82 fMatPropNames[kWLSABSLENGTH2] = "WLSABSLENGTH2";
83 fMatPropNames[kABSLENGTH] = "ABSLENGTH";
84 fMatPropNames[kPROTONSCINTILLATIONYIELD] = "PROTONSCINTILLATIONYIELD";
85 fMatPropNames[kDEUTERONSCINTILLATIONYIELD] ="DEUTERONSCINTILLATIONYIELD";
86 fMatPropNames[kTRITONSCINTILLATIONYIELD] = "TRITONSCINTILLATIONYIELD";
87 fMatPropNames[kALPHASCINTILLATIONYIELD] = "ALPHASCINTILLATIONYIELD";
88 fMatPropNames[kIONSCINTILLATIONYIELD] = "IONSCINTILLATIONYIELD";
89 fMatPropNames[kELECTRONSCINTILLATIONYIELD] ="ELECTRONSCINTILLATIONYIELD";
90 fMatPropNames[kSCINTILLATIONCOMPONENT1] = "SCINTILLATIONCOMPONENT1";
91 fMatPropNames[kSCINTILLATIONCOMPONENT2] = "SCINTILLATIONCOMPONENT2";
92 fMatPropNames[kSCINTILLATIONCOMPONENT3] = "SCINTILLATIONCOMPONENT3";
93 fMatPropNames[kCOATEDRINDEX] = "COATEDRINDEX";
94
95 fMP.assign(kNumberOfPropertyIndex, nullptr);
96
97 fMatConstPropNames.assign(kNumberOfConstPropertyIndex, "");
98 fMatConstPropNames[kSURFACEROUGHNESS] = "SURFACEROUGHNESS";
99 fMatConstPropNames[kISOTHERMAL_COMPRESSIBILITY] = "ISOTHERMAL_COMPRESSIBILITY";
100 fMatConstPropNames[kRS_SCALE_FACTOR] = "RS_SCALE_FACTOR";
101 fMatConstPropNames[kWLSMEANNUMBERPHOTONS] = "WLSMEANNUMBERPHOTONS";
102 fMatConstPropNames[kWLSTIMECONSTANT] = "WLSTIMECONSTANT";
103 fMatConstPropNames[kWLSMEANNUMBERPHOTONS2] = "WLSMEANNUMBERPHOTONS2";
104 fMatConstPropNames[kWLSTIMECONSTANT2] = "WLSTIMECONSTANT2";
105 fMatConstPropNames[kMIEHG_FORWARD] = "MIEHG_FORWARD";
106 fMatConstPropNames[kMIEHG_BACKWARD] = "MIEHG_BACKWARD";
107 fMatConstPropNames[kMIEHG_FORWARD_RATIO] = "MIEHG_FORWARD_RATIO";
108 fMatConstPropNames[kSCINTILLATIONYIELD] = "SCINTILLATIONYIELD";
109 fMatConstPropNames[kRESOLUTIONSCALE] = "RESOLUTIONSCALE";
110 fMatConstPropNames[kFERMIPOT] = "FERMIPOT";
111 fMatConstPropNames[kDIFFUSION] = "DIFFUSION";
112 fMatConstPropNames[kSPINFLIP] = "SPINFLIP";
113 fMatConstPropNames[kLOSS] = "LOSS";
114 fMatConstPropNames[kLOSSCS] = "LOSSCS";
115 fMatConstPropNames[kABSCS] = "ABSCS";
116 fMatConstPropNames[kSCATCS] = "SCATCS";
117 fMatConstPropNames[kMR_NBTHETA] = "MR_NBTHETA";
118 fMatConstPropNames[kMR_NBE] = "MR_NBE";
119 fMatConstPropNames[kMR_RRMS] = "MR_RRMS";
120 fMatConstPropNames[kMR_CORRLEN] = "MR_CORRLEN";
121 fMatConstPropNames[kMR_THETAMIN] = "MR_THETAMIN";
122 fMatConstPropNames[kMR_THETAMAX] = "MR_THETAMAX";
123 fMatConstPropNames[kMR_EMIN] = "MR_EMIN";
124 fMatConstPropNames[kMR_EMAX] = "MR_EMAX";
125 fMatConstPropNames[kMR_ANGNOTHETA] = "MR_ANGNOTHETA";
126 fMatConstPropNames[kMR_ANGNOPHI] = "MR_ANGNOPHI";
127 fMatConstPropNames[kMR_ANGCUT] = "MR_ANGCUT";
128 fMatConstPropNames[kSCINTILLATIONTIMECONSTANT1] = "SCINTILLATIONTIMECONSTANT1";
129 fMatConstPropNames[kSCINTILLATIONTIMECONSTANT2] = "SCINTILLATIONTIMECONSTANT2";
130 fMatConstPropNames[kSCINTILLATIONTIMECONSTANT3] = "SCINTILLATIONTIMECONSTANT3";
131 fMatConstPropNames[kSCINTILLATIONRISETIME1] = "SCINTILLATIONRISETIME1";
132 fMatConstPropNames[kSCINTILLATIONRISETIME2] = "SCINTILLATIONRISETIME2";
133 fMatConstPropNames[kSCINTILLATIONRISETIME3] = "SCINTILLATIONRISETIME3";
134 fMatConstPropNames[kSCINTILLATIONYIELD1] = "SCINTILLATIONYIELD1";
135 fMatConstPropNames[kSCINTILLATIONYIELD2] = "SCINTILLATIONYIELD2";
136 fMatConstPropNames[kSCINTILLATIONYIELD3] = "SCINTILLATIONYIELD3";
137 fMatConstPropNames[kPROTONSCINTILLATIONYIELD1] = "PROTONSCINTILLATIONYIELD1";
138 fMatConstPropNames[kPROTONSCINTILLATIONYIELD2] = "PROTONSCINTILLATIONYIELD2";
139 fMatConstPropNames[kPROTONSCINTILLATIONYIELD3] = "PROTONSCINTILLATIONYIELD3";
140 fMatConstPropNames[kDEUTERONSCINTILLATIONYIELD1] = "DEUTERONSCINTILLATIONYIELD1";
141 fMatConstPropNames[kDEUTERONSCINTILLATIONYIELD2] = "DEUTERONSCINTILLATIONYIELD2";
142 fMatConstPropNames[kDEUTERONSCINTILLATIONYIELD3] = "DEUTERONSCINTILLATIONYIELD3";
143 fMatConstPropNames[kTRITONSCINTILLATIONYIELD1] = "TRITONSCINTILLATIONYIELD1";
144 fMatConstPropNames[kTRITONSCINTILLATIONYIELD2] = "TRITONSCINTILLATIONYIELD2";
145 fMatConstPropNames[kTRITONSCINTILLATIONYIELD3] = "TRITONSCINTILLATIONYIELD3";
146 fMatConstPropNames[kALPHASCINTILLATIONYIELD1] = "ALPHASCINTILLATIONYIELD1";
147 fMatConstPropNames[kALPHASCINTILLATIONYIELD2] = "ALPHASCINTILLATIONYIELD2";
148 fMatConstPropNames[kALPHASCINTILLATIONYIELD3] = "ALPHASCINTILLATIONYIELD3";
149 fMatConstPropNames[kIONSCINTILLATIONYIELD1] = "IONSCINTILLATIONYIELD1";
150 fMatConstPropNames[kIONSCINTILLATIONYIELD2] = "IONSCINTILLATIONYIELD2";
151 fMatConstPropNames[kIONSCINTILLATIONYIELD3] = "IONSCINTILLATIONYIELD3";
152 fMatConstPropNames[kELECTRONSCINTILLATIONYIELD1] = "ELECTRONSCINTILLATIONYIELD1";
153 fMatConstPropNames[kELECTRONSCINTILLATIONYIELD2] = "ELECTRONSCINTILLATIONYIELD2";
154 fMatConstPropNames[kELECTRONSCINTILLATIONYIELD3] = "ELECTRONSCINTILLATIONYIELD3";
155 fMatConstPropNames[kCOATEDTHICKNESS] = "COATEDTHICKNESS";
156 fMatConstPropNames[kCOATEDFRUSTRATEDTRANSMISSION] = "COATEDFRUSTRATEDTRANSMISSION";
157 fMatConstPropNames[kPROTONSCINTILLATIONTIMECONSTANT1] = "PROTONSCINTILLATIONTIMECONSTANT1";
158 fMatConstPropNames[kPROTONSCINTILLATIONTIMECONSTANT2] = "PROTONSCINTILLATIONTIMECONSTANT2";
159 fMatConstPropNames[kPROTONSCINTILLATIONTIMECONSTANT3] = "PROTONSCINTILLATIONTIMECONSTANT3";
160 fMatConstPropNames[kDEUTERONSCINTILLATIONTIMECONSTANT1] = "DEUTERONSCINTILLATIONTIMECONSTANT1";
161 fMatConstPropNames[kDEUTERONSCINTILLATIONTIMECONSTANT2] = "DEUTERONSCINTILLATIONTIMECONSTANT2";
162 fMatConstPropNames[kDEUTERONSCINTILLATIONTIMECONSTANT3] = "DEUTERONSCINTILLATIONTIMECONSTANT3";
163 fMatConstPropNames[kTRITONSCINTILLATIONTIMECONSTANT1] = "TRITONSCINTILLATIONTIMECONSTANT1";
164 fMatConstPropNames[kTRITONSCINTILLATIONTIMECONSTANT2] = "TRITONSCINTILLATIONTIMECONSTANT2";
165 fMatConstPropNames[kTRITONSCINTILLATIONTIMECONSTANT3] = "TRITONSCINTILLATIONTIMECONSTANT3";
166 fMatConstPropNames[kALPHASCINTILLATIONTIMECONSTANT1] = "ALPHASCINTILLATIONTIMECONSTANT1";
167 fMatConstPropNames[kALPHASCINTILLATIONTIMECONSTANT2] = "ALPHASCINTILLATIONTIMECONSTANT2";
168 fMatConstPropNames[kALPHASCINTILLATIONTIMECONSTANT3] = "ALPHASCINTILLATIONTIMECONSTANT3";
169 fMatConstPropNames[kIONSCINTILLATIONTIMECONSTANT1] = "IONSCINTILLATIONTIMECONSTANT1";
170 fMatConstPropNames[kIONSCINTILLATIONTIMECONSTANT2] = "IONSCINTILLATIONTIMECONSTANT2";
171 fMatConstPropNames[kIONSCINTILLATIONTIMECONSTANT3] = "IONSCINTILLATIONTIMECONSTANT3";
172 fMatConstPropNames[kELECTRONSCINTILLATIONTIMECONSTANT1] = "ELECTRONSCINTILLATIONTIMECONSTANT1";
173 fMatConstPropNames[kELECTRONSCINTILLATIONTIMECONSTANT2] = "ELECTRONSCINTILLATIONTIMECONSTANT2";
174 fMatConstPropNames[kELECTRONSCINTILLATIONTIMECONSTANT3] = "ELECTRONSCINTILLATIONTIMECONSTANT3";
175 fMCP.assign(kNumberOfConstPropertyIndex, {0., false});
176 }
177
178 G4MaterialPropertiesTable::~G4MaterialPropertiesTable()
179 {
180 for (auto prop : fMP) {
181 delete (prop);
182 }
183 }
184
185 G4int G4MaterialPropertiesTable::GetConstPropertyIndex(const G4String& key) const
186 {
187 // Returns the constant material property index corresponding to a key
188
189 std::size_t index = std::distance(fMatConstPropNames.cbegin(),
190 std::find(fMatConstPropNames.cbegin(), fMatConstPropNames.cend(), key));
191 if (index < fMatConstPropNames.size()) {
192 return (G4int)index;
193 }
194
195 G4ExceptionDescription ed;
196 ed << "Constant Material Property Index for key " << key << " not found.";
197 G4Exception("G4MaterialPropertiesTable::GetConstPropertyIndex()", "mat200", FatalException, ed);
198 return 0;
199 }
200
201 G4int G4MaterialPropertiesTable::GetPropertyIndex(const G4String& key) const
202 {
203 // Returns the material property index corresponding to a key
204 std::size_t index = std::distance(
205 fMatPropNames.cbegin(), std::find(fMatPropNames.cbegin(), fMatPropNames.cend(), key));
206 if (index < fMatPropNames.size()) {
207 return (G4int)index;
208 }
209 G4ExceptionDescription ed;
210 ed << "Material Property Index for key " << key << " not found.";
211 G4Exception("G4MaterialPropertiesTable::GetPropertyIndex()", "mat201", FatalException, ed);
212 return 0;
213 }
214
215 G4double G4MaterialPropertiesTable::GetConstProperty(const G4int index) const
216 {
217 // Returns the constant material property corresponding to an index
218 // fatal exception if property not found
219
220 if (index < (G4int)fMCP.size() && fMCP[index].second) {
221 return fMCP[index].first;
222 }
223 G4ExceptionDescription ed;
224 ed << "Constant Material Property " << fMatConstPropNames[index] << " not found.";
225 G4Exception("G4MaterialPropertiesTable::GetConstProperty()", "mat202", FatalException, ed);
226 return 0.;
227 }
228
229 G4double G4MaterialPropertiesTable::GetConstProperty(const G4String& key) const
230 {
231 // Returns the constant material property corresponding to a key
232 // fatal exception if property not found
233
234 return GetConstProperty(GetConstPropertyIndex(key));
235 }
236
237 G4double G4MaterialPropertiesTable::GetConstProperty(const char* key) const
238 {
239 return GetConstProperty(GetConstPropertyIndex(G4String(key)));
240 }
241
242 G4bool G4MaterialPropertiesTable::ConstPropertyExists(const G4int index) const
243 {
244 // Returns true if a const property corresponding to 'index' exists
245
246 return index >= 0 && index < (G4int)fMCP.size() && fMCP[index].second;
247 }
248
249 G4bool G4MaterialPropertiesTable::ConstPropertyExists(const G4String& key) const
250 {
251 // Returns true if a const property 'key' exists
252 std::size_t index = std::distance(fMatConstPropNames.cbegin(),
253 std::find(fMatConstPropNames.cbegin(), fMatConstPropNames.cend(), key));
254 if (index < fMatConstPropNames.size()) { // index is type std::size_t so >= 0
255 return ConstPropertyExists((G4int)index);
256 }
257 return false;
258 }
259
260 G4bool G4MaterialPropertiesTable::ConstPropertyExists(const char* key) const
261 {
262 std::size_t index = std::distance(fMatConstPropNames.cbegin(),
263 std::find(fMatConstPropNames.cbegin(), fMatConstPropNames.cend(), key));
264 if (index < fMatConstPropNames.size()) { // index is type std::size_t so >= 0
265 return ConstPropertyExists((G4int)index);
266 }
267 return false;
268 }
269
270 G4MaterialPropertyVector* G4MaterialPropertiesTable::GetProperty(const G4String& key) const
271 {
272 // Returns a Material Property Vector corresponding to a key
273 if (std::find(fMatPropNames.cbegin(), fMatPropNames.cend(), key) != fMatPropNames.cend()) {
274 const G4int index = GetPropertyIndex(G4String(key));
275 return GetProperty(index);
276 }
277 return nullptr;
278 }
279
280 G4MaterialPropertyVector* G4MaterialPropertiesTable::GetProperty(const char* key) const
281 {
282 if (std::find(fMatPropNames.cbegin(), fMatPropNames.cend(), key) != fMatPropNames.cend()) {
283 const G4int index = GetPropertyIndex(G4String(key));
284 return GetProperty(index);
285 }
286 return nullptr;
287 }
288
289 G4MaterialPropertyVector* G4MaterialPropertiesTable::GetProperty(const G4int index) const
290 {
291 // Returns a Material Property Vector corresponding to an index
292 // returns nullptr if the property has not been defined by user
293 if (index >= 0 && index < (G4int)fMP.size()) {
294 return fMP[index];
295 }
296 return nullptr;
297 }
298
299 G4MaterialPropertyVector* G4MaterialPropertiesTable::AddProperty(const G4String& key,
300 const std::vector<G4double>& photonEnergies, const std::vector<G4double>& propertyValues,
301 G4bool createNewKey, G4bool spline)
302 {
303 if (photonEnergies.size() != propertyValues.size()) {
304 G4ExceptionDescription ed;
305 ed << "AddProperty error. Number of property values must be equal to the number of\n"
306 << "energy values. Property name: " << key;
307 G4Exception("G4MaterialPropertiesTable::AddProperty()", "mat204", FatalException, ed);
308 }
309
310 if (photonEnergies.size() == 1) {
311 G4ExceptionDescription ed;
312 ed << "AddProperty warning. A material property vector must have more than one value.\n"
313 << "Unless you will later add an entry, this is an error.\n"
314 << "Property name: " << key;
315 G4Exception("G4MaterialPropertiesTable::AddProperty()", "mat218", JustWarning, ed);
316 }
317
318 // G4PhysicsVector assumes energies are in increasing order
319 for (std::size_t i = 0; i < photonEnergies.size() - 1; ++i) {
320 if (photonEnergies.at(i + 1) < photonEnergies.at(i)) {
321 G4ExceptionDescription ed;
322 ed << "Energies in material property vector must be in increasing "
323 << "order. Key: " << key << " Energy: " << photonEnergies.at(i + 1);
324 G4Exception("G4MaterialPropertiesTable::AddProperty()", "mat215", FatalException, ed);
325 }
326 }
327
328 // if the key doesn't exist, add it if requested
329 if (std::find(fMatPropNames.cbegin(), fMatPropNames.cend(), key) == fMatPropNames.cend()) {
330 if (createNewKey) {
331 fMatPropNames.push_back(key);
332 fMP.push_back(nullptr);
333 }
334 else {
335 G4ExceptionDescription ed;
336 ed << "Attempting to create a new material property vector " << key << " without setting\n"
337 << "createNewKey parameter of AddProperty to true.";
338 G4Exception("G4MaterialPropertiesTable::AddProperty()", "mat205", FatalException, ed);
339 }
340 }
341
342 auto* mpv = new G4MaterialPropertyVector(photonEnergies, propertyValues, spline);
343 mpv->SetVerboseLevel(1);
344 if (spline) {
345 mpv->FillSecondDerivatives();
346 }
347 G4int index = GetPropertyIndex(key);
348 fMP[index] = mpv;
349
350 // if key is RINDEX, we calculate GROUPVEL -
351 // contribution from Tao Lin (IHEP, the JUNO experiment)
352 if (key == "RINDEX") {
353 CalculateGROUPVEL();
354 }
355
356 return mpv;
357 }
358
359 G4MaterialPropertyVector* G4MaterialPropertiesTable::AddProperty(const char* key,
360 G4double* photonEnergies, G4double* propertyValues, G4int numEntries, G4bool createNewKey,
361 G4bool spline)
362 {
363 // Provides a way of adding a property to the Material Properties
364 // Table given a pair of arrays and a key
365 G4String k(key);
366
367 std::vector<G4double> energies(photonEnergies, photonEnergies + numEntries);
368 std::vector<G4double> values(propertyValues, propertyValues + numEntries);
369 return AddProperty(k, energies, values, createNewKey, spline);
370 }
371
372 void G4MaterialPropertiesTable::AddProperty(
373 const G4String& key, G4MaterialPropertyVector* mpv, G4bool createNewKey)
374 {
375 // Provides a way of adding a property to the Material Properties
376 // Table given an G4MaterialPropertyVector Reference and a key
377
378 // G4PhysicsVector assumes energies are in increasing order
379 if (mpv->GetVectorLength() > 1) {
380 for (std::size_t i = 0; i < mpv->GetVectorLength() - 1; ++i) {
381 if (mpv->Energy(i + 1) < mpv->Energy(i)) {
382 G4ExceptionDescription ed;
383 ed << "Energies in material property vector must be in increasing "
384 << "order. Key: " << key << " Energy: " << mpv->Energy(i + 1);
385 G4Exception("G4MaterialPropertiesTable::AddProperty()", "mat216", FatalException, ed);
386 }
387 }
388 }
389
390 if (mpv->GetVectorLength() <= 1) {
391 G4ExceptionDescription ed;
392 ed << "AddProperty warning. A material property vector must have more than one value.\n"
393 << "Unless you will later add an entry, this is an error.\n"
394 << "Property name: " << key;
395 G4Exception("G4MaterialPropertiesTable::AddProperty()", "mat219", JustWarning, ed);
396 }
397
398 // if the key doesn't exist, add it
399 if (std::find(fMatPropNames.cbegin(), fMatPropNames.cend(), key) == fMatPropNames.cend()) {
400 if (createNewKey) {
401 fMatPropNames.push_back(key);
402 fMP.push_back(nullptr);
403 }
404 else {
405 G4ExceptionDescription ed;
406 ed << "Attempting to create a new material property key " << key << " without setting\n"
407 << "createNewKey parameter of AddProperty to true.";
408 G4Exception("G4MaterialPropertiesTable::AddProperty()", "mat206", FatalException, ed);
409 }
410 }
411 G4int index = GetPropertyIndex(key);
412 fMP[index] = mpv;
413
414 // if key is RINDEX, we calculate GROUPVEL -
415 // contribution from Tao Lin (IHEP, the JUNO experiment)
416 if (key == "RINDEX") {
417 CalculateGROUPVEL();
418 }
419 }
420
421 void G4MaterialPropertiesTable::AddProperty(
422 const char* key, G4MaterialPropertyVector* mpv, G4bool createNewKey)
423 {
424 AddProperty(G4String(key), mpv, createNewKey);
425 }
426
427 void G4MaterialPropertiesTable::AddProperty(const G4String& key, const G4String& mat)
428 {
429 // load a material property vector defined in Geant4 source
430 G4MaterialPropertyVector* v = G4OpticalMaterialProperties::GetProperty(key, mat);
431 AddProperty(key, v);
432 }
433
434 void G4MaterialPropertiesTable::AddConstProperty(
435 const G4String& key, G4double propertyValue, G4bool createNewKey)
436 {
437 // Provides a way of adding a constant property to the Material Properties
438 // Table given a key
439 if (std::find(fMatConstPropNames.cbegin(), fMatConstPropNames.cend(), key) ==
440 fMatConstPropNames.cend())
441 {
442 if (createNewKey) {
443 fMatConstPropNames.push_back(key);
444 fMCP.emplace_back(0., true);
445 }
446 else {
447 G4ExceptionDescription ed;
448 ed << "Attempting to create a new material constant property key " << key
449 << " without setting\n"
450 << "createNewKey parameter of AddProperty to true.";
451 G4Exception("G4MaterialPropertiesTable::AddProperty()", "mat207", FatalException, ed);
452 }
453 }
454 G4int index = GetConstPropertyIndex(key);
455
456 fMCP[index] = std::pair<G4double, G4bool>{propertyValue, true};
457 }
458
459 void G4MaterialPropertiesTable::AddConstProperty(
460 const char* key, G4double propertyValue, G4bool createNewKey)
461 {
462 // Provides a way of adding a constant property to the Material Properties
463 // Table given a key
464 AddConstProperty(G4String(key), propertyValue, createNewKey);
465 }
466
467 void G4MaterialPropertiesTable::RemoveConstProperty(const G4String& key)
468 {
469 G4int index = GetConstPropertyIndex(key);
470 if (index < (G4int)fMCP.size()) {
471 fMCP[index] = std::pair<G4double, G4bool>{0., false};
472 }
473 }
474
475 void G4MaterialPropertiesTable::RemoveConstProperty(const char* key)
476 {
477 RemoveConstProperty(G4String(key));
478 }
479
480 void G4MaterialPropertiesTable::RemoveProperty(const G4String& key)
481 {
482 G4int index = GetPropertyIndex(key);
483 delete fMP[index];
484 fMP[index] = nullptr;
485 }
486
487 void G4MaterialPropertiesTable::RemoveProperty(const char* key) { RemoveProperty(G4String(key)); }
488
489 void G4MaterialPropertiesTable::AddEntry(
490 const G4String& key, G4double aPhotonEnergy, G4double aPropertyValue)
491 {
492 // Allows to add an entry pair directly to the Material Property Vector
493 // given a key.
494 if (std::find(fMatPropNames.cbegin(), fMatPropNames.cend(), key) == fMatPropNames.cend()) {
495 G4ExceptionDescription ed;
496 ed << "Material Property Vector " << key << " not found.";
497 G4Exception("G4MaterialPropertiesTable::AddEntry()", "mat214", FatalException, ed);
498 }
499 G4int index = GetPropertyIndex(key);
500
501 G4MaterialPropertyVector* targetVector = fMP[index];
502 if (targetVector != nullptr) {
503 // do not allow duplicate energies
504 for (std::size_t i = 0; i < targetVector->GetVectorLength(); ++i) {
505 if (aPhotonEnergy == targetVector->Energy(i)) {
506 G4ExceptionDescription ed;
507 ed << "Energy values in material property vector must be unique. "
508 << "Key: " << key;
509 G4Exception("G4MaterialPropertiesTable::AddEntry()", "mat217", FatalException, ed);
510 }
511 }
512
513 targetVector->InsertValues(aPhotonEnergy, aPropertyValue);
514 }
515 else {
516 G4ExceptionDescription ed;
517 ed << "Material Property Vector " << key << " not found.";
518 G4Exception("G4MaterialPropertiesTable::AddEntry()", "mat208", FatalException, ed);
519 }
520 if (key == "RINDEX") {
521 CalculateGROUPVEL();
522 }
523 }
524
525 void G4MaterialPropertiesTable::AddEntry(
526 const char* key, G4double aPhotonEnergy, G4double aPropertyValue)
527 {
528 AddEntry(G4String(key), aPhotonEnergy, aPropertyValue);
529 }
530
531 void G4MaterialPropertiesTable::DumpTable() const
532 {
533 // material properties
534 G4int j = 0;
535 for (const auto& prop : fMP) {
536 if (prop != nullptr) {
537 G4cout << j << ": " << fMatPropNames[j] << G4endl;
538 prop->DumpValues();
539 }
540 ++j;
541 }
542 // material constant properties
543 j = 0;
544 for (const auto& cprop : fMCP) {
545 if (cprop.second) {
546 G4cout << j << ": " << fMatConstPropNames[j] << " " << cprop.first << G4endl;
547 }
548 ++j;
549 }
550 }
551
552 G4MaterialPropertyVector* G4MaterialPropertiesTable::CalculateGROUPVEL()
553 {
554 #ifdef G4MULTITHREADED
555 G4AutoLock mptm(&materialPropertyTableMutex);
556 #endif
557
558 // check if "GROUPVEL" already exists. If so, remove it.
559 if (fMP[kGROUPVEL] != nullptr) {
560 this->RemoveProperty("GROUPVEL");
561 }
562
563 // fetch RINDEX data, give up if unavailable
564 G4MaterialPropertyVector* rindex = this->GetProperty(kRINDEX);
565 if (rindex == nullptr) {
566 return nullptr;
567 }
568
569 // RINDEX exists but has no entries, give up
570 if (rindex->GetVectorLength() == 0) {
571 return nullptr;
572 }
573
574 // add GROUPVEL vector
575 auto* groupvel = new G4MaterialPropertyVector();
576 groupvel->SetVerboseLevel(1);
577
578 // fill GROUPVEL vector using RINDEX values
579 // rindex built-in "iterator" was advanced to first entry above
580 G4double E0 = rindex->Energy(0);
581 G4double n0 = (*rindex)[0];
582
583 if (E0 <= 0.) {
584 G4Exception("G4MaterialPropertiesTable::CalculateGROUPVEL()", "mat211", FatalException,
585 "Optical Photon Energy <= 0");
586 }
587
588 if (rindex->GetVectorLength() >= 2) {
589 // good, we have at least two entries in RINDEX
590 // get next energy/value pair
591
592 G4double E1 = rindex->Energy(1);
593 G4double n1 = (*rindex)[1];
594
595 if (E1 <= 0.) {
596 G4Exception("G4MaterialPropertiesTable::CalculateGROUPVEL()", "mat212", FatalException,
597 "Optical Photon Energy <= 0");
598 }
599
600 G4double vg;
601
602 // add entry at first photon energy
603 vg = c_light / (n0 + (n1 - n0) / G4Log(E1 / E0));
604
605 // allow only for 'normal dispersion' -> dn/d(logE) > 0
606 if ((vg < 0) || (vg > c_light / n0)) {
607 vg = c_light / n0;
608 }
609
610 groupvel->InsertValues(E0, vg);
611
612 // add entries at midpoints between remaining photon energies
613 for (std::size_t i = 2; i < rindex->GetVectorLength(); ++i) {
614 vg = c_light / (0.5 * (n0 + n1) + (n1 - n0) / G4Log(E1 / E0));
615
616 // allow only for 'normal dispersion' -> dn/d(logE) > 0
617 if ((vg < 0) || (vg > c_light / (0.5 * (n0 + n1)))) {
618 vg = c_light / (0.5 * (n0 + n1));
619 }
620 groupvel->InsertValues(0.5 * (E0 + E1), vg);
621
622 // get next energy/value pair, or exit loop
623 E0 = E1;
624 n0 = n1;
625 E1 = rindex->Energy(i);
626 n1 = (*rindex)[i];
627
628 if (E1 <= 0.) {
629 G4Exception("G4MaterialPropertiesTable::CalculateGROUPVEL()", "mat213", FatalException,
630 "Optical Photon Energy <= 0");
631 }
632 }
633
634 // add entry at last photon energy
635 vg = c_light / (n1 + (n1 - n0) / G4Log(E1 / E0));
636
637 // allow only for 'normal dispersion' -> dn/d(logE) > 0
638 if ((vg < 0) || (vg > c_light / n1)) {
639 vg = c_light / n1;
640 }
641 groupvel->InsertValues(E1, vg);
642 }
643 else // only one entry in RINDEX -- weird!
644 {
645 groupvel->InsertValues(E0, c_light / n0);
646 }
647
648 this->AddProperty("GROUPVEL", groupvel);
649
650 return groupvel;
651 }
652