/// \file DetectorConstruction.cc
/// \brief Implementation of the DetectorConstruction class
#include “DetectorConstruction.hh”
#include “DetectorMessenger.hh”
#include “G4Material.hh”
#include “G4NistManager.hh”
#include “G4GeometryManager.hh”
#include “G4PhysicalVolumeStore.hh”
#include “G4LogicalVolumeStore.hh”
#include “G4SolidStore.hh”
#ifdef G4MULTITHREADED
#include “G4MTRunManager.hh”
#else
#include “G4RunManager.hh”
#endif
#include “G4Box.hh”
#include “G4Cons.hh”
#include “G4Orb.hh”
#include “G4Sphere.hh”
#include “G4Trd.hh”
#include “G4Tubs.hh”
#include “G4LogicalVolume.hh”
#include “G4PVPlacement.hh”
#include “G4UnionSolid.hh”
#include “G4SystemOfUnits.hh”
#include “G4PhysicalConstants.hh”
#include “G4UnitsTable.hh”
#include “G4Isotope.hh”
#include “G4Element.hh”
//includes for sensitive detectors
#include “G4SDManager.hh”
#include “SensitiveDetector.hh”
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//have to define it as a null pointer that is just for the local thread or else multithreading won’t work
G4ThreadLocal SensitiveDetector* DetectorConstruction::SD = nullptr;
DetectorConstruction::DetectorConstruction()
: G4VUserDetectorConstruction(),fDetectorMessenger(nullptr),emitter_material(nullptr),collector_material(nullptr),wire_material(nullptr),world_material(nullptr),insulator_material(nullptr),newMat(nullptr),insulator_newMat(nullptr),collector_newMat(nullptr),wire_newMat(nullptr),world_newMat(nullptr),rhodium(nullptr),inconel(nullptr),SS304L(nullptr),MgO(nullptr),Al2O3(nullptr),air(nullptr),world_logical(nullptr),insulator_logic(nullptr),emitter_logic(nullptr),wire_logic(nullptr),collector_logic(nullptr),world_physical(nullptr),vanadium(nullptr),platinum(nullptr),cobalt(nullptr),HfO2(nullptr)
{
////////// INITIALIZE EMIITER, INSULATOR AND COLLECTOR WITH SOME MATERIALS //////////////////
//THESE CAN BE CHANGED IN MACRO FILES ON THE FLY
materialName = "rhodium"; //emitter material
insulator_materialName = "Al2O3"; //insulator material
collector_materialName = "SS304L"; //colletor material
wire_materialName = "SS304L";//wire material
world_materialName = "G4_AIR";//world material
//////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////// INITIALIZE DETECTOR DIMENSIONS //////////////////////////////////////////
//emitter dimensions
emitter_radius=0.23*mm; //Radius of emitter
emitter_length= 2*mm;//400*mm; //Length of emitter
//SPND/collector dimensions
//G4double detector_radius=0.97*mm; //Radius from center of detector to outside the inconel
collector_thickness=0.25*mm; //thickness of the inconel collector on the sides
//wire dimensions
wire_radius=0.0125*mm; //Radius of wire connected to emitter
// scale factor for how big want the diameter of the bottom part of the SPND to be relative to the detector
scale=0.5;
//insulator dimensions
insulator_thickness=0.5*mm; //thickness insulator on the top and bottom of emitter (in between emitter and endcap)
// define how much bigger that the detector the world volume should be
world_buffer=10*mm; //space in between detector and end of "world" on each side
//G4double world_radius=63.5*mm; //assume pipe is 2.5"
detector_radius=emitter_radius+insulator_thickness+collector_thickness; //Radius from center of detector to outside the inconel
world_length = emitter_length + (insulator_thickness*2) + detector_radius; // world is the length of the emitter + length of the insulator or either side of emitter + length of end cap + buffer, are assuming detector sits on bottom of world
world_radius = detector_radius + world_buffer; //world radius is the radius of the detector + buffer
//INSULATOR
//Define dimensions of insulator
//inner radii
insulator_bottom_radiusIn=wire_radius; //inner radius below emitter is same as wire
insulator_middle_radiusIn=emitter_radius; //inner radius around emitter is same as emitter
//outer radius
insulator_radiusOut=detector_radius-collector_thickness; //outer radius of insulator
//lengths of insulator - these must be half lengths because counts from middle
insulator_bottom_length=insulator_thickness/2;
insulator_top_length=insulator_thickness/2;
insulator_middle_length=emitter_length/2;
//COLLECTOR
collector_side_length=insulator_middle_length+insulator_bottom_length+insulator_top_length; // length of middle part of collector
collector_inner_radius=detector_radius-collector_thickness; // inner radius of collector
collector_bottom_length = detector_radius/2; //half length of cone on bottom of collector
cone_inner_radius =collector_inner_radius*scale; // inner radius of code on bottom of collector
cone_outer_radius =detector_radius*scale; //outer radius of bottom of cone
//WIRE
wire_length=insulator_bottom_length+collector_bottom_length;
/////////////////////////////////////////////////////////
//These are the dimensions that need to scale with what input for the definitions
//INSULATOR
//Define dimensions of insulator
//inner radii
insulator_bottom_radiusIn=wire_radius; //inner radius below emitter is same as wire
insulator_middle_radiusIn=emitter_radius; //inner radius around emitter is same as emitter
//outer radius
insulator_radiusOut=detector_radius-collector_thickness; //outer radius of insulator
//lengths of insulator - these must be half lengths because counts from middle
insulator_bottom_length=insulator_thickness/2;
insulator_top_length=insulator_thickness/2;
insulator_middle_length=emitter_length/2;
//COLLECTOR
collector_side_length=insulator_middle_length+insulator_bottom_length+insulator_top_length; // length of middle part of collector
collector_inner_radius=detector_radius-collector_thickness; // inner radius of collector
collector_bottom_length = detector_radius/2; //half length of cone on bottom of collector
cone_inner_radius =collector_inner_radius*scale; // inner radius of code on bottom of collector
cone_outer_radius =detector_radius*scale; //outer radius of bottom of cone
//WIRE
wire_length=insulator_bottom_length+collector_bottom_length;
////////////////////////////////////////////////////////////////////////////////////////////
// Create commands for interactive definition of the geometry
fDetectorMessenger = new DetectorMessenger(this);
}
//////////////////////////////////////////////////////////////////////////////////////
DetectorConstruction::~DetectorConstruction()
{delete fDetectorMessenger;
//delete SD;
//SD = nullptr;
}
///////////////////////////////////////////////////////////////////////////////////////
G4VPhysicalVolume* DetectorConstruction::Construct()
{
return ConstructVolumes();
}
///////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////
///////////// DEFINE SOME USEFUL METHODS ///////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////// CONSTRUCT VOLUMES METHOD ////////////////////////////////////////
G4VPhysicalVolume* DetectorConstruction::ConstructVolumes()
{
// Cleanup old geometry
G4GeometryManager::GetInstance()->OpenGeometry();
G4PhysicalVolumeStore::GetInstance()->Clean();
G4LogicalVolumeStore::GetInstance()->Clean();
G4SolidStore::GetInstance()->Clean();
// Get nist material manager
//G4NistManager* nist = G4NistManager::Instance();
//
// Option to switch on/off checking of volumes overlaps
//
checkOverlaps = true;
/////////////////////////// DEFINE MATIERALS ///////////////////////////////////
Most materials left out because of space limit but here is an example:
//////////////////////// PLATINUM //////////////////////////
G4Isotope* Pt190 =
new G4Isotope(“190Pt”, //name
Z=78, // proton number (Z)
A=190, // nucleon number (A)
atomic_weight=189.96*g/mole // atomic weight
);
G4Isotope* Pt192 =
new G4Isotope("192Pt", //name
Z=78, // proton number (Z)
A=192, // nucleon number (A)
atomic_weight=191.96*g/mole // atomic weight
);
G4Isotope* Pt194 =
new G4Isotope("194Pt", //name
Z=78, // proton number (Z)
A=194, // nucleon number (A)
atomic_weight=193.96*g/mole // atomic weight
);
G4Isotope* Pt195 =
new G4Isotope("195Pt", //name
Z=78, // proton number (Z)
A=195, // nucleon number (A)
atomic_weight=194.96*g/mole // atomic weight
);
G4Isotope* Pt196 =
new G4Isotope("196Pt", //name
Z=78, // proton number (Z)
A=196, // nucleon number (A)
atomic_weight=195.96*g/mole // atomic weight
);
G4Isotope* Pt198 =
new G4Isotope("198Pt", //name
Z=78, // proton number (Z)
A=198, // nucleon number (A)
atomic_weight=197.96*g/mole // atomic weight
);
//Define element and add isotopes
G4Element* Pt_natural =
new G4Element("Pt_natural", //element name
"Pt", // element symbol
nIso=6 // # of isotopes
);
Pt_natural -> AddIsotope(Pt190,abundance=0.014*perCent);
Pt_natural -> AddIsotope(Pt192,abundance=0.782*perCent);
Pt_natural -> AddIsotope(Pt194,abundance=32.967*perCent);
Pt_natural -> AddIsotope(Pt195,abundance=33.832*perCent);
Pt_natural -> AddIsotope(Pt196,abundance=25.242*perCent);
Pt_natural -> AddIsotope(Pt198,abundance=7.163*perCent);
G4double density_platinum = 21.45*g/cm3;
//Build emitter matieral
platinum =
new G4Material("platinum", //name
density_platinum, //density
1, //number of components
kStateSolid, //state
temperature //temperature
);
platinum -> AddElement(Pt_natural,mass_fraction=100.*perCent);
/* Mass Number Natural Abundance Atomic Weight
190 0.014% 189.959931655
192 0.782% 191.961038005
194 32.967% 193.962680253
195 33.832% 194.964791134
196 25.242% 195.964951521
198 7.163% 197.96789279
*/
// Print materials table on startup
G4cout << *(G4Material::GetMaterialTable()) << G4endl;
///////////////////////////////// DETECTOR MATERIALS /////////////////////////////////////////
//Actually assign which materials want to use for emitter, collector, wire, insulator and world here
SetEmitterMaterial(materialName);
SetInsulatorMaterial(insulator_materialName);
SetCollectorMaterial(collector_materialName);
SetWireMaterial(wire_materialName);
SetWorldMaterial(world_materialName);
/////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////// WORLD VOLUME ///////////////////////////////////////////////
//Make world solid (where define dimensions of world)
G4VSolid* world_solid =
new G4Tubs("World", //name
0 * cm, //inner radius, world is not hollow obviously
world_radius, //outer radius
world_length, //height of cylinder
0.0 * deg, 360.0 * deg); // segment angles
//Make logical volume for world, define what world is made of here
world_logical= new G4LogicalVolume(
world_solid, //its solid defining dimensions of world
world_material, //what the world is made of
"World"); //its name
//Make physical volume for world, actually place it and make it the "mother" volume
world_physical= new G4PVPlacement(
0, // no rotation
G4ThreeVector(), // at (0,0,0)
world_logical, // its logical volume
"World", // its name
0, // it is the mother volume
false, // no boolean operations
0, // copy number
checkOverlaps); // checking overlaps
///////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////INSULATOR ///////////////////////////////////////////////
// Insulator is in four parts, top, which is solid, the middle which is hollow to make room for the emitter, the bottom which is also hollow but has a smaller diameter to allow the wire leading to the emitter to pass through and cone which is the most bottom part in the cone, also is hollow to allow wire to pass through
//BOTTOM
G4ThreeVector insulator_bottom_pos = G4ThreeVector(0, 0,-(insulator_middle_length+insulator_bottom_length)); //where the bottom is located relative to middle
//describe the solid
G4VSolid* insulator_bottom_solid =
new G4Tubs("insulator_bottom", //name
insulator_bottom_radiusIn, //inner radius, is hollow cylinder
insulator_radiusOut, //outer radius
insulator_bottom_length, //height of cylinder, is a half length
0.0*deg, 360.0*deg // segment angles
);
// MIDDLE
//describe the solid
G4VSolid* insulator_middle_solid =
new G4Tubs("insulator_middle", //name
insulator_middle_radiusIn, //inner radius, is hollow cylinder
insulator_radiusOut, //outer radius
insulator_middle_length, //height of cylinder
0.0*deg, 360.0*deg // segment angles
);
//TOP
//where the top is located with respect to the middle (need this to build the it)
G4ThreeVector insulator_top_pos = G4ThreeVector(0, 0,insulator_middle_length+insulator_bottom_length);
//Describe the solid
G4VSolid* insulator_top_solid =
new G4Tubs("insulator_top", //name
0*mm, //inner radius
insulator_radiusOut, //outer radius
insulator_top_length, //height of cylinder
0.0*deg, 360.0*deg // segment angles
);
//CONE
//where the cone is located with respect to the middle (need this to build the it)
G4ThreeVector insulator_cone_pos = G4ThreeVector(0, 0,-(insulator_middle_length+2*insulator_bottom_length+collector_bottom_length));
//Describe the solid
G4VSolid* insulator_cone_solid =
new G4Tubs("insulator_cone", //name
insulator_bottom_radiusIn, //inner radius
cone_inner_radius, //outer radius
collector_bottom_length, //height of cylinder
0.0*deg, 360.0*deg // segment angles
);
// PUT THE INSULATOR TOGETHER
// Put two parts together
G4VSolid* insulator_temp
= new G4UnionSolid(
"insulator_temp", //name
insulator_middle_solid, //solid1
insulator_top_solid, //solid2
0, //rotation
insulator_top_pos //position of solid 2 relative to 1
);
//Put two more together
G4VSolid* insulator_temp2
= new G4UnionSolid(
"insulator_temp2", //name
insulator_temp, //solid1
insulator_bottom_solid, //solid2
0, //rotation
insulator_bottom_pos //position of solid 2 relative to 1
);
//Finish building insulator
G4VSolid* insulator_solid
= new G4UnionSolid(
"insulator", //name
insulator_temp2, //solid1
insulator_cone_solid, //solid2
0, //rotation
insulator_cone_pos //position of solid 2 relative to 1
);
// Describe what the insulator is
insulator_logic = new G4LogicalVolume(
insulator_solid, //its solid
insulator_material, //its material
"insulator"); //its name
// PLACE THE INSULATOR IN THE WORLD VOLUME
// Location want the insulator
G4ThreeVector insulator_pos = G4ThreeVector(0, 0, 0);
G4VPhysicalVolume* insulator_physical =
new G4PVPlacement(0, //no rotation
insulator_pos, //at position
insulator_logic, //its logical volume
"insulator", //its name
world_logical, //its mother volume
false, //no boolean operation
0, //copy number
checkOverlaps); //overlaps checking
////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////// EMITTER ////////////////////////////////////
//describe the solid
G4VSolid* emitter_solid =
new G4Tubs("emitter", //name
0, //inner radius, is hollow cylinder
emitter_radius, //outer radius
emitter_length/2, //height of cylinder, is a half length
0.0*deg, 360.0*deg // segment angles
);
// Describe what the emitter is
emitter_logic = new G4LogicalVolume(
emitter_solid, //its solid
emitter_material, //its material
"emitter"); //its name
// Location want the emitter
G4ThreeVector emitter_pos = G4ThreeVector(0, 0, 0);
//Place emitter in world volume
G4VPhysicalVolume* emitter_physical =
new G4PVPlacement(0, //no rotation
emitter_pos, //at position
emitter_logic, //its logical volume
"emitter", //its name
world_logical, //its mother volume
false, //no boolean operation
0, //copy number
checkOverlaps); //overlaps checking
/////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////// WIRE //////////////////////////////////////////////////////
//describe the solid
G4VSolid* wire_solid =
new G4Tubs("wire", //name
0, //inner radius, is hollow cylinder
wire_radius, //outer radius
wire_length, //height of cylinder
0.0*deg, 360.0*deg // segment angles
);
// Describe what the wire is
wire_logic = new G4LogicalVolume(
wire_solid, //its solid
wire_material, //its material
"wire"); //its name
// Location want the wire
G4ThreeVector wire_pos = G4ThreeVector(0, 0,-(insulator_middle_length+wire_length));
//Place wire in world volume
G4VPhysicalVolume* wire_physical =
new G4PVPlacement(0, //no rotation
wire_pos, //at position
wire_logic, //its logical volume
"wire", //its name
world_logical, //its mother volume
false, //no boolean operation
0, //copy number
checkOverlaps); //overlaps checking
/////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////// COLLECTOR ////////////////////////////////////////////
//Collector is composed of 3 parts, the endcap which is a hemisphere, the "side" which is the part that runs down the outside of the emitter and the cone which is a cone at the bottom of the collector
// Location want the endcap relative to collector side
G4ThreeVector endcap_pos = G4ThreeVector(0,0,collector_side_length);
//Top of collector is a hemispherical cap
G4VSolid* endcap_solid =
new G4Sphere("endcap", //name
0*mm, //inner radius, =0 because is solid
detector_radius, // outer radius
0.0*deg, 360*deg, // segment angles phi
0.0*deg, 90*deg // segment angles theta
);
//Side of collector
G4VSolid* collector_side_solid =
new G4Tubs("collector_side", //name
collector_inner_radius, //inner radius
detector_radius, //outer radius
collector_side_length, //height of cylinder
0.0*deg, 360.0*deg // segment angles
);
// Location of the bottom cone relative to collector side
G4ThreeVector collector_bottom_pos = G4ThreeVector(0,0,-(collector_side_length+collector_bottom_length));
// Generate a rotation matrix because need to flip cone 180 degrees
G4RotationMatrix* collector_bottom_rotation = new G4RotationMatrix;
collector_bottom_rotation -> rotateX(180*deg); //rotate 180 degrees in x, leave y and z alone
//Bottom of collector is a cone for now
G4VSolid* collector_bottom_solid =
new G4Cons( "collector_bottom",
cone_inner_radius, //inner radius top
detector_radius, //outer radius top
cone_inner_radius, //inner radius bottom
cone_outer_radius, //outer radius bottom
collector_bottom_length, // height of cone
0.0* deg, 360*deg // segment angles
);
// PUT THE COLLECTOR TOGETHER
// Put two parts together
G4VSolid* collector_temp
= new G4UnionSolid(
"collector_temp", //name
collector_side_solid, //solid1
endcap_solid, //solid2
0, //rotation
endcap_pos //position of solid 2 relative to 1
);
// Put rest together
G4VSolid* collector_solid
= new G4UnionSolid(
"collector", //name
collector_temp, //solid1
collector_bottom_solid, //solid2
collector_bottom_rotation, //rotation
collector_bottom_pos //position of solid 2 relative to 1
);
// Describe what the collector is
collector_logic = new G4LogicalVolume(
collector_solid, //its solid
collector_material, //its material
"collector"); //its name
// Location want the collector
G4ThreeVector collector_pos = G4ThreeVector(0,0,0);
//Place collector in world volume
G4VPhysicalVolume* collector_physical =
new G4PVPlacement(0, //no rotation
collector_pos, //at position
collector_logic, //its logical volume
"collector", //its name
world_logical, //its mother volume
false, //no boolean operation
0, //copy number
checkOverlaps); //overlaps checking
/////////////////////////////////////////////////////////////////////////////////////////////////
// Print what emitter material is on start up
G4cout << "\n--> Emitter Material Set to : "
<< emitter_logic->GetMaterial()->GetName() << G4endl;
//
//always return the physical World
//
return world_physical;
}
//////////////////// END CONSTRUCT VOLUMES METHOD ////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////// SET SENSITIVE DETECTOR METHOD ///////////////////////////////////
void DetectorConstruction::ConstructSDandField()
{
//need to check and see if the sensitive detector has already been delcared (for when re-initialize geometry)
if (SD == nullptr)
{
SD = new SensitiveDetector(“SD”, “DetectorHitsCollection”);// analysis);
G4SDManager::GetSDMpointer()->AddNewDetector(SD);
SetSensitiveDetector("emitter", SD);
SetSensitiveDetector("insulator", SD);
SetSensitiveDetector("wire", SD);
SetSensitiveDetector("collector", SD);
SetSensitiveDetector("World", SD);
}
/*
SetSensitiveDetector("emitter", SD);
SetSensitiveDetector("insulator", SD);
SetSensitiveDetector("wire", SD);
SetSensitiveDetector("collector", SD);
SetSensitiveDetector("World", SD);*/
}
////////////////////////////// END SET SENSITIVE DETECTOR METHOD ////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////
Set Detector parameter methods with some removed for because of space limit, but they are all virtually identical
/////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////// SET EMITTER MATERIAL METHOD //////////////////////////////
void DetectorConstruction::SetEmitterMaterial(G4String materialChoice)
{
//DetectorConstruction::SD = nullptr;
materialName = materialChoice;
newMat = G4Material::GetMaterial(materialChoice);
if(newMat)
{
emitter_material = newMat;
if(emitter_logic) { emitter_logic->SetMaterial(emitter_material); }
#ifdef G4MULTITHREADED
G4MTRunManager::GetRunManager()->ReinitializeGeometry();
G4MTRunManager::GetRunManager()->PhysicsHasBeenModified();
#else
G4RunManager::GetRunManager()->ReinitializeGeometry();
G4RunManager::GetRunManager()->PhysicsHasBeenModified();
#endif
}
else
{
G4cout << “\n–> warning from DetectorConstruction::SetMaterial : "
<< materialChoice << " not found” << G4endl;
}
G4cout << "\n--> Emitter Material Set to: "
<< emitter_material->GetName() << G4endl;
}
/////////////////////////// END SET EMITTER MATERIAL METHOD ////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////// SET EMITTER RADIUS METHOD //////////////////////////////
void DetectorConstruction::SetEmitterRadius(G4double value)
{
//DetectorConstruction::SD = nullptr;
// Set new emitter radius
emitter_radius = value;
//Recalculate every calculated detector parameter
detector_radius=emitter_radius+insulator_thickness+collector_thickness; //Radius from center of detector to outside the inconel
world_length = emitter_length + (insulator_thickness*2) + detector_radius; // world is the length of the emitter + length of the insulator or either side of emitter + length of end cap + buffer, are assuming detector sits on bottom of world
world_radius = detector_radius + world_buffer; //world radius is the radius of the detector + buffer
insulator_bottom_radiusIn=wire_radius; //inner radius below emitter is same as wire
insulator_middle_radiusIn=emitter_radius; //inner radius around emitter is same as emitter
insulator_radiusOut=detector_radius-collector_thickness; //outer radius of insulator
insulator_bottom_length=insulator_thickness/2;
insulator_top_length=insulator_thickness/2;
insulator_middle_length=emitter_length/2;
collector_side_length=insulator_middle_length+insulator_bottom_length+insulator_top_length; // length of middle part of collector
collector_inner_radius=detector_radius-collector_thickness; // inner radius of collector
collector_bottom_length = detector_radius/2; //half length of cone on bottom of collector
cone_inner_radius =collector_inner_radius*scale; // inner radius of code on bottom of collector
cone_outer_radius =detector_radius*scale; //outer radius of bottom of cone
wire_length=insulator_bottom_length+collector_bottom_length;
//Re-initialize the geometry
#ifdef G4MULTITHREADED
G4MTRunManager::GetRunManager()->ReinitializeGeometry();
#else
G4RunManager::GetRunManager()->ReinitializeGeometry();
#endif
//Write to the output confirming changes
G4cout << "\n Emitter radius set to: " << G4BestUnit(emitter_radius,"Length")
<< "and is made of " << emitter_material->GetName()
<< "\n \n" << emitter_material << G4endl;
}
/////////////////////////// END SET EMITTER RADIUS METHOD ////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////