EIC software

Tutorial Part 2: Modifying and Adding Detectors

Setup

Note all these commands assume you are in an eic-shell singularity session.

Clone part2 repo

git clone https://eicweb.phy.anl.gov/EIC/tutorials/ip6_tutorial_2.git part2
cd part2

Compile the detector library

cmake -B build -S . -DCMAKE_INSTALL_PREFIX=$ATHENA_PREFIX -DCMAKE_CXX_STANDARD=17
cmake --build build -j4 -- install

How to build a detector from scratch

Building a new (generic) detector using dd4hep is rather straight forward if we make the following assumptions.

  1. We will use the built-in sensitive detectors
  2. We will use the dd4pod data model (i.e. use npsim)

Compiling a new detector

For this tutorial we will build a more detailed GEM detector. We will discuss the detector built in the source file src/TrapEndcapTracker_geo.cpp.

In the CMakeLists.txt it adds all cpp files, otherwise we would do

dd4hep_add_plugin(${a_lib_name} SOURCES
src/TrapEndcapTracker_geo.cpp
...
)

Building the geometry (cpp + xml)

The work of constructing the detector is done in a static function (here called build_detector) that is registered using the DD4hep plugin macro DECLARE_DETELEMENT.

static Ref_t build_detector(Detector& dtor, xml_h e, SensitiveDetector sens)
{
xml_det_t x_det = e;
Material air = dtor.air();
...
}
DECLARE_DETELEMENT(MyGEMTrackerEndcap, build_detector)

The argument signature of the build_detector is:

  • Detector& dtor: This handle provides the main hook to the detector tree (dd4hep::Detector).
  • xml_h e: Handle to the XML <detector> tag associated with the detector in the "compact" detector description file (more on this later). This provides the mechanism for feeding in the run-time construction parameters.
  • SensitiveDetector sens: The sensitive detector to be assigned to the sensitive volumes/elements of the detector.

The DD4hep plugin macro DECLARE_DETELEMENT(MyGEMTrackerEndcap, build_detector) stamps out the necessary boiler plate code to register a new detector called MyGEMTrackerEndcap which is build by calling build_detector.

Compact detector description entry element

The <detector> tag defines a new instance of a detector and requires the attributes "id", "name", and "type". For example from part 1:

<detector id="2" name="GEMTracker" vis="RedVis" type="my_GEMTracker" readout="GEMTrackerHits" >
<layer id="1" z="-100 *cm" inner_r="40*cm" outer_r="120*cm" phi0_offset="0.0*deg" />
<layer id="2" z="-80 *cm" inner_r="30*cm" outer_r="90*cm" phi0_offset="0.0*deg" />
<layer id="3" z="-60 *cm" inner_r="20*cm" outer_r="70*cm" phi0_offset="0.0*deg" />
<layer id="4" z="-40 *cm" inner_r="10*cm" outer_r="20.0*cm" phi0_offset="0.0*deg" />
<layer id="5" z=" 40 *cm" inner_r="10*cm" outer_r="20.0*cm" phi0_offset="0.0*deg" />
<layer id="6" z=" 60 *cm" inner_r="25*cm" outer_r="70.0*cm" phi0_offset="0.0*deg" />
<layer id="7" z=" 80 *cm" inner_r="30*cm" outer_r="90.0*cm" phi0_offset="0.0*deg" />
<layer id="8" z="100 *cm" inner_r="40*cm" outer_r="100.0*cm" phi0_offset="0.0*deg" />
</detector>

This defines an instance of the detector named "GEMTracker" of type "my_GEMTracker" (i.e. the type-name given in the first argument of the DD4hep DECLARE_DETELEMENT macro) and with id=2. Each <detector> must have a unique id. The additional attributes (vis, readout, zoffset) are optional.

The detector tag is provided as the second argument in the build_detector function. It can be parsed how ever you want in order to extract detector information. The allowed attributes are listed in UnicodeValues.h where it is clear how to add new attributes.

Geometry Construction

If you are familiar with Geant4 or TGeo geometry construction then this will be easy. DD4hep has TGeo under hood but there are a few naming tweaks. The following table will help orient the user.

DD4hepGeant4*TGeo
SolidG4VSolidTGeoShape
VolumeG4LogicalVolumeTGeoVolume
PlacedVolumeG4PVPlacementTGeoNode
ElementG4ElementTGeoElement
MaterialG4MaterialTGeoMaterial/TGeoMedium
XML Parsing Tip : Provide good default values

If you have a detector parameter which we later will tweak (while optimizing the design) try to get the value from the xml element but provide a good default value. For example:

double radius = ( x_det.hasAttr(_Unicode(radius)) ) ? x_det.attr<double>(_Unicode(radius)) : 5.0*dd4hep::cm;
//or
double radius = dd4hep::getAttrOrDefault(x_det, _Unicode(radius), 5.0*dd4hep::cm)

