Made it work, more or less... Includes still need to be done more sensibly (823d5e91) · Commits · XRFLab109 / Allpix Squared

tools/root_analysis_macros/etaCorrectionResiduals.C

+321 −113

Original line number	Diff line number	Diff line
		@@ -2,19 +2,71 @@
		// SPDX-License-Identifier: MIT

		#include <TCanvas.h>
		#include <TF1.h>
		#include <TFile.h>
		#include <TH1D.h>
		#include <TH2D.h>
		#include <TProfile.h>
		#include <TSystem.h>
		#include <TTree.h>
		#include <vector>

		#include <memory>

		#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"
		//#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"
		//#include "../../src/core/utils/unit.h"

		#include "/home/hawk/apVersions/ap2p4p0/src/objects/MCParticle.hpp"
		#include "/home/hawk/apVersions/ap2p4p0/src/objects/PixelCharge.hpp"
		#include "/home/hawk/apVersions/ap2p4p0/src/objects/PixelHit.hpp"
		#include "/home/hawk/apVersions/ap2p4p0/src/objects/PropagatedCharge.hpp"

		#include "/home/hawk/apVersions/ap2p4p0/src/modules/DetectorHistogrammer/Cluster.hpp"

		#include "/home/hawk/apVersions/ap2p4p0/src/tools/units.h"

		// Global things, for simplicity
		auto pitch_x = 25 / 1000.0;
		auto pitch_y = 25 / 1000.0;

		// Initialise histograms
		std::string mod_axes_x = "in-2pixel x_{cluster} [mm];in-2pixel x_{track} [mm];";
		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);
		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);
		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);
		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);

		std::vector<allpix::Cluster> doClustering(const std::vector<allpix::PixelHit*>& input_hits) {
		std::vector<allpix::Cluster> clusters;
		@@ -69,24 +121,39 @@ std::vector<allpix::Cluster> doClustering(const std::vector<allpix::PixelHit*>&
		return clusters;
		}

		void EtaCalculationModule::calculate_eta(const MCParticle* track, Cluster* cluster) {
		std::vector<const allpix::MCParticle> getPrimaryMCParticles(const std::vector<allpix::MCParticle> allMCparticles) {
		std::vector<const allpix::MCParticle*> primaries;

		// Loop over all MCParticles available
		for(const auto& mc_particle : allMCparticles) {
		// Check for possible parents:
		const auto* parent = mc_particle->getParent();
		if(parent != nullptr) {
		continue;
		}
		// This particle has no parent particles in the regarded sensor, return it.
		primaries.push_back(mc_particle);
		}
		return primaries;
		}

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

		auto model = m_detector->getModel();
		auto localIntercept = track->getLocalReferencePoint();
		auto [size_x, size_y] = cluster->getSizeXY();

		if(size_x == 2) {
		LOG(DEBUG) << "Size in X is 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()) {
		auto idx = pixel->getIndex();
		auto center = model->getPixelCenter(idx.x(), idx.y());
		auto center = pixel->getPixel().getLocalCenter();

		if(center.x() > max_center) {
		max_center = center.x();
		@@ -96,22 +163,21 @@ void EtaCalculationModule::calculate_eta(const MCParticle* track, Cluster* clust
		}
		}
		auto reference_X = max_center - (max_center - min_center) / 2.;
		LOG(DEBUG) << "min = " << min_center;
		LOG(DEBUG) << "max = " << max_center;
		LOG(DEBUG) << "ref = " << reference_X;
		// 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;
		LOG(DEBUG) << "xmod_cluster = " << xmod_cluster << " xmod_track = " << xmod_track;
		m_etaDistributionX->Fill(xmod_cluster, xmod_track);
		m_etaDistributionXprofile->Fill(xmod_cluster, xmod_track);
		// 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) {
		LOG(DEBUG) << "Size in Y is 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()) {
		auto idx = pixel->getIndex();
		auto center = model->getPixelCenter(idx.x(), idx.y());
		auto center = pixel->getPixel().getLocalCenter();

		if(center.y() > max_center) {
		max_center = center.y();
		@@ -123,37 +189,81 @@ void EtaCalculationModule::calculate_eta(const MCParticle* track, Cluster* clust
		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;
		LOG(DEBUG) << "ymod_cluster = " << ymod_cluster << " ymod_track = " << ymod_track;
		m_etaDistributionY->Fill(ymod_cluster, ymod_track);
		m_etaDistributionYprofile->Fill(ymod_cluster, ymod_track);
		// 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 fucntion, the fit range (min and max), and the TProfile to fit
		std::string EtaCalculationModule::fit(const std::string& fname, double minRange, double maxRange, TProfile* profile) {

		auto formula = config_.get<std::string>(fname);
		auto function = new TF1(fname.c_str(), formula.c_str(), minRange, maxRange);
		// 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) {

		if(!function->IsValid()) {
		throw InvalidValueError(config_, fname, "The provided function is not a valid ROOT::TFormula expression");
		}
		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());

		std::stringstream parameters;
		for(int i = 0; i < fit->GetNumberFreeParameters(); i++) {
		parameters << " " << fit->GetParameter(i);
		return fit;
		}
		return parameters.str();

		// Applies correction, and 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();
		// Get cluster size
		auto sizeXY = cluster.getSizeXY();
		// Ignore clusters where neither directional size is 2; return old position
		if(sizeXY.first != 2 && sizeXY.second != 2) {
		return {old_position.x(), old_position.y(), old_position.z()};
		}
		double new_x = old_position.x();
		double new_y = old_position.y();

		if(sizeXY.first == 2) {
		auto min_center = std::numeric_limits<double>::max(), max_center = std::numeric_limits<double>::lowest();
		for(const auto& pixel_hit : cluster.getPixelHits()) {
		auto center = pixel_hit->getPixel().getLocalCenter();
		if(center.x() > max_center) {
		max_center = center.x();
		}
		if(center.x() < min_center) {
		min_center = center.x();
		}
		}
		auto reference_X = max_center - (max_center - min_center) / 2.0;
		auto xmod_cluster = old_position.x() - reference_X;
		new_x = eta_corrector_x->Eval(xmod_cluster) + reference_X;
		}

		if(sizeXY.second == 2) {
		auto min_center = std::numeric_limits<double>::max(), max_center = std::numeric_limits<double>::lowest();
		for(const auto& pixel_hit : cluster.getPixelHits()) {
		auto center = pixel_hit->getPixel().getLocalCenter();
		if(center.y() > max_center) {
		max_center = center.y();
		}
		if(center.y() < min_center) {
		min_center = center.y();
		}
		}
		auto reference_Y = max_center - (max_center - min_center) / 2.0;
		auto ymod_cluster = old_position.y() - reference_Y;
		new_y = eta_corrector_y->Eval(ymod_cluster) + reference_Y;
		}

		ROOT::Math::XYZPoint new_position(new_x, new_y, old_position.z());
		return new_position;
		}

		/**
		* An example of how to perform an eta correction of residuals, by fitting of a fifth order polynomial
		*/
		void etaCorrectionResiduals(TFile* file, std::string detector) {
		// void etaCorrectionResiduals(TFile* file, std::string detector) {
		void etaCorrectionResiduals(std::string fileName, std::string detector) {

		auto file = new TFile(fileName.c_str());

		// Initialise reading of the PixelHit TTrees
		TTree* pixel_hit_tree = static_cast<TTree*>(file->Get("PixelHit"));
		@@ -187,43 +297,101 @@ void etaCorrectionResiduals(TFile* file, std::string detector) {
		std::vector<allpix::MCParticle*> input_particles;
		mc_particle_branch->SetObject(&input_particles);

		// Initialise histograms
		// Hitmap with arbitrary field of view
		TH2D* hitmap = new TH2D("hitmap", "Hitmap; x [mm]; y [mm]; hits", 200, 0, 20, 200, 0, 20);

		// Residuals:
		// first for all hits with the mean position of all MCParticles,
		// then only using the MCParticles that are part of the PixelHit history
		TH1D* residual_x = new TH1D("residual_x", "residual x; x_{MC} - x_{hit} [mm]; hits", 200, -5, 5);
		TH1D* residual_x_related =
		new TH1D("residual_x_related", "residual X, related hits; x_{MC} - x_{hit} [mm]; hits", 200, -5, 5);
		TH1D* residual_y = new TH1D("residual_y", "residual y; y_{MC} - y_{hit} [mm]; hits", 200, -5, 5);
		TH1D* residual_y_related =
		new TH1D("residual_y_related", "residual Y, related hits; y_{MC} - y_{hit} [mm]; hits", 200, -5, 5);

		// Spectrum of the PixelHit signal
		TH1D* spectrum = new TH1D("spectrum", "PixelHit signal spectrum; signal; hits", 200, 0, 100000);

		// Iterate over all events
		// Define histograms
		int noOfResBinsXY = 200;
		double resHalfRange = 50.0 / 1000.0;
		// 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_corrected = new TH1D("residual_x",
		"Residual x, corrected; x_{MC} - x_{cluster} [#mum]; hits",
		noOfResBinsXY,
		-resHalfRange * 1000,
		resHalfRange * 1000);
		TH1D* residual_y_corrected = new TH1D("residual_y",
		"Residual y, corrected; y_{MC} - y_{cluster} [#mum]; hits",
		noOfResBinsXY,
		-resHalfRange * 1000,
		resHalfRange * 1000);
		// 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);
		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);
		// 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);
		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);
		// 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,
		noOfResBinsXY,
		0,
		pitch_y * 1000);
		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,
		noOfResBinsXY,
		0,
		pitch_y * 1000);
		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,
		noOfResBinsXY,
		0,
		pitch_y * 1000);

		// Iterate over all events, for eta calculation
		for(int i = 0; i < pixel_hit_tree->GetEntries(); ++i) {
		if(i % 100 == 0) {
		std::cout << "Processing event " << i << std::endl;
		if(i % 1000 == 0) {
		std::cout << "Processing event " << i << " for eta calculation" << std::endl;
		}

		// Access next event. Pushes information into input_*
		pixel_hit_tree->GetEntry(i);
		mc_particle_tree->GetEntry(i);

		auto pixels_message = messenger_->fetchMessage<PixelHitMessage>(this, event);
		auto mcparticle_message = messenger_->fetchMessage<MCParticleMessage>(this, event);

		// Perform clustering
		std::vector<Cluster> clusters = doClustering(pixels_message);
		std::vector<allpix::Cluster> clusters = doClustering(input_hits);

		// Loop over all tracks 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:
		for(const auto& m : mcp) {
		if(m->isPrimary()) {
		@@ -232,66 +400,106 @@ void etaCorrectionResiduals(TFile* file, std::string detector) {
		}
		}
		}

		// Calculate the mean position of all MCParticles in this event
		double position_mcparts_x = 0.;
		double position_mcparts_y = 0.;
		for(auto& mc_part : input_particles) {
		position_mcparts_x += mc_part->getLocalReferencePoint().x();
		position_mcparts_y += mc_part->getLocalReferencePoint().y();
		}
		position_mcparts_x /= input_particles.size();
		position_mcparts_y /= input_particles.size();

		// Iterate over all PixelHits
		for(auto& hit : input_hits) {
		// Retrieve information using Allpix Squared methods
		double position_hit_x = hit->getPixel().getLocalCenter().x();
		double position_hit_y = hit->getPixel().getLocalCenter().y();
		double charge = hit->getSignal();
		// 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);

		// Access history of the PixelHit
		auto parts = hit->getMCParticles();

		// Calculate the mean position of all MCParticles related to this hit
		double position_mcparts_related_x = 0.;
		double position_mcparts_related_y = 0.;
		for(auto& part : parts) {
		position_mcparts_related_x += part->getLocalReferencePoint().x();
		position_mcparts_related_y += part->getLocalReferencePoint().y();
		// Iterate over all events, to draw the rest of the owl
		for(int i = 0; i < pixel_hit_tree->GetEntries(); ++i) {
		if(i % 1000 == 0) {
		std::cout << "Processing event " << i << " for analysis" << std::endl;
		}
		position_mcparts_related_x /= parts.size();
		position_mcparts_related_y /= parts.size();

		// Fill histograms
		hitmap->Fill(position_hit_x, position_hit_y);

		residual_x->Fill(position_mcparts_x - position_hit_x);
		residual_x_related->Fill(position_mcparts_related_x - position_hit_x);
		residual_y->Fill(position_mcparts_y - position_hit_y);
		residual_y_related->Fill(position_mcparts_related_y - position_hit_y);
		// Access next event. Pushes information into input_*
		pixel_hit_tree->GetEntry(i);
		mc_particle_tree->GetEntry(i);

		spectrum->Fill(charge);
		// Perform clustering
		std::vector<allpix::Cluster> clusters = doClustering(input_hits);

		// 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
		std::set<const allpix::MCParticle*> cluster_particles = cluster.getMCParticles();
		std::vector<const allpix::MCParticle*> intersection;
		std::set_intersection(primaryMCparticles.begin(),
		primaryMCparticles.end(),
		cluster_particles.begin(),
		cluster_particles.end(),
		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
		const allpix::MCParticle* clusterOriginParticle = intersection.at(0);

		ROOT::Math::XYZPoint truthParticlePosition =
		clusterOriginParticle->getLocalReferencePoint(); // Gets the extrapolated position in the middle of the
		// sensor plane, in local coords

		// In um, in pixel coordinate system
		auto inPixelPos = ROOT::Math::XYVector(std::fmod(truthParticlePosition.x() + pitch_x / 2, pitch_x) * 1000,
		std::fmod(truthParticlePosition.y() + pitch_y / 2, pitch_y) * 1000);

		// Units? Cluster position is now in local coordinates, in mm.
		// Then let's also do the whole thing in um. So, another factor of 1000.
		double residualXClustering = (truthParticlePosition.x() - cluster.getPosition().x()) * 1000;
		double residualYClustering = (truthParticlePosition.y() - cluster.getPosition().y()) * 1000;
		residual_x->Fill(residualXClustering);
		residual_y->Fill(residualYClustering);

		residual_x_vs_x->Fill(inPixelPos.x(), fabs(residualXClustering));
		residual_y_vs_y->Fill(inPixelPos.y(), fabs(residualYClustering));
		residual_x_vs_y->Fill(inPixelPos.y(), fabs(residualXClustering));
		residual_y_vs_x->Fill(inPixelPos.x(), fabs(residualYClustering));
		// Residual maps
		residual_map_full->Fill(
		inPixelPos.x(),
		inPixelPos.y(),
		sqrt(residualXClustering * residualXClustering + residualYClustering * residualYClustering));
		residual_map_x->Fill(inPixelPos.x(), inPixelPos.y(), fabs(residualXClustering));
		residual_map_y->Fill(inPixelPos.x(), inPixelPos.y(), fabs(residualYClustering));

		// Apply the eta correction; retrieve the updated cluster position
		auto updatedClusterPosition = applyEtaCorrection(cluster, eta_corrector_x, eta_corrector_y);
		double residualXClustering_corrected = (truthParticlePosition.x() - updatedClusterPosition.x()) * 1000;
		double residualYClustering_corrected = (truthParticlePosition.y() - updatedClusterPosition.y()) * 1000;
		residual_x_corrected->Fill(residualXClustering_corrected);
		residual_y_corrected->Fill(residualYClustering_corrected);
		}
		}

		// Draw histograms
		TCanvas* c0 = new TCanvas("c0", "Hitmap", 600, 400);
		hitmap->Draw("colz");

		TCanvas* c1 = new TCanvas("c1", "Residuals", 1200, 800);
		c1->Divide(2, 2);

		c1->cd(1);
		TCanvas* etaCorrectionCanvas = new TCanvas("etaCorrectionCanvas", "Eta correction", 1200, 800);
		etaCorrectionCanvas->Divide(2, 2);
		etaCorrectionCanvas->cd(1);
		etaDistributionX->Draw("colz");
		etaCorrectionCanvas->cd(2);
		etaDistributionXprofile->Draw();
		etaCorrectionCanvas->cd(3);
		etaDistributionY->Draw("colz");
		etaCorrectionCanvas->cd(4);
		etaDistributionYprofile->Draw();

		TCanvas* residualsCanvas = new TCanvas("residualsCanvas", "Residuals", 1200, 800);
		residualsCanvas->Divide(2, 2);
		residualsCanvas->cd(1);
		residual_x->Draw();
		c1->cd(2);
		residualsCanvas->cd(2);
		residual_x_corrected->Draw();
		residualsCanvas->cd(3);
		residual_y->Draw();

		c1->cd(3);
		residual_x_related->Draw();
		c1->cd(4);
		residual_y_related->Draw();

		TCanvas* c2 = new TCanvas("c2", "Signal spectrum", 600, 400);
		spectrum->Draw();
		residualsCanvas->cd(4);
		residual_y_corrected->Draw();
		}