// The includes to draw from other folders in ChoozSimulation
#include "DetectorConstruction.hh"
#include "DetectorMessenger.hh"

#include "G4Material.hh"
#include "G4NistManager.hh"
#include "G4Box.hh"
#include "G4Cons.hh"
#include "G4Tubs.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4PVReplica.hh"
#include "G4GeometryManager.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4SolidStore.hh"
#include "G4VisAttributes.hh"
#include "G4Colour.hh"
#include "G4SystemOfUnits.hh"
#include "G4RunManager.hh"
#include "G4RotationMatrix.hh"
#include "G4Element.hh"
#include "G4SDManager.hh"
#include "CrystalSensitiveDetector.hh"

#include <cstdlib>
#include <cmath>

using namespace std;

// Macro Command calling things to be built

DetectorConstruction::~DetectorConstruction()
{ delete fDetectorMessenger;}

G4VPhysicalVolume*  DetectorConstruction::Construct()
{
	InitializeWorld();
	ConstructPit();
	ConstructWaterSheilding();
	ConstructDespereaux();
	phystest = ConstructDetectors();
  
	return phystest;
}

void DetectorConstruction::DefineMaterials()
{
		// use G4-NIST materials data base
		G4NistManager* man = G4NistManager::Instance();
		fMaterialGalactic = man->FindOrBuildMaterial("G4_Galactic");
		fMaterialGe = man->FindOrBuildMaterial("G4_Ge");
		fMaterialSi = man->FindOrBuildMaterial("G4_Si");
	fMaterialAir = man->FindOrBuildMaterial("G4_AIR");
	fMaterialPoly = man->FindOrBuildMaterial("G4_POLYETHYLENE");
	fMaterialConcrete = man->FindOrBuildMaterial("G4_CONCRETE");
	fMaterialWater = man->FindOrBuildMaterial("G4_WATER");
	fMaterialLead = man->FindOrBuildMaterial("G4_Pb");
	fMaterialOs = man->FindOrBuildMaterial("G4_Os");
	fMaterialAl = man->FindOrBuildMaterial("G4_Al");
	fMaterialW = man->FindOrBuildMaterial("G4_W");
	fMaterialCu = man->FindOrBuildMaterial("G4_Cu");
	fMaterialPb = man->FindOrBuildMaterial("G4_Pb");
	fMaterialZn = man->FindOrBuildMaterial("G4_Zn");
	fMaterialZr = man->FindOrBuildMaterial("G4_Zr");  
	fMaterialGd = man->FindOrBuildMaterial("G4_Gd"); 

	elH = man->FindOrBuildElement("H");
	elO = man->FindOrBuildElement("O");
	elSi = man->FindOrBuildElement("Si");
	elMg = man->FindOrBuildElement("Mg");
	elAl = man->FindOrBuildElement("Al");
	elFe = man->FindOrBuildElement("Fe");
	elK = man->FindOrBuildElement("K");
	elB = man->FindOrBuildElement("B");
	elC = man->FindOrBuildElement("C");
  

	fMaterialScintillator = new G4Material("Scintillator",1.032*g/cm3,2);
	fMaterialScintillator->AddElement(elC,9);
	fMaterialScintillator->AddElement(elH,10);
  
	Chrysotile = new G4Material("Chrysotile",2.53*g/cm3,4);
	Chrysotile->AddElement(elMg,3);
	Chrysotile->AddElement(elSi,2);
	Chrysotile->AddElement(elO,9);
	Chrysotile->AddElement(elH,4);  
	fMaterialPoly_Borated = new G4Material("BoratedPoly",1.07*g/cm3,2);
	fMaterialPoly_Borated->AddMaterial(fMaterialPoly,95.0*perCent);
	fMaterialPoly_Borated->AddElement(elB,5.0*perCent);
  
	fMaterialQuartz = new G4Material("Quartz",2.81*g/cm3,2);
	fMaterialQuartz->AddElement(elSi,1);
	fMaterialQuartz->AddElement(elO,2);
  
	fMaterialCorundum = new G4Material("Corundum",2.81*g/cm3,2);
	fMaterialCorundum->AddElement(elAl,2);
	fMaterialCorundum->AddElement(elO,3);
  
	fMaterialIronOxide = new G4Material("IronOxide",2.81*g/cm3,2);
	fMaterialIronOxide->AddElement(elFe,1);
	fMaterialIronOxide->AddElement(elO,1);
  
	fMaterialMgOxide = new G4Material("MgOxide",2.81*g/cm3,2);
	fMaterialMgOxide->AddElement(elMg,1);
	fMaterialMgOxide->AddElement(elO,1);
  
	fMaterialKOxide = new G4Material("KOxide",2.81*g/cm3,2);
	fMaterialKOxide->AddElement(elK,2);
	fMaterialKOxide->AddElement(elO,1);

	fMaterialChoozRock = new G4Material("ChoozRock",2.81*g/cm3,5);
	fMaterialChoozRock->AddMaterial(fMaterialQuartz,58.0*perCent);
	fMaterialChoozRock->AddMaterial(fMaterialCorundum,19.0*perCent);
	fMaterialChoozRock->AddMaterial(fMaterialIronOxide,17.0*perCent);
	fMaterialChoozRock->AddMaterial(fMaterialMgOxide,4.0*perCent);
	fMaterialChoozRock->AddMaterial(fMaterialKOxide,2.0*perCent);
  
	mineraloil = new G4Material("MineralOil",0.86*g/cm3,2);
	mineraloil->AddElement(elC,25);
	mineraloil->AddElement(elH,43);  

	fMaterialSteel = new G4Material("Steel",7.9*g/cm3,2);
	fMaterialSteel->AddElement(elFe,98.0*perCent);
	fMaterialSteel->AddElement(elC,2.0*perCent);  

	fMaterialDodecane = new G4Material("Dodecane",0.75*g/cm3,2);
	fMaterialDodecane->AddElement(elC,12);
	fMaterialDodecane->AddElement(elH,24);  

	fMaterialPXE = new G4Material("PXE",0.988*g/cm3,2);
	fMaterialPXE->AddElement(elC,16);
	fMaterialPXE->AddElement(elH,18);

	fMaterialPPO = new G4Material("PPO",1.06*g/cm3,3);
	fMaterialPPO->AddElement(elC,3);
	fMaterialPPO->AddElement(elH,6);
	fMaterialPPO->AddElement(elO,1);
  
	fMaterialNuTarget = new G4Material("NuTargetLiquid",804.59*kg/m3,4);
	fMaterialNuTarget->AddMaterial(fMaterialDodecane,74.5*perCent);
	fMaterialNuTarget->AddMaterial(fMaterialPXE,24.6*perCent);
	fMaterialNuTarget->AddMaterial(fMaterialPPO,0.8*perCent);
	fMaterialNuTarget->AddMaterial(fMaterialGd,0.1*perCent);

	fMaterialGammaCatcher = new G4Material("GammaCatcherLiquid",800.6*kg/m3,3); 
	fMaterialGammaCatcher->AddMaterial(fMaterialDodecane,74.94*perCent);
	fMaterialGammaCatcher->AddMaterial(fMaterialPXE,24.68*perCent);
	fMaterialGammaCatcher->AddMaterial(fMaterialPPO,0.38*perCent);
  
}


