G4StationConstruction.cc

Go to the documentation of this file.

#include "G4StationConstruction.h"

#include <fwk/RunController.h>
#include <fwk/CentralConfig.h>

#include <utl/CoordinateSystem.h>
#include <utl/Point.h>
#include <utl/TabulatedFunction.h>
#include <utl/Reader.h>
#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/Math.h>

#include <G4PhysicalConstants.hh>
#include <G4SystemOfUnits.hh>
#include <G4PVReplica.hh>

#include <vector>
#include <sstream>
#include <algorithm>

using sdet::SDetector;
using utl::CoordinateSystem;
using utl::CoordinateSystemPtr;
using utl::Point;
using utl::TabulatedFunction;
using utl::Branch;

using namespace fwk;

using namespace std;


namespace G4StationSimulatorOG {

  // declare variables that static members above refer to.
  G4ThreeVector G4StationConstruction::fgTankCenter;

  double G4StationConstruction::fgTankRadius = 0;
  double G4StationConstruction::fgTankHalfHeight = 0;
  double G4StationConstruction::fgTankThickness = 0;

  double G4StationConstruction::fgHeightz = 0;
  double G4StationConstruction::fgPmtRmax = 0;
  double G4StationConstruction::fgPmtRzmax = 0;
  double G4StationConstruction::fgInterfaceRmax = 0;
  double G4StationConstruction::fgInterfaceRzmax = 0;
  double G4StationConstruction::fgDomeRmax = 0;
  double G4StationConstruction::fgDomeRzmax = 0;

  double G4StationConstruction::fPmtRmin_SPMT = 0;
  double G4StationConstruction::fPmtRmax_SPMT = 0;
  double G4StationConstruction::fInterfaceRadius_SPMT = 0;
  double G4StationConstruction::fDomeRadius_SPMT = 0;
  double G4StationConstruction::fGlassThickness_SPMT = 0;
  double G4StationConstruction::fInterfaceThickness_SPMT = 0;
  double G4StationConstruction::fDomeThickness_SPMT = 0;

  TabulatedFunction G4StationConstruction::fgInterfaceRINDEX;
  TabulatedFunction G4StationConstruction::fgPmtdomeRINDEX;
  TabulatedFunction G4StationConstruction::fgInterfaceABSORPTION;
  TabulatedFunction G4StationConstruction::fgPmtdomeABSORPTION;

  double G4StationConstruction::fgSIGMA_ALPHA = 0;

  TabulatedFunction G4StationConstruction::fgLinerREFLECTIVITY;
  TabulatedFunction G4StationConstruction::fgLinerSPECULARLOBECONSTANT;
  TabulatedFunction G4StationConstruction::fgLinerSPECULARSPIKECONSTANT;
  TabulatedFunction G4StationConstruction::fgLinerBACKSCATTERCONSTANT;
  TabulatedFunction G4StationConstruction::fgWaterABSORPTION;

  //

  const double kDistanceToG4 = CLHEP::m / utl::m;
  const double kEnergyToG4 = CLHEP::eV / utl::eV;


  // helpers

  template<class T>
  inline
  void
  Scale(T& v, const double factor)
  {
    for (auto& x : v)
      x *= factor;
  }


  template<class T>
  inline
  void
  Scale(const T& b, const T& e, const double factor)
  {
    for (auto it = b; it != e; ++it)
      *it *= factor;
  }


  template<class T>
  inline
  void
  Renew(T* (&p), T* const val)
  {
    delete p;
    p = val;
  }


  // TODO this function needs to be removed; suitable member pointers replaced with local vars
  // Needed for replication of identical logical volumes places in
  // different postions and, therefore, generating different physical volumes
  template<class T>
  inline
  void
  NotRenew(T* (&p), T* const val)
  {
    p = val;
  }


  inline
  G4ThreeVector
  Convert(const utl::Point p, const utl::CoordinateSystemPtr cs)
  {
    const utl::Point::Triple t = p.GetCoordinates(cs);
    return kDistanceToG4 * G4ThreeVector(t.get<0>(), t.get<1>(), t.get<2>());
  }


  template<class Material>
  inline
  void
  AddProperty(Material* const m, const char* const name, const TabulatedFunction& f)
  {
    // const-correctness of G4 is broken, AddProperty actually does not modify the arrays
    m->AddProperty(name, (double*)&f.XFront(), (double*)&f.YFront(), f.GetNPoints());
  }


  G4StationConstruction::G4StationConstruction(const double simRadius, const bool umdEnabled,
                                               const double scintYield, const bool martaEnabled) :
    fG4StationSimulator(
      dynamic_cast<G4StationSimulator&>(RunController::GetInstance().GetModule("G4StationSimulatorOG"))
    )
  {
    /* Set hardware version. This is used in determining whether or not
      to set xmlparameters and construct and register hardware-version-dependent volumes (e.g.
      those relating to the small PMT and the scintillator) */

    const auto& dStation = G4StationSimulator::fgCurrent.GetDetectorStation();
    fHasSmallPMT = dStation.HasSmallPMT();
    fHasScintillator = dStation.HasScintillator();
    // Needed to be consistent with CachedShowerRegenerator area of injection (NB: variable loaded with units)
    fSimulationRadius = simRadius;
    fUMDEnabled = umdEnabled;
    fScintYield = scintYield;
    fMARTAEnabled = martaEnabled;

    SetRequiredParameters();
    SetXMLParameters();
    SetDetectorParameters();

    // disable solar panel, electronics box for now:

    fGroundEnable = fUMDEnabled;

    fNistManager = G4NistManager::Instance();
  }


  void
  G4StationConstruction::SetRequiredParameters()
  {
    // The "experimental hall" in which the detector is contained:
    const double exp_hall_side_side = fSimulationRadius;

    // The ground beneath the tank
    fGroundThickness = 3*CLHEP::m;

    fExpHall_x = exp_hall_side_side;
    fExpHall_y = exp_hall_side_side;
    fExpHall_z = fGroundThickness;

    // The position of the center of the tank in the local Geant4 coordinate system
    fTankPos_x = 0;
    fTankPos_y = 0;
    fTankPos_z = 0;
  }


  void
  G4StationConstruction::SetDetectorParameters()
  {
    const auto& dStation = G4StationSimulator::fgCurrent.GetDetectorStation();
    // Set tank properties
    fgTankRadius = dStation.GetRadius() * kDistanceToG4;
    fgTankHalfHeight = 0.5 * dStation.GetHeight() * kDistanceToG4;
    fgTankThickness = dStation.GetThickness() * kDistanceToG4;

    fgTankCenter =
      G4ThreeVector(fTankPos_x, fTankPos_y, fTankPos_z) +
      G4ThreeVector(0, 0, fgTankHalfHeight + fgTankThickness);

    // Set pmt position (relative to the center of the tank)

    const auto cs = dStation.GetLocalCoordinateSystem();

    for (auto& p : fPMTs)
      p.Clear();
    fPMTs.clear();
    for (const auto& dpmt : dStation.PMTsRange()) {
      PMTConstruction pmt;
      pmt.fPMTPosition = Convert(dpmt.GetPosition(), cs);
      pmt.fId = dpmt.GetId();
      pmt.fPMTType = dpmt.GetType();
      fPMTs.push_back(pmt);
    }

    // Set pmt optical properties
    // For now, assume that all pmts have the same properties

    fPmtfaceRINDEX = dStation.GetPMT(1).GetFaceRefractionIndex();
    Scale(fPmtfaceRINDEX.XBegin(), fPmtfaceRINDEX.XEnd(), kEnergyToG4);

    fgInterfaceRINDEX = dStation.GetPMT(1).GetRTVRefractionIndex();
    Scale(fgInterfaceRINDEX.XBegin(), fgInterfaceRINDEX.XEnd(), kEnergyToG4);

    fgPmtdomeRINDEX = dStation.GetPMT(1).GetDomeRefractionIndex();
    Scale(fgPmtdomeRINDEX.XBegin(), fgPmtdomeRINDEX.XEnd(), kEnergyToG4);

    // These parameters specify the pmt geometry (already converted to G4 units)

    fPmtRmin = fFaceRadius - fGlassThickness;
    fgPmtRmax = fFaceRadius;
    fPmtRzmin = fFaceRadiusz - fGlassThickness;
    fgPmtRzmax = fFaceRadiusz;
    fInterfaceRmin = fFaceRadius;
    fInterfaceRzmin = fFaceRadiusz;
    fgInterfaceRmax = fFaceRadius + fInterfaceThickness;
    fgInterfaceRzmax = fFaceRadiusz + fInterfaceThickness;
    fDomeRmin = fFaceRadius + fInterfaceThickness;
    fDomeRzmin = fFaceRadiusz + fInterfaceThickness;
    fgDomeRmax = fFaceRadius + fInterfaceThickness + fDomeThickness;
    fgDomeRzmax = fFaceRadiusz + fInterfaceThickness + fDomeThickness;

    G4double alpha = 0;
    // Taking into account that the radius of the photocathode is fFaceActiveRadius
    if (fFaceActiveRadius >= 114*mm) {
      fFaceActiveRadius = 113.926112814*mm;
      fgHeightz = 0;
    } else if (fFaceActiveRadius > 97.27*mm) {
      alpha = acos(1/62.*(fFaceActiveRadius/mm - 133*sin(47*deg) + 62*cos(43*deg)));
      fgHeightz = 62*mm*sin(alpha);
    } else {
      alpha = asin(fFaceActiveRadius/mm/133);
      fgHeightz = (133*cos(alpha) - (133 - 62)*cos(47*deg))*mm;
    }

    // Set liner properties

    fgSIGMA_ALPHA = dStation.GetLinerSigmaAlpha();

    fgLinerREFLECTIVITY = dStation.GetLinerReflectivity();
    Scale(fgLinerREFLECTIVITY.XBegin(), fgLinerREFLECTIVITY.XEnd(), kEnergyToG4);

    fgLinerSPECULARLOBECONSTANT = dStation.GetLinerSpecularLobe();
    Scale(fgLinerSPECULARLOBECONSTANT.XBegin(), fgLinerSPECULARLOBECONSTANT.XEnd(), kEnergyToG4);

    fgLinerSPECULARSPIKECONSTANT = dStation.GetLinerSpecularSpike();
    Scale(fgLinerSPECULARSPIKECONSTANT.XBegin(), fgLinerSPECULARSPIKECONSTANT.XEnd(), kEnergyToG4);

    // Set water properties

    fgWaterABSORPTION = dStation.GetWaterAbsorptionLength();
    Scale(fgWaterABSORPTION.XBegin(), fgWaterABSORPTION.XEnd(), kEnergyToG4);
    Scale(fgWaterABSORPTION.YBegin(), fgWaterABSORPTION.YEnd(), kDistanceToG4);

    fWaterRINDEX = dStation.GetWaterRefractionIndex();
    Scale(fWaterRINDEX.XBegin(), fWaterRINDEX.XEnd(), kEnergyToG4);

    // hardware-version-specific detector parameters

    if (dStation.HasScintillator()) {
      const sdet::Scintillator& scintillator = dStation.GetScintillator();
      // scintillator bars
      fScintillatorBarLength = scintillator.GetBarLength() * kDistanceToG4;
      fScintillatorBarWidth = scintillator.GetBarWidth() * kDistanceToG4;
      fScintillatorBarThickness = scintillator.GetBarThickness() * kDistanceToG4;
      fNScintillatorBars = scintillator.GetNBars();
      fScintillatorGap = scintillator.GetGap() * kDistanceToG4;
      // polystyrene housing
      fScintillatorHousingLength = scintillator.GetHousingLength() * kDistanceToG4;
      fScintillatorHousingWidth = scintillator.GetHousingWidth() * kDistanceToG4;
      fScintillatorHousingThickness = scintillator.GetHousingThickness() * kDistanceToG4;
      // aluminum casing
      fScintillatorCasingThickness = scintillator.GetCasingThickness() * kDistanceToG4;
      // sandwich panel
      fScintillatorSandwichPanelThickness = scintillator.GetSandwichPanelThickness() * kDistanceToG4;
      // roof
      fScintillatorRoofOffset = scintillator.GetRoofOffset() * kDistanceToG4;
      fScintillatorRoofThickness = scintillator.GetRoofThickness() * kDistanceToG4;
      // positioning
      fScintillatorPosition = Convert(scintillator.GetPosition(), cs);
    }

    // These parameters specify the small PMT geometry ( already converted to G4 units )

    fPmtRmin_SPMT = fFaceActiveRadius_SPMT;
    fPmtRmax_SPMT = fFaceRadius_SPMT;
    fInterfaceRadius_SPMT = fWindowRadius_SPMT;
    fDomeRadius_SPMT = fWindowRadius_SPMT;

    fPmtZ_SPMT = fDomeThickness_SPMT + fInterfaceThickness_SPMT + fGlassThickness_SPMT/2;
    fInterfaceZ_SPMT = fDomeThickness_SPMT + fInterfaceThickness_SPMT/2;
    fDomeZ_SPMT = fDomeThickness_SPMT/2;

    // Creates Amiga underground primitives
    if (fUMDEnabled) {
      if (!dStation.ExistsAssociatedCounter()) {
        ostringstream err;
        err << "Trying to simulate the UMD in a WCD (Id "<< dStation.GetId()
            << ") that has not an underground AMIGA scintillator associated!";
        //throw utl::NonExistentComponentException(err.str());
        WARNING(err);
      }
      SetUMDParameters();
    }

    // Creates MARTA underground primitives
    if (fMARTAEnabled) {
      SetMARTAParameters();

      // Raises the tank by the height of the precast
      const double martaSupportHeight =
        fTankSupportTopSlabDimensions[2] +
        fTankSupportOuterFootDimensions[2] +
        fTankSupportOuterFootBaseDimensions[2];

      fTankPos_z += martaSupportHeight;

      fgTankCenter =
        G4ThreeVector(fTankPos_x, fTankPos_y, fTankPos_z) +
        G4ThreeVector(0, 0, fgTankHalfHeight + fgTankThickness);
    }
  }


