Possible to dump all available readout elements?

[alexandertuna]

Hi,

@headunderheels and I are simulating an electromagnetic calorimeter with dd4hep and a compact detector description. Thanks for the very nice code! We have a small question: is it possible to dump all readout elements somehow?

For example, if we simulate a particle passing through the ecal, we can get a collection of readout elements where the particle deposited energy. But is it possible to get a list of all readout elements, without simulating a particle? It would be useful for us to have this info.

Thanks for your help, we're new to this!

---

Our ecal, for example: https://github.com/madbaron/detector-simulation/blob/KITP_10TeV/geometries/MuColl_10TeV_v0A/ECalEndcap_o2_v01_02.xml

<lccdd>

	<readouts>
		<readout name="ECalEndcapCollection">
			<segmentation type="CartesianGridXY" grid_size_x="ECal_cell_size" grid_size_y="ECal_cell_size" />
			<id>${GlobalCalorimeterReadoutID}</id>
		</readout>
	</readouts>

	<detectors>
		<detector name="ECalEndcap" type="GenericCalEndcap_o1_v01" id="DetID_ECal_Endcap" readout="ECalEndcapCollection" vis="ECALVis" >

        and so on

[atolosadelgado]

Hi

The following snippet shows how to loop over the detector element map and sensitive detector map:

#include "DD4hep/Detector.h"
#include "DD4hep/DD4hepUnits.h"
#include "DD4hep/BitFieldCoder.h"

R__LOAD_LIBRARY(libDDCore) 

auto lcdd = &(dd4hep::Detector::getInstance());
lcdd->fromCompact( path_to_main_XML )

// map of detectors, geometry
for( auto & [name, det] : lcdd->detectors() ) std::cout << name.c_str() << std::endl;

// map of sensitive detectors, readout
for( auto & [name, sd] : lcdd->sensitiveDetectors() ) std::cout << name.c_str() << std::endl;

Is this snippet what you were looking for?

Best,
Alvaro

[alexandertuna]

Thanks for the fast response @atolosadelgado ! We’ll try this asap and let you know.

[alexandertuna]

Hi again @atolosadelgado , I think this is almost what we want. We're interested in dumping a full list of readout elements, not just readouts. In our example, we have:

<detector> -> <readout> -> <segmentation> -> CartesianGridXY

We're interested in dumping each element of the CartesianGridXY. I'll look for methods within dd4hep for finding this, but feel free to post any suggestions if you have them.

[vvolkl]

Yes, this is not that straightforward, as it depends on your segmentation. Previous discussion here: #580 , along with a partial solution

[atolosadelgado]

Hi again @atolosadelgado , I think this is almost what we want. We're interested in dumping a full list of readout elements, not just readouts. In our example, we have:

<detector> -> <readout> -> <segmentation> -> CartesianGridXY

We're interested in dumping each element of the CartesianGridXY. I'll look for methods within dd4hep for finding this, but feel free to post any suggestions if you have them.

Hi,

I think the last line of the snippet now provides what you are looking for :)

#include "DD4hep/Detector.h"
#include "DD4hep/DD4hepUnits.h"
#include "DD4hep/BitFieldCoder.h"

R__LOAD_LIBRARY(libDDCore) 

auto lcdd = &(dd4hep::Detector::getInstance());
lcdd->fromCompact( path_to_main_XML )

// map of detectors, geometry
for( auto & [name, det] : lcdd->detectors() ) std::cout << name.c_str() << std::endl;

// map of sensitive detectors, readout
for( auto & [name, sd] : lcdd->sensitiveDetectors() ) std::cout << name.c_str() << std::endl;

// map of readouts
for( auto & [name, ro_b] : lcdd->readouts() )
{ 
  dd4hep::Readout r = ro_b; 
  std::cout << r.segmentation().type()<< std::endl;
}

[BrieucF]

Hi @alexandertuna ,

It looks indeed that, since the segmentation does not know, a priori, what is the extent of the volume to be segmented there is no magic function to get what you want. I think you have to retrieve the segmentation and the dimension of the volume at the same time to be able to extract extremas.

See for instance: https://github.com/key4hep/k4geo/blob/main/detectorCommon/src/DetUtils_k4geo.cpp#L382 .Which is used (for a different segmentation) e.g. here: https://github.com/HEP-FCC/k4RecCalorimeter/blob/main/RecFCChhCalorimeter/src/components/CreateFCChhCaloNeighbours.cpp#L92

Cheers,
Brieuc

[atolosadelgado]

Hi again @atolosadelgado , I think this is almost what we want. We're interested in dumping a full list of readout elements, not just readouts. In our example, we have:

