Deprecated/UpgradeASCIITests/G4TankSimulatorASCII/G4TankConstruction.cc

Go to the documentation of this file.

#include "G4TankConstruction.h"

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

#include <sdet/SDetector.h>

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

#include <vector>
//#include <iostream>
#include <sstream>

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

using namespace fwk;

using std::vector;
//using std::cout;
//using std::cerr;
//using std::endl;

using std::ostringstream;

using namespace G4TankSimulatorASCII;

// declare variables that static members above refer to.
double G4TankConstruction::fExpHall_x, G4TankConstruction::fExpHall_y;
double G4TankConstruction::fExpHall_z;
  
double G4TankConstruction::fTankRadius;
double G4TankConstruction::fTankHalfHeight;
double G4TankConstruction::fTankThickness;
double G4TankConstruction::fTankPos_x;
double G4TankConstruction::fTankPos_y;
double G4TankConstruction::fTankPos_z;

utl::Point G4TankConstruction::fPmt1;
utl::Point G4TankConstruction::fPmt2;
utl::Point G4TankConstruction::fPmt3;
  
double G4TankConstruction::fFaceRadius, G4TankConstruction::fFaceRadiusz;
double G4TankConstruction::fFaceActiveRadius, G4TankConstruction::fHeightz;
double G4TankConstruction::fPmtRmin, G4TankConstruction::fPmtRmax;
double G4TankConstruction::fPmtRzmin, G4TankConstruction::fPmtRzmax;
double G4TankConstruction::fInterfaceRmin, G4TankConstruction::fInterfaceRzmin;
double G4TankConstruction::fInterfaceRmax, G4TankConstruction::fInterfaceRzmax;
double G4TankConstruction::fInterfaceThickness, G4TankConstruction::fGlassThickness;
double G4TankConstruction::fDomeThickness;
double G4TankConstruction::fDomeRmin, G4TankConstruction::fDomeRmax;
double G4TankConstruction::fDomeRzmin, G4TankConstruction::fDomeRzmax;
double G4TankConstruction::fMinPhi, G4TankConstruction::fMaxPhi;
double G4TankConstruction::fMinTheta, G4TankConstruction::fMaxTheta;
double G4TankConstruction::fSolarPanelLength, G4TankConstruction::fSolarPanelWidth;
double G4TankConstruction::fSolarPanelThickness;
double G4TankConstruction::fSolarPanelX, G4TankConstruction::fSolarPanelY;
double G4TankConstruction::fSolarPanelZ, G4TankConstruction::fSolarPanelTiltAngle;
double G4TankConstruction::fElecBoxLength, G4TankConstruction::fElecBoxWidth;
double G4TankConstruction::fElecBoxThickness;
double G4TankConstruction::fElecBoxX, G4TankConstruction::fElecBoxY;
double G4TankConstruction::fElecBoxZ, G4TankConstruction::fElecBoxTiltAngle;

TabulatedFunction G4TankConstruction::fPmtfaceRINDEX;
TabulatedFunction G4TankConstruction::fInterfaceRINDEX;
TabulatedFunction G4TankConstruction::fPmtdomeRINDEX;
TabulatedFunction G4TankConstruction::fPmtfaceABSORPTION;
TabulatedFunction G4TankConstruction::fInterfaceABSORPTION;
TabulatedFunction G4TankConstruction::fPmtdomeABSORPTION;
double G4TankConstruction::fSIGMA_ALPHA;  
TabulatedFunction G4TankConstruction::fLinerREFLECTIVITY;
TabulatedFunction G4TankConstruction::fLinerSPECULARLOBECONSTANT;
TabulatedFunction G4TankConstruction::fLinerSPECULARSPIKECONSTANT;
TabulatedFunction G4TankConstruction::fLinerBACKSCATTERCONSTANT;
TabulatedFunction G4TankConstruction::fLinerTYVEK_RINDEX;
TabulatedFunction G4TankConstruction::fLinerABSORPTION;
TabulatedFunction G4TankConstruction::fWaterABSORPTION;
TabulatedFunction G4TankConstruction::fWaterRINDEX;

// Older versions of Geant4 do not have the G4Ellipsoid defined

G4TankConstruction::G4TankConstruction() : expHall_phys(0)
{    
  SetRequiredParameters();
  SetXMLParameters();
  SetDetectorParameters();
  
  // disable solar panel, electronics box, and ground for now:
  
  fSolarPanelEnable = false; 
  fElecBoxEnable    = false;
  fGroundEnable     = false;  
  fASCIIEnable      = true;

  Water = NULL;
  Pyrex = NULL;
  Pyrex1 = NULL;
  Lucite = NULL;
  Interface = NULL;
  HDPE = NULL;
  linerOpticalMPT = NULL;
  Polystyrene = NULL;
  theIron=NULL;
  Plastic=NULL;

  return;
}

G4TankConstruction::~G4TankConstruction() {}

void G4TankConstruction::SetRequiredParameters(void) 
{         
  // The "experimental hall" in which the detector 
  // is contained:
  
  fExpHall_x = 5.0*m;
  fExpHall_y = 5.0*m;
  fExpHall_z = 5.0*m;
    
  // The position of the center of the tank 
  // in the local Geant4 coordinate system
  
  fTankPos_x = 0.0*m;
  fTankPos_y = 0.0*m;
  fTankPos_z = 0.0*m;
}
    
