DetectorConstruction.cc

// ********************************************************************
// * 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 /src/DetectorConstruction.cc
/// \brief Implementation of the DetectorConstruction class
//
//

// USER //
#include "DetectorConstruction.hh"
#include "DetectorMessenger.hh"
#include "Materials.hh"
#include "PMTSD.hh"

#include "CADMesh.hh"

#include <math.h>

// GEANT4 //
#include "G4Element.hh"
#include "G4SDManager.hh"
#include "G4RunManager.hh"
#include "G4GeometryManager.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4SolidStore.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4LogicalSkinSurface.hh"
#include "G4OpticalSurface.hh"
#include "G4ThreeVector.hh"
#include "G4PVPlacement.hh"
#include "G4SystemOfUnits.hh"
#include "G4Tubs.hh"
#include "G4Trd.hh"
#include "G4SubtractionSolid.hh"
#include "G4NistManager.hh"
#include "G4VisAttributes.hh"
#include "G4VisExtent.hh"
#include "G4Colour.hh"
#include "G4UImanager.hh"

/*
 *
 */
DetectorConstruction::DetectorConstruction()
  : G4VUserDetectorConstruction(),
    m_worldDim( new G4ThreeVector(1.0*m, 1.0*m, 1.0*m) ),
    m_LGpos( new G4ThreeVector() ),
    m_pmtPos( new G4ThreeVector(28.5*mm, 35.0*mm, -7.5*mm) ),
    m_rotation( new G4RotationMatrix() ),
    m_pmtDia(10.5*mm),
    m_PMTthickness(0.001*mm),
    m_ConstructionHasBeenDone(false),
    m_UsingCADmodel(true),
    m_logicWorld(0),
    m_physLG(0),
    m_physPMT(0),
    m_DetectorMessenger(nullptr)
{
  materials = Materials::getInstance();
  materials->UseOpticalMaterials(true);
  materials->DefineOpticalProperties();

  m_DetectorMessenger = new DetectorMessenger(this);
  m_runMan = G4RunManager::GetRunManager();

}


/*
 *
 */
DetectorConstruction::~DetectorConstruction(){
  delete m_DetectorMessenger;
}

/*
 * Execute geometry.mac here and construct the world according
 * to what it contains
 */
G4VPhysicalVolume* DetectorConstruction::Construct(){
  if ( m_ConstructionHasBeenDone ) {
    G4GeometryManager::GetInstance()->OpenGeometry();
    G4PhysicalVolumeStore::GetInstance()->Clean();
    G4LogicalVolumeStore::GetInstance()->Clean();
    G4SolidStore::GetInstance()->Clean();
    G4LogicalSkinSurface::CleanSurfaceTable();
    G4LogicalBorderSurface::CleanSurfaceTable();
  } else {
    G4UImanager* UImanager = G4UImanager::GetUIpointer();
    UImanager->ApplyCommand("/control/execute geometry.mac");
    BuildPMT();
  }

  if(!m_logicWorld) BuildWorld();
  if(!m_physLG) BuildDefaultLG();
  PlaceGeometry();

  m_ConstructionHasBeenDone = true;

  return m_physWorld;

}

/*
 *
 */
void DetectorConstruction::BuildWorld(){
  bool checkOverlaps = false;

  //----------------- Define the world volume -----------------//
      m_solidWorld =
	new G4Box("World",       //name
		  m_worldDim->x(),  //sizeX
		  m_worldDim->y(),  //sizeY
		  m_worldDim->z()); //sizeZ

            m_logicWorld =
	      new G4LogicalVolume(m_solidWorld, //solid
				  materials->Air,          //material
				  "World");     //name

	            m_physWorld =
		      new G4PVPlacement(0,                  //no rotation
					G4ThreeVector(),    //at (0,0,0)
					m_logicWorld,       //logical volume
					"World",            //name
					0,                  //mother  volume
					false,              //no boolean operation
					0,                  //copy number
					checkOverlaps);     //overlaps checking

		   // G4VisAttributes* boxVisAtt_world= new G4VisAttributes();     ***UNCOMMENT FOR VISIBLE WORLD VOLUME***
		   G4VisAttributes* boxVisAtt_world= new G4VisAttributes(G4VisAttributes::Invisible);
		   m_logicWorld ->SetVisAttributes(boxVisAtt_world);
}

