// Written by Peter Stromberger
// based on work of Mirco Deckenhoff
// modified by Bastian Rössler


#include "DetectorConstruction.hh"
#include "G4Material.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4LogicalSkinSurface.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4EllipticalTube.hh"
#include "G4SubtractionSolid.hh"
#include "G4LogicalVolume.hh"
#include "G4RotationMatrix.hh"
#include "G4PVPlacement.hh"
#include "G4OpBoundaryProcess.hh"
#include "G4MaterialPropertyVector.hh"
#include "G4NistManager.hh"
#include "G4VisAttributes.hh"
#include "G4Colour.hh"
#include "G4PVParameterised.hh"
#include "SensitiveDetector.hh"
#include "G4SDManager.hh"
#include "Analysis.hh"
#include "Parameters.hh"
#include "TF1.h"
#include "TF2.h"
#include "Randomize.hh"
#include <fstream>
#include <sstream>
#include "G4SystemOfUnits.hh"
#include "G4VPVParameterisation.hh"
#include "Convert.hh"
#include <math.h>
#include <stdlib.h> 
#include <time.h>  
#include <iostream>
#include <fstream>

DetectorConstruction::DetectorConstruction()
: G4VUserDetectorConstruction(),
  fScoringVolume(0)
{
    // Size of experimental hall and fibre length //
    expHall_x = 35*mm/2.+0.5*m;
    expHall_y = 10*mm/2.+0.5*m;
    expHall_z = (Parameters::GetInstance()->FibreLength()+0.5)/2.*m+0.5*m;
    scint_z = Parameters::GetInstance()->FibreLength()/2.*m;

    // mat/detector dimensions
    Nj = 8;
    Nk = 128;
  
    xDist = 0.350 * mm;
    yDist = 0.185 * mm;
    stripWidth = 0.25 * mm;
    stripHeight = 1.62 * mm;

    airGap = 0.001 * mm;
    pixelDimX = 0.0625 * mm;
    pixelDimY = 0.0675 * mm;
    epoxy_strip_width = 0.1*mm;
	
    Nx = 4;
    Ny = 24;

    
    fMirrorToggle = true;
    fMirrorPolish = 1.;
    fMirrorReflectivity = 0.8;

   
    fMirrorZ = 0.1*mm;
    fMirrorRmax = 0.250*mm;
}

DetectorConstruction::~DetectorConstruction(){}

void DetectorConstruction::DefineMaterials() 
{
    // Get nist material manager
    G4NistManager* man = G4NistManager::Instance();
    // Option to switch on/off checking of volumes overlaps
    //
    G4bool checkOverlaps = true;
    

    // Elements to cunstruct inner cladding material (PMMA)
    G4double densityPMMA = 1190*kg/m3;
    std::vector<G4String> PMMA_elm;
    std::vector<G4int> PMMA_nbAtoms;
    PMMA_elm.push_back("H"); PMMA_nbAtoms.push_back(8);
    PMMA_elm.push_back("C"); PMMA_nbAtoms.push_back(5);
    PMMA_elm.push_back("O"); PMMA_nbAtoms.push_back(2);
    PMMA = man->ConstructNewMaterial("PMMA", PMMA_elm, PMMA_nbAtoms, densityPMMA);
    
  	// Elements to cunstruct outer cladding material (PTFEMA)
    G4double densityPMMA2 = 1430*kg/m3;
    std::vector<G4String> PMMA2_elm;
    std::vector<G4int> PMMA2_nbAtoms;
    PMMA2_elm.push_back("H"); PMMA2_nbAtoms.push_back(7);
    PMMA2_elm.push_back("C"); PMMA2_nbAtoms.push_back(6);
    PMMA2_elm.push_back("O"); PMMA2_nbAtoms.push_back(2);
    PMMA2_elm.push_back("F"); PMMA2_nbAtoms.push_back(3);
    PMMA2 = man->ConstructNewMaterial("PMMA2", PMMA2_elm, PMMA2_nbAtoms, densityPMMA2=1430*kg/m3);
  
    // Glue (Epo-Tek 301)
    G4double Glue_density = 1.15*g/cm3;
    std::vector<G4String> Glue_elm;
    std::vector<G4int>	Glue_nbAtoms;
    Glue_elm.push_back("C"); Glue_nbAtoms.push_back(19);
    Glue_elm.push_back("H"); Glue_nbAtoms.push_back(27);
    Glue_elm.push_back("O"); Glue_nbAtoms.push_back(3);
    Glue = man->ConstructNewMaterial("Glue",Glue_elm, Glue_nbAtoms, Glue_density);
  
    // TiO2
    G4double TiO2_density = 4.26*g/cm3;
    std::vector<G4String> TiO2_elm;
    std::vector<G4int>	TiO2_nbAtoms;
    TiO2_elm.push_back("Ti"); TiO2_nbAtoms.push_back(1);
    TiO2_elm.push_back("Ti"); TiO2_nbAtoms.push_back(2);
    TiO2 = man->ConstructNewMaterial("TiO2",TiO2_elm, TiO2_nbAtoms, TiO2_density);
  
    // Abs plastic
    G4double Abs_density = 1.07*g/cm3;
    std::vector<G4String> Abs_elm;
    std::vector<G4int> Abs_nbAtoms;
    Abs_elm.push_back("H"); Abs_nbAtoms.push_back(17);
    Abs_elm.push_back("C"); Abs_nbAtoms.push_back(13);
    Abs_elm.push_back("N"); Abs_nbAtoms.push_back(1);
    Abs_plastic = man->ConstructNewMaterial("Abs_plastic", Abs_elm, Abs_nbAtoms, Abs_density);

    // Epoxy
    G4double Epoxy_density = 1.5*g/cm3;
    Epoxy = new G4Material("Epoxy", Epoxy_density, 2);
    Epoxy->AddMaterial(TiO2, 25*perCent);
    Epoxy->AddMaterial(Glue, 75*perCent);

    // Environment
    Air     = man->FindOrBuildMaterial("G4_AIR");
    Vacuum  = man->FindOrBuildMaterial("G4_Galactic");

    // Polystyrene   G4_POLYSTYRENE
    Pstyrene = man->FindOrBuildMaterial("G4_POLYSTYRENE");

    // Aluminium
    alu = man->FindOrBuildMaterial("G4_Al");

    // Core sections
    G4String materialNameCore = "ScintCoreMaterial";
    scintCoreMaterial = new G4Material(materialNameCore,Pstyrene->GetDensity(),1);
    scintCoreMaterial->AddMaterial(Pstyrene,1.);

    // Inner cladding sections
    G4String materialNameCladding1 = "InnerCladdingMaterial";
    innerCladdingMaterial = new G4Material(materialNameCladding1,PMMA->GetDensity(),1);
    innerCladdingMaterial->AddMaterial(PMMA,1.);

    // Outer cladding sections
    G4String materialNameCladding2 = "OuterCladdingMaterial";
    outerCladdingMaterial = new G4Material(materialNameCladding2,PMMA2->GetDensity(),1);
    outerCladdingMaterial->AddMaterial(PMMA2,1.);
}




