Commit f5c0e4db authored by Paul Schütze's avatar Paul Schütze
Browse files

Merge branch 'better_pulse' into 'master' 

Improve Pulse Object

See merge request allpix-squared/allpix-squared!565
parents a88c7c0b 58434684

src/modules/CSADigitizer/CSADigitizerModule.cpp

+7 −8
Original line number Diff line number Diff line
@@ -221,7 +221,6 @@ void CSADigitizerModule::run(Event* event) { 
            throw ModuleError("No pulse information available.");

        }



        auto pulse_vec = pulse.getPulse(); // the vector of the charges

        auto timestep = pulse.getBinning();

        LOG(DEBUG) << "Timestep: " << timestep << " integration_time: " << integration_time_;

        auto ntimepoints = static_cast<size_t>(ceil(integration_time_ / timestep));
@@ -255,18 +254,18 @@ void CSADigitizerModule::run(Event* event) { 
        });



        std::vector<double> amplified_pulse_vec(ntimepoints);

        auto input_length = pulse_vec.size();

        LOG(TRACE) << "Preparing pulse for pixel " << pixel_index << ", " << pulse_vec.size() << " bins of "

        LOG(TRACE) << "Preparing pulse for pixel " << pixel_index << ", " << pulse.size() << " bins of "

                   << Units::display(timestep, {"ps", "ns"}) << ", total charge: " << Units::display(pulse.getCharge(), "e");

        // convolution of the pulse (size input_length) with the impulse response (size ntimepoints)



        // Convolution of the input pulse with the impulse response (size ntimepoints)

        for(size_t k = 0; k < ntimepoints; ++k) {

            double outsum{};

            // convolution: multiply pulse_vec.at(k - i) * impulse_response_function_.at(i), when (k - i) < input_length

            // Convolution: multiply pulse.at(k - i) * impulse_response_function_.at(i), when (k - i) < input length

            // -> no point to start i at 0, start from jmin:

            size_t jmin = (k >= input_length - 1) ? k - (input_length - 1) : 0;

            size_t jmin = (k >= pulse.size() - 1) ? k - (pulse.size() - 1) : 0;

            for(size_t i = jmin; i <= k; ++i) {

                if((k - i) < input_length) {

                    outsum += pulse_vec.at(k - i) * impulse_response_function_.at(i);

                if((k - i) < pulse.size()) {

                    outsum += pulse.at(k - i) * impulse_response_function_.at(i);

                }

            }

            amplified_pulse_vec.at(k) = outsum;

src/modules/DefaultDigitizer/DefaultDigitizerModule.cpp

+3 −4
Original line number Diff line number Diff line
@@ -329,16 +329,15 @@ double DefaultDigitizerModule::time_of_arrival(const PixelCharge& pixel_charge, 
    // If this PixelCharge has a pulse, we can find out when it crossed the threshold:

    const auto& pulse = pixel_charge.getPulse();

    if(pulse.isInitialized()) {

        auto charges = pulse.getPulse();

        std::vector<double>::iterator bin;

        auto bin = pulse.begin();

        double integrated_charge = 0;

        for(bin = charges.begin(); bin != charges.end(); bin++) {

        for(; bin != pulse.end(); bin++) {

            integrated_charge += *bin;

            if(integrated_charge >= threshold) {

                break;

            }

        }

        return pulse.getBinning() * static_cast<double>(std::distance(charges.begin(), bin));

        return pulse.getBinning() * static_cast<double>(std::distance(pulse.begin(), bin));

    } else {

        LOG_ONCE(INFO) << "Simulation chain does not allow for time-of-arrival calculation";

        return 0;

src/modules/PulseTransfer/PulseTransferModule.cpp

+15 −11
Original line number Diff line number Diff line
@@ -145,7 +145,13 @@ void PulseTransferModule::run(Event* event) { 


            // Generate pseudo-pulse:

            Pulse pulse(timestep_);

            try {

                pulse.addCharge(propagated_charge.getCharge(), propagated_charge.getLocalTime());

            } catch(const PulseBadAllocException& e) {

                LOG(ERROR) << e.what() << std::endl

                           << "Ignoring pulse contribution at time "

                           << Units::display(propagated_charge.getLocalTime(), {"ms", "us", "ns"});

            }

            pixel_pulse_map[pixel_index] += pulse;



            auto px = pixel_charge_map[pixel_index];
@@ -183,11 +189,10 @@ void PulseTransferModule::run(Event* event) { 
            h_induced_pixel_charge_->Fill(pulse.getCharge() / 1e3);



            auto step = pulse.getBinning();

            auto pulse_vec = pulse.getPulse();

            double charge = 0;



            for(auto bin = pulse_vec.begin(); bin != pulse_vec.end(); ++bin) {

                auto time = step * static_cast<double>(std::distance(pulse_vec.begin(), bin));

            for(auto bin = pulse.begin(); bin != pulse.end(); ++bin) {

                auto time = step * static_cast<double>(std::distance(pulse.begin(), bin));

                h_induced_pulses_->Fill(time, *bin);

                p_induced_pulses_->Fill(time, *bin);
@@ -255,19 +260,18 @@ void PulseTransferModule::finalize() { 


void PulseTransferModule::create_pulsegraphs(uint64_t event_num, const Pixel::Index& index, const Pulse& pulse) const {

    auto step = pulse.getBinning();

    auto pulse_vec = pulse.getPulse();

    LOG(TRACE) << "Preparing pulse for pixel " << index << ", " << pulse_vec.size() << " bins of "

    LOG(TRACE) << "Preparing pulse for pixel " << index << ", " << pulse.size() << " bins of "

               << Units::display(step, {"ps", "ns"}) << ", total charge: " << Units::display(pulse.getCharge(), "e");



    // Generate x-axis:

    std::vector<double> time(pulse_vec.size());

    std::vector<double> time(pulse.size());

    // clang-format off

    std::generate(time.begin(), time.end(), [n = 0.0, step]() mutable { auto now = n; n += step; return now; });

    // clang-format on



    std::string name =

        "pulse_ev" + std::to_string(event_num) + "_px" + std::to_string(index.x()) + "-" + std::to_string(index.y());

    auto* pulse_graph = new TGraph(static_cast<int>(pulse_vec.size()), &time[0], &pulse_vec[0]);

    auto* pulse_graph = new TGraph(static_cast<int>(pulse.size()), &time[0], &pulse[0]);

    pulse_graph->GetXaxis()->SetTitle("t [ns]");

    pulse_graph->GetYaxis()->SetTitle("Q_{ind} [e]");

    pulse_graph->SetTitle(("Induced charge per unit step time in pixel (" + std::to_string(index.x()) + "," +
@@ -276,7 +280,7 @@ void PulseTransferModule::create_pulsegraphs(uint64_t event_num, const Pixel::In 
                              .c_str());

    getROOTDirectory()->WriteTObject(pulse_graph, name.c_str());



    std::vector<double> abs_vec = pulse_vec;

    std::vector<double> abs_vec = pulse;

    std::for_each(abs_vec.begin(), abs_vec.end(), [](auto& bin) { bin = std::fabs(bin); });

    name = "pulse_abs_ev" + std::to_string(event_num) + "_px" + std::to_string(index.x()) + "-" + std::to_string(index.y());

    auto* pulse_abs_graph = new TGraph(static_cast<int>(abs_vec.size()), &time[0], &abs_vec[0]);
@@ -289,7 +293,7 @@ void PulseTransferModule::create_pulsegraphs(uint64_t event_num, const Pixel::In 
                                  .c_str());

    getROOTDirectory()->WriteTObject(pulse_abs_graph, name.c_str());



    std::vector<double> current_vec = pulse_vec;

    std::vector<double> current_vec = pulse;

    // Convert charge bins to current in uA

    std::for_each(current_vec.begin(), current_vec.end(), [step](auto& bin) {

        bin = static_cast<double>(Units::convert(bin, "fC") / Units::convert(step, "ns"));
@@ -309,7 +313,7 @@ void PulseTransferModule::create_pulsegraphs(uint64_t event_num, const Pixel::In 
    // Generate graphs of integrated charge over time:

    std::vector<double> charge_vec;

    double charge = 0;

    for(const auto& bin : pulse_vec) {

    for(const auto& bin : pulse) {

        charge += bin;

        charge_vec.push_back(charge);

    }

src/modules/TransientPropagation/TransientPropagationModule.cpp

+8 −2
Original line number Diff line number Diff line
@@ -402,8 +402,14 @@ TransientPropagationModule::propagate(Event* event, 
                       << " q = " << Units::display(induced, "e");



            // Create pulse if it doesn't exist. Store induced charge in the returned pulse iterator

            auto pixel_map_iterator = pixel_map.emplace(pixel_index, Pulse(timestep_));

            auto pixel_map_iterator = pixel_map.emplace(pixel_index, Pulse(timestep_, integration_time_));

            try {

                pixel_map_iterator.first->second.addCharge(induced, initial_time + runge_kutta.getTime());

            } catch(const PulseBadAllocException& e) {

                LOG(ERROR) << e.what() << std::endl

                           << "Ignoring pulse contribution at time "

                           << Units::display(initial_time + runge_kutta.getTime(), {"ms", "us", "ns"});

            }



            if(output_plots_) {

                potential_difference_->Fill(std::fabs(ramo - last_ramo));

src/objects/Pulse.cpp

+23 −17
Original line number Diff line number Diff line
@@ -8,35 +8,43 @@ 
 */



#include "Pulse.hpp"

#include "objects/exceptions.h"



#include <cmath>

#include <numeric>



using namespace allpix;



Pulse::Pulse(double time_bin) : bin_(time_bin), initialized_(true) {}

Pulse::Pulse(double time_bin) noexcept : bin_(time_bin), initialized_(true) {}



Pulse::Pulse(double time_bin, double total_time) : bin_(time_bin), initialized_(true) {

    auto bins = static_cast<size_t>(std::lround(total_time / bin_));

    try {

        this->reserve(bins);

    } catch(const std::bad_alloc& e) {

        PulseBadAllocException(bins, total_time, e.what());

    }

}



void Pulse::addCharge(double charge, double time) {

    // For uninitialized pulses, store all charge in the first bin:

    auto bin = (initialized_ ? static_cast<size_t>(std::lround(time / bin_)) : 0);



    try {

        // Adapt pulse storage vector:

    if(bin >= pulse_.size()) {

        pulse_.resize(bin + 1);

        if(bin >= this->size()) {

            this->resize(bin + 1);

        }

        this->at(bin) += charge;

    } catch(const std::bad_alloc& e) {

        PulseBadAllocException(bin + 1, time, e.what());

    }

    pulse_.at(bin) += charge;

}



int Pulse::getCharge() const {

    double charge = std::accumulate(pulse_.begin(), pulse_.end(), 0.0);

    double charge = std::accumulate(this->begin(), this->end(), 0.0);

    return static_cast<int>(std::round(charge));

}



const std::vector<double>& Pulse::getPulse() const {

    return pulse_;

}



double Pulse::getBinning() const {

    return bin_;

}
@@ -46,8 +54,6 @@ bool Pulse::isInitialized() const { 
}



Pulse& Pulse::operator+=(const Pulse& rhs) {

    auto rhs_pulse = rhs.getPulse();



    // Allow to initialize uninitialized pulse

    if(!this->initialized_) {

        this->bin_ = rhs.getBinning();
@@ -60,13 +66,13 @@ Pulse& Pulse::operator+=(const Pulse& rhs) { 
    }



    // If new pulse is longer, extend:

    if(this->pulse_.size() < rhs_pulse.size()) {

        this->pulse_.resize(rhs_pulse.size());

    if(this->size() < rhs.size()) {

        this->resize(rhs.size());

    }



    // Add up the individual bins:

    for(size_t bin = 0; bin < rhs_pulse.size(); bin++) {

        this->pulse_.at(bin) += rhs_pulse.at(bin);

    for(size_t bin = 0; bin < rhs.size(); bin++) {

        this->at(bin) += rhs.at(bin);

    }



    return *this;