#include "DetectorConstruction.hh"

#include "G4NistManager.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4RotationMatrix.hh"
#include "G4Transform3D.hh"
#include "G4SDManager.hh"
#include "G4MultiFunctionalDetector.hh"
#include "G4VPrimitiveScorer.hh"
#include "G4PSEnergyDeposit.hh"
#include "G4VisAttributes.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"

#include <sstream>

#include "G4UniformElectricField.hh"
#include "G4UniformMagField.hh"
#include "G4MagneticField.hh"
#include "G4FieldManager.hh"
#include "G4TransportationManager.hh"
#include "G4EquationOfMotion.hh"
#include "G4EqMagElectricField.hh"
#include "G4Mag_UsualEqRhs.hh"
#include "G4MagIntegratorStepper.hh"
#include "G4MagIntegratorDriver.hh"
#include "G4ChordFinder.hh"

#include "G4ExplicitEuler.hh"
#include "G4ImplicitEuler.hh"
#include "G4SimpleRunge.hh"
#include "G4SimpleHeum.hh"
#include "G4ClassicalRK4.hh"
#include "G4HelixExplicitEuler.hh"
#include "G4HelixImplicitEuler.hh"
#include "G4HelixSimpleRunge.hh"
#include "G4CashKarpRKF45.hh"
#include "G4RKG3_Stepper.hh"

#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"

using namespace std;

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------


DetectorConstruction::DetectorConstruction()
: G4VUserDetectorConstruction(),
  fElectriclogical_1(0),
  fElectriclogical_2(0),
  fEMfield_1(0),
  fEquation_1(0),
  fStepper_1(0),
  fFieldMgr_1(0),
  fMinStep_1(0),
  fChordFinder_1(0),
  fEMfield_2(0),
  fEquation_2(0),
  fStepper_2(0),
  fFieldMgr_2(0),
  fMinStep_2(0),
  fChordFinder_2(0),
  fEMfield_WORLD(0),
  fEquation_WORLD(0),
  fStepper_WORLD(0),
  fFieldMgr_WORLD(0),
  fMinStep_WORLD(0),
  fChordFinder_WORLD(0),
  fCheckOverlaps(true)
{
  //Constructor
}

//-----------------------------------------------------------------------------

DetectorConstruction::~DetectorConstruction()
{
  //Destructor
}


//-----------------------------------------------------------------------------

G4VPhysicalVolume* DetectorConstruction::Construct()
{
  //World dimensions
  G4double worldsizeX = 5. * m;
  G4double worldsizeY = 5. * m;
  G4double worldsizeZ = 5. * m;

  //Box dimensions
  G4double box_X = 300.* cm;
  G4double box_Y = 20.* cm;
  G4double box_Z = 5.* cm;

  auto nistManager = G4NistManager::Instance();
  //-----------------------------------------------------------
  G4double atomicNumber = 1.;
  G4double massOfMole = 1.008*g/mole;
  G4double density = 1.e-25*g/cm3;
  G4double temperature = 2.73*kelvin;
  G4double pressure = 3.e-18*pascal;
  auto Vacuum = new G4Material("interGalactic", atomicNumber,massOfMole, density, kStateGas,temperature, pressure);
  //-----------------------------------------------------------
  auto BoxMaterial = nistManager->FindOrBuildMaterial("G4_STAINLESS-STEEL");    //Box

  //Let's create the "World"
  G4Box* WorldBox = new G4Box("World",worldsizeX / 2, worldsizeY / 2, worldsizeZ / 2);
  G4LogicalVolume* logicWorld = new G4LogicalVolume(WorldBox,Vacuum,"World");
  G4VPhysicalVolume* physWorld = new G4PVPlacement(0,G4ThreeVector(),logicWorld,"World",0,false,0,fCheckOverlaps);                           // checking overlaps

   //Box_1 inside the world
   G4Box* Box_1 = new G4Box("Box_1",box_X / 2, box_Y / 2, box_Z / 2);

   //This logical volume has a local electric field. Check ConstructSDandField function below.
   fElectriclogical_1 = new G4LogicalVolume(Box_1,BoxMaterial,"BoxLV_1");
   new G4PVPlacement(0,G4ThreeVector(0.*cm, 0.*cm, -5.*cm),fElectriclogical_1,"BoxLV_1",logicWorld,false,0,fCheckOverlaps);

   //Box_2 inside the world
   G4Box* Box_2 = new G4Box("Box_2",box_X / 2, box_Y / 2, box_Z / 2);

   //This logical volume has a local electric field. Check ConstructSDandField function below.
   fElectriclogical_2 = new G4LogicalVolume(Box_2,BoxMaterial,"BoxLV_2");
   new G4PVPlacement(0,G4ThreeVector(0.*cm, 0.*cm, 5.*cm),fElectriclogical_2,"BoxLV_2",logicWorld,false,0,fCheckOverlaps);

    //world
    //logicWorld->SetVisAttributes(G4VisAttributes::GetInvisible());

    //Box_1
    auto g_color = new G4VisAttributes(G4Colour::Blue());
    g_color->SetVisibility(true);
    //g_color->SetForceSolid(true);
    fElectriclogical_1->SetVisAttributes(g_color);
    fElectriclogical_2->SetVisAttributes(g_color);

    //Let's print the materials Table
    G4cout << *(G4Material::GetMaterialTable()) << G4endl;

    //Always return the physical world
    return physWorld;
}

