#include "DetectorConstruction.hh"

#include "G4RunManager.hh"
#include "G4NistManager.hh"
#include "G4Box.hh"
#include "G4Cons.hh"
#include "G4Orb.hh"
#include "G4Sphere.hh"
#include "G4Trd.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4PVReplica.hh"
#include "G4SubtractionSolid.hh"
#include "G4SystemOfUnits.hh"
#include "G4ElementVector.hh"
#include "G4Material.hh"
#include "G4Isotope.hh"
#include "G4Element.hh"
#include "G4UniformMagField.hh"
#include "G4Mag_UsualEqRhs.hh"
#include "G4ClassicalRK4.hh"
#include "G4ChordFinder.hh"
#include "G4FieldManager.hh"
#include "G4TransportationManager.hh"
#include "G4Tubs.hh"
#include "CalorimeterSD.hh"
#include "G4SDManager.hh"
#include "TFile.h"
#include "G4VisAttributes.hh"
namespace Lohengrin 
{


G4VPhysicalVolume* DetectorConstruction::Construct()
{
  
  G4NistManager* nist = G4NistManager::Instance();

// here the world volume is created as a box with dimensions world_x, world_y and world_z
  G4double world_x = 5.0*m;//1.*m; // the world volume extends in all directions from -1m to +1m
  G4double world_y = 5.0*m;//1.*m;
  G4double world_z = 9.0*m;//1.*m;
  G4bool checkOverlaps = true; // all volumes may overlapp with the world volume

// here the material of the world volume, in this case 'vacuum', is defined with actomic number, mass in g/mole, etc.
  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;
  G4Material *Vacuum = new G4Material("vacuum", atomicNumber, massOfMole, density, kStateGas, temperature, pressure ); // a new G4Material is defined


  G4Box* worldSolid = new G4Box("World", world_x, world_y, world_z ); // create the world box with the dimensions defined above
  G4LogicalVolume* worldLogical = new G4LogicalVolume(worldSolid, Vacuum, "World" ); // create the logical volume filled with the material Vacuum and called World
  G4PVPlacement* worldPhysical = new G4PVPlacement( nullptr, G4ThreeVector(), worldLogical, "World", nullptr, false, 0, checkOverlaps ); // here we place the world box with a null pointer, this will be out relative coordinate system (center in center of world box)
 // G4PVPlacement::G4PVPlacement(G4RotationMatrix *	pRot, const G4ThreeVector &	tlate, G4LogicalVolume * pCurrentLogical, const G4String & pName, G4LogicalVolume * pMotherLogical, G4bool pMany, G4int	pCopyNo, G4bool pSurfChk = false)	 



// here the target is created
  G4double radlen_W = 0.3504*cm;
  G4double target_x = 13.*cm;//50.*cm; // dimensions of the target (again actually going from -50cm to +50cm)
  G4double target_y = 13.*cm;//50.*cm;
  G4double target_z = 0.03504*cm;//0.14*mm;////0.1*radlen_W;//50.*cm;


  // Define tungsten material of only W-184 
  G4int Z_W = 74; // atomic number of tungsten
  G4int N_W = 184; // mass number of tungsten-184
  G4double atomicmass_W = 183.84*g/mole; // atomic mass of W-184
  G4Isotope* W184 = new G4Isotope("W184", Z_W, N_W, atomicmass_W); // defining the new isotope
  G4Element* elW184 = new G4Element("Tungsten184", "W184", 1); // define element consisting of only that isotope
  elW184->AddIsotope(W184, 100.0*perCent); 

  G4int Z_Pb = 82; // atomic number 
  G4int N_Pb = 208; // mass number
  G4double atomicmass_Pb = 207.976652*g/mole; // atomic mass 
  G4Isotope* Pb208 = new G4Isotope("Pb208", Z_Pb, N_Pb, atomicmass_Pb); 
  G4Element* elPb208 = new G4Element("Lead208", "Pb208", 1); 
  elPb208->AddIsotope(Pb208, 100.0*perCent); 

  G4double density_Pb = 19.3*g/cm3;
  G4Material* lead208 = new G4Material("Lead208", density_Pb, 1); 
  lead208->AddElement(elPb208, 1);

  G4double density_W = 11.3*g/cm3;
  G4Material* tungsten184 = new G4Material("Tungsten184", density_W, 1); // define material of only W184
  tungsten184->AddElement(elW184, 1);

  // Target volume is created
  G4Box* targetSolid = new G4Box("Target", target_x/2, target_y/2, target_z/2 ); // target box is created with respective name and dimensions
  G4LogicalVolume *targetLogical = new G4LogicalVolume( targetSolid, tungsten184, "Target"); // create a logical target volume
  G4VPhysicalVolume *targetPhysical = new G4PVPlacement( 0, G4ThreeVector(0., 0., 0.*cm), targetLogical, "Target", worldLogical, false, 0 ); // creates the physical target volume and places it in the world box -> it is placed 1cm off the origin in z- direction

  

// Let's build a pixel detector 
  // G4double zPos1 = 1.0*cm; // z positions of the different pixel detectors behind the target
  // G4double zPos2 = 3.0*cm;
  // G4double zPos3 = 4.5*cm;
  // G4double zPos4 = 7.0*cm;
  G4Material* air = nist->FindOrBuildMaterial("G4_AIR");
  G4Material* silicon = nist->FindOrBuildMaterial("G4_Si");

  // G4double pixelSizeXY = 0.003304*cm;
  // G4double pixelSizeZ = 0.01*cm; 
  // G4Box* pixelBox = new G4Box("pixelBox", pixelSizeXY/2, pixelSizeXY/2, pixelSizeZ/2);
  // G4LogicalVolume* pixelLogical = new G4LogicalVolume(pixelBox, silicon, "PixelLogical");

  // G4double numPixelsX = 1024;
  // G4double numPixelsY = 1024;
  // G4double detectorSizeXY = numPixelsX*pixelSizeXY*cm;
  // G4double detectorSizeZ = pixelSizeZ*cm;
  // G4Box* pixelDetectorBox = new G4Box("PixelDetectorBox", detectorSizeXY/2, detectorSizeXY/2, detectorSizeZ/2);
  // G4LogicalVolume* pixelDetectorLogical1 = new G4LogicalVolume(pixelDetectorBox, air, "pixelDetectorLogical1");
  // G4LogicalVolume* pixelDetectorLogical2 = new G4LogicalVolume(pixelDetectorBox, air, "pixelDetectorLogical2");
  // G4LogicalVolume* pixelDetectorLogical3 = new G4LogicalVolume(pixelDetectorBox, air, "pixelDetectorLogical3");
  // G4LogicalVolume* pixelDetectorLogical4 = new G4LogicalVolume(pixelDetectorBox, air, "pixelDetectorLogical4");

  // // // Place pixel detector volumes in world volume
  // G4PVPlacement* pixelDetectorPhysical1 = new G4PVPlacement(0, G4ThreeVector(0,0,zPos1), pixelDetectorLogical1, "pixelDetectorPhysical1", worldLogical, false, 0);
  // G4PVPlacement* pixelDetectorPhysical2 = new G4PVPlacement(0, G4ThreeVector(0,0,zPos2), pixelDetectorLogical2, "pixelDetectorPhysical2", worldLogical, false, 0);
  // G4PVPlacement* pixelDetectorPhysical3 = new G4PVPlacement(0, G4ThreeVector(0,0,zPos3), pixelDetectorLogical3, "pixelDetectorPhysical3", worldLogical, false, 0);
  // G4PVPlacement* pixelDetectorPhysical4 = new G4PVPlacement(0, G4ThreeVector(0,0,zPos4), pixelDetectorLogical4, "pixelDetectorPhysical4", worldLogical, false, 0);


  // for(G4int i = 0; i < numPixelsX; ++i){ // placement of the individual pixels
  //   for(G4int j = 0; j < numPixelsY; ++j){
  //     G4double xPos = (i - numPixelsX/2 + pixelSizeXY/2)*pixelSizeXY;
  //     G4double yPos = (j - numPixelsY/2 + pixelSizeXY/2)*pixelSizeXY;
  //     // The pixels are placed in the coordinate system of the pixelDetectorPhysical and thus z=0
  //     new G4PVPlacement(0, G4ThreeVector(xPos, yPos, 0), pixelLogical, "PixelPhysical1", pixelDetectorLogical1, false, i*numPixelsY+j);
  //     new G4PVPlacement(0, G4ThreeVector(xPos, yPos, 0), pixelLogical, "PixelPhysical2", pixelDetectorLogical2, false, i*numPixelsY+j);
  //     new G4PVPlacement(0, G4ThreeVector(xPos, yPos, 0), pixelLogical, "PixelPhysical3", pixelDetectorLogical3, false, i*numPixelsY+j);
  //     new G4PVPlacement(0, G4ThreeVector(xPos, yPos, 0), pixelLogical, "PixelPhysical4", pixelDetectorLogical4, false, i*numPixelsY+j);
  //     }
  // }

 
// Magnetic Field here we go
  G4ThreeVector magFieldValue = G4ThreeVector(0, 0.9*tesla, 0); // defines the field strength; 
  //important: to have deflection in the y direction (which I assume we want) we need the magnetic field along the x-axis (cross porduct of v,B gives F)
  G4UniformMagField* magField = new G4UniformMagField(magFieldValue); // created magnetic field
  G4FieldManager* fieldManager = G4TransportationManager::GetTransportationManager()->GetFieldManager(); // create a field manager 
  fieldManager->SetDetectorField(magField); // and set field here
  G4Mag_UsualEqRhs* equation = new G4Mag_UsualEqRhs(magField); // create EoM
  G4ChordFinder* chordFinder = new G4ChordFinder(magField, 1.0e-2*mm, new G4ClassicalRK4(equation)); // create chord finder
  fieldManager->SetChordFinder(chordFinder);
  // to propagate a particle through a magn field its EoM is solved and then the curved path is divided into chord segments to come as close to EoM as possible 
  // each chord segment is checked if it crosses a material/volume boundary


// Okkk let's go with the magnet
  G4double innerR = 50*cm;
  G4double outerR = 90*cm;
  G4double height =  1.0*m; // magnet will be in total this 'long'
  G4double startAngle = 2.86*deg;
  G4double spanningAngle = (360-2.86)*deg;
  G4Tubs* solidMagnet = new G4Tubs("magnet", innerR, outerR, height/2, startAngle, spanningAngle);
  //G4NistManager* nist = G4NistManager::Instance();
  G4Material* iron = nist->FindOrBuildMaterial("G4_Fe");
  G4LogicalVolume* magnetLogical= new G4LogicalVolume(solidMagnet, iron, "LogicalMagnet");
  G4ThreeVector magnPosition = G4ThreeVector(0,0,height/2);
  new G4PVPlacement(0, magnPosition, magnetLogical, "magnet", worldLogical, false, 0, true);
  magnetLogical->SetFieldManager(fieldManager, true);
  


 // Time for the calorimeters
 // Calorimeter parameters:
  G4double thickSi = 0.5*cm;
  G4double thickAb = 2*cm;
  G4double numlayers = 30;
  G4double thickecal = numlayers*(thickSi+thickAb); 
  G4double lengthX = 48*cm;
  G4double lengthY = 48*cm;
  G4double slitecalX = 48*cm/2;
  G4double thickAbHCal = 2*cm;
  G4double numlayersHCal = 145; 
  G4double thickhcal = numlayersHCal*(thickSi+thickAbHCal); 
  G4double lengthXhcal = 270*cm;
  G4double lengthYhcal = 250*cm;

  G4double slitX = lengthXhcal/2; // needs to be half of the width of the hcal, as it only covers one side
  G4double slitY = 5*cm;

  G4Material* absorberHCal = nist->FindOrBuildMaterial("G4_Fe");
  G4Material* absorber = nist->FindOrBuildMaterial("G4_W");

  G4ThreeVector hcalPos = G4ThreeVector(0., 0., 350*cm+thickhcal/2);
  G4ThreeVector ecalPos = G4ThreeVector(0., 0., 350*cm+thickecal/2);//387.5*cm);
  G4ThreeVector slitPosHCal(lengthXhcal/4, 0., 0.);
  G4ThreeVector holePosition(0., 0., -thickhcal/2+thickecal/2);
  G4ThreeVector holePositionLayer(0., 0., 0.); // this position is relative to the Absorber box itself and hole is centered in it
  G4ThreeVector slitPosECal(lengthX/4, 0., 0.);


  // HCal detector Box
  G4Box* hcalDetectorBox = new G4Box("HcalDetectorBox", lengthXhcal/2, lengthYhcal/2, thickhcal/2);
  G4Box* ecalHoleBox = new G4Box("EcalHoleBox", lengthX/2, lengthY/2, thickecal/2);
  G4Box* slitBox = new G4Box("SlitBox", slitX/2, slitY/2, thickhcal/2);

  G4VSolid* hcalDetectorBoxWithHole = new G4SubtractionSolid("hcalDetectorBoxWithHole", hcalDetectorBox, ecalHoleBox, 0, holePosition);
  G4VSolid* hcalDetectorBoxWithHoles = new G4SubtractionSolid("hcalDetectorBoxWithHoles", hcalDetectorBoxWithHole, slitBox, 0, slitPosHCal);

  G4LogicalVolume* hcalDetectorLogical = new G4LogicalVolume(hcalDetectorBoxWithHoles, air, "hcalDetectorLogical");
  G4PVPlacement* hcalDetectorBoxPhysical = new G4PVPlacement(nullptr, hcalPos, hcalDetectorLogical, "hcalDetectorPhysical", worldLogical, false, 0, false);

//   Layers inside HCal
  G4Box* layers = new G4Box("Layer", lengthXhcal/2, lengthYhcal/2, thickSi+thickAbHCal);
  G4LogicalVolume* layerLogical = new G4LogicalVolume(layers, air, "Layer");
  
  G4double holeX = 48*cm;
  G4double holeY = 48*cm;

//   Absorber per layer  
  G4Box* AbsorberBoxHCal = new G4Box("AbsorberBoxHCal", lengthXhcal/2, lengthYhcal/2 , thickAbHCal/2);
  G4Box* HoleBoxHCalAb = new G4Box("HoleBoxAbsorber", holeX/2, holeY/2, thickAbHCal/2);
  G4Box* SlitBoxHCalAb = new G4Box("Hole2BoxAbsorber", slitX/2, slitY/2, thickAbHCal/2);
  G4VSolid* AbsorberBoxHCalWithHole = new G4SubtractionSolid("AbsorberBoxHCalWithHole", AbsorberBoxHCal, HoleBoxHCalAb, 0, holePositionLayer);
  G4VSolid* AbsorberBoxHCalWithHoles = new G4SubtractionSolid("AbsorberBoxHCalWithHoles", AbsorberBoxHCalWithHole, SlitBoxHCalAb, 0, slitPosHCal);
  G4LogicalVolume* AbsorberLogicalHCal = new G4LogicalVolume(AbsorberBoxHCalWithHoles, absorberHCal, "AbsorberLogicalHCal");
  new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -thickAbHCal/2), AbsorberLogicalHCal, "AbsorberHCalPhysical", layerLogical, false, 0, false); 