void DetectorConstruction::InitializeWorld()
{

	fWorldSize = 20.*m;
	fSolidWorld = new G4Box("fSolidWorld", fWorldSize/2., fWorldSize/2., fWorldSize/2.);
	fLogicWorld = new G4LogicalVolume(fSolidWorld, fMaterialGalactic, "world");
	fLogicWorld->SetVisAttributes(G4VisAttributes::GetInvisible());
	fPhysiWorld = new G4PVPlacement(0, G4ThreeVector(), fLogicWorld, "world", 0, false, 0);

}  


void DetectorConstruction::ConstructPit()
{

	if(!fPhysiWorld)
	{
		G4cout << "ERROR in DetectorConstruction::ConstructPit" << G4endl;
		G4cout << "World volume does not exist!!" << G4endl;
		return;
	}  
	G4ThreeVector vec_offset(0*mm,0*mm,750*mm);
	//G4ThreeVector concrete_offset(0*mm,0*mm,100*mm);
	zeroradius = 0.*cm;
	startAngle = 0.*deg;
	spanningAngleFull = 360.*deg;
	G4Box *solidRock = new G4Box("Rock",6.*m,6.*m,5.*m);
	G4Tubs *PitTube = new G4Tubs("Pittube", zeroradius, 4350*mm, 8500*mm/2.0, startAngle, spanningAngleFull);
	G4Tubs *Jacket2Tube = new G4Tubs("JacketShell21", zeroradius, 4250*mm, 8400*mm/2.0,startAngle, spanningAngleFull);
	G4Tubs *Jacket21Tube = new G4Tubs("JacketShell21", zeroradius, 4350*mm, 8500*mm/2.0,startAngle, spanningAngleFull);
	this->PitSolid = new G4SubtractionSolid("ChoozPit",solidRock,PitTube,0,vec_offset);
	this->PitLog = new G4LogicalVolume(PitSolid, fMaterialChoozRock, "logicalRock");
	this->PitPhys =  new G4PVPlacement(0, vec_offset, PitLog, "PVRock", fLogicWorld, false, 0);
	this->PitJacketSolid = new G4SubtractionSolid("ChoozPit_ConcreteJacket",Jacket21Tube,Jacket2Tube,0,concrete_offset);
	this->PitJacketLog = new G4LogicalVolume(PitJacketSolid, fMaterialConcrete, "logicalConcrete");
	this->PitJacketPhys =  new G4PVPlacement(0,G4ThreeVector(0), PitJacketLog, "ConcreteJacket", fLogicWorld, false, 0);

	// Visualization
	G4VisAttributes visRock(G4Colour(0.398,0.199,0.));
	visRock.SetForceWireframe(true);
	PitLog->SetVisAttributes(visRock);
	G4VisAttributes visJacket(G4Colour(0.625,0.625,0.625));
	visJacket.SetForceWireframe(true);
	PitJacketLog->SetVisAttributes(visJacket);

}