This provides a default radius of 5 cm when x_det does not have a "radius" attribute defined. We will return to this later.

Critical parts of build_detector

We will now look at the function in src/GEMTrackerDisc_geo.cpp. Note that this is not a good example to copy, rather it is a simple example to walk through.

static Ref_t build_detector(Detector& lcdd, xml_h e, SensitiveDetector sens)
{
typedef vector<PlacedVolume> Placements;
xml_det_t x_det = e;
Material air = lcdd.air();
Material carbon = lcdd.material("CarbonFiber");
Material silicon = lcdd.material("SiliconOxide");
int det_id = x_det.id();
string det_name = x_det.nameStr();
PlacedVolume pv;
DetElement sdet(det_name, det_id);
Assembly assembly(det_name+"_assembly");
sens.setType("tracker");
string module_name = "GEM";
double thickness = 0.01*dd4hep::cm;
int N_layers = 0;
for(xml_coll_t lay( x_det, _U(layer) ); lay; ++lay, ++N_layers) {
xml_comp_t x_layer = lay;
double inner_r = x_layer.attr<double>( _Unicode(inner_r) ) ;
double outer_r = x_layer.attr<double>( _Unicode(outer_r) ) ;
double phi0_offset = x_layer.attr<double>( _Unicode(phi0_offset) ) ;
double z = x_layer.attr<double>( _Unicode(z) ) ;
int layer_id = x_layer.id();//attr<double>( _Unicode(z) ) ;
string layer_name = std::string("gem_layer") + std::to_string(layer_id);
Tube gem_layer(inner_r, outer_r, thickness / 2.0);
Volume gem_layer_vol("gem_layer_vol", gem_layer, carbon);
gem_layer_vol.setSensitiveDetector(sens);
DetElement layer_DE( sdet, _toString(layer_id,"layer%d"), layer_id );
pv = assembly.placeVolume( gem_layer_vol, Transform3D(RotationZ(phi0_offset),Position(0.0,0.0,z)) );
pv.addPhysVolID( "layer", layer_id );
layer_DE.setPlacement(pv);
}
sdet.setAttributes(lcdd, assembly,x_det.regionStr(),x_det.limitsStr(),x_det.visStr());
pv = lcdd.pickMotherVolume(sdet).placeVolume(assembly);
pv.addPhysVolID("system", det_id); // Set the subdetector system ID.
sdet.setPlacement(pv);
assembly->GetShape()->ComputeBBox() ;
return sdet;
}
DECLARE_DETELEMENT(my_GEMTracker, build_detector)

Here is the XML again:

<detector id="2" name="GEMTracker" vis="RedVis" type="my_GEMTracker" readout="GEMTrackerHits" >
<layer id="1" z="-100 *cm" inner_r="40*cm" outer_r="120*cm" phi0_offset="0.0*deg" />
<layer id="2" z="-80 *cm" inner_r="30*cm" outer_r="90*cm" phi0_offset="0.0*deg" />
<layer id="3" z="-60 *cm" inner_r="20*cm" outer_r="70*cm" phi0_offset="0.0*deg" />
<layer id="4" z="-40 *cm" inner_r="10*cm" outer_r="20.0*cm" phi0_offset="0.0*deg" />
<layer id="5" z=" 40 *cm" inner_r="10*cm" outer_r="20.0*cm" phi0_offset="0.0*deg" />
<layer id="6" z=" 60 *cm" inner_r="25*cm" outer_r="70.0*cm" phi0_offset="0.0*deg" />
<layer id="7" z=" 80 *cm" inner_r="30*cm" outer_r="90.0*cm" phi0_offset="0.0*deg" />
<layer id="8" z="100 *cm" inner_r="40*cm" outer_r="100.0*cm" phi0_offset="0.0*deg" />
</detector>

What problems do you see with this construction?

  • rigid
  • hard coded
  • no defaults
  • single material
  • lots more...

What materials exist?

To browse the list of materials and elements (other than looking in the compact files), go to geo viewer and scroll through the materials. browse_materials

Here we are grabbing the materials that are assumed to be already defined. Also we are getting the detector id and name defined in the detector tag.

Next we define an Assembly volume. Here we also stumble upon the important class dd4hep::DetElement. It is a means of providing the detector hierarchy/tree, but doesn't necessarily have to map exactly to detector geometry. However, it typically will typically parallel the geometry (and probably should).

A more detailed GEM Tracker

Uncomment the include line in gem_tracker.xml

<include ref="compact/gem_tracker_endcap.xml"/>

And have a look at the detector part of compact/gem_tracker_endcap.xml.

