diff --git a/CMakeLists.txt b/CMakeLists.txt index e367b4c..e6fe12e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/exampleB4a.cc b/exampleB4a.cc index bc988bd..ec929cd 100644 --- a/exampleB4a.cc +++ b/exampleB4a.cc @@ -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; diff --git a/include/ConstructionWrapper.hh b/include/ConstructionWrapper.hh new file mode 100644 index 0000000..b93aae5 --- /dev/null +++ b/include/ConstructionWrapper.hh @@ -0,0 +1,48 @@ + +#include +#include +#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& getLayers() ; + const std::vector& getLayers() const; + + double getXYWidth() const{ + return xywidth; + } +private: + double xywidth; + std::vector layers; +}; \ No newline at end of file diff --git a/include/DetectorConstruction.hh b/include/DetectorConstruction.hh index 681b295..9db1e56 100644 --- a/include/DetectorConstruction.hh +++ b/include/DetectorConstruction.hh @@ -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 diff --git a/src/ConstructionWrapper.cc b/src/ConstructionWrapper.cc new file mode 100644 index 0000000..4273c75 --- /dev/null +++ b/src/ConstructionWrapper.cc @@ -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 & ConstructionWrapper::getLayers() const{ + return layers; +} + +std::vector & ConstructionWrapper::getLayers(){ + return layers; +} diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc index 3d19716..44c437d 100644 --- a/src/DetectorConstruction.cc +++ b/src/DetectorConstruction.cc @@ -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 - - 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 - - // - // 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 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"); + //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.)); + + 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()