<detector> -> <readout> -> <segmentation> -> CartesianGridXY

We're interested in dumping each element of the CartesianGridXY. I'll look for methods within dd4hep for finding this, but feel free to post any suggestions if you have them.

Hi,

sorry, I may have misunderstood the original question. If the question is how to get the size of each pixel/segment, please find below an example that shows the case of cartesian XY segmentation of TGeoBBox. i do not know exactly how to extend it in an elegant way because segmentation can be fully custom

// foo.cxx
// Root script that walks over the tree checking for sensitive volumes and looking for the corresponding segmentation
// To run it: root foo.cxx

#include "DD4hep/Detector.h"
#include "DD4hep/DD4hepUnits.h"
#include "DD4hep/BitFieldCoder.h"
R__LOAD_LIBRARY(libDDCore) 

constexpr int myverbose = 0;

// Function to recursively check sensitive volumes
void checkSensitiveVolumes(dd4hep::VolumeManager& vm, dd4hep::PlacedVolume pv, int level = 0) {

    dd4hep::Volume volume = pv.volume();
    if( 4 < myverbose ) std::cout << volume.name() << std::endl;
    
    // Check if the volume is sensitive
    dd4hep::SensitiveDetector sd = volume.sensitiveDetector();
    if (sd.isValid()) 
    {
    
        if( 3 < myverbose ) std::cout << volume.name() << '\t' << sd.name() << std::endl;
        
        // Get the readout associated with the sensitive volume
        dd4hep::Readout readout = sd.readout();

        // Print the name of the readout
        if( 2 < myverbose ) 
        std::cout << std::string(level, ' ') 
        	  << "Sensitive Volume: " << volume.name()
                  << ", Readout: "        << readout.name() 
                  << ", Segmentation: "   << readout.segmentation().type() 
                  << std::endl;
        dd4hep::Solid s = volume.solid();

        // Print using ROOT Print method
        // s.access()->Print();
        TGeoShape * ss = s.access();
        if( 0 == std::string("TGeoBBox").compare( ss->ClassName() ) )
        {
        auto box_s = static_cast<TGeoBBox*>(ss);
        std::cout << std::string(level, ' ') 
        	  << "Sensitive Volume: "            << volume.name()
                  << ", Readout: "                   << readout.name() 
                  << ", Segmentation type: "         << readout.segmentation().type() 
                  << ", cell x-size / cm: "          << readout.segmentation().cellDimensions(0)[0] /dd4hep::cm
                  << ", Box x-side/cm: "             << 2*box_s->GetDX()/dd4hep::cm
                  << ", cell y-size / cm: "          << readout.segmentation().cellDimensions(0)[1] /dd4hep::cm
                  << ", Box y-side/cm: "             << 2*box_s->GetDY()/dd4hep::cm
                  << std::endl;
        }
        
    }

    // Recursively check child volumes
    int ndaug = pv.num_daughters();
    for (int idaug = 0; idaug<ndaug; ++idaug) {
        checkSensitiveVolumes(vm, pv.daughter(idaug), level + 1);
    }
}

void foo()
{


auto lcdd = &(dd4hep::Detector::getInstance());
lcdd->fromCompact( "lcgeoTests/compact/ARC_standalone_o1_v01.xml" );
dd4hep::VolumeManager vm = lcdd->volumeManager();
vm = vm.getVolumeManager( *lcdd); 
dd4hep::PlacedVolume world_pv = vm.detector().placement();
checkSensitiveVolumes(vm, world_pv);
}

The output would look like this:

     Sensitive Volume: ARCENDCAP_sensor, Readout: ArcCollection, Segmentation type: CartesianGridXY, cell x-size / cm: 0.08, Box x-side/cm: 8, cell y-size / cm: 0.08, Box y-side/cm: 8

[alexandertuna]

Thanks everyone. I believe @vvolkl and @BrieucF understand my request - it is the same as #580 , if I understand correctly.

Thanks for these links, @BrieucF . It seems there is no magic way to list all the readout channels without e.g. calculating cellsEta and minEtaID with ceil/floor. That's good to know. We'll need to think about if we want to repeat this logic ourselves.

[MarkusFrankATcernch]

As you have already noted issue #580 already states that that enumerating all sensitive elements within a sensitive volume in general is impossible. One can enumerate if the volume is of a very concrete shape (such as a box, a cylinder, etc.), but otherwise is very time consuming and not possible in a only half way straight forward way. At the end one would always to sort of scan through an approximation and for each cell check if the edges are "inside" or "outside" the concrete sensitive volume. This is very expensive.

If you can come up with a sort of generic strategy, please let me know! I am happy to implement.