<detector
id="GEMTrackerEndcap_ID"
name="GEMTrackerEndcap"
type="MyGEMTrackerEndcap"
readout="GEMTrackerEndcapHits"
vis="BlueVis"
reflect="false">
<module name="GEMModule1" vis="GreenVis">
<trd x1="GEMTrackerEndcapFoilX1/2.0" x2="GEMTrackerEndcapFoilX2/2.0" z="GEMTrackerEndcapFoilY/2"/>
<comment> Going from HV side to readout side</comment>
<module_component thickness="0.127 * mm" material="Mylar"/>
<module_component thickness="50.0*um" material="Kapton" name="entrance_window"/>
<module_component thickness=" 3.0*mm" material="Ar10CO2" name="entrance region" />
<module_component thickness="50.0*um" material="Kapton"/>
<module_component thickness=" 3.0*um" material="Copper"/>
<module_component thickness=" 3.0*mm" material="Ar10CO2" name="drift region"/>
<module_component thickness="30.0*um" material="Kapton" name="gem_foil"/>
<module_component thickness=" 3.0*um" material="Copper" name="gem_foil_Cu"/>
<module_component thickness=" 2.0*mm" material="Ar10CO2" name="transfer region I"/>
<module_component thickness="30.0*um" material="Kapton" name="gem_foil"/>
<module_component thickness=" 3.0*um" material="Copper" name="gem_foil_Cu"/>
<module_component thickness=" 2.0*mm" material="Ar10CO2" name="transfer region II"/>
<module_component thickness="30.0*um" material="Kapton" name="gem_foil"/>
<module_component thickness=" 3.0*um" material="Copper" name="gem_foil_Cu"/>
<module_component thickness=" 2.0*mm" material="Ar10CO2" name="induction region"/>
<module_component thickness="30.0*um" material="Kapton" name="readout" sensitive="true"/>
<module_component thickness=" 3.0*um" material="Copper" name="readout_Cu"/>
<module_component thickness="127.0*um" material="Mylar"/>
</module>
<module name="GEMSupportModule1" vis="OrangeVis">
<trd x1="GEMTrackerEndcapFoilX2/2.0" x2="GEMTrackerEndcapFoilX1/2.0" z="GEMTrackerEndcapFrameBotEdge_width"/>
<module_component thickness="GEMTrackerEndcapFrame_thickness" material="Mylar"/>
</module>
<comment>
GEMSupportModule2 is the support in between gem foils.
</comment>
<module name="GEMSupportModule2" vis="OrangeVis">
<trd x1="GEMTrackerEndcapFrameSideEdge_width" x2="GEMTrackerEndcapFrameSideEdge_width" z="GEMTrackerEndcapFoilY/2"/>
<module_component thickness="4.0*mm" material="Mylar"/>
</module>
<layer id="1" >
<ring vis="PurpleVis"
r="GEMTrackerEndcapFoil_rmin+GEMTrackerEndcapFoilY/2.0"
zstart="GEMTrackerEndcap_zmin + 0.5*GEMTrackerEndcapLayer_thickness"
nmodules="12" dz="10 * mm" module="GEMModule1" />
<ring vis="PurpleVis" phi0="15.0*degree"
r="GEMTrackerEndcapFoil_rmin+GEMTrackerEndcapFoilY/2.0"
zstart="GEMTrackerEndcap_zmin + 0.5*GEMTrackerEndcapLayer_thickness"
nmodules="12" dz="0 * mm" module="GEMSupportModule2" />
</layer>
</detector>

Here you see a different strategy: build a module, then build layers constructed from the module.

The corresponding detector constructor function is in src/TrapEndcapTracker_geo.cpp