//   Si per layer
  G4Box* SiBoxHCal = new G4Box("SiBoxHCal", lengthXhcal/2, lengthYhcal/2 , thickSi/2);
  G4Box* HoleBoxHCalSi = new G4Box("HoleBoxSi", holeX/2, holeY/2, thickSi/2);
  G4Box* SlitBoxHCalSi = new G4Box("Hole2BoxSi", slitX/2, slitY/2, thickSi/2);
  G4VSolid* SiBoxHCalWithHole = new G4SubtractionSolid("SiBoxHCalWithHole", SiBoxHCal, HoleBoxHCalSi, 0, holePositionLayer);
  G4VSolid* SiBoxHCalWithHoles = new G4SubtractionSolid("SiBoxHCalWithHoles", SiBoxHCalWithHole, SlitBoxHCalSi, 0, slitPosHCal);
  G4LogicalVolume* SiLogicalHCal = new G4LogicalVolume(SiBoxHCalWithHoles, silicon, "SiLogicalHCal");
  new G4PVPlacement(nullptr, G4ThreeVector(0., 0., thickSi/2), SiLogicalHCal, "SiHCalPhysical", layerLogical, false, 0, false); 

  new G4PVReplica("Layer", layerLogical, hcalDetectorLogical, kZAxis, numlayersHCal, thickSi+thickAbHCal);

  // ECal detector box
  G4Box* ecalDetectorBox = new G4Box("EcalDetectorBox", lengthX/2, lengthY/2, thickecal/2); // defined this above for creating hole in hcal
  // G4Box* slitEcalBox = new G4Box("SlitBox", slitecalX/2, slitY/2, thickecal/2);
  // G4VSolid* ecalDetectorBoxWithHole = new G4SubtractionSolid("ecalDetectorBoxWithHole", ecalDetectorBox, slitEcalBox, 0, slitPosECal);
  G4LogicalVolume* ecalDetectorLogical = new G4LogicalVolume(ecalDetectorBox, air, "ecalDetectorLogical");
  G4PVPlacement* ecalDetectorBoxPhysical = new G4PVPlacement(nullptr, ecalPos, ecalDetectorLogical, "ecalDetectorPhysical", worldLogical, false, 0, true);

  // Layers
  G4Box* layersECal = new G4Box("LayerECal", lengthX/2, lengthY/2, thickSi+thickAb);
  G4LogicalVolume* layerECalLogical = new G4LogicalVolume(layersECal, air, "LayerECal");
  
 
  // Absorber per layer
  G4Box* AbsorberBox = new G4Box("AbsorberBox", lengthX/2, lengthY/2 , thickAb/2);
  // G4Box* HoleEcalBoxAbsorber = new G4Box("Hole2BoxAbsorber", slitecalX/2, slitY/2, thickAb/2);
  // G4VSolid* AbsorberBoxECalWithHole = new G4SubtractionSolid("AbsorberBoxECalWithHole", AbsorberBox, HoleEcalBoxAbsorber, 0, slitPosECal);
  G4LogicalVolume* AbsorberLogical = new G4LogicalVolume(AbsorberBox, absorber, "AbsorberLogical");
  G4PVPlacement* AbsorberPhysical = new G4PVPlacement(nullptr, G4ThreeVector(0., 0., -thickAb/2), AbsorberLogical, "AbsorberPhysical", layerECalLogical, false, 0, false);

  // Si per layer
  G4Box* SiBox = new G4Box("SiBox", lengthX/2, lengthY/2 , thickSi/2);
  // G4Box* HoleEcalBoxSi = new G4Box("Hole2BoxSi", slitecalX/2, slitY/2, thickSi/2);
  // G4VSolid* SiBoxECalWithHole = new G4SubtractionSolid("SiBoxECalWithHole", SiBox, HoleEcalBoxSi, 0, slitPosECal);
  G4LogicalVolume* SiLogical = new G4LogicalVolume(SiBox, silicon, "SiLogical");
  G4PVPlacement* SiPhysical = new G4PVPlacement(nullptr, G4ThreeVector(0., 0., thickSi/2), SiLogical, "SiPhysical", layerECalLogical, false, 0, false);

  new G4PVReplica("LayerECal", layerECalLogical, ecalDetectorLogical, kZAxis, numlayers, thickSi+thickAb);

  // G4VisAttributes* worldVisAtt = new G4VisAttributes(G4Colour(1.0, 1.0, 1.0)); // white for world
  // worldVisAtt->SetVisibility(true);
  // worldLogical->SetVisAttributes(worldVisAtt);

  // // Set visualization attributes for HCal
  // G4VisAttributes* HCalVisAtt = new G4VisAttributes(G4Colour(0., 0., 1.)); // Blue for HCal Ab
  // HCalVisAtt->SetVisibility(true);
  // hcalDetectorLogical->SetVisAttributes(HCalVisAtt);

  // G4VisAttributes* EcalVisAtt = new G4VisAttributes(G4Colour(0., 1., 0.)); // Green for ECal Si
  // EcalVisAtt->SetVisibility(true);
  // ecalDetectorLogical->SetVisAttributes(EcalVisAtt);

  // // Set visualization attributes for target
  // G4VisAttributes* TargetVisAtt = new G4VisAttributes(G4Colour(1.0, 0.0, 1.)); // magenta for HCal Ab
  // TargetVisAtt->SetVisibility(true);
  // targetLogical->SetVisAttributes(TargetVisAtt);

  // G4VisAttributes* MagnetVisAtt = new G4VisAttributes(G4Colour(1.0, 0., 0.)); // red for HCal Ab
  // MagnetVisAtt->SetVisibility(true);
  // magnetLogical->SetVisAttributes(MagnetVisAtt);




  std::cout << "Checking overlaps..." << std::endl;
  ecalDetectorBoxPhysical->CheckOverlaps();
  AbsorberPhysical->CheckOverlaps();
  SiPhysical->CheckOverlaps();

  G4SDManager* sdManager = G4SDManager::GetSDMpointer();

  CalorimeterSD* ecalSD = new CalorimeterSD("ecalSD");
  sdManager->AddNewDetector(ecalSD);
  SiLogical->SetSensitiveDetector(ecalSD);

  CalorimeterSD* hcalSD = new CalorimeterSD("hcalSD");
  sdManager->AddNewDetector(hcalSD);
  SiLogicalHCal->SetSensitiveDetector(hcalSD);


  // TFile *file = TFile::Open("geometry.root", "RECREATE");
  // file->WriteObject(worldLogical);
  // file->Close();

  return worldPhysical;
}


}