void G4TankConstruction::SetDetectorParameters(void) 
{
  vector<double>::iterator tfIt;
  G4double alpha;
  // Get the pointer to the current detector
  G4TankSimulator* const theTankSimulator = 
    dynamic_cast<G4TankSimulator*>(&RunController::GetInstance().GetModule("G4TankSimulatorASCII"));
  
  const sdet::Station* theCurrentDetectorStation = 
    theTankSimulator->GetCurrentDetectorStation();
  if (!theCurrentDetectorStation) return; // no point getting data yet
  
  // Set tank properties
  
  fTankRadius = 
    theCurrentDetectorStation->GetRadius()*m;
  fTankHalfHeight = 
    theCurrentDetectorStation->GetHeight()*0.5*m;
  fTankThickness = 
    theCurrentDetectorStation->GetThickness()*m;
  
  // Set pmt position (relative to the center of the tank)
  
  fPmt1 = theCurrentDetectorStation->GetPMT(1).GetPosition(); // converted to G4 units later
  fPmt2 = theCurrentDetectorStation->GetPMT(2).GetPosition(); // converted to G4 units later
  fPmt3 = theCurrentDetectorStation->GetPMT(3).GetPosition(); // converted to G4 units later
  
  // Set pmt optical properties
  // For now, assume that all three pmts have the same properties
  
  fPmtfaceRINDEX = 
    theCurrentDetectorStation->GetPMT(1).GetFaceRefractionIndex();
  for (tfIt = fPmtfaceRINDEX.XBegin(); tfIt != fPmtfaceRINDEX.XEnd(); tfIt++)
    (*tfIt) *= eV;
  fInterfaceRINDEX = 
    theCurrentDetectorStation->GetPMT(1).GetRTVRefractionIndex();
  for (tfIt = fInterfaceRINDEX.XBegin(); tfIt != fInterfaceRINDEX.XEnd();
       tfIt++) (*tfIt) *= eV;
  fPmtdomeRINDEX = 
    theCurrentDetectorStation->GetPMT(1).GetDomeRefractionIndex();
  for (tfIt = fPmtdomeRINDEX.XBegin(); tfIt != fPmtdomeRINDEX.XEnd(); tfIt++)
    (*tfIt) *= eV;
  
  // These parameters specify the pmt geometry ( already converted to G4 units )

  fPmtRmin =        fFaceRadius-fGlassThickness;
  fPmtRmax =        fFaceRadius;
  fPmtRzmin =       fFaceRadiusz-fGlassThickness;
  fPmtRzmax =       fFaceRadiusz;
  fInterfaceRmin =  fFaceRadius;
  fInterfaceRzmin = fFaceRadiusz;
  fInterfaceRmax =  fFaceRadius + fInterfaceThickness;
  fInterfaceRzmax = fFaceRadiusz + fInterfaceThickness;
  fDomeRmin =       fFaceRadius + fInterfaceThickness;
  fDomeRzmin =      fFaceRadiusz + fInterfaceThickness;
  fDomeRmax =       fFaceRadius + fInterfaceThickness + fDomeThickness;
  fDomeRzmax =      fFaceRadiusz + fInterfaceThickness + fDomeThickness;

  // Taking into account that the radius of the photocathode is fFaceActiveRadius
  if ( fFaceActiveRadius/mm >= 114. ) 
    {
      fFaceActiveRadius = 113.926112814*mm;
      fHeightz = 0.;
    }
  else if ( fFaceActiveRadius/mm > 97.27 ) 
    {
      alpha = acos(1./62.*(fFaceActiveRadius/mm-133.*sin(47.*deg)+62.*cos(43.*deg)));
      fHeightz = 62.*sin(alpha)*mm;
    }
  else
    {
      alpha = asin(fFaceActiveRadius/mm/133.);
      fHeightz = (133.*cos(alpha)-(133.-62.)*cos(47.*deg))*mm;
    }
  
  // Set liner properties
  
  fSIGMA_ALPHA = theCurrentDetectorStation->GetLinerSigmaAlpha();  
  
  fLinerREFLECTIVITY = theCurrentDetectorStation->GetLinerReflectivity();
  for (tfIt = fLinerREFLECTIVITY.XBegin(); tfIt != fLinerREFLECTIVITY.XEnd();
       tfIt++) (*tfIt) *= eV;
  fLinerSPECULARLOBECONSTANT = theCurrentDetectorStation->GetLinerSpecularLobe();   
  for (tfIt = fLinerSPECULARLOBECONSTANT.XBegin();
       tfIt != fLinerSPECULARLOBECONSTANT.XEnd(); tfIt++) (*tfIt) *= eV;
  fLinerSPECULARSPIKECONSTANT = theCurrentDetectorStation->GetLinerSpecularSpike(); 
  for (tfIt = fLinerSPECULARSPIKECONSTANT.XBegin();
       tfIt != fLinerSPECULARSPIKECONSTANT.XEnd(); tfIt++) (*tfIt) *= eV;
  
  // Set water properties
  
  fWaterABSORPTION = theCurrentDetectorStation->GetWaterAbsorptionLength(); 
  for (tfIt = fWaterABSORPTION.XBegin(); tfIt != fWaterABSORPTION.XEnd();
       tfIt++) (*tfIt) *= eV;
  for (tfIt = fWaterABSORPTION.YBegin(); tfIt != fWaterABSORPTION.YEnd();
       tfIt++) (*tfIt) *= m;
  fWaterRINDEX = theCurrentDetectorStation->GetWaterRefractionIndex();
  for (tfIt = fWaterRINDEX.XBegin(); tfIt != fWaterRINDEX.XEnd(); tfIt++)
    (*tfIt) *= eV;
  return;
}

void G4TankConstruction::SetXMLParameters() 
{
  CentralConfig* theConfig = CentralConfig::GetInstance();
  
  Branch TopBranch = theConfig->GetTopBranch("G4TankSimulatorASCII");
  

  Branch ASCIIEnableB = TopBranch.GetChild("ASCIIEnable");
  ASCIIEnableB.GetData(fASCIIEnable);

  //cout<<"Getting parameters from XML file..."<<endl;
  
  unsigned int size;

  vector<double> pmt_pp;
  vector<double> pmt_face_abs;
  vector<double> pmt_interface_abs;
  vector<double> pmt_dome_abs;
    
  TopBranch.GetChild("pmt_momentum_bins").GetData(pmt_pp);
  TopBranch.GetChild("pmtface_AbsLength").GetData(pmt_face_abs); 
  TopBranch.GetChild("interface_AbsLength").GetData(pmt_interface_abs);
  TopBranch.GetChild("dome_AbsLength").GetData(pmt_dome_abs);
  TopBranch.GetChild("interface_thickness").GetData(fInterfaceThickness);
  TopBranch.GetChild("pmtface_thickness").GetData(fGlassThickness); 
  fInterfaceThickness*=m;
  fGlassThickness*=m;
  
  size = pmt_pp.size();
    
  if ( pmt_face_abs.size()      != size || 
       pmt_interface_abs.size() != size || 
       pmt_dome_abs.size()      != size )
    {
      ostringstream errMsg;
      errMsg<<"TabulatedFunction size mismatch for pmt properties";     
      ERROR(errMsg);

      throw utl::OutOfBoundException(errMsg.str());
    }
  
  for (unsigned int i = 0; i < size; ++i)
    {
      fPmtfaceABSORPTION.PushBack( pmt_pp[i]*eV, pmt_face_abs[i]*m );
      fInterfaceABSORPTION.PushBack( pmt_pp[i]*eV, pmt_interface_abs[i]*m ); 
      fPmtdomeABSORPTION.PushBack( pmt_pp[i]*eV, pmt_dome_abs[i]*m );
    }

  TopBranch.GetChild("face_radius").GetData(fFaceRadius); 
  TopBranch.GetChild("face_radius_z").GetData(fFaceRadiusz); 
  TopBranch.GetChild("face_active_radius").GetData(fFaceActiveRadius); 
  fFaceRadius*=m;
  fFaceRadiusz*=m;
  fFaceActiveRadius*=m;
    
  vector<double> liner_pp;
  vector<double> liner_rindex;
  vector<double> liner_abs;
  vector<double> liner_backscatter;
    
  TopBranch.GetChild("liner_momentum_bins").GetData(liner_pp);
  TopBranch.GetChild("liner_Rindex").GetData(liner_rindex);
  TopBranch.GetChild("liner_AbsLength").GetData(liner_abs); 
  TopBranch.GetChild("backScatter").GetData(liner_backscatter);
  TopBranch.GetChild("dome_thickness").GetData(fDomeThickness);
  fDomeThickness *= m; 

  size = liner_pp.size();
    
  if ( liner_rindex.size() != size ||
       liner_abs.size()    != size ||
       liner_backscatter.size() != size )
    {
      ostringstream errMsg;
      errMsg<<"Tabulated function size mismatch for liner properties";  
      ERROR(errMsg);
        
      throw utl::OutOfBoundException(errMsg.str());
    }

  for (unsigned int i = 0; i < size; ++i)
    {
      fLinerTYVEK_RINDEX.PushBack( liner_pp[i]*eV, liner_rindex[i] );
      fLinerABSORPTION.PushBack( liner_pp[i]*eV, liner_abs[i]*m );
      fLinerBACKSCATTERCONSTANT.PushBack( liner_pp[i]*eV, liner_backscatter[i] );
    }
  return;
}