  // Defines area kinds for UMD modules
  void
  G4StationConstruction::SetUMDParameters()
  {
    const auto& dStation = G4StationSimulator::fgCurrent.GetDetectorStation();
    // Get the muon detector
    const mdet::MDetector& mDetector = det::Detector::GetInstance().GetMDetector();
    const mdet::Counter& cmDet = mDetector.GetCounter(dStation.GetAssociatedCounterId());

    fUMDThicknessCasing = 1*CLHEP::cm;
    // Largest length on manifold is 940mm, shortest is 530mm
    // Therefore the largest relative difference is 410mm
    fUMDManifoldL = 1*CLHEP::m;

    UMDScintiBox sbox;
    UMDFiberTube ftub;
    // Loop over muon modules (may only differs in bar lengths)
    for (auto mmDet = cmDet.ModulesBegin(); mmDet != cmDet.ModulesEnd(); ++mmDet) {

       const mdet::Scintillator& s = *mmDet->ScintillatorsBegin();
       const mdet::Fiber& f = mmDet->GetFiberFor(s);

       double l = 0;
       f.VisitShape(ftub).GetDimensions(l, fUMDFiberRadius);
       s.VisitShape(sbox).GetDimensions(l, fUMDScintsW, fUMDScintsH);

       l *= kDistanceToG4;
       fUMDScintsW *= kDistanceToG4;
       fUMDScintsH *= kDistanceToG4;
       fUMDScints = mmDet->GetNumberOfScintillators();
       fUMDFiberRadius *= kDistanceToG4;

       mdet::Module::AreaKind area = mmDet->GetAreaKind();
       if (!fModAreaLenghts.count(area))
         fModAreaLenghts[area] = l;
    }
  }


  void
  G4StationConstruction::SetMARTAParameters()
  {
    const auto& dStation = G4StationSimulator::fgCurrent.GetDetectorStation();
    const int sId = dStation.GetId();
    const cdet::Station& station = det::Detector::GetInstance().GetCDetector().GetStation(sId);

    {
      vector<double> t;

      t = station.GetTankSupportTopSlabDimensions();
      fTankSupportTopSlabDimensions = G4ThreeVector(t[0], t[1], t[2]) * kDistanceToG4;

      t = station.GetTankSupportCentralFootDimensions();
      fTankSupportCentralFootDimensions = G4ThreeVector(t[0], t[1], t[2]) * kDistanceToG4;

      t = station.GetTankSupportCentralFootBaseDimensions();
      fTankSupportCentralFootBaseDimensions = G4ThreeVector(t[0], t[1], t[2]) * kDistanceToG4;

      t = station.GetTankSupportOuterFootDimensions();
      fTankSupportOuterFootDimensions = G4ThreeVector(t[0], t[1], t[2]) * kDistanceToG4;

      t = station.GetTankSupportOuterFootBaseDimensions();
      fTankSupportOuterFootBaseDimensions = G4ThreeVector(t[0], t[1], t[2]) * kDistanceToG4;

      t = station.GetRPCHousingInnerDimensions();
      fAlBoxInnerDimensions = G4ThreeVector(t[0], t[1], t[2]) * kDistanceToG4;

      t = station.GetRPCHousingThickness();
      fAlBoxThickness = G4ThreeVector(t[0], t[1], t[2]) * kDistanceToG4;
    }

    fTankSupportOuterFootDistanceToCenter = station.GetTankSupportOuterFootDistanceToCenter() * kDistanceToG4;

    fPCBThickness = 1*CLHEP::mm;
    fSpacerThickness = 2*CLHEP::cm;

    fRPCsizeX = station.GetRPCSizeX() * kDistanceToG4;
    fRPCsizeY = station.GetRPCSizeY() * kDistanceToG4;

    const unsigned int nRPC = station.GetNumberRPCChambers();

    {
      ostringstream msg;
      msg << "Number of RPC chambers " << nRPC;
      INFO(msg);
    }

    const CoordinateSystemPtr cs = dStation.GetLocalCoordinateSystem();

    for (unsigned int i = 1; i <= nRPC; ++i) {
      const G4ThreeVector pG4 = Convert(station.GetRPCPosition(i), cs);
      fRPCPositions.push_back(pG4);
      const G4double rG4 = station.GetRPCRotation(i) * CLHEP::degree / utl::degree;
      fRPCRotations.push_back(rG4);
    }

    INFO("RPC positions and rotation in SD station coordinate system:");
    for (unsigned int i = 0; i < fRPCPositions.size(); ++i)
      G4cerr << fRPCPositions[i] / CLHEP::m << " m "
             << fRPCRotations[i] / CLHEP::degree << " deg" << G4endl;

    // Do not build the tank support if any dimension is zero
    fMakeTankSupport =
      fTankSupportTopSlabDimensions[0] &&
      fTankSupportTopSlabDimensions[1] &&
      fTankSupportTopSlabDimensions[2] &&
      fTankSupportOuterFootDimensions[0] &&
      fTankSupportOuterFootDimensions[1] &&
      fTankSupportOuterFootDimensions[2] &&
      fTankSupportOuterFootBaseDimensions[0] &&
      fTankSupportOuterFootBaseDimensions[1] &&
      fTankSupportOuterFootBaseDimensions[2] &&
      fTankSupportCentralFootDimensions[0] &&
      fTankSupportCentralFootDimensions[1] &&
      fTankSupportCentralFootDimensions[2] &&
      fTankSupportCentralFootBaseDimensions[0] &&
      fTankSupportCentralFootBaseDimensions[1] &&
      fTankSupportCentralFootBaseDimensions[2];

    if (!fMakeTankSupport)
      WARNING("Tank support will not be build!");
    else {
      // Check validity of tank support dimensions
      // Feet must have same height
      if (fTankSupportCentralFootDimensions[2] != fTankSupportOuterFootDimensions[2] ||
          fTankSupportCentralFootBaseDimensions[2] != fTankSupportOuterFootBaseDimensions[2]) {
        const string err = "Feet of tank support have different heights! "
                           "Refuse to simulate this configuration!";
        throw utl::InvalidConfigurationException(err);
      }

      // Feet must have same length (i.e. along Y)
      if (fTankSupportCentralFootDimensions[1] != fTankSupportOuterFootDimensions[1] ||
          fTankSupportCentralFootBaseDimensions[1] != fTankSupportOuterFootBaseDimensions[1]) {
        const string err = "Feet of tank support have different lengths! "
                           "Refuse to simulate this configuration!";
        throw utl::InvalidConfigurationException(err);
      }

      // Outer feet can not stick out of top slab along X or Y
      if (fTankSupportOuterFootDistanceToCenter > 0.5*(fTankSupportTopSlabDimensions[0] - fTankSupportOuterFootDimensions[0]) ||
          fTankSupportTopSlabDimensions[1] < fTankSupportCentralFootDimensions[1]) {
        const string err = "Outer feet of tank support stick out of top slab! "
                           "Refuse to simulate this configuration!";
        throw utl::InvalidConfigurationException(err);
      }

      // Feet base can not overlap
      if (fTankSupportOuterFootDistanceToCenter < 0.5*(fTankSupportCentralFootDimensions[0] + fTankSupportOuterFootDimensions[0])) {
        const string err = "Central foot overlaps outer feet! "
                           "Refuse to simulate this configuration !";
        throw utl::InvalidConfigurationException(err);
      }

      if (0.5*fTankSupportCentralFootBaseDimensions[0] > fTankSupportOuterFootDistanceToCenter - 0.5*fTankSupportOuterFootBaseDimensions[0]) {
        const string err = "Base of central foot overlap base of outer feet! "
                           "Refuse to simulate this configuration !";
        throw utl::InvalidConfigurationException(err);
      }

      const G4ThreeVector alBoxDimensions(fAlBoxInnerDimensions[0] + 2*fAlBoxThickness[0],
                                          fAlBoxInnerDimensions[1] + 2*fAlBoxThickness[1],
                                          fAlBoxInnerDimensions[2] + 2*fAlBoxThickness[2]);

      const G4ThreeVector precastHoleDimensions(
        fTankSupportOuterFootDistanceToCenter -
          0.5*(fTankSupportCentralFootDimensions[0] + fTankSupportOuterFootDimensions[0]),
        0.5*fTankSupportTopSlabDimensions[1],
        fTankSupportCentralFootDimensions[2]
      );

      if (alBoxDimensions[0] >= precastHoleDimensions[0] ||
          alBoxDimensions[1] >= precastHoleDimensions[1] ||
          alBoxDimensions[2] >= precastHoleDimensions[2]) {
        ostringstream err;
        err << "RPC housing larger than inner size of concrete tank support!\n"
               "x,y,z (m) : ("
            << alBoxDimensions[0]/CLHEP::m << ", "
            << alBoxDimensions[1]/CLHEP::m << ", "
            << alBoxDimensions[2]/CLHEP::m << ")   >   ("
            << precastHoleDimensions[0]/CLHEP::m << ", "
            << precastHoleDimensions[1]/CLHEP::m << ", "
            << precastHoleDimensions[2]/CLHEP::m << ')';
        throw utl::InvalidConfigurationException(err.str());
      }
    }

    if (fRPCsizeX > fAlBoxInnerDimensions[0] || fRPCsizeY > fAlBoxInnerDimensions[1]) {
      ostringstream err;
      err << "Dimensions of RPC chamber greater than inner size of Al box!\n"
             "x,y (m) : ("
          << fRPCsizeX / CLHEP::m << ", "
          << fRPCsizeY / CLHEP::m << ")   >   ("
          << fAlBoxInnerDimensions[0] / CLHEP::m << ", "
          << fAlBoxInnerDimensions[1] / CLHEP::m << ')';
      throw utl::InvalidConfigurationException(err.str());
    }
  }


  void
  G4StationConstruction::SetXMLParameters()
  {
    CentralConfig& cc = *CentralConfig::GetInstance();

    Branch topBranch = cc.GetTopBranch("G4StationSimulator");

    vector<double> pmt_pp;
    topBranch.GetChild("pmt_momentum_bins").GetData(pmt_pp);

    vector<double> pmt_face_abs;
    topBranch.GetChild("pmtface_AbsLength").GetData(pmt_face_abs);

    vector<double> pmt_interface_abs;
    topBranch.GetChild("interface_AbsLength").GetData(pmt_interface_abs);

    vector<double> pmt_dome_abs;
    topBranch.GetChild("dome_AbsLength").GetData(pmt_dome_abs);

    topBranch.GetChild("interface_thickness").GetData(fInterfaceThickness);
    fInterfaceThickness *= kDistanceToG4;

    topBranch.GetChild("pmtface_thickness").GetData(fGlassThickness);
    fGlassThickness *= kDistanceToG4;

    {
      const unsigned int size = pmt_pp.size();

      if (pmt_face_abs.size() != size ||
          pmt_interface_abs.size() != size ||
          pmt_dome_abs.size() != size) {
        const string err = "TabulatedFunction size mismatch for pmt properties";
        ERROR(err);
        throw utl::OutOfBoundException(err);
      }

      for (unsigned int i = 0; i < size; ++i) {
        fPmtfaceABSORPTION.PushBack(pmt_pp[i] * kEnergyToG4, pmt_face_abs[i] * kDistanceToG4);
        fgInterfaceABSORPTION.PushBack(pmt_pp[i] * kEnergyToG4, pmt_interface_abs[i] * kDistanceToG4);
        fgPmtdomeABSORPTION.PushBack( pmt_pp[i] * kEnergyToG4, pmt_dome_abs[i] * kDistanceToG4);
      }
    }

    topBranch.GetChild("face_radius").GetData(fFaceRadius);
    fFaceRadius *= kDistanceToG4;

    topBranch.GetChild("face_radius_z").GetData(fFaceRadiusz);
    fFaceRadiusz *= kDistanceToG4;

    topBranch.GetChild("face_active_radius").GetData(fFaceActiveRadius);
    fFaceActiveRadius *= kDistanceToG4;

    // hardware-version-specific detector parameters
    if (fHasSmallPMT) {
      topBranch.GetChild("interface_thickness_spmt").GetData(fInterfaceThickness_SPMT);
      fInterfaceThickness_SPMT *= kDistanceToG4;

      topBranch.GetChild("pmtface_thickness_spmt").GetData(fGlassThickness_SPMT);
      fGlassThickness_SPMT *= kDistanceToG4;

      topBranch.GetChild("dome_thickness_spmt").GetData(fDomeThickness_SPMT);
      fDomeThickness_SPMT *= kDistanceToG4;

      topBranch.GetChild("face_radius_spmt").GetData(fFaceRadius_SPMT);
      fFaceRadius_SPMT *= kDistanceToG4;

      topBranch.GetChild("face_active_radius_spmt").GetData(fFaceActiveRadius_SPMT);
      fFaceActiveRadius_SPMT *= kDistanceToG4;

      topBranch.GetChild("window_radius_spmt").GetData(fWindowRadius_SPMT);
      fWindowRadius_SPMT *= kDistanceToG4;

    }

    vector<double> liner_pp;
    topBranch.GetChild("liner_momentum_bins").GetData(liner_pp);

    vector<double> liner_rindex;
    topBranch.GetChild("liner_Rindex").GetData(liner_rindex);

    vector<double> liner_abs;
    topBranch.GetChild("liner_AbsLength").GetData(liner_abs);

    vector<double> liner_backscatter;
    topBranch.GetChild("backScatter").GetData(liner_backscatter);

    topBranch.GetChild("dome_thickness").GetData(fDomeThickness);
    fDomeThickness *= kDistanceToG4;

    {
      const unsigned int size = liner_pp.size();

      if (liner_rindex.size() != size ||
          liner_abs.size() != size ||
          liner_backscatter.size() != size) {
        const string err = "Tabulated function size mismatch for liner properties";
        ERROR(err);
        throw utl::OutOfBoundException(err);
      }

      for (unsigned int i = 0; i < size; ++i) {
        fLinerTYVEK_RINDEX.PushBack(liner_pp[i] * kEnergyToG4, liner_rindex[i]);
        fLinerABSORPTION.PushBack(liner_pp[i] * kEnergyToG4, liner_abs[i] * kDistanceToG4);
        fgLinerBACKSCATTERCONSTANT.PushBack(liner_pp[i] * kEnergyToG4, liner_backscatter[i]);
      }
    }
  }


  G4VPhysicalVolume*
  G4StationConstruction::Construct()
  {
    if (!expHall_phys) {
      CreateElements();
      CreateMaterials();
      return CreateTank();
    }
    return expHall_phys;
  }


