MaterialPropertyLoader.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file MaterialPropertyLoader.cxx
//
/// \author  bjpjones@mit.edu
////////////////////////////////////////////////////////////////////////
// Class to set material properties for different materials in
// the detector. Currently mainly used to set optical properties
// for LAr and other optical components
//

// TODO verify the inclusion list
#include "larsim/LegacyLArG4/MaterialPropertyLoader.h"
#include "larcore/CoreUtils/ServiceUtil.h"
#include "lardata/DetectorInfoServices/LArPropertiesService.h"
#include "lardataalg/DetectorInfo/DetectorPropertiesData.h"

#include "Geant4/G4LogicalSkinSurface.hh"
#include "Geant4/G4LogicalVolume.hh"
#include "Geant4/G4LogicalVolumeStore.hh"
#include "Geant4/G4Material.hh"
#include "Geant4/G4MaterialPropertiesTable.hh"
#include "Geant4/G4OpticalSurface.hh"

#include "messagefacility/MessageLogger/MessageLogger.h"

namespace larg4 {

  //----------------------------------------------
  void
  MaterialPropertyLoader::SetMaterialProperty(std::string Material,
                                              std::string Property,
                                              std::map<double, double> PropertyVector,
                                              double Unit)
  {
    std::map<double, double> PropVectorWithUnit;
    for (std::map<double, double>::const_iterator it = PropertyVector.begin();
         it != PropertyVector.end();
         it++) {
      PropVectorWithUnit[it->first * CLHEP::eV] = it->second * Unit;
    }
    fPropertyList[Material][Property] = PropVectorWithUnit;
    // replace with MF_LOGDEBUG()
    mf::LogInfo("MaterialPropertyLoader") << "Added property " << Material << "  " << Property;
  }

  //----------------------------------------------
  void
  MaterialPropertyLoader::SetMaterialConstProperty(std::string Material,
                                                   std::string Property,
                                                   double PropertyValue,
                                                   double Unit)
  {
    fConstPropertyList[Material][Property] = PropertyValue * Unit;
    // replace with MF_LOGDEBUG()
    mf::LogInfo("MaterialPropertyLoader")
      << "Added const property " << Material << "  " << Property << " = " << PropertyValue;
  }

  //----------------------------------------------
  void
  MaterialPropertyLoader::SetBirksConstant(std::string Material, double PropertyValue, double Unit)
  {
    fBirksConstants[Material] = PropertyValue * Unit;
    // replace with MF_LOGDEBUG()
    mf::LogInfo("MaterialPropertyLoader") << "Set Birks constant " << Material;
  }