G4VPhysicalVolume *G4TankConstruction::Construct(void) 
{ 
  if(!expHall_phys){
    CreateElements();
    CreateMaterials();
    return CreateTank();
  }else{
    return expHall_phys;
  }
}
void G4TankConstruction::CreateElements(void) 
{
  G4double a;
  G4String name, symbol;
  G4int z;
  
  // Elements
  //---------

  a = 14.01*g/mole;
  elN = new  G4Element(name = "Nitrogen", symbol = "N", z = 7, a);
  a = 16.00*g/mole;
  elO = new  G4Element(name = "Oxygen", symbol = "O", z = 8, a);
  a = 29.98*g/mole;
  elAl = new G4Element(name = "Aluminum", symbol = "Al", z = 13, a);     
  a = 55.85*g/mole;
  elFe = new G4Element(name = "Iron", symbol = "Fe", z = 26, a);
  a = 1.01*g/mole;
  elH = new  G4Element(name = "Hydrogen", symbol = "H", z = 1, a);
  a = 28.09*g/mole;
  elSi = new G4Element(name = "Silicon", symbol = "Si", z = 14, a);
  a = 10.811*g/mole;
  elB = new  G4Element(name = "Boron", symbol = "B", z = 5, a);
  a = 22.98977*g/mole;
  elNa = new G4Element(name = "Sodium", symbol = "Na", z = 11, a);
  a = 12.0107*g/mole;
  elC = new  G4Element(name = "Carbon", symbol = "C", z = 6, a);
}
  
void G4TankConstruction::CreateMaterials(void) 
{
  CreateAir();
  CreateWater();
  CreateVacuum();
  CreatePyrex();
  CreatePyrex1();
  CreateLucite();
  CreateInterface();
  CreateHDPE();
  CreateLiner();
  CreateAluminium();
  CreateDirt();

  CreatePolystyrene();
  CreateIron();
  CreatePlastic();

  return;
}

G4VPhysicalVolume *G4TankConstruction::CreateTank(void)
{
  CreatePrimitives();
  CreateHall();
  AssembleTank();
  return expHall_phys;
}

void G4TankConstruction::CreateAir(void)
{
  // Air
  //----
  G4String name;
  G4double density = 1.29e-03*g/cm3;
  G4int nel;
  Air = new G4Material(name = "Air_G4Tank", density, nel = 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.0*m, 10000.0*m };
  
  airMPT = new G4MaterialPropertiesTable();
  airMPT->AddProperty("RINDEX", airPP, airRINDEX, 2);
  airMPT->AddProperty("ABSLENGTH", airPP, airABSLENGTH, 2);

  Air->SetMaterialPropertiesTable(airMPT);
  return;
}

void G4TankConstruction::CreateWater(void)
{

  // Water
  //------
  G4String name;
  G4double density = 1.0*g/cm3;
  G4int nel;

  if (Water != NULL) {
    delete Water;
    delete waterMPT;
  }
  Water = new G4Material(name = "Water", density, nel = 2);
  Water->AddElement(elH, 2);
  Water->AddElement(elO, 1);
      
  // Set water material properties

  waterMPT = new G4MaterialPropertiesTable();
  
  // water rindex

  int index;
  int num = fWaterRINDEX.GetNPoints();

  double* array_ppWATER     = new double[num+1];
  double* array_waterRINDEX = new double[num+1];
  
  for ( index = 0; index < num; index++ ) 
    {
      array_ppWATER[index]     = fWaterRINDEX[index].X();       
      array_waterRINDEX[index] = fWaterRINDEX[index].Y();
    }
  
  waterMPT->AddProperty("RINDEX", array_ppWATER, array_waterRINDEX, num);    
  
  // water abs. length
  
  num = fWaterABSORPTION.GetNPoints();
  
  double* array_ppwaterABSORPTION = new double[num+1];
  double* array_waterABSORPTION   = new double[num+1];
  
  for ( index = 0; index < num; index++ ) 
    {
      array_ppwaterABSORPTION[index] = fWaterABSORPTION[index].X();
      array_waterABSORPTION[index]   = fWaterABSORPTION[index].Y();
    }
  
  waterMPT->AddProperty("ABSLENGTH", array_ppwaterABSORPTION, array_waterABSORPTION, num);
  
  //cout<<"Water Material Properties:"<<endl;
  //waterMPT->DumpTable();
  
  Water->SetMaterialPropertiesTable(waterMPT);
  
  delete[] array_ppWATER;
  delete[] array_waterRINDEX;
  delete[] array_ppwaterABSORPTION;
  delete[] array_waterABSORPTION;
  return;
}

void G4TankConstruction::CreateVacuum(void)
{
  // Vacuum
  //-------
  G4String name;
  G4double a = 1.* g / mole; 
  G4double density = universe_mean_density;
  G4int z;

  // ATTENTION: change from STP to avoid the warning...
  Vacuum =  new G4Material(name = "Vacuum", z = 1, a, density,
                                       kStateUndefined, STP_Temperature, 
                                       STP_Pressure);
  return;
}

void G4TankConstruction::CreatePyrex(void)
{  
  // 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?
  
  G4String name;
  G4double density = 2.23*g/cm3;
  G4int nel, natoms, ncomponents;
  
  if (Pyrex != NULL) {
    delete SiO2;
    delete B2O2;
    delete Na2O;
    delete Pyrex;
    delete pmtfaceMPT;
  }
  
  SiO2 = new G4Material(name = "SiO2", density, nel = 2);
  SiO2->AddElement(elSi, natoms = 1);
  SiO2->AddElement(elO, natoms = 2);
  
  B2O2 = new G4Material(name = "B2O2", density, nel = 2);
  B2O2->AddElement(elB, natoms = 2);
  B2O2->AddElement(elO, natoms = 2);
  
  Na2O = new G4Material(name = "Na2O", density, nel = 2);
  Na2O->AddElement(elNa, natoms = 2);
  Na2O->AddElement(elO, natoms = 1);
  
  Pyrex = new G4Material(name = "Pyrex", density, ncomponents = 3);
  Pyrex->AddMaterial(SiO2, 0.80);
  Pyrex->AddMaterial(B2O2, 0.13);
  Pyrex->AddMaterial(Na2O, 0.07);
  
  // Fill material properties table for pmt face
  
  pmtfaceMPT = new G4MaterialPropertiesTable();
  
  // pmt face rindex

  int index;
  int num = fPmtfaceRINDEX.GetNPoints();
    
  double* array_pmtfaceRINDEX_PP = new double[num+1];
  double* array_pmtfaceRINDEX    = new double[num+1];

  // debug
  //fout<<"PMT face rindex:"<<endl;
    
  for ( index = 0; index < num; index++ ) 
    {            
      array_pmtfaceRINDEX_PP[index] = fPmtfaceRINDEX[index].X();
      array_pmtfaceRINDEX[index]    = fPmtfaceRINDEX[index].Y();
    }
  
  pmtfaceMPT->AddProperty("RINDEX", array_pmtfaceRINDEX_PP, array_pmtfaceRINDEX, num);

  // Fill and add pmt face absorption length

  num = fPmtfaceABSORPTION.GetNPoints();
    
  double* array_pmtfaceABSORPTION_PP = new double[num+1];
  double* array_pmtfaceABSORPTION    = new double[num+1];

  // debug 
  //fout<<"PMT face abs. length:"<<endl;
    
  for ( index = 0; index < num; index++ ) 
    {
      // debug
      //fout<< fPmtfaceABSORPTION[index].X() <<"    "<< fPmtfaceABSORPTION[index].Y() <<endl;
        
      array_pmtfaceABSORPTION_PP[index] = fPmtfaceABSORPTION[index].X();
      array_pmtfaceABSORPTION[index]    = fPmtfaceABSORPTION[index].Y();
    }
    
  pmtfaceMPT->AddProperty("ABSLENGTH", array_pmtfaceABSORPTION_PP, array_pmtfaceABSORPTION, num);
  

  //cout<<"PMT Face Material Properties:"<<endl;
  //pmtfaceMPT->DumpTable();
  
  Pyrex->SetMaterialPropertiesTable(pmtfaceMPT);

  delete[] array_pmtfaceRINDEX_PP;
  delete[] array_pmtfaceRINDEX;
  delete[] array_pmtfaceABSORPTION_PP;
  delete[] array_pmtfaceABSORPTION;
  return;
}

