Working towards MIPs along an arbitrary direction. Added calculation of... (b8cfbca1) · Commits · XRFLab109 / Allpix Squared

src/modules/DepositionPointCharge/DepositionPointChargeModule.cpp

+119 −67

Original line number	Diff line number	Diff line
		@@ -62,7 +62,7 @@ DepositionPointChargeModule::DepositionPointChargeModule(Configuration& config,

		void DepositionPointChargeModule::initialize() {

		auto model = detector_->getModel();
		detector_model_ = detector_->getModel();
		output_plots_ = config_.get<bool>("output_plots");
		output_plots_bins_per_um_ = config_.get<int>("output_plots_bins_per_um");

		@@ -73,7 +73,8 @@ void DepositionPointChargeModule::initialize() {

		// Calculate voxel size and ensure granularity is not zero:
		auto granularity = std::max(config_.get<unsigned int>("number_of_steps"), 1u);
		step_size_z_ = model->getSensorSize().z() / granularity;
		step_size_z_ = detector_model_->getSensorSize().z() / granularity;
		// NOTE: probably want to change this. To go along the axis we deposit along... But how?

		// We should deposit the equivalent of about 80 e/h pairs per micro meter (80`000 per mm):
		auto eh_per_um = config_.get<unsigned int>("number_of_charges");
		@@ -81,6 +82,8 @@ void DepositionPointChargeModule::initialize() {
		LOG(INFO) << "Step size for MIP energy deposition: " << Units::display(step_size_z_, {"um", "mm"}) << ", depositing "
		<< carriers_ << " e/h pairs per step (" << Units::display(eh_per_um, "/um") << ")";

		mip_direction_ = config_.get<ROOT::Math::XYZVector>("mip_direction", ROOT::Math::XYZVector(0.0, 0.0, 1.0)).Unit();

		// Check if the number of charge carriers is larger than zero
		if(carriers_ == 0) {
		throw InvalidValueError(config_,
		@@ -105,7 +108,32 @@ void DepositionPointChargeModule::initialize() {
		scan_y_ = std::find(scan_coordinates_.begin(), scan_coordinates_.end(), "y") != scan_coordinates_.end();
		scan_z_ = std::find(scan_coordinates_.begin(), scan_coordinates_.end(), "z") != scan_coordinates_.end();

		if(scan_coordinates_.size() > 3 \|\| (scan_coordinates_.size() == 3 && !(scan_x_ && scan_y_ && scan_z_))) {
		size_t no_of_coordinates = scan_coordinates_.size();

		// If MIP, and along one of the cardinal directions: don't do a scan in that direction, as the MIP goes along it
		// anyway
		if(scan_x_ && mip_direction_ == ROOT::Math::XYZVector(1.0, 0.0, 0.0)) {
		scan_x_ = false;
		no_of_coordinates--;
		LOG(WARNING) << "MIP shot in the x-direction; don't scan in x";
		}
		if(scan_y_ && mip_direction_ == ROOT::Math::XYZVector(0.0, 1.0, 0.0)) {
		scan_y_ = false;
		no_of_coordinates--;
		LOG(WARNING) << "MIP shot in the y-direction; don't scan in y";
		}
		if(scan_z_ && mip_direction_ == ROOT::Math::XYZVector(0.0, 0.0, 1.0)) {
		scan_z_ = false;
		no_of_coordinates--;
		LOG(WARNING) << "MIP shot in the z-direction; don't scan in z";
		}

		if(no_of_coordinates < 1) {
		LOG(WARNING) << "A scan will not be performed, but all MIPs will be along the same line.";
		}

		if(no_of_coordinates > 3 \|\| !(scan_x_ \|\| scan_y_ \|\| scan_z_) \|\|
		(no_of_coordinates == 3 && !(scan_x_ && scan_y_ && scan_z_))) {
		throw InvalidValueError(
		config_,
		"scan_coordinates",
		@@ -113,90 +141,73 @@ void DepositionPointChargeModule::initialize() {
		}

		// Scan with points required 3D scanning, scan with MIPs only 2D:
		if(scan_coordinates_.size() == 1) {
		root_ = events;
		// Calculate voxel size:
		voxel_ = ROOT::Math::XYZVector(model->getPixelSize().x() / (scan_x_ ? root_ : 1.0),
		model->getPixelSize().y() / (scan_y_ ? root_ : 1.0),
		model->getSensorSize().z() / (scan_z_ ? root_ : 1.0));
		if(!(scan_x_ \|\| scan_y_ \|\| scan_z_)) {
		throw InvalidValueError(config_, "scan_coordinates", "The coordinates must be x, y, or z.");
		}
		} else if(scan_coordinates_.size() == 2) {
		if(no_of_coordinates == 2) {
		// Throw if we don't have a valid combination. Need 2 valid entries; x y, x z, or y z
		root_ = static_cast<unsigned int>(std::lround(std::sqrt(events)));
		if(events != root_ * root_) {
		LOG(WARNING) << "Number of events is not a square, pixel cell volume cannot fully be covered in scan. "
		<< "Closest square is " << root_ * root_;
		}
		voxel_ = ROOT::Math::XYZVector(model->getPixelSize().x() / (scan_x_ ? root_ : 1.0),
		model->getPixelSize().y() / (scan_y_ ? root_ : 1.0),
		model->getSensorSize().z() / (scan_z_ ? root_ : 1.0));
		if(!((scan_x_ && scan_y_) \|\| (scan_x_ && scan_z_) \|\| (scan_y_ && scan_z_))) {
		throw InvalidValueError(config_,
		"scan_coordinates",
		"The coordinates must be x, y, or z, and a coordinate must not be repeated");
		}
		} else if(scan_coordinates_.size() == 3 && type_ == SourceType::MIP) {
		root_ = static_cast<unsigned int>(std::lround(std::sqrt(events)));
		if(events != root_ * root_) {
		LOG(WARNING) << "Number of events is not a square, pixel cell volume cannot fully be covered in scan. "
		<< "Closest square is " << root_ * root_;
		}
		// Calculate voxel size:
		voxel_ = ROOT::Math::XYZVector(
		model->getPixelSize().x() / root_, model->getPixelSize().y() / root_, model->getSensorSize().z());
		} else {
		} else if(no_of_coordinates == 3) {
		root_ = static_cast<unsigned int>(std::lround(std::cbrt(events)));
		if(events != root_ * root_ * root_) {
		LOG(WARNING) << "Number of events is not a cube, pixel cell volume cannot fully be covered in scan. "
		<< "Closest cube is " << root_ * root_ * root_;
		}
		// Calculate voxel size:
		voxel_ = ROOT::Math::XYZVector(
		model->getPixelSize().x() / root_, model->getPixelSize().y() / root_, model->getSensorSize().z() / root_);
		}
		// Calculate voxel size:
		voxel_ = ROOT::Math::XYZVector(detector_model_->getPixelSize().x() / (scan_x_ ? root_ : 1.0),
		detector_model_->getPixelSize().y() / (scan_y_ ? root_ : 1.0),
		detector_model_->getSensorSize().z() / (scan_z_ ? root_ : 1.0));
		LOG(INFO) << "Voxel size for scan of pixel volume: " << Units::display(voxel_, {"um", "mm"});
		}

		if(output_plots_) {
		auto bins_x = static_cast<int>(output_plots_bins_per_um_ * Units::convert(model->getPixelSize().x(), "um"));
		auto bins_y = static_cast<int>(output_plots_bins_per_um_ * Units::convert(model->getPixelSize().y(), "um"));
		auto bins_z = static_cast<int>(output_plots_bins_per_um_ * Units::convert(model->getSensorSize().z(), "um"));
		auto bins_x =
		static_cast<int>(output_plots_bins_per_um_ * Units::convert(detector_model_->getPixelSize().x(), "um"));
		auto bins_y =
		static_cast<int>(output_plots_bins_per_um_ * Units::convert(detector_model_->getPixelSize().y(), "um"));
		auto bins_z =
		static_cast<int>(output_plots_bins_per_um_ * Units::convert(detector_model_->getSensorSize().z(), "um"));
		deposition_position_xy =
		CreateHistogram<TH2D>("deposition_position_xy",
		"Deposition position, x-y plane;x [#mum];y [#mum]",
		bins_x,
		-static_cast<double>(Units::convert(model->getPixelSize().x() / 2, "um")),
		static_cast<double>(Units::convert(model->getPixelSize().x() / 2, "um")),
		-static_cast<double>(Units::convert(detector_model_->getPixelSize().x() / 2, "um")),
		static_cast<double>(Units::convert(detector_model_->getPixelSize().x() / 2, "um")),
		bins_y,
		-static_cast<double>(Units::convert(model->getPixelSize().y() / 2, "um")),
		static_cast<double>(Units::convert(model->getPixelSize().y() / 2, "um")));
		-static_cast<double>(Units::convert(detector_model_->getPixelSize().y() / 2, "um")),
		static_cast<double>(Units::convert(detector_model_->getPixelSize().y() / 2, "um")));
		deposition_position_xz =
		CreateHistogram<TH2D>("deposition_position_xz",
		"Deposition position, x-z plane;x [#mum];z [#mum]",
		bins_x,
		-static_cast<double>(Units::convert(model->getPixelSize().x() / 2, "um")),
		static_cast<double>(Units::convert(model->getPixelSize().x() / 2, "um")),
		-static_cast<double>(Units::convert(detector_model_->getPixelSize().x() / 2, "um")),
		static_cast<double>(Units::convert(detector_model_->getPixelSize().x() / 2, "um")),
		bins_z,
		-static_cast<double>(Units::convert(model->getSensorSize().z() / 2, "um")),
		static_cast<double>(Units::convert(model->getSensorSize().z() / 2, "um")));
		-static_cast<double>(Units::convert(detector_model_->getSensorSize().z() / 2, "um")),
		static_cast<double>(Units::convert(detector_model_->getSensorSize().z() / 2, "um")));
		deposition_position_yz =
		CreateHistogram<TH2D>("deposition_position_yz",
		"Deposition position, y-z plane;y [#mum];z [#mum]",
		bins_y,
		-static_cast<double>(Units::convert(model->getPixelSize().y() / 2, "um")),
		static_cast<double>(Units::convert(model->getPixelSize().y() / 2, "um")),
		-static_cast<double>(Units::convert(detector_model_->getPixelSize().y() / 2, "um")),
		static_cast<double>(Units::convert(detector_model_->getPixelSize().y() / 2, "um")),
		bins_z,
		-static_cast<double>(Units::convert(model->getSensorSize().z() / 2, "um")),
		static_cast<double>(Units::convert(model->getSensorSize().z() / 2, "um")));
		-static_cast<double>(Units::convert(detector_model_->getSensorSize().z() / 2, "um")),
		static_cast<double>(Units::convert(detector_model_->getSensorSize().z() / 2, "um")));
		}
		}

		void DepositionPointChargeModule::run(Event* event) {

		ROOT::Math::XYZPoint position;
		auto model = detector_->getModel();

		if(model_ == DepositionModel::FIXED) {
		// Fixed position as read from the configuration:
		@@ -204,9 +215,10 @@ void DepositionPointChargeModule::run(Event* event) {
		} else if(model_ == DepositionModel::SCAN) {
		// Center the volume to be scanned in the center of the sensor,
		// reference point is lower left corner of one pixel volume
		auto ref = position_ + model->getMatrixSize() / 2.0 + voxel_ / 2.0 -
		ROOT::Math::XYZVector(
		model->getPixelSize().x() / 2.0, model->getPixelSize().y() / 2.0, model->getSensorSize().z() / 2.0);
		auto ref = position_ + detector_model_->getMatrixSize() / 2.0 + voxel_ / 2.0 -
		ROOT::Math::XYZVector(detector_model_->getPixelSize().x() / 2.0,
		detector_model_->getPixelSize().y() / 2.0,
		detector_model_->getSensorSize().z() / 2.0);
		LOG(DEBUG) << "Reference: " << Units::display(ref, {"um", "mm"});
		if(scan_coordinates_.size() == 3) {
		position =
		@@ -251,11 +263,9 @@ void DepositionPointChargeModule::run(Event* event) {
		DepositLine(event, position);
		} else {
		DepositPoint(event, position);
		}

		if(output_plots_) {
		auto [xpixel, ypixel] = model->getPixelIndex(position);
		auto inPixelPos = position - model->getPixelCenter(xpixel, ypixel);
		auto [xpixel, ypixel] = detector_model_->getPixelIndex(position);
		auto inPixelPos = position - detector_model_->getPixelCenter(xpixel, ypixel);
		auto in_pixel_um_x = static_cast<double>(Units::convert(inPixelPos.x(), "um"));
		auto in_pixel_um_y = static_cast<double>(Units::convert(inPixelPos.y(), "um"));
		auto in_pixel_um_z = static_cast<double>(Units::convert(position.z(), "um"));
		@@ -264,6 +274,7 @@ void DepositionPointChargeModule::run(Event* event) {
		deposition_position_yz->Fill(in_pixel_um_y, in_pixel_um_z);
		}
		}
		}

		void DepositionPointChargeModule::finalize() {
		if(output_plots_) {
		@@ -284,7 +295,7 @@ void DepositionPointChargeModule::DepositPoint(Event* event, const ROOT::Math::X

		LOG(DEBUG) << "Position (local coordinates): " << Units::display(position, {"um", "mm"});
		// Cross-check calculated position to be within sensor:
		if(!detector_->getModel()->isWithinSensor(position)) {
		if(!detector_model_->isWithinSensor(position)) {
		LOG(DEBUG) << "Requested position is outside active sensor volume.";
		return;
		}
		@@ -312,23 +323,22 @@ void DepositionPointChargeModule::DepositPoint(Event* event, const ROOT::Math::X
		}

		void DepositionPointChargeModule::DepositLine(Event* event, const ROOT::Math::XYZPoint& position) {
		auto model = detector_->getModel();

		// Vector of deposited charges and their "MCParticle"
		std::vector<DepositedCharge> charges;
		std::vector<MCParticle> mcparticles;

		// Cross-check calculated position to be within sensor:
		if(!detector_->getModel()->isWithinSensor(ROOT::Math::XYZPoint(position.x(), position.y(), 0))) {
		if(!detector_model_->isWithinSensor(position)) {
		LOG(DEBUG) << "Requested position is outside active sensor volume.";
		return;
		}

		// Start and end position of MCParticle:
		auto start_local = ROOT::Math::XYZPoint(position.x(), position.y(), -model->getSensorSize().z() / 2.0);
		auto end_local = ROOT::Math::XYZPoint(position.x(), position.y(), model->getSensorSize().z() / 2.0);
		auto start_global = detector_->getGlobalPosition(start_local);
		auto end_global = detector_->getGlobalPosition(end_local);
		// End point? It's the intersection along the line to the box. Start position is also that, but int he other direction.
		// The given position is a point the line intersects. Need to extrapolate to surfaces using this
		auto [start_local, end_local] = SensorIntersection(position);
		ROOT::Math::XYZPoint start_global = detector_->getGlobalPosition(start_local);
		ROOT::Math::XYZPoint end_global = detector_->getGlobalPosition(end_local);

		// Total number of carriers will be:
		auto charge = static_cast<unsigned int>(carriers_ * (end_local.z() - start_local.z()) / step_size_z_);
		@@ -350,6 +360,17 @@ void DepositionPointChargeModule::DepositLine(Event* event, const ROOT::Math::XY
		charges.emplace_back(position_local, position_global, CarrierType::HOLE, carriers_, 0., 0., &(mcparticles.back()));
		LOG(TRACE) << "Deposited " << carriers_ << " charge carriers of both types at global position "
		<< Units::display(position_global, {"um", "mm"}) << " in detector " << detector_->getName();

		if(output_plots_) {
		auto [xpixel, ypixel] = detector_model_->getPixelIndex(position_local);
		auto inPixelPos = position_local - detector_model_->getPixelCenter(xpixel, ypixel);
		auto in_pixel_um_x = static_cast<double>(Units::convert(inPixelPos.x(), "um"));
		auto in_pixel_um_y = static_cast<double>(Units::convert(inPixelPos.y(), "um"));
		auto in_pixel_um_z = static_cast<double>(Units::convert(position.z(), "um"));
		deposition_position_xy->Fill(in_pixel_um_x, in_pixel_um_y);
		deposition_position_xz->Fill(in_pixel_um_x, in_pixel_um_z);
		deposition_position_yz->Fill(in_pixel_um_y, in_pixel_um_z);
		}
		}

		// Dispatch the messages to the framework
		@@ -359,3 +380,34 @@ void DepositionPointChargeModule::DepositLine(Event* event, const ROOT::Math::XY
		auto deposit_message = std::make_shared<DepositedChargeMessage>(std::move(charges), detector_);
		messenger_->dispatchMessage(this, std::move(deposit_message), event);
		}

		std::tuple<ROOT::Math::XYZPoint, ROOT::Math::XYZPoint>
		DepositionPointChargeModule::SensorIntersection(const ROOT::Math::XYZPoint& line_origin) {

		double sensorLowerEdges[3] = {-detector_model_->getSensorSize().x() / 2.0,
		-detector_model_->getSensorSize().y() / 2.0,
		-detector_model_->getSensorSize().z() / 2.0};
		double sensorUpperEdges[3] = {detector_model_->getSensorSize().x() / 2.0,
		detector_model_->getSensorSize().y() / 2.0,
		detector_model_->getSensorSize().z() / 2.0};
		double lineOriginPoint[3] = {line_origin.x(), line_origin.y(), line_origin.z()};
		double lineDirection[3] = {mip_direction_.x(), mip_direction_.y(), mip_direction_.z()};

		double tMin = -INFINITY;
		double tMax = INFINITY;
		// Loop over the three coordinates
		for(int i = 0; i < 3; i++) {
		double t1 = (sensorLowerEdges[i] - lineOriginPoint[i]) / lineDirection[i];
		double t2 = (sensorUpperEdges[i] - lineOriginPoint[i]) / lineDirection[i];

		tMin = std::max(tMin, std::min(t1, t2));
		tMax = std::min(tMax, std::max(t1, t2));
		}

		ROOT::Math::XYZPoint lowerIntersectPos = line_origin + tMin * mip_direction_;
		ROOT::Math::XYZPoint upperIntersectPos = line_origin + tMax * mip_direction_;
		LOG(WARNING) << "lowerIntersectPos: " << lowerIntersectPos << ", upperIntersectPos: " << upperIntersectPos;
		LOG(WARNING) << "MIP direction: " << mip_direction_;

		return {lowerIntersectPos, upperIntersectPos};
		}

src/modules/DepositionPointCharge/DepositionPointChargeModule.hpp

+8 −0

Original line number	Diff line number	Diff line
		@@ -79,9 +79,16 @@ namespace allpix {
		*/
		void DepositLine(Event*, const ROOT::Math::XYZPoint& position);

		/**
		* @brief Finds and returns the points where a line with mip_direction through a given point intersects the sensor
		* @param line_origin Point the line goes through
		*/
		std::tuple<ROOT::Math::XYZPoint, ROOT::Math::XYZPoint> SensorIntersection(const ROOT::Math::XYZPoint& line_origin);

		Messenger* messenger_;

		std::shared_ptr<Detector> detector_;
		std::shared_ptr<allpix::DetectorModel> detector_model_;

		DepositionModel model_;
		SourceType type_;
		@@ -90,6 +97,7 @@ namespace allpix {
		double step_size_z_{};
		unsigned int root_{}, carriers_{};
		ROOT::Math::XYZVector position_{};
		ROOT::Math::XYZVector mip_direction_{};
		std::vector<std::string> scan_coordinates_{};

		bool scan_x_;

src/modules/DepositionPointCharge/README.md

+1 −0

Original line number	Diff line number	Diff line
		@@ -34,6 +34,7 @@ All charge carriers are deposited at time zero, i.e. at the beginning of the eve
		* `position`: Position in local coordinates of the sensor, where charge carriers should be deposited. Expects three values for local-x, local-y and local-z position in the sensor volume and defaults to `0um 0um 0um`, i.e. the center of first (lower left) pixel. When using source type `mip`, providing a 2D position is sufficient since it only uses the x and y coordinates. If used in scan mode, it allows you to shift the origin of each deposited charge by adding this value. If the scan is only performed in one or two dimensions, the remaining coordinate will constantly have the value given by `position`.
		* `spot_size`: Width of the Gaussian distribution used to smear the position in the `spot` model. Only one value is taken and used for all three dimensions.
		* `scan_coordinates`: Coordinates to scan over, a combination of x, y, z. Defaults to `x y z`, i.e. all three spatial coordinates. The `position` parameter is used to determine the value of the coordinates that are not scanned over if a partial scan is requested, and the start offset of the scan for the other coordinates.
		* `MIP_direction`: Unit vector giving the direction of the line along which deposits are made when the `mip` source type is used. Defaults to `0 0 1`, i.e. along the z-axis. The `position` keyword gives a point that the line will cross through with this direction.

		### Plotting parameters