LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
MaterialPropertyLoader.cxx
Go to the documentation of this file.
1 //
6 // Class to set material properties for different materials in
7 // the detector. Currently mainly used to set optical properties
8 // for LAr and other optical components
9 //
10 
11 // TODO verify the inclusion list
16 
17 #include "Geant4/G4LogicalSkinSurface.hh"
18 #include "Geant4/G4LogicalVolume.hh"
19 #include "Geant4/G4LogicalVolumeStore.hh"
20 #include "Geant4/G4Material.hh"
21 #include "Geant4/G4MaterialPropertiesTable.hh"
22 #include "Geant4/G4OpticalSurface.hh"
23 
25 
26 namespace larg4 {
27 
28  //----------------------------------------------
30  std::string Property,
31  std::map<double, double> PropertyVector,
32  double Unit)
33  {
34  std::map<double, double> PropVectorWithUnit;
35  for (std::map<double, double>::const_iterator it = PropertyVector.begin();
36  it != PropertyVector.end();
37  it++) {
38  PropVectorWithUnit[it->first * CLHEP::eV] = it->second * Unit;
39  }
40  fPropertyList[Material][Property] = PropVectorWithUnit;
41  // replace with MF_LOGDEBUG()
42  mf::LogInfo("MaterialPropertyLoader") << "Added property " << Material << " " << Property;
43  }
44 
45  //----------------------------------------------
47  std::string Property,
48  double PropertyValue,
49  double Unit)
50  {
51  fConstPropertyList[Material][Property] = PropertyValue * Unit;
52  // replace with MF_LOGDEBUG()
53  mf::LogInfo("MaterialPropertyLoader")
54  << "Added const property " << Material << " " << Property << " = " << PropertyValue;
55  }
56 
57  //----------------------------------------------
58  void MaterialPropertyLoader::SetBirksConstant(std::string Material,
59  double PropertyValue,
60  double Unit)
61  {
62  fBirksConstants[Material] = PropertyValue * Unit;
63  // replace with MF_LOGDEBUG()
64  mf::LogInfo("MaterialPropertyLoader") << "Set Birks constant " << Material;
65  }
66 
67  //----------------------------------------------
68  void MaterialPropertyLoader::UpdateGeometry(G4LogicalVolumeStore* lvs)
69  {
70  std::map<std::string, G4MaterialPropertiesTable*> MaterialTables;
71  std::map<std::string, bool> MaterialsSet;
72 
73  // TODO replace console output with messagefacility output
74  mf::LogInfo("MaterialPropertyLoader") << "UPDATING GEOMETRY";
75 
76  // Loop over each material with a property vector and create a new material table for it
77  for (std::map<std::string, std::map<std::string, std::map<double, double>>>::const_iterator i =
78  fPropertyList.begin();
79  i != fPropertyList.end();
80  i++) {
81  std::string Material = i->first;
82  MaterialsSet[Material] = true;
83  MaterialTables[Material] = new G4MaterialPropertiesTable;
84  }
85 
86  // Loop over each material with a const property,
87  // if material table does not exist, create one
88  for (std::map<std::string, std::map<std::string, double>>::const_iterator i =
89  fConstPropertyList.begin();
90  i != fConstPropertyList.end();
91  i++) {
92  std::string Material = i->first;
93  if (!MaterialsSet[Material]) {
94  MaterialsSet[Material] = true;
95  MaterialTables[Material] = new G4MaterialPropertiesTable;
96  }
97  }
98 
99  // For each property vector, convert to an array of g4doubles and
100  // feed to materials table Lots of firsts and seconds! See annotation
101  // in MaterialPropertyLoader.h to follow what each element is
102 
103  for (std::map<std::string, std::map<std::string, std::map<double, double>>>::const_iterator i =
104  fPropertyList.begin();
105  i != fPropertyList.end();
106  i++) {
107  std::string Material = i->first;
108  for (std::map<std::string, std::map<double, double>>::const_iterator j = i->second.begin();
109  j != i->second.end();
110  j++) {
111  std::string Property = j->first;
112  std::vector<G4double> g4MomentumVector;
113  std::vector<G4double> g4PropertyVector;
114 
115  for (std::map<double, double>::const_iterator k = j->second.begin(); k != j->second.end();
116  k++) {
117  g4MomentumVector.push_back(k->first);
118  g4PropertyVector.push_back(k->second);
119  }
120  int NoOfElements = g4MomentumVector.size();
121  MaterialTables[Material]->AddProperty(
122  Property.c_str(), &g4MomentumVector[0], &g4PropertyVector[0], NoOfElements);
123  // replace with mf::LogVerbatim()
124  mf::LogInfo("MaterialPropertyLoader")
125  << "Added property " << Property << " to material table " << Material;
126  }
127  }
128 
129  //Add each const property element
130  for (std::map<std::string, std::map<std::string, double>>::const_iterator i =
131  fConstPropertyList.begin();
132  i != fConstPropertyList.end();
133  i++) {
134  std::string Material = i->first;
135  for (std::map<std::string, double>::const_iterator j = i->second.begin();
136  j != i->second.end();
137  j++) {
138  std::string Property = j->first;
139  G4double PropertyValue = j->second;
140  MaterialTables[Material]->AddConstProperty(Property.c_str(), PropertyValue);
141  // replace with mf::LogVerbatim()
142  mf::LogInfo("MaterialPropertyLoader")
143  << "Added const property " << Property << " to material table " << Material;
144  }
145  }
146 
147  //Loop through geometry elements and apply relevant material table where materials match
148  for (G4LogicalVolumeStore::iterator i = lvs->begin(); i != lvs->end(); ++i) {
149  G4LogicalVolume* volume = (*i);
150  G4Material* TheMaterial = volume->GetMaterial();
151  std::string Material = TheMaterial->GetName();
152 
153  //
154  // create reflective surfaces corresponding to the volumes made of some
155  // selected materials
156  //
157 
158  //--------------------------> FIXME <-----------------(parameters from fcl files(?))
159  G4MaterialPropertyVector* PropertyPointer = 0;
160  if (MaterialTables[Material])
161  PropertyPointer = MaterialTables[Material]->GetProperty("REFLECTIVITY");
162 
163  if (Material == "Copper") {
164  std::cout << "copper foil surface set " << volume->GetName() << std::endl;
165  if (PropertyPointer) {
166  std::cout << "defining Copper optical boundary " << std::endl;
167  G4OpticalSurface* refl_opsurfc =
168  new G4OpticalSurface("Surface copper", glisur, ground, dielectric_metal);
169  refl_opsurfc->SetMaterialPropertiesTable(MaterialTables[Material]);
170  refl_opsurfc->SetPolish(0.2);
171  new G4LogicalSkinSurface("refl_surfacec", volume, refl_opsurfc);
172  }
173  else
174  std::cout << "Warning: Copper surface in the geometry without REFLECTIVITY assigned"
175  << std::endl;
176  }
177 
178  if (Material == "G10") {
179  std::cout << "G10 surface set " << volume->GetName() << std::endl;
180  if (PropertyPointer) {
181  std::cout << "defining G10 optical boundary " << std::endl;
182  G4OpticalSurface* refl_opsurfg =
183  new G4OpticalSurface("g10 Surface", glisur, ground, dielectric_metal);
184  refl_opsurfg->SetMaterialPropertiesTable(MaterialTables[Material]);
185  refl_opsurfg->SetPolish(0.1);
186  new G4LogicalSkinSurface("refl_surfaceg", volume, refl_opsurfg);
187  }
188  else
189  std::cout << "Warning: G10 surface in the geometry without REFLECTIVITY assigned"
190  << std::endl;
191  }
192 
193  if (Material == "vm2000") {
194  std::cout << "vm2000 surface set " << volume->GetName() << std::endl;
195  if (PropertyPointer) {
196  std::cout << "defining vm2000 optical boundary " << std::endl;
197  G4OpticalSurface* refl_opsurf = new G4OpticalSurface(
198  "Reflector Surface", unified, groundfrontpainted, dielectric_dielectric);
199  refl_opsurf->SetMaterialPropertiesTable(MaterialTables[Material]);
200  G4double sigma_alpha = 0.8;
201  refl_opsurf->SetSigmaAlpha(sigma_alpha);
202  new G4LogicalSkinSurface("refl_surface", volume, refl_opsurf);
203  }
204  else
205  std::cout << "Warning: vm2000 surface in the geometry without REFLECTIVITY assigned"
206  << std::endl;
207  }
208  if (Material == "ALUMINUM_Al") {
209  std::cout << "ALUMINUM_Al surface set " << volume->GetName() << std::endl;
210  if (PropertyPointer) {
211  std::cout << "defining ALUMINUM_Al optical boundary " << std::endl;
212  G4OpticalSurface* refl_opsurfs =
213  new G4OpticalSurface("Surface Aluminum", glisur, ground, dielectric_metal);
214  refl_opsurfs->SetMaterialPropertiesTable(MaterialTables[Material]);
215  refl_opsurfs->SetPolish(0.5);
216  new G4LogicalSkinSurface("refl_surfaces", volume, refl_opsurfs);
217  }
218  else
219  std::cout << "Warning: ALUMINUM_Al surface in the geometry without REFLECTIVITY assigned"
220  << std::endl;
221  }
222  if (Material == "STEEL_STAINLESS_Fe7Cr2Ni") {
223  std::cout << "STEEL_STAINLESS_Fe7Cr2Ni surface set " << volume->GetName() << std::endl;
224  if (PropertyPointer) {
225  std::cout << "defining STEEL_STAINLESS_Fe7Cr2Ni optical boundary " << std::endl;
226  G4OpticalSurface* refl_opsurfs =
227  new G4OpticalSurface("Surface Steel", glisur, ground, dielectric_metal);
228  refl_opsurfs->SetMaterialPropertiesTable(MaterialTables[Material]);
229  refl_opsurfs->SetPolish(0.5);
230  new G4LogicalSkinSurface("refl_surfaces", volume, refl_opsurfs);
231  }
232  else
233  std::cout << "Warning: STEEL_STAINLESS_Fe7Cr2Ni surface in the geometry without "
234  "REFLECTIVITY assigned"
235  << std::endl;
236  }
237  //-----------------------------------------------------------------------------
238 
239  //
240  // apply the remaining material properties
241  //
243  MaterialTables.begin();
244  j != MaterialTables.end();
245  j++) {
246  if (Material == j->first) {
247  TheMaterial->SetMaterialPropertiesTable(j->second);
248  //Birks Constant, for some reason, must be set separately
249  if (fBirksConstants[Material] != 0)
250  TheMaterial->GetIonisation()->SetBirksConstant(fBirksConstants[Material]);
251  volume->SetMaterial(TheMaterial);
252  }
253  }
254  }
255  }
256 
258  std::string /*Material*/,
259  std::map<std::string, std::map<double, double>> Reflectances,
260  std::map<std::string, std::map<double, double>> DiffuseFractions)
261  {
262  std::map<double, double> ReflectanceToStore;
263  std::map<double, double> DiffuseToStore;
264 
265  for (std::map<std::string, std::map<double, double>>::const_iterator itMat =
266  Reflectances.begin();
267  itMat != Reflectances.end();
268  ++itMat) {
269  std::string ReflectancePropName = std::string("REFLECTANCE_") + itMat->first;
270  ReflectanceToStore.clear();
271  for (std::map<double, double>::const_iterator itEn = itMat->second.begin();
272  itEn != itMat->second.end();
273  ++itEn) {
274  ReflectanceToStore[itEn->first] = itEn->second;
275  }
276  SetMaterialProperty("LAr", ReflectancePropName, ReflectanceToStore, 1);
277  }
278 
279  for (std::map<std::string, std::map<double, double>>::const_iterator itMat =
280  DiffuseFractions.begin();
281  itMat != DiffuseFractions.end();
282  ++itMat) {
283  std::string DiffusePropName = std::string("DIFFUSE_REFLECTANCE_FRACTION_") + itMat->first;
284  DiffuseToStore.clear();
285  for (std::map<double, double>::const_iterator itEn = itMat->second.begin();
286  itEn != itMat->second.end();
287  ++itEn) {
288  DiffuseToStore[itEn->first] = itEn->second;
289  }
290  SetMaterialProperty("LAr", DiffusePropName, DiffuseToStore, 1);
291  }
292  }
293 
295  std::map<std::string, std::map<double, double>> Reflectances)
296  {
297  std::map<double, double> ReflectanceToStore;
298 
299  for (std::map<std::string, std::map<double, double>>::const_iterator itMat =
300  Reflectances.begin();
301  itMat != Reflectances.end();
302  ++itMat) {
303  ReflectanceToStore.clear();
304  for (std::map<double, double>::const_iterator itEn = itMat->second.begin();
305  itEn != itMat->second.end();
306  ++itEn) {
307  ReflectanceToStore[itEn->first] = itEn->second;
308  }
309  SetMaterialProperty(itMat->first, "REFLECTIVITY", ReflectanceToStore, 1);
310  }
311  }
312 
314  detinfo::DetectorPropertiesData const& detProp)
315  {
316  const detinfo::LArProperties* LarProp = lar::providerFrom<detinfo::LArPropertiesService>();
317 
318  // wavelength dependent quantities
319 
320  SetMaterialProperty("LAr", "FASTCOMPONENT", LarProp->FastScintSpectrum(), 1);
321  SetMaterialProperty("LAr", "SLOWCOMPONENT", LarProp->SlowScintSpectrum(), 1);
322  SetMaterialProperty("LAr", "RINDEX", LarProp->RIndexSpectrum(), 1);
323  SetMaterialProperty("LAr", "ABSLENGTH", LarProp->AbsLengthSpectrum(), CLHEP::cm);
324  SetMaterialProperty("LAr", "RAYLEIGH", LarProp->RayleighSpectrum(), CLHEP::cm);
325 
326  // scalar properties
327 
329  "SCINTILLATIONYIELD",
330  LarProp->ScintYield(true),
331  1 / CLHEP::MeV); // true = scaled down by prescale in larproperties
332  SetMaterialConstProperty("LAr", "RESOLUTIONSCALE", LarProp->ScintResolutionScale(), 1);
333  SetMaterialConstProperty("LAr", "FASTTIMECONSTANT", LarProp->ScintFastTimeConst(), CLHEP::ns);
334  SetMaterialConstProperty("LAr", "SLOWTIMECONSTANT", LarProp->ScintSlowTimeConst(), CLHEP::ns);
335  SetMaterialConstProperty("LAr", "YIELDRATIO", LarProp->ScintYieldRatio(), 1);
336  SetMaterialConstProperty("LAr", "ELECTRICFIELD", detProp.Efield(), CLHEP::kilovolt / CLHEP::cm);
337 
338  SetBirksConstant("LAr", LarProp->ScintBirksConstant(), CLHEP::cm / CLHEP::MeV);
339  if (detProp.SimpleBoundary())
341  "LAr", LarProp->SurfaceReflectances(), LarProp->SurfaceReflectanceDiffuseFractions());
342  else
344 
345  // If we are using scint by particle type, load these
346 
347  if (LarProp->ScintByParticleType()) {
348  // true = scaled down by prescale in larproperties
350  "LAr", "PROTONSCINTILLATIONYIELD", LarProp->ProtonScintYield(true), 1. / CLHEP::MeV);
351  SetMaterialConstProperty("LAr", "PROTONYIELDRATIO", LarProp->ProtonScintYieldRatio(), 1.);
353  "LAr", "MUONSCINTILLATIONYIELD", LarProp->MuonScintYield(true), 1. / CLHEP::MeV);
354  SetMaterialConstProperty("LAr", "MUONYIELDRATIO", LarProp->MuonScintYieldRatio(), 1.);
356  "LAr", "KAONSCINTILLATIONYIELD", LarProp->KaonScintYield(true), 1. / CLHEP::MeV);
357  SetMaterialConstProperty("LAr", "KAONYIELDRATIO", LarProp->KaonScintYieldRatio(), 1.);
359  "LAr", "PIONSCINTILLATIONYIELD", LarProp->PionScintYield(true), 1. / CLHEP::MeV);
360  SetMaterialConstProperty("LAr", "PIONYIELDRATIO", LarProp->PionScintYieldRatio(), 1.);
362  "LAr", "ELECTRONSCINTILLATIONYIELD", LarProp->ElectronScintYield(true), 1. / CLHEP::MeV);
363  SetMaterialConstProperty("LAr", "ELECTRONYIELDRATIO", LarProp->ElectronScintYieldRatio(), 1.);
365  "LAr", "ALPHASCINTILLATIONYIELD", LarProp->AlphaScintYield(true), 1. / CLHEP::MeV);
366  SetMaterialConstProperty("LAr", "ALPHAYIELDRATIO", LarProp->AlphaScintYieldRatio(), 1.);
367  }
368 
369  // If we are simulating the TPB load this
370 
371  if (LarProp->ExtraMatProperties()) {
372  SetMaterialProperty("TPB", "RINDEX", LarProp->RIndexSpectrum(), 1);
373  SetMaterialProperty("TPB", "WLSABSLENGTH", LarProp->TpbAbs(), CLHEP::m);
374  SetMaterialProperty("TPB", "WLSCOMPONENT", LarProp->TpbEm(), 1);
375 
376  SetMaterialConstProperty("TPB", "WLSTIMECONSTANT", LarProp->TpbTimeConstant(), CLHEP::ns);
377  }
378  }
379 
380 }
intermediate_table::iterator iterator
virtual double ElectronScintYieldRatio() const =0
virtual double ScintYieldRatio() const =0
virtual double AlphaScintYield(bool prescale=false) const =0
virtual std::map< double, double > AbsLengthSpectrum() const =0
Utilities related to art service access.
void SetMaterialProperty(std::string Material, std::string Property, std::map< double, double > Values, double Unit)
Stores the specified emergy-dependent material property.
virtual std::map< double, double > RayleighSpectrum() const =0
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
virtual std::map< double, double > RIndexSpectrum() const =0
virtual double AlphaScintYieldRatio() const =0
void GetPropertiesFromServices(detinfo::DetectorPropertiesData const &detProp)
Imports properties from LArSoft services.
virtual double ScintSlowTimeConst() const =0
std::map< std::string, std::map< std::string, std::map< double, double > > > fPropertyList
virtual double ScintFastTimeConst() const =0
Geant4 interface.
intermediate_table::const_iterator const_iterator
virtual double ProtonScintYieldRatio() const =0
void SetMaterialConstProperty(std::string Material, std::string Property, double Value, double Unit)
Stores the specified material property.
kilovolt_as<> kilovolt
Type of potential stored in kilovolt, in double precision.
virtual double ScintResolutionScale() const =0
double Efield(unsigned int planegap=0) const
kV/cm
virtual double ProtonScintYield(bool prescale=false) const =0
virtual double PionScintYieldRatio() const =0
std::map< std::string, std::map< std::string, double > > fConstPropertyList
virtual bool ExtraMatProperties() const =0
virtual double ScintBirksConstant() const =0
virtual double ElectronScintYield(bool prescale=false) const =0
virtual double PionScintYield(bool prescale=false) const =0
virtual std::map< double, double > TpbEm() const =0
virtual double TpbTimeConstant() const =0
virtual std::map< std::string, std::map< double, double > > SurfaceReflectances() const =0
void UpdateGeometry(G4LogicalVolumeStore *lvs)
Updates the material properties with the collected values.
virtual std::map< double, double > SlowScintSpectrum() const =0
virtual double MuonScintYield(bool prescale=false) const =0
virtual double KaonScintYield(bool prescale=false) const =0
virtual std::map< double, double > TpbAbs() const =0
virtual double ScintYield(bool prescale=false) const =0
virtual double MuonScintYieldRatio() const =0
void SetReflectances(std::string, std::map< std::string, std::map< double, double >>, std::map< std::string, std::map< double, double >>)
virtual std::map< std::string, std::map< double, double > > SurfaceReflectanceDiffuseFractions() const =0
void SetBirksConstant(std::string, double, double)
virtual double KaonScintYieldRatio() const =0
std::map< std::string, double > fBirksConstants
virtual bool ScintByParticleType() const =0
virtual std::map< double, double > FastScintSpectrum() const =0