Tutorial Part 2: Modifying and Adding Detectors

Setup

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

Clone part2 repo

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

Compile the detector library

Copy
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

Copy
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.

Copy
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:

Copy
<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.

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:

Copy
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.

Copy
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:

Copy
<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.

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

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

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

Copy
<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

Copy
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)

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?