  void
  G4StationConstruction::CreateElements()
  {
    // G4Element(name, symbol, Z, n(molar mass))
    Renew(elN, new G4Element("Nitrogen", "N", 7, 14.01 * g/mole));
    Renew(elO, new G4Element("Oxygen", "O", 8, 16.00 * g/mole));
    Renew(elAl, new G4Element("Aluminum", "Al", 13, 29.98 * g/mole));
    Renew(elFe, new G4Element("Iron", "Fe", 26, 55.85 * g/mole));
    Renew(elH, new G4Element("Hydrogen", "H", 1, 1.01 * g/mole));
    Renew(elSi, new G4Element("Silicon", "Si", 14, 28.09 * g/mole));
    Renew(elB, new G4Element("Boron", "B", 5, 10.811 * g/mole));
    Renew(elNa, new G4Element("Sodium", "Na", 11, 22.98977 * g/mole));
    Renew(elC, new G4Element("Carbon", "C", 6, 12.0107 * g/mole));
    Renew(elCl, new G4Element("Chlorine", "Cl", 17, 35.453 * g/mole));
    Renew(elTi, new G4Element("Titanium", "Ti", 22, 47.867 * g/mole));

    // Re-used material for Pyrex
    const G4double density = 2.23*g/cm3;

    // G4Material(name, density, nel)
    Renew(SiO2, new G4Material("SiO2", density, 2));
    // AddElement(ptr, natoms)
    SiO2->AddElement(elSi, 1);
    SiO2->AddElement(elO, 2);

    Renew(B2O2, new G4Material("B2O2", density, 2));
    B2O2->AddElement(elB, 2);
    B2O2->AddElement(elO, 2);

    Renew(Na2O, new G4Material("Na2O", density, 2));
    Na2O->AddElement(elNa, 2);
    Na2O->AddElement(elO, 1);

    Renew(TiO2, new G4Material( "TiO2", 1.060 * g/cm3, 2));
    TiO2->AddElement(elTi, 1);
    TiO2->AddElement(elO, 2);
  }


  void
  G4StationConstruction::CreateMaterials()
  {
    CreateAir();
    CreateWater();
    CreateVacuum();
    CreatePyrex();
    CreatePyrex1();
    CreateLucite();
    CreateInterface();
    CreateHDPE();
    CreateLiner();
    CreateAluminium();
    CreateDirt();
    CreatePolystyrene();
    CreateExpandedPolystyreneFoam();
    CreateExtrudedPolystyreneFoam();
    CreatePVC();
    CreateScintillator();
    CreateWLS();

    CreateConcrete();
    CreateAcrylic();
    CreateSodaLimeGlass();
    CreateR134A();
    CreateBakelite();
  }


  G4VPhysicalVolume*
  G4StationConstruction::CreateTank()
  {
    CreatePrimitives();
    CreateHall();
    AssembleStation();
    return expHall_phys;
  }


  void
  G4StationConstruction::CreateAir()
  {
    Renew(Air, new G4Material("Air", 1.29e-3 * g/cm3, 2));
    Air->AddElement(elN, 0.7);
    Air->AddElement(elO, 0.3);

    // Add material property table for Air to simulate empty tank
    double airPP[2] = { 2.08*eV, 4.20*eV };
    double airRINDEX[2] = { 1.000273, 1.000273 };
    double airABSLENGTH[2] = { 10000*m, 10000*m };

    Renew(airMPT, new G4MaterialPropertiesTable());
    airMPT->AddProperty("RINDEX", airPP, airRINDEX, 2);
    airMPT->AddProperty("ABSLENGTH", airPP, airABSLENGTH, 2);

    Air->SetMaterialPropertiesTable(airMPT);
  }


  void
  G4StationConstruction::CreateWater()
  {
    Renew(Water, new G4Material("Water", 1 * g/cm3, 2));
    Water->AddElement(elH, 2);
    Water->AddElement(elO, 1);

    // Set water material properties
    Renew(waterMPT, new G4MaterialPropertiesTable());

    // water rindex
    AddProperty(waterMPT, "RINDEX", fWaterRINDEX);

    // water abs. length
    AddProperty(waterMPT, "ABSLENGTH", fgWaterABSORPTION);

    Water->SetMaterialPropertiesTable(waterMPT);
  }


  void
  G4StationConstruction::CreateVacuum()
  {
    // ATTENTION: change from STP to avoid the warning...
    Renew(Vacuum, new G4Material("Vacuum", 1, 1 * g/mole, universe_mean_density,
                                 kStateUndefined, STP_Temperature, STP_Pressure));
  }


  void
  G4StationConstruction::CreatePyrex()
  {
    // Borosilicate glass (Pyrex) for the active pmt faces
    //----------------------------------------------------
    /*
      NOTE : The "face" is actually an ellipsoidal surface
      The active area (photocathode) covers 90% of it.
      The refractive index is the same as the dome
      refractive index. The absorption length is set
      so that no photons propagate through the
      photocathode.
    */

    // from PDG:
    // 80% SiO2 + 13% B2O2 + 7% Na2O
    // by fractional mass?

    Renew(Pyrex, new G4Material("Pyrex", 2.23 * g/cm3, 3));
    Pyrex->AddMaterial(SiO2, 0.80);
    Pyrex->AddMaterial(B2O2, 0.13);
    Pyrex->AddMaterial(Na2O, 0.07);

    // Fill material properties table for pmt face
    Renew(pmtfaceMPT, new G4MaterialPropertiesTable());

    // pmt face rindex
    AddProperty(pmtfaceMPT, "RINDEX", fPmtfaceRINDEX);

    // Fill and add pmt face absorption length
    AddProperty(pmtfaceMPT, "ABSLENGTH", fPmtfaceABSORPTION);

    Pyrex->SetMaterialPropertiesTable(pmtfaceMPT);
  }


  void
  G4StationConstruction::CreatePyrex1()
  {
    // Borosilicate glass (Pyrex1) for the inactive part of the pmt faces
    /* NOTE : The inactive part of the "face" is about 10% of
      the total ellipsoidal surface seen by the photons.
      The refractive index is the same as the dome
      refractive index. The absorption length is set
      so that photons can propagate through this part.
      To simplify things we use the same absorption length
      as for the dome. */
    const G4double density = 2.23 * g/cm3;

    Renew(Pyrex1, new G4Material("Pyrex1", density, 3));
    Pyrex1->AddMaterial(SiO2, 0.80);
    Pyrex1->AddMaterial(B2O2, 0.13);
    Pyrex1->AddMaterial(Na2O, 0.07);

    Renew(pmtfaceMPT1, new G4MaterialPropertiesTable());

    // Fill material properties table for non sensitive pmt face

    // non-sensitive pmt face rindex = sensitive pmt face rindex
    AddProperty(pmtfaceMPT1, "RINDEX", fPmtfaceRINDEX);

    // Fill and add non-sensitive pmt face absorption length = dome absorption length
    AddProperty(pmtfaceMPT1, "ABSLENGTH", fgPmtdomeABSORPTION);

    Pyrex1->SetMaterialPropertiesTable(pmtfaceMPT1);
  }


  void
  G4StationConstruction::CreateLucite()
  {
    // Lucite (for the pmt domes)
    const G4double density = 1.20 * g/cm3;

    Renew(CH3, new G4Material("CH3", density, 2));
    CH3->AddElement(elC, 1);
    CH3->AddElement(elH, 3);

    Renew(CH2, new G4Material("CH2", density, 2));
    CH2->AddElement(elC, 1);
    CH2->AddElement(elH, 2);

    Renew(C, new G4Material("C", density, 1));
    C->AddElement(elC, 1);

    Renew(CO2, new G4Material("CO2", density, 2));
    CO2->AddElement(elC, 1);
    CO2->AddElement(elO, 2);

    /* NOTE : Do not know the fractional mass comp.
      of lucite, only the molecular composition.
      Perhaps it is a standard material defined in g4.
      for now define by number of atoms.
      It really doesn't matter at this point
      anyway, since we are more concerned about
      the optical properties. */

    Renew(Lucite, new G4Material("Lucite", density, 3));
    Lucite->AddElement(elC, 4);
    Lucite->AddElement(elH, 8);
    Lucite->AddElement(elO, 2);

    Renew(pmtdomeMPT, new G4MaterialPropertiesTable());

    // pmt dome rindex
    AddProperty(pmtdomeMPT, "RINDEX", fgPmtdomeRINDEX);

    // pmt dome absorption length
    AddProperty(pmtdomeMPT, "ABSLENGTH", fgPmtdomeABSORPTION);

    Lucite->SetMaterialPropertiesTable(pmtdomeMPT);
  }


  void
  G4StationConstruction::CreateInterface()
  {
    // Interface for production tanks is Wacker SilGel 612 (Registered trademark of Wacker Group
    const G4double density = 0.97 * g/cm3; // guess - not explicitly set by Troy

    Renew(Interface, new G4Material("Interface", density, 3));
    Interface->AddElement(elC, 4);
    Interface->AddElement(elH, 8);
    Interface->AddElement(elO, 2);

    Renew(interfaceMPT, new G4MaterialPropertiesTable());

    // interface rindex
    AddProperty(interfaceMPT, "RINDEX", fgInterfaceRINDEX);

    // interface abs. length
    AddProperty(interfaceMPT, "ABSLENGTH", fgInterfaceABSORPTION);

    Interface->SetMaterialPropertiesTable(interfaceMPT);
  }


  void
  G4StationConstruction::CreateHDPE()
  {
    // HDPE (for tank walls)
    const G4double density = 0.94 * g/cm3;

    Renew(HDPE, new G4Material("HDPE", density, 2));
    HDPE->AddElement(elC, 2);
    HDPE->AddElement(elH, 4);

    Renew(linerMPT, new G4MaterialPropertiesTable());

    // liner abs. length
    AddProperty(linerMPT, "ABSLENGTH", fLinerABSORPTION);

    HDPE->SetMaterialPropertiesTable(linerMPT);
  }


  void
  G4StationConstruction::CreateLiner()
  {
    // Optical surface properties

    // Define the optical model and set optical properties of the surfaces

    /* Note that diffuse reflection is implicitly defined, since:
      SPECULARSPIKECONSTANT+SPECULARLOBECONSTANT+BACKSCATTERCONSTANT=
      1-DIFFUSEL0BECONSTANT */

    Renew(linerOpticalMPT, new G4MaterialPropertiesTable());
    AddProperty(linerOpticalMPT, "SPECULARLOBECONSTANT", fgLinerSPECULARLOBECONSTANT);
    AddProperty(linerOpticalMPT, "SPECULARSPIKECONSTANT", fgLinerSPECULARSPIKECONSTANT);
    AddProperty(linerOpticalMPT, "REFLECTIVITY", fgLinerREFLECTIVITY);
    AddProperty(linerOpticalMPT, "BACKSCATTERCONSTANT", fgLinerBACKSCATTERCONSTANT);
    AddProperty(linerOpticalMPT, "RINDEX", fLinerTYVEK_RINDEX);

    Renew(OpLinerSurface, new G4OpticalSurface("WallSurface"));
    OpLinerSurface->SetModel(unified);
    OpLinerSurface->SetType(dielectric_metal);
    OpLinerSurface->SetFinish(ground);
    OpLinerSurface->SetSigmaAlpha(fgSIGMA_ALPHA);
    OpLinerSurface->SetMaterialPropertiesTable(linerOpticalMPT);
  }


  void
  G4StationConstruction::CreateAluminium()
  {
    // Aluminum for electronics boxes and (possibly) solar panel
    const G4double density = 2.700 * g/cm3;
    const G4double a = 29.98 * g/mole;
    Renew(Al, new G4Material("Aluminum", /*z*/13, a, density));
  }


  void
  G4StationConstruction::CreateDirt()
  {
    // Dirt for the ground
    /* Soil density samples (see Informe_Preliminar__Intemin_1_11_07
      Samples at:
         A.- Hilda           tank
         B.- Lety            tank
         C.- TierraDelFuego  tank */

    const G4double density = (
      2.43 + 2.39 + 2.41 +  // Hilda
      2.39 + 2.45 + 2.32 +  // Lety
      2.35 + 2.40 + 2.28    // TierraDelFuego
    ) / 9 * g/cm3;
    Renew(Dirt, new G4Material("Dirt", density, 4));
    Dirt->AddElement(elO, 51*perCent);
    Dirt->AddElement(elSi, 35.2*perCent);
    Dirt->AddElement(elAl, 8.2*perCent);
    Dirt->AddElement(elFe, 5.6*perCent);
  }


  void
  G4StationConstruction::CreatePolystyrene()
  {
    // Polystyrene for the scintillator bars Dow STYRON 663 W
    const G4double density = 1.04 * g/cm3;
    Renew(Polystyrene, new G4Material("Polystyrene", density, 2));
    Polystyrene->AddElement(elC, 8);
    Polystyrene->AddElement(elH, 8);
  }


  void
  G4StationConstruction::CreateScintillator()
  {
    /* Scintillator made of Polystyrene ( C6H5-CH-CH2 ) + PPO + POPOP
      Plastic:
      Is a replication on an aromatic ring (C6H5) with one out of the six H
      replaced by a CH bond with another CH2.

                 CH2
                //
                CH
                |
               / \
       C6H5:  |   |
               \ /
    */

    {
      const G4double density = 1.060 * g/cm3;
      // PPO:
      Renew(PPO, new G4Material("PPO", density, 4));
      PPO->AddElement(elC, 15);
      PPO->AddElement(elH, 11);
      PPO->AddElement(elN, 1);
      PPO->AddElement(elO, 1);
      // POPOP:
      Renew(POPOP, new G4Material("POPOP", density, 4));
      POPOP->AddElement(elC, 24);
      POPOP->AddElement(elH, 16);
      POPOP->AddElement(elN, 2);
      POPOP->AddElement(elO, 2);
      // Scintillator material:
      Renew(ScintMat, new G4Material("ScintMat", density, 3));
      ScintMat->AddMaterial(Polystyrene, 98.7*perCent);
      ScintMat->AddMaterial(PPO, 1*perCent);
      ScintMat->AddMaterial(POPOP, 0.3*perCent);
    }

    // Scintillator Coating - 15% TiO2 and 85% polystyrene by weight.
    const G4double density = 1.52 * g/cm3;
    Renew(ScintCoating, new G4Material("ScintCoating", density, 2));
    ScintCoating->AddMaterial(TiO2, 15*perCent);
    ScintCoating->AddMaterial(Polystyrene, 85*perCent);

    // Do not add scintillation properties to material for fast mode (PE photons will be created
    // from parametriezed curve in (see G4UMDScintStripAction.cc)
    if (!fG4StationSimulator.fUMDFastMode)
      SetScintPropertyTable();
  }