void DetectorConstruction::ConstructWaterShielding()
{

	zeroradius = 0.*cm;
	startAngle = 0.*deg;
	spanningAngleFull = 360.*deg;
	G4ThreeVector vec_zero(0*mm,0*mm,50*mm); // Offset include to have top of water shield line up with top of pit
	G4ThreeVector vec_offset(0*mm,0*mm,250*mm);
	this->TempTube_inner = new G4Tubs("Temp_inner", zeroradius, 1150.0*mm, 3200.0*mm,startAngle, spanningAngleFull);
	this->TempTube_outer = new G4Tubs("Temp_outer", zeroradius, 4250.0*mm, 4200*mm,startAngle, spanningAngleFull);
	this->WaterTubeSolid = new G4SubtractionSolid("WaterShielding",TempTube_outer,TempTube_inner,0,vec_offset);
	this->WaterTubeLog = new G4LogicalVolume(WaterTubeSolid, fMaterialWater, "WaterShielding");
	this->WaterTubePhys = new G4PVPlacement(0,vec_zero, WaterTubeLog, "WaterShielding",fLogicWorld, false,0);
	
	// Visualization 
	G4VisAttributes visWaterTube(G4Colour(0,0,1));
	visWaterTube.SetForceWireframe(true);
	visWaterTube.SetForceAuxEdgeVisible(true);
	WaterTubeLog->SetVisAttributes(visWaterTube); 

}

