Merge branch 'dethist_pixelcenter' into 'master' (80c5c100) · Commits · XRFLab109 / Allpix Squared

src/modules/DetectorHistogrammer/DetectorHistogrammerModule.cpp

+83 −46

Original line number	Diff line number	Diff line
		@@ -89,30 +89,36 @@ void DetectorHistogrammerModule::initialize() {
		// Create histogram of cluster map
		std::string cluster_size_map_title = "Cluster size as function of in-pixel impact position for " + detector_->getName() +
		";x%pitch [#mum];y%pitch [#mum]";
		cluster_size_map = CreateHistogram<TProfile2D>(
		"cluster_size_map", cluster_size_map_title.c_str(), inpixel_bins.x(), 0., pitch_x, inpixel_bins.y(), 0., pitch_y);
		cluster_size_map = CreateHistogram<TProfile2D>("cluster_size_map",
		cluster_size_map_title.c_str(),
		inpixel_bins.x(),
		-pitch_x / 2,
		pitch_x / 2,
		inpixel_bins.y(),
		-pitch_y / 2,
		pitch_y / 2);

		std::string cluster_size_x_map_title = "Cluster size in X as function of in-pixel impact position for " +
		detector_->getName() + ";x%pitch [#mum];y%pitch [#mum]";
		cluster_size_x_map = CreateHistogram<TProfile2D>("cluster_size_x_map",
		cluster_size_x_map_title.c_str(),
		inpixel_bins.x(),
		0.,
		pitch_x,
		-pitch_x / 2,
		pitch_x / 2,
		inpixel_bins.y(),
		0.,
		pitch_y);
		-pitch_y / 2,
		pitch_y / 2);

		std::string cluster_size_y_map_title = "Cluster size in Y as function of in-pixel impact position for " +
		detector_->getName() + ";x%pitch [#mum];y%pitch [#mum]";
		cluster_size_y_map = CreateHistogram<TProfile2D>("cluster_size_y_map",
		cluster_size_y_map_title.c_str(),
		inpixel_bins.x(),
		0.,
		pitch_x,
		-pitch_x / 2,
		pitch_x / 2,
		inpixel_bins.y(),
		0.,
		pitch_y);
		-pitch_y / 2,
		pitch_y / 2);

		// Charge maps:
		std::string cluster_charge_map_title = "Cluster charge as function of in-pixel impact position for " +
		@@ -120,15 +126,21 @@ void DetectorHistogrammerModule::initialize() {
		cluster_charge_map = CreateHistogram<TProfile2D>("cluster_charge_map",
		cluster_charge_map_title.c_str(),
		inpixel_bins.x(),
		0.,
		pitch_x,
		-pitch_x / 2,
		pitch_x / 2,
		inpixel_bins.y(),
		0.,
		pitch_y);
		-pitch_y / 2,
		pitch_y / 2);
		std::string seed_charge_map_title = "Seed pixel charge as function of in-pixel impact position for " +
		detector_->getName() + ";x%pitch [#mum];y%pitch [#mum];<seed pixel charge> [ke]";
		seed_charge_map = CreateHistogram<TProfile2D>(
		"seed_charge_map", seed_charge_map_title.c_str(), inpixel_bins.x(), 0., pitch_x, inpixel_bins.y(), 0., pitch_y);
		seed_charge_map = CreateHistogram<TProfile2D>("seed_charge_map",
		seed_charge_map_title.c_str(),
		inpixel_bins.x(),
		-pitch_x / 2,
		pitch_x / 2,
		inpixel_bins.y(),
		-pitch_y / 2,
		pitch_y / 2);

		// Create cluster size plots, preventing a zero-bin histogram by scaling with integer ceiling: (x + y - 1) / y
		std::string cluster_size_title = "Cluster size for " + detector_->getName() + ";cluster size [px];clusters";
		@@ -157,27 +169,33 @@ void DetectorHistogrammerModule::initialize() {
		// Residual projections
		std::string residual_x_vs_x_title = "Mean absolute deviation of residual in X as function of in-pixel X position for " +
		detector_->getName() + ";x%pitch [#mum];MAD(#Deltax) [#mum]";
		residual_x_vs_x =
		CreateHistogram<TProfile>("residual_x_vs_x", residual_x_vs_x_title.c_str(), inpixel_bins.x(), 0., pitch_x);
		residual_x_vs_x = CreateHistogram<TProfile>(
		"residual_x_vs_x", residual_x_vs_x_title.c_str(), inpixel_bins.x(), -pitch_x / 2, pitch_x / 2);
		std::string residual_y_vs_y_title = "Mean absolute deviation of residual in Y as function of in-pixel Y position for " +
		detector_->getName() + ";y%pitch [#mum];MAD(#Deltay) [#mum]";
		residual_y_vs_y =
		CreateHistogram<TProfile>("residual_y_vs_y", residual_y_vs_y_title.c_str(), inpixel_bins.y(), 0., pitch_y);
		residual_y_vs_y = CreateHistogram<TProfile>(
		"residual_y_vs_y", residual_y_vs_y_title.c_str(), inpixel_bins.y(), -pitch_y / 2, pitch_y / 2);
		std::string residual_x_vs_y_title = "Mean absolute deviation of residual in X as function of in-pixel Y position for " +
		detector_->getName() + ";y%pitch [#mum];MAD(#Deltax) [#mum]";
		residual_x_vs_y =
		CreateHistogram<TProfile>("residual_x_vs_y", residual_x_vs_y_title.c_str(), inpixel_bins.y(), 0., pitch_y);
		residual_x_vs_y = CreateHistogram<TProfile>(
		"residual_x_vs_y", residual_x_vs_y_title.c_str(), inpixel_bins.y(), -pitch_y / 2, pitch_y / 2);
		std::string residual_y_vs_x_title = "Mean absolute deviation of residual in Y as function of in-pixel X position for " +
		detector_->getName() + ";x%pitch [#mum];MAD(#Deltay) [#mum]";
		residual_y_vs_x =
		CreateHistogram<TProfile>("residual_y_vs_x", residual_y_vs_x_title.c_str(), inpixel_bins.x(), 0., pitch_x);
		residual_y_vs_x = CreateHistogram<TProfile>(
		"residual_y_vs_x", residual_y_vs_x_title.c_str(), inpixel_bins.x(), -pitch_x / 2, pitch_x / 2);

		// Residual maps
		std::string residual_map_title = "Mean absolute deviation of residual as function of in-pixel impact position for " +
		detector_->getName() +
		";x%pitch [#mum];y%pitch [#mum];MAD(#sqrt{#Deltax^{2}+#Deltay^{2}}) [#mum]";
		residual_map = CreateHistogram<TProfile2D>(
		"residual_map", residual_map_title.c_str(), inpixel_bins.x(), 0., pitch_x, inpixel_bins.y(), 0., pitch_y);
		residual_map = CreateHistogram<TProfile2D>("residual_map",
		residual_map_title.c_str(),
		inpixel_bins.x(),
		-pitch_x / 2,
		pitch_x / 2,
		inpixel_bins.y(),
		-pitch_y / 2,
		pitch_y / 2);
		std::string residual_detector_title = "Mean absolute deviation of residual of " + detector_->getName() +
		";x (pixels);y (pixels);MAD(#sqrt{#Deltax^{2}+#Deltay^{2}}) [#mum]";
		residual_detector = CreateHistogram<TProfile2D>(
		@@ -186,8 +204,14 @@ void DetectorHistogrammerModule::initialize() {
		std::string residual_x_map_title =
		"Mean absolute deviation of residual in X as function of in-pixel impact position for " + detector_->getName() +
		";x%pitch [#mum];y%pitch [#mum];MAD(#Deltax) [#mum]";
		residual_x_map = CreateHistogram<TProfile2D>(
		"residual_x_map", residual_x_map_title.c_str(), inpixel_bins.x(), 0., pitch_x, inpixel_bins.y(), 0., pitch_y);
		residual_x_map = CreateHistogram<TProfile2D>("residual_x_map",
		residual_x_map_title.c_str(),
		inpixel_bins.x(),
		-pitch_x / 2,
		pitch_x / 2,
		inpixel_bins.y(),
		-pitch_y / 2,
		pitch_y / 2);
		std::string residual_x_detector_title =
		"Mean absolute deviation of residual in X of " + detector_->getName() + ";x (pixels);y (pixels);MAD(#Deltax) [#mum]";
		residual_x_detector = CreateHistogram<TProfile2D>("residual_x_detector",
		@@ -202,8 +226,14 @@ void DetectorHistogrammerModule::initialize() {
		std::string residual_y_map_title =
		"Mean absolute deviation of residual in Y as function of in-pixel impact position for " + detector_->getName() +
		";x%pitch [#mum];y%pitch [#mum];MAD(#Deltay) [#mum]";
		residual_y_map = CreateHistogram<TProfile2D>(
		"residual_y_map", residual_y_map_title.c_str(), inpixel_bins.x(), 0., pitch_x, inpixel_bins.y(), 0., pitch_y);
		residual_y_map = CreateHistogram<TProfile2D>("residual_y_map",
		residual_y_map_title.c_str(),
		inpixel_bins.x(),
		-pitch_x / 2,
		pitch_x / 2,
		inpixel_bins.y(),
		-pitch_y / 2,
		pitch_y / 2);
		std::string residual_y_detector_title =
		"Mean absolute deviation of residual in Y of " + detector_->getName() + ";x (pixels);y (pixels);MAD(#Deltay) [#mum]";
		residual_y_detector = CreateHistogram<TProfile2D>("residual_y_detector",
		@@ -218,8 +248,16 @@ void DetectorHistogrammerModule::initialize() {
		// Efficiency maps:
		std::string efficiency_map_title = "Efficiency as function of in-pixel impact position for " + detector_->getName() +
		";x%pitch [#mum];y%pitch [#mum];efficiency";
		efficiency_map = CreateHistogram<TProfile2D>(
		"efficiency_map", efficiency_map_title.c_str(), inpixel_bins.x(), 0, pitch_x, inpixel_bins.y(), 0, pitch_y, 0, 1);
		efficiency_map = CreateHistogram<TProfile2D>("efficiency_map",
		efficiency_map_title.c_str(),
		inpixel_bins.x(),
		-pitch_x / 2,
		pitch_x / 2,
		inpixel_bins.y(),
		-pitch_y / 2,
		pitch_y / 2,
		0,
		1);
		std::string efficiency_detector_title = "Efficiency of " + detector_->getName() + ";x (pixels);y (pixels);efficiency";
		efficiency_detector = CreateHistogram<TProfile2D>("efficiency_detector",
		efficiency_detector_title.c_str(),
		@@ -234,12 +272,12 @@ void DetectorHistogrammerModule::initialize() {
		// Efficiency projections
		std::string efficiency_vs_x_title =
		"Efficiency as function of in-pixel X position for " + detector_->getName() + ";x%pitch [#mum];efficiency";
		efficiency_vs_x =
		CreateHistogram<TProfile>("efficiency_vs_x", efficiency_vs_x_title.c_str(), inpixel_bins.x(), 0., pitch_x, 0, 1);
		efficiency_vs_x = CreateHistogram<TProfile>(
		"efficiency_vs_x", efficiency_vs_x_title.c_str(), inpixel_bins.x(), -pitch_x / 2, pitch_x / 2, 0, 1);
		std::string efficiency_vs_y_title =
		"Efficiency as function of in-pixel Y position for " + detector_->getName() + ";y%pitch [#mum];efficiency";
		efficiency_vs_y =
		CreateHistogram<TProfile>("efficiency_vs_y", efficiency_vs_y_title.c_str(), inpixel_bins.y(), 0., pitch_y, 0, 1);
		efficiency_vs_y = CreateHistogram<TProfile>(
		"efficiency_vs_y", efficiency_vs_y_title.c_str(), inpixel_bins.y(), -pitch_y / 2, pitch_y / 2, 0, 1);

		// Create number of clusters plot
		std::string n_cluster_title = "Number of clusters for " + detector_->getName() + ";clusters;events";
		@@ -333,10 +371,11 @@ void DetectorHistogrammerModule::run(Event* event) {
		auto particlePos = particle->getLocalReferencePoint() + track_smearing(track_resolution_);
		LOG(DEBUG) << "MCParticle at " << Units::display(particlePos, {"mm", "um"});

		auto inPixelPos = XYVector(std::fmod(particlePos.x() + pitch.x() / 2, pitch.x()),
		std::fmod(particlePos.y() + pitch.y() / 2, pitch.y()));
		LOG(TRACE) << "MCParticle in pixel at " << Units::display(inPixelPos, {"mm", "um"});
		// Find the nearest pixel
		auto [xpixel, ypixel] = detector_->getModel()->getPixelIndex(particlePos);

		auto inPixelPos = particlePos - detector_->getModel()->getPixelCenter(static_cast<unsigned int>(xpixel),
		static_cast<unsigned int>(ypixel));
		auto inPixel_um_x = static_cast<double>(Units::convert(inPixelPos.x(), "um"));
		auto inPixel_um_y = static_cast<double>(Units::convert(inPixelPos.y(), "um"));

		@@ -348,9 +387,6 @@ void DetectorHistogrammerModule::run(Event* event) {
		cluster_charge_map->Fill(
		inPixel_um_x, inPixel_um_y, static_cast<double>(Units::convert(clus.getCharge(), "ke")));

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

		// Retrieve the pixel to which this MCParticle points:
		const auto* pixel = clus.getPixelHit(static_cast<unsigned int>(xpixel), static_cast<unsigned int>(ypixel));
		if(pixel != nullptr) {
		@@ -387,14 +423,15 @@ void DetectorHistogrammerModule::run(Event* event) {

		// Calculate 2D local position of particle:
		auto particlePos = particle->getLocalReferencePoint() + track_smearing(track_resolution_);
		auto inPixelPos = XYVector(std::fmod(particlePos.x() + pitch.x() / 2, pitch.x()),
		std::fmod(particlePos.y() + pitch.y() / 2, pitch.y()));
		auto inPixel_um_x = static_cast<double>(Units::convert(inPixelPos.x(), "um"));
		auto inPixel_um_y = static_cast<double>(Units::convert(inPixelPos.y(), "um"));

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

		auto inPixelPos = particlePos - detector_->getModel()->getPixelCenter(static_cast<unsigned int>(xpixel),
		static_cast<unsigned int>(ypixel));
		auto inPixel_um_x = static_cast<double>(Units::convert(inPixelPos.x(), "um"));
		auto inPixel_um_y = static_cast<double>(Units::convert(inPixelPos.y(), "um"));

		auto matched_cluster =
		std::find_if(clusters.begin(), clusters.end(), [this, &particlePos, &pitch](const Cluster& clus) {
		return (std::fabs(clus.getPosition().x() * pitch.x() - particlePos.x()) < matching_cut_.x()) &&