  void
  G4StationConstruction::CreateWLS()
  {
    // WLSfiber PMMA
    Renew(PMMA, new G4Material("PMMA", 1.190 * g/cm3, 3));
    PMMA->AddElement(elC, 5);
    PMMA->AddElement(elH, 8);
    PMMA->AddElement(elO, 2);

    // Cladding (polyethylene)
    Renew(Pethylene, new G4Material("Pethylene", 1.200 * g/cm3, 2));
    Pethylene->AddElement(elC, 2);
    Pethylene->AddElement(elH, 4);

    // Double Cladding (fluorinated polyethylene)
    Renew(FPethylene, new G4Material("FPethylene", 1.400*g/cm3, 2));
    FPethylene->AddElement(elC, 2);
    FPethylene->AddElement(elH, 4);

    G4double photonEnergy[] = {
      2.00, 2.03, 2.06, 2.09, 2.12,
      2.15, 2.18, 2.21, 2.24, 2.27,
      2.30, 2.33, 2.36, 2.39, 2.42,
      2.45, 2.48, 2.51, 2.54, 2.57,
      2.60, 2.63, 2.66, 2.69, 2.72,
      2.75, 2.78, 2.81, 2.84, 2.87,
      2.90, 2.93, 2.96, 2.99, 3.02,
      3.05, 3.08, 3.11, 3.14, 3.17,
      3.20, 3.23, 3.26, 3.29, 3.32,
      3.35, 3.38, 3.41, 3.44, 3.47
    };
    Scale(photonEnergy, CLHEP::eV);

    const auto n = utl::Length(photonEnergy);

    // PMMA for WLSfibers
    G4double refractiveIndexWLSfiber[] = {
      1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60,
      1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60,
      1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60,
      1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60,
      1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60
    };
    assert(utl::Length(refractiveIndexWLSfiber) == n);

    G4double absWLSfiber[] = {
     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00,
     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00,
     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 1.0,
     0.0001, 0.0001, 0.0001, 0.0001, 0.0001,0.0001, 0.0001, 0.0001, 0.0001, 0.0001,
     0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001
    };
    Scale(absWLSfiber, CLHEP::m);
    assert(utl::Length(absWLSfiber) == n);

    G4double emissionFib[] = {
      0.05, 0.10, 0.30, 0.50, 0.75, 1.00, 1.50, 1.85, 2.30, 2.75,
      3.25, 3.80, 4.50, 5.20, 6.00, 7.00, 8.50, 9.50, 11.1, 12.4,
      12.9, 13.0, 12.8, 12.3, 11.1, 11.0, 12.0, 11.0, 17.0, 16.9,
      15.0, 9.00, 2.50, 1.00, 0.05, 0.00, 0.00, 0.00, 0.00, 0.00,
      0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00
    };
    assert(utl::Length(emissionFib) == n);

    // Add entries into properties table
    G4MaterialPropertiesTable* const mptWLSfiber = new G4MaterialPropertiesTable();
    mptWLSfiber->AddProperty("RINDEX", photonEnergy, refractiveIndexWLSfiber, n);
    mptWLSfiber->AddProperty("WLSABSLENGTH", photonEnergy, absWLSfiber, n);
    mptWLSfiber->AddProperty("WLSCOMPONENT", photonEnergy, emissionFib, n);
    mptWLSfiber->AddConstProperty("WLSTIMECONSTANT", 20*ns);

    PMMA->SetMaterialPropertiesTable(mptWLSfiber);

    //  Polyethylene

    G4double refractiveIndexClad1[] = {
      1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49,
      1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49,
      1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49,
      1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49,
      1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49, 1.49
    };

    assert(utl::Length(refractiveIndexClad1) == n);

    G4double absClad[] = {
      20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0,
      20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0,
      20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0,
      20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0,
      20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0
    };
    Scale(absClad, CLHEP::m);
    assert(utl::Length(absClad) == n);

    // Add entries into properties table
    G4MaterialPropertiesTable* const mptClad1 = new G4MaterialPropertiesTable();
    mptClad1->AddProperty("RINDEX", photonEnergy, refractiveIndexClad1, n);
    mptClad1->AddProperty("ABSLENGTH", photonEnergy, absClad, n);

    Pethylene->SetMaterialPropertiesTable(mptClad1);

    // Fluorinated Polyethylene

    G4double refractiveIndexClad2[] = {
      1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42,
      1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42,
      1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42,
      1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42,
      1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42, 1.42
    };

    assert(utl::Length(refractiveIndexClad2) == n);

    // Add entries into properties table
    G4MaterialPropertiesTable* const mptClad2 = new G4MaterialPropertiesTable();
    mptClad2->AddProperty("RINDEX", photonEnergy, refractiveIndexClad2, n);
    mptClad2->AddProperty("ABSLENGTH", photonEnergy, absClad, n);

    FPethylene->SetMaterialPropertiesTable(mptClad2);
  }


  void
  G4StationConstruction::SetScintPropertyTable()
  {
    // ------------ Generate & Add Material Properties Table ------------

    G4double photonEnergy[] = {
      2.00, 2.03, 2.06, 2.09, 2.12,
      2.15, 2.18, 2.21, 2.24, 2.27,
      2.30, 2.33, 2.36, 2.39, 2.42,
      2.45, 2.48, 2.51, 2.54, 2.57,
      2.60, 2.63, 2.66, 2.69, 2.72,
      2.75, 2.78, 2.81, 2.84, 2.87,
      2.90, 2.93, 2.96, 2.99, 3.02,
      3.05, 3.08, 3.11, 3.14, 3.17,
      3.20, 3.23, 3.26, 3.29, 3.32,
      3.35, 3.38, 3.41, 3.44, 3.47
    };
    Scale(photonEnergy, CLHEP::eV);
    const auto n = utl::Length(photonEnergy);

    //  Scintillator (Polystyrene based)

    G4double refractiveIndexPS[] = {
      1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50,
      1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50,
      1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50,
      1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50,
      1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50
    };
    assert(utl::Length(refractiveIndexPS) == n);

    G4double absPS[] = {
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0,
      20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0
    };
    Scale(absPS, CLHEP::cm);
    assert(utl::Length(absPS) == n);

    G4double scintilFast[] = {
      0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
      0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
      0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
      1, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0
    };
    assert(utl::Length(scintilFast) == n);

    G4double scintilSlow[] = {
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
    };
    assert(utl::Length(scintilSlow) == n);

    // Add entries into properties table
    G4MaterialPropertiesTable* const mptPolystyrene = new G4MaterialPropertiesTable();
    mptPolystyrene->AddProperty("RINDEX", photonEnergy, refractiveIndexPS, n);
    mptPolystyrene->AddProperty("ABSLENGTH", photonEnergy, absPS, n);
    mptPolystyrene->AddProperty("FASTCOMPONENT", photonEnergy, scintilFast, n);
    mptPolystyrene->AddProperty("SLOWCOMPONENT", photonEnergy, scintilSlow, n);

    mptPolystyrene->AddConstProperty("SCINTILLATIONYIELD", fScintYield / CLHEP::keV);
    mptPolystyrene->AddConstProperty("RESOLUTIONSCALE", 1);
    mptPolystyrene->AddConstProperty("FASTTIMECONSTANT", 5.*CLHEP::ns);
    mptPolystyrene->AddConstProperty("SLOWTIMECONSTANT", 50.*CLHEP::ns);
    mptPolystyrene->AddConstProperty("YIELDRATIO", 1);

    ScintMat->SetMaterialPropertiesTable(mptPolystyrene);

    // Set the Birks Constant for the Polystyrene scintillator
    //ScintMat->GetIonisation()->SetBirksConstant(0.126*CLHEP::mm/CLHEP::MeV);

    // Optical Skin properties
    // Strip coating properties
    const double extrusion_polish = 1;
    Renew(fUMDScintSkinSurf, new G4OpticalSurface("StripSkin", glisur, ground, dielectric_metal, extrusion_polish));
    Renew(fUMDScintSkinSurfBack, new G4OpticalSurface("StripSkinBack", glisur, ground, dielectric_metal, extrusion_polish));

    const double refl = 1;
    G4double skin_refl[] = {
      refl, refl, refl, refl, refl, refl, refl, refl, refl, refl,
      refl, refl, refl, refl, refl, refl, refl, refl, refl, refl,
      refl, refl, refl, refl, refl, refl, refl, refl, refl, refl,
      refl, refl, refl, refl, refl, refl, refl, refl, refl, refl,
      refl, refl, refl, refl, refl, refl, refl, refl, refl, refl
    };
    assert(utl::Length(skin_refl) == n);

    G4double skin_effi[] = {
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0
    };
    assert(utl::Length(skin_effi) == n);

    G4MaterialPropertiesTable* const scint_skin_MPT = new G4MaterialPropertiesTable();
    G4MaterialPropertiesTable* const scint_skin_back_MPT = new G4MaterialPropertiesTable();
    scint_skin_MPT->AddProperty("REFLECTIVITY", photonEnergy, skin_refl, n);
    scint_skin_MPT->AddProperty("EFFICIENCY", photonEnergy, skin_effi, n);
    // black painted (total absorption)
    scint_skin_back_MPT->AddProperty("REFLECTIVITY", photonEnergy, skin_effi, n);
    scint_skin_back_MPT->AddProperty("EFFICIENCY", photonEnergy, skin_effi, n);

    fUMDScintSkinSurf->SetMaterialPropertiesTable(scint_skin_MPT);
    fUMDScintSkinSurfBack->SetMaterialPropertiesTable(scint_skin_back_MPT);
  }


  void
  G4StationConstruction::SetScintCoatingProperties()
  {
    new G4LogicalSkinSurface("SkinTopBot_large", umd_top_bot_coat_log_large, fUMDScintSkinSurf);
    new G4LogicalSkinSurface("SkinSide_large", umd_side_coat_log_large, fUMDScintSkinSurf);
    new G4LogicalSkinSurface("SkinTopBot_small", umd_top_bot_coat_log_small, fUMDScintSkinSurf);
    new G4LogicalSkinSurface("SkinSide_small", umd_side_coat_log_small, fUMDScintSkinSurf);
    new G4LogicalSkinSurface("SkinBack", umd_back_side_coat_log, fUMDScintSkinSurfBack);
  }


  void
  G4StationConstruction::CreateExpandedPolystyreneFoam()
  {
    // Styrofoam used to fill upper portion of SSD
    Renew(ExpandedPolystyreneFoam, new G4Material("ExpandedPolystyreneFoam", 0.01414 * g/cm3, 2));
    ExpandedPolystyreneFoam->AddElement(elC, 8);
    ExpandedPolystyreneFoam->AddElement(elH, 8);
  }


  void
  G4StationConstruction::CreateExtrudedPolystyreneFoam()
  {
    // Foam in-between aluminum sheets of the sandwich panel
    Renew(ExtrudedPolystyreneFoam, new G4Material("ExtrudedPolystyreneFoam", 0.032 * g/cm3, 2));
    ExtrudedPolystyreneFoam->AddElement(elC, 8);
    ExtrudedPolystyreneFoam->AddElement(elH, 8);
  }


  void
  G4StationConstruction::CreatePVC()
  {
    // PVC for the casing
    // flexible PVC
    Renew(PVC, new G4Material("PVC", 1.35 * g/cm3, 3));
    PVC->AddElement(elC, 2);
    PVC->AddElement(elH, 3);
    PVC->AddElement(elCl, 1);
  }


  void
  G4StationConstruction::CreateConcrete()
  {
    fG4_Concrete = fNistManager->FindOrBuildMaterial("G4_CONCRETE");
  }


  void
  G4StationConstruction::CreateAcrylic()
  {
    fG4_Acrylic = fNistManager->FindOrBuildMaterial("G4_POLYACRYLONITRILE");
  }


  void
  G4StationConstruction::CreateSodaLimeGlass()
  {
    G4Material* const G4_CaO = fNistManager->FindOrBuildMaterial("G4_CALCIUM_OXIDE");
    G4Material* const G4_MgO = fNistManager->FindOrBuildMaterial("G4_MAGNESIUM_OXIDE");

    Renew(fSoda_lime_glass, new G4Material("Soda-lime_glass", 2.52 * g/cm3, 4));
    fSoda_lime_glass->AddMaterial(SiO2, 0.73);
    fSoda_lime_glass->AddMaterial(Na2O, 0.14);
    fSoda_lime_glass->AddMaterial(G4_CaO, 0.09);
    fSoda_lime_glass->AddMaterial(G4_MgO, 0.04);
  }


  void
  G4StationConstruction::CreateR134A()
  {
    const vector<G4String> elements({ "H", "C", "F" });
    const vector<G4int> nbAtoms({ 2, 2, 4 });
    fNistManager->ConstructNewMaterial(
      "TetraFluoroethane_STP",
      elements, nbAtoms, 0.00425 * g/cm3,
      true, kStateGas
    );
    fR134a = fNistManager->ConstructNewGasMaterial(
      "TetraFluoroethane", "TetraFluoroethane_STP",
      273.*CLHEP::kelvin,
      0.8*CLHEP::atmosphere
    );
  }


  void
  G4StationConstruction::CreateBakelite()
  {
    fG4_Bakelite = fNistManager->FindOrBuildMaterial("G4_BAKELITE");
  }


