#include "GeneDetectorConstruction.hh"

#include "G4Box.hh"
#include "G4Colour.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4LogicalVolume.hh"
#include "G4Material.hh"
#include "G4NistManager.hh"
#include "G4OpticalSurface.hh"
#include "G4PVPlacement.hh"
#include "G4PhysicalConstants.hh"
#include "G4RotationMatrix.hh"
#include "G4SDManager.hh"
#include "G4SubtractionSolid.hh"
#include "G4SystemOfUnits.hh"
#include "G4ThreeVector.hh"
#include "G4Tubs.hh"
#include "G4VisAttributes.hh"
#include "GeneConst.hh"
#include "GeneHe3detSD.hh"
#include "globals.hh"

GeneDetectorConstruction::GeneDetectorConstruction() {}

GeneDetectorConstruction::~GeneDetectorConstruction() {}

G4VPhysicalVolume* GeneDetectorConstruction::Construct() {
    // ********** Materials **********
    G4NistManager* manager = G4NistManager::Instance();
    G4double z, n, a, density, temperature, pressure, fractionmass;
    G4String name, symbol;
    G4int ncomponents;

    // Hydrogen 氢
    G4Element* ele_H = manager->FindOrBuildElement("H");
    // Carbon 碳
    G4Element* ele_C = manager->FindOrBuildElement("C");
    // Boron 硼
    G4Material* mat_B = manager->FindOrBuildMaterial("G4_B");
    // Stainless stell 不锈钢
    G4Material* mat_stell = manager->FindOrBuildMaterial("G4_STAINLESS-STEEL");
    // Air 空气
    G4Material* mat_Air = manager->FindOrBuildMaterial("G4_AIR");
    // Polyethylene 聚乙烯
    G4Material* mat_PE = manager->FindOrBuildMaterial("G4_POLYETHYLENE");

    // shielding 屏蔽层 PE(B) 含硼聚乙烯
    density = 0.938 * g / cm3;
    G4Material* shield_PE = new G4Material("BPE", density, ncomponents = 2);
    shield_PE->AddMaterial(mat_PE, fractionmass = 93 * perCent);
    shield_PE->AddMaterial(mat_B, fractionmass = 7 * perCent);

    // Vacuum 真空
    temperature = 77 * kelvin;
    pressure = 1.0e-5 * pascal;
    density = 1.29 * kg / m3 * pressure / atmosphere;
    G4Material* Vacuum = new G4Material("Vacuum", density, ncomponents = 1, kStateGas, temperature, pressure);
    Vacuum->AddMaterial(mat_Air, 100.0 * perCent);

    // He3 gas He3气体
    G4Isotope* iso_He3 = new G4Isotope(name = "He3", z = 2.0, n = 3.0, a = 3.016 * g / mole);
    G4Element* ele_He = new G4Element(name = "iso_He3", symbol = "He3", ncomponents = 1);
    ele_He->AddIsotope(iso_He3, 100.0 * perCent);

    // Work gas 工作气体
    temperature = 300. * kelvin;
    pressure = 10. * atmosphere;
    density = 1.22685 * mg / cm3;  // 0.49074(3He, 4atm) Calculate by LISE++
    G4Material* work_gas = new G4Material("Workgas", density, ncomponents = 1, kStateGas, temperature, pressure);
    work_gas->AddElement(ele_He, 100. * perCent);

    // Ne22 Target Ne22靶
    temperature = 300 * kelvin;
    pressure = 100 * pascal;
    density = (0.981146 * kg / m3) * (pressure / atmosphere) * (273.15 * kelvin / temperature);
    G4Isotope* iso_Ne22 = new G4Isotope(name = "Ne22", z = 10, n = 12, a = 21.99138 * g / mole);
    G4Element* ele_Ne22 = new G4Element(name = "Neon-22", symbol = "Ne22", ncomponents = 1);
    ele_Ne22->AddIsotope(iso_Ne22, 100 * perCent);
    G4Material* target_Ne = new G4Material("targetNe", density, ncomponents = 1, kStateGas, temperature, pressure);
    target_Ne->AddElement(ele_Ne22, 100 * perCent);

    // ********** Sciintillator **********
    // Plastic scintillator 塑料闪烁体 EJ-200
    density = 1.023 * g / cm3;
    G4Material* ej_200 = new G4Material("ej_200", density, ncomponents = 2);
    ej_200->AddElement(ele_H, 8.5 * perCent);
    ej_200->AddElement(ele_C, 91.5 * perCent);

    /*
        Material properties table 物质光学属性表
        Ref:
            1. Geant4 Book For Application Developers, TrackingAndPhysics, physicsProcess, optical-photon-processes
            2. Eljen-Catalog: Ej-200
    */
    G4MaterialPropertiesTable* ej_200_MPT = new G4MaterialPropertiesTable();
    // Optical properties 光学属性
    const G4int nEntries = 21;
    G4double PhotonEnergy[nEntries] = {2.411 * eV, 2.483 * eV, 2.533 * eV, 2.586 * eV, 2.628 * eV, 2.666 * eV,
                                       2.699 * eV, 2.779 * eV, 2.821 * eV, 2.865 * eV, 2.903 * eV, 2.921 * eV,
                                       2.936 * eV, 2.946 * eV, 2.964 * eV, 3.008 * eV, 3.022 * eV, 3.039 * eV,
                                       3.079 * eV, 3.104 * eV, 3.143 * eV};  // 光子能量 395 nm -> 515 nm
    G4double AbsorptionLength[nEntries] = {380. * cm, 380. * cm, 380. * cm, 380. * cm, 380. * cm, 380. * cm,
                                           380. * cm, 380. * cm, 380. * cm, 380. * cm, 380. * cm, 380. * cm,
                                           380. * cm, 380. * cm, 380. * cm, 380. * cm, 380. * cm, 380. * cm,
                                           380. * cm, 380. * cm, 380. * cm};  // 吸收长度
    G4double RefractiveIndex[nEntries] = {1.58, 1.58, 1.58, 1.58, 1.58, 1.58, 1.58, 1.58, 1.58, 1.58, 1.58,
                                          1.58, 1.58, 1.58, 1.58, 1.58, 1.58, 1.58, 1.58, 1.58, 1.58};  // 折射率
    G4double FastComponent[nEntries] = {0.00000, 0.05885, 0.10053, 0.17197, 0.26426, 0.37440, 0.43394,
                                        0.60064, 0.74651, 0.89833, 0.98168, 1.00000, 0.97870, 0.93703,
                                        0.79712, 0.40119, 0.30444, 0.19876, 0.03801, 0.00926, 0.00000};  // 快成分能谱
    ej_200_MPT->AddProperty("ABSLENGTH", PhotonEnergy, AbsorptionLength, nEntries);
    ej_200_MPT->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex, nEntries);
    ej_200_MPT->AddProperty("FASTCOMPONENT", PhotonEnergy, FastComponent, nEntries);

    // Intrinsic properties 本征属性
    ej_200_MPT->AddConstProperty("RESOLUTIONSCALE", 10.0);                // 本征分辨率
    ej_200_MPT->AddConstProperty("SCINTILLATIONYIELD", 10000. / MeV);     // 本征光产额
    ej_200_MPT->AddConstProperty("YIELDRATIO", 1.0);                      // 快成分占比
    ej_200_MPT->AddConstProperty("FASTSCINTILLATIONRISETIME", 0.9 * ns);  // 快成分上升时间
    ej_200_MPT->AddConstProperty("FASTTIMECONSTANT", 2.1 * ns);           // 快成分衰减时间
    ej_200->SetMaterialPropertiesTable(ej_200_MPT);

    // Boundary Optical Properties 边界光学属性
    G4OpticalSurface* ej_200_OpticalSurface = new G4OpticalSurface("ej_200_OpticalSurface");
    ej_200_OpticalSurface->SetType(dielectric_metal);
    ej_200_OpticalSurface->SetFinish(polished);
    ej_200_OpticalSurface->SetModel(glisur);

    // ********** Geometry **********
    // World
    G4double world_xyz = 2. * m;
    G4Box* solidWorld = new G4Box("World", 0.5 * world_xyz, 0.5 * world_xyz, 0.5 * world_xyz);
    logicalWorld = new G4LogicalVolume(solidWorld, mat_Air, "World", 0, 0, 0);
    physicalWorld = new G4PVPlacement(0, G4ThreeVector(), "World", logicalWorld,
                                      NULL,   // mother volume
                                      false,  // no boolean operation
                                      0);     // copy number
    logicalWorld->SetVisAttributes(G4VisAttributes::GetInvisible);

    // shielding PE
    G4double BPE_xyz = 60.0 * cm;
    G4Box* BPE_box = new G4Box("BPE_box", 0.5 * BPE_xyz, 0.5 * BPE_xyz, 0.5 * BPE_xyz);
    // Plastic scintillator
    G4double PS_xyz = 50.0 * cm;
    G4Box* PS_box = new G4Box("PS_box", 0.5 * PS_xyz, 0.5 * PS_xyz, 0.5 * PS_xyz);

    // Beam hole
    G4double rMax = 5.0 * cm;
    G4Tubs* PS_hole = new G4Tubs("PS_hole", 0, rMax, 0.5 * BPE_xyz, 0, 360. * deg);

    // He3 Tube holes
    G4double active_len = 30 * cm;  // length
    G4double dead_len = 2 * cm;
    G4double tube_r = 2.54 / 2 * cm;  // radius
    G4double tube_t = 0.5 * mm;       // thickness

    G4double gas_region_len = active_len + 2.0 * dead_len;  // 34cm
    G4double gas_region_r = tube_r - tube_t;                // active region radius and dead region radius is the same
    G4double tube_len = gas_region_len + 2 * tube_t;        // 32.1cm

    G4double tube_hole_len = 0.5 * BPE_xyz + 0.5 * active_len + dead_len + 1 * cm;  // 48cm
    G4double tube_hole_r = 0.5 * 2.745 * cm;
    G4double tube_hole_z = -0.5 * BPE_xyz + 0.5 * tube_hole_len;
    G4Tubs* He3tube_Hole = new G4Tubs("He3tube_Hole", 0, tube_hole_r, 0.5 * tube_hole_len, 0, 360. * deg);

    G4Tubs* He3tub = new G4Tubs("He3tub", 0, tube_r, 0.5 * tube_len, 0, 360. * deg);
    G4Tubs* He3gas = new G4Tubs("He3gas", 0, gas_region_r, 0.5 * gas_region_len, 0, 360. * deg);
    G4Tubs* He3active = new G4Tubs("He3active", 0, gas_region_r, 0.5 * active_len, 0, 360. * deg);
    G4SubtractionSolid* He3tub_sub = new G4SubtractionSolid("He3tub_sub", He3tub, He3gas);
    G4RotationMatrix rot0;
    G4SubtractionSolid* He3gas_sub =
        new G4SubtractionSolid("He3gas_sub", He3gas, He3active, G4Transform3D(rot0, G4ThreeVector(0, 0, 0)));

    G4LogicalVolume* He3tub_log = new G4LogicalVolume(He3tub_sub, mat_stell, "He3tub_log");
    G4LogicalVolume* He3gas_log = new G4LogicalVolume(He3gas_sub, work_gas, "He3gas_log");
    G4LogicalVolume* He3active_log = new G4LogicalVolume(He3active, work_gas, "He3active_log");

    G4VisAttributes* He3active_vis = new G4VisAttributes(G4Colour(1.0, 0.0, 0.0));
    He3active_vis->SetForceSolid(true);
    He3active_log->SetVisAttributes(He3active_vis);
    He3tub_log->SetVisAttributes(G4VisAttributes::GetInvisible);

    G4SubtractionSolid* solid_PS = new G4SubtractionSolid("solid_PS", PS_box, PS_hole);
    G4SubtractionSolid* sub_BPE = new G4SubtractionSolid("sub_BPE", BPE_box, PS_box);
    G4SubtractionSolid* solid_BPE = new G4SubtractionSolid("solid_BPE", sub_BPE, PS_hole);

    char tmp[50];
    G4double posx, posy;
    for (G4int i = 0; i < 12; i++) {
        posx = innerR * std::cos(i * 30.0 * deg);
        posy = innerR * std::sin(i * 30.0 * deg);
        sprintf(tmp, "solid_PS%d", i);
        solid_PS = new G4SubtractionSolid(tmp, solid_PS, He3tube_Hole,
                                          G4Transform3D(rot0, G4ThreeVector(posx, posy, tube_hole_z)));
        sprintf(tmp, "solid_BPE%d", i);
        solid_BPE = new G4SubtractionSolid(tmp, solid_BPE, He3tube_Hole,
                                           G4Transform3D(rot0, G4ThreeVector(posx, posy, tube_hole_z)));
        posx = outerR * std::cos(i * 30.0 * deg + 15.0 * deg);
        posy = outerR * std::sin(i * 30.0 * deg + 15.0 * deg);
        sprintf(tmp, "solid_PS%d", i + 12);
        solid_PS = new G4SubtractionSolid(tmp, solid_PS, He3tube_Hole,
                                          G4Transform3D(rot0, G4ThreeVector(posx, posy, tube_hole_z)));
        sprintf(tmp, "solid_BPE%d", i + 12);
        solid_BPE = new G4SubtractionSolid(tmp, solid_BPE, He3tube_Hole,
                                           G4Transform3D(rot0, G4ThreeVector(posx, posy, tube_hole_z)));
    }

    G4LogicalVolume* PS_log = new G4LogicalVolume(solid_PS, ej_200, "logic_PS");
    new G4PVPlacement(0, G4ThreeVector(0, 0, 0), "PS_phys", PS_log, physicalWorld, false, 1, true);

    G4LogicalVolume* BPE_log = new G4LogicalVolume(solid_BPE, shield_PE, "logic_BPE");
    new G4PVPlacement(0, G4ThreeVector(0, 0, 0), "BPE_phys", BPE_log, physicalWorld, false, 2, true);

    G4VisAttributes* PS_vis = new G4VisAttributes(G4Colour(0.0, 0.8, 0.8, 0.5));
    PS_vis->SetForceSolid(true);
    // PS_vis->SetForceWireframe(true);
    PS_log->SetVisAttributes(PS_vis);

    G4VisAttributes* BPE_vis = new G4VisAttributes(G4Colour(0.1, 0.3, 1, 0.1));
    BPE_log->SetVisAttributes(BPE_vis);
    // BPE_log->SetVisAttributes(G4VisAttributes::GetInvisible);

    G4double posz = tube_hole_len - 0.5 * BPE_xyz - 0.5 * tube_len;
    for (G4int i = 0; i < 12; i++) {
        posx = innerR * std::cos(i * 30.0 * deg);
        posy = innerR * std::sin(i * 30.0 * deg);
        sprintf(tmp, "He3tub_phys%d", i);
        new G4PVPlacement(0, G4ThreeVector(posx, posy, posz + 0), tmp, He3tub_log, physicalWorld, false, 100 + i, true);
        sprintf(tmp, "He3gas_phys%d", i);
        new G4PVPlacement(0, G4ThreeVector(posx, posy, posz + 0), tmp, He3gas_log, physicalWorld, false, 200 + i, true);
        sprintf(tmp, "He3active_phys%d", i);
        new G4PVPlacement(0, G4ThreeVector(posx, posy, posz + 0 + 0), tmp, He3active_log, physicalWorld, false, 300 + i,
                          true);
        posx = outerR * std::cos(i * 30.0 * deg + 15.0 * deg);
        posy = outerR * std::sin(i * 30.0 * deg + 15.0 * deg);
        sprintf(tmp, "He3tub_phys%d", i + 12);
        new G4PVPlacement(0, G4ThreeVector(posx, posy, posz + 0), tmp, He3tub_log, physicalWorld, false, 100 + i + 12,
                          true);
        sprintf(tmp, "He3gas_phys%d", i + 12);
        new G4PVPlacement(0, G4ThreeVector(posx, posy, posz + 0), tmp, He3gas_log, physicalWorld, false, 200 + i + 12,
                          true);
        sprintf(tmp, "He3active_phys%d", i + 12);
        new G4PVPlacement(0, G4ThreeVector(posx, posy, posz + 0 + 0), tmp, He3active_log, physicalWorld, false,
                          300 + i + 12, true);
    }

    G4VisAttributes* tube_vis = new G4VisAttributes(G4Colour(0.8, 0.8, 0.8, 1.0));
    tube_vis->SetForceSolid(true);
    G4VisAttributes* Ne22_vis = new G4VisAttributes(G4Colour(1.0, 1.0, 1.0, 1.0));
    Ne22_vis->SetForceSolid(true);

    // Target Ne22
    G4double rBeam = rMax - 5 * mm;
    G4double beam_len = 60.0 * cm;
    G4double target_len = beam_len - 1 * cm;
    G4Tubs* Ne22tub = new G4Tubs("Ne22tub", 0, rBeam, 0.5 * target_len, 0, 360. * deg);
    G4LogicalVolume* Ne22log = new G4LogicalVolume(Ne22tub, target_Ne, "Ne22log");
    new G4PVPlacement(0, G4ThreeVector(), "Ne22phys", Ne22log, physicalWorld, false, 3, true);
    Ne22log->SetVisAttributes(Ne22_vis);

    // Beam Line Tube Stainlesss 不锈钢管
    G4Tubs* beamtub = new G4Tubs("Beam_tub", 0, rMax, 0.5 * beam_len, 0. * deg, 360. * deg);
    G4SubtractionSolid* beamtub_sub = new G4SubtractionSolid("Beam_sub", beamtub, Ne22tub);
    G4LogicalVolume* beam_log = new G4LogicalVolume(beamtub_sub, mat_stell, "beam_log");
    new G4PVPlacement(0, G4ThreeVector(0, 0, 0), "beam_phys", beam_log, physicalWorld, false, 10, true);
    beam_log->SetVisAttributes(tube_vis);

    G4SDManager* SDman = G4SDManager::GetSDMpointer();
    GeneHe3detSD* detSD = new GeneHe3detSD("/GeneHe3det/");
    SDman->AddNewDetector(detSD);
    He3active_log->SetSensitiveDetector(detSD);

    return physicalWorld;
}