#include <iostream>
#include <cstdlib>
#define THICKNESS 0.001		// define the thickness of air and aluminum in each forward calorimeter

#include "LGDetectorConstruction.hh"

#include "G4Element.hh"
#include "G4ElementTable.hh"
#include "G4Material.hh"
#include "G4MaterialTable.hh"
#include "G4NistManager.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4Trd.hh"
#include "G4Cons.hh"
#include "G4Torus.hh"
#include "G4LogicalVolume.hh"
#include "G4ThreeVector.hh"
#include "G4RotationMatrix.hh"
#include "G4Transform3D.hh"
#include "G4PVPlacement.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4OpBoundaryProcess.hh"
#include "G4SubtractionSolid.hh"
#include "G4ios.hh"

LGDetectorConstruction::LGDetectorConstruction(G4double zLG, G4double Lround, 
					       G4double BendR, G4double BendA, G4int BendFirst): 
  zTapperedSection(zLG), zLGRoundSection(Lround), LGBendRadius(BendR), LGBendAngle(BendA), FirstDoBend(BendFirst)
{

//--- 1. LightGuide: You must be carefule because all the length means half of the actual length.

  x_LGrec = 2.54*cm/2.;   //thicknes light guide rectangular cross section (paddle thickness)
  y_LGrec = 6.0*cm/2.;     //hight light guide rectangular cross section (paddle width)
  z_LGrec = 3.0*cm/2.; //length or rectangular part

  x_Scin  = x_LGrec;   // part of scintillator to be used to generate light distribution at light guide entrance
  y_Scin  = y_LGrec;
  z_Scin  = 30.*cm;

  r_LGcyl = 2.5*cm;                    //radius of strait cylindrical part of light guide
  z_LGcyl = zLGRoundSection*cm/2.;     //length of strait cylindrical part of light guide

  x1_LGtrap =  x_LGrec;
  y1_LGtrap =  y_LGrec;
  x2_LGtrap =  r_LGcyl;
  y2_LGtrap =  x2_LGtrap;
  z_LGtrap  = zTapperedSection*cm/2.;

  r1m_LGcone = sqrt(x1_LGtrap*x1_LGtrap + y1_LGtrap*y1_LGtrap);
  r1M_LGcone = y_LGrec+2.*cm;
  r2m_LGcone = r_LGcyl;
  r2M_LGcone = r_LGcyl+2.*cm;
  z_LGcone   = z_LGtrap;

  //			-------- 4. World 
  x_World = 200*cm;
  y_World = 200*cm;
  z_World = 400*cm;		
  
}

LGDetectorConstruction::~LGDetectorConstruction() {;}