  void
  G4StationConstruction::CreatePrimitives()
  {
    Renew(expHall_box, new G4Box("World", fExpHall_x, fExpHall_y, fExpHall_z));
    Renew(ground_solid, new G4Tubs("ground_solid", 0, fSimulationRadius, fGroundThickness/2, 0, 360*deg));
    // tank casing
    Renew(tank_solid, new G4Tubs("tank_solid", 0, fgTankRadius, fgTankHalfHeight, 0, 360*deg));
    Renew(top_solid, new G4Tubs("top_solid", 0, fgTankRadius + fgTankThickness, fgTankThickness/2, 0, 360*deg));
    Renew(side_solid, new G4Tubs("side_solid", fgTankRadius, fgTankRadius + fgTankThickness, fgTankHalfHeight, 0, 360*deg));
    Renew(inner_solid, new G4Ellipsoid("inner_solid", fPmtRmin, fPmtRmin, fPmtRzmin, -fPmtRzmin, 0));
    // Active photocathode covered area
    Renew(pmt_aux, new G4Ellipsoid("pmt_aux", fgPmtRmax, fgPmtRmax, fgPmtRzmax, -fgPmtRzmax, -fgHeightz));
    Renew(pmt_solid, new G4SubtractionSolid("pmt_solid", pmt_aux, inner_solid));
    Renew(pmt_aux1, new G4Ellipsoid("pmt_aux1", fgPmtRmax, fgPmtRmax, fgPmtRzmax, -fgHeightz, 0));
    Renew(pmt_solid1, new G4SubtractionSolid("pmt_solid1", pmt_aux1, inner_solid));

    Renew(interface_in_aux, new G4Ellipsoid("interface_in_aux", fInterfaceRmin, fInterfaceRmin, fInterfaceRzmin, -fInterfaceRzmin, 0));
    Renew(interface_out_aux, new G4Ellipsoid("interface_out_aux", fgInterfaceRmax, fgInterfaceRmax, fgInterfaceRzmax, -fgInterfaceRzmax, 0));
    Renew(interface_solid, new G4SubtractionSolid("interface_solid", interface_out_aux, interface_in_aux));

    Renew(dome_in_aux, new G4Ellipsoid("dome_in_aux", fDomeRmin, fDomeRmin, fDomeRzmin, -fDomeRzmin, 0));
    Renew(dome_out_aux, new G4Ellipsoid("dome_out_aux", fgDomeRmax, fgDomeRmax, fgDomeRzmax, -fgDomeRzmax, 0));
    Renew(dome_solid, new G4SubtractionSolid("dome_solid", dome_out_aux, dome_in_aux));

    if (fHasSmallPMT) {
      // SPMT AT LINER LEVEL
      Renew(s_pmt_solid, new G4Tubs("s_pmt_solid", 0, fPmtRmin_SPMT, fGlassThickness_SPMT/2, 0,360*deg));
      Renew(s_pmt_aux, new G4Tubs("s_pmt_aux", 0, fPmtRmin_SPMT, fGlassThickness_SPMT/2, 0, 360*deg));
      Renew(s_pmt_aux1, new G4Tubs("s_pmt_aux1", 0, fPmtRmax_SPMT, fGlassThickness_SPMT/2, 0, 360*deg));
      Renew(s_pmt_solid1, new G4SubtractionSolid("s_pmt_solid1", s_pmt_aux1, s_pmt_aux));
      Renew(s_interface_solid, new G4Tubs("s_interface_solid", 0, fInterfaceRadius_SPMT, fInterfaceThickness_SPMT/2, 0, 360*deg));
      Renew(s_dome_solid, new G4Tubs("s_dome_solid", 0, fDomeRadius_SPMT, fDomeThickness_SPMT/2, 0, 360*deg));
    }

    // Creates Amiga underground primitives
    if (fUMDEnabled) {
      const auto& dStation = G4StationSimulator::fgCurrent.GetDetectorStation();
      if (!dStation.ExistsAssociatedCounter()) {
        ostringstream err;
        err << "Trying to simulate the UMD in a WCD (Id "<< dStation.GetId() << ") "
               "that has not an underground AMIGA scintillator associated!";
        WARNING(err);
      }
      CreatePrimitivesUMD();
    }

    // Creates MARTA primitives
    if (fMARTAEnabled)
      CreatePrimitivesMARTA();
  }


  // Creates casings for scintillators
  void
  G4StationConstruction::CreatePrimitivesUMD()
  {
    fCoatingThickness = 0.25*CLHEP::mm;
    fCladingThickness = 0.10*CLHEP::mm;
    fPixelL = 0.5*CLHEP::cm;

    INFO("UMD areas:");
    for (const auto& ml : fModAreaLenghts) {
      const auto module = ml.first;
      const auto len  = ml.second;
      const double x = len + fUMDManifoldL/2;
      const double y = fUMDScintsW/2 * fUMDScints/2;
      const double z = fUMDScintsH/2;
      const double t = fUMDThicknessCasing;
      if (module == mdet::Module::AreaKind::eLarge) {

        Renew(umd_casing_solid_large, new G4Box("umd_casing_large", x+t, y+t, z+t));
        Renew(umd_strip_solid_large, new G4Box("umd_strip_solid_large", len/2 - fCoatingThickness, fUMDScintsW/2 - fCoatingThickness, fUMDScintsH/2 - fCoatingThickness));
        Renew(umd_top_bot_coat_solid_large, new G4Box("umd_top_bot_coat_solid_large", len/2 - fCoatingThickness, fUMDScintsW/2 - fCoatingThickness, fCoatingThickness/2));
        Renew(umd_side_coat_solid_large, new G4Box("umd_side_coat_solid_large", len/2 - fCoatingThickness, fCoatingThickness/2, fUMDScintsH/2));

        Renew(umd_fiber_core_solid_large, new G4Tubs("umd_fiber_core_solid_large", 0, fUMDFiberRadius - 2*fCladingThickness, len/2 - fCoatingThickness, 0, 360*CLHEP::deg));
        Renew(umd_fiber_clad1_solid_large, new G4Tubs("umd_fiber_clad1_solid_large", 0, fUMDFiberRadius - fCladingThickness, len/2 - fCoatingThickness, 0, 360*CLHEP::deg));
        Renew(umd_fiber_clad2_solid_large, new G4Tubs("umd_fiber_clad2_solid_large", 0, fUMDFiberRadius, len/2 - fCoatingThickness, 0, 360*CLHEP::deg));

      } else {

        Renew(umd_casing_solid_small, new G4Box("umd_casing_small", x+t, y+t, z+t));
        Renew(umd_strip_solid_small, new G4Box("umd_strip_solid_small", len/2 - fCoatingThickness, fUMDScintsW/2 - fCoatingThickness, fUMDScintsH/2 - fCoatingThickness));
        Renew(umd_top_bot_coat_solid_small, new G4Box("umd_top_bot_coat_solid_small", len/2 - fCoatingThickness, fUMDScintsW/2 - fCoatingThickness, fCoatingThickness/2));
        Renew(umd_side_coat_solid_small, new G4Box("umd_side_coat_solid_small", len/2 - fCoatingThickness, fCoatingThickness/2, fUMDScintsH/2));

        Renew(umd_fiber_core_solid_small, new G4Tubs("umd_fiber_core_solid_small", 0, fUMDFiberRadius - 2*fCladingThickness, len/2 - fCoatingThickness, 0, 360*CLHEP::deg));
        Renew(umd_fiber_clad1_solid_small, new G4Tubs("umd_fiber_clad1_solid_small", 0, fUMDFiberRadius - fCladingThickness, len/2 - fCoatingThickness, 0, 360*CLHEP::deg));
        Renew(umd_fiber_clad2_solid_small, new G4Tubs("umd_fiber_clad2_solid_small", 0, fUMDFiberRadius, len/2 - fCoatingThickness, 0, 360*CLHEP::deg));

      }

      // Same for all areas
      Renew(umd_back_side_coat_solid, new G4Box("umd_back_side_coat", fCoatingThickness/2, fUMDScintsW/2, fUMDScintsH/2));
      Renew(pixel_solid, new G4Tubs("pixel_solid", 0, fUMDFiberRadius, fPixelL/2, 0, 360*CLHEP::deg));
    }
  }


  void
  G4StationConstruction::CreatePrimitivesMARTA()
  {
    // Al box
    const G4ThreeVector alBoxSize = fAlBoxInnerDimensions + 2*fAlBoxThickness;

    // We sum 0.1 mm of the thickness to the size of the inner part of the tank support to guarantee that
    // if a wall thickness of 0 is chosen the solid subtraction below will be correctly performed.
    const G4ThreeVector alBoxHoleSize = fAlBoxInnerDimensions + 2*0.1*CLHEP::mm * G4ThreeVector(1,1,1);

    G4Box* const alBox_hole = new G4Box("AlBox_hole", alBoxHoleSize[0]/2, alBoxHoleSize[1]/2, alBoxHoleSize[2]/2);

    G4Box* const alBox_out = new G4Box("AlBox_out", alBoxSize[0]/2, alBoxSize[1]/2, alBoxSize[2]/2);

    Renew(AlBox_solid, new G4SubtractionSolid("AlBox_solid", alBox_out, alBox_hole));

    // Double 1mm gap RPC
    const double rpcAcrylicThickness = 0.80*CLHEP::mm;
    const double rpcGlassThickness = 1.85*CLHEP::mm;
    const double rpcGasThickness = 1.0*CLHEP::mm;
    const double rpcHeight = 2*rpcAcrylicThickness + 3*rpcGlassThickness + 2*rpcGasThickness;

    const double rpc_halfX = fRPCsizeX/2;
    const double rpc_halfY = fRPCsizeY/2;
    const double rpc_halfZ = rpcHeight/2;

    const double glass_halfX = rpc_halfX - 0.5*CLHEP::mm;
    const double glass_halfY = rpc_halfY - 0.5*CLHEP::mm;
    const double glass_halfZ = (3*rpcGlassThickness + 2*rpcGasThickness)/2;

    const double gas_halfX = glass_halfX - 0.5*CLHEP::mm;
    const double gas_halfY = glass_halfY - 0.5*CLHEP::mm;
    const double gas_halfZ = rpcGasThickness/2;

    const double pcb_halfX = rpc_halfX;
    const double pcb_halfY = rpc_halfY;
    const double pcb_halfZ = fPCBThickness/2;

    const double spacer_thickness = fSpacerThickness;
    const double spacer_halfX = rpc_halfX;
    const double spacer_halfY = rpc_halfY;
    const double spacer_halfZ = spacer_thickness/2;

    if (rpcHeight + 2*fPCBThickness + spacer_thickness > fAlBoxInnerDimensions[2] + 5*CLHEP::mm) {
      ostringstream errMsg;
      errMsg << "Total height of RPC + 2 x PCB + Spacer = "
             << rpcHeight/CLHEP::mm << " + "
             << 2*fPCBThickness/CLHEP::mm << " + "
             << spacer_thickness/CLHEP::mm << " mm "
             << " greater than inner height of Al box ("
             << fAlBoxInnerDimensions[2]/CLHEP::mm << " mm)!";
      ERROR(errMsg);
      throw utl::OutOfBoundException(errMsg.str());
    }

    Renew(rpc_solid, new G4Box("rpc_solid", rpc_halfX, rpc_halfY, rpc_halfZ));
    Renew(glass_solid, new G4Box("rpc_glass_solid", glass_halfX, glass_halfY, glass_halfZ));
    Renew(gas_solid, new G4Box("rpc_gas_solid", gas_halfX, gas_halfY, gas_halfZ));
    Renew(pcb_solid, new G4Box("pcb_solid", pcb_halfX, pcb_halfY, pcb_halfZ));
    Renew(spacer_solid, new G4Box("spacer_solid", spacer_halfX, spacer_halfY, spacer_halfZ));
  }


  void
  G4StationConstruction::CreateHall()
  {
    // Experimental hall
    Renew(expHall_log, new G4LogicalVolume(expHall_box, Air, "World", 0, 0, 0));
    Renew(expHall_phys, new G4PVPlacement(nullptr, G4ThreeVector(), "World", expHall_log, 0, false, 0));
    expHall_log->SetVisAttributes(G4VisAttributes::Invisible);

    // Ground around the detector
    if (fGroundEnable) {
      // Create sensitive soil to avoid simulation of low energy particles (or other restrictions to speed up cpu time)
      Renew(ground_log, new G4LogicalVolume(ground_solid, Dirt, "ground_solid"));
      Renew(ground_phys, new G4PVPlacement(nullptr, G4ThreeVector(0, 0, -fGroundThickness/2), "ground", ground_log, expHall_phys, false, 0));
      G4VisAttributes gVisatt(G4Colour::Brown());
      if (fG4StationSimulator.fVerbosity<2) ground_log->SetVisAttributes(G4VisAttributes::Invisible);
      // registration
      G4SoilAction* const soilSD = new G4SoilAction("/soil");
      G4SDManager* const sdMan = G4SDManager::GetSDMpointer();
      sdMan->AddNewDetector(soilSD);
      ground_log->SetSensitiveDetector(soilSD);
    }
  }