void G4TankConstruction::CreatePyrex1(void)
{  
  // 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.
  */

  G4String name;
  G4double density = 2.23 * g / cm3;
  G4int nel, natoms;
  
  if (Pyrex1 != NULL) {
    delete SiO2;
    delete B2O2;
    delete Na2O;
    delete Pyrex1;
    delete pmtfaceMPT1;
  }
  SiO2 = new G4Material(name = "SiO2", density, nel = 2);
  SiO2->AddElement(elSi, natoms = 1);
  SiO2->AddElement(elO, natoms = 2);
  
  B2O2 = new G4Material(name = "B2O2", density, nel = 2);
  B2O2->AddElement(elB, natoms = 2);
  B2O2->AddElement(elO, natoms = 2);
    
  Na2O = new G4Material(name = "Na2O", density, nel = 2);
  Na2O->AddElement(elNa, natoms = 2);
  Na2O->AddElement(elO, natoms = 1);
  
  int ncomponents;
  Pyrex1 = new G4Material(name = "Pyrex1", density, ncomponents = 3);
  Pyrex1->AddMaterial(SiO2, 0.80);
  Pyrex1->AddMaterial(B2O2, 0.13);
  Pyrex1->AddMaterial(Na2O, 0.07);
  
  pmtfaceMPT1 = new G4MaterialPropertiesTable();

  // Fill material properties table for non sensitive pmt face
  
  // non-sensitive pmt face rindex = sensitive pmt face rindex

  int index;
  int num = fPmtfaceRINDEX.GetNPoints();
    
  double* array_pmtface1RINDEX_PP = new double[num+1];
  double* array_pmtface1RINDEX    = new double[num+1];

  for ( index = 0; index < num; index++ ) 
    {            
      array_pmtface1RINDEX_PP[index] = fPmtfaceRINDEX[index].X();
      array_pmtface1RINDEX[index]    = fPmtfaceRINDEX[index].Y();
    }
  
  pmtfaceMPT1->AddProperty("RINDEX", array_pmtface1RINDEX_PP, array_pmtface1RINDEX, num);


  // Fill and add non-sensitive pmt face absorption length = dome absorption length

  num = fPmtdomeABSORPTION.GetNPoints();

  double* array_pmtface1ABSORPTION_PP = new double[num+1];
  double* array_pmtface1ABSORPTION    = new double[num+1];

  for ( index = 0; index < num; index++ ) 
    {
      array_pmtface1ABSORPTION_PP[index] = fPmtdomeABSORPTION[index].X();
      array_pmtface1ABSORPTION[index]    = fPmtdomeABSORPTION[index].Y();
    }
    
  pmtfaceMPT1->AddProperty("ABSLENGTH", array_pmtface1ABSORPTION_PP, array_pmtface1ABSORPTION, num);

  Pyrex1->SetMaterialPropertiesTable(pmtfaceMPT1);
  
  delete[] array_pmtface1RINDEX_PP;
  delete[] array_pmtface1RINDEX;
  delete[] array_pmtface1ABSORPTION_PP;
  delete[] array_pmtface1ABSORPTION;
  return;
}
  