  //----------------------------------------------
  void
  MaterialPropertyLoader::UpdateGeometry(G4LogicalVolumeStore* lvs)
  {
    std::map<std::string, G4MaterialPropertiesTable*> MaterialTables;
    std::map<std::string, bool> MaterialsSet;

    // TODO replace console output with messagefacility output
    mf::LogInfo("MaterialPropertyLoader") << "UPDATING GEOMETRY";

    // Loop over each material with a property vector and create a new material table for it
    for (std::map<std::string, std::map<std::string, std::map<double, double>>>::const_iterator i =
           fPropertyList.begin();
         i != fPropertyList.end();
         i++) {
      std::string Material = i->first;
      MaterialsSet[Material] = true;
      MaterialTables[Material] = new G4MaterialPropertiesTable;
    }

    // Loop over each material with a const property,
    // if material table does not exist, create one
    for (std::map<std::string, std::map<std::string, double>>::const_iterator i =
           fConstPropertyList.begin();
         i != fConstPropertyList.end();
         i++) {
      std::string Material = i->first;
      if (!MaterialsSet[Material]) {
        MaterialsSet[Material] = true;
        MaterialTables[Material] = new G4MaterialPropertiesTable;
      }
    }

    // For each property vector, convert to an array of g4doubles and
    // feed to materials table Lots of firsts and seconds!  See annotation
    // in MaterialPropertyLoader.h to follow what each element is

    for (std::map<std::string, std::map<std::string, std::map<double, double>>>::const_iterator i =
           fPropertyList.begin();
         i != fPropertyList.end();
         i++) {
      std::string Material = i->first;
      for (std::map<std::string, std::map<double, double>>::const_iterator j = i->second.begin();
           j != i->second.end();
           j++) {
        std::string Property = j->first;
        std::vector<G4double> g4MomentumVector;
        std::vector<G4double> g4PropertyVector;

        for (std::map<double, double>::const_iterator k = j->second.begin(); k != j->second.end();
             k++) {
          g4MomentumVector.push_back(k->first);
          g4PropertyVector.push_back(k->second);
        }
        int NoOfElements = g4MomentumVector.size();
        MaterialTables[Material]->AddProperty(
          Property.c_str(), &g4MomentumVector[0], &g4PropertyVector[0], NoOfElements);
        // replace with mf::LogVerbatim()
        mf::LogInfo("MaterialPropertyLoader")
          << "Added property " << Property << " to material table " << Material;
      }
    }

    //Add each const property element
    for (std::map<std::string, std::map<std::string, double>>::const_iterator i =
           fConstPropertyList.begin();
         i != fConstPropertyList.end();
         i++) {
      std::string Material = i->first;
      for (std::map<std::string, double>::const_iterator j = i->second.begin();
           j != i->second.end();
           j++) {
        std::string Property = j->first;
        G4double PropertyValue = j->second;
        MaterialTables[Material]->AddConstProperty(Property.c_str(), PropertyValue);
        // replace with mf::LogVerbatim()
        mf::LogInfo("MaterialPropertyLoader")
          << "Added const property " << Property << " to material table " << Material;
      }
    }

    //Loop through geometry elements and apply relevant material table where materials match
    for (G4LogicalVolumeStore::iterator i = lvs->begin(); i != lvs->end(); ++i) {
      G4LogicalVolume* volume = (*i);
      G4Material* TheMaterial = volume->GetMaterial();
      std::string Material = TheMaterial->GetName();

      //
      // create reflective surfaces corresponding to the volumes made of some
      // selected materials
      //

      //--------------------------> FIXME <-----------------(parameters from fcl files(?))
      G4MaterialPropertyVector* PropertyPointer = 0;
      if (MaterialTables[Material])
        PropertyPointer = MaterialTables[Material]->GetProperty("REFLECTIVITY");

      if (Material == "Copper") {
        std::cout << "copper foil surface set " << volume->GetName() << std::endl;
        if (PropertyPointer) {
          std::cout << "defining Copper optical boundary " << std::endl;
          G4OpticalSurface* refl_opsurfc =
            new G4OpticalSurface("Surface copper", glisur, ground, dielectric_metal);
          refl_opsurfc->SetMaterialPropertiesTable(MaterialTables[Material]);
          refl_opsurfc->SetPolish(0.2);
          new G4LogicalSkinSurface("refl_surfacec", volume, refl_opsurfc);
        }
        else
          std::cout << "Warning: Copper surface in the geometry without REFLECTIVITY assigned"
                    << std::endl;
      }

      if (Material == "G10") {
        std::cout << "G10 surface set " << volume->GetName() << std::endl;
        if (PropertyPointer) {
          std::cout << "defining G10 optical boundary " << std::endl;
          G4OpticalSurface* refl_opsurfg =
            new G4OpticalSurface("g10 Surface", glisur, ground, dielectric_metal);
          refl_opsurfg->SetMaterialPropertiesTable(MaterialTables[Material]);
          refl_opsurfg->SetPolish(0.1);
          new G4LogicalSkinSurface("refl_surfaceg", volume, refl_opsurfg);
        }
        else
          std::cout << "Warning: G10 surface in the geometry without REFLECTIVITY assigned"
                    << std::endl;
      }

      if (Material == "vm2000") {
        std::cout << "vm2000 surface set " << volume->GetName() << std::endl;
        if (PropertyPointer) {
          std::cout << "defining vm2000 optical boundary " << std::endl;
          G4OpticalSurface* refl_opsurf = new G4OpticalSurface(
            "Reflector Surface", unified, groundfrontpainted, dielectric_dielectric);
          refl_opsurf->SetMaterialPropertiesTable(MaterialTables[Material]);
          G4double sigma_alpha = 0.8;
          refl_opsurf->SetSigmaAlpha(sigma_alpha);
          new G4LogicalSkinSurface("refl_surface", volume, refl_opsurf);
        }
        else
          std::cout << "Warning: vm2000 surface in the geometry without REFLECTIVITY assigned"
                    << std::endl;
      }
      if (Material == "ALUMINUM_Al") {
        std::cout << "ALUMINUM_Al surface set " << volume->GetName() << std::endl;
        if (PropertyPointer) {
          std::cout << "defining ALUMINUM_Al optical boundary " << std::endl;
          G4OpticalSurface* refl_opsurfs =
            new G4OpticalSurface("Surface Aluminum", glisur, ground, dielectric_metal);
          refl_opsurfs->SetMaterialPropertiesTable(MaterialTables[Material]);
          refl_opsurfs->SetPolish(0.5);
          new G4LogicalSkinSurface("refl_surfaces", volume, refl_opsurfs);
        }
        else
          std::cout << "Warning: ALUMINUM_Al surface in the geometry without REFLECTIVITY assigned"
                    << std::endl;
      }
      if (Material == "STEEL_STAINLESS_Fe7Cr2Ni") {
        std::cout << "STEEL_STAINLESS_Fe7Cr2Ni surface set " << volume->GetName() << std::endl;
        if (PropertyPointer) {
          std::cout << "defining STEEL_STAINLESS_Fe7Cr2Ni optical boundary " << std::endl;
          G4OpticalSurface* refl_opsurfs =
            new G4OpticalSurface("Surface Steel", glisur, ground, dielectric_metal);
          refl_opsurfs->SetMaterialPropertiesTable(MaterialTables[Material]);
          refl_opsurfs->SetPolish(0.5);
          new G4LogicalSkinSurface("refl_surfaces", volume, refl_opsurfs);
        }
        else
          std::cout << "Warning: STEEL_STAINLESS_Fe7Cr2Ni surface in the geometry without "
                       "REFLECTIVITY assigned"
                    << std::endl;
      }
      //-----------------------------------------------------------------------------

      //
      // apply the remaining material properties
      //
      for (std::map<std::string, G4MaterialPropertiesTable*>::const_iterator j =
             MaterialTables.begin();
           j != MaterialTables.end();
           j++) {
        if (Material == j->first) {
          TheMaterial->SetMaterialPropertiesTable(j->second);
          //Birks Constant, for some reason, must be set separately
          if (fBirksConstants[Material] != 0)
            TheMaterial->GetIonisation()->SetBirksConstant(fBirksConstants[Material]);
          volume->SetMaterial(TheMaterial);
        }
      }
    }
  }