void DetectorConstruction::DefineMaterialProperties() 
{
    // Initialise considered ENERGIES and EMISSION SPECTRA for scintillation and wls

    const G4int numInterpolPoints = Parameters::GetInstance()->NumberOfInterpolatedPoints();

    // Set scintillation emmision spectrum
    G4int E_NUMENTRIES = Parameters::GetInstance()->NumberOfEnergies();
    G4double* Energy = new G4double[E_NUMENTRIES];
    for(int i=0; i<E_NUMENTRIES; i++)
        Energy[i] = Parameters::GetInstance()->Energy[i]*eV;

    G4double* ScintilEnergyDist = new G4double[E_NUMENTRIES];
    for(int i=0; i<E_NUMENTRIES; i++)
        ScintilEnergyDist[i] = Parameters::GetInstance()->Intensity[i];

    G4MaterialPropertyVector* scintSpecVector = new G4MaterialPropertyVector(Energy, ScintilEnergyDist, E_NUMENTRIES);
    scintSpecVector->SetSpline(true);

    // Set the Birks Constant for the Polystyrene scintillator
    Pstyrene->GetIonisation()->SetBirksConstant(Parameters::GetInstance()->BirksConstant()*mm/MeV);

    // Set wls emmision spectrum
    G4int WLS_E_NUMENTRIES = Parameters::GetInstance()->NumberOfWlsEmissionEnergies();
    G4double* WlsEnergy = new G4double[WLS_E_NUMENTRIES];
    for(int i=0; i<WLS_E_NUMENTRIES; i++)
        WlsEnergy[i] = Parameters::GetInstance()->WlsEmissionEnergy[i]*eV;

    G4double* WlsEnergyDist = new G4double[WLS_E_NUMENTRIES];
    for(int i=0; i<WLS_E_NUMENTRIES; i++)
        WlsEnergyDist[i] = Parameters::GetInstance()->WlsEmissionIntensity[i];

    G4MaterialPropertyVector* wlsSpecVector = new G4MaterialPropertyVector(WlsEnergy, WlsEnergyDist, WLS_E_NUMENTRIES);
  	wlsSpecVector->SetSpline(true);

    if(numInterpolPoints>0)
    {
        // Interpolate scintillation spectrum

        const G4int E_NUMENTRIES_New = (E_NUMENTRIES-1)*(numInterpolPoints+1)+1;
        G4double* NewEnergy = new G4double[E_NUMENTRIES_New];
        G4double* NewValue = new G4double[E_NUMENTRIES_New];
        G4double* scintInterpolValue = new G4double[numInterpolPoints];
    
        for(int j=0; j<E_NUMENTRIES-1; j++)
        {
            G4double interpolDist = (Energy[j+1]-Energy[j])/(numInterpolPoints+1);
            NewEnergy[j*(numInterpolPoints+1)] = Energy[j];
            NewValue[j*(numInterpolPoints+1)] = scintSpecVector->Value(Energy[j]);

            for(int k=0; k<numInterpolPoints; k++)
            {
                scintInterpolValue[k] = scintSpecVector->Value(Energy[j]+(k+1)*interpolDist);
                NewEnergy[j*(numInterpolPoints+1)+(k+1)] = Energy[j]+(k+1)*interpolDist;

                if(scintInterpolValue[k]>=0)
                {
                    NewValue[j*(numInterpolPoints+1)+(k+1)] = scintInterpolValue[k];
                }
                else
                {
                    NewValue[j*(numInterpolPoints+1)+(k+1)] = 0;
                    G4cout << "Warning: Intensity of emission spectrum set to 0 for energy "
                        << NewEnergy[j*(numInterpolPoints+1)+(k+1)]*1e6
                        << " eV.\nSpline interpolation led to negative value!"  << G4endl;
				}
            }
        }
        
        NewEnergy[E_NUMENTRIES_New-1] = Energy[E_NUMENTRIES-1];
        NewValue[E_NUMENTRIES_New-1] = scintSpecVector->Value(Energy[E_NUMENTRIES-1]);
    
        delete[] Energy;
        delete[] ScintilEnergyDist;
        delete scintSpecVector;

        E_NUMENTRIES = E_NUMENTRIES_New;
        Energy = new G4double[E_NUMENTRIES];
        ScintilEnergyDist = new G4double[E_NUMENTRIES];

        for(int i=0; i<E_NUMENTRIES; i++)
        {
            Energy[i] = NewEnergy[i];
            ScintilEnergyDist[i] = NewValue[i];
    	}

        scintSpecVector = new G4MaterialPropertyVector(Energy, ScintilEnergyDist, E_NUMENTRIES);
        scintSpecVector->SetSpline(true);
   
        delete[] NewEnergy;
        delete[] NewValue;
        delete[] scintInterpolValue;


        // Interpolate wls spectrum
    
        const G4int WLS_E_NUMENTRIES_New = (WLS_E_NUMENTRIES-1)*(numInterpolPoints+1)+1;
        NewEnergy = new G4double[WLS_E_NUMENTRIES_New];
        NewValue = new G4double[WLS_E_NUMENTRIES_New];

        G4double* wlsInterpolValue = new G4double[numInterpolPoints];
    
        for(int j=0; j<WLS_E_NUMENTRIES-1; j++)
        {
            G4double interpolDist = (WlsEnergy[j+1]-WlsEnergy[j])/(numInterpolPoints+1);
            NewEnergy[j*(numInterpolPoints+1)] = WlsEnergy[j];
            NewValue[j*(numInterpolPoints+1)] = wlsSpecVector->Value(WlsEnergy[j]);

            for(int k=0; k<numInterpolPoints; k++)
            {
                wlsInterpolValue[k] = wlsSpecVector->Value(WlsEnergy[j]+(k+1)*interpolDist);
                NewEnergy[j*(numInterpolPoints+1)+(k+1)] = WlsEnergy[j]+(k+1)*interpolDist;
                if(wlsInterpolValue[k]>=0)
                {
                    NewValue[j*(numInterpolPoints+1)+(k+1)] = wlsInterpolValue[k];
                }
                else
                {
                    NewValue[j*(numInterpolPoints+1)+(k+1)] = 0;
                    G4cout << "Warning: Intensity of WLS emission spectrum set to 0 for energy "
                           << NewEnergy[j*(numInterpolPoints+1)+(k+1)]*1e6
                           << " eV.\nSpline interpolation led to negative value!" << G4endl;
                }
            }
        }
        
        NewEnergy[WLS_E_NUMENTRIES_New-1] = WlsEnergy[WLS_E_NUMENTRIES-1];
        NewValue[WLS_E_NUMENTRIES_New-1] = wlsSpecVector->Value(WlsEnergy[WLS_E_NUMENTRIES-1]);
    
        delete[] WlsEnergy;
        delete[] WlsEnergyDist;
        delete wlsSpecVector;

        WLS_E_NUMENTRIES = WLS_E_NUMENTRIES_New;
        WlsEnergy = new G4double[WLS_E_NUMENTRIES];
        WlsEnergyDist = new G4double[WLS_E_NUMENTRIES];

        for(int i=0; i<WLS_E_NUMENTRIES; i++)
        {
            WlsEnergy[i] = NewEnergy[i];
            WlsEnergyDist[i] = NewValue[i];
        }

        wlsSpecVector = new G4MaterialPropertyVector(WlsEnergy, WlsEnergyDist, WLS_E_NUMENTRIES);
        wlsSpecVector->SetSpline(true);

        delete[] NewEnergy;
        delete[] NewValue;
        delete[] wlsInterpolValue;
    }

    // Save spectra to parameter file
    std::ofstream parameterOutputFile;
    parameterOutputFile.open(Parameters::GetInstance()->ParameterOutputFileName(), std::ios_base::app);
    
    parameterOutputFile << G4endl << "Scintillation emission spectrum:" << G4endl;
    parameterOutputFile << "Energy/eV\tWavelength/nm\tIntensity" << "\n";
    for(int i=0; i<E_NUMENTRIES; i++)
    {
        parameterOutputFile << Energy[i]*1e6 << "\t" << Parameters::hcPERe/Energy[i]*1e3
                            << "\t" << ScintilEnergyDist[i]<< "\n";
	}

    parameterOutputFile << G4endl << "WLS emission spectrum:" << G4endl;
    parameterOutputFile << "Energy/eV\tWavelength/nm\tIntensity" << "\n";
    for(int i=0; i<WLS_E_NUMENTRIES; i++)
    {
        parameterOutputFile << WlsEnergy[i]*1e6 << "\t" << Parameters::hcPERe/WlsEnergy[i]*1e3
                            << "\t" << WlsEnergyDist[i]<< "\n";
	}
    
    parameterOutputFile.close();

    // WLS ABSORPTION
    const G4int WLS_ABS_ENTRIES = Parameters::GetInstance()->NumberOfWlsAbsEnergies();
    G4double* WlsAbsEnergy = new G4double[WLS_ABS_ENTRIES];
    G4double* WlsAbsLength = new G4double[WLS_ABS_ENTRIES];

    for(int j=0; j<WLS_ABS_ENTRIES; j++)
    {
        WlsAbsEnergy[j] = Parameters::GetInstance()->WlsAbsEnergy[j]*eV;
        WlsAbsLength[j] = Parameters::GetInstance()->WlsAbsLength[j]*m;
    }

    // Save WLS absorption lengths to parameter file
    parameterOutputFile.open(Parameters::GetInstance()->ParameterOutputFileName(), std::ios_base::app);
    parameterOutputFile << G4endl << "WLS absorption length / m:" << G4endl;
    parameterOutputFile << "Energy/eV\tWavelength/nm\tAbsorption length" << "\n";
    for(int i=0; i<WLS_ABS_ENTRIES; i++)
    {
        parameterOutputFile << WlsAbsEnergy[i]*1e6 << "\t" << Parameters::hcPERe/WlsAbsEnergy[i]*1e3
                            << "\t" << WlsAbsLength[i]*1e-3 << "\n";
	}
    parameterOutputFile.close();

    // REFRACTIVE INDICES

    G4double* Vacuum_RIND = new G4double[E_NUMENTRIES];
    G4double* Pstyrene_RIND = new G4double[E_NUMENTRIES];
    G4double* PMMA_RIND = new G4double[E_NUMENTRIES];
    G4double* PMMA2_RIND = new G4double[E_NUMENTRIES];
    G4double* Epoxy_RIND = new G4double[E_NUMENTRIES];
	
    // Functions are saved in parameter-file
    // todo: write this functions in c++ code
    TF1 vacuumRind("vacuumRind",Parameters::GetInstance()->RefractiveIndexVacuum(),300,800);
    TF1 coreRind("coreRind",Parameters::GetInstance()->RefractiveIndexCore(),300,800);
    TF1 clad1Rind("clad1Rind",Parameters::GetInstance()->RefractiveIndexClad1(),300,800);
    TF1 clad2Rind("clad2Rind",Parameters::GetInstance()->RefractiveIndexClad2(),300,800);

    for(int i=0; i<E_NUMENTRIES; i++)
    {
        double wavelengthNanometer = Parameters::hcPERe/Energy[i]*1e3;
        Vacuum_RIND[i] = vacuumRind.Eval(wavelengthNanometer);
        Pstyrene_RIND[i] = coreRind.Eval(wavelengthNanometer);
        PMMA_RIND[i] = clad1Rind.Eval(wavelengthNanometer);
        PMMA2_RIND[i] = clad2Rind.Eval(wavelengthNanometer);
        Epoxy_RIND[i] = 1.59;
    }

    // Save refractive indices to parameter file

    parameterOutputFile.open(Parameters::GetInstance()->ParameterOutputFileName(), std::ios_base::app);
    parameterOutputFile << G4endl << "Refractive indices:" << G4endl;
    parameterOutputFile << "Energy/eV\tWavelength/nm\tVacuum\tClad2\tClad1\tCore" << "\n";
    for(int i=0; i<E_NUMENTRIES; i++)
    {
        parameterOutputFile << Energy[i]*1e6 << "\t" << Parameters::hcPERe/Energy[i]*1e3
                            << "\t" << Vacuum_RIND[i]<< "\t" << PMMA2_RIND[i] << "\t"
                            << PMMA_RIND[i] << "\t" << Pstyrene_RIND[i] << "\n";
    }
    parameterOutputFile.close();


    // Set material properties table of Vacuum

    G4double* Vacuum_ABS = new G4double[E_NUMENTRIES];
    G4double* Epoxy_ABS = new G4double[E_NUMENTRIES];
    for(int i=0; i<E_NUMENTRIES; i++){
        Vacuum_ABS[i] = 5e4*m; // absorption length in vacuum...
        Epoxy_ABS[i] = 1*m;
    }

    G4MaterialPropertiesTable *Vacuum_mt = new G4MaterialPropertiesTable();
    Vacuum_mt->AddProperty("RINDEX", Energy, Vacuum_RIND,E_NUMENTRIES);
    Vacuum_mt->AddProperty("ABSLENGTH",Energy,Vacuum_ABS,E_NUMENTRIES);
    Vacuum->SetMaterialPropertiesTable(Vacuum_mt);

    G4MaterialPropertiesTable *Epoxy_mt = new G4MaterialPropertiesTable();
    Epoxy_mt->AddProperty("RINDEX", Energy, Epoxy_RIND, E_NUMENTRIES);
    Epoxy_mt->AddProperty("ABSLENGTH", Energy, Epoxy_ABS, E_NUMENTRIES);
    Glue->SetMaterialPropertiesTable(Epoxy_mt);

    // give air same material properties as vacuum
    Air->SetMaterialPropertiesTable(Vacuum_mt);

    // RAYLEIGH SCATTERING

    G4double* Pstyrene_RAYLEIGH = new G4double[E_NUMENTRIES];
    G4double* PMMA_RAYLEIGH = new G4double[E_NUMENTRIES];
    G4double* PMMA2_RAYLEIGH = new G4double[E_NUMENTRIES];

    TF1 coreRayleigh("coreRayleigh",Parameters::GetInstance()->RayleighCore(),300,800);
    TF1 clad1Rayleigh("clad1Rayleigh",Parameters::GetInstance()->RayleighClad1(),300,800);
    TF1 clad2Rayleigh("clad2Rayleigh",Parameters::GetInstance()->RayleighClad2(),300,800);
  
    // Calculate scattering lengths
  
    for(int j=0; j<E_NUMENTRIES; j++)
    {
        double wavelengthNanometer = Parameters::hcPERe/Energy[j]*1e3;
        // todo : write this functions in c++
        Pstyrene_RAYLEIGH[j] = 1./coreRayleigh.Eval(wavelengthNanometer)*m;
        PMMA_RAYLEIGH[j] 		= 1./clad1Rayleigh.Eval(wavelengthNanometer)*m;
        PMMA2_RAYLEIGH[j] 	= 1./clad2Rayleigh.Eval(wavelengthNanometer)*m;
  	}

    // Save Rayleigh scattering lengths to parameter file

    parameterOutputFile.open(Parameters::GetInstance()->ParameterOutputFileName(), std::ios_base::app);
    parameterOutputFile << G4endl << "Rayleigh scattering length / m:" << G4endl;
    parameterOutputFile << "Energy/eV\tWavelength/nm\tClad2\tClad1\tCore" << "\n";
    for(int i=0; i<E_NUMENTRIES; i++)
    {
    parameterOutputFile << Energy[i]*1e6 << "\t" << Parameters::hcPERe/Energy[i]*1e3
                        << "\t" << PMMA2_RAYLEIGH[i]*1e-3 << "\t" << PMMA_RAYLEIGH[i]*1e-3
                        << "\t" << Pstyrene_RAYLEIGH[i]*1e-3 << "\n";
    }
    parameterOutputFile.close();



    // Material properties table

    G4MaterialPropertiesTable* scintCoreMaterialProperties = new G4MaterialPropertiesTable();
    G4MaterialPropertiesTable* innerCladMaterialProperties = new G4MaterialPropertiesTable();
    G4MaterialPropertiesTable* outerCladMaterialProperties = new G4MaterialPropertiesTable();

    // Set refractive indices

    scintCoreMaterialProperties->AddProperty("RINDEX",Energy,Pstyrene_RIND,E_NUMENTRIES)->SetSpline(true);
    innerCladMaterialProperties->AddProperty("RINDEX",Energy,PMMA_RIND,E_NUMENTRIES)->SetSpline(true);
    outerCladMaterialProperties->AddProperty("RINDEX",Energy,PMMA2_RIND,E_NUMENTRIES)->SetSpline(true);

    // Set absorption

    // todo : write functions in c++ code
    TF1 coreAbs("coreAbs",Parameters::GetInstance()->AbsorptionCore(),300,800);
    TF1 clad1Abs("clad1Abs",Parameters::GetInstance()->AbsorptionClad1(),300,800);
    TF1 clad2Abs("clad2Abs",Parameters::GetInstance()->AbsorptionClad2(),300,800);

    G4double* Pstyrene_ABSLENGTH = new G4double[E_NUMENTRIES];
    G4double* PMMA_ABSLENGTH = new G4double[E_NUMENTRIES];
    G4double* PMMA2_ABSLENGTH = new G4double[E_NUMENTRIES];

    // Calculate absorption lengths

    for(int j=0; j<E_NUMENTRIES; j++)
    {
        double wavelengthNanometer = Parameters::hcPERe/Energy[j]*1e3;

        Pstyrene_ABSLENGTH[j] = 1./coreAbs.Eval(wavelengthNanometer)*m;
        PMMA_ABSLENGTH[j] 	 = 1./clad1Abs.Eval(wavelengthNanometer)*m;
        PMMA2_ABSLENGTH[j] 	 = 1./clad2Abs.Eval(wavelengthNanometer)*m;
    }
    
    scintCoreMaterialProperties->AddProperty("ABSLENGTH",Energy,Pstyrene_ABSLENGTH,E_NUMENTRIES)->SetSpline(true);
    innerCladMaterialProperties->AddProperty("ABSLENGTH",Energy,PMMA_ABSLENGTH,E_NUMENTRIES)->SetSpline(true);
    outerCladMaterialProperties->AddProperty("ABSLENGTH",Energy,PMMA2_ABSLENGTH,E_NUMENTRIES)->SetSpline(true);
    
    delete[] Pstyrene_ABSLENGTH;
    delete[] PMMA_ABSLENGTH;
    delete[] PMMA2_ABSLENGTH;

    // Set Rayleigh scattering

    scintCoreMaterialProperties->AddProperty("RAYLEIGH",Energy,Pstyrene_RAYLEIGH,E_NUMENTRIES)->SetSpline(true);
    innerCladMaterialProperties->AddProperty("RAYLEIGH",Energy,PMMA_RAYLEIGH,E_NUMENTRIES)->SetSpline(true);
    outerCladMaterialProperties->AddProperty("RAYLEIGH",Energy,PMMA2_RAYLEIGH,E_NUMENTRIES)->SetSpline(true);

    // Set scintillation and WLS properties

    scintCoreMaterialProperties->AddProperty("FASTCOMPONENT",scintSpecVector);
    scintCoreMaterialProperties->AddProperty("SLOWCOMPONENT",scintSpecVector);
    scintCoreMaterialProperties->AddProperty("WLSCOMPONENT",wlsSpecVector);


    scintCoreMaterialProperties->AddConstProperty("SCINTILLATIONYIELD",
                                                  Parameters::GetInstance()->ScintillationYield()/keV);
    scintCoreMaterialProperties->AddConstProperty("RESOLUTIONSCALE",
                                                  Parameters::GetInstance()->ResolutionScale());
    scintCoreMaterialProperties->AddConstProperty("FASTTIMECONSTANT",
                                                  Parameters::GetInstance()->DecayTimeFast()*ns);
    scintCoreMaterialProperties->AddConstProperty("SLOWTIMECONSTANT",
                                                  Parameters::GetInstance()->DecayTimeSlow()*ns);
    scintCoreMaterialProperties->AddConstProperty("YIELDRATIO", Parameters::GetInstance()->YieldRatio());
    // Is set in "PhysicsList.hh" as well due to inconsistency in Geant4 (G4OpticalPhsysics default value)

    scintCoreMaterialProperties->AddProperty("WLSABSLENGTH",
                                             WlsAbsEnergy,WlsAbsLength,WLS_ABS_ENTRIES)->SetSpline(true);
	
    scintCoreMaterialProperties->AddConstProperty("WLSTIMECONSTANT", Parameters::GetInstance()->WlsDecayTime()*ns);

    // Assign material properties tables

    scintCoreMaterial->SetMaterialPropertiesTable(scintCoreMaterialProperties);
    innerCladdingMaterial->SetMaterialPropertiesTable(innerCladMaterialProperties);
    outerCladdingMaterial->SetMaterialPropertiesTable(outerCladMaterialProperties);

    // Set the Birks Constant for the Polystyrene scintillator
    scintCoreMaterial->GetIonisation()->SetBirksConstant(Parameters::GetInstance()->BirksConstant()*mm/MeV);

    delete[] Energy;
    delete[] ScintilEnergyDist;

    delete[] WlsEnergy;
    delete[] WlsEnergyDist;

    delete[] WlsAbsEnergy;
    delete[] WlsAbsLength;

    delete[] Vacuum_RIND;
    delete[] Vacuum_ABS;

    delete[] Pstyrene_RIND;
    delete[] PMMA_RIND;
    delete[] PMMA2_RIND;

    delete[] Pstyrene_RAYLEIGH;
    delete[] PMMA_RAYLEIGH;
    delete[] PMMA2_RAYLEIGH;
}



