BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
BDSMaterials.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2022.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "BDSDebug.hh"
20#include "BDSException.hh"
21#include "BDSMaterials.hh"
22#include "BDSParser.hh"
23#include "BDSUtilities.hh"
24#include "BDSWarning.hh"
25
26#include "G4MaterialTable.hh"
27#include "G4String.hh"
28#include "G4NistManager.hh"
29#include "G4Version.hh"
30
31#include <iomanip>
32#include <list>
33#include <map>
34#include <vector>
35
36BDSMaterials* BDSMaterials::instance = nullptr;
37
38BDSMaterials* BDSMaterials::Instance()
39{
40 if (!instance)
41 {instance = new BDSMaterials();}
42 return instance;
43}
44
45BDSMaterials::BDSMaterials()
46{
47#ifdef BDSDEBUG
48 G4cout << "BDSMaterials: G4 predefined units: " << G4endl;
49 G4cout << "g= " << CLHEP::g << G4endl;
50 G4cout << "m= " << CLHEP::m << G4endl;
51 G4cout << "mole= " << CLHEP::mole << G4endl;
52 G4cout << "kelvin= " << CLHEP::kelvin << G4endl;
53#endif
54
55 DefineMetals();
56 DefineSuperconductors();
57 DefineNonMetalSolids();
58 DefineScintillators();
59 DefineLHCComponents();
60 DefineLiquids();
61 DefineGases();
62 DefinePlasmas();
63 DefineVacuums();
64}
65
66void BDSMaterials::DefineMetals()
67{
68 //default Geant4 temperature = 293.15 K (NTP_Temperature)
69 //default Geant4 pressure = 1atm
70#if G4VERSION_NUMBER < 1011
71 G4double NTP_Temperature = 293.15;
72#endif
73
74 // solid materials
75 // metals
76 // standard single element metals
77 AddExistingMaterialAlias("Al", "aluminium");
78 AddExistingMaterialAlias("Be", "beryllium");
79 AddExistingMaterialAlias("C", "carbon");
80 AddExistingMaterialAlias("Cr", "chromium");
81 AddExistingMaterialAlias("Cu", "copper");
82 AddExistingMaterialAlias("Fe", "iron");
83 AddExistingMaterialAlias("Pb", "lead");
84 AddExistingMaterialAlias("Mg", "magnesium");
85 AddExistingMaterialAlias("Ni", "nickel");
86 AddExistingMaterialAlias("Si", "silicon");
87 AddExistingMaterialAlias("Ti", "titanium");
88 AddExistingMaterialAlias("W", "tungsten");
89 AddExistingMaterialAlias("U", "uranium");
90 AddExistingMaterialAlias("V", "vanadium");
91 AddExistingMaterialAlias("Zn", "zinc");
92
93 // special forms
94 std::list<int> singleElement = {1};
95 AddMaterial("graphite", 2.265, kStateSolid, NTP_Temperature, 1, {"C"}, singleElement);
96 AddMaterial("graphitefoam", 0.61, kStateSolid, NTP_Temperature, 1, {"C"}, singleElement);
97 AddMaterial("solidhydrogen",8.96, kStateSolid, NTP_Temperature, 1, {"H"}, singleElement);
98 AddMaterial("solidnitrogen",8.96, kStateSolid, NTP_Temperature, 1, {"N"}, singleElement);
99 AddMaterial("solidoxygen", 8.96, kStateSolid, NTP_Temperature, 1, {"O"}, singleElement);
100 AddMaterial("weightiron", 7.87, kStateSolid, NTP_Temperature, 1, {"Fe"}, singleElement);
101
102 // composites and alloys
103 AddMaterial("leadtungstate" ,
104 8.27,
105 kStateSolid, NTP_Temperature, 1,
106 {"Pb","W" ,"O" },
107 std::list<int>{1, 1, 4});
108
109 AddMaterial("smco",
110 8.4,
111 kStateSolid, 300, 1,
112 {"Sm","Co"},
113 std::list<double>{0.338, 0.662});
114
115 // Titanium alloy (BDS collimator material SLAC-TN-03-069 p25) deacon 15 Jun 2007
116 AddMaterial("titaniumalloy",
117 4.48,
118 kStateSolid, 300, 1,
119 {"V" ,"Al","Ti"},
120 std::list<double>{0.025, 0.03, 0.945});
121
122 // Carbon Steel (shell of cryomodule). LDeacon 21 Feb 2006
123 AddMaterial("carbonsteel",
124 7.87,
125 kStateSolid, 100, 1,
126 {"C","Mn","P","S","Fe"},
127 std::list<double>{0.0017, 0.0045, 0.0004, 0.0005, 0.9929});
128
129 // Copper Alloy 17410 (C17410) "Beryllium Copper"
130 AddMaterial("berylliumcopper",
131 8.8,
132 kStateSolid, 300, 1,
133 {"Cu", "Be", "Co", "Al", "Fe", "Ni"},
134 std::list<double>{0.991, 0.0031, 0.00500, 0.0004, 0.0003, 0.0002});
135
136 // Stainless Steel 316L
137 AddMaterial("stainlesssteel",
138 8.0,
139 kStateSolid, 295, 1,
140 {"C","Mn","Si","P","S","Cr","Mo","Ni","N","Fe"},
141 std::list<double>{0.0003, 0.02, 0.0075, 0.00045, 0.0003, 0.17, 0.025, 0.12, 0.001, 0.65545});
142
143 // Stainless Steel AISI code 304L (low-carbon) @ 300K
144 AddMaterial("stainless_steel_304L",
145 8.02,
146 kStateSolid, 300, 1,
147 {"Fe", "Cr", "Ni", "Mn", "Si", "P", "S", "C"},
148 std::list<double>{0.67145, 0.185, 0.1125, 0.02, 0.01, 0.00045, 0.0003, 0.0003});
149
150 // Stainless Steel AISI code 304L (low-carbon) @ 87K
151 AddMaterial("stainless_steel_304L_87K",
152 8.02,
153 kStateSolid, 87, 1,
154 {"Fe", "Cr", "Ni", "Mn", "Si", "P", "S", "C"},
155 std::list<double>{0.67145, 0.185, 0.1125, 0.02, 0.01, 0.00045, 0.0003, 0.0003});
156
157 // Stainless Steel AISI code 304L (low-carbon) @ 2K
158 AddMaterial("stainless_steel_304L_2K",
159 8.02,
160 kStateSolid, 2, 1,
161 {"Fe", "Cr", "Ni", "Mn", "Si", "P", "S", "C"},
162 std::list<double>{0.67145, 0.185, 0.1125, 0.02, 0.01, 0.00045, 0.0003, 0.0003});
163
164 // Stainless Steel AISI code 316LN
165 // (Type 316, low carbon, nitrogen-enhanced) @ 300K
166 AddMaterial("stainless_steel_316LN",
167 8.03,
168 kStateSolid, 300, 1,
169 {"Fe", "Cr", "Ni", "Mo", "Mn", "Si", "Ti", "N", "Nb",
170 "Cu", "Co", "P", "C", "S", "Ta", "B"},
171 std::list<double>{0.65093, 0.1700, 0.12000, 0.02500, 0.0200, 0.00750,
172 0.00150, 0.0014, 0.00100, 0.00100, 0.0005, 0.00045,
173 0.00030, 0.0003, 0.00010, 0.00002});
174
175 // Stainless Steel AISI code 316LN
176 // (Type 316, low-carbon nitrogen-enhanced) @ 87K
177 AddMaterial("stainless_steel_316LN_87K",
178 8.03,
179 kStateSolid, 87, 1,
180 {"Fe", "Cr", "Ni", "Mo", "Mn", "Si", "Ti", "N",
181 "Nb", "Cu", "Co", "P", "C", "S", "Ta", "B"},
182 std::list<double>{0.65093, 0.1700, 0.12000, 0.02500, 0.0200,
183 0.00750, 0.00150, 0.0014, 0.00100, 0.00100,
184 0.0005, 0.00045, 0.00030, 0.0003, 0.00010,
185 0.00002});
186
187 // Stainless Steel AISI code 316LN
188 // (Type 316, low-carbon nitrogen-enhanced) @ 2K
189 AddMaterial("stainless_steel_316LN_2K",
190 8.03,
191 kStateSolid, 2, 1,
192 {"Fe", "Cr", "Ni", "Mo", "Mn", "Si", "Ti", "N",
193 "Nb", "Cu", "Co", "P", "C", "S", "Ta", "B"},
194 std::list<double>{0.65093, 0.1700, 0.12000, 0.02500, 0.0200,
195 0.00750, 0.00150, 0.0014, 0.00100, 0.00100,
196 0.0005, 0.00045, 0.00030, 0.0003, 0.00010,
197 0.00002});
198
199
200 // Mild Steel
201 AddMaterial("mild_steel", 8.000, kStateSolid, 295, 1,
202 {"C", "Mn", "Si", "Fe"},
203 std::list<double>{0.002, 0.005, 0.0015, 0.99150});
204
205 // Pure tungsten is not typically used, but instead part of "heavy
206 // alloy." I think "heavy" in the sense that the tungsten makes up almost
207 // all of the composition, and tungsten is a very dense metal.
208 AddMaterial("tungsten_heavy_alloy",
209 18.5,
210 kStateSolid, 87, 1,
211 {"W", "Ni", "Fe"},
212 std::list<double>{0.97, 0.02, 0.01});
213
214 // For HL-LHC Collimation - Inermet170
215 AddMaterial("inermet170",
216 17.0,
217 kStateSolid, 300, 1,
218 {"W", "Ni", "Cu"},
219 std::list<double>{0.90, 0.05, 0.05});
220
221 // For HL-LHC Collimation - Inermet176
222 AddMaterial("inermet176",
223 17.6,
224 kStateSolid, 300, 1,
225 {"W", "Ni", "Cu"},
226 std::list<double>{0.925, 0.0375, 0.0375});
227
228 // For HL-LHC Collimation - Inermet180
229 AddMaterial("inermet180",
230 18.0,
231 kStateSolid, 300, 1,
232 {"W", "Ni", "Cu"},
233 std::list<double>{0.95, 0.025, 0.025});
234}
235
236void BDSMaterials::DefineSuperconductors()
237{
238 // niobium at 2K
239 AddMaterial("niobium_2k", 8.57 , kStateSolid, 2, 1, {"Nb"}, std::list<int>{1});
240 AddExistingMaterialAlias("niobium_2k", "nb_2k"); // alias
241 // niobium titanium at 4K
242 AddMaterial("nbti_4k", 5.6 , kStateSolid, 4, 1, {"Nb","Ti"}, std::list<int>{1,1});
243}
244
245void BDSMaterials::DefineNonMetalSolids()
246{
247 // boron Nickel (absorber)
248 AddMaterial("bn5000",
249 1.925,
250 kStateSolid, 300, 1,
251 {"B","Ni","O","Ca","Si"},
252 std::list<double>{0.383249242, 0.472071387, 0.0366276887, 0.0228923054, 0.0851593762});
253
254 // calcium carbonate (calcite)
255 AddMaterial("calciumCarbonate",
256 2.711,
257 kStateSolid, 300, 1,
258 {"Ca","O","C"},
259 std::list<int>{1,3,1});
260
261 // clay
262 AddMaterial("clay",
263 1.746,
264 kStateSolid, 300, 1,
265 {"Al","O","Si","H"},
266 std::list<int>{1,9,2,4});
267
268 AddMaterial("concrete",
269 2.3,
270 kStateSolid, 300, 1, {"Si","O","H","Ca","Al","Fe"},
271 std::list<double>{0.227915, 0.60541, 0.09972, 0.04986, 0.014245, 0.00285});
272
273 // we make the material manually here so we have a pointer to it. note density units must
274 // be used here.
275 G4Material* tmpMaterial = new G4Material("fusedsilica",
276 1.032*CLHEP::g/CLHEP::cm3, // density
277 2,
278 kStateSolid);
279 tmpMaterial->AddElement(GetElement("O") , 2);
280 tmpMaterial->AddElement(GetElement("Si"), 1);
281 const G4int FusedSilica_NUMENTRIES = 3; //Number of entries in the material properties table
282 G4double FusedSilica_RIND[FusedSilica_NUMENTRIES]={1.49,1.49,1.49};
283 G4double FusedSilica_AbsLength[FusedSilica_NUMENTRIES]={420.*CLHEP::cm,420.*CLHEP::cm,420.*CLHEP::cm};
284 G4double FusedSilica_Energy[FusedSilica_NUMENTRIES] = {2.0*CLHEP::eV,7.0*CLHEP::eV,7.14*CLHEP::eV};
285 G4MaterialPropertiesTable* fsMaterialPropertiesTable = CreatePropertiesTable();
286 fsMaterialPropertiesTable->AddProperty("ABSLENGTH",FusedSilica_Energy,FusedSilica_AbsLength,FusedSilica_NUMENTRIES);
287 fsMaterialPropertiesTable->AddProperty("RINDEX",FusedSilica_Energy,FusedSilica_RIND,FusedSilica_NUMENTRIES);
288 tmpMaterial->SetMaterialPropertiesTable(fsMaterialPropertiesTable);
289 AddMaterial(tmpMaterial, "fusedsilica");
290
291 // perspex
292 AddMaterial("perspex",
293 1.18,
294 kStateSolid,
295 300,
296 1,
297 {"C","O","H"},
298 std::list<double>{0.59984,0.31961,0.08055});
299
300 // invar - Temperature 2 kelvin. LDeacon 6th Feburary 2006
301 AddMaterial("invar" ,
302 8.1 ,
303 kStateSolid, 2, 1,
304 {"Ni","Fe"},
305 std::list<double>{0.35,0.65});
306
307 // kapton polyimide film
308 AddMaterial("kapton",
309 1.42,
310 kStateSolid, 295, 1,
311 {"H","C","N","O"},
312 std::list<double>{0.026362,0.691133,0.073270,0.209235});
313
314 // n-bk7
315 tmpMaterial = new G4Material("n-bk7", 1.032*CLHEP::g/CLHEP::cm3, 2, kStateSolid);
316 tmpMaterial->AddElement(GetElement("O") , 2);
317 tmpMaterial->AddElement(GetElement("Si"), 1);
318 const G4int N_Bk7_NUMENTRIES = 3; //Number of entries in the material properties table
319 G4double N_Bk7_RIND[N_Bk7_NUMENTRIES]={1.51680,1.51680,1.51680};
320 G4double N_Bk7_AbsLength[N_Bk7_NUMENTRIES]={420.*CLHEP::cm,420.*CLHEP::cm,420.*CLHEP::cm};
321 G4double N_Bk7_Energy[N_Bk7_NUMENTRIES] = {2.0*CLHEP::eV,7.0*CLHEP::eV,7.14*CLHEP::eV};
322 G4MaterialPropertiesTable* nbk7MaterialPropertiesTable = CreatePropertiesTable();
323 nbk7MaterialPropertiesTable->AddProperty("ABSLENGTH",N_Bk7_Energy,N_Bk7_AbsLength,N_Bk7_NUMENTRIES);
324 nbk7MaterialPropertiesTable->AddProperty("RINDEX",N_Bk7_Energy,N_Bk7_RIND,N_Bk7_NUMENTRIES);
325 tmpMaterial->SetMaterialPropertiesTable(nbk7MaterialPropertiesTable);
326 AddMaterial(tmpMaterial, "n-bk7");
327
328 // quartz
329 AddMaterial("quartz", 2.655, kStateSolid, 300, 1, {"Si","O"}, std::list<int>{1,2});
330
331 // marl
332 AddMaterial("marl",
333 1.474,
334 kStateSolid, 300, 1,
335 {"clay","calciumCarbonate"},
336 std::list<double>{0.5,0.5});
337
338 // clayousMarl
339 AddMaterial("clayousMarl",
340 1.555,
341 kStateSolid, 300, 1,
342 {"clay","calciumCarbonate"},
343 std::list<double>{0.65,0.35});
344
345 // limousMarl
346 AddMaterial("limousMarl",
347 1.392,
348 kStateSolid, 300, 1,
349 {"clay", "calciumCarbonate"},
350 std::list<double>{0.35,0.65});
351
352 // "standard" soil (dry)
353 AddMaterial("soil",
354 1.9,
355 kStateSolid, 300, 1,
356 {"Si","O","H","Al"},
357 std::list<double>{0.33377483443708611, 0.57218543046357617, 0.022516556291390728, 0.071523178807947022});
358
359 // epoxy Resin components
360 // the main component of epoxy resin commonly used to insulate magnet coils.
361 AddMaterial("aralditef", 1.175, kStateSolid, 300, 1, {"C","H","O"},std::list<int>{12,18,4});
362
363 // a hardener for the epoxy resin
364 AddMaterial("hy906", 1.225, kStateSolid, 300, 1, {"C","H","O"},std::list<int>{10,5,3});
365
366 // an accelerator for epoxy resin
367 AddMaterial("dy061", 1.025, kStateSolid, 300, 1, {"C","H","O","N"},std::list<int>{15,25,1,3});
368
369 // material type 3 from CERN 81-05, "The Selection and Properties of Epoxide Resins
370 // Used for the Insulation of Magnet Systems in Radiation Environments".
371 AddMaterial("epoxyresin3",
372 1.20,
373 kStateSolid, 300, 1,
374 {"aralditef","hy906","dy061"},
375 std::list<double>{0.497512,0.497512,0.004976});
376
377 // cellulose
378 tmpMaterial = G4NistManager::Instance()->FindOrBuildMaterial("G4_CELLULOSE_CELLOPHANE",true,true);
379 const G4int Cellulose_NUMENTRIES = 3; //Number of entries in the material properties table
380 G4double Cellulose_RIND[Cellulose_NUMENTRIES] = {1.532,1.532,1.532};//Assume constant refractive index.
381 G4double Cellulose_Energy[Cellulose_NUMENTRIES] = {2.0*CLHEP::eV,7.0*CLHEP::eV,7.14*CLHEP::eV}; //The energies.
382 G4MaterialPropertiesTable* celluloseMaterialPropertiesTable = CreatePropertiesTable();
383 celluloseMaterialPropertiesTable->AddProperty("RINDEX",Cellulose_Energy, Cellulose_RIND, Cellulose_NUMENTRIES);
384 tmpMaterial->SetMaterialPropertiesTable(celluloseMaterialPropertiesTable);
385 AddMaterial(tmpMaterial, "cellulose");
386
387 // polyurethane
388 AddMaterial("polyurethane",
389 1.05,
390 kStateSolid, NTP_Temperature, 1,
391 {"C","H","N","O"},
392 std::list<int>{6,10,2,4});
393
394 // RCH 1000 - Ultra high molecular weight polyethylene [PE-UHMW]
395 // at 4K for LHC dipoles
396 AddMaterial("rch1000_4k",
397 0.925,
398 kStateSolid, 4, 1,
399 {"C", "H"},
400 std::list<int>{2,4});
401}
402
403void BDSMaterials::DefineScintillators()
404{
405 //YAG
406 G4Material* tmpMaterial = new G4Material("yag", 4.56*CLHEP::g/CLHEP::cm3, 3);
407 tmpMaterial->AddElement(GetElement("Y") , 3);
408 tmpMaterial->AddElement(GetElement("Al"), 5);
409 tmpMaterial->AddElement(GetElement("O") , 12);
410 G4double birks = 0.08*CLHEP::mm/CLHEP::MeV;
411 tmpMaterial->GetIonisation()->SetBirksConstant(birks);
412 G4MaterialPropertiesTable* mpt_YAG = CreatePropertiesTable();
413 // const int nEntries3=60;
414 const G4int nEntries = 9;
415 G4double PhotonEnergyYAG[nEntries];
416 G4double dNEntries2=(G4double)nEntries;
417 G4double energyMin=1.91*CLHEP::eV;
418 G4double energyMax=2.76*CLHEP::eV;
419 G4double deltaEnergy=(energyMax-energyMin)/(dNEntries2-1.0);
420 G4double energy=energyMin;
421 for (G4int i=0; i<nEntries; energy += deltaEnergy, i++)
422 {PhotonEnergyYAG[i]=energy;}
423
424 G4double RefractiveIndexYAG[nEntries] = //Approximately correct, but check for different wavelengths
425 { 1.82, 1.82, 1.82, 1.82, 1.82, 1.82, 1.82,
426 1.82, 1.82 };
427
428 mpt_YAG->AddProperty("RINDEX",PhotonEnergyYAG, RefractiveIndexYAG, nEntries);
429#if G4VERSION_NUMBER < 1079
430 G4double scintFastYAG[nEntries] = //Approximately correct
431 { 0, 0.25, 2.0, 14.0, 13.0, 7.0, 4.0, 2.0, 0.0 };
432 // All of these parameters are deprecated in V11 onwards and will cause a crash
433 mpt_YAG->AddProperty("FASTCOMPONENT",PhotonEnergyYAG, scintFastYAG, nEntries)->SetSpline(true);
434 mpt_YAG->AddConstProperty("FASTTIMECONSTANT",70.*CLHEP::ns); //Approximately correct
435 mpt_YAG->AddConstProperty("YIELDRATIO",1.0);
436#endif
437 mpt_YAG->AddConstProperty("SCINTILLATIONYIELD",8000./CLHEP::MeV); //Approximately correct
438 mpt_YAG->AddConstProperty("RESOLUTIONSCALE",2.0); //Check this
439 tmpMaterial->SetMaterialPropertiesTable(mpt_YAG);
440 AddMaterial(tmpMaterial, "yag");
441
442 //UPS-923A - see http://www.amcrys-h.com/
443 //Define the material properties (copy from NIST table of materials).
444 G4NistManager* nistManager = G4NistManager::Instance();
445 G4Material* polystyrene = nistManager->FindOrBuildMaterial("G4_POLYSTYRENE",true,true);
446 tmpMaterial = new G4Material("ups923a",polystyrene->GetDensity(),1);
447 tmpMaterial->AddMaterial(polystyrene,1);
448 tmpMaterial->SetName("ups923a");
449 std::vector<G4double> ups923a_PhotonEnergy = {
450 3.35, 3.31, 3.28, 3.26, 3.25, 3.23, 3.23,
451 3.22, 3.21, 3.19, 3.18, 3.17, 3.16, 3.15,
452 3.14, 3.14, 3.13, 3.11, 3.1, 3.09, 3.09,
453 3.08, 3.07, 3.04, 3.02, 3.02, 3.01, 2.99,
454 2.98, 2.97, 2.97, 2.95, 2.95, 2.93, 2.93,
455 2.92, 2.92, 2.91, 2.89, 2.88, 2.87, 2.86,
456 2.85, 2.83, 2.81, 2.8, 2.79, 2.78, 2.76,
457 2.74, 2.72, 2.71, 2.68, 2.66, 2.64, 2.62,
458 2.61, 2.58, 2.55, 2.53, 2.5, 2.48, 2.46,
459 2.44, 2.41, 2.38, 2.35 };
460 std::reverse(ups923a_PhotonEnergy.begin(), ups923a_PhotonEnergy.end());
461
462 // AUG 21 - these were previously just one number as a const property but they should be non-const
463 // which requires vs energy numbers - so just use arrays the same shape
464 std::vector<G4double> ups923a_RINDEX(ups923a_PhotonEnergy.size(), 1.52);
465 std::vector<G4double> ups923a_ABSLENGTH(ups923a_PhotonEnergy.size(), 1*CLHEP::m);
466
467 G4MaterialPropertiesTable* ups923a_mt = CreatePropertiesTable();
468#if G4VERSION_NUMBER < 1070
469 ups923a_mt->AddProperty("RINDEX", ups923a_PhotonEnergy.data(), ups923a_RINDEX.data(), (int)ups923a_PhotonEnergy.size());
470 ups923a_mt->AddProperty("ABSLENGTH", ups923a_PhotonEnergy.data(), ups923a_ABSLENGTH.data(), (int)ups923a_PhotonEnergy.size());
471#else
472 ups923a_mt->AddProperty("RINDEX", ups923a_PhotonEnergy, ups923a_RINDEX);
473 ups923a_mt->AddProperty("ABSLENGTH", ups923a_PhotonEnergy, ups923a_ABSLENGTH);
474#endif
475 //Birk's constant
476 birks = (0.014/1.06)*CLHEP::cm/CLHEP::MeV;
477 tmpMaterial->GetIonisation()->SetBirksConstant(birks);
478#if G4VERSION_NUMBER < 1079
479 const G4int ups923a_numentries = 67;
480 G4double ups923a_emission[ups923a_numentries] = {
481 0, 0.04, 0.11, 0.2, 0.3, 0.4, 0.52,
482 0.62, 0.67, 0.68, 0.67, 0.62, 0.53, 0.48,
483 0.44, 0.42, 0.4, 0.41, 0.42, 0.51, 0.46,
484 0.57, 0.67, 0.78, 0.91, 0.93, 0.95, 0.96,
485 0.94, 0.91, 0.85, 0.76, 0.67, 0.61, 0.57,
486 0.55, 0.52, 0.51, 0.52, 0.54, 0.57, 0.58,
487 0.6, 0.6, 0.59, 0.58, 0.55, 0.48, 0.42,
488 0.37, 0.33, 0.31, 0.29, 0.28, 0.26, 0.24,
489 0.2, 0.17, 0.12, 0.09, 0.08, 0.07,
490 0.06, 0.04, 0.02, 0.01, 0.01 };
491 ups923a_mt->AddConstProperty("FASTTIMECONSTANT",3.3*CLHEP::ns);
492 ups923a_mt->AddProperty("FASTCOMPONENT",ups923a_PhotonEnergy.data(), ups923a_emission, ups923a_numentries)->SetSpline(true);
493 ups923a_mt->AddConstProperty("YIELDRATIO",1.0);
494#endif
495 ups923a_mt->AddConstProperty("RESOLUTIONSCALE",2.0); //Check this
496 G4double scintYieldAnthracene=14200; //Anthracene yield per 1 CLHEP::MeV
497 G4double scintYieldUPS923A=scintYieldAnthracene*0.60;//60% of anthracene
498 ups923a_mt->AddConstProperty("SCINTILLATIONYIELD",scintYieldUPS923A/CLHEP::MeV);
499
500 tmpMaterial->SetMaterialPropertiesTable(ups923a_mt);
501 AddMaterial(tmpMaterial, "ups923a");
502
503 // PET
504 G4double pet_density=1.4*CLHEP::g/CLHEP::cm3;
505 G4int pet_nelements=3;
506 G4State pet_state=kStateSolid;
507 tmpMaterial= new G4Material("pet",
508 pet_density,
509 pet_nelements,
510 pet_state);
511 tmpMaterial->AddElement(nistManager->FindOrBuildElement("C",true),10);
512 tmpMaterial->AddElement(nistManager->FindOrBuildElement("H",true),8);
513 tmpMaterial->AddElement(nistManager->FindOrBuildElement("O",true),4);
514 const G4int Pet_NUMENTRIES = 3; //Number of entries in the material properties table
515 G4double Pet_RIND[Pet_NUMENTRIES] = {1.570,1.570,1.570};//Assume constant refractive index.
516 G4double Pet_Energy[Pet_NUMENTRIES] = {2.0*CLHEP::eV,7.0*CLHEP::eV,7.14*CLHEP::eV}; //The energies.
517 G4MaterialPropertiesTable* petMaterialPropertiesTable = CreatePropertiesTable();
518 petMaterialPropertiesTable->AddProperty("RINDEX",Pet_Energy, Pet_RIND, Pet_NUMENTRIES);
519 tmpMaterial->SetMaterialPropertiesTable(petMaterialPropertiesTable);
520 AddMaterial(tmpMaterial, "pet");
521
522 // Opaque PET (Dacron)
523 tmpMaterial= new G4Material("pet_opaque",
524 pet_density,
525 pet_nelements,
526 pet_state);
527 tmpMaterial->AddElement(nistManager->FindOrBuildElement("C",true),10);
528 tmpMaterial->AddElement(nistManager->FindOrBuildElement("H",true),8);
529 tmpMaterial->AddElement(nistManager->FindOrBuildElement("O",true),4);
530 const G4int Pet_Opaque_NUMENTRIES = 3; //Number of entries in the material properties table
531 G4double Pet_Opaque_RIND[Pet_Opaque_NUMENTRIES] = {1.570,1.570,1.570};//Assume constant refractive index.
532 G4double Pet_Opaque_Energy[Pet_Opaque_NUMENTRIES] = {2.0*CLHEP::eV,7.0*CLHEP::eV,7.14*CLHEP::eV}; //The energies.
533 G4double Pet_Opaque_abslen[]={1*CLHEP::um, 1*CLHEP::um, 1*CLHEP::um};
534 G4MaterialPropertiesTable* pet_opaqueMaterialPropertiesTable = CreatePropertiesTable();
535 pet_opaqueMaterialPropertiesTable->AddProperty("RINDEX",Pet_Opaque_Energy, Pet_Opaque_RIND, Pet_Opaque_NUMENTRIES);
536 pet_opaqueMaterialPropertiesTable->AddProperty("ABSLENGTH",Pet_Opaque_Energy, Pet_Opaque_abslen, Pet_Opaque_NUMENTRIES);
537 tmpMaterial->SetMaterialPropertiesTable(pet_opaqueMaterialPropertiesTable);
538 AddMaterial(tmpMaterial, "pet_opaque");
539
540 // Gadolinium oxysulphate Gd_2 O_2 S
541 G4Material* GOS = nistManager->FindOrBuildMaterial("G4_GADOLINIUM_OXYSULFIDE",true,true);
542
543 // Ganolinium oxysulphate in a polyurethane elastomer (lanex)
544 G4double fill_factor=0.5;
545 G4double lanex_density=fill_factor*GOS->GetDensity()+(1-fill_factor)*GetMaterial("polyurethane")->GetDensity();
546 G4double gos_fraction_by_mass=fill_factor*GOS->GetDensity()/lanex_density;
547 G4double pur_fraction_by_mass=1-gos_fraction_by_mass;
548
549 tmpMaterial = new G4Material("lanex", lanex_density, 2);
550 tmpMaterial->AddMaterial(GOS, gos_fraction_by_mass);
551 tmpMaterial->AddMaterial(GetMaterial("polyurethane"), pur_fraction_by_mass);
552
553 G4MaterialPropertiesTable* mptLanex = CreatePropertiesTable();
554 const G4int nentLanex=2;
555 G4double rindex=1.50;//(1.82+1.50)/2.0;
556 G4double energytab[]={2.239*CLHEP::eV, 2.241*CLHEP::eV};
557 G4double rindextab[]={rindex, rindex};
558 G4double abslen[]={7*CLHEP::mm, 7*CLHEP::mm};
559 mptLanex->AddProperty("RINDEX",energytab, rindextab, nentLanex); //Average refractive index of bulk material
560 mptLanex->AddProperty("ABSLENGTH", energytab, abslen, nentLanex);
561#if G4VERSION_NUMBER < 1079
562 mptLanex->AddConstProperty("MIEHG", 60.3e-3*CLHEP::mm); // interface changed in V11
563 G4double emitspec[]={1.0, 1.0};
564 mptLanex->AddProperty("FASTCOMPONENT",energytab, emitspec, nentLanex);
565 mptLanex->AddConstProperty("FASTTIMECONSTANT", 1.*CLHEP::ns);
566#endif
567 G4double scintScalingFactor=1;
568 mptLanex->AddConstProperty("SCINTILLATIONYIELD",7.8e4/CLHEP::MeV);
569 mptLanex->AddConstProperty("RESOLUTIONSCALE",1.0);
570 mptLanex->AddConstProperty("MIEHG_FORWARD", 0.91);
571 mptLanex->AddConstProperty("MIEHG_BACKWARD", 0.91);
572 mptLanex->AddConstProperty("MIEHG_FORWARD_RATIO", 1.0);
573 tmpMaterial->SetMaterialPropertiesTable(mptLanex);
574 AddMaterial(tmpMaterial, "lanex");
575
576 //Ganolinium oxysulphate in a polyurethane elastomer (lanex) - version 2
577 tmpMaterial = new G4Material("lanex2", lanex_density, 2);
578 tmpMaterial->AddMaterial(GOS, gos_fraction_by_mass);
579 tmpMaterial->AddMaterial(GetMaterial("polyurethane"), pur_fraction_by_mass);
580 G4MaterialPropertiesTable* mptLanex2 = CreatePropertiesTable();
581 mptLanex2->AddProperty("RINDEX",energytab, rindextab, nentLanex); //Average refractive index of bulk material
582 mptLanex2->AddProperty("ABSLENGTH", energytab, abslen, nentLanex);
583#if G4VERSION_NUMBER < 1079
584 mptLanex2->AddConstProperty("MIEHG", 60.3e-3*CLHEP::mm);
585 mptLanex2->AddProperty("FASTCOMPONENT",energytab, emitspec, nentLanex);
586 mptLanex2->AddConstProperty("FASTTIMECONSTANT", 1.*CLHEP::ns);
587#endif
588 mptLanex2->AddConstProperty("SCINTILLATIONYIELD",8.9e4/CLHEP::MeV);
589 mptLanex2->AddConstProperty("RESOLUTIONSCALE",1.0);
590 mptLanex2->AddConstProperty("MIEHG_FORWARD", 0.91);
591 mptLanex2->AddConstProperty("MIEHG_BACKWARD", 0.91);
592 mptLanex2->AddConstProperty("MIEHG_FORWARD_RATIO", 0.5);
593 tmpMaterial->SetMaterialPropertiesTable(mptLanex);
594 AddMaterial(tmpMaterial, "lanex2");
595
596 //gos_lanex - GOS with the bulk optical transport properties of lanex particles suspended in an elastomer but the atomic, density and scintillation properties of GOS
597 G4double gos_lanex_density=GOS->GetDensity();
598 tmpMaterial = new G4Material("gos_lanex", gos_lanex_density, 1);
599 tmpMaterial->AddMaterial(GOS, 1.0);
600 G4MaterialPropertiesTable* mptGOSLanex = CreatePropertiesTable();
601 const G4int nentGOSLanex=2;
602 G4double rindexGOSLanex=1.50;
603 G4double energyGOSLanexTab[]={2.239*CLHEP::eV, 2.241*CLHEP::eV};
604 G4double rindexGOSLanexTab[]={rindexGOSLanex, rindexGOSLanex};
605 G4double abslenGOSLanex[]={7*CLHEP::mm, 7*CLHEP::mm};
606 G4double gosLanexMiehgForward=0.911;
607 G4double gosLanexMiehgBackward=0.911;
608 G4double gosLanexMiehgForwardRatio=0.5;
609#if G4VERSION_NUMBER < 1079
610 G4double mieHgTimeConst=1.0*CLHEP::ns;
611 G4double emitspecGOSLanex[]={1.0, 1.0};
612 G4double mieScatteringLengthGOSLanex=60.3*CLHEP::um;
613 mptGOSLanex->AddProperty("FASTCOMPONENT",energyGOSLanexTab, emitspecGOSLanex, nentGOSLanex);
614 mptGOSLanex->AddConstProperty("FASTTIMECONSTANT", mieHgTimeConst);
615 mptGOSLanex->AddConstProperty("YIELDRATIO", 1.0);
616 mptGOSLanex->AddConstProperty("MIEHG", mieScatteringLengthGOSLanex);
617#endif
618 mptGOSLanex->AddConstProperty("SCINTILLATIONYIELD",8.9e4/CLHEP::MeV); //Intrinisic scintilation yield of GOS
619 mptGOSLanex->AddConstProperty("RESOLUTIONSCALE", 1.0);
620 mptGOSLanex->AddConstProperty("MIEHG_FORWARD", gosLanexMiehgForward);
621 mptGOSLanex->AddConstProperty("MIEHG_BACKWARD", gosLanexMiehgBackward);
622 mptGOSLanex->AddConstProperty("MIEHG_FORWARD_RATIO", gosLanexMiehgForwardRatio);
623 mptGOSLanex->AddProperty("RINDEX",energyGOSLanexTab, rindexGOSLanexTab, nentGOSLanex); //Average refractive index of bulk material
624 mptGOSLanex->AddProperty("ABSLENGTH", energyGOSLanexTab, abslenGOSLanex, nentGOSLanex);
625 tmpMaterial->SetMaterialPropertiesTable(mptGOSLanex);
626 AddMaterial(tmpMaterial, "gos_lanex");
627
628 //Same as gos_lanex but refractive index = 1
629 tmpMaterial = new G4Material("gos_ri1", gos_lanex_density, 1);
630 tmpMaterial->AddMaterial(GOS, 1.0);
631 G4MaterialPropertiesTable* mptGOSLanexRi1 = CreatePropertiesTable();
632 G4double rindexGOSLanexRi1Tab[]={1.0, 1.0};
633#if G4VERSION_NUMBER < 1079
634 mptGOSLanexRi1->AddProperty("FASTCOMPONENT",energyGOSLanexTab, emitspecGOSLanex, nentGOSLanex);
635 mptGOSLanexRi1->AddConstProperty("FASTTIMECONSTANT", mieHgTimeConst);
636 mptGOSLanexRi1->AddConstProperty("YIELDRATIO", 1.0);
637 mptGOSLanexRi1->AddConstProperty("MIEHG", mieScatteringLengthGOSLanex);
638#endif
639 mptGOSLanexRi1->AddConstProperty("SCINTILLATIONYIELD",8.9e4/CLHEP::MeV); //Intrinisic scintilation yield of GOS
640 mptGOSLanexRi1->AddConstProperty("RESOLUTIONSCALE", 1.0);
641 mptGOSLanexRi1->AddConstProperty("MIEHG_FORWARD", gosLanexMiehgForward);
642 mptGOSLanexRi1->AddConstProperty("MIEHG_BACKWARD", gosLanexMiehgBackward);
643 mptGOSLanexRi1->AddConstProperty("MIEHG_FORWARD_RATIO", gosLanexMiehgForwardRatio);
644 mptGOSLanexRi1->AddProperty("RINDEX",energyGOSLanexTab, rindexGOSLanexRi1Tab, nentGOSLanex); //Average refractive index of bulk material
645 mptGOSLanexRi1->AddProperty("ABSLENGTH", energyGOSLanexTab, abslenGOSLanex, nentGOSLanex);
646 tmpMaterial->SetMaterialPropertiesTable(mptGOSLanexRi1);
647 AddMaterial(tmpMaterial, "gos_ri1");
648
649 //pet_lanex - PET with the bulk optical transport properties of lanex particles suspended in an elastomer but the atomic, density and scintillation properties of PET
650 G4double pet_lanex_density=GetMaterial("polyurethane")->GetDensity();
651 tmpMaterial = new G4Material("pet_lanex", pet_lanex_density, 1);
652 tmpMaterial->AddMaterial(GetMaterial("polyurethane"), 1.0);
653 G4MaterialPropertiesTable* mptPETLanex = CreatePropertiesTable();
654 mptPETLanex->AddConstProperty("MIEHG_FORWARD", gosLanexMiehgForward);
655 mptPETLanex->AddConstProperty("MIEHG_BACKWARD", gosLanexMiehgBackward);
656 mptPETLanex->AddConstProperty("MIEHG_FORWARD_RATIO", gosLanexMiehgForwardRatio);
657#if G4VERSION_NUMBER < 1079
658 mptPETLanex->AddConstProperty("MIEHG", mieScatteringLengthGOSLanex);
659#endif
660 mptPETLanex->AddProperty("RINDEX",energyGOSLanexTab, rindexGOSLanexTab, nentGOSLanex); //Average refractive index of bulk material
661 mptPETLanex->AddProperty("ABSLENGTH", energyGOSLanexTab, abslenGOSLanex, nentGOSLanex);
662 tmpMaterial->SetMaterialPropertiesTable(mptPETLanex);
663 AddMaterial(tmpMaterial, "pet_lanex");
664
665 //Medex (larger grained lanex)
666 // G4double medex_fill_factor=0.5;
667 G4double medex_density=fill_factor*GOS->GetDensity()+(1-fill_factor)*GetMaterial("polyurethane")->GetDensity();
668 G4double medex_gos_fraction_by_mass=fill_factor*GOS->GetDensity()/medex_density;
669 G4double medex_pur_fraction_by_mass=1-medex_gos_fraction_by_mass;
670 tmpMaterial = new G4Material("medex", medex_density, 2);
671 tmpMaterial->AddMaterial(GOS, medex_gos_fraction_by_mass);
672 tmpMaterial->AddMaterial(GetMaterial("polyurethane"), medex_pur_fraction_by_mass);
673 G4MaterialPropertiesTable* mptMedex = CreatePropertiesTable();
674 const G4int nentMedex=2;
675 // G4double medexRindex=(1.82+1.50)/2.0;
676 // G4double medexEnergytab[]={2.239*CLHEP::eV, 2.241*CLHEP::eV};
677 G4double medexRindextab[]={rindex, rindex};
678
679 G4double medexAbslen[]={7*CLHEP::mm, 7*CLHEP::mm};
680 mptMedex->AddProperty("RINDEX",energytab, medexRindextab, nentMedex); //Average refractive index of bulk material
681 mptMedex->AddProperty("ABSLENGTH", energytab, medexAbslen, nentMedex);
682#if G4VERSION_NUMBER < 1079
683 mptMedex->AddConstProperty("MIEHG", 230e-3*CLHEP::mm);
684 G4double medexEmitspec[]={1.0, 1.0};
685 mptMedex->AddProperty("FASTCOMPONENT",energytab, medexEmitspec, nentMedex);
686 mptMedex->AddConstProperty("FASTTIMECONSTANT", 1.*CLHEP::ns);
687#endif
688 mptMedex->AddConstProperty("SCINTILLATIONYIELD",scintScalingFactor*2.94e4/CLHEP::MeV);
689 mptMedex->AddConstProperty("RESOLUTIONSCALE",1.0);
690 mptMedex->AddConstProperty("MIEHG_FORWARD", 0.93);
691 mptMedex->AddConstProperty("MIEHG_BACKWARD", 0.93);
692 mptMedex->AddConstProperty("MIEHG_FORWARD_RATIO", 1.0);
693 tmpMaterial->SetMaterialPropertiesTable(mptMedex);
694 AddMaterial(tmpMaterial, "medex");
695
696 // carbon fiber
697 G4double frac_graph=0.5;
698 G4double frac_poly=0.5;
699 G4Material* graph=GetMaterial("G4_GRAPHITE");
700 G4Material* poly=GetMaterial("G4_POLYACRYLONITRILE");
701 G4double dens_graph=graph->GetDensity();
702 G4double dens_poly=poly->GetDensity();
703
704 G4double dens_cf =frac_graph*dens_graph+frac_poly*dens_poly;
705 G4double frac_graph_bw=frac_graph*dens_graph/dens_cf;
706 G4double frac_poly_bw=frac_poly*dens_poly/dens_cf;
707 tmpMaterial = new G4Material("carbonfiber", dens_cf, 2);
708 tmpMaterial->AddMaterial(graph, frac_graph_bw);
709 tmpMaterial->AddMaterial(poly, frac_poly_bw);
710 G4MaterialPropertiesTable* mptCarbonfiber = CreatePropertiesTable();
711 const G4int nentCarbonfiber=2;
712 G4double energytab_cf[]={2.239*CLHEP::eV, 2.241*CLHEP::eV};
713 G4double carbonfiberRindextab[]={2.6, 2.6};
714 G4double carbonfiberAbslen[]={2.1*CLHEP::um, 2.1*CLHEP::um};
715 mptCarbonfiber->AddProperty("RINDEX",energytab_cf, carbonfiberRindextab, nentCarbonfiber); //Average refractive index of bulk material
716 mptCarbonfiber->AddProperty("ABSLENGTH", energytab_cf, carbonfiberAbslen, nentCarbonfiber);
717 tmpMaterial->SetMaterialPropertiesTable(mptCarbonfiber);
718 AddMaterial(tmpMaterial, "carbonfiber");
719}
720
721void BDSMaterials::DefineLHCComponents()
722{
723 AddMaterial("lhcconcrete",
724 2.42,
725 kStateSolid, 300, 1,
726 {"H","C","O","Na","Mg","Al","Si","K","Ca","Fe","P","S",
727 "Ti","Mn","Zn","Zr","Ba","Pb","Sr","Eu"},
728 std::list<double>{0.59785345499811 *CLHEP::perCent,
729 5.59989402848226 *CLHEP::perCent,
730 49.1111702720319 *CLHEP::perCent,
731 0.45137935852357 *CLHEP::perCent,
732 0.66062806777291 *CLHEP::perCent,
733 2.05561946276849 *CLHEP::perCent,
734 18.7995018924154 *CLHEP::perCent,
735 0.65365311079793 *CLHEP::perCent,
736 20.0191229406116 *CLHEP::perCent,
737 1.11400027114647 *CLHEP::perCent,
738 0.04782827639985 *CLHEP::perCent,
739 0.01195706909996 *CLHEP::perCent,
740 0.3457585814739 *CLHEP::perCent,
741 0.03856154784738 *CLHEP::perCent,
742 0.02401378044242 *CLHEP::perCent,
743 0.00737352594498 *CLHEP::perCent,
744 0.01783596140744 *CLHEP::perCent,
745 0.04623400051985 *CLHEP::perCent,
746 0.39757254757374 *CLHEP::perCent,
747 4.184974185E-05 *CLHEP::perCent});
748
749 // LHC Rock
750 AddMaterial("lhc_rock",
751 2,
752 kStateSolid, 300, 1,
753 {"H", "Na", "Si", "Fe", "C", "Mg", "K", "O", "Al", "Ca"},
754 std::list<double>{0.006, 0.01, 0.2, 0.014, 0.03, 0.005, 0.01, 0.5, 0.03, 0.195});
755
756 // Superconducting components of LHC magnet elements
757 // Definitions taken from FLUKA
758
759 std::list<int> singleElement = {1};
760 // Liquid helium at 1.9K
761 AddMaterial("lhe_1.9k",
762 0.1472,
763 kStateLiquid, 1.9, 1,
764 {"He"}, singleElement);
765
766 // Niobium @ 87K
767 AddMaterial("nb_87k",
768 8.902 ,
769 kStateSolid , 87 , 1,
770 {"Nb"}, singleElement);
771
772 // Titanium @ 87K
773 AddMaterial("ti_87k",
774 4.54,
775 kStateSolid, 87 , 1,
776 {"Ti"}, singleElement);
777
778 // superconductor NbTi with Ti = 47% by weight
779 AddMaterial("nbti_87k",
780 6.0471,
781 kStateSolid , 87 , 1,
782 {"nb_87k","ti_87k"},
783 std::list<double>{0.53,0.47});
784
785 // copper at 4 Kelvin
786 AddMaterial("cu_4k",
787 8.96,
788 kStateSolid, 4, 1,
789 {"Cu"}, singleElement);
790
791 // copper at 2 Kelvin
792 AddMaterial("cu_2k",
793 8.96,
794 kStateSolid, 2, 1,
795 {"Cu"}, singleElement);
796
797 // naked superconductor NbTi wire with Cu/SC volume ratio (>= 4.0 and <4.8)
798 AddMaterial("nbti.1",
799 8.4206,
800 kStateSolid, 4, 1,
801 {"nbti_87k","cu_4k"},
802 std::list<double>{1.0/5.4, 4.4/5.4});
803
804 // For HL-LHC Collimation - Copper Diamond - TBC - component fractions
805 AddMaterial("copperdiamond",
806 5.3,
807 kStateSolid, 300, 1,
808 {"Cu", "C"},
809 std::list<double>{0.32, 0.68});
810 AddExistingMaterialAlias("copperdiamond", "cucd");
811
812 // For HL-LHC Collimation - Copper Diamond
813 AddMaterial("molybdenumcarbide",
814 8.9,
815 kStateSolid, 300, 1,
816 {"Mo", "C"},
817 std::list<int>{1, 1});
818 AddExistingMaterialAlias("molybdenumcarbide", "mogr");
819}
820
821void BDSMaterials::DefineLiquids()
822{
823 AddMaterial("liquidhelium", 0.12498, kStateLiquid, 4.15, 1, {"He"}, std::list<int>{1});
824
825 // water
826 AddExistingMaterialAlias("G4_WATER", "water"); // use NIST water
827 G4Material* water = GetMaterial("water"); // get the material pointer back
828
829 // add refractive index to G4_WATER material for cherenkov radiation to work
830 const G4int nEntries = 9;
831 G4MaterialPropertiesTable* waterProperties = CreatePropertiesTable();
832 G4double photonEnergy[nEntries];
833 G4double dNEntries = (G4double)nEntries;
834 G4double energyMin = 1.*CLHEP::eV;
835 G4double energyMax = 3.*CLHEP::eV;
836 G4double deltaEnergy = (energyMax-energyMin)/(dNEntries-1.0);
837 G4double energy = energyMin;
838 for (G4int i = 0; i < nEntries; energy += deltaEnergy, i++)
839 {photonEnergy[i] = energy;}
840 G4double refractiveIndex[nEntries] = {1.325, 1.325, 1.326, 1.327, 1.328, 1.33, 1.333, 1.336, 1.343};
841 waterProperties->AddProperty("RINDEX", photonEnergy, refractiveIndex, nEntries);
842 water->SetMaterialPropertiesTable(waterProperties);
843}
844
845void BDSMaterials::DefineGases()
846{
847 AddExistingMaterialAlias("N", "nitrogen");
848
849 // air
850 AddExistingMaterialAlias("G4_AIR", "air");
851 G4Material* g4AirMaterial = GetMaterial("air");
852 const G4int airNEntries = 3;
853 // source: NPL Tables of Physical & Chemical Constants. Refractive indices at
854 // different energies.
855 G4double airRefractiveIndex[airNEntries] = {1.000292,1.000292,1.000292};
856 G4double airEnergy[airNEntries] = {2.0*CLHEP::eV,7.0*CLHEP::eV,7.14*CLHEP::eV};
857 G4MaterialPropertiesTable* airProperties = CreatePropertiesTable();
858 airProperties->AddProperty("RINDEX", airEnergy, airRefractiveIndex, airNEntries);
859 g4AirMaterial->SetMaterialPropertiesTable(airProperties);
860
861 // old way - manually constructed air
862 // room temperature temperature and standard pressure
863 G4double temperature = 300*CLHEP::kelvin;
864 G4double pressure = 1.0*CLHEP::atmosphere;
865 G4double airDensity = 0.001225; // g/cm3
866 G4Material* tmpMaterial = new G4Material("airbdsim",
867 airDensity*CLHEP::g/CLHEP::cm3,
868 2,
869 kStateGas,
870 temperature,
871 pressure);
872 tmpMaterial->AddElement(GetElement("O"), 0.2);
873 tmpMaterial->AddElement(GetElement("N"), 0.8);
874 tmpMaterial->SetMaterialPropertiesTable(airProperties);
875 AddMaterial(tmpMaterial, "airbdsim");
876
877 // Carbon monoxide
878 G4double coDensity = 0.001145; // g/cm3
879 AddMaterial("carbonmonoxide",
880 coDensity,
881 kStateGas,
882 temperature,
883 pressure,
884 {"C","O"},
885 std::list<int>{1,1});
886
887 // Carbon monoxide beam pipe gas
888 G4double bp_pressure = 0.0133e-9 * (CLHEP::bar/CLHEP::atmosphere); //10 nTorr pressure
889 AddMaterial("bp_carbonmonoxide",
890 coDensity,
891 kStateGas,
892 temperature,
893 bp_pressure,
894 {"C","O"},
895 std::list<int>{1,1});
896}
897
898void BDSMaterials::DefinePlasmas()
899{
900 //Awake plasma - rubidium at density of 7e14 atoms/cm3
901 // G4double numberDensity = 7.0e14/CLHEP::cm3;
902 G4double a = 85.4678*CLHEP::g/CLHEP::mole;
903 G4double density = 1e-7 * CLHEP::g/CLHEP::cm3;
904 G4Material* tmpMaterial = new G4Material("awakeplasma", 37., a, density);
905 AddMaterial(tmpMaterial, "awakeplasma");
906}
907
908void BDSMaterials::DefineVacuums()
909{
910 // Default vacuum (same composition as residual vacuum in warm sections of LHC).
911 G4double vacpressure;
912 // check if parser is initialised (in case list materials is asked by command line)
913 if (BDSParser::IsInitialised())
914 {vacpressure=BDSParser::Instance()->GetOptions().vacuumPressure*CLHEP::bar;}
915 else
916 {vacpressure=1e-12*CLHEP::bar;}
917 G4double temperature = 300*CLHEP::kelvin;
918 G4double density = (CLHEP::STP_Temperature/temperature) * (vacpressure/(1.*CLHEP::atmosphere)) * 29*CLHEP::g/(22.4*1.e-3*CLHEP::m3);
919
920 G4Material* tmpMaterial = new G4Material("vacuum",
921 density,
922 3, kStateGas,
923 temperature, vacpressure);
924 tmpMaterial->AddElement(GetElement("H"), 0.482);
925 tmpMaterial->AddElement(GetElement("C"), 0.221);
926 tmpMaterial->AddElement(GetElement("O"), 0.297);
927 AddMaterial(tmpMaterial, "vacuum");
928
929 const G4int Vac_NUMENTRIES = 3;
930 // Assume refractive index = 1 in a vacuum.
931 G4double Vac_RIND[Vac_NUMENTRIES] = {1.000,1.000,1.000};
932 G4double Vac_Energy[Vac_NUMENTRIES] = {2.0*CLHEP::eV,7.0*CLHEP::eV,7.14*CLHEP::eV};
933 G4MaterialPropertiesTable* vacMaterialPropertiesTable = CreatePropertiesTable();
934 vacMaterialPropertiesTable->AddProperty("RINDEX",Vac_Energy, Vac_RIND, Vac_NUMENTRIES);
935 tmpMaterial->SetMaterialPropertiesTable(vacMaterialPropertiesTable);
936
937 // copy regular vacuum as named 'laservac' to detect where laserwire is.
938 // G4Material doesn't have a copy constructor but it has a sort of copy
939 // constructor.
940 G4Material* regularVacuum = GetMaterial("vacuum");
941 G4Material* laservac = new G4Material("laservac",
942 regularVacuum->GetDensity(),
943 regularVacuum,
944 kStateGas,
945 regularVacuum->GetTemperature(),
946 regularVacuum->GetPressure());
947 AddMaterial(laservac, "laservac");
948}
949
950void BDSMaterials::AddMaterial(G4Material* material, G4String name)
951{
952 name = BDS::LowerCase(name);
953 if (materials.insert(make_pair(name, material)).second)
954 {
955#ifdef BDSDEBUG
956 G4cout << "New material : " << name << G4endl;
957#endif
958 }
959 else
960 {throw BDSException(__METHOD_NAME__, "Material \"" + name + "\" already exists");}
961}
962
963void BDSMaterials::AddExistingMaterialAlias(const G4String &existingMaterialName,
964 G4String alias)
965{
966 alias = BDS::LowerCase(alias);
967 G4Material* material = GetMaterial(existingMaterialName);
968 aliases[alias] = material; // store in lower case as that's how we search
969}
970
971void BDSMaterials::AddMaterial(G4String name,
972 G4double Z,
973 G4double A,
974 G4double density,
975 G4State state,
976 G4double temperature,
977 G4double pressure)
978{
979 // convention: material name in small letters (to be able to find materials regardless of capitalisation)
980 name = BDS::LowerCase(name);
981 DensityCheck(density, name);
982
983 G4Material* tmpMaterial = new G4Material(name,
984 Z,
985 A*CLHEP::g/CLHEP::mole,
986 density*CLHEP::g/CLHEP::cm3,
987 state,
988 temperature*CLHEP::kelvin,
989 pressure*CLHEP::atmosphere);
990 AddMaterial(tmpMaterial, name);
991}
992
993template <typename Type>
994void BDSMaterials::AddMaterial(G4String name,
995 G4double density,
996 G4State state,
997 G4double temperature,
998 G4double pressure,
999 const std::list<G4String>& components,
1000 const std::list<Type>& componentFractions)
1001{
1002 name = BDS::LowerCase(name);
1003 DensityCheck(density, name);
1004
1005 G4Material* tmpMaterial = new G4Material(name,
1006 density*CLHEP::g/CLHEP::cm3,
1007 (G4int)components.size(),
1008 state,
1009 temperature*CLHEP::kelvin,
1010 pressure*CLHEP::atmosphere);
1011 std::list<G4String>::const_iterator sIter;
1012 typename std::list<Type>::const_iterator dIter;
1013 for (sIter = components.begin(), dIter = componentFractions.begin();
1014 sIter != components.end();
1015 ++sIter, ++dIter)
1016 {
1017#ifdef BDSDEBUG
1018 G4cout << "BDSMaterials::AddMaterial: " << *sIter << G4endl;
1019#endif
1020 G4Element* element = CheckElement(*sIter);
1021 if (element)
1022 {tmpMaterial->AddElement(element, (*dIter));}
1023 else
1024 {tmpMaterial->AddMaterial(GetMaterial(*sIter), (*dIter));}
1025 }
1026 AddMaterial(tmpMaterial, name);
1027}
1028
1029G4Material* BDSMaterials::GetMaterial(G4String material) const
1030{
1031 if (material.empty())
1032 {throw BDSException(__METHOD_NAME__, "empty material name");}
1033 G4String materialOriginal = material;
1034 // for short names we assume they're elements so we prefix with G4_ and
1035 // get them from NIST
1036 G4String nistString ("G4_");
1037 if (material.length() <= 2)
1038#if G4VERSION_NUMBER > 1099
1039 {material = nistString + material;}
1040#else
1041 {material.prepend(nistString);}
1042#endif
1043
1044 G4String start = material.substr(0,3);
1045 if (nistString == start)
1046 {
1047#ifdef BDSDEBUG
1048 G4cout << "Using NIST material " << material << G4endl;
1049#endif
1050 G4Material* mat = G4NistManager::Instance()->FindOrBuildMaterial(material, true, true);
1051 if (!mat)
1052 {throw BDSException(__METHOD_NAME__, "\"" + material + "\" could not be found by NIST.");}
1053 return mat;
1054 }
1055 else
1056 {
1057 // find material regardless of capitalisation
1058 material = BDS::LowerCase(material);
1059 auto search = possibleDuplicates.find(material);
1060 if (search != possibleDuplicates.end())
1061 {
1062 if (search->second > 1)
1063 {
1064 throw BDSException(__METHOD_NAME__, "material \"" + materialOriginal
1065 + "\" has been loaded from multiple GDML files and is ambiguous.\n"
1066 + "Please prepend with the BDSIM element used to load the file to be explicit.");
1067 }
1068 }
1069 auto iter = materials.find(material);
1070 if (iter != materials.end())
1071 {return (*iter).second;}
1072 else
1073 {
1074 // search aliases
1075 auto iter2 = aliases.find(material);
1076 if (iter2 != aliases.end())
1077 {return iter2->second;}
1078 else
1079 {// can't find it -> warn and exit
1080 ListMaterials();
1081 throw BDSException(__METHOD_NAME__, "\"" + materialOriginal + "\" is unknown.");
1082 }
1083 }
1084 }
1085}
1086
1087void BDSMaterials::CacheMaterialsFromGDML(const std::map<G4String, G4Material*>& materialsGDML,
1088 const G4String& prepend,
1089 G4bool prependWasUsed)
1090{
1091 // Register in "aliases" member so we don't try to delete - geant4 will
1092 // do this for ones loaded in GDML. Therefore, avoid double deletion.
1093 for (const auto& kv : materialsGDML)
1094 {
1095 G4String nameLower = BDS::LowerCase(kv.first);
1096 //G4bool startsWithPrepend = prependExists ? BDS::StartsWith(kv.first, prepend) : false;
1097 if (BDS::StartsWith(nameLower, "g4_") || materials.find(nameLower) != materials.end())
1098 {continue;} // a Geant4 material or a BDSIM one
1099 aliases[nameLower] = kv.second;
1100
1101 if (prependWasUsed)
1102 {// cache without prefix
1103 G4String nameCopy = kv.first;
1104 nameCopy.erase(0, prepend.size() + 1);
1105 nameCopy = BDS::LowerCase(nameCopy);
1106 aliases[nameCopy] = kv.second;
1107 possibleDuplicates[nameCopy]++;
1108 }
1109 }
1110}
1111
1112void BDSMaterials::AddElement(G4Element* element, const G4String& symbol)
1113{
1114 if (CheckElement(symbol) != nullptr)
1115 {throw BDSException(__METHOD_NAME__, "Element \"" + symbol + "\" already exists.");}
1116
1117 elements.insert(make_pair(symbol, element));
1118#ifdef BDSDEBUG
1119 G4cout << "New element : " << symbol << G4endl;
1120#endif
1121}
1122
1123void BDSMaterials::AddElement(const G4String& name,
1124 const G4String& symbol,
1125 G4double Z,
1126 G4double A)
1127{
1128 G4Element* tmpElement = new G4Element(name, symbol, Z, A*CLHEP::g/CLHEP::mole);
1129 AddElement(tmpElement, symbol);
1130}
1131
1132void BDSMaterials::DensityCheck(G4double density,
1133 const G4String& materialName) const
1134{
1135 if (density > 1e2)
1136 {// so greater than 100g / cm3, the densest natural material is around 23g/cm3
1137 G4String msg = "material \"" + materialName + "\"has a density higher than 100g/cm3! Perhaps check this!\n";
1138 msg += "Density: " + std::to_string(density) + " g/cm3... proceeding";
1139 BDS::Warning(__METHOD_NAME__, msg);
1140 }
1141}
1142
1143G4Element* BDSMaterials::CheckElement(const G4String& symbol) const
1144{
1145 // first look in defined element list
1146 auto iter = elements.find(symbol);
1147 if(iter != elements.end())
1148 {return (*iter).second;}
1149 else
1150 {
1151 G4Element* element = G4NistManager::Instance()->FindOrBuildElement(symbol, true);
1152 return element;
1153 }
1154}
1155
1156G4Element* BDSMaterials::GetElement(const G4String& symbol) const
1157{
1158 G4Element* element = CheckElement(symbol);
1159 if (!element)
1160 {throw BDSException(__METHOD_NAME__, "Element \"" + symbol + "\" could not be found.");}
1161 return element;
1162}
1163
1164void BDSMaterials::PrintBasicMaterialMassFraction(G4Material* material) const
1165{
1166 // simpler version of geant4 code
1167 for (G4int i = 0; i < (G4int)material->GetNumberOfElements(); i++)
1168 {
1169 G4cout << (*(material->GetElementVector()))[i]->GetName() << "\t "
1170 << (material->GetFractionVector()[i])/CLHEP::perCent << " %" << G4endl;
1171 }
1172}
1173
1174void BDSMaterials::ListMaterials() const
1175{
1176 auto flagsCache(G4cout.flags());
1177 G4cout << __METHOD_NAME__ << G4endl;
1178 // always print out vacuum composition
1179 G4Material* vacuum = GetMaterial("vacuum");
1180 G4cout<< "\"vacuum\" composition: " << G4endl;
1181 PrintBasicMaterialMassFraction(vacuum);
1182 G4cout<< "pressure = " << vacuum->GetPressure()/CLHEP::bar << " bar" << G4endl;
1183 G4cout<< "temperature = " << vacuum->GetTemperature()/CLHEP::kelvin << " K" << G4endl;
1184 G4cout<< "density = " << vacuum->GetDensity()/(CLHEP::g/CLHEP::cm3)<< " g/cm^3" << G4endl << G4endl;
1185
1186 G4cout << "All elements are available with their 1 or 2 letter chemical symbol. ie C or G4_C" << G4endl << G4endl;
1187
1188 if (!elements.empty())
1189 {
1190 G4cout << "Extra defined elements are:" << G4endl;
1191 for (const auto& element : elements)
1192 {G4cout << std::left << std::setw(12) << element.second->GetName() << " - " << element.second->GetSymbol() << G4endl;}
1193 G4cout << G4endl;
1194 }
1195
1196 G4cout << "Defined materials are:" << G4endl;
1197 for (const auto& material : materials)
1198 {
1199 G4cout << material.first;
1200 G4String realName = material.second->GetName();
1201 if (realName != material.first)
1202 {G4cout << " (" << material.second->GetName() << ")" << G4endl;}
1203 else
1204 {G4cout << G4endl;}
1205 }
1206 G4cout << G4endl;
1207 G4cout << "Available NIST materials are:" << G4endl;
1208 G4NistManager::Instance()->ListMaterials("all");
1209 G4cout.flags(flagsCache);
1210}
1211
1212BDSMaterials::~BDSMaterials()
1213{
1214 for (auto& material : materials)
1215 {delete material.second;}
1216 materials.clear();
1217
1218 for (auto& element : elements)
1219 {delete element.second;}
1220 elements.clear();
1221
1222 for (G4MaterialPropertiesTable* table : propertiesTables)
1223 {delete table;}
1224
1225 instance = nullptr;
1226}
1227
1228void BDSMaterials::PrepareRequiredMaterials(G4bool verbose)
1229{
1230 // This function uses the list from the parser and prepares
1231 // the necessary materials for this run.
1232#ifdef BDSDEBUG
1233 G4bool debug = true;
1234#else
1235 G4bool debug = false;
1236#endif
1237
1238 // convert the parsed atom list to list of Geant4 G4Elements
1239 if (verbose || debug)
1240 {G4cout << __METHOD_NAME__ << "parsing the atom list..." << G4endl;}
1241 for (const auto& it : BDSParser::Instance()->GetAtoms())
1242 {AddElement(it.name,it.symbol,it.Z,it.A);}
1243 if (verbose || debug)
1244 {G4cout << "size of atom list: "<< BDSParser::Instance()->GetAtoms().size() << G4endl;}
1245
1246 // convert the parsed material list to list of Geant4 G4Materials
1247 if (verbose || debug)
1248 {G4cout << "parsing the material list..."<< G4endl;}
1249 for (auto it : BDSParser::Instance()->GetMaterials())
1250 {
1251 G4State itsState;
1252 if (it.state=="solid")
1253 {itsState = kStateSolid;}
1254 else if (it.state=="liquid")
1255 {itsState = kStateLiquid;}
1256 else if (it.state=="gas")
1257 {itsState = kStateGas;}
1258 else
1259 {
1260 G4cout << "Unknown material state "<< it.state
1261 << ", setting it to default (solid)"
1262 << G4endl;
1263 it.state="solid";
1264 itsState = kStateSolid;
1265 }
1266
1267#ifdef BDSDEBUG
1268 G4cout << "---->adding Material, ";
1269 it.print();
1270#endif
1271
1272 if(it.Z != 0)
1273 {
1274 AddMaterial(it.name,
1275 it.Z,
1276 it.A,
1277 it.density,
1278 itsState,
1279 it.temper,
1280 it.pressure);
1281 }
1282 else if(!(it.components.empty()))
1283 {
1284 std::list<G4String> tempComponents;
1285 for (const auto& jt : it.components)
1286 {tempComponents.emplace_back(G4String(jt));}
1287
1288 if(it.componentsWeights.size()==it.components.size())
1289 {
1290 AddMaterial(it.name,
1291 it.density,
1292 itsState,
1293 it.temper,
1294 it.pressure,
1295 tempComponents,
1296 it.componentsWeights);
1297 }
1298 else if(it.componentsFractions.size()==it.components.size())
1299 {
1300 AddMaterial(it.name,
1301 it.density,
1302 itsState,
1303 it.temper,
1304 it.pressure,
1305 tempComponents,
1306 it.componentsFractions);
1307 }
1308 else
1309 {throw BDSException(__METHOD_NAME__, "Badly defined material - number of components is not equal to number of weights or mass fractions!");}
1310 }
1311 else
1312 {throw BDSException(__METHOD_NAME__, "Badly defined material - need more information!");}
1313 }
1314 if (verbose || debug)
1315 {G4cout << "size of material list: "<< BDSParser::Instance()->GetMaterials().size() << G4endl;}
1316}
1317
1318G4MaterialPropertiesTable* BDSMaterials::CreatePropertiesTable()
1319{
1320 G4MaterialPropertiesTable* table = new G4MaterialPropertiesTable();
1321 propertiesTables.push_back(table);
1322 return table;
1323}
BDSException
General exception with possible name of object and message.
Definition: BDSException.hh:35
BDSMaterials
A class for all material definitions known to BDSIM. Additional materials can be added in the parser ...
Definition: BDSMaterials.hh:38
BDSMaterials::DefineMetals
void DefineMetals()
Methods called by constructor.
Definition: BDSMaterials.cc:66
BDSMaterials::DefineLHCComponents
void DefineLHCComponents()
Methods called by constructor.
Definition: BDSMaterials.cc:721
BDSMaterials::DefineGases
void DefineGases()
Methods called by constructor.
Definition: BDSMaterials.cc:845
BDSMaterials::DefineNonMetalSolids
void DefineNonMetalSolids()
Methods called by constructor.
Definition: BDSMaterials.cc:245
BDSMaterials::DensityCheck
void DensityCheck(G4double density, const G4String &materialName) const
Print warning if density suspiciously high. Should be in g/cm3 in G4 units already.
Definition: BDSMaterials.cc:1132
BDSMaterials::DefineScintillators
void DefineScintillators()
Methods called by constructor.
Definition: BDSMaterials.cc:403
BDSMaterials::DefineVacuums
void DefineVacuums()
Methods called by constructor.
Definition: BDSMaterials.cc:908
BDSMaterials::Instance
static BDSMaterials * Instance()
Singleton pattern access.
Definition: BDSMaterials.cc:38
BDSMaterials::CreatePropertiesTable
G4MaterialPropertiesTable * CreatePropertiesTable()
Create new properties table and store in vector.
Definition: BDSMaterials.cc:1318
BDSMaterials::aliases
std::map< G4String, G4Material * > aliases
Maps of other names to existing materials. To avoid double deletion. Also in lower case.
Definition: BDSMaterials.hh:139
BDSMaterials::instance
static BDSMaterials * instance
Singleton instance.
Definition: BDSMaterials.hh:110
BDSMaterials::materials
std::map< G4String, G4Material * > materials
<ap of materials, convention name lowercase.
Definition: BDSMaterials.hh:137
BDSMaterials::elements
std::map< G4String, G4Element * > elements
Map of elements, no lowercase convention.
Definition: BDSMaterials.hh:144
BDSMaterials::DefineLiquids
void DefineLiquids()
Methods called by constructor.
Definition: BDSMaterials.cc:821
BDSMaterials::CacheMaterialsFromGDML
void CacheMaterialsFromGDML(const std::map< G4String, G4Material * > &materialsGDML, const G4String &prepend, G4bool prependWasUsed)
Definition: BDSMaterials.cc:1087
BDSMaterials::GetMaterial
G4Material * GetMaterial(G4String material) const
Get material by name.
Definition: BDSMaterials.cc:1029
BDSMaterials::PrepareRequiredMaterials
void PrepareRequiredMaterials(G4bool verbose=false)
converts parser material list
Definition: BDSMaterials.cc:1228
BDSMaterials::ListMaterials
void ListMaterials() const
output available materials
Definition: BDSMaterials.cc:1174
BDSMaterials::DefinePlasmas
void DefinePlasmas()
Methods called by constructor.
Definition: BDSMaterials.cc:898
BDSMaterials::AddElement
void AddElement(G4Element *element, const G4String &symbol)
Add a G4Element.
Definition: BDSMaterials.cc:1112
BDSMaterials::AddMaterial
void AddMaterial(G4Material *material, G4String name)
Add G4Material.
Definition: BDSMaterials.cc:950
BDSMaterials::DefineSuperconductors
void DefineSuperconductors()
Methods called by constructor.
Definition: BDSMaterials.cc:236
BDSMaterials::CheckElement
G4Element * CheckElement(const G4String &symbol) const
Return element if defined (nullptr if not)
Definition: BDSMaterials.cc:1143
BDSMaterials::AddExistingMaterialAlias
void AddExistingMaterialAlias(const G4String &existingMaterialName, G4String alias)
Add alias to a material.
Definition: BDSMaterials.cc:963
BDSMaterials::possibleDuplicates
std::map< G4String, G4int > possibleDuplicates
Definition: BDSMaterials.hh:142
BDSMaterials::GetElement
G4Element * GetElement(const G4String &symbol) const
Get element by name.
Definition: BDSMaterials.cc:1156
BDSMaterials::PrintBasicMaterialMassFraction
void PrintBasicMaterialMassFraction(G4Material *material) const
Print mass fractions of consituents of a given material.
Definition: BDSMaterials.cc:1164
BDSMaterials::propertiesTables
std::vector< G4MaterialPropertiesTable * > propertiesTables
Material tables for storing pointers.
Definition: BDSMaterials.hh:146
BDSParser::Instance
static BDSParser * Instance()
Access method.
Definition: BDSParser.cc:28
BDSParser::GetAtoms
std::vector< GMAD::Atom > GetAtoms() const
Return the vector of atom objects.
Definition: BDSParser.hh:88
BDSParser::IsInitialised
static bool IsInitialised()
Returns if parser is initialised.
Definition: BDSParser.cc:49
BDSParser::GetOptions
const GMAD::Options & GetOptions() const
Return options.
Definition: BDSParser.hh:50
BDSParser::GetMaterials
std::vector< GMAD::Material > GetMaterials() const
Return material list.
Definition: BDSParser.hh:111
BDS::LowerCase
G4String LowerCase(const G4String &str)
Utility function to simplify lots of syntax changes for pedantic g4 changes.
BDS::StartsWith
G4bool StartsWith(const std::string &expression, const std::string &prefix)
Return true if a string "expression" starts with "prefix".
Definition: BDSUtilities.hh:266
Generated on Tue Sep 6 2022 10:41:16 for BDSIM by doxygen 1.9.3