#include "DetectorConstruction.hh"
#include "G4VUserDetectorConstruction.hh"
#include "G4VPhysicalVolume.hh"
#include "G4Material.hh"
#include "G4Isotope.hh"
#include "G4Element.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4SystemOfUnits.hh"
#include "G4NistManager.hh"
#include "G4SDManager.hh"
#include "G4MultiFunctionalDetector.hh"
#include "G4PSEnergyDeposit.hh"
#include "G4PSNofSecondary.hh"
#include "G4VPrimitiveScorer.hh"
#include "MySensitiveDetector.hh"
#include "G4SubtractionSolid.hh"
#include "G4MTRunManager.hh"

using namespace CLHEP;

DetectorConstruction::DetectorConstruction() : G4VUserDetectorConstruction() {}

DetectorConstruction::~DetectorConstruction() {}

// Construct method
G4VPhysicalVolume* DetectorConstruction::Construct() {

    // Get nist material manager
    G4NistManager* nist = G4NistManager::Instance();
    nist->SetVerbose(1);

    // Individual isotopes, elements, and materials
    G4Material *Air = nist->FindOrBuildMaterial("G4_AIR");
    G4Material *Hg = nist->FindOrBuildMaterial("G4_Hg");
    
    G4int iz, n;                       //iz=number of protons  in an isotope;
                                     // n=number of nucleons in an isotope;
    G4int ncomponents;                              
    G4double a;
    G4double li6density;       
    G4double abundance;                      
    G4int natoms;               

    G4Isotope* Li6 = new G4Isotope("Li6", iz=3, n=6, a=6.01512*g/mole);
    G4Element* Li  = new G4Element("Pure Lithium 6", "Li", ncomponents=1);
         Li->AddIsotope(Li6, abundance= 100.*perCent);
    G4Material* Lithium6 = new G4Material("Lithium 6", li6density= 1.000*g/cm3, ncomponents=1);
         Lithium6->AddElement(Li, natoms=1);
    
    G4bool checkOverlaps = true;

    // World
    G4double worldSize = 20 * m;
    G4Box* worldSolid = new G4Box("World", 0.5 * worldSize, 0.5 * worldSize, 0.5 * worldSize);
    G4LogicalVolume* worldLogical = new G4LogicalVolume(worldSolid, G4Material::GetMaterial("G4_AIR"), "World");
    G4VPhysicalVolume* worldPhysical = new G4PVPlacement(0, G4ThreeVector(), worldLogical, "World", 0, false, 0, true);

    // Mercury Target
    G4double targetRadius = 10 * cm;
    G4double targetLength = 10 * cm;

    G4Tubs* targetSolid = new G4Tubs("Target", 0, targetRadius, 0.5 * targetLength, 0, 2 * pi);
    G4LogicalVolume* targetLogical = new G4LogicalVolume(targetSolid, G4Material::GetMaterial("G4_Hg"), "Target");
    new G4PVPlacement(0, G4ThreeVector(), targetLogical, "Target", worldLogical, false, 0, true);

    // Variables for additional radius (to compensate for space)
    G4double extraRadius = 1.0 * cm;

    // Variables for initial outer detector shell
    G4double neutronDetectorThickness = 1.0 * cm;
    G4double neutronDetectorRadius = targetRadius + neutronDetectorThickness + extraRadius;
    G4double neutronDetectorLength = targetLength + (4 * neutronDetectorThickness);

    // Variables for cavity
    G4double neutronDetectorCavityRadius = 11.0 * cm;
    G4double neutronDetectorCavityLength = 12.0 * cm;

    // Initial detector geometry (cylinder without any holes)
    G4Tubs* initialDetectorSolid = new G4Tubs("Detector", 0.0, neutronDetectorRadius, neutronDetectorLength / 2.0, 0.0, 2.0 * M_PI);
    G4Material* detectorMaterial = G4Material::GetMaterial("G4_POLYETHYLENE");

    G4LogicalVolume* initialDetectorLogical = new G4LogicalVolume(initialDetectorSolid, detectorMaterial, "InitialDetectorLogical");

    // Subtract the target cavity
    G4Tubs* cavitySolid = new G4Tubs("Cavity", 0, neutronDetectorCavityRadius, 0.5 * neutronDetectorCavityLength, 0, 2 * pi);
    G4SubtractionSolid* cavitatedDetectorSolid = new G4SubtractionSolid("CavitatedDetector", initialDetectorSolid, cavitySolid);
    neutronCavitatedDetectorLogical = new G4LogicalVolume(cavitatedDetectorSolid, Air, "neutronCavitatedDetectorLogical");
    new G4PVPlacement(0, G4ThreeVector(), neutronCavitatedDetectorLogical, "CavitatedDetectorPhysical", worldLogical, false, 0, true);

    // Calculate the dimensions of the LiF geometry
    G4double LIFRadius = 15.0 * cm;
    G4double LIFLength = neutronDetectorLength + 60;
    G4double LIFCavityRadius = 13.0 * cm;
    G4double LIFCavityLength = 16.0 * cm;

    // LiF geometry (scaled neutron detector)
    G4Tubs* LIFSolid = new G4Tubs("LIFSolid", 0.0, LIFRadius, LIFLength / 2.0, 0.0, 2.0 * M_PI);
    G4LogicalVolume* LIFLogical = new G4LogicalVolume(LIFSolid, Lithium6, "LIFLogical");

    // Subtract the target and detector cavity
    G4Tubs* LIFCavitySolid = new G4Tubs("LIFCavity", 0, LIFCavityRadius, 0.5 * LIFCavityLength, 0, 2 * pi);
    G4SubtractionSolid* cavitatedLIFSolid = new G4SubtractionSolid("CavitatedLIFSolid", LIFSolid, LIFCavitySolid);
    G4LogicalVolume* cavitatedLIFLogical = new G4LogicalVolume(cavitatedLIFSolid, Lithium6, "CavitatedLIFLogical");
    new G4PVPlacement(0, G4ThreeVector(), cavitatedLIFLogical, "CavitatedLIFPhysical", worldLogical, false, 0, true);

    // Variables for additional radius (to compensate for space)
    G4double extraAlphaDetectorRadius = 1.0 * cm;

    // Variables for initial outer detector shell
    G4double alphaDetectorThickness = 1.0 * cm;
    G4double alphaDetectorRadius = LIFRadius + alphaDetectorThickness + extraAlphaDetectorRadius;
    G4double alphaDetectorLength = LIFLength + 60;
    
    G4double alphaDetectorCavityRadius = 16.0 * cm;
    G4double alphaDetectorCavityLength =  LIFLength + 20;
    
    // Initial detector geometry (cylinder without any holes)
    G4Tubs* initialAlphaDetectorSolid = new G4Tubs("AlphaDetector", 0.0, alphaDetectorRadius, alphaDetectorLength / 2.0, 0.0, 2.0 * M_PI);
    G4LogicalVolume* initialALphaDetectorLogical = new G4LogicalVolume(initialAlphaDetectorSolid, detectorMaterial, "InitialAlphaDetectorLogical");

    // Subtract the target cavity
    G4Tubs* alphaCavitySolid = new G4Tubs("AlphaCavity", 0, alphaDetectorCavityRadius, 0.5 * alphaDetectorCavityLength, 0, 2 * pi);
    G4SubtractionSolid* cavitatedAlphaDetectorSolid = new G4SubtractionSolid("CavitatedAlphaDetector", initialAlphaDetectorSolid, alphaCavitySolid);
    alphaCavitatedDetectorLogical = new G4LogicalVolume(cavitatedAlphaDetectorSolid, Air, "alphaCavitatedDetectorLogical");
    new G4PVPlacement(0, G4ThreeVector(), alphaCavitatedDetectorLogical, "CavitatedAlphaDetectorPhysical", worldLogical, false, 0, true);

    // Return the physical world
    return worldPhysical;
}

void DetectorConstruction::ConstructSDandField()
{
    // Create a sensitive detector
    MySensitiveDetector* sensitiveDetector = new MySensitiveDetector("MyNeutronSensitiveDetector");
    MySensitiveDetector* alphaSensitiveDetector = new MySensitiveDetector("MyAlphaSensitiveDetector");
    
    // Set sensitive detector
    neutronCavitatedDetectorLogical->SetSensitiveDetector(sensitiveDetector);
    alphaCavitatedDetectorLogical->SetSensitiveDetector(alphaSensitiveDetector);
}