G4VPhysicalVolume* LGDetectorConstruction::Construct()
{
  //			----------------------------------------------------------------------
  //			-------- 1. Elements and Materials 
  //			----------------------------------------------------------------------
  G4double a, z, density;			// a: mass of a mole, z: mean number of protons
  G4int n_Elements;			// n_Elements: the number of kinds of the elements
  G4int ncomponents, natoms;
  G4String name, symbol; 
 
  G4NistManager* nist_Manager = G4NistManager::Instance();
  
  //			-------- 1.2. Materials for various parts of the detector
  G4Material* air		= nist_Manager->FindOrBuildMaterial("G4_AIR");	  // Air
  G4Material* Al          = nist_Manager->FindOrBuildMaterial("G4_Al");     // Aluminum
  
  G4Material* PlexiGlass	= nist_Manager->FindOrBuildMaterial("G4_PLEXIGLASS");	// Plexiglass
  
  //G4Material* pmtGlass	= nist_Manager->FindOrBuildMaterial("G4_SILICON_DIOXIDE");// PMTWindow
  // for pmtGlass see lower

  // define Elements
  G4Element *elH, *elC;
  a = 1.01*g/mole;
  elH  = new G4Element(name="Hydrogen",symbol="H" , z= 1., a);
  
  a = 12.01*g/mole;
  elC  = new G4Element(name="Carbon"  ,symbol="C" , z= 6., a);

  //G4Material* Hydrogen    = nist_Manager->FindOrBuildMaterial("G4_H");	
  //G4Material* Carbon      = nist_Manager->FindOrBuildMaterial("G4_C");	
  G4Material* Oxigen      = nist_Manager->FindOrBuildMaterial("G4_O");	
  G4Material* Silicon     = nist_Manager->FindOrBuildMaterial("G4_Si");	
  G4Material* Potassium   = nist_Manager->FindOrBuildMaterial("G4_K");	
  G4Material* Sodium      = nist_Manager->FindOrBuildMaterial("G4_Na");	
  //G4Material* Lead        = nist_Manager->FindOrBuildMaterial("G4_Pb");	
  G4Material* Boron       = nist_Manager->FindOrBuildMaterial("G4_B");	
  G4Material* Aluminum    = nist_Manager->FindOrBuildMaterial("G4_Al");	
  
  // data for ELJEN 200 scintillator
  density = 1.023*g/cm3;
  G4Material* Scintillator = new G4Material(name="Scintillator", density, ncomponents=2);
  Scintillator->AddElement(elC, natoms=469);
  Scintillator->AddElement(elH, natoms=517);
   
  //Fraction by weight for Borosicilcat PMT window glass 
  //source http://en.wikipedia.org/wiki/Borosilicate_glass
  //Element    Atomicnumber 	Fraction
  //----------------------------------------
  //B 		5 		0.040064
  //O 		8 		0.539562
  //Na 	 	11 	 	0.028191
  //Al 	 	13 	 	0.011644
  //Si 	 	14 	 	0.377220
  //K 		19	 	0.003321
  G4Material* pmtGlass = new G4Material("Borosicilate",density=2.23*g/cm3,n_Elements=6);
  pmtGlass->AddMaterial(Oxigen,     0.539562); 
  pmtGlass->AddMaterial(Boron,      0.040064);
  pmtGlass->AddMaterial(Sodium,     0.028191);
  pmtGlass->AddMaterial(Aluminum,   0.011644);
  pmtGlass->AddMaterial(Silicon,    0.377220);
  pmtGlass->AddMaterial(Potassium,  0.003321);
  
  //			----------------------------------------------------------------------
  //			-------- 2. Generate and Add Material Property Tables
  //			----------------------------------------------------------------------
  
  const G4int num = 4;
  
  G4double eRange_Photon[num]		= {0.5*eV,1.24*eV, 2.27*eV, 3.76*eV};	
  
  G4double absLength_PlexiGlass[num]	= {500.*cm,500.*cm, 500.*cm, 500.*cm};	
  G4double absLength_PMTGlass[num]	= {200.*cm,200.*cm, 200.*cm, 200.*cm};	
  G4double index_Air[num]               = {1.,1., 1., 1.};	       // Index of refraction of air

  G4double index_Scintillator[num]               = {1.58,1.58, 1.58, 1.58};	      
  G4double absLength_Scintillator[num]  = {400.*cm,400.*cm, 400.*cm, 400.*cm};
	 
  // source http://en.wikipedia.org/wiki/Borosilicate_glass
  //G4double index_PMTGlass[num]		= {1.52,1.52, 1.52, 1.52};
  G4double index_PMTGlass[num]		= {1.48,1.48, 1.48, 1.48}; // in original code
  
  // Index of refraction of the PlexiGlass light guide
  G4double index_PlexiGlass[num]		= {1.51,1.51, 1.51, 1.51};  
  
  G4double eff_Cathode[num]		= {1.,1., 1., 1.};	       // eff = efficiency	
  G4double ref_Cathode[num]		= {0.,0., 0., 0.};	       // ref = reflectivity
  
  //			-------- 2.1. Table - Air MPT= Material Property Table
  G4MaterialPropertiesTable* air_MPT = new G4MaterialPropertiesTable();
  air_MPT->AddProperty("RINDEX", eRange_Photon, index_Air, num);
  air->SetMaterialPropertiesTable(air_MPT);
  
  G4MaterialPropertiesTable* Scintillator_MPT = new G4MaterialPropertiesTable();
  Scintillator_MPT->AddProperty("RINDEX", eRange_Photon, index_Scintillator, num);
  Scintillator_MPT->AddProperty("ABSLENGTH", eRange_Photon, absLength_Scintillator, num);
  Scintillator->SetMaterialPropertiesTable(Scintillator_MPT);

  //			-------- 2.4. Table - plexiglass
  G4MaterialPropertiesTable* PlexiGlass_MPT = new G4MaterialPropertiesTable();
  PlexiGlass_MPT->AddProperty("RINDEX", eRange_Photon, index_PlexiGlass, num);
  PlexiGlass_MPT->AddProperty("ABSLENGTH", eRange_Photon, absLength_PlexiGlass, num);
  PlexiGlass->SetMaterialPropertiesTable(PlexiGlass_MPT);
  
  //			-------- 2.5. Table - PMTglass
  G4MaterialPropertiesTable* pmtGlass_MPT = new G4MaterialPropertiesTable();
  pmtGlass_MPT->AddProperty("RINDEX", eRange_Photon, index_PMTGlass, num);
  pmtGlass_MPT->AddProperty("ABSLENGTH", eRange_Photon, absLength_PMTGlass, num);
  pmtGlass->SetMaterialPropertiesTable(pmtGlass_MPT);
  
  //			-------- 2.6. Table - PMTcathode
  G4MaterialPropertiesTable* cathode_MPT = new G4MaterialPropertiesTable();
  cathode_MPT->AddProperty("EFFICIENCY", eRange_Photon, eff_Cathode, num);
  cathode_MPT->AddProperty("REFLECTIVITY", eRange_Photon, ref_Cathode, num);
  Al->SetMaterialPropertiesTable(cathode_MPT);
  
  
  //			----------------------------------------------------------------------
  //			-------- 3. Volumes
  //			----------------------------------------------------------------------
  //			-------- 3.1. World Volume
  //			----------------------------------------------------------------------
  
  // Definiton for the world volume
  G4Box* world_box			= new G4Box("World", x_World, y_World, z_World);
  
  // Definiton of the logical volume of the world volume filled with the air
  G4LogicalVolume* world_log		= new G4LogicalVolume(world_box, air, "World");
  
  // Placement of the world volume
  G4VPhysicalVolume* world_phys		= new G4PVPlacement(0, G4ThreeVector(), world_log, "World", 0, false, 0);
  
  
  //			----------------------------------------------------------------------
  //			-------- 3.2. Scintillator with light guide
  //			----------------------------------------------------------------------

  G4Box* LGBox_box			= new G4Box("LGBox_box", x_LGrec, y_LGrec, z_LGrec);
  G4LogicalVolume* LGBox_log		= new G4LogicalVolume(LGBox_box, PlexiGlass, "LGBox_log");
  G4VPhysicalVolume* LGBox_phys = NULL;

  G4double zpos = 0;
  LGBox_phys = new G4PVPlacement(0, G4ThreeVector(0., 0., zpos), 
				 LGBox_log, "LGBox", world_log, false, 0);

  G4Box* ScinBox_box                    = new G4Box("ScinBox_box", x_Scin, y_Scin, z_Scin);
  G4LogicalVolume* ScinBox_log          = new G4LogicalVolume(ScinBox_box, PlexiGlass,"ScinBox_log");
  G4VPhysicalVolume* ScinBox_phys = NULL;
  ScinBox_phys = new  G4PVPlacement(0, G4ThreeVector(0., 0., -z_Scin-z_LGrec), 
				    ScinBox_log, "ScinBox", world_log, false, 0);

  G4Trd* LGTrd_trd		 = new G4Trd("LGTrd", x1_LGtrap, x2_LGtrap, y1_LGtrap, y2_LGtrap, z_LGtrap);
  
  G4Cons* LGCon_con              = new G4Cons("LGCon",  r1m_LGcone, r1M_LGcone, r2m_LGcone, r2M_LGcone, z_LGcone, 
					      0.*deg, 360.*deg);

  G4VSolid* TapperedSection          = new G4SubtractionSolid("TapperedSection", LGTrd_trd, LGCon_con);
  G4LogicalVolume* LGTappered_log    = new G4LogicalVolume(TapperedSection, PlexiGlass, "LGTappered_log"); 
  
  G4double LGBendAng = LGBendAngle*deg;
  G4double LGBendRad = LGBendRadius*cm;

  G4VPhysicalVolume* LGTappered_phys    = NULL;
  G4LogicalVolume* LGBendSection_log    = NULL;
  G4VPhysicalVolume* LGBendSection_phys = NULL;

  G4RotationMatrix* RotX = new G4RotationMatrix();
  RotX->rotateX(-90.*deg);
  G4RotationMatrix* RotY = new G4RotationMatrix();
  RotY->rotateY(LGBendAng);
  G4double transz1 = 0;
  G4double transz2 = 0;
  G4double transx1 = 0;
  G4double transx2 = 0;
  G4double dx = 0;
  G4double dz = 0;
  
  G4Torus* LGBendSection_tor = NULL;
  G4Tubs*  LGBendSection_tub = NULL;

  if (!FirstDoBend){ // first set tapered section then add bending section as torus
    zpos += (z_LGrec + z_LGcone); 
    LGTappered_phys = new G4PVPlacement(0, G4ThreeVector(0.,0.,zpos), LGTappered_log,
					"LGTappered", world_log, false, 0);
    
    LGBendSection_tor          = new G4Torus("LGBendSection", 0., r_LGcyl, LGBendRad, 0., LGBendAng);
    LGBendSection_log  = new G4LogicalVolume(LGBendSection_tor, PlexiGlass, "LGBendSection_log");
    
    zpos += z_LGtrap; 
    LGBendSection_phys= new G4PVPlacement(RotX, G4ThreeVector(-LGBendRad, 0., zpos),
					  LGBendSection_log, "LGBendSection", world_log, false, 0);

    G4cout<<"Length of bent section is " << (LGBendRad)*LGBendAng/10.<<"cm"<<G4endl;

    transz1 = LGBendRad*sin(LGBendAng);
    //transz2 = z_LGcyl*sin(M_PI/2.-LGBendAng);
    transx1 = LGBendRad - LGBendRad*cos(LGBendAng);
    //transx2 = z_LGcyl*cos(M_PI/2.-LGBendAng);
    dx = transx1 ;//+ transx2;
    dz = transz1 ;//+ transz2;

  } else { // first add bended square section and then strait tappered section
    zpos += z_LGrec;
    LGBendSection_tub              = new G4Tubs("LGBendSection", LGBendRad, LGBendRad+x_LGrec*2.,
						    y_LGrec, 0., LGBendAng);
    LGBendSection_log = new G4LogicalVolume(LGBendSection_tub, PlexiGlass, "LGBendSection_log");
    LGBendSection_phys = new G4PVPlacement(RotX, G4ThreeVector( -LGBendRad-x_LGrec, 0., zpos), 
							       LGBendSection_log, "LGBendSection", world_log, false, 0);

    G4cout<<"Length of bent section is " << (LGBendRad+x_LGrec)*LGBendAng/10.<<"cm"<<G4endl;

    transz1 = (LGBendRad+x_LGrec)*sin(LGBendAng);
    transz2 = z_LGcone*sin(M_PI/2.-LGBendAng);
    transx1 = LGBendRad + x_LGrec - (LGBendRad+x_LGrec)*cos(LGBendAng);
    transx2 = z_LGcone*cos(M_PI/2.-LGBendAng);
    dx = transx1 + transx2;
    dz = transz1 + transz2;
    zpos += dz;
    LGTappered_phys = new G4PVPlacement(RotY, G4ThreeVector(-dx, 0., zpos), LGTappered_log,
					"LGTappered", world_log, false, 0);
    dx += transx2;
    dz = transz2;
  }
  
  G4Tubs* LGTub_tubs		 = new G4Tubs("LGTub_tubs", 0., r_LGcyl, z_LGcyl, 
					      0.*deg, 360.*deg);
  G4LogicalVolume* LGTub_log	 = new G4LogicalVolume(LGTub_tubs, PlexiGlass, "LGTub_log");

  //dx += transx2;
  //dz += transz2;
  transz2 = z_LGcyl*sin(M_PI/2.-LGBendAng);
  transx2 = z_LGcyl*cos(M_PI/2.-LGBendAng);
  dx += transx2;
  dz += transz2;
  
  zpos  += dz;
  G4VPhysicalVolume* LGTub_phys	 = new G4PVPlacement(RotY, G4ThreeVector(-dx, 0., zpos), 
						     LGTub_log, "LGTub", world_log, false, 0);
  
  
  // pmtglass 
  G4double PMTGlassThickness    = 1.0*mm;
  G4Tubs* pmtGlass_tubs		= new G4Tubs("PMTGlass_tubs", 0., r_LGcyl, 
					     PMTGlassThickness, 0.*deg, 360.*deg); 
  G4LogicalVolume* pmtGlass_log	= new G4LogicalVolume(pmtGlass_tubs, pmtGlass, "PMTGlass_log");
  
  dx += transx2;
  transx1 = PMTGlassThickness * cos(M_PI/2.-LGBendAng);
  dx += transx1;
  zpos += transz2;
  transz1 = PMTGlassThickness * sin(M_PI/2.-LGBendAng);
  zpos += transz1;
  G4VPhysicalVolume* pmtGlass_phys = NULL;
  pmtGlass_phys                    = new  G4PVPlacement(RotY, G4ThreeVector(-dx, 0, zpos), 
							pmtGlass_log, "PMTGlass", world_log, false, 0);
  
  // pmt cathode
  G4double PMTCathodeThickness    = 0.001*mm;
  G4Tubs* pmtCathode_tubs         = new G4Tubs("PMTCathode_tubs", 0., r_LGcyl*0.99, 
					       PMTCathodeThickness, 0.*deg, 360.*deg);
  G4LogicalVolume* pmtCathode_log = new G4LogicalVolume(pmtCathode_tubs, Aluminum, "PMTCathode_log");
  
  dx += transx1;
  transx2 = PMTCathodeThickness * cos(M_PI/2.-LGBendAng);
  dx += transx2;
  zpos += transz1;
  transz2 = PMTCathodeThickness * sin(M_PI/2.-LGBendAng);
  zpos += transz2;
  
  G4VPhysicalVolume* pmtCathode_phys = NULL;
  pmtCathode_phys = new 
    G4PVPlacement(RotY, G4ThreeVector(-dx, 0, zpos), 
		  pmtCathode_log, "PMTCathode", world_log, false, 0);
  
  
  
  
  //			----------------------------------------------------------------------
  //			-------- 4. Surfaces
  //			----------------------------------------------------------------------
  //			-------- 4.1. Generate and Add Surface Property Tables
  //			----------------------------------------------------------------------
  
  //G4double ref_PlexiGlass[num]	= {0.95,0.95, 0.95, 0.95}; // to be used for Alu foil air interface		
  //G4double eff_PlexiGlass[num]	= {0.05,0.05, 0.05, 0.05};		
     
  G4double specular_Lobe[num]	= {0.045,0.045, 0.045, 0.045};
  G4double specular_Spike[num]	= {0.95,.95, 0.95, 0.95};
  G4double back_Scatter[num]	= {0.005,0.005, 0.005, 0.005};	
  
  // surface properties from Chris Cude simulation
  //G4double specular_Lobe[num]	= {0.18,0.18, 0.18, 0.18};
  //G4double specular_Spike[num]	= {0.7954,.7954, 0.7954, 0.7954};
  //G4double back_Scatter[num]	= {0.0246,0.0246, 0.0246, 0.0246};	
  
  // specular_lobe + specular_spike + back_scatter + defuse_lobe = 1
  
  
  //			----------------------------------------------------------------------
  //			-------- 4.2. Surfaces
  //			----------------------------------------------------------------------
  //			-------- 4.2.1. Surface of the leadglass in contact with air
  G4OpticalSurface* PlexiGlassAir_Surf = new G4OpticalSurface("PlexiGlassAirSurface");
  PlexiGlassAir_Surf->SetType(dielectric_dielectric);
  //PlexiGlassAir_Surf->SetFinish(polished);
  PlexiGlassAir_Surf->SetFinish(ground);
  PlexiGlassAir_Surf->SetModel(unified);
  
  // becaus of option "ground" use a sigma for the angle alpha needs 
  // to be set. A sigma of 10 would indicate a very rough surface.
  PlexiGlassAir_Surf->SetSigmaAlpha(0.1);
	
  //			----------------------------------------------------------------------
  //			-------- 4.3. Property Tables
  //			----------------------------------------------------------------------

  // plexiglass to air interface
  
  G4LogicalBorderSurface* plexiGlassAir_Surf_log1 = NULL;
  G4LogicalBorderSurface* plexiGlassAir_Surf_log2 = NULL;
  G4LogicalBorderSurface* plexiGlassAir_Surf_log3 = NULL;
  G4LogicalBorderSurface* plexiGlassAir_Surf_log4 = NULL;
 
  plexiGlassAir_Surf_log1        = new 
    G4LogicalBorderSurface("PlexiGlassAirSurface1",  LGBox_phys, world_phys, PlexiGlassAir_Surf);
  plexiGlassAir_Surf_log2        = new 
    G4LogicalBorderSurface("PlexiGlassAirSurface2",  LGTappered_phys, world_phys, PlexiGlassAir_Surf);
  plexiGlassAir_Surf_log3        = new 
    G4LogicalBorderSurface("PlexiGlassAirSurface3",  LGBendSection_phys, world_phys, PlexiGlassAir_Surf);
  plexiGlassAir_Surf_log4        = new 
    G4LogicalBorderSurface("PlexiGlassAirSurface4",  LGTub_phys, world_phys, PlexiGlassAir_Surf);


  //			-------- 4.1.4. LeadGlass(LeadGlass_Air)
  G4MaterialPropertiesTable* PlexiGlass_Air_MPT = new G4MaterialPropertiesTable();
  PlexiGlass_Air_MPT->AddProperty("SPECULARLOBECONSTANT", eRange_Photon, specular_Lobe, num);
  PlexiGlass_Air_MPT->AddProperty("SPECULARSPIKECONSTANT", eRange_Photon, specular_Spike, num);
  PlexiGlass_Air_MPT->AddProperty("BACKSCATTERCONSTANT", eRange_Photon, back_Scatter, num);
  PlexiGlassAir_Surf->SetMaterialPropertiesTable(PlexiGlass_Air_MPT);

  return world_phys;  // return the physical volume of the world volume
}