  void
  MaterialPropertyLoader::SetReflectances(
    std::string /*Material*/,
    std::map<std::string, std::map<double, double>> Reflectances,
    std::map<std::string, std::map<double, double>> DiffuseFractions)
  {
    std::map<double, double> ReflectanceToStore;
    std::map<double, double> DiffuseToStore;

    for (std::map<std::string, std::map<double, double>>::const_iterator itMat =
           Reflectances.begin();
         itMat != Reflectances.end();
         ++itMat) {
      std::string ReflectancePropName = std::string("REFLECTANCE_") + itMat->first;
      ReflectanceToStore.clear();
      for (std::map<double, double>::const_iterator itEn = itMat->second.begin();
           itEn != itMat->second.end();
           ++itEn) {
        ReflectanceToStore[itEn->first] = itEn->second;
      }
      SetMaterialProperty("LAr", ReflectancePropName, ReflectanceToStore, 1);
    }

    for (std::map<std::string, std::map<double, double>>::const_iterator itMat =
           DiffuseFractions.begin();
         itMat != DiffuseFractions.end();
         ++itMat) {
      std::string DiffusePropName = std::string("DIFFUSE_REFLECTANCE_FRACTION_") + itMat->first;
      DiffuseToStore.clear();
      for (std::map<double, double>::const_iterator itEn = itMat->second.begin();
           itEn != itMat->second.end();
           ++itEn) {
        DiffuseToStore[itEn->first] = itEn->second;
      }
      SetMaterialProperty("LAr", DiffusePropName, DiffuseToStore, 1);
    }
  }

  void
  MaterialPropertyLoader::SetReflectances(
    std::map<std::string, std::map<double, double>> Reflectances)
  {
    std::map<double, double> ReflectanceToStore;

    for (std::map<std::string, std::map<double, double>>::const_iterator itMat =
           Reflectances.begin();
         itMat != Reflectances.end();
         ++itMat) {
      ReflectanceToStore.clear();
      for (std::map<double, double>::const_iterator itEn = itMat->second.begin();
           itEn != itMat->second.end();
           ++itEn) {
        ReflectanceToStore[itEn->first] = itEn->second;
      }
      SetMaterialProperty(itMat->first, "REFLECTIVITY", ReflectanceToStore, 1);
    }
  }