//-----------------------------------------------------------------------------


void DetectorConstruction::ConstructSDandField()
{
  G4SDManager::GetSDMpointer()->SetVerboseLevel(1);

  //---------------------------------------------------------
  //This is the first local electric field(Box_1)
  fEMfield_1 = new G4UniformElectricField(G4ThreeVector(4500.0*kilovolt/m, 0.0, 0.0));
  // Create an equation of motion for this field
  fEquation_1 = new G4EqMagElectricField(fEMfield_1);

  G4int nvar1 = 8;
  fStepper_1 = new G4ClassicalRK4(fEquation_1,nvar1);

  auto localFieldMgr_1 = new G4FieldManager(fEMfield_1);
  fElectriclogical_1->SetFieldManager(localFieldMgr_1,true);

  fMinStep_1 = 0.00010*mm ; // minimal step of 10 microns

  auto fIntgrDriver_1 = new G4MagInt_Driver(fMinStep_1,fStepper_1,fStepper_1->GetNumberOfVariables());

  fChordFinder_1 = new G4ChordFinder(fIntgrDriver_1);
  localFieldMgr_1->SetChordFinder(fChordFinder_1);
  //----------------------------------------------------------
  //----------------------------------------------------------
  //----------------------------------------------------------
  //---------------------------------------------------------

  //This is the second local electric field(Box_2)
  fEMfield_2 = new G4UniformElectricField(G4ThreeVector(-4500.0*kilovolt/m, 0.0, 0.0));
  // Create an equation of motion for this field
  fEquation_2 = new G4EqMagElectricField(fEMfield_2);

  G4int nvar2 = 8;
  fStepper_2 = new G4ClassicalRK4(fEquation_2,nvar2);

  auto localFieldMgr_2 = new G4FieldManager(fEMfield_2);
  fElectriclogical_2->SetFieldManager(localFieldMgr_2,true);

  fMinStep_2 = 0.00010*mm ; // minimal step of 10 microns

  auto fIntgrDriver_2 = new G4MagInt_Driver(fMinStep_2,fStepper_2,fStepper_2->GetNumberOfVariables());

  fChordFinder_2 = new G4ChordFinder(fIntgrDriver_2);
  localFieldMgr_2->SetChordFinder(fChordFinder_2);
  //----------------------------------------------------------
  //----------------------------------------------------------
  //----------------------------------------------------------
  //----------------------------------------------------------

  //This is the world's electric field.Currently, it is set to 0
  fEMfield_WORLD= new G4UniformElectricField(G4ThreeVector(0.0, 0.0,0.0*kilovolt/m));

  // Create an equation of motion for this field
  fEquation_WORLD = new G4EqMagElectricField(fEMfield_WORLD);

  G4int nvar3 = 8;
  fStepper_WORLD= new G4ClassicalRK4(fEquation_WORLD,nvar3);

  // Get the global field manager
  auto fFieldManager= G4TransportationManager::GetTransportationManager()->GetFieldManager();

  // Set this field to the global field manager
  fFieldManager->SetDetectorField(fEMfield_WORLD);

  fMinStep_WORLD = 0.010*mm ; // minimal step of 10 microns

  auto fIntgrDriver_WORLD = new G4MagInt_Driver(fMinStep_WORLD,fStepper_WORLD,fStepper_WORLD->GetNumberOfVariables());

  fChordFinder_WORLD = new G4ChordFinder(fIntgrDriver_WORLD);
  fFieldManager->SetChordFinder(fChordFinder_WORLD);
  //----------------------------------------------------------
  //----------------------------------------------------------
  //----------------------------------------------------------
  //----------------------------------------------------------


  //We need to declare the box as sensitive
  //in order to get information about the energy
  G4MultiFunctionalDetector* boxDetector = new G4MultiFunctionalDetector("Box_1");
  G4SDManager::GetSDMpointer()->AddNewDetector(boxDetector);
  G4VPrimitiveScorer* primitive = new G4PSEnergyDeposit("Edep");
  boxDetector->RegisterPrimitive(primitive);
  SetSensitiveDetector("BoxLV_1",boxDetector);

}

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------