Skip to content
Snippets Groups Projects
Commit f9c8b6e1 authored by Sylvester Joosten's avatar Sylvester Joosten
Browse files

Example to access HCal segmentation from DD4hep in a ROOT script

parents
Branches
No related tags found
No related merge requests found
{
gROOT->ProcessLine(".include(\"/usr/local/include/DD4hep\")")
}
# HCal segmentation example Example
To run the example, please follow these steps:
1. Download the latest `eic-shell` container in your location of preference:
```bash
curl https://eicweb.phy.anl.gov/containers/eic_container/-/raw/master/install.sh | bash
```
2. Start eic-shell from this location:
```bash
eic-shell
```
3. Load the canyonlands v2.1 detector:
```bash
source /opt/detector/athena-canyonlands-v2.1/setup.sh
```
4. Go to your directory where you unpacked this example, and execute the root macro in
compiled mode (it will NOT work in normal interpreted mode):
```bash
root -q -b hcal_segmentation.cxx+
```
5. That's all!
dis.root 0 → 100644
File added
R__LOAD_LIBRARY(libDDRec)
R__LOAD_LIBRARY(libfmt)
// Includes necessary for compilation
#include <DD4hep/Detector.h>
#include <DDRec/CellIDPositionConverter.h>
#include <DDSegmentation/BitFieldCoder.h>
#include <iostream>
#include <vector>
#include <ROOT/RDataFrame.hxx>
#include <dd4pod/CalorimeterHitData.h>
#include <fmt/format.h>
int hcal_segmentation() {
std::cout << "hello world" << std::endl;
// Get a handle to the global detector
dd4hep::Detector* det = &(dd4hep::Detector::getInstance());
// Load the detector from the compact file. Use the correct compact file for
// your simulation here. Also make sure your environment variables are setup
// correctly Example:
// - source /opt/detector/athena-canyonlands-v2.1/setup.sh
// - use "/opt/detector/athena-canyonlands-v2.1/share/athena/athena.xml"
det->fromCompact(
"/opt/detector/athena-canyonlands-v2.1/share/athena/athena.xml");
det->volumeManager();
det->apply("DD4hepVolumeManager", 0, 0);
// Create our cellID converter. This is useful to get cell coordinates. Not
// used below.
// auto cellIDConverter = std::make_shared<const
// dd4hep::rec::CellIDPositionConverter>(*det);
// Get our readout ID spec for HcalBarrelHits
auto idSpecBarrel = det->readout("HcalBarrelHits").idSpec();
auto idSpecEndcapN = det->readout("HcalEndcapNHits").idSpec();
auto idSpecEndcapP = det->readout("HcalEndcapPHits").idSpec();
// Get our cell ID decoder
auto idDecBarrel = idSpecBarrel.decoder();
auto idDecEndcapN = idSpecEndcapN.decoder();
auto idDecEndcapP = idSpecEndcapP.decoder();
// Get a bit mask for module, layer, x and y
// module: Which stave are we on?
// layer: layer within the module (stave)
// (x,): pad coordinates within the layer
//
const int moduleIdxBarrel = idDecBarrel->index("module");
const int layerIdxBarrel = idDecBarrel->index("layer");
const int xIdxBarrel = idDecBarrel->index("x");
const int yIdxBarrel = idDecBarrel->index("y");
// same for negative and positive endcaps (they don't use module)
const int layerIdxEndcapN = idDecEndcapN->index("layer");
const int xIdxEndcapN = idDecEndcapN->index("x");
const int yIdxEndcapN = idDecEndcapN->index("y");
const int layerIdxEndcapP = idDecEndcapP->index("layer");
const int xIdxEndcapP = idDecEndcapP->index("x");
const int yIdxEndcapP = idDecEndcapP->index("y");
// Simple lambda function that prints out the hit segmentation for the barrel
auto HcalBarrelPrinter =
[&](int event, const std::vector<dd4pod::CalorimeterHitData>& h) {
for (size_t i = 0; i < h.size(); ++i) {
// Only do hits above .1 MeV
if (h[i].energyDeposit < 1e-4) {
continue;
}
const int module = idDecBarrel->get(h[i].cellID, moduleIdxBarrel);
const int layer = idDecBarrel->get(h[i].cellID, layerIdxBarrel);
const int x = idDecBarrel->get(h[i].cellID, xIdxBarrel);
const int y = idDecBarrel->get(h[i].cellID, yIdxBarrel);
fmt::print(
"HcalBarrel event: {:3d} hit: {:3d} module: {:2d} layer: {:2d} x: {:4d} "
"y {:4d} energy: {:1.3e}\n",
event, i, module, layer, x, y, h[i].energyDeposit);
}
};
// Open the ROOT file from full simulation, change this by the name of
// your ROOT file
ROOT::RDataFrame d{"events", "dis.root"};
// For this demonstration we only use the first 50 events, not something
// you would want to do in general. For production, remove this line and
// change all occurences of d50 with d
auto d50 = d.Range(50);
// Print out hit segmentation data for the barren HCal
int event = 0;
d50.Define("event", [&]{return event++;}, {})
.Foreach(HcalBarrelPrinter, {"event", "HcalBarrelHits"});
// Let's do the same for the negative endcap
auto HcalEndcapNPrinter =
[&](int event, const std::vector<dd4pod::CalorimeterHitData>& h) {
for (size_t i = 0; i < h.size(); ++i) {
// Only do hits above .1 MeV
if (h[i].energyDeposit < 1e-4) {
continue;
}
const int layer = idDecEndcapN->get(h[i].cellID, layerIdxEndcapN);
const int x = idDecEndcapN->get(h[i].cellID, xIdxEndcapN);
const int y = idDecEndcapN->get(h[i].cellID, yIdxEndcapN);
fmt::print(
"HcalEndcapN event: {:3d} hit: {:3d} layer: {:2d} x: {:4d} "
"y {:4d} energy: {:1.3e}\n",
event, i, layer, x, y, h[i].energyDeposit);
}
};
event = 0;
d50.Define("event", [&]{return event++;}, {})
.Foreach(HcalEndcapNPrinter, {"event", "HcalEndcapNHits"});
return 0;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment