#include "DetectorConstruction.hh"

#include "DetectorSensitiveDetector.hh"

#include "G4RunManager.hh"
#include "G4SDManager.hh"
#include "G4NistManager.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4SystemOfUnits.hh"
#include "G4RotationMatrix.hh"
#include "G4ThreeVector.hh"
#include "G4VisAttributes.hh"

DetectorConstruction::DetectorConstruction()
: G4VUserDetectorConstruction(),
  fScoringVolume(nullptr),
  fAnodeLogical(nullptr),
  fCathodeLogical(nullptr),
  fSX3Logical(nullptr),
  fQQQLogical(nullptr)
{}

DetectorConstruction::~DetectorConstruction()
{}

G4VPhysicalVolume* DetectorConstruction::Construct()
{
  // Get nist material manager
  G4NistManager* nist = G4NistManager::Instance();

  // Materials
  G4Material* vacuum = nist->FindOrBuildMaterial("G4_Galactic");
  G4Material* al = nist->FindOrBuildMaterial("G4_Al");

  // World
  G4double worldx = 200.*mm;
  G4double worldy = 200.*mm;
  G4double worldz = 200.*mm;
  G4Box* solidWorld = new G4Box("World", worldx, worldy, worldz);
  G4LogicalVolume* logicWorld = new G4LogicalVolume(solidWorld, vacuum, "World");
  G4VPhysicalVolume* physWorld = new G4PVPlacement(0, G4ThreeVector(), logicWorld, "World", 0, false, 0, true);

  // Helium target

  G4double pressure = 379.0 * 133.322368 * pascal;
  G4double temperature = 293.15 * kelvin;
  G4double density = 8.29e-5 * g/cm3;

  G4Material* HeGas =
      new G4Material("HeGas",
                    density,
                    1,
                    kStateGas,
                    temperature,
                    pressure);

  HeGas->AddElement(nist->FindOrBuildElement("He"), 1);
  const G4double targetRadius = 37.0*mm;
  const G4double targetHalfLength = 175.*mm;
  G4Tubs* solidTarget = new G4Tubs("Target", 0, targetRadius, targetHalfLength, 0, 360*deg);
  G4LogicalVolume* logicTarget =
    new G4LogicalVolume(solidTarget, HeGas, "Target");
  new G4PVPlacement(0, G4ThreeVector(), logicTarget, "Target", logicWorld, false, 0, true);
  logicTarget->SetVisAttributes(G4VisAttributes::GetInvisible());
  fScoringVolume = logicTarget;

  // Axes (optional, for visualization)
  G4Tubs* solidAxis = new G4Tubs("Axis", 0, 0.1*mm, 5.*mm, 0, 360*deg);
  G4LogicalVolume* logicAxisX = new G4LogicalVolume(solidAxis, al, "AxisX");
  G4LogicalVolume* logicAxisY = new G4LogicalVolume(solidAxis, al, "AxisY");
  G4LogicalVolume* logicAxisZ = new G4LogicalVolume(solidAxis, al, "AxisZ");

  G4RotationMatrix* rotX = new G4RotationMatrix();
  rotX->rotateY(90*deg);
  rotX->rotateZ(90*deg);
  new G4PVPlacement(rotX, G4ThreeVector(5*mm, 0, 0), logicAxisX, "AxisX", logicWorld, false, 0);

  G4RotationMatrix* rotY = new G4RotationMatrix();
  rotY->rotateX(90*deg);
  new G4PVPlacement(rotY, G4ThreeVector(0, 5*mm, 0), logicAxisY, "AxisY", logicWorld, false, 0);

  new G4PVPlacement(0, G4ThreeVector(0, 0, 5*mm), logicAxisZ, "AxisZ", logicWorld, false, 0);

  // ANASEN geometry
  const int nWire = 24;
  const int wireShift = 4;
  const G4double zLen = 350.*mm;
  const G4double radiusA = 38.*mm;
  const G4double radiusC = 43.*mm;

  G4double dAngle = wireShift * 2 * CLHEP::pi / nWire;
  G4double radiusAnew = radiusA * cos(dAngle / 2.);
  G4double wireALength = sqrt(zLen*zLen + pow(2 * radiusA * sin(dAngle/2.), 2));
  G4double wireATheta = atan2(2 * radiusA * sin(dAngle / 2.), zLen);

  G4Tubs* solidPC_A = new G4Tubs("PC_A", 0, .1*mm, wireALength/2., 0, 360*deg);
  G4LogicalVolume* logicPC_A = new G4LogicalVolume(solidPC_A, al, "PC_A");
  fAnodeLogical = logicPC_A;

  for(int i = 0; i < nWire; i++){
    G4double phi = 2 * CLHEP::pi / nWire * i + dAngle / 2.;
    G4ThreeVector pos(radiusAnew * cos(phi), radiusAnew * sin(phi), 0);
    // Euler angles: phi (Z), theta (Y), psi (Z)
    G4double phi_euler = 360./nWire * (i + wireShift/2.) * deg;
    G4double theta_euler = wireATheta;
    G4double psi_euler = 0;
    G4RotationMatrix* rot = new G4RotationMatrix(phi_euler, theta_euler, psi_euler);
    new G4PVPlacement(rot, pos, logicPC_A, "PC_A", logicWorld, false, i);
  }

  G4double radiusCnew = radiusC * cos(dAngle / 2.);
  G4double wireCLength = sqrt(zLen*zLen + pow(2 * radiusC * sin(dAngle/2.), 2));
  G4double wireCTheta = atan2(2 * radiusC * sin(dAngle / 2.), zLen);

  G4Tubs* solidPC_C = new G4Tubs("PC_C", 0, .1*mm, wireCLength/2., 0, 360*deg);
  G4LogicalVolume* logicPC_C = new G4LogicalVolume(solidPC_C, al, "PC_C");
  fCathodeLogical = logicPC_C;

  for(int i = 0; i < nWire; i++){
    G4double phi = 2 * CLHEP::pi / nWire * i - dAngle/2.;
    G4ThreeVector pos(radiusCnew * cos(phi), radiusCnew * sin(phi), 0);
    // Euler angles: phi (Z), theta (Y), psi (Z)
    G4double phi_euler = 360./nWire * (i - wireShift/2.) * deg;
    G4double theta_euler = -wireCTheta;
    G4double psi_euler = 0;
    G4RotationMatrix* rot = new G4RotationMatrix(phi_euler, theta_euler, psi_euler);
    new G4PVPlacement(rot, pos, logicPC_C, "PC_C", logicWorld, false, i);
  }

  const int nSX3 = 12;
  const G4double sx3Radius = 88.*mm;
  const G4double sx3Width = 40.*mm;
  const G4double sx3Length = 75.*mm;
  const G4double sx3Gap = 5.*mm;

  G4Box* solidSX3 = new G4Box("SX3", 0.1*mm, sx3Width/2., sx3Length/2.);
  G4LogicalVolume* logicSX3 = new G4LogicalVolume(solidSX3, al, "SX3");
  fSX3Logical = logicSX3;
  fScoringVolume = logicSX3;

  for(int i = 0; i < nSX3; i++){
    G4double phi = 2 * CLHEP::pi / nSX3 * (i + 0.5);
    G4ThreeVector pos1(sx3Radius * cos(phi), sx3Radius * sin(phi), sx3Length/2. + sx3Gap);
    // Euler angles: phi (Z), theta (Y), psi (Z)
    G4double phi_euler = 360./nSX3 * (i + 0.5) * deg;
    G4RotationMatrix* rot1 = new G4RotationMatrix(phi_euler, 0, 0);
    new G4PVPlacement(rot1, pos1, logicSX3, "SX3_front", logicWorld, false, 2*i);

    G4ThreeVector pos2(sx3Radius * cos(phi), sx3Radius * sin(phi), -sx3Length/2. - sx3Gap);
    G4RotationMatrix* rot2 = new G4RotationMatrix(phi_euler, 0, 0);
    new G4PVPlacement(rot2, pos2, logicSX3, "SX3_back", logicWorld, false, 2*i+1);
  }

  const G4double qqqR1 = 50.*mm;
  const G4double qqqR2 = 100.*mm;
  G4Tubs* solidQQQ = new G4Tubs("QQQ", qqqR1, qqqR2, 0.5*mm, 5*deg, 85*deg);
  G4LogicalVolume* logicQQQ = new G4LogicalVolume(solidQQQ, al, "QQQ");
  fQQQLogical = logicQQQ;

  for(int i = 0; i < 4; i++){
    G4ThreeVector pos(0, 0, 100.*mm);
    // Euler angles: phi (Z), theta (Y), psi (Z)
    G4double phi_euler = 360./4 * i * deg;
    G4RotationMatrix* rot = new G4RotationMatrix(phi_euler, 0, 0);
    new G4PVPlacement(rot, pos, logicQQQ, "QQQ", logicWorld, false, i);
  }

  return physWorld;
}

void DetectorConstruction::ConstructSDandField()
{
  auto* detectorSD = new DetectorSensitiveDetector("AnasenSensitiveDetector");
  G4SDManager::GetSDMpointer()->AddNewDetector(detectorSD);

  if (fSX3Logical != nullptr) {
    SetSensitiveDetector(fSX3Logical, detectorSD);
  }
  if (fQQQLogical != nullptr) {
    SetSensitiveDetector(fQQQLogical, detectorSD);
  }
}