DecayFit.h

Go to the documentation of this file.

#ifndef TTTRLIB_DECAYFIT_H
#define TTTRLIB_DECAYFIT_H
 
#include <iostream>
#include <cmath>
#include <algorithm>
#include <string>
#include <sstream>
 
#include "i_lbfgs.h"
#include "LvArrays.h"
#include "DecayConvolution.h"
#include "DecayStatistics.h"
 
 
struct DecayFitCorrections{
 
    double gamma = 0.0;
    double g = 1.0;
    double l1 = 0.0;
    double l2 = 0.0;
    double period = 1000;
    int convolution_stop = 0;
 
    void set_gamma(double v){
        if (v < 0.)
            gamma = 0.;        // 0 < gamma < 0.999
        else if (v > 0.999)
            gamma = 0.999;
        else
            gamma = v;
    }
 
    std::string str(){
        auto s = std::stringstream();
        s << "-- Correction factors:\n";
        s << "-- g-factor: " << g << std::endl;
        s << "-- l1, l2: " << l1 << ", " << l2 << std::endl;
        s << "-- period: " << period << std::endl;
        s << "-- convolution_stop: " << convolution_stop << std::endl;
        return s.str();
    }
 
    explicit DecayFitCorrections(
            double gamma = 0.0,
            double g = 1.0,
            double l1 = 0.0,
            double l2 = 0.0,
            double period = 1000,
            int convolution_stop = 0
    ) {
        this->gamma = gamma;
        this->g = g;
        this->l1 = l1;
        this->l2 = l2;
        this->period = period;
        this->convolution_stop = convolution_stop;
    }
 
};
 
 
struct DecayFitSettings{
 
    int fixedrho = 0;
    int softbifl = 0;
    int p2s_twoIstar = 0;
    int firstcall = 1;
    double penalty = 0.0;
 
    std::string str(){
        auto s = std::stringstream();
        s << "DECAYFITSETTINGS: " << std::endl;
        s << "-- fixedrho: " << fixedrho << std::endl;
        s << "-- softbifl: " << softbifl << std::endl;
        s << "-- p2s_twoIstar: " << p2s_twoIstar << std::endl;
        s << "-- firstcall: " << firstcall << std::endl;
        s << "-- penalty: " << penalty << std::endl;
        return s.str();
    }
};
 
 
struct DecayFitIntegrateSignals{
 
    DecayFitCorrections* corrections = nullptr;
 
    double Sp = 0.0;
 
    double Ss = 0.0;
 
    double Bp = 0.0;
 
    double Bs = 0.0;
 
    double B = 0.0;
 
    double Bexpected = 0.0;
 
 
    double Fp(){
        double g = 1.0;
        if(corrections != nullptr){
            g = corrections->g;
        }
        if(g == 1.0){
            return (Sp - Bp);
        } else{
            return (Sp - g * Bp) / (1. - g);
        }
    }
 
    double Fs(){
        double g = 1.0;
        if(corrections != nullptr){
            g = corrections->g;
        }
        if(g == 1.0){
            return (Ss - Bs);
        } else{
            return (Ss - g * Bs) / (1. - g);
        }
    }
 
    double r(){
        double fp = Fp();
        double fs = Fs();
        double g = 1.0, l1 = 0.0, l2 = 0.0;
        if(corrections != nullptr){
            g = corrections->g;
            l1 = corrections->l1;
            l2 = corrections->l2;
        }
        double nom = (fp - g * fs);
        double denom = (fp * (1. - 3. * l2) + (2. - 3. * l1) * g * fs);
        return nom / denom;
    }
 
    double rho(double tau, double r0){
        double rh = tau / (r0 / r() - 1.);        // rho = tau/(r0/r-1)
        return std::max(rh, 1.e-4);
    }
 
    double rs(){
        double g = 1.0, l1 = 0.0, l2 = 0.0;
        if(corrections != nullptr){
            g = corrections->g;
            l1 = corrections->l1;
            l2 = corrections->l2;
        }
        return (Sp - g * Ss) / (Sp * (1. - 3. * l2) + (2. - 3. * l1) * g * Ss);
    }
 
    void compute_signal_and_background(MParam *p);
 
    void normM(double *M, int Nchannels);
 
    void normM(double *M, double s, int Nchannels);
 
 
    void normM_p2s(double *M, int Nchannels);
 
 
    std::string str(){
        auto s = std::stringstream();
        s << "-- Signals: " << std::endl;
        s << "-- Bp, Bs: " << Bp << ", " << Bs << std::endl;
        s << "-- Sp, Ss: " << Sp << ", " << Ss << std::endl;
        s << "-- Fp, Fs: " << Fp() << ", " << Fs() << std::endl;
        s << "-- r: " << r() << std::endl;
        return s.str();
    }
 
    DecayFitIntegrateSignals(DecayFitCorrections *corrections = nullptr){
        this->corrections = corrections;
    }
 
};
 
 
class DecayFit{
 
 
public:
 
    static int modelf(double *param,
                       double *irf,
                       double *bg,
                       int Nchannels,
                       double dt,
                       double *corrections,
                       double *mfunction
    ) {
        return 0;
    };
 
    static double targetf(double *x, void *pv){
        return 0.0;
    };
 
 
    static double fit(double* x, short* fixed, MParam* p){
        return 0.0;
    };
 
 
    static void correct_input(double *x, double *xm, LVDoubleArray *corrections, int return_r){};
 
};
 
 
#endif //TTTRLIB_DECAYFIT_H