Tidied up in includes, print-outs, and comments (139e75ac) · Commits · XRFLab109 / Allpix Squared

tools/root_analysis_macros/etaCorrectionResiduals.C

+72 −89

Original line number	Diff line number	Diff line
		@@ -14,16 +14,14 @@

		#include <memory>

		#include "../../src/core/utils/unit.h"
		#include "../../src/modules/DetectorHistogrammer/Cluster.hpp"
		#include "../../src/objects/MCParticle.hpp"
		#include "../../src/objects/PixelCharge.hpp"
		#include "../../src/objects/PixelHit.hpp"
		#include "../../src/objects/PropagatedCharge.hpp"

		// Global things, for simplicity
		auto pitch_x = 35 / 1000.0;
		auto pitch_y = 35 / 1000.0;
		// Pixel pitch in µm
		double pitch_x = 35;
		double pitch_y = 35;

		// Initialise histograms
		std::string mod_axes_x = "in-2pixel x_{cluster} [mm];in-2pixel x_{track} [mm];";
		@@ -31,35 +29,38 @@ std::string mod_axes_y = "in-2pixel y_{cluster} [mm];in-2pixel y_{track} [mm];";

		auto etaDistributionX = new TH2F("etaDistributionX",
		std::string("2D #eta distribution X;" + mod_axes_x + "No. entries").c_str(),
		pitch_x * 1000 * 2,
		-pitch_x / 2,
		pitch_x / 2,
		pitch_x * 1000 * 2,
		-pitch_x / 2,
		pitch_x / 2);
		pitch_x * 2,
		-(pitch_x / 1000) / 2,
		(pitch_x / 1000) / 2,
		pitch_x * 2,
		-(pitch_x / 1000) / 2,
		(pitch_x / 1000) / 2);
		auto etaDistributionY = new TH2F("etaDistributionY",
		std::string("2D #eta distribution Y;" + mod_axes_y + "No. entries").c_str(),
		pitch_y * 1000 * 2,
		-pitch_y / 2,
		pitch_y / 2,
		pitch_y * 1000 * 2,
		-pitch_y / 2,
		pitch_y / 2);
		pitch_y * 2,
		-(pitch_y / 1000) / 2,
		(pitch_y / 1000) / 2,
		pitch_y * 2,
		-(pitch_y / 1000) / 2,
		(pitch_y / 1000) / 2);
		auto etaDistributionXprofile = new TProfile("etaDistributionXprofile",
		std::string("#eta distribution X;" + mod_axes_x).c_str(),
		pitch_x * 1000 * 2,
		-pitch_x / 2,
		pitch_x / 2,
		-pitch_x / 2,
		pitch_x / 2);
		pitch_x * 2,
		-(pitch_x / 1000) / 2,
		(pitch_x / 1000) / 2,
		-(pitch_x / 1000) / 2,
		(pitch_x / 1000) / 2);
		auto etaDistributionYprofile = new TProfile("etaDistributionYprofile",
		std::string("#eta distribution Y;" + mod_axes_y).c_str(),
		pitch_y * 1000 * 2,
		-pitch_y / 2,
		pitch_y / 2,
		-pitch_y / 2,
		pitch_y / 2);
		pitch_y * 2,
		-(pitch_y / 1000) / 2,
		(pitch_y / 1000) / 2,
		-(pitch_y / 1000) / 2,
		(pitch_y / 1000) / 2);

		/**
		* Function to perform clustering, returning a vector of Cluster objects
		*/
		std::vector<allpix::Cluster> doClustering(const std::vector<allpix::PixelHit*>& input_hits) {
		std::vector<allpix::Cluster> clusters;
		std::map<const allpix::PixelHit*, bool> usedPixel;
		@@ -76,8 +77,6 @@ std::vector<allpix::Cluster> doClustering(const std::vector<allpix::PixelHit*>&
		// Create new cluster
		allpix::Cluster cluster(pixel_hit);
		usedPixel[pixel_hit] = true;
		std::cout << "Creating new cluster with seed: " << pixel_hit->getPixel().getIndex().X() << ", "
		<< pixel_hit->getPixel().getIndex().Y() << std::endl;

		auto touching = [&](const allpix::PixelHit* pixel) {
		auto pxi1 = pixel->getIndex();
		@@ -103,8 +102,6 @@ std::vector<allpix::Cluster> doClustering(const std::vector<allpix::PixelHit*>&
		}

		cluster.addPixelHit(neighbour);
		std::cout << "Adding pixel: " << neighbour->getPixel().getIndex().X() << ", "
		<< neighbour->getPixel().getIndex().Y() << std::endl;
		usedPixel[neighbour] = true;
		other_pixel = pixel_it;
		}
		@@ -113,6 +110,9 @@ std::vector<allpix::Cluster> doClustering(const std::vector<allpix::PixelHit*>&
		return clusters;
		}

		/**
		* Function to get the primary particles of an input vector of MCParticle pointers
		*/
		std::vector<const allpix::MCParticle> getPrimaryMCParticles(const std::vector<allpix::MCParticle> allMCparticles) {
		std::vector<const allpix::MCParticle*> primaries;

		@@ -129,20 +129,19 @@ std::vector<const allpix::MCParticle*> getPrimaryMCParticles(const std::vector<a
		return primaries;
		}

		/**
		* Fill up the eta distributions, for each cluster
		*/
		void calculate_eta(const allpix::MCParticle* track, allpix::Cluster* cluster) {
		// Ignore single pixel clusters
		if(cluster->getSize() == 1) {
		// std::cout << "Cluster size 1, skipping" << std::endl;
		return;
		}
		std::cout << "Cluster size larger than 1, doing eta correction" << std::endl;

		auto localIntercept = track->getLocalReferencePoint();
		auto [size_x, size_y] = cluster->getSizeXY();
		std::cout << "Size x: " << size_x << ", size y: " << size_y << std::endl;

		if(size_x == 2) {
		// std::cout << "Size in X is 2" << std::endl;
		auto min_center = std::numeric_limits<double>::max(), max_center = std::numeric_limits<double>::lowest();

		for(const auto& pixel : cluster->getPixelHits()) {
		@@ -156,17 +155,12 @@ void calculate_eta(const allpix::MCParticle* track, allpix::Cluster* cluster) {
		}
		}
		auto reference_X = max_center - (max_center - min_center) / 2.;
		// std::cout << "min = " << min_center << std::endl;
		// std::cout << "max = " << max_center << std::endl;
		// std::cout << "ref = " << reference_X << std::endl;
		auto xmod_cluster = cluster->getPosition().X() - reference_X;
		auto xmod_track = localIntercept.X() - reference_X;
		// std::cout << "xmod_cluster = " << xmod_cluster << " xmod_track = " << xmod_track << std::endl;
		etaDistributionX->Fill(xmod_cluster, xmod_track);
		etaDistributionXprofile->Fill(xmod_cluster, xmod_track);
		}
		if(size_y == 2) {
		// std::cout << "Size in Y is 2" << std::endl;
		auto min_center = std::numeric_limits<double>::max(), max_center = std::numeric_limits<double>::lowest();

		for(const auto& pixel : cluster->getPixelHits()) {
		@@ -182,26 +176,31 @@ void calculate_eta(const allpix::MCParticle* track, allpix::Cluster* cluster) {
		auto reference_Y = max_center - (max_center - min_center) / 2.;
		auto ymod_cluster = cluster->getPosition().Y() - reference_Y;
		auto ymod_track = localIntercept.Y() - reference_Y;
		// std::cout << "ymod_cluster = " << ymod_cluster << " ymod_track = " << ymod_track << std::endl;
		etaDistributionY->Fill(ymod_cluster, ymod_track);
		etaDistributionYprofile->Fill(ymod_cluster, ymod_track);
		}
		}

		// Give the function, the fit range (min and max), and the TProfile to fit
		TF1* fit(const std::string& fname, double minRange, double maxRange, TProfile* profile) {
		/**
		* Fit a fifth order polynomial to the eta distributions.
		* Takes the function name, the fit range (min and max), and the TProfile to be fit as input
		*/
		TF1* fitEta(const std::string& fname, double minRange, double maxRange, TProfile* profile) {

		auto formula = "[0] + [1]x + [2]x^2 + [3]x^3 + [4]x^4 + [5]*x^5";
		auto function = new TF1(fname.c_str(), formula, minRange, maxRange);

		// Get the eta distribution profiles and fit them to extract the correction parameters
		profile->Fit(function, "q R");
		TF1* fit = profile->GetFunction(fname.c_str());
		TF1* fitResult = profile->GetFunction(fname.c_str());

		return fit;
		return fitResult;
		}

		// Applies correction, and returns the updated cluster position
		/**
		* Apply the eta correction. Update cluster position.
		* Returns the updated cluster position.
		*/
		ROOT::Math::XYZPoint applyEtaCorrection(const allpix::Cluster& cluster, TF1* eta_corrector_x, TF1* eta_corrector_y) {
		// Store old cluster position
		auto old_position = cluster.getPosition();
		@@ -251,12 +250,10 @@ ROOT::Math::XYZPoint applyEtaCorrection(const allpix::Cluster& cluster, TF1* eta
		}

		/**
		* An example of how to perform an eta correction of residuals, by fitting of a fifth order polynomial
		* An example of how to perform an eta correction of residuals, by fitting of a fifth order polynomial.
		* Takes a TFile (containing MCParticle objects and PixelHit objects) and a detector name as input.
		*/
		// void etaCorrectionResiduals(TFile* file, std::string detector) {
		void etaCorrectionResiduals(std::string fileName, std::string detector) {

		auto file = new TFile(fileName.c_str());
		void etaCorrectionResiduals(TFile* file, std::string detector) {

		// Initialise reading of the PixelHit TTrees
		TTree* pixel_hit_tree = static_cast<TTree*>(file->Get("PixelHit"));
		@@ -292,82 +289,77 @@ void etaCorrectionResiduals(std::string fileName, std::string detector) {

		// Define histograms
		int noOfResBinsXY = 200;
		double resHalfRange = 50.0 / 1000.0;
		double resHalfRangeX = pitch_x;
		double resHalfRangeY = pitch_y;
		// Residuals from clustering; in micrometres!
		TH1D* residual_x = new TH1D("residual_x",
		"Residual x; x_{MC} - x_{cluster} [#mum]; hits",
		noOfResBinsXY,
		-resHalfRange * 1000,
		resHalfRange * 1000);
		TH1D* residual_y = new TH1D("residual_y",
		"Residual y; y_{MC} - y_{cluster} [#mum]; hits",
		noOfResBinsXY,
		-resHalfRange * 1000,
		resHalfRange * 1000);
		TH1D* residual_x = new TH1D(
		"residual_x", "Residual x; x_{MC} - x_{cluster} [#mum]; hits", noOfResBinsXY, -resHalfRangeX, resHalfRangeX);
		TH1D* residual_y = new TH1D(
		"residual_y", "Residual y; y_{MC} - y_{cluster} [#mum]; hits", noOfResBinsXY, -resHalfRangeY, resHalfRangeY);
		TH1D* residual_x_corrected = new TH1D("residual_x_corrected",
		"Residual x, corrected; x_{MC} - x_{cluster} [#mum]; hits",
		noOfResBinsXY,
		-resHalfRange * 1000,
		resHalfRange * 1000);
		-resHalfRangeX,
		resHalfRangeX);
		TH1D* residual_y_corrected = new TH1D("residual_y_corrected",
		"Residual y, corrected; y_{MC} - y_{cluster} [#mum]; hits",
		noOfResBinsXY,
		-resHalfRange * 1000,
		resHalfRange * 1000);
		-resHalfRangeY,
		resHalfRangeY);
		// Mean absolute deviation, i.e. absolute value of residual, vs pixel position.
		TProfile* residual_x_vs_x =
		new TProfile("residual_x_vs_x",
		"Mean absolute deviation in x, vs in-pixel position in x; x [#mum]; \|x_{MC} - x_{cluster}\| [#mum]",
		noOfResBinsXY,
		0,
		pitch_x * 1000);
		pitch_x);
		TProfile* residual_y_vs_y =
		new TProfile("residual_y_vs_y",
		"Mean absolute deviation in y, vs in-pixel position in y; y [#mum]; \|y_{MC} - y_{cluster}\| [#mum]",
		noOfResBinsXY,
		0,
		pitch_y * 1000);
		pitch_y);
		// Mean absolute deviation, i.e. absolute value of residual, vs pixel position.
		TProfile* residual_x_vs_y =
		new TProfile("residual_x_vs_y",
		"Mean absolute deviation in x, vs in-pixel position in y; y [#mum]; \|x_{MC} - x_{cluster}\| [#mum]",
		noOfResBinsXY,
		0,
		pitch_y * 1000);
		pitch_y);
		TProfile* residual_y_vs_x =
		new TProfile("residual_y_vs_x",
		"Mean absolute deviation in y, vs in-pixel position in x; x [#mum]; \|y_{MC} - y_{cluster}\| [#mum]",
		noOfResBinsXY,
		0,
		pitch_x * 1000);
		pitch_x);
		// In 2D map form
		TProfile2D* residual_map_full =
		new TProfile2D("residual_map_full",
		"Mean 2D residual vs particle hit position; x [#mum]; y [#mum]; Residual [#mum]",
		noOfResBinsXY,
		0,
		pitch_x * 1000,
		pitch_x,
		noOfResBinsXY,
		0,
		pitch_y * 1000);
		pitch_y);
		TProfile2D* residual_map_x =
		new TProfile2D("residual_map_x",
		"Mean absolute deviation in x vs particle hit position; x [#mum]; y [#mum]; MAD_x [#mum]",
		noOfResBinsXY,
		0,
		pitch_x * 1000,
		pitch_x,
		noOfResBinsXY,
		0,
		pitch_y * 1000);
		pitch_y);
		TProfile2D* residual_map_y =
		new TProfile2D("residual_map_y",
		"Mean absolute deviation in y vs particle hit position; x [#mum]; y [#mum]; MAD_y [#mum]",
		noOfResBinsXY,
		0,
		pitch_x * 1000,
		pitch_x,
		noOfResBinsXY,
		0,
		pitch_y * 1000);
		pitch_y);

		// Iterate over all events, for eta calculation
		for(int i = 0; i < pixel_hit_tree->GetEntries(); ++i) {
		@@ -379,14 +371,10 @@ void etaCorrectionResiduals(std::string fileName, std::string detector) {
		pixel_hit_tree->GetEntry(i);
		mc_particle_tree->GetEntry(i);

		if(input_hits.size() > 1) {
		std::cout << "Pixel index: " << input_hits[1]->getPixel().getIndex().x() << ", "
		<< input_hits[0]->getPixel().getIndex().y() << std::endl;
		}
		// Perform clustering
		std::vector<allpix::Cluster> clusters = doClustering(input_hits);

		// Loop over all tracks and look at the associated clusters to plot the eta distribution
		// Loop over all MC particles and look at the associated clusters to plot the eta distribution
		for(auto& cluster : clusters) {
		auto mcp = cluster.getMCParticles();
		// Calculate eta for first primary particle:
		@@ -400,8 +388,8 @@ void etaCorrectionResiduals(std::string fileName, std::string detector) {
		}

		// Fit and save the eta correction functions for x and y
		auto eta_corrector_x = fit("eta_corrector_x", -pitch_x / 2, pitch_x / 2, etaDistributionXprofile);
		auto eta_corrector_y = fit("eta_corrector_y", -pitch_y / 2, pitch_y / 2, etaDistributionYprofile);
		auto eta_corrector_x = fitEta("eta_corrector_x", -pitch_x / 2, pitch_x / 2, etaDistributionXprofile);
		auto eta_corrector_y = fitEta("eta_corrector_y", -pitch_y / 2, pitch_y / 2, etaDistributionYprofile);

		// Iterate over all events, to draw the rest of the owl
		for(int i = 0; i < pixel_hit_tree->GetEntries(); ++i) {
		@@ -418,14 +406,10 @@ void etaCorrectionResiduals(std::string fileName, std::string detector) {

		// Get all the primary MC particles for the current event
		std::vector<const allpix::MCParticle*> primaryMCparticles = getPrimaryMCParticles(input_particles);
		// std::cout << "Number of MC particles: " << input_particles.size() << " Number of primaries: " <<
		// primaryMCparticles.size() << std::endl;

		for(const auto& cluster : clusters) {

		// Find primaries conencted with cluster (should be the only one we have for these single-particle events, but
		// still) Might be zero though, if the cluster has been created by a secondary particle. So; not all clusters are
		// connected to primary particles. Delta rays, in this case
		// Find primaries associated with cluster
		std::set<const allpix::MCParticle*> cluster_particles = cluster.getMCParticles();
		std::vector<const allpix::MCParticle*> intersection;
		std::set_intersection(primaryMCparticles.begin(),
		@@ -435,7 +419,6 @@ void etaCorrectionResiduals(std::string fileName, std::string detector) {
		std::back_inserter(intersection));

		if(intersection.size() == 0) {
		// std::cout << "Cluster from secondary MC particle!" << std::endl;
		continue;
		}
		// NOTE: anything below this will only have clusters from primary particles