  void
  G4StationConstruction::AssembleStation()
  {
    G4SDManager* const sdMan = G4SDManager::GetSDMpointer();

    // Water tank
    Renew(tank_log, new G4LogicalVolume(tank_solid, Water, "tank_log", 0, 0, 0));
    Renew(tank_phys, new G4PVPlacement(nullptr, G4ThreeVector(fTankPos_x, fTankPos_y, fTankPos_z + fgTankHalfHeight + fgTankThickness), "tank", tank_log, expHall_phys, false, 0));

    // Top, bottom, side walls
    Renew(top_log, new G4LogicalVolume(top_solid, HDPE, "top_log", 0, 0, 0));
    Renew(bottom_log, new G4LogicalVolume(top_solid, HDPE, "bottom_log", 0, 0, 0));
    Renew(side_log, new G4LogicalVolume(side_solid, HDPE, "side_log", 0, 0, 0));

    Renew(top_phys, new G4PVPlacement(nullptr, G4ThreeVector(fTankPos_x, fTankPos_y, fTankPos_z + 2*fgTankHalfHeight + 1.5*fgTankThickness), "top", top_log, expHall_phys, false, 0));
    Renew(bottom_phys, new G4PVPlacement(nullptr, G4ThreeVector(fTankPos_x, fTankPos_y, fTankPos_z + 0.5*fgTankThickness), "bottom", bottom_log, expHall_phys, false, 0));
    Renew(side_phys, new G4PVPlacement(nullptr, G4ThreeVector(fTankPos_x, fTankPos_y, fTankPos_z + fgTankHalfHeight + fgTankThickness), "side", side_log, expHall_phys, false, 0));

    // the tank liner
    Renew(topsurface, new G4LogicalBorderSurface("topsurface", tank_phys, top_phys, OpLinerSurface));
    Renew(bottomsurface, new G4LogicalBorderSurface("bottomsurface", tank_phys, bottom_phys, OpLinerSurface));
    Renew(sidesurface, new G4LogicalBorderSurface("sidesurface", tank_phys, side_phys, OpLinerSurface));

    // PMTs and such

    // the interior volume of the pmt
    Renew(inner_log, new G4LogicalVolume(inner_solid, Vacuum, "inner_log", 0, 0, 0));  // standard pmt

    for (auto& pmt : fPMTs) {
      const unsigned int id = pmt.fId;
      const int pmtType = pmt.fPMTType;
      if (pmtType)
        continue;

      const auto& pos = pmt.fPMTPosition;
      const double x = pos.x();
      const double y = pos.y();
      const double z = fgTankHalfHeight;

      const G4ThreeVector p(x, y, z);

      ostringstream name;
      name << "inner" << id;
      Renew(pmt.fInner_phys, new G4PVPlacement(nullptr, p, inner_log, name.str().c_str(), tank_log, false, 0));

      // Sensitive surface of pmt faces
      name.str("");
      name << "pmt" << id << "_log";
      Renew(pmt.fPMT_log, new G4LogicalVolume(pmt_solid, Pyrex, name.str().c_str(), 0, 0, 0));

      name.str("");
      name << "pmt" << id;
      Renew(pmt.fPMT_phys, new G4PVPlacement(nullptr, p, pmt.fPMT_log, name.str().c_str(), tank_log, false, 0));

      // Non-sensitive surface of pmt faces
      name.str("");
      name << "pmt" << id << "_log1";
      Renew(pmt.fPMT_log1, new G4LogicalVolume(pmt_solid1, Pyrex1, name.str().c_str(), 0, 0, 0));

      name.str("");
      name << "pmt" << id << "1";
      Renew(pmt.fPMT_phys1, new G4PVPlacement(nullptr, p, pmt.fPMT_log1, name.str().c_str(), tank_log, false, 0));

      // The dome/face interface
      name.str("");
      name << "interface" << id << "_log";
      Renew(pmt.fInterface_log, new G4LogicalVolume(interface_solid, Interface, name.str().c_str(), 0, 0, 0));

      name.str("");
      name << "interface" << id;
      Renew(pmt.fInterface_phys, new G4PVPlacement(nullptr, p, pmt.fInterface_log, name.str().c_str(), tank_log, false, 0));

      // PMT domes
      name.str("");
      name << "dome" << id << "_log";
      Renew(pmt.fDome_log, new G4LogicalVolume(dome_solid, Lucite, name.str().c_str(), 0, 0, 0));

      name.str("");
      name << "dome" << id;
      Renew(pmt.fDome_phys, new G4PVPlacement(nullptr, p, pmt.fDome_log, name.str().c_str(), tank_log, false, 0));

      // Register the PMT's as sensitive detectors
      name.str("");
      name << "/Tank/pmt" << id;
      G4TankPMTAction* const pmtSD = new G4TankPMTAction(name.str().c_str(), pmt.fId);
      sdMan->AddNewDetector(pmtSD);
      pmt.fPMT_log->SetSensitiveDetector(pmtSD);
    }

    // Solar panel
    if (fSolarPanelEnable) {
      Renew(solarPanel_solid, new G4Box("SolarPanel", fSolarPanelThickness/2, fSolarPanelWidth/2, fSolarPanelLength/2));
      Renew(solarPanel_log, new G4LogicalVolume(solarPanel_solid, Al, "solarPanel_log"));  //  MAY BE CHANGED

      G4RotationMatrix* const solarPanelRot = new G4RotationMatrix();
      solarPanelRot->rotate(M_PI/2, G4ThreeVector(0, 0, 1));
      solarPanelRot->rotate(fSolarPanelTiltAngle, G4ThreeVector(0, 1, 0));

      Renew(solarPanel_phys, new G4PVPlacement(solarPanelRot, G4ThreeVector(fSolarPanelX, fSolarPanelY, fSolarPanelZ), "solarPanel", solarPanel_log, expHall_phys, false, 0));
    }

    // Chunk of Aluminum representing electronics box
    if (fElecBoxEnable) {
      Renew(elecBox_solid, new G4Box("ElecBox", fElecBoxThickness/2, fElecBoxWidth/2, fElecBoxLength/2));
      Renew(elecBox_log, new G4LogicalVolume(elecBox_solid, Al, "elecBox_log"));  // MAY BE CHANGED

      G4RotationMatrix* const elecBoxRot = new G4RotationMatrix();
      elecBoxRot->rotate(M_PI/2, G4ThreeVector(0, 0, 1));
      elecBoxRot->rotate(fElecBoxTiltAngle, G4ThreeVector(0, 1, 0));

      Renew(elecBox_phys, new G4PVPlacement(elecBoxRot, G4ThreeVector(fElecBoxX, fElecBoxY, fElecBoxZ), "elecBox", elecBox_log, expHall_phys, false, 0));
    }

    // Elements specific to non-zero hardware versions
    if (fHasSmallPMT) {

      // Half the tank's height is used for the z-coordinate since the PMT volumes
      // are placed in the logical volume of the water, not the world volume

      // only for spmt but loop could be used to get scint pmt too
      for (auto& pmt : fPMTs) {

        const unsigned int id = pmt.fId;
        const int pmtType = pmt.fPMTType;
        if (pmtType != 1)  // when scintillator we will put a || 2
          continue;

        const auto& pos = pmt.fPMTPosition;
        const double x = pos.x();
        const double y = pos.y();
        const double z = -fgTankThickness / 2;
        ostringstream name;

        // Sensitive surface of spmt face
        name.str("");
        name << "pmt" << id << "_log";
        Renew(pmt.fPMT_log, new G4LogicalVolume(s_pmt_solid, Pyrex, name.str().c_str(), 0, 0, 0));

        name.str("");
        name << "pmt" << id;
        Renew(pmt.fPMT_phys, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z + fPmtZ_SPMT), pmt.fPMT_log, name.str().c_str(), top_log, false, 0));

        // Non-sensitive surface of spmt face
        name.str("");
        name << "pmt" << id << "_log1";
        Renew(pmt.fPMT_log1, new G4LogicalVolume(s_pmt_solid1, Pyrex1,  name.str().c_str(), 0, 0, 0));

        name.str("");
        name << "pmt" << id << "1";
        Renew(pmt.fPMT_phys1, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z + fPmtZ_SPMT), pmt.fPMT_log1, name.str().c_str(), top_log, false, 0));

        // The dome/face interface
        name.str("");
        name << "interface" << id << "_log";
        Renew(pmt.fInterface_log, new G4LogicalVolume(s_interface_solid, Interface, name.str().c_str(), 0, 0, 0));

        name.str("");
        name << "interface" << id;
        Renew(pmt.fInterface_phys, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z + fInterfaceZ_SPMT), pmt.fInterface_log, name.str().c_str(), top_log, false, 0));

        // sPMT dome
        name.str("");
        name << "dome" << id << "_log";
        Renew(pmt.fDome_log, new G4LogicalVolume(s_dome_solid, Lucite, name.str().c_str(), 0, 0, 0));

        name.str("");
        name << "dome" << id;
        Renew(pmt.fDome_phys, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z + fDomeZ_SPMT), pmt.fDome_log, name.str().c_str(), top_log, false, 0));

        // Register the sPMT as sensitive detector
        name.str("");
        name << "/Tank/pmt" << id;
        G4TankPMTAction* const pmtSD = new G4TankPMTAction(name.str().c_str(), pmt.fId);
        sdMan->AddNewDetector(pmtSD);
        pmt.fPMT_log->SetSensitiveDetector(pmtSD);

      }
    }

    if (fHasScintillator) {
      /* The scintillator origin is directly between two scintillator halves
        and mid-way through the thickness of the polystyrene scintillator material.
        The locations of all volumes are defined here are defined relative to this
        position. See SModelConfig for a schematic and more details. */
      const double x = fScintillatorPosition.x();
      const double y = fScintillatorPosition.y();
      const double z = fScintillatorPosition.z();
      // aluminum casing
      const double casing_z = z + fScintillatorHousingThickness/2 - 2*fScintillatorCasingThickness - fScintillatorSandwichPanelThickness - fScintillatorBarThickness/2;
      Renew(scin_casing_solid, new G4Box("scin_casing_solid", fScintillatorCasingThickness + fScintillatorHousingLength/2, fScintillatorCasingThickness + fScintillatorHousingWidth/2, fScintillatorCasingThickness + fScintillatorHousingThickness/2));
      Renew(scin_casing_log, new G4LogicalVolume(scin_casing_solid, Al, "scin_casing_log"));
      Renew(scin_casing_phys, new G4PVPlacement(nullptr, G4ThreeVector(x, y, casing_z), "scin_casing", scin_casing_log, expHall_phys, false, 0));
      // sandwich panel (mother volume: casing)
      const double sandwich_z = fScintillatorSandwichPanelThickness/2 - fScintillatorHousingThickness/2;
      Renew(scin_sandwich_solid, new G4Box("scin_sandwich_solid", fScintillatorHousingLength/2, fScintillatorHousingWidth/2, fScintillatorSandwichPanelThickness/2));
      Renew(scin_sandwich_log, new G4LogicalVolume(scin_sandwich_solid, ExtrudedPolystyreneFoam, "scin_sandwich_log"));
      Renew(scin_sandwich_phys, new G4PVPlacement(nullptr, G4ThreeVector(0, 0, sandwich_z), scin_sandwich_log, "scin_sandwich", scin_casing_log, true, 0));
      // styrofoam (mother volume: casing)
      const double styro_thickness = fScintillatorHousingThickness - fScintillatorCasingThickness - fScintillatorSandwichPanelThickness;
      const double styro_z = fScintillatorHousingThickness/2 - styro_thickness/2;
      Renew(scin_styro_solid, new G4Box("scin_sandwich_solid", fScintillatorHousingLength/2, fScintillatorHousingWidth/2, styro_thickness/2));
      Renew(scin_styro_log, new G4LogicalVolume(scin_styro_solid, ExpandedPolystyreneFoam, "scin_styro_log"));
      Renew(scin_styro_phys, new G4PVPlacement(nullptr, G4ThreeVector(0, 0, styro_z), scin_styro_log, "scin_styro", scin_casing_log, true, 0));
      // scintillator (mother volume: styrofoam)
      const double bars_z = fScintillatorCasingThickness + fScintillatorBarThickness/2 - styro_thickness/2;
      Renew(scin_solid, new G4Box("scin_solid", fScintillatorBarLength/2, fScintillatorBarWidth*fNScintillatorBars/4, fScintillatorBarThickness/2));
      Renew(scin_log, new G4LogicalVolume(scin_solid, Polystyrene, "scin_log"));
      Renew(scin_phys1, new G4PVPlacement(nullptr, G4ThreeVector(-(fScintillatorGap + fScintillatorBarLength)/2, 0, bars_z), scin_log, "scin_bars", scin_styro_log, true, 0));
      Renew(scin_phys2, new G4PVPlacement(nullptr, G4ThreeVector((fScintillatorGap + fScintillatorBarLength)/2, 0, bars_z), scin_log, "scin_bars", scin_styro_log, true, 1));
      // registration
      G4ScintillatorAction* const scintillatorSD = new G4ScintillatorAction("/Scintillator/scin_bars");
      sdMan->AddNewDetector(scintillatorSD);
      scin_log->SetSensitiveDetector(scintillatorSD);
      // roof
      const double roof_z = z - 0.5*fScintillatorBarThickness - fScintillatorCasingThickness - fScintillatorSandwichPanelThickness + fScintillatorHousingThickness + fScintillatorRoofOffset + 0.5*fScintillatorRoofThickness;
      Renew(scin_roof_solid, new G4Box("scin_roof_solid", fScintillatorCasingThickness + fScintillatorHousingLength/2, fScintillatorCasingThickness + fScintillatorHousingWidth/2, fScintillatorRoofThickness/2));
      Renew(scin_roof_log, new G4LogicalVolume(scin_roof_solid, Al, "scin_roof_log"));
      Renew(scin_roof_phys, new G4PVPlacement(nullptr, G4ThreeVector(x, y, roof_z), "scin_roof", scin_roof_log, expHall_phys, false, 0));
    }

    // Creates Amiga underground primitives
    if (fUMDEnabled) {
      const auto& dStation = G4StationSimulator::fgCurrent.GetDetectorStation();
      if (!dStation.ExistsAssociatedCounter()) {
        ostringstream err;
        err << "Trying to simulate the UMD in a WCD (Id " << dStation.GetId() << ") that has not an underground AMIGA scintillator associated!";
        //throw utl::NonExistentComponentException(err.str());
        WARNING(err);
      }
      AssembleUMD();
    }

    // Assembles MARTA RPC
    if (fMARTAEnabled)
      AssembleMARTA();
  }


  void
  G4StationConstruction::AssembleUMD()
  {
    G4SDManager* const sdMan = G4SDManager::GetSDMpointer();

    // Casings
    Renew(umd_casing_log_large, new G4LogicalVolume(umd_casing_solid_large, PVC, "umd_casing_log_large", 0, 0, 0));
    Renew(umd_casing_log_small, new G4LogicalVolume(umd_casing_solid_small, PVC, "umd_casing_log_small", 0, 0, 0));

    // Scintillator Strips
    Renew(umd_strip_log_large, new G4LogicalVolume(umd_strip_solid_large, ScintMat, "umd_strip_log_large", 0, 0, 0));
    const G4VisAttributes blue(G4Colour::Blue());
    umd_strip_log_large->SetVisAttributes(blue);

    // registration
    G4UMDScintStripAction* const umdScintStripLargeSD = new G4UMDScintStripAction("/UMDStrip/umd_strip_large");
    sdMan->AddNewDetector(umdScintStripLargeSD);
    umd_strip_log_large->SetSensitiveDetector(umdScintStripLargeSD);

    Renew(umd_strip_log_small, new G4LogicalVolume(umd_strip_solid_small, ScintMat, "umd_strip_log_small", 0, 0, 0));
    umd_strip_log_small->SetVisAttributes(blue);
    // registration
    G4UMDScintStripAction* const umdScintStripSmallSD = new G4UMDScintStripAction("/UMDStrip/umd_strip_small");
    sdMan->AddNewDetector(umdScintStripSmallSD);
    umd_strip_log_small->SetSensitiveDetector(umdScintStripSmallSD);

    // Reflective TiO2 coating (surface)
    Renew(umd_top_bot_coat_log_large, new G4LogicalVolume(umd_top_bot_coat_solid_large, ScintCoating, "umd_top_bot_coat_log_large", 0, 0, 0));
    const G4VisAttributes black(G4Colour::Black());
    umd_top_bot_coat_log_large->SetVisAttributes(black);

    Renew(umd_side_coat_log_large, new G4LogicalVolume(umd_side_coat_solid_large, ScintCoating, "umd_side_coat_log_large", 0, 0, 0));
    umd_side_coat_log_large->SetVisAttributes(black);

    Renew(umd_top_bot_coat_log_small, new G4LogicalVolume(umd_top_bot_coat_solid_small, ScintCoating, "umd_top_bot_coat_log_small", 0, 0, 0));
    umd_top_bot_coat_log_small->SetVisAttributes(black);

    Renew(umd_side_coat_log_small, new G4LogicalVolume(umd_side_coat_solid_small, ScintCoating, "umd_side_coat_log_small", 0, 0, 0));
    umd_side_coat_log_small->SetVisAttributes(black);

    Renew(umd_back_side_coat_log, new G4LogicalVolume(umd_back_side_coat_solid, ScintCoating, "umd_back_side_coat_log", 0, 0, 0));
    const G4VisAttributes brown(G4Colour::Brown());
    umd_back_side_coat_log->SetVisAttributes(brown);

    const G4VisAttributes green(G4Colour::Green());
    if (!fG4StationSimulator.fUMDFastMode || fG4StationSimulator.fGeoVisOn || fG4StationSimulator.fTrajVisOn) {
      // if detector needed for full simulation mode (optical photons)
      // Fibers (neglects OnManifold differences for the moment)
      // Core
      Renew(umd_fiber_core_log_large, new G4LogicalVolume(umd_fiber_core_solid_large, PMMA, "umd_fiber_core_log_large", 0, 0, 0));
      umd_fiber_core_log_large->SetVisAttributes(green);
      // Clad1 (internal)
      Renew(umd_fiber_clad1_log_large, new G4LogicalVolume(umd_fiber_clad1_solid_large, Pethylene, "umd_fiber_clad1_log_large", 0, 0, 0));
      umd_fiber_clad1_log_large->SetVisAttributes(blue);
      // Clad2 (external)
      Renew(umd_fiber_clad2_log_large, new G4LogicalVolume(umd_fiber_clad2_solid_large, FPethylene, "umd_fiber_clad2_log_large", 0, 0, 0));
      umd_fiber_clad2_log_large->SetVisAttributes(blue);
      // Core
      Renew(umd_fiber_core_log_small, new G4LogicalVolume(umd_fiber_core_solid_small, PMMA, "umd_fiber_core_log_small", 0, 0, 0));
      umd_fiber_core_log_small->SetVisAttributes(green);
      // Clad1 (internal)
      Renew(umd_fiber_clad1_log_small, new G4LogicalVolume(umd_fiber_clad1_solid_small, Pethylene, "umd_fiber_clad1_log_small", 0, 0, 0));
      umd_fiber_clad1_log_small->SetVisAttributes(blue);
      // Clad2 (external)
      Renew(umd_fiber_clad2_log_small, new G4LogicalVolume(umd_fiber_clad2_solid_small, FPethylene, "umd_fiber_clad2_log_small", 0, 0, 0));
      umd_fiber_clad2_log_small->SetVisAttributes(blue);
      // Pixels (for sensitive detector)
      Renew(pixel_log, new G4LogicalVolume(pixel_solid, Pyrex, "pixel_log", 0, 0, 0));
      pixel_log->SetVisAttributes(black);
      // registration
      G4UMDPixelAction* const umdPixelSD = new G4UMDPixelAction("/UMDPixel/umd_pixel");
      sdMan->AddNewDetector(umdPixelSD);
      pixel_log->SetSensitiveDetector(umdPixelSD);
      // Skin Properties of coating
      SetScintCoatingProperties();
    }

    // Get the muon detector
    const mdet::MDetector& mDet = det::Detector::GetInstance().GetMDetector();
    const auto& dStation = G4StationSimulator::fgCurrent.GetDetectorStation();
    const mdet::Counter& dCounter = mDet.GetCounter(dStation.GetAssociatedCounterId());

    ostringstream name;

    // Loop over muon modules (may only differs in bar lengths)
    for (auto mIt = dCounter.ModulesBegin(); mIt != dCounter.ModulesEnd(); ++mIt) {

      const mdet::Module::AreaKind area = mIt->GetAreaKind();
      const mdet::Module& mod = *mIt;

      std::cout << (area == mdet::Module::AreaKind::eLarge ? "10m2" : "5m2") << std::endl;

      const auto cs = dCounter.GetLocalCoordinateSystem();
      const auto& modPos = mod.GetPosition();
      const double x = modPos.GetX(cs) * kDistanceToG4;
      const double y = modPos.GetY(cs) * kDistanceToG4;
      const double z = modPos.GetZ(cs) * kDistanceToG4 + fGroundThickness/2;
      // zero here is measured from the center of ground cylinder
      const G4ThreeVector pos(x, y, z);

      // Casing rotation
      G4RotationMatrix* const rotMod = new G4RotationMatrix();
      rotMod->rotateZ(-mod.GetPhi0()); // The minus sign is to take into account that positive angle means clockwise rotation for G4
      // while it is counterclockwise for Offline.

      name.str("");
      name << "module_" << mod.GetId() << '_' << dCounter.GetId();

      if (area == mdet::Module::AreaKind::eLarge)
        NotRenew(umd_casing_phys_large, new G4PVPlacement(rotMod, pos, name.str() + "_large", umd_casing_log_large, ground_phys, false, mod.GetId()));
      else
        NotRenew(umd_casing_phys_small, new G4PVPlacement(rotMod, pos, name.str() + "_small", umd_casing_log_small, ground_phys, false, mod.GetId()));

    }

    // As mother volumes (mod_phys in this case) has multiple copies, child volumes
    // (umd_strip_phys) do not need to be replicated
    // Loop only over area kinds
    for (const auto& ml : fModAreaLenghts) {
      const auto area = ml.first;
      std::cout << "Creating G4 physcal volume for " << (area == mdet::Module::AreaKind::eLarge ? "10m2" : "5m2") << " module" << std::endl;

      auto mIt = dCounter.ModulesBegin();
      for ( ; mIt != dCounter.ModulesEnd(); ++mIt)
        if (mIt->GetAreaKind() == area)
          break;
      if (mIt == dCounter.ModulesEnd())
        continue;

      const mdet::Module& mod = *mIt;
      const double shortestFiber = mod.GetShortestFiber();

      // Scintillators
      Renew(umd_strip_phys_large, static_cast<decltype(umd_strip_phys_large)>(nullptr));
      Renew(umd_strip_phys_small, static_cast<decltype(umd_strip_phys_small)>(nullptr));

      for (auto scIt = mod.ScintillatorsBegin(); scIt != mod.ScintillatorsEnd(); ++scIt) {

        const mdet::Scintillator& sc = *scIt;
        const mdet::Fiber& fiber = mod.GetFiberFor(sc);

        const auto cs = mod.GetLocalCoordinateSystem();
        const double x = sc.GetPosition().GetX(cs) * kDistanceToG4;
        const double y = sc.GetPosition().GetY(cs) * kDistanceToG4;
        const double z = sc.GetPosition().GetZ(cs) * kDistanceToG4;

        // This cannot be <= 0
        // removes the shortest length. Only relative difference is important
        const double onManiFold = fiber.GetOnManifoldLength() * kDistanceToG4 - shortestFiber * kDistanceToG4 + 2*fPixelL;
        const double dx = (x > 0 ? -1 : 1) * (fModAreaLenghts[area]/2 + onManiFold/2 - fCoatingThickness);

        NotRenew(extra_fiber_solid, new G4Tubs("extra_fiber_solid", 0, fUMDFiberRadius - 2*fCladingThickness, onManiFold/2, 0, 360*CLHEP::deg));
        NotRenew(extra_fiber_log, new G4LogicalVolume(extra_fiber_solid, PMMA, "extra_fiber_log", 0, 0, 0));
        extra_fiber_log->SetVisAttributes(green);

        NotRenew(extra_clad1_solid, new G4Tubs("extra_clad1_solid", 0, fUMDFiberRadius - fCladingThickness, onManiFold/2, 0, 360*CLHEP::deg));
        NotRenew(extra_clad1_log, new G4LogicalVolume(extra_clad1_solid, Pethylene, "extra_clad1_log", 0, 0, 0));
        extra_clad1_log->SetVisAttributes(blue);

        NotRenew(extra_clad2_solid, new G4Tubs("extra_clad2_solid", 0, fUMDFiberRadius, onManiFold/2, 0, 360*CLHEP::deg));
        NotRenew(extra_clad2_log, new G4LogicalVolume(extra_clad2_solid, FPethylene, "extra_clad2_log", 0, 0, 0));
        extra_clad2_log->SetVisAttributes(blue);

        if (area == mdet::Module::AreaKind::eLarge) {

          name.str("");
          name << "umd_strip_" << sc.GetId();
          NotRenew(umd_strip_phys_large, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z), name.str() + "_large", umd_strip_log_large, umd_casing_phys_large, false, sc.GetId()));

          if (!fG4StationSimulator.fUMDFastMode || fG4StationSimulator.fGeoVisOn || fG4StationSimulator.fTrajVisOn) {

            string nm = "umd_top_coat";
            NotRenew(umd_top_coat_phys_large, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z + fUMDScintsH/2 - fCoatingThickness/2), nm + "_large", umd_top_bot_coat_log_large, umd_casing_phys_large, false,  0));

            nm = "umd_bot_coat";
            NotRenew(umd_bot_coat_phys_large, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z - fUMDScintsH/2 + fCoatingThickness/2), nm + "_large", umd_top_bot_coat_log_large, umd_casing_phys_large, false,  0));

            nm = "umd_side1_coat";
            NotRenew(umd_side1_coat_phys_large, new G4PVPlacement(nullptr, G4ThreeVector(x, y + fUMDScintsW/2 - fCoatingThickness/2, z), nm + "_large", umd_side_coat_log_large, umd_casing_phys_large, false, 0));

            nm = "umd_side2_coat";
            NotRenew(umd_side2_coat_phys_large, new G4PVPlacement(nullptr, G4ThreeVector(x, y - fUMDScintsW/2 + fCoatingThickness/2, z), nm + "_large", umd_side_coat_log_large, umd_casing_phys_large, false, 0));

            nm = "umd_back_side_coat";
            const double xx = (x > 0 ? 1 : -1) * (fModAreaLenghts[area]/2 - fCoatingThickness/2);
            NotRenew(umd_back_side_coat_phys, new G4PVPlacement(nullptr, G4ThreeVector(x + xx, y, z), nm + "_large", umd_back_side_coat_log, umd_casing_phys_large, false, 0));

            G4RotationMatrix* const erot = new G4RotationMatrix();
            erot->rotateY(M_PI/2);
            const double z = fUMDScintsH/2 - fCoatingThickness - fUMDFiberRadius - fCladingThickness;  // at top of the strip

            NotRenew(extra_clad2_phy, new G4PVPlacement(erot, G4ThreeVector(x+dx, y, z), "extra_clad2", extra_clad2_log, umd_casing_phys_large, false, 0));

            NotRenew(extra_clad1_phy, new G4PVPlacement(nullptr, G4ThreeVector(), "extra_clad1", extra_clad1_log, extra_clad2_phy, false, 0));

            NotRenew(extra_fiber_phy, new G4PVPlacement(nullptr, G4ThreeVector(), "extra_fiber", extra_fiber_log, extra_clad1_phy, false, 0));

            G4RotationMatrix* const prot = new G4RotationMatrix();
            prot->rotateY(M_PI/2);

            name.str("");
            name << "mod_pix_" << sc.GetId();
            NotRenew(pixel_phy, new G4PVPlacement(prot, G4ThreeVector(x + dx + (x > 0 ? -1 : 1)*(onManiFold/2 + fPixelL/2), y, z), name.str(), pixel_log, umd_casing_phys_large, false, sc.GetId()));
          }

        } else {

          name.str("");
          name << "umd_strip_" << sc.GetId();
          NotRenew(umd_strip_phys_small, new G4PVPlacement(nullptr, G4ThreeVector(x,y,z), name.str() + "_small", umd_strip_log_small, umd_casing_phys_small, false, sc.GetId()));

          if (!fG4StationSimulator.fUMDFastMode || fG4StationSimulator.fGeoVisOn || fG4StationSimulator.fTrajVisOn) {
            NotRenew(umd_top_coat_phys_small, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z + fUMDScintsH/2 - fCoatingThickness/2), "umd_top_coat_small", umd_top_bot_coat_log_small, umd_casing_phys_small, false, 0));

            NotRenew(umd_bot_coat_phys_small, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z - fUMDScintsH/2 + fCoatingThickness/2), "umd_bot_coat_small", umd_top_bot_coat_log_small, umd_casing_phys_small, false, 0));

            NotRenew(umd_side1_coat_phys_small, new G4PVPlacement(nullptr, G4ThreeVector(x, y + fUMDScintsW/2 - fCoatingThickness/2, z), "umd_side1_coat_small", umd_side_coat_log_small, umd_casing_phys_small, false, 0));

            NotRenew(umd_side2_coat_phys_small, new G4PVPlacement(nullptr, G4ThreeVector(x, y - fUMDScintsW/2 + fCoatingThickness/2, z), "umd_side2_coat_small", umd_side_coat_log_small, umd_casing_phys_small, false, 0));

            const double xx = (x > 0 ? 1 : -1) * (fModAreaLenghts[area]/2 - fCoatingThickness/2);
            NotRenew(umd_back_side_coat_phys, new G4PVPlacement(nullptr, G4ThreeVector(x + xx, y, z), "umd_back_side_coat_small", umd_back_side_coat_log, umd_casing_phys_small, false, 0));

            G4RotationMatrix* const erot = new G4RotationMatrix();
            erot->rotateY(M_PI/2);
            const double z = fUMDScintsH/2 - fCoatingThickness - fUMDFiberRadius - fCladingThickness;  // at top of the strip

            NotRenew(extra_clad2_phy, new G4PVPlacement(erot, G4ThreeVector(x + dx, y, z), "extra_clad2", extra_clad2_log, umd_casing_phys_small, false, 0));

            NotRenew(extra_clad1_phy, new G4PVPlacement(nullptr, G4ThreeVector(), "extra_clad1", extra_clad1_log, extra_clad2_phy, false, 0));

            NotRenew(extra_fiber_phy, new G4PVPlacement(nullptr, G4ThreeVector(), "extra_fiber", extra_fiber_log, extra_clad1_phy, false, 0));

            G4RotationMatrix* const prot = new G4RotationMatrix();
            prot->rotateY(M_PI/2);

            name.str("");
            name << "mod_pix_" << sc.GetId();
            NotRenew(pixel_phy, new G4PVPlacement(prot, G4ThreeVector(x + dx + (x > 0 ? -1 : 1)*(onManiFold/2 + fPixelL/2), y, z), name.str(), pixel_log, umd_casing_phys_small, false, sc.GetId()));
          }

        }

      }

    }

    if (!fG4StationSimulator.fUMDFastMode || fG4StationSimulator.fGeoVisOn || fG4StationSimulator.fTrajVisOn) {
      // As mother volume (umd_strip_phys) has multiple copies, child volumes (umd_strip_log) do not need to be replicated
      if (umd_strip_phys_large) {

        double x = 0;
        double y = 0;
        double z = fUMDScintsH/2 - fCoatingThickness - fUMDFiberRadius - fCladingThickness;  // at top of the strip
        G4RotationMatrix* const frotMod_large = new G4RotationMatrix();
        frotMod_large->rotateY(M_PI/2);
        NotRenew(umd_fiber_clad2_phys_large, new G4PVPlacement(frotMod_large, G4ThreeVector(x, y, z), "umd_fiber_clad2_large", umd_fiber_clad2_log_large, umd_strip_phys_large, false, 0));

        x = 0;
        y = 0;
        z = 0;
        NotRenew(umd_fiber_clad1_phys_large, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z), "umd_fiber_clad1_large", umd_fiber_clad1_log_large, umd_fiber_clad2_phys_large, false,  0));

        x = 0;
        y = 0;
        z = 0;
        NotRenew(umd_fiber_core_phys_large, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z), "umd_fiber_core_large", umd_fiber_core_log_large, umd_fiber_clad1_phys_large, false, 0));

      }

      if (umd_strip_phys_small) {

        double x = 0;
        double y = 0;
        double z = fUMDScintsH/2 - fCoatingThickness - fUMDFiberRadius - fCladingThickness;  // at top of the strip

        G4RotationMatrix* const frotMod_small = new G4RotationMatrix();
        frotMod_small->rotateY(M_PI/2);
        NotRenew(umd_fiber_clad2_phys_small, new G4PVPlacement(frotMod_small, G4ThreeVector(x, y, z), "umd_fiber_clad2_small", umd_fiber_clad2_log_small, umd_strip_phys_small, false, 0));

        x = 0;
        y = 0;
        z = 0;
        NotRenew(umd_fiber_clad1_phys_small, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z), "umd_fiber_clad1_small", umd_fiber_clad1_log_small, umd_fiber_clad2_phys_small, false,  0));

        x = 0;
        y = 0;
        z = 0;
        NotRenew(umd_fiber_core_phys_small, new G4PVPlacement(nullptr, G4ThreeVector(x, y, z), "umd_fiber_core_small", umd_fiber_core_log_small, umd_fiber_clad1_phys_small, false, 0));

      }
    }
  }


  void
  G4StationConstruction::AssembleMARTA()
  {
    if (fMakeTankSupport)
      AssembleTankSupport();

    // Aluminium box
    const G4ThreeVector alBoxSize = fAlBoxInnerDimensions + 2*fAlBoxThickness;

    Renew(al_box_log, new G4LogicalVolume(AlBox_solid, Al, "Al_box_log", 0, 0, 0));

    // RPC
    const double rpcAcrylicThickness = 0.80*CLHEP::mm;
    const double rpcGlassThickness = 1.85*CLHEP::mm;
    const double rpcGasThickness = 1.0*CLHEP::mm;
    const double rpcHeight = 2.0*rpcAcrylicThickness + 3.0*rpcGlassThickness + 2.0*rpcGasThickness;

    const double pcb_halfZ = fPCBThickness/2;

    const double spacer_thickness = fSpacerThickness;
    const double spacer_halfZ = spacer_thickness/2;

    Renew(pcb_log, new G4LogicalVolume(pcb_solid, fG4_Bakelite, "pcb_log", 0, 0, 0));
    Renew(spacer_log, new G4LogicalVolume(spacer_solid, ExpandedPolystyreneFoam, "spacer_log", 0, 0, 0));
    Renew(rpc_log, new G4LogicalVolume(rpc_solid, fG4_Acrylic, "rpc_log", 0, 0, 0));
    Renew(glass_log, new G4LogicalVolume(glass_solid, fSoda_lime_glass, "rpc_glass_log", 0, 0, 0));
    Renew(gas_log, new G4LogicalVolume(gas_solid, fR134a, "rpc_gas_log", 0, 0, 0));

    // Position the two gas gaps inside the glass

    const G4ThreeVector gasPos1(0, 0, (rpcGlassThickness + rpcGasThickness)/2);
    const G4ThreeVector gasPos2 = -gasPos1;

    new G4PVPlacement(nullptr, gasPos1, gas_log, "rpc_gas", glass_log, false, 0);
    new G4PVPlacement(nullptr, gasPos2, gas_log, "rpc_gas", glass_log, false, 1);

    // Position the glass inside the acrylic

    const G4ThreeVector glassPos;
    new G4PVPlacement(nullptr, glassPos, glass_log, "rpc_glass", rpc_log, false, 0);

    // end RPC

    // Position tank support in the world volume;
    // in the G4 coordinate system z=0 is the outer surface of the tank bottom

    // Position different components in the world volume
    for (unsigned int i = 0 ; i < fRPCPositions.size() ; ++i) {

      // Position Al box
      // RPC position are given in the SD coordinate system

      const G4ThreeVector alBox_pos =
        fRPCPositions[i] + fgTankCenter + G4ThreeVector(0, 0, alBoxSize[2]/2) + G4ThreeVector(0, 0, 0.1*CLHEP::cm);  // for tolerance

      // Position RPC inside concrete precast

      const G4ThreeVector rpc_pos =
        alBox_pos - G4ThreeVector(0, 0, fAlBoxInnerDimensions[2]/2) + G4ThreeVector(0, 0, rpcHeight/2) + G4ThreeVector(0, 0, 0.1*CLHEP::cm);  // for tolerance

      // Position PCB + Spacer

      const G4ThreeVector pcb_1_pos =
        rpc_pos + G4ThreeVector(0, 0, (rpcHeight/2 + pcb_halfZ)) + G4ThreeVector(0, 0, 0.1*CLHEP::mm);  // for tolerance

      const G4ThreeVector spacer_pos =
        pcb_1_pos + G4ThreeVector(0, 0, (pcb_halfZ + spacer_halfZ)) + G4ThreeVector(0, 0, 0.1*CLHEP::mm);  // for tolerance

      const G4ThreeVector pcb_2_pos =
        spacer_pos + G4ThreeVector(0, 0, (spacer_halfZ + pcb_halfZ)) + G4ThreeVector(0, 0, 0.1*CLHEP::mm);  // for tolerance

      // Rotation of RPC chamber
      G4RotationMatrix* const r = new G4RotationMatrix();
      r->rotate(fRPCRotations[i], G4ThreeVector(0, 0, 1));

      new G4PVPlacement(r, alBox_pos, "Al_box",al_box_log, expHall_phys, false, i+1);
      new G4PVPlacement(r, rpc_pos, "rpc", rpc_log, expHall_phys, false, i+1);
      new G4PVPlacement(r, pcb_1_pos, "pcb1", pcb_log, expHall_phys, false, i+1);
      new G4PVPlacement(r, pcb_2_pos, "pcb2", pcb_log, expHall_phys, false, i+1);
      new G4PVPlacement(r, spacer_pos, "spacer", spacer_log, expHall_phys, false, i+1);
    }

    // Visualization attributes

    G4VisAttributes* const rpcVisAtt = new G4VisAttributes(G4Colour::Red());
    rpcVisAtt->SetForceWireframe(true);
    rpc_log->SetVisAttributes(rpcVisAtt);

    G4VisAttributes* const glassVisAtt = new G4VisAttributes(G4Colour::Yellow());
    glassVisAtt->SetForceWireframe(true);
    glass_log->SetVisAttributes(glassVisAtt);

    G4VisAttributes* const gasVisAtt = new G4VisAttributes(G4Colour::Cyan());
    gasVisAtt->SetForceWireframe(true);
    gas_log->SetVisAttributes(gasVisAtt);

    G4VisAttributes* const alBoxVisAtt = new G4VisAttributes(G4Colour::Grey());
    alBoxVisAtt->SetForceSolid(true);
    alBoxVisAtt->SetVisibility(true);
    alBoxVisAtt->SetForceWireframe(true);
    al_box_log->SetVisAttributes(alBoxVisAtt);

    G4VisAttributes* const pcbVisAtt = new G4VisAttributes(G4Colour::Green());
    pcbVisAtt->SetForceSolid(true);
    //pcbVisAtt->SetForceWireframe(true);
    pcb_log->SetVisAttributes(pcbVisAtt);

    spacer_log->SetVisAttributes(new G4VisAttributes(G4Colour::White()));
  }


  // Creates the MARTA RPC concrete support
  void
  G4StationConstruction::AssembleTankSupport()
  {
    const G4ThreeVector tankSupport_pos =
      G4ThreeVector(0, 0, fTankPos_z - (fTankSupportTopSlabDimensions[2]/2 + 0.1*CLHEP::cm));  // for tolerance

    G4Box* const topSlab = new G4Box("topSlab", fTankSupportTopSlabDimensions[0]/2, fTankSupportTopSlabDimensions[1]/2, fTankSupportTopSlabDimensions[2]/2);

    G4Box* const centralFoot = new G4Box("centralFoot", fTankSupportCentralFootDimensions[0]/2, fTankSupportCentralFootDimensions[1]/2, fTankSupportCentralFootDimensions[2]/2);

    G4Box* const centralFootBase = new G4Box("centralFootBase", fTankSupportCentralFootBaseDimensions[0]/2, fTankSupportCentralFootBaseDimensions[1]/2, fTankSupportCentralFootBaseDimensions[2]/2);

    const G4ThreeVector trans1(0, 0, -(fTankSupportCentralFootDimensions[2]/2 + fTankSupportCentralFootBaseDimensions[2]/2 - 0.1*CLHEP::mm));

    G4BooleanSolid* const centralFoot_solid = new G4UnionSolid("centralFoot_solid", centralFoot, centralFootBase, 0, trans1);

    // Build now the outer foot
    G4Box* const outerFoot = new G4Box("outerFoot", fTankSupportOuterFootDimensions[0]/2, fTankSupportOuterFootDimensions[1]/2, fTankSupportOuterFootDimensions[2]/2);

    G4Box* const outerFootBase = new G4Box("outerFootBase", fTankSupportOuterFootBaseDimensions[0]/2, fTankSupportOuterFootBaseDimensions[1]/2, fTankSupportOuterFootBaseDimensions[2]/2);

    const G4ThreeVector trans2(0, 0, -(fTankSupportOuterFootDimensions[2]/2 + fTankSupportOuterFootBaseDimensions[2]/2 - 0.1*CLHEP::mm));

    G4BooleanSolid* const outerFoot_solid = new G4UnionSolid("outerFoot_solid", outerFoot, outerFootBase, 0, trans2);

    // Now join everything
    const G4ThreeVector trans3(0, 0, -(fTankSupportTopSlabDimensions[2]/2 + fTankSupportCentralFootDimensions[2]/2 - 0.1*CLHEP::mm));

    G4BooleanSolid* const tankSupport_solid_0 = new G4UnionSolid("tank_support_0", topSlab, centralFoot_solid, 0, trans3);

    const G4ThreeVector trans4(-fTankSupportOuterFootDistanceToCenter, 0, -(fTankSupportTopSlabDimensions[2]/2 + fTankSupportOuterFootDimensions[2]/2 - 0.1*CLHEP::mm));

    G4BooleanSolid* const tankSupport_solid_1 = new G4UnionSolid("tank_support_1", tankSupport_solid_0, outerFoot_solid, 0, trans4);

    const G4ThreeVector trans5(fTankSupportOuterFootDistanceToCenter, 0, -(fTankSupportTopSlabDimensions[2]/2 + fTankSupportOuterFootDimensions[2]/2 - 0.1*CLHEP::mm));

    G4BooleanSolid* const tankSupport_solid = new G4UnionSolid("tank_support", tankSupport_solid_1, outerFoot_solid, 0, trans5);

    // tank support is made of concrete
    G4LogicalVolume* const tankSupport_log = new G4LogicalVolume(tankSupport_solid, fG4_Concrete, "tank_support_log", 0, 0, 0);

    new G4PVPlacement(nullptr, tankSupport_pos, "tank_support", tankSupport_log, expHall_phys, false, 0);

    G4VisAttributes* const tankSupportVisAtt = new G4VisAttributes(G4Colour::Grey());
    tankSupportVisAtt->SetForceSolid(true);
    tankSupportVisAtt->SetVisibility(true);
    tankSupportVisAtt->SetForceWireframe(true);
    if (tankSupport_log)
      tankSupport_log->SetVisAttributes(tankSupportVisAtt);
  }


  void
  G4StationConstruction::GetDataForThisTank()
  {
    SetDetectorParameters();  // updates G4StationConstruction data tables

    // clean up old geometry list
    G4GeometryManager::GetInstance()->OpenGeometry();
    G4PhysicalVolumeStore::GetInstance()->Clean();
    G4LogicalVolumeStore::GetInstance()->Clean();
    G4SolidStore::GetInstance()->Clean();

    // update material properties (deleting out of date materials and
    // replacing them with new ones).
    CreateWater();
    CreatePyrex();
    CreateLucite();
    CreateInterface();
    CreateHDPE();
    CreateLiner();

    // assemble objects with new sizes and materials
    G4RunManager::GetRunManager()->DefineWorldVolume(CreateTank(), true);
  }

}