/*
 * Make the default light guide
 */
void DetectorConstruction::BuildDefaultLG( ){

  // Build the light guide here up until the logical volume
  // Maybe build the aluminum portion into a single CAD solid and add
  // fibers and air here

  UseCADModel("../models/LED-GMS-light-guide.stl");

  bool CHECK_OVERLAPS = true;

  //Change these variables to what they really should be
  double fiberLength = 270*mm;
  double core_diam = 0.600*mm;
  double cladding_diam = 0.660*mm;
  double buffer_diam = 0.710*mm;

  //Core
  m_fiberCore = new G4Tubs("fiber_core", // name
			   0.0*mm,       // inner radius
			   core_diam/2.0,// outter radius
			   fiberLength,  // length
			   0.0*deg,      // sweep starting angle
			   360.0*deg);   // sweep stopping angle

  m_logicCore = new G4LogicalVolume( m_fiberCore,        // solid to be created from
				     materials->pQuartz, // material
				     "fiber_logicCore"); // name

  //Cladding
  m_fiberCladding = new G4Tubs("fiber_cladding",  // name
			       core_diam/2.0,     // inner radius
			       cladding_diam/2.0, // outter radius
			       fiberLength,       // length
			       0.0*deg,           // sweep starting angle
			       360.0*deg);        // sweep stopping angle

  m_logicCladding = new G4LogicalVolume( m_fiberCladding,    // solid to be created from
					 materials->pQuartz, // material
					 "fiber_logicCladding"); // name

  //Buffer
  m_fiberBuffer = new G4Tubs("fiber_buffer",   // name
			     cladding_diam/2.0,// inner radius
			     buffer_diam/2.0,  // outter radius
			     fiberLength,      // length
			     0.0*deg,          // sweep starting angle
			     360.0*deg);       // sweep stopping angle

  m_logicBuffer = new G4LogicalVolume( m_fiberBuffer,      // solid to be created from
				       materials->pQuartz, // material
				       "fiber_logicBuffer"); // name

  // Place the fibers
  G4RotationMatrix* fiberRotation = new G4RotationMatrix();
  fiberRotation->rotateX(90.*deg);
  char name[32];
  for(int i = 0; i < 32; i++){

    // calculate cooridnates for center of each fiber (use CAD model for placement):
    double x = 8.5 + 1.20*i;  // x-offset +  fiber spacing
    double y = 345;
    double z = -6.00;
    //Core
    sprintf(name,"fiberCore_%d", i);
    m_physCore.push_back(
			 new G4PVPlacement(fiberRotation,                  //Rotation
					   G4ThreeVector( x , y , z ), //Position
					   m_logicCore,                    //Logical volume this is built from
					   name,                           //Name
					   m_logicWorld,                   //Mother volume
					   false,                          //If there are boolean operations
					   i,                              //Copy number
					   CHECK_OVERLAPS) );              //Check for overlaps with other geometry


    //Cladding
    sprintf(name,"fiberCladding_%d", i);
    m_physCladding.push_back(
			     new G4PVPlacement(fiberRotation,                  //Rotation
					       G4ThreeVector( x , y , z ), //Position
					       m_logicCladding,                //Logical volume this is built from
					       name,                           //Name
					       m_logicWorld,                   //Mother volume
					       false,                          //If there are boolean operations
					       i,                              //Copy number
					       CHECK_OVERLAPS) );              //Check for overlaps with other geometry


    //Core
    sprintf(name,"fiberBuffer_%d", i);
    m_physBuffer.push_back(
			   new G4PVPlacement(fiberRotation,                  //Rotation
					     G4ThreeVector( x , y, z ), //Position
					     m_logicBuffer,                  //Logical volume this is built from
					     name,                           //Name
					     m_logicWorld,                   //Mother volume
					     false,                          //If there are boolean operations
					     i,                              //Copy number
					     CHECK_OVERLAPS) );              //Check for overlaps with other geometry

  }
}