static Ref_t build_detector(Detector& description, xml_h e, SensitiveDetector sens) {
typedef vector<PlacedVolume> Placements;
xml_det_t x_det = e;
Material vacuum = description.vacuum();
int det_id = x_det.id();
string det_name = x_det.nameStr();
bool reflect = x_det.reflect(false);
DetElement sdet(det_name, det_id);
Assembly assembly(det_name);
Volume motherVol = description.pickMotherVolume(sdet);
int m_id = 0, c_id = 0, n_sensor = 0;
map<string, Volume> modules;
map<string, Placements> sensitives;
PlacedVolume pv;
assembly.setVisAttributes(description.invisible());
sens.setType("tracker");
// Loop over and build modules
for (xml_coll_t mi(x_det, _U(module)); mi; ++mi, ++m_id) {
xml_comp_t x_mod = mi;
string m_nam = x_mod.nameStr();
xml_comp_t trd = x_mod.trd();
double posY;
double x1 = trd.x1();
double x2 = trd.x2();
double z = trd.z();
double y1, y2, total_thickness = 0.;
xml_coll_t ci(x_mod, _U(module_component));
for (ci.reset(), total_thickness = 0.0; ci; ++ci)
total_thickness += xml_comp_t(ci).thickness();
y1 = y2 = total_thickness / 2;
Volume m_volume(m_nam, Trapezoid(x1, x2, y1, y2, z), vacuum);
m_volume.setVisAttributes(description.visAttributes(x_mod.visStr()));
for (ci.reset(), n_sensor = 1, c_id = 0, posY = -y1; ci; ++ci, ++c_id) {
xml_comp_t c = ci;
double c_thick = c.thickness();
auto comp_x1 = getAttrOrDefault(c, _Unicode(x1), x1);
auto comp_x2 = getAttrOrDefault(c, _Unicode(x2), x2);
auto comp_height = getAttrOrDefault(c, _Unicode(height), z);
Material c_mat = description.material(c.materialStr());
string c_name = _toString(c_id, "component%d");
Volume c_vol(c_name, Trapezoid(comp_x1, comp_x2, c_thick / 2e0, c_thick / 2e0, comp_height), c_mat);
// use the module vis attributes if not set for component.
auto comp_vis = x_mod.visStr();
if(( c.visStr().size() >0 ) ) {
comp_vis = c.visStr();
}
c_vol.setVisAttributes(description.visAttributes(comp_vis));
pv = m_volume.placeVolume(c_vol, Position(0, posY + c_thick / 2, 0));
if (c.isSensitive()) {
sdet.check(n_sensor > 2, "SiTrackerEndcap2::fromCompact: " + c_name + " Max of 2 modules allowed!");
pv.addPhysVolID("sensor", n_sensor);
c_vol.setSensitiveDetector(sens);
sensitives[m_nam].push_back(pv);
++n_sensor;
}
posY += c_thick;
}
// Modules are stored here
modules[m_nam] = m_volume;
}
// Construct each layer
for (xml_coll_t li(x_det, _U(layer)); li; ++li) {
xml_comp_t x_layer(li);
int l_id = x_layer.id();
int mod_num = 1;
for (xml_coll_t ri(x_layer, _U(ring)); ri; ++ri) {
xml_comp_t x_ring = ri;
double r = x_ring.r();
double phi0 = x_ring.phi0(0);
double zstart = x_ring.zstart();
double dz = x_ring.dz(0);
int nmodules = x_ring.nmodules();
string m_nam = x_ring.moduleStr();
Volume m_vol = modules[m_nam];
double iphi = 2 * M_PI / nmodules;
double phi = phi0;
Placements& sensVols = sensitives[m_nam];
for (int k = 0; k < nmodules; ++k) {
string m_base = _toString(l_id, "layer%d") + _toString(mod_num, "_module%d");
double x = -r * std::cos(phi);
double y = -r * std::sin(phi);
DetElement module(sdet, m_base + "_pos", det_id);
pv = assembly.placeVolume(m_vol,
Transform3D(RotationZYX(0, -M_PI / 2 - phi, -M_PI / 2), Position(x, y, zstart + dz)));
pv.addPhysVolID("barrel", 1).addPhysVolID("layer", l_id).addPhysVolID("module", mod_num);
module.setPlacement(pv);
for (size_t ic = 0; ic < sensVols.size(); ++ic) {
PlacedVolume sens_pv = sensVols[ic];
DetElement comp_elt(module, sens_pv.volume().name(), mod_num);
comp_elt.setPlacement(sens_pv);
}
if (reflect) {
pv = assembly.placeVolume(
m_vol, Transform3D(RotationZYX(M_PI, -M_PI / 2 - phi, -M_PI / 2), Position(x, y, -zstart - dz)));
pv.addPhysVolID("barrel", 2).addPhysVolID("layer", l_id).addPhysVolID("module", mod_num);
DetElement r_module(sdet, m_base + "_neg", det_id);
r_module.setPlacement(pv);
for (size_t ic = 0; ic < sensVols.size(); ++ic) {
PlacedVolume sens_pv = sensVols[ic];
DetElement comp_elt(r_module, sens_pv.volume().name(), mod_num);
comp_elt.setPlacement(sens_pv);
}
}
dz = -dz;
phi += iphi;
++mod_num;
}
}
}
pv = motherVol.placeVolume(assembly);
pv.addPhysVolID("system", det_id);
sdet.setPlacement(pv);
return sdet;
}
// clang-format off
DECLARE_DETELEMENT(TrapEndcapTracker, build_detector)
DECLARE_DETELEMENT(MyGEMTrackerEndcap, build_detector)

gem tracker

Hands on challenges

Update gem tracker

Swap out the gem tracker from part 1 with the detailed construction.

Add top and bottom support for the GEM foil frame.

There is a side currently there

Make the modules overlap

Currently the modules do not overlap much. Can you come up with a way to modify to get a nice overlap?

What bug exists with the negative endcap?

If you uncomment the negative endcap, can you find out what is wrong with the construction?

Edit this page on eicweb