G4VPhysicalVolume* DetectorConstruction::Construct()
{
    DefineMaterials();
    DefineMaterialProperties();

    // Placing the world (experimental hall)
    G4Box* expHall_box = new G4Box("World",expHall_x,expHall_y,expHall_z);
    expHall_log = new G4LogicalVolume(expHall_box,Air,"World",0,0,0);
    expHall_phys = new G4PVPlacement(0,G4ThreeVector(),expHall_log,"World",0,false,0);


    // Epoxy
    G4double xEpoxy = 0.5*(Nk*xDist+2*Parameters::GetInstance()->SemiAxisZ()*mm)+5*mm;
    G4double yEpoxy = 0.5*(Nj*yDist+2*Parameters::GetInstance()->SemiAxisZ()*mm)+0.2*mm;
    G4Box* epoxyBox = new G4Box("EpoxyBox",xEpoxy, yEpoxy, scint_z);
    epoxyLog = new G4LogicalVolume(epoxyBox, Epoxy, "EpoxyBox", 0, 0, 0);
    epoxyPhy = new G4PVPlacement(0, G4ThreeVector() , epoxyLog, "EpoxyBox", expHall_log, false, 0);



   //--------------------------------------------------
    // Mirror for reflection at one of the end
    //--------------------------------------------------
  
    // Place the mirror only if the user wants the mirror
    if (fMirrorToggle) {  
  
       G4VSolid* solidMirror = new G4Box("Mirror",
                                         xEpoxy,
                                         yEpoxy,
                                         fMirrorZ);
   
       G4LogicalVolume* logicMirror = new G4LogicalVolume(solidMirror,
                                                          alu,
                                                          "Mirror");
 
     G4OpticalSurface* mirrorSurface = new G4OpticalSurface("MirrorSurface",
                                                              glisur,
                                                              ground,
                                                              dielectric_metal,
                                                              fMirrorPolish);
  
       G4MaterialPropertiesTable* mirrorSurfaceProperty =
                                                new G4MaterialPropertiesTable();
  
       G4double p_mirror[] = {2.00*eV, 3.47*eV};
       const G4int nbins = sizeof(p_mirror)/sizeof(G4double);
       G4double refl_mirror[] = {fMirrorReflectivity,fMirrorReflectivity};
       assert(sizeof(refl_mirror) == sizeof(p_mirror));
       G4double effi_mirror[] = {0, 0};
       assert(sizeof(effi_mirror) == sizeof(effi_mirror));
  
       mirrorSurfaceProperty->
                         AddProperty("REFLECTIVITY",p_mirror,refl_mirror,nbins);
       mirrorSurfaceProperty->
                         AddProperty("EFFICIENCY",p_mirror,effi_mirror,nbins);
  
       mirrorSurface -> SetMaterialPropertiesTable(mirrorSurfaceProperty);
  
       new G4PVPlacement(0,
                         G4ThreeVector(0.0,0.0,-scint_z-2*fMirrorZ),
                         logicMirror,
                         "Mirror",
                         expHall_log,
                         false,
                         0);
  
       new G4LogicalSkinSurface("MirrorSurface",logicMirror,mirrorSurface);
    }



    // ABS plastic
    G4double xABS = xEpoxy;
    G4double yABS = 2.5*mm;
    G4double zABS = scint_z;
    G4Box* absBox = new G4Box("AbsBox", xABS, yABS/2., zABS);
    absLog = new G4LogicalVolume(absBox, Abs_plastic, "AbsBox", 0, 0,0);
    absPhy = new G4PVPlacement(0, G4ThreeVector(0., -(yEpoxy+yABS/2.), 0.), absLog, "AbsBox", expHall_log, false, 0);

     ConstructFiber();
    // ConstructFiberSheet();

    /* ++ Construction and placement of the detector strips ++ */

    // Epoxy layer infront of det strip
    G4VSolid* epoxy_strip = new G4Box("EpoxyStrip",  stripWidth/2., stripHeight/2., epoxy_strip_width/2.);
    G4LogicalVolume* epoxy_strip_log = new G4LogicalVolume(epoxy_strip, Glue, "EpoxyStrip", 0, 0, 0);
    
    
    G4VSolid* pixelS = new G4Box("Pixel", pixelDimX/2., pixelDimY/2., stripWidth/2.);
    G4LogicalVolume* pixelL = new G4LogicalVolume(pixelS, Glue, "Pixel", 0, 0, 0);

    SensitiveDetector* sensitive = new SensitiveDetector("/Sensitive");
    G4SDManager* sdman = G4SDManager::GetSDMpointer();
    sdman->AddNewDetector(sensitive);
    pixelL->SetSensitiveDetector(sensitive);

    G4int Nk_s = ceil(((G4double)Nk)*(xDist/stripWidth)); // Determ. autom. nb. of detector strips.
    for(int k = 1; k < Nk_s-1; k++)
    {
        new G4PVPlacement(0, objectPos(Nj/2, k, scint_z+epoxy_strip_width/2.+airGap), epoxy_strip_log,
                          "EpoxyStrip", expHall_log, false, 0);

        for(int i = 0; i < Nx; i++)
        {
            for(int j = 0; j < Ny; j++){
                new G4PVPlacement(0, objectPos(Nj/2, i, j, k, scint_z+stripWidth/2.+epoxy_strip_width+airGap),
                                  pixelL, "SensitiveDetector"+C::c1(k)+C::c2(i)+C::c3(j), expHall_log, false, 0);
            }
        }
    }
	/* ++ End of Constr. and placem. of det. strips ++ */

    /*  Construction and placement of trigger
     *  Trigger should have the same x,y-dimensions as epoxyBox
     */
  
    G4double xTrigger = xEpoxy;//Parameters::GetInstance()->TriggerX()*mm;
    G4double yTrigger = Parameters::GetInstance()->TriggerY()*mm;
    G4double zTrigger = scint_z*2.0;//;Parameters::GetInstance()->TriggerZ()*mm;
  
    G4double triggerXPos = 0.0; //Parameters::GetInstance()->TriggerXPos()*mm;
    G4double triggerZPos = 0.0; //Parameters::GetInstance()->TriggerZPos()*mm;

    G4VSolid* triggerS = new G4Box("Trigger", xTrigger/2., yTrigger/2., zTrigger/2.);
    G4LogicalVolume* triggerL = new G4LogicalVolume(triggerS, Air, "Trigger", 0, 0, 0);
    new G4PVPlacement(0 , G4ThreeVector(triggerXPos, -(yEpoxy+yTrigger/2.+yABS), triggerZPos),
                      triggerL, "Trigger", expHall_log, false, 0);
	
    return expHall_phys;
}


