lars/minicalosim
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

construction wrapper seems to work

Browse Source
This commit is contained in:
Jan Kieseler
2023-09-07 23:31:46 +02:00
parent 3521cfd79e
commit 0adb64f763
6 changed files with 237 additions and 109 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CMakeLists.txt
+10 -5
View File
@@ -4,6 +4,8 @@
cmake_minimum_required(VERSION 3.16...3.21)
project(B4a)
#----------------------------------------------------------------------------
# Find Geant4 package, activating all available UI and Vis drivers by default
# You can set WITH_GEANT4_UIVIS to OFF via the command line or ccmake/cmake-gui
@@ -36,6 +38,13 @@ file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.hh)
add_executable(exampleB4a exampleB4a.cc ${sources} ${headers})
target_link_libraries(exampleB4a ${Geant4_LIBRARIES})
## needs to be added later
#file(GLOB_RECURSE SOURCES "${PROJECT_SOURCE_DIR}/src/*.cc" )
#add_subdirectory(lib/pybind11)
#pybind11_add_module(compiled ${SOURCES} "${PROJECT_SOURCE_DIR}/bind/bindings.cpp")
#----------------------------------------------------------------------------
# Copy all scripts to the build directory, i.e. the directory in which we
# build B4a. This is so that we can run the executable directly because it
@@ -61,11 +70,7 @@ foreach(_script ${EXAMPLEB4A_SCRIPTS})
)
endforeach()
file(GLOB_RECURSE SOURCES "${PROJECT_SOURCE_DIR}/src/*.cc" )
add_subdirectory(lib/pybind11)
pybind11_add_module(compiled ${SOURCES} "${PROJECT_SOURCE_DIR}/bind/bindings.cpp")
#----------------------------------------------------------------------------
# Install the executable to 'bin' directory under CMAKE_INSTALL_PREFIX
#
install(TARGETS exampleB4a DESTINATION bin)
#install(TARGETS exampleB4a DESTINATION bin)
exampleB4a.cc
+11 -1
View File
@@ -112,9 +112,19 @@ int main(int argc,char** argv)
}
#endif
ConstructionWrapper cw;
cw.addLayer(1, "G4_Pb", false, 2, 2);
cw.addLayer(10, "G4_Si", true, 2, 2);
cw.addLayer(1, "G4_Pb", false, 2, 2);
cw.addLayer(20, "G4_Si", true, 2, 2);
cw.addLayer(4, "G4_Pb", false, 2, 2);
cw.addLayer(3, "G4_Si", true, 2, 2);
// Set mandatory initialization classes
//
auto detConstruction = new B4::DetectorConstruction();
auto detConstruction = new B4::DetectorConstruction(cw);
// FOR LATER: assign constructionwrapper here (will be passed from other function)
runManager->SetUserInitialization(detConstruction);
auto physicsList = new FTFP_BERT;
include/ConstructionWrapper.hh
+48
View File
@@ -0,0 +1,48 @@
#include <vector>
#include <string>
#include "G4VPhysicalVolume.hh"
class Layer{
public:
Layer();
~Layer(){};
void setThickness(double thickness_cm);
void setMaterial(std::string material);
void setNx(int nx);
void setNy(int ny);
void setIsActive(bool isActive);
void assignPhysicalVolume(G4VPhysicalVolume* physicalVolume);
double thickness;
double sens_xwidth;
double sens_ywidth;
std::string material;
int nx;
int ny;
bool isActive;
G4VPhysicalVolume* physicalVolume;
std::string name;
};
class ConstructionWrapper{
public:
ConstructionWrapper(double xy_width=50);
~ConstructionWrapper() {};
void addLayer(double thickness_cm, std::string material, bool isActive=true, int nx=1, int ny=1);
std::vector<Layer>& getLayers() ;
const std::vector<Layer>& getLayers() const;
double getXYWidth() const{
return xywidth;
}
private:
double xywidth;
std::vector<Layer> layers;
};
include/DetectorConstruction.hh
+8 -1
View File
@@ -32,6 +32,7 @@
#include "G4VUserDetectorConstruction.hh"
#include "globals.hh"
#include "ConstructionWrapper.hh"
class G4VPhysicalVolume;
class G4GlobalMagFieldMessenger;
@@ -53,10 +54,14 @@ namespace B4
/// In addition a transverse uniform magnetic field is defined
/// via G4GlobalMagFieldMessenger class.
class DetectorConstruction : public G4VUserDetectorConstruction
{
public:
DetectorConstruction() = default;
DetectorConstruction(const ConstructionWrapper& cw){
this->cw = cw;
}
~DetectorConstruction() override = default;
public:
@@ -83,6 +88,8 @@ class DetectorConstruction : public G4VUserDetectorConstruction
G4VPhysicalVolume* fGapPV = nullptr; // the gap physical volume
G4bool fCheckOverlaps = true; // option to activate checking of volumes overlaps
ConstructionWrapper cw;
};
// inline functions
src/ConstructionWrapper.cc
+62
View File
@@ -0,0 +1,62 @@
#include "ConstructionWrapper.hh"
Layer::Layer(){
thickness = 0;
material = "";
nx = 1;
ny = 1;
isActive = false;
physicalVolume = nullptr;
}
void Layer::setThickness(double thickness){
this->thickness = thickness;
}
void Layer::setMaterial(std::string material){
this->material = material;
}
void Layer::setNx(int nx){
this->nx = nx;
}
void Layer::setNy(int ny){
this->ny = ny;
}
void Layer::setIsActive(bool isActive){
this->isActive = isActive;
}
void Layer::assignPhysicalVolume(G4VPhysicalVolume* physicalVolume){
this->physicalVolume = physicalVolume;
}
ConstructionWrapper::ConstructionWrapper(double xy_width){
xywidth = xy_width;
}
void ConstructionWrapper::addLayer(double thickness, std::string material, bool isActive, int nx, int ny){
Layer layer;
layer.setThickness(thickness);
layer.setMaterial(material);
layer.setIsActive(isActive);
if(!isActive){
nx = 1;
ny = 1;
}
layer.setNx(nx);
layer.setNy(ny);
layer.sens_xwidth = xywidth/(float)nx;
layer.sens_ywidth = xywidth/(float)ny;
layers.push_back(layer);
}
const std::vector<Layer> & ConstructionWrapper::getLayers() const{
return layers;
}
std::vector<Layer> & ConstructionWrapper::getLayers(){
return layers;
}
src/DetectorConstruction.cc
+96 -100
View File
@@ -49,10 +49,48 @@
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "G4PVParameterised.hh"
#include "G4VPVParameterisation.hh"
namespace B4
{
//helper
class LayerParametrisation: public G4VPVParameterisation{
public:
LayerParametrisation(Layer layer, G4double position = 0): G4VPVParameterisation(), layer(layer), position(position){
}
~LayerParametrisation() = default;
void ComputeTransformation(const G4int copyNo, G4VPhysicalVolume* physVol) const{
G4double x = 0;
G4double y = 0;
G4double z = 0;
if(layer.nx > 1){
x = (copyNo % layer.nx) * layer.sens_xwidth*cm;
}
if(layer.ny > 1){
y = (copyNo / layer.nx) * layer.sens_ywidth*cm;
}
G4ThreeVector origin(x, y, z);
origin -= G4ThreeVector(layer.sens_xwidth * ((float)layer.nx-1) / 2.*cm, layer.sens_ywidth * ((float)layer.ny-1) / 2.*cm, 0);
origin += G4ThreeVector(0., 0., position);
physVol->SetTranslation(origin);
}
void ComputeDimensions(G4Box& box, const G4int copyNo, const G4VPhysicalVolume* physVol) const{
box.SetXHalfLength(layer.sens_xwidth/2*cm);
box.SetYHalfLength(layer.sens_ywidth/2*cm);
box.SetZHalfLength(layer.thickness/2*cm);
}
private:
Layer layer;
G4double position;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4ThreadLocal
@@ -75,18 +113,11 @@ void DetectorConstruction::DefineMaterials()
{
// Lead material defined using NIST Manager
auto nistManager = G4NistManager::Instance();
nistManager->FindOrBuildMaterial("G4_Pb");
// Liquid argon material
G4double a; // mass of a mole;
G4double z; // z=mean number of protons;
G4double density;
new G4Material("liquidArgon", z=18., a= 39.95*g/mole, density= 1.390*g/cm3);
// The argon by NIST Manager is a gas with a different density
// Vacuum
new G4Material("Galactic", z=1., a=1.01*g/mole,density= universe_mean_density,
kStateGas, 2.73*kelvin, 3.e-18*pascal);
nistManager->FindOrBuildMaterial("G4_AIR");
auto cwLayers = cw.getLayers();
for(auto layer : cwLayers){
nistManager->FindOrBuildMaterial(layer.material);
}
// Print materials
G4cout << *(G4Material::GetMaterialTable()) << G4endl;
@@ -97,28 +128,21 @@ void DetectorConstruction::DefineMaterials()
G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
{
// Geometry parameters
G4int nofLayers = 10;
G4double absoThickness = 10.*mm;
G4double gapThickness = 5.*mm;
G4double calorSizeXY = 10.*cm;
auto & cwLayers = cw.getLayers();
G4int nofLayers = cwLayers.size();
G4double caloLength = 0;
for(auto layer : cwLayers){
caloLength += layer.thickness * cm;
}
G4double calorSizeXY = cw.getXYWidth() * cm;
auto layerThickness = absoThickness + gapThickness;
auto calorThickness = nofLayers * layerThickness;
auto worldSizeXY = 1.2 * calorSizeXY;
auto worldSizeZ = 1.2 * calorThickness;
auto worldSizeZ = 1.2 * caloLength;
// Get materials
auto defaultMaterial = G4Material::GetMaterial("Galactic");
auto absorberMaterial = G4Material::GetMaterial("G4_Pb");
auto gapMaterial = G4Material::GetMaterial("liquidArgon");
if ( ! defaultMaterial || ! absorberMaterial || ! gapMaterial ) {
G4ExceptionDescription msg;
msg << "Cannot retrieve materials already defined.";
G4Exception("DetectorConstruction::DefineVolumes()",
"MyCode0001", FatalException, msg);
}
auto defaultMaterial = G4Material::GetMaterial("G4_AIR");
//
// World
//
@@ -146,7 +170,7 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
//
auto calorimeterS
= new G4Box("Calorimeter", // its name
calorSizeXY/2, calorSizeXY/2, calorThickness/2); // its size
calorSizeXY/2, calorSizeXY/2, caloLength/2); // its size
auto calorLV
= new G4LogicalVolume(
@@ -155,7 +179,7 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
"Calorimeter"); // its name
new G4PVPlacement(nullptr, // no rotation
G4ThreeVector(), // at (0,0,0)
G4ThreeVector(0,0,0), // at (0,0,0)
calorLV, // its logical volume
"Calorimeter", // its name
worldLV, // its mother volume
@@ -164,82 +188,52 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
fCheckOverlaps); // checking overlaps
//
// Layer
// construct layers here; this is where the layers are added to the calorimeter
// they will be flagged active or inactive based on the isActive flag later in
// the ActionInitialization by passing the ConstructioWrapper to the EventAction class.
//
auto layerS
= new G4Box("Layer", // its name
calorSizeXY/2, calorSizeXY/2, layerThickness/2); // its size
auto layerLV
= new G4LogicalVolume(
layerS, // its solid
defaultMaterial, // its material
"Layer"); // its name
G4double position = -caloLength/2;
int layerNumber = 0;
for(auto& layer : cwLayers){
layer.name = "Layer_"+std::to_string(layerNumber);
layerNumber++;
new G4PVReplica(
"Layer", // its name
layerLV, // its logical volume
calorLV, // its mother
kZAxis, // axis of replication
nofLayers, // number of replica
layerThickness); // witdth of replica
position += layer.thickness / 2 *cm;
auto layerS
= new G4Box(layer.name, // its name
calorSizeXY/2, calorSizeXY/2, layer.thickness/2 *cm); // its size
//
// Absorber
//
auto absorberS
= new G4Box("Abso", // its name
calorSizeXY/2, calorSizeXY/2, absoThickness/2); // its size
auto layerLV
= new G4LogicalVolume(
layerS, // its solid
defaultMaterial, //G4Material::GetMaterial(layer.material), // its material
layer.name); // its name
auto absorberLV
= new G4LogicalVolume(
absorberS, // its solid
absorberMaterial, // its material
"Abso"); // its name
auto sensorLV = new G4LogicalVolume(
new G4Box("sensor", layer.sens_xwidth/2*cm, layer.sens_ywidth/2*cm, layer.thickness/2*cm),
G4Material::GetMaterial(layer.material),
"sensor");
fAbsorberPV = new G4PVPlacement(nullptr, // no rotation
G4ThreeVector(0., 0., -gapThickness / 2), // its position
absorberLV, // its logical volume
"Abso", // its name
layerLV, // its mother volume
false, // no boolean operation
0, // copy number
fCheckOverlaps); // checking overlaps
//needs RepeatPlacement for xy granularity
//use G4PVParameterised to create a grid of sensitive detectors
auto ppv = new G4PVParameterised(layer.name,
sensorLV,
layerLV, kUndefined,
layer.nx*layer.ny, new LayerParametrisation(layer,0.));
//
// Gap
//
auto gapS
= new G4Box("Gap", // its name
calorSizeXY/2, calorSizeXY/2, gapThickness/2); // its size
auto gapLV
= new G4LogicalVolume(
gapS, // its solid
gapMaterial, // its material
"Gap"); // its name
fGapPV = new G4PVPlacement(nullptr, // no rotation
G4ThreeVector(0., 0., absoThickness / 2), // its position
gapLV, // its logical volume
"Gap", // its name
layerLV, // its mother volume
false, // no boolean operation
0, // copy number
fCheckOverlaps); // checking overlaps
//
// print parameters
//
G4cout
<< G4endl
<< "------------------------------------------------------------" << G4endl
<< "---> The calorimeter is " << nofLayers << " layers of: [ "
<< absoThickness/mm << "mm of " << absorberMaterial->GetName()
<< " + "
<< gapThickness/mm << "mm of " << gapMaterial->GetName() << " ] " << G4endl
<< "------------------------------------------------------------" << G4endl;
auto pv = new G4PVPlacement(nullptr, // no rotation
G4ThreeVector(0,0,position), // at (0,0,0)
layerLV, // its logical volume
layer.name, // its name
calorLV, // its mother volume
false, // no boolean operation
0, // copy number
fCheckOverlaps); // checking overlaps
layer.assignPhysicalVolume(pv);
position += layer.thickness / 2 *cm; //assign the physical volume to the layer
}
//
// Visualization attributes
//
@@ -255,6 +249,8 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
return worldPV;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void DetectorConstruction::ConstructSDandField()