  void
  MaterialPropertyLoader::GetPropertiesFromServices(detinfo::DetectorPropertiesData const& detProp)
  {
    const detinfo::LArProperties* LarProp = lar::providerFrom<detinfo::LArPropertiesService>();

    // wavelength dependent quantities

    SetMaterialProperty("LAr", "FASTCOMPONENT", LarProp->FastScintSpectrum(), 1);
    SetMaterialProperty("LAr", "SLOWCOMPONENT", LarProp->SlowScintSpectrum(), 1);
    SetMaterialProperty("LAr", "RINDEX", LarProp->RIndexSpectrum(), 1);
    SetMaterialProperty("LAr", "ABSLENGTH", LarProp->AbsLengthSpectrum(), CLHEP::cm);
    SetMaterialProperty("LAr", "RAYLEIGH", LarProp->RayleighSpectrum(), CLHEP::cm);

    // scalar properties

    SetMaterialConstProperty("LAr",
                             "SCINTILLATIONYIELD",
                             LarProp->ScintYield(true),
                             1 / CLHEP::MeV); // true = scaled down by prescale in larproperties
    SetMaterialConstProperty("LAr", "RESOLUTIONSCALE", LarProp->ScintResolutionScale(), 1);
    SetMaterialConstProperty("LAr", "FASTTIMECONSTANT", LarProp->ScintFastTimeConst(), CLHEP::ns);
    SetMaterialConstProperty("LAr", "SLOWTIMECONSTANT", LarProp->ScintSlowTimeConst(), CLHEP::ns);
    SetMaterialConstProperty("LAr", "YIELDRATIO", LarProp->ScintYieldRatio(), 1);
    SetMaterialConstProperty("LAr", "ELECTRICFIELD", detProp.Efield(), CLHEP::kilovolt / CLHEP::cm);

    SetBirksConstant("LAr", LarProp->ScintBirksConstant(), CLHEP::cm / CLHEP::MeV);
    if (detProp.SimpleBoundary())
      SetReflectances(
        "LAr", LarProp->SurfaceReflectances(), LarProp->SurfaceReflectanceDiffuseFractions());
    else
      SetReflectances(LarProp->SurfaceReflectances());

    // If we are using scint by particle type, load these

    if (LarProp->ScintByParticleType()) {
      // true = scaled down by prescale in larproperties
      SetMaterialConstProperty(
        "LAr", "PROTONSCINTILLATIONYIELD", LarProp->ProtonScintYield(true), 1. / CLHEP::MeV);
      SetMaterialConstProperty("LAr", "PROTONYIELDRATIO", LarProp->ProtonScintYieldRatio(), 1.);
      SetMaterialConstProperty(
        "LAr", "MUONSCINTILLATIONYIELD", LarProp->MuonScintYield(true), 1. / CLHEP::MeV);
      SetMaterialConstProperty("LAr", "MUONYIELDRATIO", LarProp->MuonScintYieldRatio(), 1.);
      SetMaterialConstProperty(
        "LAr", "KAONSCINTILLATIONYIELD", LarProp->KaonScintYield(true), 1. / CLHEP::MeV);
      SetMaterialConstProperty("LAr", "KAONYIELDRATIO", LarProp->KaonScintYieldRatio(), 1.);
      SetMaterialConstProperty(
        "LAr", "PIONSCINTILLATIONYIELD", LarProp->PionScintYield(true), 1. / CLHEP::MeV);
      SetMaterialConstProperty("LAr", "PIONYIELDRATIO", LarProp->PionScintYieldRatio(), 1.);
      SetMaterialConstProperty(
        "LAr", "ELECTRONSCINTILLATIONYIELD", LarProp->ElectronScintYield(true), 1. / CLHEP::MeV);
      SetMaterialConstProperty("LAr", "ELECTRONYIELDRATIO", LarProp->ElectronScintYieldRatio(), 1.);
      SetMaterialConstProperty(
        "LAr", "ALPHASCINTILLATIONYIELD", LarProp->AlphaScintYield(true), 1. / CLHEP::MeV);
      SetMaterialConstProperty("LAr", "ALPHAYIELDRATIO", LarProp->AlphaScintYieldRatio(), 1.);
    }

    // If we are simulating the TPB load this

    if (LarProp->ExtraMatProperties()) {
      SetMaterialProperty("TPB", "RINDEX", LarProp->RIndexSpectrum(), 1);
      SetMaterialProperty("TPB", "WLSABSLENGTH", LarProp->TpbAbs(), CLHEP::m);
      SetMaterialProperty("TPB", "WLSCOMPONENT", LarProp->TpbEm(), 1);

      SetMaterialConstProperty("TPB", "WLSTIMECONSTANT", LarProp->TpbTimeConstant(), CLHEP::ns);
    }
  }

}