FDsimG4Materials.cc

Go to the documentation of this file.

#include "FDsimG4Materials.hh"
#include "FDsimG4DetectorConstruction.hh"
#include "G4Material.hh"
#include "G4MaterialTable.hh"
#include "G4Element.hh"
#include "G4ElementTable.hh"
#include "G4UserLimits.hh"
#include "G4RunManager.hh"

#include "geomdefs.hh"
#include <unconfig.h>
#include <utl/AugerUnits.h>
#include <utl/Point.h>
#include <utl/ReferenceEllipsoid.h>
#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Corrector.h>
#include <fdet/Filter.h>
#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>

using namespace TelescopeSimulatorLX ;
using namespace std;

FDsimG4Materials* FDsimG4Materials::fMaterials = 0;

FDsimG4Materials::FDsimG4Materials(const fdet::Telescope& tel)
{
  if(fMaterials)
    { G4Exception("FDsimG4Materials constructed twice."); }

  fMaterials = this;
  fTelescope = &tel;
  ConstructTableOfMaterials() ;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

FDsimG4Materials::~FDsimG4Materials(){;}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void FDsimG4Materials::ConstructTableOfMaterials()
{
  const G4RunManager* runManager = G4RunManager::GetRunManager();
  const FDsimG4DetectorConstruction* DetectorConstruction =  
    dynamic_cast<const FDsimG4DetectorConstruction*>(runManager->GetUserDetectorConstruction());
  G4int VerbosityLevel = DetectorConstruction->fVerbosityLevel;

  det::Detector& detector = det::Detector::GetInstance();
  const fdet::FDetector& detFD = detector.GetFDetector();
  const unsigned int eyeId = fTelescope->GetEyeId();
  const fdet::Eye& eye = detFD.GetEye(eyeId);
  const fdet::Corrector& corrector = fTelescope->GetCorrector();
  const fdet::Filter& filter = fTelescope->GetFilter();

  const utl::Point& eyePos = eye.GetPosition(); 
  const utl::ReferenceEllipsoid& wgs84 = 
    utl::ReferenceEllipsoid::Get(utl::ReferenceEllipsoid::eWGS84);

  const double minWl = detFD.GetModelMinWavelength();
  const double maxWl = detFD.GetModelMaxWavelength();

  double wvl;
  G4double energy;
  G4double rindex;
  G4double transmittance;
  G4double abslength;

  G4double a;
  G4double z; 
  G4double density;
  G4String name;
  G4String symbol;
  G4int nel;

  //    ------------- Elements -------------
  a = 1.01*g/mole;
  G4Element* elH = new G4Element(name="Hydrogen", symbol="H", z=1., a);

  a = 12.01*g/mole;
  G4Element* elC = new G4Element(name="Carbon",   symbol="C", z=6., a);

  a = 14.01*g/mole;
  G4Element* elN = new G4Element(name="Nitrogen", symbol="N", z=7., a);

  a = 16.00*g/mole;
  G4Element* elO = new G4Element(name="Oxygen",   symbol="O", z=8., a);

  //    ------------- Materials -------------
  // Air
  density = 1.29e-03*g/cm3;
  G4Material* Air = new G4Material(name="Air", density, nel=2);
  Air->AddElement(elN, .7);
  Air->AddElement(elO, .3);


  // Atmosphere 
  {
    const atm::Atmosphere& atmosphere = detector.GetAtmosphere();
    double eyeAltitude = wgs84.PointToLatitudeLongitudeHeight (eyePos).get<2> ();

    const atm::ProfileResult& temperatureVsHeight = 
      atmosphere.EvaluateTemperatureVsHeight();
    const atm::ProfileResult& pressureVsHeight = 
      atmosphere.EvaluatePressureVsHeight();
    const atm::ProfileResult& densityVsHeight = 
      atmosphere.EvaluateDensityVsHeight();
    
    G4double pressure = 
      pressureVsHeight.Y(eyeAltitude)/utl::atmosphere * CLHEP::atmosphere;
    G4double temperature = 
      temperatureVsHeight.Y(eyeAltitude)/utl::kelvin * CLHEP::kelvin;
    G4double density = 
      densityVsHeight.Y(eyeAltitude)/(utl::g/utl::cm3) * (CLHEP::g/CLHEP::cm3);


    G4State state = kStateGas;
    const G4int nElements = 2;
    G4Material* Atmosphere = 
      new G4Material("Atmosphere", density,nElements,state,temperature, pressure);
    Atmosphere->AddElement(elN,0.7);
    Atmosphere->AddElement(elO,0.3);
    
  // air refractive index

    G4MaterialPropertiesTable* atmosphereMPT = new G4MaterialPropertiesTable(); 
    atmosphereMPT->AddProperty("RINDEX", new G4MaterialPropertyVector()) ;

    const unsigned int NSTEPS = 10;
    double deltaWl = (maxWl - minWl)/float(NSTEPS-1);   

    for(unsigned int j=0 ; j<NSTEPS ; ++j){
      wvl = minWl + j * deltaWl;
      const atm::ProfileResult& refractionIndexVsHeight = 
        atmosphere.EvaluateRefractionIndexVsHeight(wvl);
      rindex = refractionIndexVsHeight.Y(eyeAltitude);
      energy = (CLHEP::h_Planck*CLHEP::c_light)/(wvl/utl::nanometer * CLHEP::nm);
      atmosphereMPT->AddEntry("RINDEX", energy, rindex);
    }
    Atmosphere->SetMaterialPropertiesTable(atmosphereMPT);

    if (VerbosityLevel > 0) {
      G4cerr << "Material Properties Table for atmosphere: " << G4endl ;
      atmosphereMPT->DumpTable();
    }
  }

  // Aluminum 
  // ---------
  a = 26.98*g/mole;
  density = 2.7*g/cm3;
  new G4Material(name="Aluminum", z=13., a, density); 


  // Plastic  ; 
  // -------------

  density = 1.032*g/cm3;
  G4Material *Plastic = new G4Material(name="MirrorPlastic",density, nel=2);
  Plastic->AddElement(elC,9) ;
  Plastic->AddElement(elH,10) ;

  {
    const G4int    NENTRIES = 2 ;
    G4double LAMBDA[NENTRIES]     = {0.0*eV,10.0*eV};
    G4double REFLECTIVITY[NENTRIES] = {0.91,0.91};

    G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
    MPT->AddProperty("REFLECTIVITY", LAMBDA, REFLECTIVITY, NENTRIES);
    Plastic->SetMaterialPropertiesTable(MPT);
  }


  // BK7 
  // -------------

  // The chemical composition is that of Acrylic. It is not relevant for the optical processes.

  density = 1.19*g/cm3;
  G4Material* BK7 = new G4Material(name="BK7", density, nel=3);
  BK7->AddElement(elC, 5);
  BK7->AddElement(elH, 8);
  BK7->AddElement(elO, 2);

  G4MaterialPropertiesTable *MPT_BK7 = new G4MaterialPropertiesTable();
  MPT_BK7->AddProperty("RINDEX", new G4MaterialPropertyVector()) ;
  MPT_BK7->AddProperty("ABSLENGTH", new G4MaterialPropertyVector()) ;

  const utl::TabulatedFunction& refractiveIndexBK7 = corrector.GetRefractiveIndex();
  const utl::TabulatedFunction& transmittanceBK7 = corrector.GetTransmittance();

  unsigned int nPoints;

  nPoints = refractiveIndexBK7.GetNPoints();
  for(unsigned int j=0 ; j<nPoints ; ++j){
    wvl = refractiveIndexBK7.GetX(j);
    rindex = refractiveIndexBK7.GetY(j);
    energy = (CLHEP::h_Planck*CLHEP::c_light)/(wvl/utl::nanometer * CLHEP::nm);

    MPT_BK7->AddEntry("RINDEX",energy,rindex);
  }

  nPoints = transmittanceBK7.GetNPoints();
  for(unsigned int j=0 ; j<nPoints ; ++j){
    wvl = transmittanceBK7.GetX(j);
    transmittance = transmittanceBK7.GetY(j); // Internal transmittance for 1cm thickness BK7 glass.
    energy = (CLHEP::h_Planck*CLHEP::c_light)/(wvl/utl::nanometer * CLHEP::nm);

    if (transmittance <= DBL_MIN) 
      abslength = DBL_MIN;
    else 
      abslength = -10.0*CLHEP::mm/log(transmittance) ;

    MPT_BK7->AddEntry("ABSLENGTH", energy, abslength);
  }
  
  BK7->SetMaterialPropertiesTable(MPT_BK7);
  
  if (VerbosityLevel > 0) {
    G4cerr << "Material properties table for lens (BK7):" << G4endl ;
    MPT_BK7->DumpTable();
  }
  
  // M-UG6 Filter glass from Schott-Desag 
  // ------------------

  // For now the chemical composition is that of Acrylic - not relevant for the optical processes anyway

  density = 1.19*g/cm3;
  G4Material* MUG6 = new G4Material(name="M-UG6", density, nel=3);
  MUG6->AddElement(elC, 5);
  MUG6->AddElement(elH, 8);
  MUG6->AddElement(elO, 2);

  G4MaterialPropertiesTable *MPT_MUG6 = new G4MaterialPropertiesTable();
  
  MPT_MUG6->AddProperty("RINDEX", new G4MaterialPropertyVector()) ;
  MPT_MUG6->AddProperty("ABSLENGTH", new G4MaterialPropertyVector()) ;

  MPT_MUG6->AddEntry("RINDEX",0.1*eV, 1.526);
  MPT_MUG6->AddEntry("RINDEX",10.0*eV, 1.526);

  // Transmittance for 3.25 mm thickness M-UG6 glass.
  const utl::TabulatedFunction& transmittanceMUG6 = filter.GetTransmittance();

  nPoints = transmittanceMUG6.GetNPoints();
  G4double intTransmittance ;
  for(unsigned int j=0 ; j<nPoints ; ++j){
    wvl = transmittanceMUG6.GetX(j);
    transmittance = transmittanceMUG6.GetY(j);
    energy = (CLHEP::h_Planck*CLHEP::c_light)/(wvl/utl::nanometer * CLHEP::nm);

    // Compute the internal transmittance (Ti) using the relation Ti = T/P,  where T is the transmittance and P is the reflection factor: P = 2n/(n^2+1) 

    rindex = MPT_MUG6->GetProperty("RINDEX")->GetProperty(energy) ;
    G4double P = 2.0*rindex/(rindex*rindex + 1.0) ;

    if (P > 0.0)
      intTransmittance = transmittance/P ;
    else
      intTransmittance =  DBL_MIN;
    
    if (intTransmittance <= DBL_MIN) 
      abslength = DBL_MIN;
    else if (intTransmittance >= 1.0) 
      abslength = DBL_MAX;
    else 
      abslength = -3.25*CLHEP::mm/(std::log(intTransmittance)) ;

    MPT_MUG6->AddEntry("ABSLENGTH", energy, abslength);
  }

  MUG6->SetMaterialPropertiesTable(MPT_MUG6);

  if (VerbosityLevel > 0) {
    G4cerr << "Material properties table for filter (MUG6):" << G4endl ;
    MPT_MUG6->DumpTable();
  }

  // Lime glass window of XP3062 photomultipliers of the FD. 
  // ------------------

  // For now the chemical composition is that of Acrylic - not relevant for the optical processes anyway !!!

  density = 1.19*g/cm3;
  G4Material* LimeGlass = new G4Material(name="LimeGlass", density, nel=3);
  LimeGlass->AddElement(elC, 5);
  LimeGlass->AddElement(elH, 8);
  LimeGlass->AddElement(elO, 2);

  G4MaterialPropertiesTable *MPT_LimeGlass = new G4MaterialPropertiesTable();

  // Data from Schott B270

  MPT_LimeGlass->AddProperty("RINDEX", new G4MaterialPropertyVector()) ;
  MPT_LimeGlass->AddProperty("ABSLENGTH", new G4MaterialPropertyVector()) ;
  

  // n(lambda) for lambda=240,260,280 nm, were obtained by fitting the Cauchy dispersion law 
  // (n(lambda) = a + b/lambda^2) to the data for lambda>= 313nm.
  {
    const G4int nEntries = 10;
    G4double lambda[nEntries] = 
      {240.0*nm,260.0*nm,280.0*nm,313.0*nm,334.0*nm,365.0*nm,405.0*nm,436.0*nm,546.0*nm,700.0*nm};
    G4double Rindex[nEntries] = 
      {1.59,1.58,1.57,1.56,1.5525,1.545,1.54,1.535,1.525,1.52};
    
    G4double energy;
    for (G4int i=0; i < nEntries; i++) {
      energy    = h_Planck*c_light/lambda[i] ;
      MPT_LimeGlass->AddEntry("RINDEX", energy, Rindex[i]);
    }
  }

  {
    const G4int nEntries = 17;
    G4double lambda[nEntries] = {
      240.0*nm,280*nm,300*nm,310*nm,320*nm,330*nm,340*nm,350*nm,360*nm,370*nm,380*nm,
      390*nm,400*nm,450*nm,500*nm,600*nm,700*nm};
    
    G4double Transmittance[nEntries] = 
      {0.0,0.0,35.1,60.0,76.0,84.2,88.0,89.8,90.6,90.8,90.9,
       89.2,90.0,90.3,91.2,91.3,91.1};
    
    G4double Thickness[nEntries] = 
      {1.0*mm, 1.0*mm, 1.0*mm, 1.0*mm, 1.0*mm, 1.0*mm, 1.0*mm, 1.0*mm, 1.0*mm, 1.0*mm, 1.0*mm, 
       10.0*mm, 10.0*mm, 10.0*mm, 10.0*mm, 10.0*mm, 10.0*mm};
    
    G4double abslength ;
    G4double energy    ;
    for (G4int i = 0; i < nEntries; i++) {
      energy    = h_Planck*c_light/lambda[i] ;
      Transmittance[i] = Transmittance[i]/100.0;
      
      // Compute the internal transmittance (Ti) using the relation Ti =
      // T/P,    where T is the transmittance 
      // and P is the reflection factor: P = 2n/(n^2+1) 
     
      G4double Rind = 0.0;
      Rind = MPT_LimeGlass->GetProperty("RINDEX")->GetProperty(energy) ;
      G4double P = 2.0*Rind/(Rind*Rind+1.0) ;
      
      G4double InternalTransmittance ;

      if (P > 0.0)
        InternalTransmittance = Transmittance[i]/P;
      else
        InternalTransmittance = DBL_MIN;
      
      if (InternalTransmittance <= DBL_MIN) 
        abslength = DBL_MIN;
      else if (InternalTransmittance >= 1.0) 
        abslength = DBL_MAX;
      else
        abslength = -Thickness[i]/(std::log(InternalTransmittance));

      MPT_LimeGlass->AddEntry("ABSLENGTH", energy, abslength);
    }
  }
  
  LimeGlass->SetMaterialPropertiesTable(MPT_LimeGlass);
  
  if (VerbosityLevel > 0) {
    G4cerr << "Material properties table for PMT window (Lime Glass):" << G4endl ;
    MPT_LimeGlass->DumpTable();
  }
}