void DetectorConstruction::ConstructFiberSheet()
{
    // dimensions
    G4double dim_z;
    G4double sphi, ephi;
    dim_z = scint_z;
    sphi = 0.00*deg;
    ephi = 360.*deg;
  G4double xEpoxy = 0.5*(Nk*xDist+2*Parameters::GetInstance()->SemiAxisZ()*mm)+5*mm;
            // Scintillating core
    //   G4Tubs* coreSection_tube = new G4Tubs("CoreSection",core_rZmin,core_rZmax,dim_z, sphi,ephi);
    //	    G4Box* coreSection_tube= new G4Box("CoreSection",  stripWidth/3., 0.85*mm, epoxy_strip_width/2.);
	    G4Box* CoreBox = new G4Box("CoreBox",xEpoxy-1*mm, 0.8294/2.0*mm, scint_z);
            G4LogicalVolume *coreSection_log = new G4LogicalVolume(CoreBox,
                                                                   scintCoreMaterial, "CoreSection",0,0,0);
            
            new G4PVPlacement(0, G4ThreeVector(), coreSection_log, "Core", epoxyLog, true, 0);

}
void DetectorConstruction::ConstructFiber()
{
    // dimensions
    G4double dim_z;
    G4double sphi, ephi;
    G4double core_rZmin,core_rZmax;
    G4double core_rYmin,core_rYmax;
    G4double clad1_rZmin,clad1_rZmax;
    G4double clad1_rYmin,clad1_rYmax;
    G4double clad2_rZmin,clad2_rZmax;
    G4double clad2_rYmin,clad2_rYmax;

    dim_z = scint_z;
    sphi = 0.00*deg;
    ephi = 360.*deg;

    core_rZmin = 0.00*cm;
    core_rYmin = 0.00*cm;

    core_rZmax = Parameters::GetInstance()->SemiAxisZ()*(88./100.)*mm;
    core_rYmax = Parameters::GetInstance()->SemiAxisY()*(88./100.)*mm;

    clad1_rZmin = core_rZmax;
    clad1_rYmin = core_rYmax;

    clad1_rZmax = core_rZmax + 3./88.*core_rZmax*2.;
    clad1_rYmax = core_rYmax + 3./88.*core_rYmax*2.;

    clad2_rZmin = clad1_rZmax;
    clad2_rYmin = clad1_rYmax;

    clad2_rZmax = Parameters::GetInstance()->SemiAxisZ()*mm;
    clad2_rYmax = Parameters::GetInstance()->SemiAxisY()*mm;

    G4ThreeVector origin;
    /* ++ Fibre Placement ++ */
    for(int j = 0; j < Nj; j++)
    {
        for(int k = 0; k < Nk; k++)
        {
            origin = objectPos(j,k);
      
            // Outer cladding
            G4Tubs* clad2Section_tube = new G4Tubs("Cladding2Section",clad2_rZmin,clad2_rZmax,dim_z, sphi,ephi);
            G4LogicalVolume *clad2Section_log = new G4LogicalVolume(clad2Section_tube,
                                                                    outerCladdingMaterial, "Cladding2Section",0,0,0);

            new G4PVPlacement(0, origin, clad2Section_log, "Cladding2"+C::c(j)+C::c(k), epoxyLog, true, 0);


            // Inner cladding
            G4Tubs* clad1Section_tube = new G4Tubs("Cladding1Section",clad1_rZmin,clad1_rZmax,dim_z, sphi,ephi);
            G4LogicalVolume *clad1Section_log = new G4LogicalVolume(clad1Section_tube,
                                                                    innerCladdingMaterial, "Cladding1Section",0,0,0);
      
            new G4PVPlacement(0, origin, clad1Section_log, "Cladding1"+C::c(j)+C::c(k), epoxyLog, true, 0);

            // Scintillating core
            G4Tubs* coreSection_tube = new G4Tubs("CoreSection",core_rZmin,core_rZmax,dim_z, sphi,ephi);
            G4LogicalVolume *coreSection_log = new G4LogicalVolume(coreSection_tube,
                                                                   scintCoreMaterial, "CoreSection",0,0,0);
            
            new G4PVPlacement(0, origin, coreSection_log, "Core"+C::c(j)+C::c(k), epoxyLog, true, 0);
        }
    }
	/* ++ End of Fibre Placement ++ */
}


