Commit b3c5ca4e authored by Simon Spannagel's avatar Simon Spannagel
Browse files

Pulse: directly inherit from std::vector

parent 15565f2b
Loading
Loading
Loading
Loading
+4 −5
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 "
        auto input_length = pulse.size();
        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)
        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;
            for(size_t i = jmin; i <= k; ++i) {
                if((k - i) < input_length) {
                    outsum += pulse_vec.at(k - i) * impulse_response_function_.at(i);
                    outsum += pulse.at(k - i) * impulse_response_function_.at(i);
                }
            }
            amplified_pulse_vec.at(k) = outsum;
+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; 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;
+8 −10
Original line number Diff line number Diff line
@@ -183,11 +183,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 +254,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 +274,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 +287,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 +307,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);
    }
+8 −14
Original line number Diff line number Diff line
@@ -22,21 +22,17 @@ void Pulse::addCharge(double charge, double time) {
    auto bin = (initialized_ ? static_cast<size_t>(std::lround(time / bin_)) : 0);

    // Adapt pulse storage vector:
    if(bin >= pulse_.size()) {
        pulse_.resize(bin + 1);
    if(bin >= this->size()) {
        this->resize(bin + 1);
    }
    pulse_.at(bin) += charge;
    this->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 +42,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 +54,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;
+2 −9
Original line number Diff line number Diff line
@@ -20,7 +20,7 @@ namespace allpix {
     * @brief Pulse holding induced charges as a function of time
     * @warning This object is special and is not meant to be written directly to a tree (not inheriting from \ref Object)
     */
    class Pulse {
    class Pulse : public std::vector<double> {
    public:
        /**
         * @brief Construct a new pulse
@@ -45,12 +45,6 @@ namespace allpix {
         */
        int getCharge() const;

        /**
         * @brief Function to retrieve the full pulse shape
         * @return Constant reference to the pulse vector
         */
        const std::vector<double>& getPulse() const;

        /**
         * @brief Function to retrieve time binning of pulse
         * @return Width of one pulse bin in nanoseconds
@@ -72,10 +66,9 @@ namespace allpix {
        /**
         * @brief Default constructor for ROOT I/O
         */
        ClassDef(Pulse, 2); // NOLINT
        ClassDef(Pulse, 3); // NOLINT

    private:
        std::vector<double> pulse_;
        double bin_{};
        bool initialized_{};
    };