3.3. Changing geometry
The description of the componentes of a particle detector can be extremely complex. DD4hep exploits the properties of xml file format to store hierarchically the geometry and materials of each component of a given detector. Thus the detector description is given in one or more xml files.
In this section, the following topics are covered:
how tow to simulate the existing detector model of the
CLDdetectormodify its XML files
check for overlaps with Geant4
get a better understanding of what DD4hep does under the hood
Write your own C++ detector constructor and run with it
This section is adapted from the talk given by Andre Sailer in the FCC Software Tutorial 2022.
3.3.1. Setup
The first step is to set up a CentOS 7 machine with CVMFS. It can be a local machine, a remote machine (e.g. lxplus.cern.ch), a virtual machine (e.g. CernVM) or a containerized system (e.g. Docker, or Apptainer).
Then, the stack with the FCC software can be sourced as
source /cvmfs/fcc.cern.ch/sw/latest/setup.sh
The following step is to create a working directory,
mkdir mydd4heptutorial
cd mydd4heptutorial
Then, copy the CLD detector description into the working directory,
cp -r $LCGEO/FCCee/compact/FCCee_o1_v05 .
In order to plot the histograms in the terminal, the histoprint python package is used and it can be installed with the following command
pip install --user histoprint
If everything goes as expected, the contents of your working directory shall look this:
ls -ltra
total 8
drwxr-xr-x. 28 sailer zf 4096 Oct 14 11:29 ..
drwxr-xr-x. 3 sailer zf 26 Oct 14 11:30 .
drwxr-xr-x. 2 sailer zf 4096 Oct 14 11:30 FCCee_o1_v05
Another ingredient is a script to plot some information after the simulations to show that the geometry really did change. The showPlots.py script plots the per hit energy deposition and position along Z axis in the ECal endcap. The histograms are printed in the terminal to avoid the use of a graphical interface. If a graphical system is available, one can also visualize the geometry with geoDisplay or the Geant4 visualizations.
Click to show the code of the python script showPlots.py
#!/bin/env python
import sys, ROOT
from histoprint import print_hist
ROOT.gROOT.SetBatch(True)
if len(sys.argv) < 2:
print("Please specify input file"); sys.exit(1)
inputFile = sys.argv[1]; print("Reading:", inputFile)
tfile = ROOT.TFile.Open(inputFile)
myTree = tfile.Get("events")
myTree.Draw("ECalEndcapCollection.position.z>>zHist(100, 2300, 2510)",
"ECalEndcapCollection.position.z > 0")
myTree.Draw("ECalEndcapCollection.energy>>eHist(30, 0, 0.002)")
print_hist(tfile.Get("zHist"), title="Hits per Layer")
print_hist(tfile.Get("eHist"), title="Energy per Hit")
In this section, the detector geometry will be changed in every step. The geometry can be visualized using different tools, as it is explained in the next section.
3.3.2. Modify an existing XML file
The detector geometry is complex, and takes a while to be built. To speed up the following steps, the Silicon Tracker can be removed. To do so,
Open
FCCee_o1_v05/FCCee_o1_v05.xmlFind the following lines in the file, and delete/comment them
<include ref="InnerTracker_o2_v06_02.xml"/>
<include ref="OuterTracker_o2_v06_02.xml"/>
Remember that a comment in xml looks like <!-- blablabla -->.
Comment out the plugins section
<plugins> ... </plugins>in the same file.
To simulate the interaction of 100 muons with the CLD detector, the following command is used inside the directory mydd4heptutorial. It takes some minutes to build the geometry and run the actual simulation.
ddsim --compactFile FCCee_o1_v05/FCCee_o1_v05.xml \
--enableGun \
--gun.distribution uniform \
--gun.energy "10*GeV" \
--gun.particle mu- \
--numberOfEvents 100 \
--outputFile Step1_edm4hep.root
And now plot the distribution with showPlots.py
python showPlots.py Step1_edm4hep.root
The simulation runs 20 times faster than the full geometry simulation and the distributions do not change.
3.3.3. Overlap checking
Whenever you change the geometry in a non-trivial way there are the possibilities of overlaps and the following things should be kept in mind
There are no trivial changes
See point 1
Run the overlap check
Geant4 can perform the overlap check of a given geometry. This action is available for ddsim. A macro file with the Geant4 instructions has to be provided to ddsim. To do so, create a file named as overlap.mac and write these commands inside
/geometry/test/run
exit
And then we run ddsim with this macro file, and dump the output to a text file for easy browsing:
ddsim --compactFile FCCee_o1_v05/FCCee_o1_v05.xml \
--runType run \
--macroFile overlap.mac > overlapDump &
With the full detector model including the tracker this would take about 30 minutes. The output shows some exceptions:
The exceptions about the VertexEndcapModule are due to a too small envelope, and could be easily fixed
The exceptions about the Boolean Volume are more vexing, let’s ignore those for now. The HOMAbsorber include can be dropped as well to silent these errors.
3.3.4. Modifying the ECAL
In this section, it is shown how to change the number of layers and silicon thicknesses of the ECal Endcap. To do so,
Open the file
FCCee_o1_v05/ECalEndcap_o2_v01_03.xmlFind this block:
<layer repeat="40" vis="ECalLayerVis">
<slice material = "TungstenDens24" thickness = "1.90*mm" vis="ECalAbsorberVis" radiator="yes"/>
<slice material = "G10" thickness = "0.15*mm" vis="InvisibleNoDaughters"/>
<slice material = "GroundOrHVMix" thickness = "0.10*mm" vis="ECalAbsorberVis"/>
<slice material = "Silicon" thickness = "0.50*mm" sensitive="yes" limits="cal_limits"
vis="ECalSensitiveVis"/>
<slice material = "Air" thickness = "0.10*mm" vis="InvisibleNoDaughters"/>
<slice material = "siPCBMix" thickness = "1.30*mm" vis="ECalAbsorberVis"/>
<slice material = "Air" thickness = "0.25*mm" vis="InvisibleNoDaughters"/>
<slice material = "G10" thickness = "0.75*mm" vis="InvisibleNoDaughters"/>
</layer>
Change the number of layers from 40 to 20, and the silicon thickness from
0.50*mmto1.00*mm.
Click to show that piece of code from FCCee_o1_v05/ECalEndcap_o2_v01_03.xml after the change
<layer repeat="20" vis="ECalLayerVis">
<slice material = "TungstenDens24" thickness = "1.90*mm" vis="ECalAbsorberVis" radiator="yes"/>
<slice material = "G10" thickness = "0.15*mm" vis="InvisibleNoDaughters"/>
<slice material = "GroundOrHVMix" thickness = "0.10*mm" vis="ECalAbsorberVis"/>
<slice material = "Silicon" thickness = "1.00*mm" sensitive="yes" limits="cal_limits"
vis="ECalSensitiveVis"/>
<slice material = "Air" thickness = "0.10*mm" vis="InvisibleNoDaughters"/>
<slice material = "siPCBMix" thickness = "1.30*mm" vis="ECalAbsorberVis"/>
<slice material = "Air" thickness = "0.25*mm" vis="InvisibleNoDaughters"/>
<slice material = "G10" thickness = "0.75*mm" vis="InvisibleNoDaughters"/>
</layer>
Simulate again with the new ECal,
ddsim --compactFile FCCee_o1_v05/FCCee_o1_v05.xml \
--enableGun \
--gun.distribution uniform \
--gun.energy "10*GeV" \
--gun.particle mu- \
--numberOfEvents 100 \
--outputFile Step2_edm4hep.root
Compare the histograms produced from Step1_edm4hep.root and Step2_edm4hep.root.
3.3.5. Geometry driver modifications
The type attribute of the detector tag tells DD4hep which detector constructor to load
<detector name="ECalEndcap"
type="GenericCalEndcap_o1_v01"
id="DetID_ECal_Endcap"
readout="ECalEndcapCollection"
vis="ECALVis" >
</detector>
DD4hep’s plugin service will look in the *.components files it finds via the LD_LIBRARY_PATH (or DD4HEP_LIBRARY_PATH on macOS because of SIP) environment variables, and load the library on-demand, and then instantiate the function.
For example, to know which library contains the GenericCalEndcap plugin
for DIR in $(echo $LD_LIBRARY_PATH | tr ":" " ") ; do
grep GenericCalEndcap $DIR/*.components 2> /dev/null;
done
Usually we would add our detector to an existing project, but because overwriting existing libraries in the environment requires some work, we start with a new project. This requires two main steps, first creating a new detector constructor, and second creating the detector description in a xml file. Start by executing the following code inside the directory mydd4heptutorial, in order to create the detector constructor project
mkdir MyFirstDetector
cd MyFirstDetector
touch CMakeLists.txt
mkdir src
touch src/MyFirstDetector.cpp
The CMakeLists.txt file has to contain the following code
cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
set(PackageName MyFirstDetector)
project(${PackageName})
find_package(DD4hep REQUIRED COMPONENTS DDG4)
# our sources
set(sources ./src/MyFirstDetector.cpp)
# create our library and make the components file
add_dd4hep_plugin(${PackageName} SHARED ${sources})
# link it with DDCore, or whatever you need
target_link_libraries(${PackageName} DD4hep::DDCore)
# Create this_package.sh file, and install
dd4hep_instantiate_package(${PackageName})
And the MyFirstDetector.cpp file has to contain the following code
#include "DD4hep/DetFactoryHelper.h"
#include "XML/Layering.h"
#include "XML/Utilities.h"
#include "DDRec/DetectorData.h"
static dd4hep::Ref_t create_detector(
dd4hep::Detector &theDetector,
xml_h entities,
dd4hep::SensitiveDetector sens)
{
// XML Detector Element (confusingly also XML::DetElement)
xml_det_t x_det = entities;
// DetElement of our detector instance, attach additional information, sub-elements...
// uses name of detector and ID number as defined in the XML detector tag
std::string detName = x_det.nameStr();
sens.setType("calorimeter");
dd4hep::DetElement sdet(detName, x_det.id());
// get the dimensions tag
xml_dim_t dim = x_det.dimensions();
// read its attributes
double rmin = dim.rmin();
double rmax = dim.rmax();
double zmax = dim.zmax();
// Make a Cylinder
dd4hep::Tube envelope(rmin, rmax, zmax);
dd4hep::Material air = theDetector.air();
dd4hep::Volume envelopeVol(detName + "_envelope",
envelope,
air);
dd4hep::PlacedVolume physvol =
theDetector.pickMotherVolume(sdet).placeVolume(envelopeVol);
// add system ID and identify as barrel (as opposed to endcap +/-1)
physvol.addPhysVolID("system", sdet.id()).addPhysVolID(_U(side), 0);
sdet.setPlacement(physvol);
return sdet;
}
DECLARE_DETELEMENT(MyFirstDetector, create_detector)
Still inside the MyFirstDetector directory, execute the following shell commands,
mkdir build install
cd build
cmake -D CMAKE_INSTALL_PREFIX=$PWD/../install ..
make install
source ../install/bin/thisMyFirstDetector.sh
Look at $PWD/../install and note the content of the bin and lib directories, have a look at the libMyFirstDetector.components file.
After creating the detector constructor project, a detector description has to be given as a xml file. In the FCCee_o1_v05 directory, create a mydetector.xml file with the following content
Click to show the code of mydetector.xml
<lccdd>
<define>
<constant name="ECal_cell_size" value="5.1*mm" />
</define>
<readouts>
<readout name="MyReadout">
<segmentation type="GridRPhiEta"
grid_size_eta="ECal_cell_size" phi_bins="360"
offset_r="ECalBarrel_inner_radius"
grid_size_r="1*cm" />
<id>system:5,side:2,module:8,stave:4,layer:9,submodule:4,r:32:10,eta:-11,phi:-11</id>
</readout>
</readouts>
<display>
<vis name="MyVis" alpha="0.1"
r="0.1" g=".5" b=".5"
showDaughters="true"
visible="false" />
</display>
<detectors>
<detector name="MyDetectorName"
type="MyFirstDetector"
id="1234"
readout="MyReadout"
vis="MyVis">
<dimensions
zmax="ECalBarrel_half_length"
rmin="ECalBarrel_inner_radius"
rmax="ECalBarrel_outer_radius" />
</detector>
</detectors>
</lccdd>
Now, modify the FCCee_o1_v05/FCCee_o1_v05.xml and replace the include for the ECalBarrel
with mydetector.xml. After every change in the geometry, the overlap check must be run inside FCCee_o1_v05 directory
ddsim --compactFile FCCee_o1_v05/FCCee_o1_v05.xml \
--runType run \
--macroFile overlap.mac > overlapDump2 &
If there are no exceptions, we can continue. Now, we add some layers to the detector. First, we add a layer component inside the <detector> ... </detector> section of FCCee_o1_v05/mydetector.xml detector description file.
Click to show the code of mydetector.xml
<lccdd>
<define>
<constant name="ECal_cell_size" value="5.1*mm" />
</define>
<readouts>
<readout name="MyReadout">
<segmentation type="GridRPhiEta"
grid_size_eta="ECal_cell_size" phi_bins="360"
offset_r="ECalBarrel_inner_radius"
grid_size_r="1*cm" />
<id>system:5,side:2,module:8,stave:4,layer:9,submodule:4,r:32:10,eta:-11,phi:-11</id>
</readout>
</readouts>
<display>
<vis name="MyVis" alpha="0.1"
r="0.1" g=".5" b=".5"
showDaughters="true"
visible="false" />
</display>
<detectors>
<detector name="MyDetectorName"
type="MyFirstDetector"
id="1234"
readout="MyReadout"
vis="MyVis">
<dimensions
zmax="ECalBarrel_half_length"
rmin="ECalBarrel_inner_radius"
rmax="ECalBarrel_outer_radius" />
<layer repeat="10" vis="MyVis">
<slice material="Iron" thickness="1*cm" />
<slice material="G10" thickness="1*mm" />
<slice material="Silicon" thickness="1*mm" sensitive="true" />
<slice material="G10" thickness="1*mm" />
</layer>
</detector>
</detectors>
</lccdd>
In order to actually build the layers that were included in the detector description, we have to add the corresponding interpretation in the detector constructor defined in MyFirstDetector.cpp, which should look like the following:
Click to show the code of MyFirstDetector.cpp
#include "DD4hep/DetFactoryHelper.h"
#include "DDRec/DetectorData.h"
#include "XML/Layering.h"
#include "XML/Utilities.h"
static dd4hep::Ref_t create_detector(dd4hep::Detector &theDetector,
xml_h entities,
dd4hep::SensitiveDetector sens) {
// XML Detector Element (confusingly also XML::DetElement)
xml_det_t x_det = entities;
// DetElement of our detector instance, attach additional information,
// sub-elements... uses name of detector and ID number as defined in the XML
// detector tag
std::string detName = x_det.nameStr();
sens.setType("calorimeter");
dd4hep::DetElement sdet(detName, x_det.id());
// get the dimensions tag
xml_dim_t dim = x_det.dimensions();
// read its attributes
double rmin = dim.rmin();
double rmax = dim.rmax();
double zmax = dim.zmax();
// Make a Cylinder
dd4hep::Tube envelope(rmin, rmax, zmax);
dd4hep::Material air = theDetector.air();
dd4hep::Volume envelopeVol(detName + "_envelope", envelope, air);
dd4hep::PlacedVolume physvol =
theDetector.pickMotherVolume(sdet).placeVolume(envelopeVol);
// add system ID and identify as barrel (as opposed to endcap +/-1)
physvol.addPhysVolID("system", sdet.id()).addPhysVolID(_U(side), 0);
sdet.setPlacement(physvol);
// Interpretation of layers:
double currentInnerRadius = rmin; // running inner radius
dd4hep::Layering layering(x_det); // convenience class
int layerNum = 0;
for (xml_coll_t c(x_det, _U(layer)); c; ++c, ++layerNum) {
xml_comp_t x_layer = c;
const dd4hep::Layer *lay =
layering.layer(layerNum); // Get the layer from the layering engine.
const double layerThickness = lay->thickness();
// loop over the number of repetitions
for (int i = 0, repeat = x_layer.repeat(); i < repeat; ++i, ++layerNum) {
std::string layerName = detName + dd4hep::_toString(layerNum, "_layer%d");
// make a volume for the layer
dd4hep::Tube layerTube(currentInnerRadius,
currentInnerRadius + layerThickness, zmax);
dd4hep::Volume layerVol(layerName, layerTube, air);
dd4hep::DetElement layerElement(sdet, layerName, layerNum);
dd4hep::PlacedVolume layerVolPlaced = envelopeVol.placeVolume(layerVol);
layerVolPlaced.addPhysVolID("layer", layerNum);
int sliceNum = 0;
for (xml_coll_t slice(x_layer, _U(slice)); slice; ++slice, ++sliceNum) {
xml_comp_t x_slice = slice;
double sliceThickness = x_slice.thickness();
dd4hep::Material sliceMat = theDetector.material(x_slice.materialStr());
std::string sliceName =
layerName + dd4hep::_toString(sliceNum, "slice%d");
dd4hep::Tube sliceTube(currentInnerRadius,
currentInnerRadius + sliceThickness, zmax);
dd4hep::Volume sliceVol(sliceName, sliceTube, sliceMat);
if (x_slice.isSensitive()) {
sliceVol.setSensitiveDetector(sens);
}
// place the slice in the layer
layerVol.placeVolume(sliceVol);
currentInnerRadius += sliceThickness;
} // slices
} // repetitions
} // layers
return sdet;
}
DECLARE_DETELEMENT(MyFirstDetector, create_detector)
To take effect, the code has to be recompiled,
cd MyFirstDetector/build
make install
After every change in the geometry, the overlap check must be run. To do so, go to the top working directory and run it as before
ddsim --compactFile FCCee_o1_v05/FCCee_o1_v05.xml \
--runType run \
--macroFile overlap.mac > overlapDump3 &
The number of checks now increased, and if everything is done properly, no overlaps are found. Now the geometry is ready to go for a simulation. As it was done previously, we use ddsim to run a simulation of 100 muons in the geometry that has just been defined, as in the previous steps
ddsim --compactFile FCCee_o1_v05/FCCee_o1_v05.xml \
--enableGun \
--gun.distribution uniform \
--gun.energy "10*GeV" \
--gun.particle mu- \
--numberOfEvents 100 \
--outputFile Step3_edm4hep.root
Modify showPlots.py to display properties from this collection (MyReadout).
Click to show the code of the python script showPlots.py
#!/bin/env python
import sys, ROOT
from histoprint import print_hist
ROOT.gROOT.SetBatch(True)
if len(sys.argv) < 2:
print("Please specify input file"); sys.exit(1)
inputFile = sys.argv[1]; print("Reading:", inputFile)
tfile = ROOT.TFile.Open(inputFile)
myTree = tfile.Get("events")
myTree.Draw("sqrt(MyReadout.position.x*MyReadout.position.x+MyReadout.position.y*MyReadout.position.y)>>rHist(100, 2150, 2352)")
myTree.Draw("MyReadout.energy>>eHist(30, 0, 0.002)")
print_hist(tfile.Get("rHist"), title="Hits per Layer")
print_hist(tfile.Get("eHist"), title="Energy per Hit")
3.3.6. Final note
Changing geometries with DD4hep can range from trivial to sophisticated. The overlap checker shall be run after every change.
Check the following links for further information and support:
Browse the DD4hep documentation and beginners guide
Post your questions or suggestions at FCC software forum or Mattermost channel