// for fibre placement
G4ThreeVector DetectorConstruction::objectPos(G4int j, G4int k)
{	
    G4double xOffset = xDist*(((G4double)Nk)/2.)-Parameters::GetInstance()->SemiAxisZ()*mm/2.;
    G4double yOffset = yDist*(((G4double)Nj)/2.)-Parameters::GetInstance()->SemiAxisZ()*mm;
	
    G4double xDispl;
    G4double xsigma = (-0.49*j*j + 7.0*j-1.8)/1000.;
    G4double xvar = G4RandGauss::shoot(0.,xsigma);
    while(fabs(xvar)>=0.020*mm) {
      xvar = G4RandGauss::shoot(0.,xsigma);
    }
	
    j % 2 == 0 ? xDispl = 0. : xDispl = xDist/2.;
	
    G4ThreeVector origin(xDist*k+xDispl-xOffset+xvar, -yDist*j+yOffset, 0.);
    
    // print detector positions into file
    std::ofstream outFile;
    outFile.open("fibrePos.txt",std::ios::app);
    outFile << origin.x() << "\t" << origin.y() << std::endl;
    outFile.close();
    
    //G4cout << "Placing Fibre[" << k << ":" << j << "]: " << origin << G4endl;
    
    //G4cout << origin.x() << "\t" << origin.y() << G4endl;
    
    return origin;
}