void G4TankConstruction::CreateLucite(void)
{
  // Lucite (for the pmt domes)
  G4double density = 1.20*g/cm3;
  G4String name;
  G4int nel, natoms;

  if (Lucite != NULL) {
    delete CH3;
    delete CH2;
    delete C;
    delete CO2;
    delete Lucite;
    delete pmtdomeMPT;
  }
  
  CH3 = new G4Material(name = "CH3", density, nel = 2);
  CH3->AddElement(elC, natoms = 1);
  CH3->AddElement(elH, natoms = 3);
  
  CH2 = new G4Material(name = "CH2", density, nel = 2);
  CH2->AddElement(elC, natoms = 1);
  CH2->AddElement(elH, natoms = 2);
    
  C = new G4Material(name = "C", density, nel = 1);
  C->AddElement(elC, natoms = 1);
  
  CO2 = new G4Material(name = "CO2", density, nel = 2);
  CO2->AddElement(elC, natoms = 1);
  CO2->AddElement(elO, natoms = 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.
  */
    
  Lucite = new G4Material(name="Lucite", density, nel=3);
  Lucite->AddElement(elC, natoms=4);
  Lucite->AddElement(elH, natoms=8);
  Lucite->AddElement(elO, natoms=2);
  
  pmtdomeMPT = new G4MaterialPropertiesTable();

  // pmt dome rindex
  
  int index;
  int num = fPmtdomeRINDEX.GetNPoints();
  
  double *array_pmtdomeRINDEX_PP = new double[num+1];
  double *array_pmtdomeRINDEX    = new double[num+1];

  for ( index = 0; index < num; index++ ) 
    {
      array_pmtdomeRINDEX_PP[index] = fPmtdomeRINDEX[index].X();
      array_pmtdomeRINDEX[index]    = fPmtdomeRINDEX[index].Y();
    }
    
  pmtdomeMPT->AddProperty("RINDEX", array_pmtdomeRINDEX_PP, array_pmtdomeRINDEX, num);

  // pmt dome absorption length

  num = fPmtdomeABSORPTION.GetNPoints();

  double *array_pmtdomeABSORPTION_PP = new double[num+1];
  double *array_pmtdomeABSORPTION    = new double[num+1];

  for ( index = 0; index < num; index++ ) 
    {
      array_pmtdomeABSORPTION_PP[index] = fPmtdomeABSORPTION[index].X();
      array_pmtdomeABSORPTION[index]    = fPmtdomeABSORPTION[index].Y();
    }
  
  pmtdomeMPT->AddProperty("ABSLENGTH", array_pmtdomeABSORPTION_PP, array_pmtdomeABSORPTION, num);
  
  //cout<<"PMT Dome Material Properties:"<<endl;
  //pmtdomeMPT->DumpTable();
  
  Lucite->SetMaterialPropertiesTable(pmtdomeMPT);

  delete[] array_pmtdomeRINDEX_PP;
  delete[] array_pmtdomeRINDEX;
  delete[] array_pmtdomeABSORPTION_PP;
  delete[] array_pmtdomeABSORPTION;
  return;
}
  
void G4TankConstruction::CreateInterface(void)
{
  // Interface for production tanks is Wacker SilGel 612 (Registered trademark
  // of Wacker Group
  G4String name;
  G4int nel, natoms;
  G4double density = 0.97*g/cm3; // guess - not explicitly set by Troy

  if (Interface != NULL) {
    delete Interface;
    delete interfaceMPT;
  }
  Interface = new G4Material(name = "Interface", density, nel = 3);
  Interface->AddElement(elC, natoms = 4);
  Interface->AddElement(elH, natoms = 8);
  Interface->AddElement(elO, natoms = 2);
  
  interfaceMPT = new G4MaterialPropertiesTable();
  
  // interface rindex 

  int index;
  int num = fInterfaceRINDEX.GetNPoints();
  
  double* array_interfaceRINDEX_PP = new double[num+1];
  double* array_interfaceRINDEX    = new double[num+1];
    
  for ( index = 0; index < num; index++ ) 
    {
      array_interfaceRINDEX_PP[index] = fInterfaceRINDEX[index].X();
      array_interfaceRINDEX[index]    = fInterfaceRINDEX[index].Y();
    }
  
  interfaceMPT->AddProperty("RINDEX", array_interfaceRINDEX_PP, array_interfaceRINDEX, num); 
    
  // interface abs. length

  num = fInterfaceABSORPTION.GetNPoints();
    
  double* array_interfaceABSORPTION_PP = new double[num+1];
  double* array_interfaceABSORPTION    = new double[num+1];

  for ( index = 0; index < num; index++ ) 
    {
      array_interfaceABSORPTION_PP[index] = fInterfaceABSORPTION[index].X();
      array_interfaceABSORPTION[index]    = fInterfaceABSORPTION[index].Y();
    }
        
  interfaceMPT->AddProperty("ABSLENGTH", array_interfaceABSORPTION_PP, array_interfaceABSORPTION, num);
  
  //cout<<"Interface Material Properties:"<<endl;
  //interfaceMPT->DumpTable();
  
  Interface->SetMaterialPropertiesTable(interfaceMPT);    

  delete[] array_interfaceRINDEX_PP;
  delete[] array_interfaceRINDEX;
  delete[] array_interfaceABSORPTION_PP;
  delete[] array_interfaceABSORPTION;
  return;
}
void G4TankConstruction::CreateHDPE(void)
{
  
  // HDPE (for tank walls)
  //----------------------

  G4double density = 0.94*g/cm3;
  G4int nel, natoms;
  G4String name;
  vector<double>::iterator It;

  if (HDPE != NULL) {
    delete HDPE;
    delete linerMPT;
  }

  HDPE = new G4Material(name = "HDPE", density, nel = 2);
  HDPE->AddElement(elC, natoms = 2);
  HDPE->AddElement(elH, natoms = 4);
  
  linerMPT = new G4MaterialPropertiesTable();
  
  // liner abs. length
  
  int index;
  int num = fLinerABSORPTION.GetNPoints();
  
  double* array_linerABSORPTION_PP = new double[num+1];
  double* array_linerABSORPTION    = new double[num+1];

  for ( index = 0; index < num; index++ )
    {
      array_linerABSORPTION_PP[index] = fLinerABSORPTION[index].X();
      array_linerABSORPTION[index]    = fLinerABSORPTION[index].Y();
    }
    
  linerMPT->AddProperty("ABSLENGTH", array_linerABSORPTION_PP, array_linerABSORPTION, num);

  //cout<<"Liner Material Properties:"<<endl;
  //linerMPT->DumpTable();
  
  HDPE->SetMaterialPropertiesTable(linerMPT);
 
  delete[] array_linerABSORPTION_PP;
  delete[] array_linerABSORPTION;
  return;
}
    
void G4TankConstruction::CreateLiner(void)
{
  // 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 
    Set index of refraction on the small dielectric interface to 0 so that 
    all the light will undergo total internal reflection (or absorption).  
    It is done in this weird way because the only surface that works (at the 
    moment) for the unified model is the ground surface, and for such a 
    surface a transmission coefficient is always calculated if total internal
    reflection does not occur.  Thus, the only way to ensure all the 
    (non-absorbed) light is reflected back is to choose groundbackpainted 
    finish with the index of refraction for the thin dielectric material 
    set to 0.  Note that groundbackpainted with some other index of 
    refraction doesn't work because light reflected off the back (painted) 
    surface undergoes only Lambertian reflection.. This will be fixed up 
    in the future..
  */

  if (linerOpticalMPT != NULL) delete linerOpticalMPT;    
  linerOpticalMPT = new G4MaterialPropertiesTable();
  
  int index;
  int num = fLinerSPECULARLOBECONSTANT.GetNPoints();
  double* array_linerSPECULARLOBECONSTANT_PP = new double[num+1];
  double* array_linerSPECULARLOBECONSTANT    = new double[num+1];

  for ( index = 0; index < num; index++ ) 
    {
      array_linerSPECULARLOBECONSTANT_PP[index] = fLinerSPECULARLOBECONSTANT[index].X();
      array_linerSPECULARLOBECONSTANT[index]    = fLinerSPECULARLOBECONSTANT[index].Y();
    }
    
  linerOpticalMPT->AddProperty("SPECULARLOBECONSTANT", 
                               array_linerSPECULARLOBECONSTANT_PP,
                               array_linerSPECULARLOBECONSTANT, num);


  num = fLinerSPECULARSPIKECONSTANT.GetNPoints();
  double* array_linerSPECULARSPIKECONSTANT_PP = new double[num+1];
  double* array_linerSPECULARSPIKECONSTANT    = new double[num+1];

  for ( index = 0; index < num; index++ ) 
    {
      array_linerSPECULARSPIKECONSTANT_PP[index] = fLinerSPECULARSPIKECONSTANT[index].X();
      array_linerSPECULARSPIKECONSTANT[index]    = fLinerSPECULARSPIKECONSTANT[index].Y();
    }

  linerOpticalMPT->AddProperty("SPECULARSPIKECONSTANT", 
                               array_linerSPECULARSPIKECONSTANT_PP,
                               array_linerSPECULARSPIKECONSTANT, num);


  num = fLinerREFLECTIVITY.GetNPoints();
  double* array_linerREFLECTIVITY_PP = new double[num+1];
  double* array_linerREFLECTIVITY    = new double[num+1];

  for ( index = 0; index < num; index++ ) 
    {  
      array_linerREFLECTIVITY_PP[index] = fLinerREFLECTIVITY[index].X();
      array_linerREFLECTIVITY[index]    = fLinerREFLECTIVITY[index].Y();
    }

    
  linerOpticalMPT->AddProperty("REFLECTIVITY",  
                               array_linerREFLECTIVITY_PP,       
                               array_linerREFLECTIVITY, num);
  

  num = fLinerBACKSCATTERCONSTANT.GetNPoints();
  double* array_linerBACKSCATTERCONSTANT_PP = new double[num+1];
  double* array_linerBACKSCATTERCONSTANT    = new double[num+1];

  for ( index = 0; index < num; index++ ) 
    {
      array_linerBACKSCATTERCONSTANT_PP[index] = fLinerBACKSCATTERCONSTANT[index].X();
      array_linerBACKSCATTERCONSTANT[index]    = fLinerBACKSCATTERCONSTANT[index].Y();
    }

  linerOpticalMPT->AddProperty("BACKSCATTERCONSTANT", 
                               array_linerBACKSCATTERCONSTANT_PP, 
                               array_linerBACKSCATTERCONSTANT, num);

  num = fLinerTYVEK_RINDEX.GetNPoints();
  double* array_TYVEK_RINDEX_PP = new double[num+1];
  double* array_TYVEK_RINDEX    = new double[num+1];
     
  for ( index = 0; index < num; index++ ) 
    {  
      array_TYVEK_RINDEX_PP[index] = fLinerTYVEK_RINDEX[index].X();
      array_TYVEK_RINDEX[index]    = fLinerTYVEK_RINDEX[index].Y();
    }

  //cout<<"Liner Optical Properties:"<<endl;
  //linerOpticalMPT->DumpTable();
  linerOpticalMPT->AddProperty("RINDEX",  
                               array_TYVEK_RINDEX_PP,              
                               array_TYVEK_RINDEX, num); 

  OpLinerSurface = new G4OpticalSurface("WallSurface");

  OpLinerSurface->SetModel(unified);
  OpLinerSurface->SetType(dielectric_dielectric);
  OpLinerSurface->SetFinish(groundbackpainted);
  OpLinerSurface->SetSigmaAlpha(fSIGMA_ALPHA);
    
  OpLinerSurface->SetMaterialPropertiesTable(linerOpticalMPT); 
  
  delete[] array_linerSPECULARLOBECONSTANT_PP;
  delete[] array_linerSPECULARLOBECONSTANT;
    
  delete[] array_linerSPECULARSPIKECONSTANT_PP;
  delete[] array_linerSPECULARSPIKECONSTANT;
    
  delete[] array_linerREFLECTIVITY_PP;
  delete[] array_linerREFLECTIVITY;
    
  delete[] array_linerBACKSCATTERCONSTANT_PP;
  delete[] array_linerBACKSCATTERCONSTANT;
    
  delete[] array_TYVEK_RINDEX_PP;
  delete[] array_TYVEK_RINDEX;

  return;

}

   
void G4TankConstruction::CreateAluminium(void)
{
  // Aluminum for electronics boxes and (possibly) solar panel
  //----------------------------------------------------------

  G4double density = 2.700*g/cm3;
  G4double a = 29.98*g/mole;
  G4String name;
  G4int z;
  Al = new G4Material(name = "Aluminum", z = 13, a, density);
  return;
}
  
void G4TankConstruction::CreateDirt(void)
{
  // Dirt for the ground
  //--------------------
  G4String name;
  G4int ncomponents;
  double fractionmass;    
  G4double density = 1.2 * g / cm3;
  Dirt = new G4Material(name = "Dirt", density, ncomponents = 4);
  Dirt->AddElement(elO, fractionmass = 51.0*perCent);
  Dirt->AddElement(elSi, fractionmass = 35.2*perCent);
  Dirt->AddElement(elAl, fractionmass = 8.2*perCent);
  Dirt->AddElement(elFe, fractionmass = 5.6*perCent);
}



void G4TankConstruction::CreatePolystyrene(void)
{
  
  G4String name;
  G4double density = 1.06*g/cm3;
  G4int nel;

  if (Polystyrene != NULL) {
    delete Polystyrene;
  }

  Polystyrene = new G4Material(name = "Polystyrene", density, nel = 2);
  Polystyrene->AddElement(elH, 1);
  Polystyrene->AddElement(elC, 1);
       
}

void G4TankConstruction::CreatePlastic(void) // To be changed to PVC
{
  
  G4String name;
  G4double density = 0.7*g/cm3;
  G4int nel;

  if (Plastic != NULL) {
    delete Plastic;
  }

  Plastic = new G4Material(name = "Plastic", density, nel = 2);
  Plastic->AddElement(elH, 1);
  Plastic->AddElement(elC, 1);
       
}
   
void G4TankConstruction::CreateIron(void) 
{
  if (theIron!=NULL){
    delete theIron;
  }
  G4double density = 7.800*g/cm3;
  G4double a = 55.845*g/mole;
  G4String name;
  G4int z;
  theIron = new G4Material(name = "Iron", z = 26, a, density);
  return;
}

void G4TankConstruction::CreatePrimitives(void)
{
  G4double top_TankThickness = (fTankThickness/2);
  G4double fGroundThickness = fExpHall_z/4; // arbitrary but probably okay

  expHall_box = new G4Box("World", fExpHall_x, fExpHall_y, fExpHall_z);
  ground_solid = new G4Tubs("ground_solid", 0, fExpHall_x, fGroundThickness,          
                   0.0*deg, 360.0*deg);
  tank_solid = new G4Tubs("tank_solid", 0.0*m, fTankRadius, fTankHalfHeight,
               0.0*deg, 360.0*deg);
  top_solid = new G4Tubs("top_solid", 
               0.0*m, fTankRadius + fTankThickness, top_TankThickness,
               0.0*deg, 360.0*deg);
  side_solid = new G4Tubs("side_solid", 
               fTankRadius, fTankRadius + fTankThickness, fTankHalfHeight,
               0.0*deg, 360.0*deg);

  inner_solid = new G4Ellipsoid("inner_solid", fPmtRmin, fPmtRmin, fPmtRzmin, -fPmtRzmin, 0.);

  // Active photocathode covered area
  pmt_aux = new G4Ellipsoid("pmt_aux", fPmtRmax, fPmtRmax, fPmtRzmax, -fPmtRzmax, -fHeightz);

  pmt_solid = new G4SubtractionSolid("pmt_solid", pmt_aux, inner_solid);
  pmt_aux1 = new G4Ellipsoid("pmt_aux1", fPmtRmax, fPmtRmax, fPmtRzmax, -fHeightz, 0.);
  pmt_solid1 = new G4SubtractionSolid("pmt_solid1", pmt_aux1, inner_solid);

  interface_in_aux = new G4Ellipsoid("interface_in_aux", fInterfaceRmin, fInterfaceRmin,        
                      fInterfaceRzmin,-fInterfaceRzmin,0.);
  interface_out_aux = new G4Ellipsoid("interface_out_aux", fInterfaceRmax, fInterfaceRmax,      
                      fInterfaceRzmax, -fInterfaceRzmax, 0.0*cm);
  interface_solid = new G4SubtractionSolid("interface_solid", interface_out_aux, interface_in_aux);
  dome_in_aux = new G4Ellipsoid("dome_in_aux", fDomeRmin, fDomeRmin, fDomeRzmin, -fDomeRzmin, 0.);
  dome_out_aux = new G4Ellipsoid("dome_out_aux", fDomeRmax, fDomeRmax, fDomeRzmax, -fDomeRzmax, 0.0*cm);
  dome_solid = new G4SubtractionSolid("dome_solid", dome_out_aux, dome_in_aux);

  return;
}
  
void G4TankConstruction::CreateHall(void)
{
  // Experimental hall
  //------------------

  expHall_log = new G4LogicalVolume(expHall_box, Air, "World", 0, 0, 0);
  
  expHall_phys = new G4PVPlacement(0, G4ThreeVector(), "World", expHall_log, 0, 
                        false, 0);
  
  expHall_log->SetVisAttributes(G4VisAttributes::Invisible);

  // Ground around the detector
  // --------------------------

  if (fGroundEnable) 
    {   
      double fGroundThickness = fExpHall_z/4; // arbitrary but probably okay
    
      ground_log = new G4LogicalVolume(ground_solid, Dirt, "ground_solid"); 
    
      ground_phys = new G4PVPlacement(0, 
                          G4ThreeVector(0, 0,- fGroundThickness), 
                          "ground", ground_log, expHall_phys, false, 0); 
    }
  return;
}

void G4TankConstruction::AssembleTank(void)
{
  G4SDManager *SDMan = G4SDManager::GetSDMpointer();
  G4String SDName;
  
  // Water tank
  //-----------

  tank_log = new G4LogicalVolume(tank_solid, Water, "tank_log", 0, 0, 0);
  
  tank_phys = new G4PVPlacement(0, G4ThreeVector(fTankPos_x, fTankPos_y, fTankPos_z + 
                                                 fTankHalfHeight + fTankThickness),
                        "tank", tank_log, expHall_phys, false, 0);
  
  // Top, bottom, side walls
  //------------------------

  top_log = new G4LogicalVolume(top_solid, HDPE, "top_log", 0, 0, 0);
  bottom_log = new G4LogicalVolume(top_solid, HDPE, "bottom_log", 0, 0, 0);
  side_log = new G4LogicalVolume(side_solid, HDPE, "side_log", 0, 0, 0);
  
  top_phys = new G4PVPlacement(0, G4ThreeVector(fTankPos_x, fTankPos_y, 
                                                fTankPos_z + 2.*fTankHalfHeight + 
                                                (fTankThickness)/2.0 + fTankThickness), 
                               "top", top_log, expHall_phys, false, 0);
  bottom_phys = new G4PVPlacement(0,G4ThreeVector(fTankPos_x, fTankPos_y, 
                                      fTankPos_z + fTankThickness - (fTankThickness)/2.0),
                                  "bottom", bottom_log, expHall_phys, false, 0);
  side_phys = new G4PVPlacement(0,G4ThreeVector(fTankPos_x, fTankPos_y, fTankPos_z + 
                                                + fTankHalfHeight + fTankThickness),
                                "side", side_log, expHall_phys, false, 0);
  
  //the tank liner  
  topsurface = 
    new G4LogicalBorderSurface("topsurface", tank_phys, top_phys,
                               OpLinerSurface);
  bottomsurface = 
    new G4LogicalBorderSurface("bottomsurface", tank_phys, bottom_phys,
                               OpLinerSurface);
  sidesurface = 
    new G4LogicalBorderSurface("sidesurface", tank_phys, side_phys,
                               OpLinerSurface);

  // PMTs and such
  //--------------
  // the interior volume of the pmt
  
  // 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
  double dome1_x = fPmt1.GetX(fPmt1.GetCoordinateSystem())*m;  
  double dome1_y = fPmt1.GetY(fPmt1.GetCoordinateSystem())*m;
  double dome1_z = fTankHalfHeight;
    
  double dome2_x = fPmt2.GetX(fPmt2.GetCoordinateSystem())*m;
  double dome2_y = fPmt2.GetY(fPmt2.GetCoordinateSystem())*m;
  double dome2_z = fTankHalfHeight;
  
  double dome3_x = fPmt3.GetX(fPmt3.GetCoordinateSystem())*m;
  double dome3_y = fPmt3.GetY(fPmt3.GetCoordinateSystem())*m;
  double dome3_z = fTankHalfHeight;
  
  inner_log = new G4LogicalVolume(inner_solid, Vacuum, "inner_log", 0, 0, 0);
  
  inner1_phys = new G4PVPlacement(0, 
                        G4ThreeVector(dome1_x, dome1_y, dome1_z), inner_log,
                        "inner1", tank_log, false, 0);
  inner2_phys = new G4PVPlacement(0, 
                        G4ThreeVector(dome2_x, dome2_y, dome2_z), inner_log,
                        "inner2", tank_log, false, 0);
  inner3_phys = new G4PVPlacement(0, 
                        G4ThreeVector(dome3_x, dome3_y, dome3_z), inner_log,
                        "inner3", tank_log, false, 0);
  /*
    Why 3 different logical volumes?
    because the sensitive detectors are registered 
    by logical volume
  */
        
  // Sensitive surface of pmt faces

  pmt1_log = new G4LogicalVolume(pmt_solid, Pyrex, "pmt1_log", 0, 0, 0);
  pmt2_log = new G4LogicalVolume(pmt_solid, Pyrex, "pmt2_log", 0, 0, 0);
  pmt3_log = new G4LogicalVolume(pmt_solid, Pyrex, "pmt3_log", 0, 0, 0);
  
  pmt1_phys = new G4PVPlacement(0, G4ThreeVector(dome1_x, dome1_y, dome1_z), 
                        pmt1_log,"pmt1", tank_log, false, 0);
  pmt2_phys = new G4PVPlacement(0, G4ThreeVector(dome2_x, dome2_y, dome2_z), 
                        pmt2_log,"pmt2", tank_log, false, 0);
  pmt3_phys = new G4PVPlacement(0, G4ThreeVector(dome3_x, dome3_y, dome3_z), 
                        pmt3_log,"pmt3", tank_log, false, 0);

  // Non-sensitive surface of pmt faces

  pmt1_log1 = new G4LogicalVolume(pmt_solid1, Pyrex1, "pmt1_log1", 0, 0, 0);
  pmt2_log1 = new G4LogicalVolume(pmt_solid1, Pyrex1, "pmt2_log1", 0, 0, 0);
  pmt3_log1 = new G4LogicalVolume(pmt_solid1, Pyrex1, "pmt3_log1", 0, 0, 0);
  
  pmt1_phys1 = new G4PVPlacement(0, G4ThreeVector(dome1_x, dome1_y, dome1_z), 
                        pmt1_log1,"pmt11", tank_log, false, 0);
  pmt2_phys1 = new G4PVPlacement(0, G4ThreeVector(dome2_x, dome2_y, dome2_z), 
                        pmt2_log1,"pmt21", tank_log, false, 0);
  pmt3_phys1 = new G4PVPlacement(0, G4ThreeVector(dome3_x, dome3_y, dome3_z), 
                        pmt3_log1,"pmt31", tank_log, false, 0);

  // The dome/face interface
    
  interface1_log = new G4LogicalVolume(interface_solid, Interface, "interface1_log", 
                          0, 0, 0);
  interface2_log = new G4LogicalVolume(interface_solid, Interface, "interface2_log", 
                          0, 0, 0);
  interface3_log = new G4LogicalVolume(interface_solid, Interface, "interface3_log", 
                          0, 0, 0);
  
  interface1_phys = new G4PVPlacement(0, G4ThreeVector(dome1_x, dome1_y, dome1_z), 
                        interface1_log,"interface1", tank_log, false, 0);
  interface2_phys = new G4PVPlacement(0, G4ThreeVector(dome2_x, dome2_y, dome2_z), 
                        interface2_log,"interface2", tank_log, false, 0);
  interface3_phys = new G4PVPlacement(0, G4ThreeVector(dome3_x, dome3_y, dome3_z), 
                        interface3_log,"interface3", tank_log, false, 0);
    
  // PMT domes  
  dome1_log = new G4LogicalVolume(dome_solid, Lucite, "dome1_log", 0, 0, 0);
  dome2_log = new G4LogicalVolume(dome_solid, Lucite, "dome2_log", 0, 0, 0);
  dome3_log = new G4LogicalVolume(dome_solid, Lucite, "dome3_log", 0, 0, 0);
    
  dome1_phys = new G4PVPlacement(0, G4ThreeVector(dome1_x, dome1_y, dome1_z), 
                        dome1_log, "dome1", tank_log, false, 0);
  dome2_phys = new G4PVPlacement(0, G4ThreeVector(dome2_x, dome2_y, dome2_z), 
                        dome2_log, "dome2", tank_log, false, 0);
  dome3_phys = new G4PVPlacement(0, G4ThreeVector(dome3_x, dome3_y, dome3_z), 
                        dome3_log, "dome3", tank_log, false, 0);

  // Register the PMT's as sensitive detectors
  /*
    Note that each PMT must be given a unique index. If programmer
    screws up and gives the same index to 2 (or more) PMT's
    then photoelectrons will be double counted in that PMT.
    Perhaps there is a better way to do this.
  */
  
  SDName = "/Tank/pmt1";
  G4TankPMTAction *PMT1SD = new G4TankPMTAction(SDName, 1);
  
  SDName= "/Tank/pmt2";                 
  G4TankPMTAction *PMT2SD = new G4TankPMTAction(SDName, 2);
  
  SDName= "/Tank/pmt3";                 
  G4TankPMTAction *PMT3SD = new G4TankPMTAction(SDName, 3);
  
  SDMan->AddNewDetector(PMT1SD);
  SDMan->AddNewDetector(PMT2SD);
  SDMan->AddNewDetector(PMT3SD);
  
  pmt1_log->SetSensitiveDetector(PMT1SD);
  pmt2_log->SetSensitiveDetector(PMT2SD);
  pmt3_log->SetSensitiveDetector(PMT3SD);
     

  // Solar panel
  // -----------    
  
  if (fSolarPanelEnable) 
    {   
      solarPanel_solid = new G4Box("SolarPanel", fSolarPanelThickness/2, fSolarPanelWidth/2,
                    fSolarPanelLength/2);
    
      solarPanel_log = new G4LogicalVolume(solarPanel_solid, 
                              Al,               // ** MAY BE CHANGED
                              "solarPanel_log");
      G4RotationMatrix *solarPanelRot = new G4RotationMatrix();
      solarPanelRot->rotate(M_PI/2, G4ThreeVector(0, 0, 1));
      solarPanelRot->rotate(fSolarPanelTiltAngle, G4ThreeVector(0, 1, 0));
    
      solarPanel_phys = new G4PVPlacement(solarPanelRot, 
                            G4ThreeVector(fSolarPanelX, fSolarPanelY, 
                                          fSolarPanelZ),
                            "solarPanel",solarPanel_log, expHall_phys, 
                            false, 0);
    }
    
  // Chunk of Aluminum representing electronics box
  //-----------------------------------------------

  if (fElecBoxEnable) 
    {   
      elecBox_solid = new G4Box("ElecBox", fElecBoxThickness/2, fElecBoxWidth/2,
                    fElecBoxLength/2);
    
      elecBox_log = new G4LogicalVolume(elecBox_solid, 
                              Al,               // ** MAY BE CHANGED
                              "elecBox_log");
      G4RotationMatrix *elecBoxRot = new G4RotationMatrix();
      elecBoxRot->rotate(M_PI/2, G4ThreeVector(0, 0, 1));
      elecBoxRot->rotate(fElecBoxTiltAngle, G4ThreeVector(0, 1, 0));
    
      elecBox_phys = new G4PVPlacement(elecBoxRot, 
                            G4ThreeVector(fElecBoxX, fElecBoxY, fElecBoxZ),
                            "elecBox",elecBox_log, expHall_phys, false, 0); 
    }
   
  // ASCII scintillator
  //-----------------------------------------------
  if ( fASCIIEnable )
    {
      fScinHeightAboveGround=1500.*CLHEP::mm;
      fEnclosureThickness=5*CLHEP::mm;
      fShieldThickness=0.6*CLHEP::mm;
      fSpacerThickness=20*CLHEP::mm;

      fScinLength= 1800.*CLHEP::mm ;
      fScinWidth= 80*CLHEP::mm;
      fnScinBars= 27;
      fScinHeight= 10.*CLHEP::mm;
      
      fScinPos_x = 0.0*CLHEP::m;
      fScinPos_y = 0.0*CLHEP::m;  

      double Width=double(fnScinBars)*fScinWidth;
  
      G4VPhysicalVolume *pMother=expHall_phys;
      double zPosScin= fScinHeightAboveGround ;
      
      // Enclosure (PVC)
      enclosure_solid= new G4Box("enclosure_solid",fScinLength/2+fEnclosureThickness,  Width/2.+fEnclosureThickness,  fScinHeight/2.+fEnclosureThickness );
      enclosure_log = new G4LogicalVolume(enclosure_solid, Plastic, "enclosure_log", 0, 0, 0);  
      enclosure_phys = new G4PVPlacement(0, G4ThreeVector(fScinPos_x, fScinPos_y, zPosScin),
                                         "enclosure", enclosure_log,pMother, false, 0 );  

      
      // Scintillator   
      scin_solid= new G4Box("scin_solid",fScinLength/2,fScinWidth/2.,fScinHeight/2.);
      scin_log = new G4LogicalVolume(scin_solid, Polystyrene, "scin_log", 0, 0, 0);  
      for (int ibar=0;ibar<fnScinBars;ibar++)
        {    
          double Y=-Width/2.+fScinWidth*(double(ibar)+0.5);
          scin_phys = new G4PVPlacement(0, G4ThreeVector(0,Y,0),
                                        "scin", scin_log, enclosure_phys, true, ibar );
        }
      
      //Shield ( Above/Below enclosure for light tightness)
      shield_solid= new G4Box("shield_solid",fScinLength/2+fEnclosureThickness,  Width/2.+fEnclosureThickness,  fShieldThickness/2. );
      shield_log  = new G4LogicalVolume(shield_solid, theIron, "shield_log", 0, 0, 0);
  
      for (int i=0;i<2;i++)
        {
          double sign=(i==0?1.:-1.);
          double zPos= zPosScin + sign*( fScinHeight/2. + fEnclosureThickness +  fShieldThickness/2. );
          shield_phys[i] = new G4PVPlacement(0, G4ThreeVector(fScinPos_x, fScinPos_y, zPos ),
                                             (i==0?"shieldTop0":"shieldBot"), shield_log, pMother, true, i );
        }
  
      //Shield ( for air circulation)
      double zPos= zPosScin + fScinHeight/2. + fEnclosureThickness +  fShieldThickness + fSpacerThickness + fShieldThickness/2. ;
      shield_phys[2] = new G4PVPlacement(0, G4ThreeVector(fScinPos_x, fScinPos_y, zPos ),
                                         "shieldTop1", shield_log, pMother, true, 2 );
   
      // Register scintillator bars as sensitive detectors
      G4String Name = "/Tank/ASCII_scin";      
      G4ASCIIAction *ASCII_SCIN = new G4ASCIIAction(Name);  
      SDMan->AddNewDetector(ASCII_SCIN);
      scin_log->SetSensitiveDetector(ASCII_SCIN);
       
    }

  return;
}
  

void G4TankConstruction::GetDataForThisTank()
{
  SetDetectorParameters(); // updates G4TankConstruction 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();
  
  CreatePolystyrene();
  CreateIron();
  CreatePlastic();


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