/*
 * Make the PMT logical volume and Sensitvie Detector
 */
void DetectorConstruction::BuildPMT(){

  G4SDManager* SDman = G4SDManager::GetSDMpointer();
  m_SDnames.push_back("MyPMT");
  PMTSD* PMT = new PMTSD(m_SDnames.back(), false); // isIntermediate set to false. This is the final detector
  SDman->AddNewDetector( PMT );

      m_solidPMT =
	 new G4Tubs("PMT",         //name
		   0.0*mm,         //Inner radius
		   m_pmtDia/2.0,   //Outter radius
		   m_PMTthickness, //Height
		   0.0*deg,        //Rotation start
		   360.0*deg);     //Sweep


            m_logicPMT =
	      new G4LogicalVolume(m_solidPMT,       //solid
				  materials->Air,   //material
				  "PMT");           //name

	    G4VisAttributes* VisAtt_PMT = new G4VisAttributes(G4Colour(1.0,1.0,0.6,0.7));
	    m_logicPMT->SetVisAttributes(VisAtt_PMT);
	    m_logicPMT->SetSensitiveDetector( PMT );


      /***********************************
       * Intermediate Sensitive Detector 1
       **********************************/
      m_SDnames.push_back("IntSD1");
      PMT = new PMTSD(m_SDnames.back(), true); // isIntermediate set to true
      SDman->AddNewDetector( PMT );
      m_solidIntSD.push_back( new G4Box("soldIntSD1",20.00*mm, 0.001*mm, 2.00*mm) );
      m_logicIntSD.push_back( new G4LogicalVolume(m_solidIntSD.back(), materials->Air, "logicIntSD1") );
      m_logicIntSD.back()->SetVisAttributes(VisAtt_PMT);
      m_logicIntSD.back()->SetSensitiveDetector( PMT );

      /***********************************
       * Intermediate Sensitive Detector 2
       **********************************/
      m_SDnames.push_back("IntSD2");
      PMT = new PMTSD(m_SDnames.back(), true); // isIntermediate set to true
      SDman->AddNewDetector( PMT );
      m_solidIntSD.push_back( new G4Box("soldIntSD2",50.00*mm, 0.001*mm, 2.00*mm) );
      m_logicIntSD.push_back( new G4LogicalVolume(m_solidIntSD.back(), materials->Air, "logicIntSD2") );
      m_logicIntSD.back()->SetVisAttributes(VisAtt_PMT);
      m_logicIntSD.back()->SetSensitiveDetector( PMT );

}

/*
 * Place the geometry. This is done in it's own step in case
 * the user has specified a CAD model to be used
 */
void DetectorConstruction::PlaceGeometry(){

  //----------------- Place the LightGuide -----------------//

  m_physLG = new G4PVPlacement(m_rotation,
			       *m_LGpos,
			       m_logicLG,
			       "LightGuide",
			       m_logicWorld,
			       false,
			       0);

  //----------------- Place the PMT -----------------//

  m_physPMT = new G4PVPlacement(new G4RotationMatrix(),
				*m_pmtPos,
				m_logicPMT,
				"PMT",
				m_logicWorld,
				false,
				0,
				CHECK_OVERLAPS) );


  //----------------- Place Intermediate Detector 1 -----------------//

  m_physIntSD.push_back(
    new G4PVPlacement(new G4RotationMatrix(),
				              G4ThreeVector( 28.50*mm, 615.00*mm, -5.50*mm),
				              m_logicIntSD[0],
				              "IntSD1",
				              m_logicWorld,
				              false,
				              0,
					      CHECK_OVERLAPS) );
  //----------------- Place Intermediate Detector 2 -----------------//

  m_physIntSD.push_back(
    new G4PVPlacement(new G4RotationMatrix(),
                                              G4ThreeVector( 28.50*mm, 57.50*mm, -6.00*mm),
                                              m_logicIntSD[0],
                                              "IntSD2",
                                              m_logicWorld,
                                              false,
                                              0,
					      CHECK_OVERLAPS) );


  //----------------- Define Optical Borders -----------------//

  // One of these is correct. The order of the logical volumes is different
  // between the two. This is specific to the direction of light i.e. from
  // volume a to volume b or vice versa

  // m_Surface = new G4LogicalBorderSurface("AlSurface",
  //                              m_physLightGuide,
  //                              m_physWorld,
  //                              materials->AlSurface ) );
  //
  m_Surface = new G4LogicalBorderSurface("AlSurface",
					 m_physWorld,
					 m_physLG,
					 materials->AlSurface );

}