// for detector placement
G4ThreeVector DetectorConstruction::objectPos(G4int j, G4int k, G4double zPos)
{	
    G4double randomN = Parameters::GetInstance()->RandomNumber();
    G4double offset = (stripWidth/2.)*randomN;
	
    // Save random number
    std::ofstream myfile;
    myfile.open("randomN.txt");
    myfile << randomN;
    myfile.close();
	
    G4double xOffset = xDist*(((G4double)Nk)/2.)-Parameters::GetInstance()->SemiAxisZ()*mm/2.;
    G4double yOffset = yDist*(((G4double)Nj)/2.)-Parameters::GetInstance()->SemiAxisZ()*mm;
	
    G4double xDispl;
	
    j % 2 == 0 ? xDispl = 0. : xDispl = xDist/2.;
	
    G4ThreeVector origin(stripWidth*k+xDispl-xOffset-offset, -yDist*j+yOffset, zPos);

    // print detector positions into file
    std::ofstream detFile;
    detFile.open("detPos.txt",std::ios::app);
    detFile << k << " " << stripWidth*k+xDispl-xOffset-offset << std::endl;
    detFile.close();
    
    return origin;
}

// for detector placement
G4ThreeVector DetectorConstruction::objectPos(G4int j, G4int i, G4int j2, G4int k, G4double zPos)
{	
    G4double randomN = Parameters::GetInstance()->RandomNumber();
    G4double offset = (stripWidth/2.)*randomN;
	
    // Save random number
    //std::ofstream myfile;
    //myfile.open("randomN.txt");
    //myfile << randomN;
    //myfile.close();
	
    G4double xOffset = xDist*(((G4double)Nk)/2.)-Parameters::GetInstance()->SemiAxisZ()*mm/2.;
    G4double yOffset = yDist*(((G4double)Nj)/2.)-Parameters::GetInstance()->SemiAxisZ()*mm + 3/2*pixelDimY;
	
    G4double xDispl;
	
    j % 2 == 0 ? xDispl = 0. : xDispl = xDist/2.;

    // print detector positions into file
    //std::ofstream detFile;
    //detFile.open("detPos.txt",std::ios::app);
    //detFile << k << " " << stripWidth*k+xDispl-xOffset-offset << std::endl;
    //detFile.close();

    G4double pixelX = pixelDimX/2. + i*pixelDimX - stripWidth/2.;
    G4double pixelY = pixelDimY/2. + j2*pixelDimY - stripHeight/2.;

    G4ThreeVector origin(stripWidth*k+xDispl-xOffset-offset + pixelX, -yDist*j+yOffset + pixelY, zPos);

    std::ofstream detFileY;
    detFileY.open("detPosY.txt",std::ios::app);
    detFileY << k << " " << -yDist*j+yOffset + pixelY << std::endl;
    detFileY.close();
	
    return origin;
}