Merge branch 'histogrammerEfficiencyFix' into 'master' (14176e64) · Commits · XRFLab109 / Allpix Squared

src/core/geometry/DetectorModel.hpp

+10 −0

Original line number	Diff line number	Diff line
		@@ -370,6 +370,16 @@ namespace allpix {
		*/
		virtual bool isWithinMatrix(const int x, const int y) const = 0;

		/**
		* @brief Returns if a position is within the grid of pixels defined for the device
		* @param position Position in local coordinates of the detector model
		* @return True if position within the pixel grid, false otherwise
		*/
		bool isWithinMatrix(const ROOT::Math::XYZPoint& position) const {
		auto [index_x, index_y] = getPixelIndex(position);
		return isWithinMatrix(index_x, index_y);
		}

		/**
		* @brief Returns a pixel center in local coordinates
		* @param x X- (or column-) coordinate of the pixel

src/modules/DetectorHistogrammer/DetectorHistogrammerModule.cpp

+20 −9

Original line number	Diff line number	Diff line
		@@ -40,7 +40,7 @@ DetectorHistogrammerModule::DetectorHistogrammerModule(Configuration& config,

		auto model = detector_->getModel();
		config_.setDefault<XYVector>("matching_cut", model->getPixelSize() * 3);
		config_.setDefault<XYVector>("track_resolution", ROOT::Math::XYVector(Units::get(2.0, "um"), Units::get(2.0, "um")));
		config_.setDefault<XYVector>("track_resolution", ROOT::Math::XYVector(Units::get(0.0, "um"), Units::get(0.0, "um")));

		config_.setDefault<DisplacementVector2D<Cartesian2D<int>>>(
		"granularity",
		@@ -397,7 +397,10 @@ void DetectorHistogrammerModule::run(Event* event) {
		std::shared_ptr<PixelHitMessage> pixels_message{nullptr};
		auto mcparticle_message = messenger_->fetchMessage<MCParticleMessage>(this, event);

		auto radial_model = std::dynamic_pointer_cast<RadialStripDetectorModel>(detector_->getModel());
		// Fetch detector model
		auto model = detector_->getModel();

		auto radial_model = std::dynamic_pointer_cast<RadialStripDetectorModel>(model);

		// Check that we actually received pixel hits - we might have none and just received MCParticles!
		try {
		@@ -452,7 +455,7 @@ void DetectorHistogrammerModule::run(Event* event) {
		cluster_size_y->Fill(clusSizesXY.second);

		auto clusterPos = clus.getPosition();
		auto [cluster_x, cluster_y] = detector_->getModel()->getPixelIndex(clusterPos);
		auto [cluster_x, cluster_y] = model->getPixelIndex(clusterPos);
		LOG(DEBUG) << "Cluster at indices " << cluster_x << ", " << cluster_y << "(" << clusterPos
		<< " local coordinates) with charge " << Units::display(clus.getCharge(), "ke");
		cluster_map->Fill(cluster_x, cluster_y);
		@@ -480,9 +483,9 @@ void DetectorHistogrammerModule::run(Event* event) {
		LOG(DEBUG) << "MCParticle at " << Units::display(particlePos, {"mm", "um"});

		// Find the nearest pixel
		auto [xpixel, ypixel] = detector_->getModel()->getPixelIndex(particlePos);
		auto [xpixel, ypixel] = model->getPixelIndex(particlePos);

		auto inPixelPos = particlePos - detector_->getModel()->getPixelCenter(xpixel, ypixel);
		auto inPixelPos = particlePos - model->getPixelCenter(xpixel, ypixel);
		LOG(TRACE) << "MCParticle in pixel at " << Units::display(inPixelPos, {"mm", "um"});

		// Calculate residual with cluster position:
		@@ -492,7 +495,7 @@ void DetectorHistogrammerModule::run(Event* event) {
		// If model is radial_strip, calculate polar residuals and in-pixel positions
		if(radial_model != nullptr) {
		// Transform coordinates to polar representation
		auto strip_polar = radial_model->getPositionPolar(detector_->getModel()->getPixelCenter(xpixel, ypixel));
		auto strip_polar = radial_model->getPositionPolar(model->getPixelCenter(xpixel, ypixel));
		auto particle_polar = radial_model->getPositionPolar(particlePos);
		auto cluster_polar = radial_model->getPositionPolar(clusterPos);

		@@ -551,15 +554,23 @@ void DetectorHistogrammerModule::run(Event* event) {
		// Calculate 2D local position of particle:
		auto particlePos = particle->getLocalReferencePoint() + track_smearing(track_resolution_);

		// Check whether the particle position is in the sensor excess, and exclude it from the efficiency calculation if so
		if(!model->isWithinMatrix(particlePos)) {
		LOG(DEBUG) << "Particle at local coordinate " << Units::display(particlePos, {"mm", "um"}) << ", pixel index ("
		<< model->getPixelIndex(particlePos).first << "," << model->getPixelIndex(particlePos).second
		<< "), hit in the sensor excess; removing from efficiency calculation.";
		continue;
		}

		// Find the nearest pixel
		auto [xpixel, ypixel] = detector_->getModel()->getPixelIndex(particlePos);
		auto [xpixel, ypixel] = model->getPixelIndex(particlePos);

		auto inPixelPos = particlePos - detector_->getModel()->getPixelCenter(xpixel, ypixel);
		auto inPixelPos = particlePos - model->getPixelCenter(xpixel, ypixel);

		// If model is radial_strip, recalculate inPixelPos
		if(radial_model != nullptr) {
		// Transform coordinates to polar representation
		auto strip_polar = radial_model->getPositionPolar(detector_->getModel()->getPixelCenter(xpixel, ypixel));
		auto strip_polar = radial_model->getPositionPolar(model->getPixelCenter(xpixel, ypixel));
		auto particle_polar = radial_model->getPositionPolar(particlePos);

		// Overwrite inPixelPos with correct values

src/modules/DetectorHistogrammer/README.md

+3 −3

Original line number	Diff line number	Diff line
		@@ -19,8 +19,8 @@ If the PixelHit is free-standing, a new cluster is created.

		This module serves as a quick "mini-analysis" and creates the histograms listed below.
		The Monte Carlo truth position provided by the `MCParticle` objects is used as track reference position.
		An additional uncertainty can be added by configuring a track resolution, with which every cluster residual is convolved.
		For technical reasons, this offset is drawn randomly from a Gauss distribution independently for the resolution and the efficiency measurement.
		An additional uncertainty can be added by configuring a track resolution, with which every cluster residual is convolved. This makes it possible to perform a quick test beam-like analysis.
		For technical reasons, this offset is drawn randomly from a Gaussian distribution independently for the resolution and the efficiency measurement. Note: If a non-zero track resolution is used, pixel matrix edge effects may appear as particles hit the sensor excess.

		* A hitmap of all pixels in the pixel grid, displaying the number of times a pixel has been hit during the simulation run.
		* A cluster map indicating the cluster positions for the whole simulation run.
		@@ -42,7 +42,7 @@ For technical reasons, this offset is drawn randomly from a Gauss distribution i
		* `granularity`: 2D integer vector defining the number of bins along the x and y axis for in-pixel maps. Defaults to the pixel pitch in micro meters, e.g. a detector with 100um x 100um pixels would be represented in a histogram with `100 * 100 = 10000` bins.
		* `granularity_local`: 2D integer vector defining the number of bins for each pixel along the x and y axis for maps in local coordinates where particle positions are used as reference. Defaults to `1 1` corresponding to a single bin per pixel.
		* `max_cluster_charge`: Upper limit for the cluster charge histogram, defaults to `50ke`.
		* `track_resolution`: Assumed track resolution the Monte Carlo truth is smeared with. Expects two values for the resolution in local-x and local-y directions and defaults to `2um 2um`.
		* `track_resolution`: Assumed track resolution the Monte Carlo truth is smeared with. Expects two values for the resolution in local-x and local-y directions and defaults to `0um 0um`, i.e. no smearing.
		* `matching_cut`: Required maximum matching distance between cluster position and particle position for the efficiency measurement. Expected two values and defaults to three times the pixel pitch in each dimension.

		## Usage