Pda.h

Go to the documentation of this file.

#ifndef TTTRLIB_PDA_H
#define TTTRLIB_PDA_H
 
#include <vector>
#include <iostream>
#include <cmath>
 
#include "TTTR.h"
#include "PdaCallback.h"
 
 
class Pda {
 
private:
    bool _is_valid_sgsr = false;
    PdaCallback* _histogram_function;
 
    std::vector<double> _probability_ch1;
 
    std::vector<double> _amplitudes;
 
 
protected:
 
    unsigned int _n_2d_max = 300;
 
    unsigned int _n_2d_min = 3;
 
    double _bg_ch1 = 0.0;
 
    double _bg_ch2 = 0.0;
 
    std::vector<double> _S1S2;
 
    std::vector<double> pF;
 
 
public:
 
    void evaluate();
 
    Pda(
        int hist2d_nmax=300,
        int hist2d_nmin=5,
        double background_ch1=0.0,
        double background_ch2=0.0,
        std::vector<double> pF = std::vector<double>()
    ){
        set_max_number_of_photons(std::abs(hist2d_nmax));
        _n_2d_min = std::abs(hist2d_nmin);
        _histogram_function = new PdaCallback();
        _bg_ch1 = background_ch1;
        _bg_ch2 = background_ch2;
        setPF(pF.data(), pF.size());
    }
 
    ~Pda() = default;
    
    void append(double amplitude, double probability_ch1) {
        _is_valid_sgsr = false;
        _amplitudes.push_back(amplitude);
        _probability_ch1.push_back(probability_ch1);
    }
 
    void clear_probability_ch1() {
        _is_valid_sgsr = false;
        _amplitudes.clear();
        _probability_ch1.clear();
    }
 
    void get_amplitudes(double **output_view, int *n_output) {
        *n_output = _amplitudes.size();
        *output_view = _amplitudes.data();
    }
 
    void set_amplitudes(double *input, int n_input) {
        _amplitudes.clear();
        _is_valid_sgsr = false;
        for(int i = 0; i < n_input; i++)
            _amplitudes.emplace_back(input[i]);
    }
 
    void set_callback(PdaCallback* cb){
        _histogram_function = cb;
    }
 
public:
 
    void get_S1S2_matrix(double **output, int *n_output1, int *n_output2){
        if(!_is_valid_sgsr) evaluate();
        auto* t = (double *) malloc(_S1S2.size() * sizeof(double));
        for(unsigned int i = 0; i < _S1S2.size(); i++) 
            t[i] = _S1S2[i];
        *output = t;
        *n_output1 = int (_n_2d_max + 1);
        *n_output2 = int (_n_2d_max + 1);
    }
 
    void set_probability_spectrum_ch1(double *input, int n_input){
        _is_valid_sgsr = false;
        _amplitudes.clear();
        _probability_ch1.clear();
        int n_components = n_input / 2;
        for(int i=0; i < n_components; i++){
            double amplitude = input[2 * i];
            double probability_green = input[(2 * i) + 1];
            _amplitudes.emplace_back(amplitude);
            _probability_ch1.emplace_back(probability_green);
        }
    }
 
    void get_probabilities_ch1(double **output_view, int *n_output) {
        *n_output = _probability_ch1.size();
        *output_view = _probability_ch1.data();
    }
 
    void set_probabilities_ch1(double *input, int n_input) {
        _probability_ch1.clear();
        _is_valid_sgsr = false;
        for(int i = 0; i < n_input; i++)
            _probability_ch1.emplace_back(input[i]);
    }
 
 
    void get_probability_spectrum_ch1(double** output, int* n_output) {
        int n = (int)_amplitudes.size() * 2;
        auto temp = (double*) malloc(sizeof(double) * n);
        for(unsigned int i=0; i < _amplitudes.size(); i++){
            temp[2 * i] = _amplitudes[i];
            temp[2 * i + 1] = _probability_ch1[i];
        }
        *n_output = n;
        *output = temp;
    }
 
    void setPF(double *input, int n_input){
#ifdef VERBOSE_TTTRLIB
        std::clog << "-- Setting pF " << std::endl;
#endif
        _is_valid_sgsr = false;
        pF.assign(input, input + n_input);
    }
 
    void getPF(double** output_view, int* n_output){
        *n_output = pF.size();
        *output_view = pF.data();
    }
 
    void get_1dhistogram(
            double **histogram_x, int *n_histogram_x,
            double **histogram_y, int *n_histogram_y,
            double x_max=1000.0, double x_min=0.01, int n_bins=81,
            bool log_x=true,
            std::vector<double> s1s2 = std::vector<double>(),
            int n_min=-1,
            bool skip_zero_photon=true,
            std::vector<double> amplitudes = std::vector<double>(),
            std::vector<double> probabilities_ch1 = std::vector<double>()
    );
 
    static void compute_experimental_histograms(
        TTTR* tttr_data,
        double** s1s2, int* dim1, int* dim2,
        double** ps, int* dim_ps,
        int** tttr_indices, int* n_tttr_indices,
        std::vector<int> channels_1,
        std::vector<int> channels_2,
        int maximum_number_of_photons,
        int minimum_number_of_photons,
        double minimum_time_window_length
    );
 
    static void S1S2_pF(
        std::vector<double> &S1S2,
        std::vector<double> &pF,
        unsigned int Nmax,
        double background_ch1,
        double background_ch2,
        std::vector<double> &p_ch1,
        std::vector<double> &amplitudes
    );
 
    static void conv_pF(
        std::vector<double> &S1S2,
        const std::vector<double> &F1F2,
        unsigned int Nmax,
        double background_ch1,
        double background_ch2
    );
 
    static void poisson_0toN(
        std::vector<double> &return_p,
        int start_idx,
        double lam,
        int return_dim
    );
 
 
    unsigned int get_max_number_of_photons() const{
        return _n_2d_max;
    }
 
    void set_max_number_of_photons(unsigned int nmax){
        _n_2d_max = nmax;
        _S1S2.resize((_n_2d_max + 1) * (_n_2d_max + 1));
        _is_valid_sgsr = false;
    }
 
    unsigned int get_min_number_of_photons() const{
        return this->_n_2d_min;
    }
 
    void set_min_number_of_photons(unsigned int nmin){
        this->_n_2d_min = nmin;
        _is_valid_sgsr = false;
    }
 
    double get_ch1_background() const{
        return _bg_ch1;
    }
 
    void set_ch1_background(double bg){
        _is_valid_sgsr = false;
        _bg_ch1 = bg;
    }
 
    double get_ch2_background() const{
        return _bg_ch2;
    }
 
    void set_ch2_background(double br) {
        _bg_ch2 = br;
        _is_valid_sgsr = false;
    }
 
    bool is_valid_sgsr() const{
        return _is_valid_sgsr;
    }
 
    void set_valid_sgsr(bool v){
        _is_valid_sgsr = v;
    }
 
};
 
 
#endif //TTTRLIB_PDA_H