/*
 * Set the dimensions of the world volume.
 * If the volume exists modify it, otherwise build it.
 */
void DetectorConstruction::SetWorldVolume(G4ThreeVector arg){
  if(m_ConstructionHasBeenDone){
    m_solidWorld->SetXHalfLength( arg.x() );
    m_solidWorld->SetYHalfLength( arg.y() );
    m_solidWorld->SetZHalfLength( arg.z() );
    m_runMan->GeometryHasBeenModified();
    G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/rebuild");
  } else {
    delete m_worldDim;
    m_worldDim = new G4ThreeVector(arg);
    BuildWorld();
  }
}

/*
 * If the world already exists, rotate the light guide
 * otherwise set the rotation for the light guide that will be built
 */
void DetectorConstruction::SetRotation(G4ThreeVector arg){
  if(m_rotation) delete m_rotation;
  m_rotation = new G4RotationMatrix();

  m_rotation->rotateX(arg.x());
  m_rotation->rotateY(arg.y());
  m_rotation->rotateZ(arg.z());

  if(m_ConstructionHasBeenDone){
    m_physLG->SetRotation(m_rotation);
    m_runMan->GeometryHasBeenModified();
    G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/rebuild");
  }
}

/*
 * If the world already exists, move the light guide. If there is more
 * than one light guide, set the position and rebuild.
 * otherwise set the location for the light guide that will be built
 */
void DetectorConstruction::SetTranslation(G4ThreeVector arg){
  if(m_ConstructionHasBeenDone){
    m_physLG->SetTranslation(arg);
    m_runMan->GeometryHasBeenModified();
    G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/rebuild");
  } else {
    m_LGpos = new G4ThreeVector(arg);
  }
}


/*
 * If the world already exists, move the PMT. If there is more
 * than one PMT, set the position and rebuild.
 * otherwise set the location for the PMT that will be built
 */
void DetectorConstruction::SetPMTTranslation(G4ThreeVector arg){
  if(m_ConstructionHasBeenDone){
    m_physPMT->SetTranslation(arg);
    m_runMan->GeometryHasBeenModified();
    G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/rebuild");
  } else {
    m_pmtPos = new G4ThreeVector(arg);
  }
}

/*
 * If the world already exists, modify the PMT
 * otherwise set the diameter for PMT that will be built
 */
void DetectorConstruction::SetPMTDiameter(G4double arg){
  if(m_ConstructionHasBeenDone){
    m_solidPMT->SetOuterRadius(arg/2.0);
    m_runMan->GeometryHasBeenModified();
    G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/rebuild");
  } else {
    m_pmtDia = arg;
  }
}


/*
 * Builds a light guide logical volume from CAD model file
 */