void DetectorConstruction::ConstructPbSheilding(G4double innerR, G4double outerR,G4double topthickness)
{  

	// detector shielding is arbitrarily assuemd to be two concentric cylinders of
	// polyethylene and lead; poly is 10cm thick and lead is 4cm thick (again arbitrarily)
 	// Lead Tube
	zeroradius = 0.*cm;
	startAngle = 0.*deg;
	spanningAngleFull = 360.*deg;
	G4ThreeVector vec_zero(0*mm,0*mm,0*mm);
	outerRadius_shield = 1.0*outerR*cm;
	//fDRheight = 1279*mm;
	outerRadius_shield1 = 1.0*innerR*cm;
	fDRheight1 = 1287*mm;
	this->shieldTube = new G4Tubs("Polyshield", zeroradius, outerRadius_shield, fDRheight/2.0, startAngle, spanningAngleFull);
	this->shield1Tube = new G4Tubs("Polyshield", zeroradius, outerRadius_shield1, fDRheight1/2.0,startAngle, spanningAngleFull);
	G4SubtractionSolid* fullshieldTube = new G4SubtractionSolid("Shielding Shell (Pb)",shieldTube,shield1Tube);
	this->fullshieldTubeLog = new G4LogicalVolume(fullshieldTube, fMaterialPb, "Shielding Shell (Pb)");
	this->fullshieldTubePhys = new G4PVPlacement(0,vec_zero, fullshieldTubeLog, "Shielding shell (Pb)",fLogicWorld, false,0);

	// Visualization  
	G4VisAttributes visshieldTube(G4Colour(1.0,0.586,0.0));
	visshieldTube.SetForceWireframe(true);
	visshieldTube.SetForceAuxEdgeVisible(true);
	fullshieldTubeLog->SetVisAttributes(visshieldTube);

}


void DetectorConstruction::ConstructPolySheilding(G4double innerR, G4double outerR,G4double topthickness)
{  
	// detector shielding is arbitrarily assumed to be two concentric cylinders of
	// polyethylene and lead; poly is 10cm thick and lead is 4cm thick (again arbitrarily)
	// Lead Tube
	zeroradius = 0.*cm;
	startAngle = 0.*deg;
	spanningAngleFull = 360.*deg;
	G4ThreeVector vec_zero(0*mm,0*mm,0*mm);
	outerRadius_shield = 1.0*outerR*cm;

	outerRadius_shield1 = 1.0*innerR*cm;
	fDRheight= 1279*mm;
	fDRheight1 = 1289*mm;
	this->shieldTube = new G4Tubs("Polyshield", zeroradius, outerRadius_shield, fDRheight/2.0, startAngle, spanningAngleFull);
	this->shield1Tube = new G4Tubs("Polyshield", zeroradius, outerRadius_shield1, fDRheight1/2.0,startAngle, spanningAngleFull);
	G4SubtractionSolid* fullshieldTube = new G4SubtractionSolid("Borated_Poly_Shield",shieldTube,shield1Tube);
	this->fullshieldTubeLog = new G4LogicalVolume(fullshieldTube, fMaterialPoly_Borated, "Borated_Poly_Shield");
	this->fullshieldTubePhys = new G4PVPlacement(0,vec_zero, fullshieldTubeLog, "Borated_Poly_Shield",fLogicWorld, false,0);
  
	// Visualization
	G4VisAttributes visshieldTube(G4Colour(1.0,0.586,0.0));
	visshieldTube.SetForceWireframe(true);
	visshieldTube.SetForceAuxEdgeVisible(true);
	fullshieldTubeLog->SetVisAttributes(visshieldTube);

}

void DetectorConstruction::ConstructSupportiveCasing()
{

	// Suportive Casing

	// Top Hardwear

}


