#include "CloverDetectorConstruction.hh"
#include "CloverCrystalSD.hh"
#include "Detector_Parameters.hh"

#include "G4Material.hh"
#include "G4NistManager.hh"

#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4Sphere.hh"
#include "G4Trd.hh"
#include "G4UnionSolid.hh"
#include "G4IntersectionSolid.hh"
#include "G4SubtractionSolid.hh"

#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4PVReplica.hh"
#include "G4GlobalMagFieldMessenger.hh"
#include "G4AutoDelete.hh"

#include "G4RotationMatrix.hh"

#include "G4SDManager.hh"

#include "G4VisAttributes.hh"
#include "G4Colour.hh"

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

#include <vector>


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

G4ThreadLocal 
G4GlobalMagFieldMessenger* CloverDetectorConstruction::fMagFieldMessenger = 0; 

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

CloverDetectorConstruction::CloverDetectorConstruction()
 : G4VUserDetectorConstruction(),
   fNumOfCrystal(0),
   fLogicCrystal(NULL), fCrystalMaterial(NULL),
   fCheckOverlaps(true)
{
  fNumOfCrystal = Crystal_Num;  //also need to change the fNDet in CloverEventAction.cc
  fLogicCrystal = new G4LogicalVolume* [fNumOfCrystal];
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
CloverDetectorConstruction::~CloverDetectorConstruction()
{
  delete [] fLogicCrystal;
}

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

G4VPhysicalVolume* CloverDetectorConstruction::Construct()
{
  // Define materials
  DefineMaterials();
  
  // Define volumes
  return DefineVolumes();
}

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

void CloverDetectorConstruction::DefineMaterials()
{ 
  // Lead material defined using NIST Manager
  auto nistManager = G4NistManager::Instance();
  
  // Vacuum
  new G4Material("Galactic", 1., 1.01*g/mole, universe_mean_density, kStateGas, 2.73*kelvin, 3.e-18*pascal);

  // Air
  nistManager->FindOrBuildMaterial("G4_AIR");

  // Iron
  nistManager->FindOrBuildMaterial("G4_Fe");

  // Aluminium
  nistManager->FindOrBuildMaterial("G4_Al");

  // Ge
  fCrystalMaterial = nistManager->FindOrBuildMaterial("G4_Ge");

  // Print materials
  G4cout << *(G4Material::GetMaterialTable()) << G4endl;
}

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

G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
{
  // Geometry parameters
  G4double crystalLength = 8.*cm;
  G4double crystalRadius = 25.*mm;
  G4double crystalZPos = Radius*mm + crystalLength/2.; //distance to target + distance to center of crystals.;
  G4double cutXY = 46.0 * mm;

  G4double worldSizeZ  = 2*(crystalLength + crystalZPos);
  G4double worldSizeXY = worldSizeZ;
/*
  //pipe
  G4double pipeOuterRadius     = 100/2.* mm;
  G4double pipeWallThickness   =  5 * mm;
  G4double pipeLength          = 70* mm;
  G4double pipeZPos            = 0.*cm;
*/

  //spherical chamber
  G4double chamberOuterRadius     = 150/2.* mm;
  G4double chamberWallThickness   =  5 * mm;
  G4double chamberInnerRadius      = chamberOuterRadius - chamberWallThickness;
  G4double chamberZPos            = 0.*cm;


  //Clover casing
  G4double caseXYInner = 2*(cutXY + 1.0) * mm;
  G4double caseZInner  = crystalLength/2 + crystalLength + 10*mm;
  G4double caseXYWallThickness = 3.0 * mm;
  G4double caseZWallThickness = 3.0 * mm;


  //BGO Detector 
  
  G4double bgoXYInner = 10*mm + 2*(cutXY + 1.0) * mm;
  G4double bgoZInner  = crystalLength/2 + 4 * mm + crystalLength + 15*mm;
  G4double bgoXYWallThickness = 103.5 * mm; //from schematic
  G4double bgoanglecut0 = 53 * mm ;
  G4double zbgo0 = 178 * mm;
  
  G4double bgoanglecut1 = 24.935 * mm;
  G4double zbgo1 = 80*mm;

  G4double bgoanglecut2 = 10 *mm ;
  G4double zbgo2 = 20*mm;
  
  // Get materials
  auto backgroundMaterial = G4Material::GetMaterial("G4_AIR");
  //auto backgroundMaterial = G4Material::GetMaterial("Galactic");
   
  // World
  auto worldS = new G4Box("World",  worldSizeXY/2, worldSizeXY/2, worldSizeZ); 
                         
  auto worldLV = new G4LogicalVolume( worldS,  backgroundMaterial, "World"); 
                                   
  auto worldPV = new G4PVPlacement( 0,                // no rotation
                                    G4ThreeVector(),  // at (0,0,0)
                                    worldLV,          // its logical volume                         
                                    "World",          // its name
                                    0,                // its mother  volume
                                    false,            // no boolean operation
                                    0,                // copy number
                                    fCheckOverlaps);  // checking overlaps

  worldLV->SetVisAttributes (G4VisAttributes::GetInvisible());

  // vaccum Pipe
  
 // G4Tubs * pipe = new G4Tubs("pipe", pipeOuterRadius - pipeWallThickness , pipeOuterRadius, pipeLength, 0, 360*degree);
//  G4Tubs * pipe0 = new G4Tubs("pipe0", 0 , pipeOuterRadius, pipeLength, 0, 360*degree);
 
  G4Sphere *chamber = new G4Sphere("chamber",chamberInnerRadius,chamberOuterRadius,0,2*pi,0,pi);

  G4RotationMatrix * rotX = new G4RotationMatrix(); rotX->rotateX(90*degree);
  G4RotationMatrix * rotY = new G4RotationMatrix(); rotY->rotateY(90*degree);
  
  G4LogicalVolume * chamberLV = new G4LogicalVolume(chamber, G4Material::GetMaterial("G4_Fe"),"targetchamber");
  
  new G4PVPlacement( 0,                // no rotation
                     G4ThreeVector(0, 0, chamberZPos),  // at (0,0,0)
                     chamberLV,          // its logical volume                         
                     "targetchamber",          // its name
                     worldLV,          // its mother  volume
                     false,            // no boolean operation
                     0,                // copy number
                     fCheckOverlaps);  // checking overlaps
    
  chamberLV->SetVisAttributes(new G4VisAttributes(G4Colour(1.0,1.0,1.0)));
  

  // Crystals
  G4VisAttributes * crystalVisAtt= new G4VisAttributes(G4Colour(0.5,0.5,1.0));
  crystalVisAtt->SetVisibility(true);

//define crystal shapes

  G4Tubs * base = new G4Tubs("base",  0, crystalRadius, crystalLength/2., 0, 360*degree);
  G4Box * cut = new G4Box("cut", cutXY/2. , cutXY/2. , crystalLength/2. * 1.2);
  G4IntersectionSolid * crystalS = new G4IntersectionSolid("HPGe", base, cut);

  G4Box* case1 = new G4Box("case1", caseXYInner/2., caseXYInner/2., caseZInner  );
  G4Box* case2 = new G4Box("case2", (caseXYInner + caseXYWallThickness)/2. , (caseXYInner + caseXYWallThickness)/2., caseZInner + caseZWallThickness*2.);
  G4SubtractionSolid * casing = new G4SubtractionSolid("casing", case2, case1);
  
  G4Box* bgo1 = new G4Box("bgo1", bgoXYInner/2., bgoXYInner/2., bgoZInner +zbgo1);
 
  G4Trd* bgo2 = new G4Trd("bgo2", (bgoXYInner+ bgoXYWallThickness)/2.  - bgoanglecut1 ,(bgoXYInner+ bgoXYWallThickness)/2.
                                , (bgoXYInner + bgoXYWallThickness)/2. - bgoanglecut1, (bgoXYInner + bgoXYWallThickness)/2., zbgo0/2);
    
  G4Trd* bgo3 = new G4Trd("bgo3", (bgoXYInner+ bgoXYWallThickness)/2.  - bgoanglecut0 + bgoanglecut2,(bgoXYInner+ bgoXYWallThickness)/2. - bgoanglecut1
                                , (bgoXYInner + bgoXYWallThickness)/2. - bgoanglecut0 + bgoanglecut2, (bgoXYInner + bgoXYWallThickness)/2.- bgoanglecut1, zbgo1/2);
  G4ThreeVector bgo3_v = G4ThreeVector(0.,0.,-zbgo0 + 50*mm);  
   
  G4RotationMatrix * bgo3_m = new G4RotationMatrix(); 


  G4VSolid* bgo_u = new G4UnionSolid("bgo2 + bgo3",bgo2,bgo3,bgo3_m,bgo3_v);
  //G4Box* bgo2 = new G4Box("bgo2", (bgoXYInner+ bgoXYWallThickness)/2., (bgoXYInner + bgoXYWallThickness)/2., bgoZInner + bgoZWallThickness*2.);

 
  G4Trd* bgowindow = new G4Trd("bgowindow", (bgoXYInner+ bgoXYWallThickness)/2.  - bgoanglecut0 ,(bgoXYInner+ bgoXYWallThickness)/2. - bgoanglecut0 +bgoanglecut2
                                , (bgoXYInner + bgoXYWallThickness)/2. - bgoanglecut0, (bgoXYInner + bgoXYWallThickness)/2. -bgoanglecut0 +bgoanglecut2, zbgo2/2);



  G4SubtractionSolid * bgo = new G4SubtractionSolid("bgo", bgo_u, bgo1);
 
  // creating the subtraction solid that will form the actual bgo sheild.

  
  
  double dTheta,dPsi,zmod; 

  for( G4int i = 0 ; i < fNumOfCrystal ; i++){ 
    G4String name = "HPGe"+ std::to_string(i);
    G4cout << " crystal name : " << name << G4endl;
    G4double phi = 360./4. * degree;

    G4double rho = 0;
    rho = (cutXY/2  + 0.3* mm) / sin(phi/2.);

    int iClover  = i/4; //How to get each group of crystals.
    int jCrystal = i%4; //How to renumber each group as 0,1,2,3. 
    
    fLogicCrystal[i] = new G4LogicalVolume(crystalS, fCrystalMaterial, "CrystalLV");


    auto caseLV =  new G4LogicalVolume * [fNumOfCrystal/4];
    auto bgoLV  =  new G4LogicalVolume * [fNumOfCrystal/4];
    auto bgowindowLV  =  new G4LogicalVolume * [fNumOfCrystal/4];

    fLogicCrystal[i]->SetVisAttributes(crystalVisAtt);

        
    if(int j = iClover){    
        
        dTheta = Clover_Location[j][1];

        dPsi = Clover_Location[j][2];

        zmod = Clover_Location[j][3];

        }
        

    G4ThreeVector pos = G4ThreeVector(rho * cos(jCrystal* phi + phi/2.), rho * sin(jCrystal * phi + phi/2.), crystalZPos+zmod);
    
//    if(iClover != DetNum){return 0;}


    G4ThreeVector case_pos = G4ThreeVector(0,0,88*mm+crystalZPos); //The crystalZpos is the radius of the detector,
                                                                   //bgo uses an identical rotation matrix as the case, 
                                                                   //you can also add in zmod if you want to shift the shielding. 
    G4ThreeVector bgo_pos = G4ThreeVector(0,0,115*mm+crystalZPos);
    G4ThreeVector bgowindow_pos = G4ThreeVector(0,0,-63*mm + crystalZPos);
                                                                          
    //Set rotation matrix settings for clovers.
    
    G4RotationMatrix * bgowindowmat = new G4RotationMatrix(); 
    G4RotationMatrix * bgomat = new G4RotationMatrix(); 
    G4RotationMatrix * rotmat = new G4RotationMatrix(); 
    G4RotationMatrix * crystalmat= new G4RotationMatrix();

    //x = sin(theta)cos(phi) || y = sin(theta)sin(phi) || z = cos(theta) 
   
  
    for(int j = 0; j <= iClover;j++){
        if(j == iClover){
 
    // rotation matrix must undo the operations you perform with the position vector. 
    bgowindowmat->rotateZ(-dPsi*degree);
    bgowindowmat->rotateX(-dTheta*degree);

    bgomat->rotateZ(-dPsi*degree);
    bgomat->rotateX(-dTheta*degree);
   
    rotmat->rotateZ(-dPsi*degree);
    rotmat->rotateX(-dTheta*degree);
   
    crystalmat->rotateZ(-dPsi*degree);
    crystalmat->rotateX(-dTheta*degree);
   
    //rotate position vector based on iClover.

    pos.rotateX(dTheta*degree);
    pos.rotateZ(dPsi *degree);

    case_pos.rotateX(dTheta *degree); 
    case_pos.rotateZ(dPsi*degree); 
    
    bgo_pos.rotateX(dTheta *degree); 
    bgo_pos.rotateZ(dPsi*degree); 
    
    bgowindow_pos.rotateX(dTheta *degree); 
    bgowindow_pos.rotateZ(dPsi*degree); 
        }
    }
    new G4PVPlacement( crystalmat,     // orientation matrix
                       pos,              // its position vector
                       fLogicCrystal[i],        // its logical volume                         
                       name,             // its name
                       worldLV,          // its mother  volume
                       false,            // no boolean operation
                       i,                // copy number
                       fCheckOverlaps);               // checking overlaps
     
    
  if( jCrystal  == 0 ) {  
      caseLV[iClover]  = new G4LogicalVolume(casing, G4Material::GetMaterial("G4_Al"), "Case");
      bgoLV[iClover]   = new G4LogicalVolume(bgo   , G4Material::GetMaterial("G4_Al"), "BGO ");
      bgowindowLV[iClover]   = new G4LogicalVolume(bgowindow, G4Material::GetMaterial("G4_Al"), "BGO WINDOW ");
     
      new G4PVPlacement( rotmat ,                // case orientation
                     case_pos,          // case location
                     caseLV[i/4],           // its logical volume                         
                     "Casing",         // its name
                     worldLV,          // its mother  volume
                     false,            // no boolean operation
                     iClover,                // copy number
                     fCheckOverlaps);  // checking overlaps
  
      new G4PVPlacement( bgomat ,      // bgo orientation
                     bgo_pos,          // bgo location
                     bgoLV[i/4],      // its logical volume                         
                     "Bgo",         // its name
                     worldLV,          // its mother  volume
                     false,            // no boolean operation
                     iClover,                // copy number
                     fCheckOverlaps);  // checking overlaps
      
      new G4PVPlacement( bgowindowmat ,      // bgo orientation
                     bgowindow_pos,          // bgo location
                     bgowindowLV[i/4],      // its logical volume                         
                     "Bgowindow",         // its name
                     worldLV,          // its mother  volume
                     false,            // no boolean operation
                     iClover,                // copy number
                     fCheckOverlaps);  // checking overlaps
                                       
  caseLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(1.0,1.0,0.0)));
  bgoLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(0.,1.0,0.0)));
  bgowindowLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(0.0,1.0,1.0)));
    }
       //Set a i+1%4 for each 
     //

  }

  // Al casing
  
  
  
  // Always return the physical World
  return worldPV;
}

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

void CloverDetectorConstruction::ConstructSDandField()
{

  G4cout << "********************* CloverDetectorConstruction::ConstructSDandField() " << G4endl;
  G4SDManager::GetSDMpointer()->SetVerboseLevel(1);

  // Sensitive detectors
  CloverCrystalSD * crystalSD = new CloverCrystalSD("Crystal", "CrystalHitsCollection", fNumOfCrystal);
  G4SDManager::GetSDMpointer()->AddNewDetector(crystalSD);

  //Set crystalSD to all logical volumn under the same name of "CrystalLV"
  SetSensitiveDetector("CrystalLV",crystalSD, true);

  // Create global magnetic field messenger.
  // Uniform magnetic field is then created automatically if
  // the field value is not zero.
  G4ThreeVector fieldValue;
  fMagFieldMessenger = new G4GlobalMagFieldMessenger(fieldValue);
  fMagFieldMessenger->SetVerboseLevel(1);
  
  // Register the field messenger for deleting
  G4AutoDelete::Register(fMagFieldMessenger);
}

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