void DetectorConstruction::UseCADModel(G4String fileName){
  //Delete anything we are going to redefine
  if(m_ConstructionHasBeenDone){
    delete m_logicLG;
    delete m_physLG;
    delete m_physPMT;
    delete m_Surface;
  }

  float xSize=0.0, ySize=0.0, zSize=0.0;

  G4String fileType = fileName.substr( fileName.last('.') + 1, fileName.size() - fileName.last('.'));

  if(fileType == "stl"){
    auto mesh = CADMesh::TessellatedMesh::FromSTL((char*) fileName.c_str());
    // auto mesh = CADMesh::TessellatedMesh::From((char*) fileName.c_str());
    mesh->SetScale(mm);
    mesh->SetOffset( G4ThreeVector(0, 0, 0) );
    mesh->SetReverse(false);

            m_logicLG =
	      new G4LogicalVolume(mesh->GetSolid(), //solid
				  materials->Al,           //material
				  "LightGuide");           //name
  }else{
    G4cout << "File type not supported" << G4endl;
  }


  if(m_logicLG !=0 ){

    //Print the extent to the console to help the user position the light guide
    char message[64];
    G4VisExtent extent = m_logicLG->GetSolid()->GetExtent();
    xSize = extent.GetXmax() - extent.GetXmin();
    ySize = extent.GetYmax() - extent.GetYmin();
    zSize = extent.GetZmax() - extent.GetZmin();
    G4cout << "===== Light guide extent =====" << G4endl;
    sprintf(message,"Xmin = %04.2f, Xmax = %04.2f: Width  = %04.2f", extent.GetXmin(), extent.GetXmax(), xSize );
    G4cout << message << G4endl;
    sprintf(message,"Ymin = %04.2f, Ymax = %04.2f: Height = %04.2f", extent.GetYmin(), extent.GetYmax(), ySize );
    G4cout << message << G4endl;
    sprintf(message,"Zmin = %04.2f, Zmax = %04.2f: Depth  = %04.2f", extent.GetZmin(), extent.GetZmax(), zSize );
    G4cout << message << G4endl;

  }

  if(m_pmtPos == 0) m_pmtPos = new G4ThreeVector( );

  if(m_ConstructionHasBeenDone){
    PlaceGeometry();
    m_runMan->GeometryHasBeenModified();
    G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/rebuild");
  }
  G4cout << "Using " << fileName << " for CAD model" << G4endl;
}


/*
 * Set the surface model for the light guide
 */
void DetectorConstruction::SetSurfaceModel(const G4OpticalSurfaceModel model){
  materials->AlSurface->SetModel(model);
  if(m_ConstructionHasBeenDone){
    m_runMan->GeometryHasBeenModified();
  }
}

/*
 *
 */
void DetectorConstruction::SetSurfaceFinish(const G4OpticalSurfaceFinish finish){
  materials->AlSurface->SetFinish(finish);
  if(m_ConstructionHasBeenDone){
    m_runMan->GeometryHasBeenModified();
  }
}

/*
 *
 */
void DetectorConstruction::SetSurfaceType(const G4SurfaceType type){
  materials->AlSurface->SetType(type);
  if(m_ConstructionHasBeenDone){
    m_runMan->GeometryHasBeenModified();
  }
}

/*
 *
 */
void DetectorConstruction::SetSurfaceSigmaAlpha(G4double v){
  materials->AlSurface->SetSigmaAlpha(v);

  if(m_ConstructionHasBeenDone){
    G4RunManager::GetRunManager()->GeometryHasBeenModified();
  }

  G4cout << "Surface sigma alpha set to: " << materials->AlSurface->GetSigmaAlpha() << G4endl;
}

/*
 * Set the material property table of the surface
 * of the light guide
 */
void DetectorConstruction::AddSurfaceMPV(const char* c, G4MaterialPropertyVector* mpv){
  mpv->SetSpline(true);
  materials->GetMPTArray().at(1)->AddProperty(c, mpv);
  materials->AlSurface->SetMaterialPropertiesTable(materials->GetMPTArray().at(1));
  G4cout << "The MPT for the surface is now: " << G4endl;
  materials->GetMPTArray().at(1)->DumpTable();
  G4cout << "............." << G4endl;
}