void DetectorConstruction::ConstructDespereaux()
{

	// 50 K Plate
	outerRadius50 = 180*mm;
	hz50 = 3*mm;
	zeroradius = 0.*mm;
	startAngle = 0.*deg;
	startAngle90 = 90.*deg;
	spanningAngleFull = 360.*deg;
	G4ThreeVector vec_fifty(0*mm,0*mm,403.15*mm);

	this->plate50KTube = new G4Tubs("50KPlate", zeroradius, outerRadius50, hz50/2., startAngle, spanningAngleFull);
	this->plate50KLog = new G4LogicalVolume(plateTube, fMaterialCu, "50 K Plate");
	this->plate50KPhys = new G4PVPlacement(0, vec_fifty, plate50KLog, "50 K Plate", fLogicWorld, false, 0);
	// G4VisAttributes visplate(G4Colour(1.,1.,0.));
	// visplate.SetForceWireframe(true);
	// visplate.SetForceSolid(true);
	// plate50KLog->SetVisAttributes(visplate);

	// 50 K Associated Cylinder
	outerRadius50KC = 169.5*mm
	innerRadius50KC =167.5*mm
	hz50KCouter = 832.4*mm
	hz50KCinner =830.4*mm
	  G4ThreeVector vec_fiftyc(0*mm,0*mm,-59.15*mm)
	this->cylinder50KTubeOuter = new G4Tubs("50KCylinderouter", zeroradius, outerRadius50KC, hz50KCouter/2., startAngle, spanningAngleFull);
	this->cylinder50KTubeInner = new G4Tubs("50KCylinderinner", zeroradius, innerRadius50KC, hz50KCinner/2., startAngle, spanningAngleFull);
	this->G4SubtractionSolid*hollowCylinder50K = new G4SubtractionSolid("HollowCylinder50K", cylinder50KTubeOuter, cylinder50KTubeInner);
	this->hollowCylinder50KLog = new G4LogicalVolume(plateTube, fMaterialCu, "Hollow 50 K Cylinder");
	this->hollowCylinder50KPhys = new G4PVPlacement(0, vec_fiftyc, hollowCylinder50KLog, "Hollow 50 K Cylinder", FLogicWorld, false, 0);
	// Visualization
	// G4VisAttributes visplate(G4Colour(1.,1.,0.));
	// visplate.SetForceWireframe(true);
	// visplate.SetForceSolid(true);
	// hollowCylinder50KLog->SetVisAttributes(visplate);

	// 4 K Plate
	outerRadius4K = 162*mm;
	hz4K = 2.5*mm;
	G4ThreeVector vec_four(0*mm,0*mm,259.1*mm);
	
	this->plate4KTube = new G4Tubs("4KPlate", zeroradius, outerRadius4K, hz4K/2., startAngle, spanningAngleFull);
	this->plate4KLog = newG4LogicalVolume(plateTube, fMaterialCu, "4 K Plate", fLogicWorld, false, 0);
	this->plate4KPhys = new G4PVPlacement(0,vec_four, plate4KLog, "4 K Plate", fLogicWorld, false, 0);
	// Visualization Stuff
	// G4VisAttributes visplate(G4Colour(1.,1.,0.));
	// visplate.SetForceWireframe(true);
	// visplate.SetForceSolid(true);
	// plate4KLog->SetVisAttributes(visplate);

	// 4 K Associated Cylinder
	outerRadius4KC = 151*mm
	innerRadius4KC =149*mm
	hz4KCouter = 618.2*mm
	hz4KCinner =616.2*mm
	G4ThreeVector vec_fourc(0*mm,0*mm,-50*mm)
	this->cylinder4KTubeOuter = new G4Tubs("4KCylinderouter", zeroradius, outerRadius4KC, hz4KCouter/2., startAngle, spanningAngleFull);
	this->cylinder4KTubeInner = new G4Tubs("4KCylinderinner", zeroradius, innerRadius4KC, hz4KCinner/2., startAngle, spanningAngleFull);
	this->G4SubtractionSolid*hollowCylinder4K = new G4SubtractionSolid("HollowCylinder4K", cylinder4KTubeOuter, cylinder4KTubeInner);
	this->hollowCylinder4KLog = new G4LogicalVolume(plateTube, fMaterialCu, "Hollow 4 K Cylinder");
	this->hollowCylinder4KPhys = new G4PVPlacement(0, vec_fourc, hollowCylinder4KLog, "Hollow 4 K Cylinder", FLogicWorld, false, 0);
	// Visualization
	// G4VisAttributes visplate(G4Colour(1.,1.,0.));
	// visplate.SetForceWireframe(true);
	// visplate.SetForceSolid(true);
	// hollowCylinder4KLog->SetVisAttributes(visplate); 


	// 1 K Plate
	outerRadius1K = 142*mm;
	hz1K = 2.5*mm;
	G4ThreeVector vec_one(0*mm,0*mm,128.1*mm);

	this->plate1KTube = new G4Tubs("1KPlate", zeroradius, outerRadius1K, hz1K/2., startAngle, spanningAngleFull);
	this->plate1KLog = newG4LogicalVolume(plateTube, fMaterialCu, "1 K Plate", fLogicWorld, false, 0);
	this->plate1KPhys = new G4PVPlacement(0,vec_one, plate4KLog, "1 K Plate", fLogicWorld, false, 0);
	// Visualization Stuff
	// G4VisAttributes visplate(G4Colour(1.,1.,0.));
	// visplate.SetForceWireframe(true);
	// visplate.SetForceSolid(true);
	// plate1KLog->SetVisAttributes(visplate);

	// 1 K Associated Cylinder
	outerRadius1KC = 131.5*mm
	innerRadius1KC =129.5*mm
	hz1KCouter = 504.2*mm
	hz1KCinner =502.2*mm
	G4ThreeVector vec_onec(0*mm,0*mm,-294.2*mm)
	this->cylinder1KTubeOuter = new G4Tubs("1KCylinderouter", zeroradius, outerRadius1KC, hz1KCouter/2., startAngle, spanningAngleFull);
	this->cylinder1KTubeInner = new G4Tubs("1KCylinderinner", zeroradius, innerRadius1KC, hz1KCinner/2., startAngle, spanningAngleFull);
	this->G4SubtractionSolid*hollowCylinder1K = new G4SubtractionSolid("HollowCylinder1K", cylinder1KTubeOuter, cylinder1KTubeInner);
	this->hollowCylinder1KLog = new G4LogicalVolume(plateTube, fMaterialCu, "Hollow 1 K Cylinder");
	this->hollowCylinder1KPhys = new G4PVPlacement(0, vec_onec, hollowCylinder1KLog, "Hollow 1 K Cylinder", FLogicWorld, false, 0);
	// Visualization
	// G4VisAttributes visplate(G4Colour(1.,1.,0.));
	// visplate.SetForceWireframe(true);
	// visplate.SetForceSolid(true);
	// hollowCylinder1KLog->SetVisAttributes(visplate);

	// 100 mK Plate
	outerRadius100mK = 125*mm;
	hz100mK = 1.5*mm;
	G4ThreeVector vec_hundred(0*mm,0*mm,144.1*mm);

	this->plate100mKTube = new G4Tubs("100mKPlate", zeroradius, outerRadius100mK, hz100mK/2., startAngle, spanningAngleFull);
	this->plate100mKLog = newG4LogicalVolume(plateTube, fMaterialCu, "100 mK Plate", fLogicWorld, false, 0);
	this->plate100mKPhys = new G4PVPlacement(0,vec_hundred, plate4KLog, "100 mK Plate", fLogicWorld, false, 0);
	// Within 1K Associated Cylinder
	// Visualization Stuff
	// G4VisAttributes visplate(G4Colour(1.,1.,0.));
	// visplate.SetForceWireframe(true);
	// visplate.SetForceSolid(true);
	// plate100mKLog->SetVisAttributes(visplate);

	// 15 mK Plate
	outerRadius15mK = 125*mm;
	hz15mK = 2*mm;
	G4ThreeVector vec_fifteen(0*mm,0*mm,72.9*mm);

	this->plate15mKTube = new G4Tubs("15mKPlate", zeroradius, outerRadius15mK, hz15mK/2., startAngle, spanningAngleFull);
	this->plate15mKLog = newG4LogicalVolume(plateTube, fMaterialCu, "15 mK Plate", fLogicWorld, false, 0);
	this->plate15mKPhys = new G4PVPlacement(0,vec_fifteen, plate4KLog, "15 mK Plate", fLogicWorld, false, 0);
	// Within 1K Associated Cylinder
	// Visualization Stuff
	// G4VisAttributes visplate(G4Colour(1.,1.,0.));
	// visplate.SetForceWireframe(true);
	// visplate.SetForceSolid(true);
	// plate15mKLog->SetVisAttributes(visplate);

}

