Merge branch 'tformula' into 'master' (f48741fb) · Commits · XRFLab109 / Allpix Squared

src/modules/CSADigitizer/CSADigitizerModule.cpp

+15 −8

Original line number	Diff line number	Diff line
		@@ -13,7 +13,6 @@
		#include "core/utils/unit.h"
		#include "tools/ROOT.h"

		#include <TF1.h>
		#include <TFile.h>
		#include <TGraph.h>
		#include <TH1D.h>
		@@ -84,8 +83,8 @@ CSADigitizerModule::CSADigitizerModule(Configuration& config, Messenger* messeng
		auto capacitance_feedback = config_.get<double>("feedback_capacitance");
		auto resistance_feedback = tauF / capacitance_feedback;

		calculate_impulse_response_ = std::make_unique<TF1>(
		"response_function", "[0]*(TMath::Exp(-x/[1])-TMath::Exp(-x/[2]))/([1]-[2])", 0., integration_time_);
		calculate_impulse_response_ =
		std::make_unique<TFormula>("response_function", "[0]*(TMath::Exp(-x/[1])-TMath::Exp(-x/[2]))/([1]-[2])");
		calculate_impulse_response_->SetParameters(resistance_feedback, tauF, tauR);

		LOG(DEBUG) << "Parameters: cf = " << Units::display(capacitance_feedback, {"C/V", "fC/mV"})
		@@ -114,8 +113,8 @@ CSADigitizerModule::CSADigitizerModule(Configuration& config, Messenger* messeng
		auto tauF = resistance_feedback * capacitance_feedback;
		auto tauR = (capacitance_detector * capacitance_output) / (gm * capacitance_feedback);

		calculate_impulse_response_ = std::make_unique<TF1>(
		"response_function", "[0]*(TMath::Exp(-x/[1])-TMath::Exp(-x/[2]))/([1]-[2])", 0., integration_time_);
		calculate_impulse_response_ =
		std::make_unique<TFormula>("response_function", "[0]*(TMath::Exp(-x/[1])-TMath::Exp(-x/[2]))/([1]-[2])");
		calculate_impulse_response_->SetParameters(resistance_feedback, tauF, tauR);

		LOG(DEBUG) << "Parameters: rf = " << Units::display(resistance_feedback, "V*s/C")
		@@ -127,18 +126,26 @@ CSADigitizerModule::CSADigitizerModule(Configuration& config, Messenger* messeng
		<< ", tauR = " << Units::display(tauR, {"ns", "us", "ms", "s"})
		<< ", temperature = " << Units::display(config_.get<double>("temperature"), "K");
		} else if(model_ == DigitizerType::CUSTOM) {
		calculate_impulse_response_ = std::make_unique<TF1>(
		"response_function", (config_.get<std::string>("response_function")).c_str(), 0., integration_time_);
		calculate_impulse_response_ =
		std::make_unique<TFormula>("response_function", (config_.get<std::string>("response_function")).c_str());

		if(!calculate_impulse_response_->IsValid()) {
		throw InvalidValueError(
		config_, "response_function", "The response function is not a valid ROOT::TFormula expression.");
		}

		if(calculate_impulse_response_->GetNdim() != 1) {
		throw InvalidValueError(config_,
		"response_function",
		"The response function has " +
		allpix::to_string(calculate_impulse_response_->GetNdim()) +
		" dimensions, only one expected.");
		}

		auto parameters = config_.getArray<double>("response_parameters");

		// check if number of parameters match up
		if(static_cast<size_t>(calculate_impulse_response_->GetNumberFreeParameters()) != parameters.size()) {
		if(static_cast<size_t>(calculate_impulse_response_->GetNpar()) != parameters.size()) {
		throw InvalidValueError(
		config_,
		"response_parameters",

src/modules/CSADigitizer/CSADigitizerModule.hpp

+2 −2

Original line number	Diff line number	Diff line
		@@ -20,7 +20,7 @@

		#include "objects/PixelCharge.hpp"

		#include <TF1.h>
		#include <TFormula.h>
		#include <TH1D.h>
		#include <TH2D.h>

		@@ -77,7 +77,7 @@ namespace allpix {
		DigitizerType model_;

		// Function to calculate impulse response
		std::unique_ptr<TF1> calculate_impulse_response_;
		std::unique_ptr<TFormula> calculate_impulse_response_;

		// Parameters of the electronics: Noise, time-over-threshold logic
		double sigmaNoise_{}, clockToT_{}, clockToA_{}, threshold_{};

src/modules/ElectricFieldReader/ElectricFieldReaderModule.cpp

+14 −52

Original line number	Diff line number	Diff line
		@@ -18,7 +18,7 @@
		#include <utility>

		#include <Math/Vector3D.h>
		#include <TF3.h>
		#include <TFormula.h>
		#include <TH2F.h>

		#include "core/config/exceptions.h"
		@@ -119,8 +119,7 @@ void ElectricFieldReaderModule::initialize() {
		get_parabolic_field_function(thickness_domain), thickness_domain, FieldType::CUSTOM);
		} else if(field_model == ElectricField::CUSTOM) {
		LOG(TRACE) << "Adding custom electric field";
		detector_->setElectricFieldFunction(
		get_custom_field_function(thickness_domain), thickness_domain, FieldType::CUSTOM);
		detector_->setElectricFieldFunction(get_custom_field_function(), thickness_domain, FieldType::CUSTOM);
		}

		// Produce histograms if needed
		@@ -184,31 +183,18 @@ ElectricFieldReaderModule::get_parabolic_field_function(std::pair<double, double
		};
		}

		FieldFunction<ROOT::Math::XYZVector>
		ElectricFieldReaderModule::get_custom_field_function(std::pair<double, double> thickness_domain) {
		FieldFunction<ROOT::Math::XYZVector> ElectricFieldReaderModule::get_custom_field_function() {

		auto field_functions = config_.getArray<std::string>("field_function");
		auto field_parameters = config_.getArray<double>("field_parameters");

		// Derive field extent from pixel boundaries:
		auto model = detector_->getModel();
		auto min = -0.5 * model->getPixelSize();
		auto max = 0.5 * model->getPixelSize();

		// 1D field, interpret as field along z-axis:
		if(field_functions.size() == 1) {
		LOG(DEBUG) << "Found definition of 1D custom field, applying to z axis";
		auto z = std::make_shared<TF3>("z",
		field_functions.front().c_str(),
		min.x(),
		max.x(),
		min.y(),
		max.y(),
		thickness_domain.first,
		thickness_domain.second);
		auto z = std::make_shared<TFormula>("ez", field_functions.front().c_str(), false);

		// Check if number of parameters match up
		if(static_cast<size_t>(z->GetNumberFreeParameters()) != field_parameters.size()) {
		if(static_cast<size_t>(z->GetNpar()) != field_parameters.size()) {
		throw InvalidValueError(
		config_,
		"field_parameters",
		@@ -226,34 +212,12 @@ ElectricFieldReaderModule::get_custom_field_function(std::pair<double, double> t
		};
		} else if(field_functions.size() == 3) {
		LOG(DEBUG) << "Found definition of 3D custom field, applying to three Cartesian axes";
		auto x = std::make_shared<TF3>("x",
		field_functions.at(0).c_str(),
		min.x(),
		max.x(),
		min.y(),
		max.y(),
		thickness_domain.first,
		thickness_domain.second);
		auto y = std::make_shared<TF3>("y",
		field_functions.at(1).c_str(),
		min.x(),
		max.x(),
		min.y(),
		max.y(),
		thickness_domain.first,
		thickness_domain.second);
		auto z = std::make_shared<TF3>("z",
		field_functions.at(2).c_str(),
		min.x(),
		max.x(),
		min.y(),
		max.y(),
		thickness_domain.first,
		thickness_domain.second);
		auto x = std::make_shared<TFormula>("ex", field_functions.at(0).c_str(), false);
		auto y = std::make_shared<TFormula>("ey", field_functions.at(1).c_str(), false);
		auto z = std::make_shared<TFormula>("ez", field_functions.at(2).c_str(), false);

		// Check if number of parameters match up
		if(static_cast<size_t>(x->GetNumberFreeParameters() + y->GetNumberFreeParameters() + z->GetNumberFreeParameters()) !=
		field_parameters.size()) {
		if(static_cast<size_t>(x->GetNpar() + y->GetNpar() + z->GetNpar()) != field_parameters.size()) {
		throw InvalidValueError(
		config_,
		"field_parameters",
		@@ -261,16 +225,14 @@ ElectricFieldReaderModule::get_custom_field_function(std::pair<double, double> t
		}

		// Apply parameters to the functions
		for(auto n = 0; n < x->GetNumberFreeParameters(); ++n) {
		for(auto n = 0; n < x->GetNpar(); ++n) {
		x->SetParameter(n, field_parameters.at(static_cast<size_t>(n)));
		}
		for(auto n = 0; n < y->GetNumberFreeParameters(); ++n) {
		y->SetParameter(n, field_parameters.at(static_cast<size_t>(n + x->GetNumberFreeParameters())));
		for(auto n = 0; n < y->GetNpar(); ++n) {
		y->SetParameter(n, field_parameters.at(static_cast<size_t>(n + x->GetNpar())));
		}
		for(auto n = 0; n < z->GetNumberFreeParameters(); ++n) {
		z->SetParameter(
		n,
		field_parameters.at(static_cast<size_t>(n + x->GetNumberFreeParameters() + y->GetNumberFreeParameters())));
		for(auto n = 0; n < z->GetNpar(); ++n) {
		z->SetParameter(n, field_parameters.at(static_cast<size_t>(n + x->GetNpar() + y->GetNpar())));
		}

		LOG(DEBUG) << "Value of custom field at pixel center: "

src/modules/ElectricFieldReader/ElectricFieldReaderModule.hpp

+1 −1

Original line number	Diff line number	Diff line
		@@ -75,7 +75,7 @@ namespace allpix {
		* @brief Create and apply a custom field from functions
		* @param thickness_domain Domain of the thickness where the field is defined
		*/
		FieldFunction<ROOT::Math::XYZVector> get_custom_field_function(std::pair<double, double> thickness_domain);
		FieldFunction<ROOT::Math::XYZVector> get_custom_field_function();

		/**
		* @brief Read field from a file in init or apf format