//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
/// \file DetectorConstruction.cc
/// \brief Implementation of the DetectorConstruction class
//
//

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

#include "DetectorConstruction.hh"
#include "DetectorMessenger.hh"
#include "G4Material.hh"
#include "G4NistManager.hh"

#include "G4Box.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"

#include "G4GeometryManager.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4SolidStore.hh"
#include "G4RunManager.hh"

#include "G4UnitsTable.hh"
#include "G4SystemOfUnits.hh"
#include<cmath>
#include "G4SubtractionSolid.hh"
#include "G4Tubs.hh"
#include "G4VisAttributes.hh"

#include "G4Colour.hh"

#include "BiasMultiParticleChangeCrossSection.hh" //for biasing
#include "G4SDManager.hh" //for baising
#include "G4LogicalVolumeStore.hh"
#include "G4Cons.hh"
#include <cmath>
#include <set>

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

DetectorConstruction::DetectorConstruction()
:G4VUserDetectorConstruction(),
 fPBox(0), fLBox(0), fMaterial(0), fDetectorMessenger(0)
{
  fBoxSize = 1*m;
  DefineMaterials();
  SetMaterial("Water_ts");  
  fDetectorMessenger = new DetectorMessenger(this);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

DetectorConstruction::~DetectorConstruction()
{ delete fDetectorMessenger;}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4VPhysicalVolume* DetectorConstruction::Construct()
{
  return ConstructVolumes();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void DetectorConstruction::DefineMaterials()
{
  // specific element name for thermal neutronHP
  // (see G4ParticleHPThermalScatteringNames.cc)

  G4int ncomponents, natoms;

  // pressurized water
  G4Element* H  = new G4Element("TS_H_of_Water" ,"H" , 1., 1.0079*g/mole);
  G4Element* O  = new G4Element("Oxygen"        ,"O" , 8., 16.00*g/mole);
  G4Material* H2O = //new G4Material("Water_ts", 1.000*g/cm3, ncomponents=2,kStateLiquid, 593*kelvin, 150*bar); //old
                        new G4Material("Water_ts", 1.000*g/cm3, ncomponents=2,kStateLiquid, 293.15*kelvin, 1*atmosphere);
  H2O->AddElement(H, natoms=2);
  H2O->AddElement(O, natoms=1);
  H2O->GetIonisation()->SetMeanExcitationEnergy(78.0*eV);
  //germenium
   G4Material* germanium = new G4Material("Germanium", 32, 72.63*g/mole, 5.323 * g/cm3);
  // heavy water
  G4Isotope* H2 = new G4Isotope("H2",1,2);
  G4Element* D  = new G4Element("TS_D_of_Heavy_Water", "D", 1);
  D->AddIsotope(H2, 100*perCent);  
  G4Material* D2O = new G4Material("HeavyWater", 1.11*g/cm3, ncomponents=2,
                        kStateLiquid, 293.15*kelvin, 1*atmosphere);
  D2O->AddElement(D, natoms=2);
  D2O->AddElement(O, natoms=1);
  
  // graphite
  G4Isotope* C12 = new G4Isotope("C12", 6, 12);  
  G4Element* C   = new G4Element("TS_C_of_Graphite","C", ncomponents=1);
  C->AddIsotope(C12, 100.*perCent);
  G4Material* graphite = 
  new G4Material("graphite", 2.27*g/cm3, ncomponents=1,
                         kStateSolid, 293*kelvin, 1*atmosphere);
  graphite->AddElement(C, natoms=1);    
  
 ///G4cout << *(G4Material::GetMaterialTable()) << G4endl;
//******************GdCl3
G4Element * GD= new G4Element("Gadolinium", "Gd",64.,157.25*g/mole);
G4Element* CL= new G4Element("Chlorin","Cl", 17.,35.45*g/mole);
G4Material* GDCl3 =  new G4Material("GDCl3", 4.52*g/cm3, ncomponents=2);
GDCl3->AddElement(GD,natoms=1);
GDCl3->AddElement(CL,natoms=3);
double water_fraction = 0.98;
double gdcl3_fraction = 1 - water_fraction;

double total_density = water_fraction * 1.000*g/cm3 + gdcl3_fraction *4.52*g/cm3;
double water_density = water_fraction * 1.000*g/cm3 / total_density;
double gdcl3_density = gdcl3_fraction * 4.52*g/cm3 / total_density;
auto* gdcl3_water = new G4Material("GdCl3_water", total_density, 2);
gdcl3_water->AddMaterial(H2O, water_density);
gdcl3_water->AddMaterial(GDCl3, gdcl3_density);

// Glass
  fGlass = new G4Material("Glass", 1.032 * g / cm3, 2);
  fGlass->AddElement(C, 91.533 * perCent);
  fGlass->AddElement(H, 8.467 * perCent);
  
  //-----------------adding cerenkov property to water 
 std::vector<G4double> energy = {
    2.034 * eV, 2.068 * eV, 2.103 * eV, 2.139 * eV, 2.177 * eV, 2.216 * eV,
    2.256 * eV, 2.298 * eV, 2.341 * eV, 2.386 * eV, 2.433 * eV, 2.481 * eV,
    2.532 * eV, 2.585 * eV, 2.640 * eV, 2.697 * eV, 2.757 * eV, 2.820 * eV,
    2.885 * eV, 2.954 * eV, 3.026 * eV, 3.102 * eV, 3.181 * eV, 3.265 * eV,
    3.353 * eV, 3.446 * eV, 3.545 * eV, 3.649 * eV, 3.760 * eV, 3.877 * eV,
    4.002 * eV, 4.136 * eV
  };
   std::vector<G4double> rindexWater = {
    1.3435, 1.344,  1.3445, 1.345,  1.3455, 1.346,  1.3465, 1.347,
    1.3475, 1.348,  1.3485, 1.3492, 1.35,   1.3505, 1.351,  1.3518,
    1.3522, 1.3530, 1.3535, 1.354,  1.3545, 1.355,  1.3555, 1.356,
    1.3568, 1.3572, 1.358,  1.3585, 1.359,  1.3595, 1.36,   1.3608
  };

  std::vector<G4double> absorption = {
    3.448 * m,  4.082 * m,  6.329 * m,  9.174 * m,  12.346 * m, 13.889 * m,
    15.152 * m, 17.241 * m, 18.868 * m, 20.000 * m, 26.316 * m, 35.714 * m,
    45.455 * m, 47.619 * m, 52.632 * m, 52.632 * m, 55.556 * m, 52.632 * m,
    52.632 * m, 47.619 * m, 45.455 * m, 41.667 * m, 37.037 * m, 33.333 * m,
    30.000 * m, 28.500 * m, 27.000 * m, 24.500 * m, 22.000 * m, 19.500 * m,
    17.500 * m, 14.500 * m
  };
  std::vector<G4double> rindexWorld = {
    1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00,
    1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00,
    1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00
  };
  /*G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();

   myMPT1->AddProperty("RINDEX", photonEnergy, refractiveIndex1)
    ->SetSpline(true);
    //H2O->SetMaterialPropertiesTable(myMPT1);
    G4cout<<"cerenkov property added"<<G4endl;*/
     G4MaterialPropertiesTable *myMPT1 = new G4MaterialPropertiesTable();
    //mptAerogel->AddProperty("RINDEX", energy, rindexAerogel, 2);
    //Aerogel->SetMaterialPropertiesTable(mptAerogel);
    myMPT1->AddProperty("RINDEX", energy, rindexWater);
    myMPT1->AddProperty("ABSLENGTH", energy, absorption)->SetSpline(true);
    



//***********************
// the follwoing geometry has been taken from warwick_legend
//legend_geometry*********************legend_geometry********************legend_geometry****************************legend_geometry



/*
  G4NistManager* nistManager = G4NistManager::Instance();
  nistManager->FindOrBuildMaterial("G4_Galactic");
  nistManager->FindOrBuildMaterial("G4_lAr");
  nistManager->FindOrBuildMaterial("G4_AIR");
  nistManager->FindOrBuildMaterial("G4_STAINLESS-STEEL");
  nistManager->FindOrBuildMaterial("G4_Cu");
  nistManager->FindOrBuildMaterial("G4_WATER");

  auto* C  = new G4Element("Carbon", "C", 6., 12.011 * g / mole);
  auto* O  = new G4Element("Oxygen", "O", 8., 16.00 * g / mole);
  auto* Ca = new G4Element("Calcium", "Ca", 20., 40.08 * g / mole);
  auto* Mg = new G4Element("Magnesium", "Mg", 12., 24.31 * g / mole);

  // Standard Rock definition, similar to Gran Sasso rock
  // with density from PDG report
  auto* stdRock = new G4Material("StdRock", 2.65 * g / cm3, 4);
  stdRock->AddElement(O, 52.0 * perCent);
  stdRock->AddElement(Ca, 27.0 * perCent);
  stdRock->AddElement(C, 12.0 * perCent);
  stdRock->AddElement(Mg, 9.0 * perCent);

  auto* H     = new G4Element("Hydrogen", "H", 1., 1.00794 * g / mole);
  auto* N     = new G4Element("Nitrogen", "N", 7., 14.00 * g / mole);
  auto* puMat = new G4Material("polyurethane", 0.3 * g / cm3, 4);  // high density foam
  puMat->AddElement(H, 16);
  puMat->AddElement(O, 2);
  puMat->AddElement(C, 8);
  puMat->AddElement(N, 2);

  auto*    B10           = new G4Isotope("B10", 5, 10, 10 * g / mole);
  auto*    B11           = new G4Isotope("B11", 5, 11, 11 * g / mole);
  auto*    SpecialB      = new G4Element("SpecialB", "SpecialB", 2);
  G4double B_MassRatio   = 0.2;
  G4double B_NumberRatio = 1 / (1 + (1 - B_MassRatio) / B_MassRatio * 10. / 11.);
  SpecialB->AddIsotope(B10, B_NumberRatio);
  SpecialB->AddIsotope(B11, 1 - B_NumberRatio);

  // Borated Polyethylene according to GEM TN-92-172 TART Calculations of Neutron
  // Attenuation and Neutron-induced Photons on 5 % and 20 % Borated Polyethylene Slabs
  auto* BoratedPET =
    new G4Material("BoratedPET", 0.95 * g / cm3, 4);  // high density foam
  BoratedPET->AddElement(H, 0.116);
  BoratedPET->AddElement(C, 0.612);
  BoratedPET->AddElement(SpecialB, 0.05);
  BoratedPET->AddElement(O, 0.222);

  // Estimated using the number of elements per chain elements  (C_5 O_2 H_8)
  auto* PMMA = new G4Material("PMMA", 1.18 * g / cm3, 3);
  PMMA->AddElement(H, 0.08);
  PMMA->AddElement(C, 0.60);
  PMMA->AddElement(O, 0.32);

  // Estimated using the number of elements per molecule (C_2 H_4)
  auto* PolyEthylene = new G4Material("PolyEthylene", 0.95 * g / cm3, 2);
  PolyEthylene->AddElement(H, 0.142);
  PolyEthylene->AddElement(C, 0.857);

  G4Element* elGd = new G4Element("Gadolinium", "Gd", 64, 157.25 * g / mole);
  G4Element* elS  = new G4Element("Sulfur", "S", 16., 32.066 * g / mole);
  G4cout << elGd << G4endl;

  G4double density = 3.01 * g / cm3;  // https://www.sigmaaldrich.com/catalog/product/aldrich/203300?lang=de&region=DE
                                      // @room temp
  G4Material* gadoliniumSulfate =
    new G4Material("GadoliniumSulfate", density, 3);  // Gd2(SO4)3
  gadoliniumSulfate->AddElement(elGd, 2);
  gadoliniumSulfate->AddElement(elS, 3);
  gadoliniumSulfate->AddElement(O, 12);

  G4Material* purewater = G4Material::GetMaterial("G4_WATER");  // EDIT: changed water to purewater & use it to create "special" water
  auto water = new G4Material("GdLoadedWater", 1.000000 * g / cm3, 2);
  water->AddMaterial(purewater, 1. - 0.002);
  water->AddMaterial(gadoliniumSulfate, 0.002);

  // enriched Germanium from isotopes
  auto* Ge_74 = new G4Isotope("Ge74", 32, 74, 74.0 * g / mole);
  auto* Ge_76 = new G4Isotope("Ge76", 32, 76, 76.0 * g / mole);

  G4Isotope* Ge_70 = new G4Isotope("Ge70",  32, 70, 69.92*g/mole);
  G4Isotope* Ge_72 = new G4Isotope("Ge72",  32, 72, 71.92*g/mole);
  G4Isotope* Ge_73 = new G4Isotope("Ge73",  32, 73, 73.0*g/mole);

  G4Element* eGe = new G4Element("enriched Germanium", "enrGe", 2);

  eGe->AddIsotope(Ge_76, 87. * perCent);
  eGe->AddIsotope(Ge_74, 13. * perCent);
  density      = 5.323 * g / cm3;

  if(fMaGeMaterial){
    density = 5.56 * g / cm3;
    eGe = new G4Element("enriched Germanium", "enrGe", 5);  
    eGe->AddIsotope(Ge_70,0.0*perCent);
    eGe->AddIsotope(Ge_72,0.1*perCent);
    eGe->AddIsotope(Ge_73,0.2*perCent);
    eGe->AddIsotope(Ge_74,13.1*perCent);
    eGe->AddIsotope(Ge_76,86.6*perCent);
  }

  auto* roiMat = new G4Material("enrGe", density, 1);
  roiMat->AddElement(eGe, 1);

  // Edit: 2020/02/17 by Moritz Neuberger
  // Added new def. of LAr in order to add doping with Xe and He-3

  G4double dLAr  = 1.393 * g / cm3;
  G4double dLXe  = 3.02 * g / cm3;
  G4double dLHe3 = 0.059 * g / cm3;

  G4double fArConc = 1 - (fXeConc + fHe3Conc);

  G4double dComb = 1 / ((fArConc / dLAr) + (fXeConc / dLXe) + (fHe3Conc / dLHe3));

  G4cout << "___________________________________________" << G4endl;
  G4cout << "Mass ratios of cryostat:" << G4endl;
  G4cout << "LAr:   " << fArConc << G4endl;
  G4cout << "LXe:   " << fXeConc << G4endl;
  G4cout << "LHe3:   " << fHe3Conc << G4endl;
  G4cout << "dComb:   " << dComb << G4endl;
  G4cout << "___________________________________________" << G4endl;

  
  //auto* eLAr = new G4Element("LAr", "Ar", 18., 39.95 * g / mole);
  larMat = G4Material::GetMaterial("G4_lAr");
  auto*      eLXe   = new G4Element("LXe", "Xe", 54., 131.29 * g / mole);
  auto*      eHe3   = new G4Element("He3", "He3", 1);
  G4Isotope* iHe3   = new G4Isotope("He3", 2, 3);
  eHe3->AddIsotope(iHe3, 1);

  if(fMaGeMaterial){
    double densityLAr = 1.396  * g / cm3; 
    G4Element* larEl = new G4Element("LAr", "Ar", 18., 39.95 * g / mole);
    larMat = new G4Material("LAr", densityLAr, 1);
    larMat->AddElement(larEl, 1);
  }

  CombinedArXeHe3 =
    new G4Material("CombinedArXeHe3", dComb, 3, kStateLiquid, 87. * kelvin);
  CombinedArXeHe3->AddMaterial(larMat, fArConc);
  CombinedArXeHe3->AddElement(eHe3, fHe3Conc);
  CombinedArXeHe3->AddElement(eLXe, fXeConc);

  G4cout << CombinedArXeHe3 << G4endl;*/


//legend_geometry*********************legend_geometry********************legend_geometry****************************legend_geometry



}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4Material* DetectorConstruction::MaterialWithSingleIsotope( G4String name,
                           G4String symbol, G4double density, G4int Z, G4int A)
{
 // define a material from an isotope
 //
 G4int ncomponents;
 G4double abundance, massfraction;

 G4Isotope* isotope = new G4Isotope(symbol, Z, A);
 
 G4Element* element  = new G4Element(name, symbol, ncomponents=1);
 element->AddIsotope(isotope, abundance= 100.*perCent);
 
 G4Material* material = new G4Material(name, density, ncomponents=1);
 material->AddElement(element, massfraction=100.*perCent);

 return material;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4VPhysicalVolume* DetectorConstruction::ConstructVolumes()
{   
//Air



  // Cleanup old geometry
  G4GeometryManager::GetInstance()->OpenGeometry();
  G4PhysicalVolumeStore::GetInstance()->Clean();
  G4LogicalVolumeStore::GetInstance()->Clean();
  G4SolidStore::GetInstance()->Clean();

 /* G4Box*
  sBox = new G4Box("Container",                         //its name
                   fBoxSize/2,fBoxSize/2,fBoxSize/2);   //its dimensions

  fLBox = new G4LogicalVolume(sBox,                     //its shape
                             fMaterial,                 //its material
                             fMaterial->GetName());     //its name

  fPBox = new G4PVPlacement(0,                          //no rotation
                            G4ThreeVector(),            //at (0,0,0)
                            fLBox,                      //its logical volume
                            fMaterial->GetName(),       //its name
                            0,                          //its mother  volume
                            false,                      //no boolean operation
                            0);                         //copy number
     
                           
  PrintParameters();
  return fPBox; */
  
  //always return the root volume
  //
  //Geometry Josef Jochum~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~april 23

  G4double stone       = 100.0;  // Hall wall thickness 1 m
  G4double hallrad     = 600.0;  // Hall cylinder diam 12 m
  //G4double hallrad = fBoxSize;
  G4double hallhheight = 850.0;  // Hall cylinder height 17 m
  // water tank
  G4double tankwalltop = 0.6;  // water tank thickness at top 6 mm
  G4double tankwallbot = 0.8;  // water tank thickness at bottom 8 mm
  G4double tankrad     = 550;  // water tank diam 11 m
  //G4double tankrad =fBoxSize;
  G4double tankhheight = 650;  // water tank height 13 m
  // cryostat
  G4double cryowall   = 3.0;    // cryostat wall thickness from GERDA
  G4double vacgap     = 1.0;    // vacuum gap between walls
  G4double cryrad     = 350.0;  // cryostat diam 7 m
  G4double cryhheight = 350.0;  // cryostat height 7 m

    //mother
    /*G4double len = 20 *m;
    G4Box* motherV = new G4Box("motherV",len/2,len/2,len/2);
    G4Material* mothMat = G4NistManager::Instance()->FindOrBuildMaterial("G4_Galactic");
    G4LogicalVolume* motherLogical = new G4LogicalVolume(motherV, mothMat, "motherLogical");
    G4VPhysicalVolume* motherPhysical = new G4PVPlacement(0, G4ThreeVector(), motherLogical, "motherPhysical", nullptr, false,0);*/
     
    
    // World volume (hall) this worked perfectly
    
    G4double hallRadius = 6*m;
    G4double hallHeight = 12.*m;
    G4Tubs* hallSolid = new G4Tubs("hallSolid", 0, hallRadius+100.*cm, hallHeight/2+100.*cm, 0, 2*M_PI);
    G4Material* air = G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR");
    G4LogicalVolume* hallLogical = new G4LogicalVolume(hallSolid, air, "hallLogical");
    G4VPhysicalVolume* hallPhysical = new G4PVPlacement(0, G4ThreeVector(), hallLogical, "hallPhysical",nullptr, false, 0);

    // LAr cylinder
    G4double larRadius = 3.47*m;
    G4double larHeight = 8.0*m;
    G4Tubs* larSolid = new G4Tubs("larSolid", 0, larRadius, larHeight/2, 0, 2*M_PI);
    G4Material* lar = G4NistManager::Instance()->FindOrBuildMaterial("G4_lAr");
    G4LogicalVolume* larLogical = new G4LogicalVolume(larSolid, lar, "larLogical");
    G4VPhysicalVolume* larPhysical=new G4PVPlacement(0, G4ThreeVector(), larLogical, "larPhysical", hallLogical, false, 0);

    // Steel layer
    G4double steelThickness = 0.03*m;
    G4double steelRadius = larRadius + steelThickness;
    G4Tubs* steelSolid = new G4Tubs("steelSolid", larRadius, steelRadius, larHeight/2, 0, 2*M_PI);
    G4Material* steel = G4NistManager::Instance()->FindOrBuildMaterial("G4_STAINLESS-STEEL");
    G4LogicalVolume* steelLogical = new G4LogicalVolume(steelSolid, steel, "steelLogical");
    G4VPhysicalVolume* steelPhysical=new G4PVPlacement(0, G4ThreeVector(), steelLogical, "steelPhysical", hallLogical, false, 0);

    // Water layer
    //fMaterial,                 //its material
     //fMaterial->GetName()
    G4double waterThickness = fBoxSize;
    G4double waterRadius = steelRadius + waterThickness;
    //G4Tubs* waterSolid = new G4Tubs("waterSolid", steelRadius, hallRadius, hallHeight/2, 0, 2*M_PI);
    G4Tubs* waterSolid = new G4Tubs("waterSolid", steelRadius, waterRadius-10.*cm, larHeight/2, 0, 2*M_PI);
    //G4Material* water = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
    G4LogicalVolume* fWaterLogical = new G4LogicalVolume(waterSolid, fMaterial, fMaterial->GetName());
    //G4LogicalVolume*  fLBox  = new G4LogicalVolume(waterSolid, fMaterial, fMaterial->GetName());
    G4VPhysicalVolume* waterPhysical    =new G4PVPlacement(0, G4ThreeVector(),  fWaterLogical , fMaterial->GetName(), hallLogical, false, 0); // this has been changed to fLBox becuase the macro written for this will turn water as logical volume
    
    // Set visualization attributes
  G4VisAttributes* LArVisAtt = new G4VisAttributes(G4Colour::Blue());
  larLogical->SetVisAttributes(LArVisAtt);

  G4VisAttributes* steelVisAtt = new G4VisAttributes(G4Colour::Grey());
  steelLogical->SetVisAttributes(steelVisAtt);

  G4VisAttributes* waterVisAtt = new G4VisAttributes(G4Colour::Cyan());
  fWaterLogical->SetVisAttributes(waterVisAtt);
  

/*G4double pmtRadius = 2.0*cm; // PMT radius
G4double pmtSpacing = .5*cm; // spacing between PMTs
G4int numPmtsPerRing = 700;//36; // number of PMTs per ring
//G4int numRings = 200.; // number of rings
G4Box* solidDetector = new G4Box("PMT", 2.*cm, 2*cm, 2*cm);
//G4Tubs* solidDetector = new G4Tubs("PMT",0.,1.0*cm,1.0*cm,0.,2*M_PI);
G4LogicalVolume* logicDetector = new G4LogicalVolume(solidDetector, fGlass, "PMT_log");  
//G4double waterRadius = steelRadius + waterThickness; // water cylinder radius
G4int numRings = std::ceil(800.0/cm/2.*pmtSpacing);
for (G4int ringNum = 0; ringNum < numRings; ringNum++) {
 // G4double ringRadius = waterRadius-10.0*cm; /*(pmtRadius + pmtSpacing) * (ringNum + 0.5);*/
  /*G4double ringRadius = hallRadius-10.0*cm;
  G4double circumference = 2.0*M_PI*ringRadius;
  G4double pmtAngle = 360.0/numPmtsPerRing*deg;
  
  for (G4int pmtNum = 0; pmtNum < numPmtsPerRing; pmtNum++) {
    G4double pmtX = ringRadius*cos(pmtAngle*pmtNum);
    G4double pmtY = ringRadius*sin(pmtAngle*pmtNum);
    G4double pmtZ = -larHeight/2.0+ringNum*4.*cm;
    G4ThreeVector pmtPos(pmtX, pmtY, pmtZ);
    G4RotationMatrix pmtRot;
    pmtRot.rotateX(90.0*deg);
    pmtRot.rotateZ(pmtAngle*pmtNum);
    G4Transform3D pmtTrans(pmtRot, pmtPos);
    //G4String pmtName = "pmt_" + std::to_string(ringNum) + "_" + std::to_string(pmtNum);
    new G4PVPlacement(pmtTrans, logicDetector, "pmt_log", hallLogical, false, 0);
  }
}
 */

    return hallPhysical; 
    
    //======================================================================LEGEND GEOMETRY==========================================================
    
                           
    
    //construct pmts
    
   
      
    // Define dimensions and positions of PMT and water cylinder
/*G4double pmtLength = 1. * cm;
G4double pmtWidth = 1.0 * cm;
G4double pmtThickness = .5 * cm;
G4double pmtSpacing = 0.1 * cm;
G4double waterRadius = tankrad * cm;
G4double waterHeight = (tankhheight - tankwallbot) * cm;

// Define number of rings and PMTs per ring
G4int numRings = 20;
G4int numPMTsPerRing = 1000;

// Loop over each ring
for (G4int ringNum = 0; ringNum < numRings; ringNum++) {
    // Define height of current ring
    G4double ringHeight = (2 * ringNum + 1) * pmtThickness;

    // Loop over each PMT in the ring
    for (G4int pmtNum = 0; pmtNum < numPMTsPerRing; pmtNum++) {
        // Calculate the angle and position of the PMT
        G4double angle = pmtNum * 360.0 / numPMTsPerRing;
        G4RotationMatrix* rot = new G4RotationMatrix();
        rot->rotateZ(angle * deg);
        G4double xPos = (waterRadius*cm + pmtWidth/2.0+pmtSpacing) * cos(angle * deg);
        G4double yPos = (waterRadius*cm + pmtWidth/2.0+pmtSpacing) * sin(angle * deg);
        G4double zPos = (waterHeight*cm)/2.0 + ringHeight;

        // Create and place PMT volume
        G4Box* pmtSolid = new G4Box("PMT", pmtLength/2.0, pmtWidth/2.0, pmtThickness/2.0);
        G4LogicalVolume* pmtLogical = new G4LogicalVolume(pmtSolid, fGlass, "PMT_log");
        new G4PVPlacement(rot, G4ThreeVector(xPos, yPos, zPos), pmtLogical, "PMT_phys", fWaterLogical, false, 0, true);
    }
}

    */



 // SetLogicalVolume(fWaterLogical); //new line for setting water as logicalvolume
 // fScoringVolume =fWaterLogical;

    //return fWorldPhysical;
    //return fCavernPhysical; //use cavern as world
    //====================================================================LEGEND GEOMETRY ENDED =======================================================

              
  PrintParameters();

  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  
  
  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  
  //PMT's 
   // Define PMT dimensions
/*G4double pmt_width = 2.0 * cm;
G4double pmt_height = 2.0 * cm;
G4double pmt_thickness = 1.0 * cm;

// Define cylinder dimensions
G4double cylinder_radius = 7.0 * m;
G4double cylinder_height = 10.0 * m;

// Create PMT solid and logical volume
G4Box* pmt_box = new G4Box("PMT_box", 0.5 * pmt_width, 0.5 * pmt_height, 0.5 * pmt_thickness);
G4LogicalVolume* pmt_log = new G4LogicalVolume(pmt_box, pmt_mat, "PMT_log");

// Calculate the number of rows and columns needed to cover the cylinder
G4int num_rows = (G4int)(2.0 * cylinder_radius / pmt_width);
G4int num_cols = (G4int)(cylinder_height / pmt_height);

// Create a rotation matrix for the PMT placement
G4RotationMatrix* rot = new G4RotationMatrix();
rot->rotateX(90.0 * degree);

// Place the PMTs on the curved surface of the cylinder
for (G4int i = 0; i < num_rows; i++) {
  for (G4int j = 0; j < num_cols; j++) {
    G4double x = (i + 0.5) * pmt_width - cylinder_radius;
    G4double y = (j + 0.5) * pmt_height - 0.5 * cylinder_height;
    G4double z = 0.0;

    G4ThreeVector pos(x, y, z);
    G4Transform3D transform(*rot, pos);
    G4VPhysicalVolume* pmt_phys = new G4PVPlacement(transform, pmt_log, "PMT_phys", world_log, false, i * num_cols + j, true);
  }
}

  */
  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  
  
  
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
///construnting SDAndField()----------------for biaisng and sensitivedetector
void  DetectorConstruction::ConstructSDandField(){
  //G4SDManager::GetSDMpointer()->SetVerboseLevel(1);
  
 if(!fSD.Get())
  { 
     G4String crystalSDname  = "CrystalSD";
    CrystalSD* aCrystalSD = new CrystalSD(crystalSDname,
                                                  "CrystalHitsCollection");
    fSD.Put(aCrystalSD);
    G4SDManager::GetSDMpointer()->AddNewDetector(fSD.Get());

   // SetSensitiveDetector("PMT_log", fSD.Get());

  
    G4LogicalVolumeStore* volumeStore = G4LogicalVolumeStore::GetInstance();
  
    auto* biasnXS = new BiasMultiParticleChangeCrossSection();
    biasnXS->SetNeutronFactor(fNeutronBias);
    biasnXS->SetMuonFactor(fMuonBias);
    G4cout << " >>> Detector: set neutron bias to " << fNeutronBias << G4endl;
    biasnXS->AddParticle("neutron");
    //G4LogicalVolume* logicGe = volumeStore->GetVolume("Water_log"); 
    G4LogicalVolume* logicGe = volumeStore->GetVolume(fMaterial->GetName()); //biasing water volume to enahnce the nCapture
    biasnXS->AttachTo(logicGe);

    // -- Attach muon XS biasing to all required volumes consistently
    auto* biasmuXS = new BiasMultiParticleChangeCrossSection();
    biasmuXS->SetNeutronFactor(fNeutronBias);
    biasmuXS->SetMuonFactor(fMuonBias);
    G4cout << " >>> Detector: set muon bias to " << fMuonBias << G4endl;
    biasmuXS->AddParticle("mu-");
  
  // add sensitive detector for cernkov 
  //MySensitiveDetector *sensDet = new MySensitiveDetector("SensitiveDetector");
 // logicDetector->SetSensitiveDetector(sensDet);

 } 
}
//constructing SDandField() ended here
void DetectorConstruction::PrintParameters()
{
  G4cout << "\n The Box is " << G4BestUnit(fBoxSize,"Length")
         << " of " << fMaterial->GetName() 
         << "\n \n" << fMaterial << G4endl;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//the following are commented out because to set logicaldetector and dimension of detecor as default 
void DetectorConstruction::SetMaterial(G4String materialChoice)
{
  // search the material by its name
  G4Material* pttoMaterial =
     G4NistManager::Instance()->FindOrBuildMaterial(materialChoice);   
  
 if (pttoMaterial) { 
    if(fMaterial != pttoMaterial) {
      fMaterial = pttoMaterial;
      //if(fLBox) { fLBox->SetMaterial(pttoMaterial); } //old
      if(fWaterLogical) {fWaterLogical->SetMaterial(pttoMaterial); 
       fMaterial->SetMaterialPropertiesTable(myMPT1);
      }
      
      G4RunManager::GetRunManager()->PhysicsHasBeenModified();
    }
  } else {
    G4cout << "\n--> warning from DetectorConstruction::SetMaterial : "
           << materialChoice << " not found" << G4endl;
  }            
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void DetectorConstruction::SetSize(G4double value)
{
  fBoxSize = value;
  G4RunManager::GetRunManager()->ReinitializeGeometry();
}
//new function  for biasing
void  DetectorConstruction::SetNeutronBiasFactor(G4double nf) { fNeutronBias = nf; }


void DetectorConstruction::SetMuonBiasFactor(G4double mf) { fMuonBias = mf; }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......