G4VPhysicalVolume* DetectorConstruction::ConstructDetectors()
{

	DetectorSize = 5.2*cm; // Size needed for 1 Kg of detector material
  
	std::vector<std::string> CyrstalLabels = {"Crystal_1", "Crystal_2","Crystal_3", "Crystal_4","Crystal_5"};
	std::vector<G4double> Cyrstalpos_x = {0.*cm, -6.*cm, 6.*cm, 6.*cm, -6.*cm};
	std::vector<G4double> Cyrstalpos_y = {0.*cm, 6.*cm, 6.*cm, -6.*cm, -6.*cm};
	std::vector<G4double> Cyrstalpos_z = {0.*cm, 0.*cm, 0.*cm, 0.*cm, 0.*cm};

	// Visualization
	G4VisAttributes visDetector(G4Colour(0.0,1.0,0.0));
	visDetector.SetForceWireframe(true);
	visDetector.SetForceAuxEdgeVisible(true);

	G4ThreeVector vec_zero(0*mm,0*mm,0*mm);
	G4ThreeVector vec_offset(0*mm,0*mm,-192.1*mm);
	G4SDManager* SDmanager = G4SDManager::GetSDMpointer();
	for(size_t jCrys = 0; jCrys < 1; jCrys++)
	{
		v_CrystalBoxes.push_back(new G4Box(CyrstalLabels[jCrys],DetectorSize/2.0,DetectorSize/2.0,DetectorSize/2.0));
		v_CrystalBoxesLog.push_back(new G4LogicalVolume(v_CrystalBoxes[jCrys], fMaterialOs,CyrstalLabels[jCrys]));
		G4ThreeVector crys_pos(Cyrstalpos_x[jCrys],Cyrstalpos_y[jCrys],Cyrstalpos_z[jCrys]);
		v_CrystalBoxesPhys.push_back(new G4PVPlacement(0,crys_pos+vec_offset,v_CrystalBoxesLog[jCrys],CyrstalLabels[jCrys],fLogicWorld,false,0));
		v_CrystalBoxesLog[jCrys]->SetVisAttributes(visDetector);
		if (jCrys==0)
		{  
			CrystalSensitiveDetector *Crystal1_SD = new CrystalSensitiveDetector("Crystal_1SD");
			SDmanager->AddNewDetector(Crystal1_SD);
			v_CrystalBoxesLog[jCrys]->SetSensitiveDetector(Crystal1_SD);
		}  
  
	}
	return fPhysiWorld; 

}

G4VPhysicalVolume*  DetectorConstruction::ConstructSingleDetector()
{

	// arbitrarily assuming 4-inch iZIP here [AJA]
	G4ThreeVector vec_zero(0*mm,0*mm,0*mm);
	this->crystalBox_single = new G4Box("Single Crystal",DetectorSize/2.0,DetectorSize/2.0,DetectorSize/2.0);
	this->crystalLog_single = new G4LogicalVolume(crystalBox_single,fMaterialZn, "Single Crystal");
	this->crystalPhys_single = new G4PVPlacement(0, vec_zero, crystalLog_single, "Single Crystal",fLogicWorld, false,0);
	
	// Visualization
	G4VisAttributes visDetector(G4Colour(0.0,1.0,0.0));
	visDetector.SetForceWireframe(true);
	visDetector.SetForceAuxEdgeVisible(true);
	crystalLog_single->SetVisAttributes(visDetector); 
	return fPhysiWorld;  

}


void DetectorConstruction::SetCrystalMaterial(G4String Material)
{
	if (Material == "Os" || Material == "Zn" || Material == "Zr")
	{
		crystalmaterial = Material;
		if (Material == "Os")
		{
			for (int i=0;i<1;i++)
			{  
				v_CrystalBoxesLog[i]->SetMaterial(fMaterialOs);
			}  
		}
		if (Material == "Zn")
		{
			for (int i=0;i<1;i++)
			{  
				v_CrystalBoxesLog[i]->SetMaterial(fMaterialZn);
			}
		}  
		if (Material == "Zr")
		{
			for (int i=0;i<1;i++)
			{  
				v_CrystalBoxesLog[i]->SetMaterial(fMaterialZr);
			}
		}   
	}
	else
	{
		std::cerr << "ERROR in SetCrystalMaterial source type must be Os or Zn setting to Os!" << std::endl;
		crystalmaterial = "Os";
		for (int i=0;i<1;i++)
		{  
			v_CrystalBoxesLog[i]->SetMaterial(fMaterialOs);
		}     
	}
	G4RunManager::GetRunManager()->GeometryHasBeenModified();

}