//Renata Kopecna

#include <complex>
#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <math.h>
#include <TMath.h>

#include <TDecompChol.h>
#include <TDecompBK.h>
#include <TMatrixD.h>
#include <Math/GSLIntegrator.h>

#include <funcs.hh>
#include <event.hh>
#include <helpers.hh>
#include <constants.hh>
#include <integrals.hh>

#include <spdlog.h>

namespace fcnc
{


TMatrixD get_bsz_mu_invcov(){ //Used one in toystudy.hh
    TMatrixDSym bsz_cov(21);
    bsz_cov = get_bsz_cov();

    //need to do cholesky decomposition, use root functionality for this
    //TMatrixDSym mcoeffs_sym(21, &coeffs_sym[0]);
    TDecompChol chol(bsz_cov);
    bool success = chol.Decompose();
    assert(success);
    TMatrixD mU = chol.GetU();//gets upper triangular matrix

    if (spdlog_debug()){
        for (unsigned int i =0; i<21; i++){//column
            std::cout << i << std::endl;
            for (unsigned int j =0; j<21; j++){//row
                std::cout << mU(j,i) << " ";
            }
            std::cout << std::endl;
        }
    }

    return mU;
}

TMatrixDSym get_bsz_invcov()
{
    TMatrixDSym covariance_sym = TMatrixDSym(21);
    covariance_sym = get_bsz_cov();
    bool success = false;
    TDecompBK decompbk(covariance_sym);
    TMatrixDSym bsz_invcov = TMatrixDSym(decompbk.Invert(success));
    assert(success);
    return bsz_invcov;
}

TMatrixDSym get_bsz_cov()
{

    //current covariance matrix
    double covariance_raw[21*21] = {0.000837635, 0.00471065, 0.00272131, 0.00471065, 0.0659712, 0.204148, 0.00272131, 0.204148, 2.67119, 4.69403e-5, 0.000919267, 0.00714105, 0.0020347, 0.022787, 0.0846849, 0.0222044, 0.228333, 1.10882, 0.000602086, 0.00227154, 0.000674336, 0.00331857, 0.0250234, 0.079078, 0.00187104, 0.0687236, 0.713438, 3.79001e-5, 0.00134723, 0.00865921, 0.00239851, 0.0371657, 0.029583, 0.0296877, 0.312199, 0.0215381, 1.31212e-5, 0.000248351, -0.00751951, 0.00170386, 0.0218497, -0.0768858, 0.021364, 0.208973, -0.336619, 1.73162e-5, 0.000722872, 0.00072477, 0.00172356, 0.0211513, 0.0374082, 0.0213424, 0.190446, 0.359855, 0.000826815, 0.000674452, -0.022157, 0.00932974, 0.0304414, -0.109803, 0.0686369, 0.154917, -1.14562, 4.69403e-5, 0.0020347, 0.0222044, 0.000919267, 0.022787, 0.228333, 0.00714105, 0.0846849, 1.10882, 0.000694641, 0.00347885, 0.00415084, 0.00347885, 0.0353042, 0.144041, 0.00415084, 0.144041, 1.05126, 3.04401e-5, 3.57306e-5, 0.00355147, 0.000650453, 0.00688257, 0.0593046, 0.00522776, 0.0417512, 0.36039, 0.000809924, 0.00503414, -0.0144592, 0.00424412, 0.0423395, -0.0182202, 0.00664211, 0.138549, 0.397771, 0.000652147, 0.00340264, -0.0202476, 0.00325978, 0.0270724, -0.114933, 0.00381811, 0.0778571, -0.214537, 0.00065194, 0.00303435, -0.00428123, 0.00327701, 0.028583, 0.0341862, 0.00394892, 0.10177, 0.331577, 0.000942311, 0.00179228, -0.0210447, 0.0082808, 0.0190827, -0.18232, 0.0359498, 0.0574551, -0.914838, 0.000602086, 0.00331857, 0.00187104, 0.00227154, 0.0250234, 0.0687236, 0.000674336, 0.079078, 0.713438, 3.04401e-5, 0.000650453, 0.00522776, 3.57306e-5, 0.00688257, 0.0417512, 0.00355147, 0.0593046, 0.36039, 0.000432775, 0.00165482, 0.000605487, 0.00165482, 0.0165823, 0.059165, 0.000605487, 0.059165, 0.430695, 2.32114e-5, 0.000937822, 0.00645455, -6.04184e-5, 0.0117957, 0.0462642, 0.0042958, 0.0833038, 0.140056, 6.01978e-6, 0.000154238, -0.0052701, -8.61667e-5, 0.00572307, -0.0128712, 0.00312849, 0.0501177, -0.083906, 9.00023e-6, 0.000503353, 0.000573722, -7.58925e-5, 0.00611674, 0.0257193, 0.00313087, 0.0460709, 0.133667, 0.000590954, 0.000464079, -0.0159601, 0.00315256, 0.0153668, -0.0107164, 0.0167505, 0.083256, -0.0400335, 3.79001e-5, 0.00239851, 0.0296877, 0.00134723, 0.0371657, 0.312199, 0.00865921, 0.029583, 0.0215381, 0.000809924, 0.00424412, 0.00664211, 0.00503414, 0.0423395, 0.138549, -0.0144592, -0.0182202, 0.397771, 2.32114e-5, -6.04184e-5, 0.0042958, 0.000937822, 0.0117957, 0.0833038, 0.00645455, 0.0462642, 0.140056, 0.00110852, 0.00658825, -0.0202371, 0.00658825, 0.0682608, 0.0138459, -0.0202371, 0.0138459, 2.34402, 0.000825342, 0.004522, -0.0218767, 0.00509604, 0.0448809, -0.101388, -0.0149136, -0.0227422, 0.849315, 0.000823664, 0.00388965, -0.00393922, 0.00509672, 0.0397888, 0.0548511, -0.0148331, -0.0137692, 0.473698, 0.00120989, 0.00205499, -0.0264334, 0.0116269, 0.0278953, -0.204925, 0.00157102, 2.28401e-5, 0.19199, 1.31212e-5, 0.00170386, 0.021364, 0.000248351, 0.0218497, 0.208973, -0.00751951, -0.0768858, -0.336619, 0.000652147, 0.00325978, 0.00381811, 0.00340264, 0.0270724, 0.0778571, -0.0202476, -0.114933, -0.214537, 6.01978e-6, -8.61667e-5, 0.00312849, 0.000154238, 0.00572307, 0.0501177, -0.0052701, -0.0128712, -0.083906, 0.000825342, 0.00509604, -0.0149136, 0.004522, 0.0448809, -0.0227422, -0.0218767, -0.101388, 0.849315, 0.000756032, 0.00372113, -0.0213639, 0.00372113, 0.0359388, 0.00028928, -0.0213639, 0.00028928, 2.68844, 0.000756032, 0.00313686, -0.00507995, 0.0036848, 0.0277685, 0.0556352, -0.0217047, -0.0652209, 0.825691, 0.00089759, 0.00144256, -0.0188003, 0.0072161, 0.0179819, -0.0877721, -0.0281293, -0.0385387, 1.61537, 1.73162e-5, 0.00172356, 0.0213424, 0.000722872, 0.0211513, 0.190446, 0.00072477, 0.0374082, 0.359855, 0.00065194, 0.00327701, 0.00394892, 0.00303435, 0.028583, 0.10177, -0.00428123, 0.0341862, 0.331577, 9.00023e-6, -7.58925e-5, 0.00313087, 0.000503353, 0.00611674, 0.0460709, 0.000573722, 0.0257193, 0.133667, 0.000823664, 0.00509672, -0.0148331, 0.00388965, 0.0397888, -0.0137692, -0.00393922, 0.0548511, 0.473698, 0.000756032, 0.0036848, -0.0217047, 0.00313686, 0.0277685, -0.0652209, -0.00507995, 0.0556352, 0.825691, 0.000756032, 0.00316639, -0.00503577, 0.00316639, 0.0276287, 0.0511444, -0.00503577, 0.0511444, 0.646067, 0.000900667, 0.00145139, -0.0189859, 0.00693449, 0.0173072, -0.119675, 0.00493465, 0.0386434, 0.478871, 0.000826815, 0.00932974, 0.0686369, 0.000674452, 0.0304414, 0.154917, -0.022157, -0.109803, -1.14562, 0.000942311, 0.0082808, 0.0359498, 0.00179228, 0.0190827, 0.0574551, -0.0210447, -0.18232, -0.914838, 0.000590954, 0.00315256, 0.0167505, 0.000464079, 0.0153668, 0.083256, -0.0159601, -0.0107164, -0.0400335, 0.00120989, 0.0116269, 0.00157102, 0.00205499, 0.0278953, 2.28401e-5, -0.0264334, -0.204925, 0.19199, 0.00089759, 0.0072161, -0.0281293, 0.00144256, 0.0179819, -0.0385387, -0.0188003, -0.0877721, 1.61537, 0.000900667, 0.00693449, 0.00493465, 0.00145139, 0.0173072, 0.0386434, -0.0189859, -0.119675, 0.478871, 0.0040111, 0.00460369, -0.0993374, 0.00460369, 0.0493411, 0.140689, -0.0993374, 0.140689, 4.85487 };
    
    //spdlog::info("chi2 start"std::endl;

    double covariance_reordered[21*21];
    //double covariance_reordered_old[21*21];
    //want to reorder covariances into proper matrix A0a0 A0a1 A0a2 A1a0 A1a1 A1a2 A12a0 and so on
    for (unsigned int i =0; i<21; i++)//column
        for (unsigned int j =0; j<21; j++)//row
        {
            int to = 21*j+i;
            int parcol = (to%21);
            int parrow = to / 21;
            int ffcol = parcol/3;
            int ffrow = parrow/3;
            int alphacol = parcol%3;
            int alpharow = parrow%3;
            int from = 9*(7*ffrow+ffcol) + 3*alpharow + alphacol;
            covariance_reordered[to] = covariance_raw[from];
            //covariance_reordered_old[to] = covariance_old[from];
        }

    //for (unsigned int i=0; i<21; i++)
    //  spdlog::info("COVCOMP new ",sqrt(covariance_reordered[i*21+i])" old ",sqrt(covariance_reordered_old[i*21+i]) << );
    //now need to remove 100% correlated obs
    //a12_0 = (mB^2-mK*^2)/(8mBmK*)A3(0) = (mB^2-mK*^2)/(8mBmK*)A0(0)  B6 in BSZ paper
    //t2_0 = t1_0
    double covariance_final[21*21];
    for (unsigned int i =0; i<21; i++)//column
        for (unsigned int j =0; j<21; j++)//row
        {
            covariance_final[j*21+i] = covariance_reordered[j*21+i];
            if ((i==6 || j==6) || (i==15 || j==15))
                covariance_final[j*21+i] = 0.0;
            if ((i==6 && j==6) || (i==15 && j==15))
                covariance_final[j*21+i] = 0.1;
        }
    //test output
    //for (unsigned int i =0; i<21; i++)//column
    //  {
    //    std::cout << i << std::endl;
    //    for (unsigned int j =0; j<21; j++)//row
    //	std::cout << covariance_final[21*j+i] << " ";
    //    std::cout << std::endl << std::endl;
    //  }//seems to be ok
    // spdlog::info("BSZ correlations",std::endl << "    ";
    // for (unsigned int i =0; i<21; i++)
    //   std::cout << std::fixed << std::setw(5) << i << " ";
    // std::cout << std::endl;
    // for (unsigned int i =0; i<21; i++)//column
    //   {
    //     std::cout << i << (i<10 ? "   " : "  ");
    //     for (unsigned int j =0; j<21; j++)//row
    // 	  std::cout << std::fixed << std::setprecision(2) << std::setw(5) << covariance_final[21*j+i]/sqrt(covariance_final[21*i+i]*covariance_final[21*j+j]) << " ";
    //     std::cout << std::endl;
    //  }//seems to be ok


    //we should be able to invert this now
    TMatrixDSym covariance_sym(21, &covariance_final[0]);
    return covariance_sym;
}

double crystalball(double m, double mean, double sigma, double alpha, double n){

    //implementation taken from RooFit: https://root.cern.ch/doc/master/RooCBShape_8cxx_source.html#l00056

    double t = (m-mean)/sigma;
    if (alpha < 0) t = -t;
    
    double absAlpha = fabs(alpha);
    
    if (t >= -absAlpha) {
        return exp(-0.5*t*t);
    }
    else {
        double a =  TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
        double b= n/absAlpha - absAlpha;
        return a/TMath::Power(b - t, n);
    }
}

double twotailedcrystalball(double m, double mean, double sigma, double alpha1, double alpha2, double n1, double n2){

    double result = 0.0;
    if (alpha1 < 0.0 || alpha2 < 0.0 || isnan(alpha1) || isnan(alpha2)){
        spdlog::critical("sigma {0:f}\talpha1 {1:f}\talpha2 {2:f}\tn1 {3:f}\tn2 {4:f}",sigma,alpha1,alpha2,n1,n2);
        assert(0);
    }

    //implementation taken from RooFit: RooDoubleCB.cxx
    double t = (m-mean)/sigma;
    if(t>-alpha1 && t<alpha2){
        result = TMath::Exp(-0.5*t*t);
    }
    else if(t<=-alpha1){
        double alpha1invn1 = alpha1/n1;
        result = TMath::Exp(-0.5*alpha1*alpha1)*TMath::Power(1. - alpha1invn1*(alpha1+t), -n1);
    }
    else if(t>=alpha2){
        double alpha2invn2 = alpha2/n2;
        result = TMath::Exp(-0.5*alpha2*alpha2)*TMath::Power(1. - alpha2invn2*(alpha2-t), -n2);
    }
    else{
        spdlog::critical("Cannot sort t={0:f} in range with alphas: alpha1={1:f}, alpha2={2:f}",
                         t, alpha1, alpha2);
        assert(0);
    }

    if(std::isnan(result) || std::isinf(result)){
        spdlog::critical("Calculation of twotailed CB failed: {0:f}", result);
        spdlog::critical("m(B): {0:f}\tmean: {1:f}\tsigma: {2:f}\talpha1: {3:f}\talpha2: {4:f}\tn1: {5:f}\tn2: {6:f}",m, mean, sigma, alpha1, alpha2, n1, n2);
        assert(0);
    }

    return result;
}

double sqr(double x)
{
    return x*x;
}

//basic trigonometric functions
double costheta2(double costheta){
    return costheta*costheta; //cos^2(x)
}
double cos2theta(double costheta){
    return 2.0*costheta2(costheta) - 1.0; //cos(2x) = 2cos^2(x)-1
}
double sintheta2(double costheta){
    return 1.0-costheta2(costheta); //sin^2 = 1-cos^2
}
double sintheta_abs(double costheta){ //TODO: check if this is fine as sin(x) != |sin(x)|
    return sqrt(sintheta2(costheta));
}
double sin2theta(double costheta){ //TODO: check if this is fine as sin(x) != |sin(x)|
    return 2.0 * sintheta_abs(costheta) * costheta;
}


//calculates all chebyshev polynomials up to including degree n
void chebyshev(double x, int n, std::vector<double>& results)
{
    assert(n >= 0);
    results.clear();
    results.reserve(n+1);
    results.push_back(1.0);
    results.push_back(x);
    for (int i=2; i<n+1; i++){
        results.push_back(2.0*x*results.at(i-1) - results.at(i-2));
    }
}

//get polynomial coefficients from chebyshev coefficients
void chebyshev_to_poly(const std::vector<double>& chebychev, std::vector<double>& poly)
{
    poly.clear();
    poly.resize(chebychev.size(), 0.0);
    if (chebychev.size() > 0) poly.at(0) = chebychev.at(0);
    if (chebychev.size() > 1) poly.at(1) = chebychev.at(1);
    for (unsigned int n =2; n<chebychev.size(); n++)
    {
        for (unsigned int r=0; r<=n/2; r++)
        {
            poly.at(n-2*r) += chebychev.at(n)
                    *n/2.0*pow(-1.0,r)/double(n-r)*TMath::Binomial(n-r,r)*pow(2.0,n-2*r);
        }
    }
}


//correct the polynomial coefficients from -1..1 to min..max
void correct_poly(const std::vector<double>& poly, std::vector<double>& correct, double min, double max){
    correct.clear();

    double c = 2.0/(max-min); //What idiot names stuff like c, d and e without explanations
    double d = -2.0*min/(max-min)-1.0;
    double e = d/c; //e = -min ?!

    unsigned int order = poly.size();
    correct.resize(order, 0.0);
    for (unsigned int i = 0; i<order; i++){
        double c_dash = poly.at(i);
        for (unsigned int m=0; m<=i; m++){
            correct.at(m) +=
                    c_dash
                    *pow(c,int(i))
                    *TMath::Binomial(i,m)
                    *pow(e,int(i-m));
        }
    }
}

//calculates legendre polynomial of order n
double legendre(double x, int n){
    assert(n >= 0);
    if (n == 0)      return 1;
    else if (n == 1) return x;
    else             return ((2.0*n-1.0)*x*legendre(x, n-1) - (n-1.0)*legendre(x, n-2))/double(n);
}

//calculates all legendre polynomials up to including degree n
void legendre(double x, int n, std::vector<double>& results){
    //optimized
    assert(n > 0);
    //results = std::vector<double>(n+1,1.0);
    results.resize(n+1);
    results[0] = 1.0;
    results[1] = x;
    for (int i=2; i<n+1; i++){
        results[i] = ((2.0*i-1.0) * x * results[i-1] - (i-1.0) * results[i-2]) / i;
    }
}

void orthonormal_legendre(double x, int n, std::vector<double>& results){
    legendre(x, n, results);
    for (unsigned int i=0; i<results.size(); i++){
        results.at(i) *= sqrt(0.5*(2.0*i+1.0));
    }
}

void legendre_to_poly(const std::vector<double>& legendre, std::vector<double>& poly){
    poly.clear();
    poly.resize(legendre.size(), 0.0);
    if (legendre.size() > 0){//constant
        poly.at(0) = legendre.at(0);
    }
    if (legendre.size() > 1){//linear
        poly.at(1) = legendre.at(1);
    }
    for (unsigned int i=2; i<legendre.size(); i++)//quadratic and higher
        //example size = 3 -> quadratic
        //i = 2 and 2 only
        //k=0
        //k=1
        //as it should be
        //poly(2) +=
        //poly(0) +=
    {
        for (unsigned int k=0; k<=i/2; k++)
            poly.at(i-2*k) += legendre.at(i)
                    * pow(-1.0, k)/pow(2.0, i)*TMath::Binomial(i, k)*TMath::Binomial(2*i-2*k, i);
    }

}

//calculates chebyshev polynomial of order n
double chebyshev(double x, int n){
    assert(n >= 0);
    if (n == 0)         return 1;
    else if (n == 1)    return x;
    else                return 2.0*x*chebyshev(x, n-1) - chebyshev(x, n-2);
}


double legendre_poly(unsigned int l, double x){
    return ROOT::Math::legendre(l, x);
}

double assoc_legendre(unsigned int l, unsigned int m, double x){
    return ROOT::Math::assoc_legendre(l, m, x);
}

double spherical_harmonic_re(unsigned int l, unsigned int m, double cos_theta, double phi){
    return sqrt((2*l+1.0)/(4.0*TMath::Pi()))*sqrt(TMath::Factorial(l-m)/TMath::Factorial(l+m))
            *assoc_legendre(l, m, cos_theta)
            *cos(m*phi);
}

double spherical_harmonic_im(unsigned int l, unsigned int m, double cos_theta, double phi){
    return sqrt((2*l+1.0)/(4.0*TMath::Pi()))*sqrt(TMath::Factorial(l-m)/double(TMath::Factorial(l+m)))
            *assoc_legendre(l, m, cos_theta)
            *sin(m*phi);
}

void calculate_transversity_angles(const LorentzVector& mu_minus,
                                   const LorentzVector& mu_plus,
                                   const LorentzVector& k_minus,
                                   const LorentzVector& k_plus,
                                   double & cos_theta, double & angle_phi, double & cos_psi){
    //LorentzVector bs = mu_plus + mu_minus + k_plus + k_minus;
    LorentzVector jpsi = mu_plus + mu_minus;
    LorentzVector phi = k_plus + k_minus;
    LorentzBoost jpsiboost(jpsi.BoostToCM());//boosttocm actually gives beta vector
    LorentzBoost phiboost(phi.BoostToCM());//pxyzm has boost to cm pxyzm not?
    //boosted quantities in jpsi system
    //LorentzVector mu_minus_d = jpsiboost(mu_minus);
    LorentzVector mu_plus_d = jpsiboost(mu_plus);
    LorentzVector k_minus_d = jpsiboost(k_minus);
    LorentzVector k_plus_d = jpsiboost(k_plus);
    LorentzVector phi_d = jpsiboost(phi);
    //boosted quantities in phi system
    LorentzVector k_plus_dd = phiboost(k_plus);
    LorentzVector jpsi_dd = phiboost(jpsi);

    //cos theta is calculated in the cms of the jpsi
    //cos theta
    Vector3 normalxy = k_minus_d.Vect().Cross(k_plus_d.Vect());
    normalxy = normalxy.Unit();//normals should be normalized
    cos_theta = (normalxy.Dot(mu_plus_d.Vect())
                 / sqrt(normalxy.Mag2())
                 / sqrt(mu_plus_d.Vect().Mag2()));
    //phi is calculated in the cms of the jpsi
    //cos phi
    Vector3 mu_perp = normalxy * mu_plus_d.Vect().Dot(normalxy);
    Vector3 mu_parallel = mu_plus_d.Vect() - mu_perp;
    angle_phi = acos(mu_parallel.Dot(phi_d.Vect())
                     / sqrt(mu_parallel.Mag2())
                     / sqrt(phi_d.Vect().Mag2()));
    //now get orientation
    Vector3 checknormal = phi_d.Vect().Cross(mu_parallel);
    bool samedirection = (checknormal.Dot(normalxy) > 0.0);
    if (!samedirection)
        angle_phi = -angle_phi;

    //cos psi is calculated in the cms of the phi
    //cos psi
    cos_psi = k_plus_dd.Vect().Dot(-jpsi_dd.Vect())
            /sqrt(k_plus_dd.Vect().Mag2())
            /sqrt(jpsi_dd.Vect().Mag2());

}

//this is not normalized!
double convoluted_exp(double ct, double sigma_ct, double ctau){
    //use error function instead of numeric integration, should be a lot faster
    if (ct > -6.0*sigma_ct)//this safety is important for the plotting
    {
        const double sqrt2 = TMath::Sqrt(2.0);
        //  const double norm_factor = 2.0/TMath::Sqrt(TMath::Pi());
        // Erf(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between 0 and x, already taken into account?
        double exp_part = exp(-ct/ctau + sigma_ct*sigma_ct/(2.0*ctau*ctau));
        double erfc = TMath::Erfc(sigma_ct/ctau/sqrt2 - ct/sigma_ct/sqrt2);
        // / norm_factor; the normalization here is already done?
        //the expression below was normalized, but it shouldnt have been
        //double prob = 0.5/ctau * exp * erfc;//is the factor 2 correct? yes.
        double prob = 0.5 * exp_part * erfc;//is the factor 2 correct? yes.
        if (std::isnan(prob) || prob <= 0.0)
        {
            spdlog::warn("CONVEXP gives prob: {0:f}; ct/sigma_ct/ctau/exp/erfc/exparg/erfarg: {1:f} {2:f}",
                         prob, ct,sigma_ct );
            spdlog::warn("ctau/exp{0:f} {1:f}",ctau,exp_part);
            spdlog::warn("erfc/exparg/erfarg {0:f} {1:f} {2:f}",erfc,
                         -ct/ctau + sigma_ct*sigma_ct/(2.0*ctau*ctau),
                         sigma_ct/ctau/sqrt2 - ct/sigma_ct/sqrt2 );
            prob = 0.0;
        }
        return prob;
    }
    else return 0.0;//safety added
}


//may be needed later...
double expdecay(double t, double tau, bool normalized){
    if (t >= 0.0) {
        if (normalized) return 1.0/tau*exp(-t/tau);
        else            return exp(-t/tau);
    }
    else{
        return 0.0;
    }
}

double gauss(double x, double sigma, double mean){
    //what do we have libraries for...
    //this one is normalised
    return TMath::Gaus(x, mean, sigma, kTRUE);
}

double normed_gauss(double x, double sigma, double mean, double min, double max){
    //corrected error
    double norm = 1.0/2.0
            *(TMath::Erf((max-mean)/(TMath::Sqrt(2.0)*sigma))
              -TMath::Erf((min-mean)/(TMath::Sqrt(2.0)*sigma)));
    return 1.0/(TMath::Sqrt(2.0*TMath::Pi())*sigma)*exp(-(mean-x)*(mean-x)/(2*sigma*sigma))/norm;
}

double linear(double x, double min, double max, double slope){
    //have corrected the slope (parameter needs to be big enough for minuit -> divide by width = 500^2)
    double width = max - min;
    double ymin = 1.0/width - 0.5*slope/(width*width)*width;//corrected the 0.5 factor
    return ymin + slope/(width*width) * (x - min);
}


//computes erfc, because erfc(x)=exp(-x^2)*w(i*x)
std::complex<double> cErrF(const std::complex<double>& x){
    const std::complex<double> i(0.0,1.0);
    std::complex<double> z(i*x);
    //std::complex<double> z(-x.imag(), x.real());
    std::complex<double> result = exp(-x*x)*wErrF(z);
    return result;
}

//computes erfc, because erfc(x)=exp(-x^2)*w(i*x)
std::complex<double> cErrF_2(const std::complex<double>& x){
    const std::complex<double> i(0.0,1.0);
    std::complex<double> z(i*x);
    std::complex<double> result = exp(-x*x)*Faddeeva_2(z);

    if (x.real() > 20.0)  result = 0.0;
    if (x.real() < -20.0) result = 2.0;

    return result;
}

std::complex<double> ErrF_2(const std::complex<double>& x){
    return 1.0 - cErrF_2(x);
}

//computes erfc, because erfc(x)=exp(-x^2)*w(i*x)
std::complex<double> cErrF_3(const std::complex<double>& x){
    const std::complex<double> i(0.0,1.0);
    std::complex<double> z(i*x);
    std::complex<double> result = exp(-x*x)*nwwerf(z);
    return result;
    //careful! the fortran version is apparently not thread safe //TODO
}

//from matpack
//precision 1e-14
std::complex<double> Faddeeva_2 (const std::complex<double>& z)
{
    // table 1: coefficients for h = 0.5
    static double n1[12] =
    { 0.25, 1.0, 2.25, 4.0, 6.25, 9.0, 12.25, 16.0,
      20.25, 25.0, 30.25, 36.0 };
    static double e1[12] =
    { 0.7788007830714049,    0.3678794411714423,
      1.053992245618643e-1,  1.831563888873418e-2,
      1.930454136227709e-3,  1.234098040866795e-4,
      4.785117392129009e-6,  1.125351747192591e-7,
      1.605228055185612e-9,  1.388794386496402e-11,
      7.287724095819692e-14, 2.319522830243569e-16 };

    // table 2: coefficients for h = 0.53
    static double n2[12] =
    { 0.2809, 1.1236, 2.5281, 4.4944, 7.0225, 10.1124,
      13.7641, 17.9776, 22.7529, 28.09, 33.9889, 40.4496 };
    static double e2[12] =
    { 0.7551038420890235,    0.3251072991205958,
      7.981051630007964e-2,  1.117138143353082e-2,
      0.891593719995219e-3,  4.057331392320188e-5,
      1.052755021528803e-6,  1.557498087816203e-8,
      1.313835773243312e-10, 6.319285885175346e-13,
      1.733038792213266e-15, 2.709954036083074e-18 };
    
    // tables for Pade approximation
    static double C[7] =
    { 65536.0, -2885792.0, 69973904.0, -791494704.0,
      8962513560.0, -32794651890.0, 175685635125.0 };
    static double D[7] =
    { 192192.0, 8648640.0, 183783600.0, 2329725600.0,
      18332414100.0, 84329104860.0, 175685635125.0 };

    double *n,*e,t,u,r,s,d,f,g,h;
    std::complex<double> c,d2,v,w,zz;
    int i;
    
    // use Pade approximation
    s = norm(z);
    if (s < 1e-7) {
        zz = z*z;
        v  = exp(zz);
        c  = C[0];
        d2 = D[0];
        for (i = 1; i <= 6; i++) {
            c  = c  * zz + C[i];
            d2 = d2 * zz + D[i];
        }
        w = 1.0 / v + std::complex<double>(0.0,m_2_SQRTPI) * c/d2 * z * v;
        return w;

        // use trapezoid rule
    }
    else {
        // select default table 1
        n = n1;
        e = e1;
        r = M_1_PI * 0.5;

        // if z is too close to a pole select table 2
        if (fabs(imag(z)) < 0.01 && fabs(real(z)) < 6.01) {
            h = modf(2*fabs(real(z)),&g);
            if (h < 0.02 || h > 0.98) {
                n = n2;
                e = e2;
                r = M_1_PI * 0.53;
            }
        }
        
        d = (imag(z) - real(z)) * (imag(z) + real(z));
        f = 4 * real(z) * real(z) * imag(z) * imag(z);

        g = h = 0.0;
        for (i = 0; i < 12; i++) {
            t = d + n[i];
            u = e[i] / (t * t + f);
            g += (s + n[i]) * u;
            h += (s - n[i]) * u;
        }
        u = 1 / s;
        c = r * std::complex<double>(imag(z) * (u + 2.0 * g),
                                     real(z) * (u + 2.0 * h) );
        
        if (imag(z) < M_2PI) {
            s = 2.0 / r;
            t = s * real(z);
            u = s * imag(z);
            s = sin(t);
            h = cos(t);
            f = exp(- u) - h;
            g = 2.0 * exp(d-u) / (s * s + f * f);
            u = 2.0 * real(z) * imag(z);
            h = cos(u);
            t = sin(u);
            c += g * std::complex<double>( (h * f - t * s), -(h * s + t * f));
        }
        return c;
    }
}

//this is the complex error function w(z)
std::complex<double> wErrF(const std::complex<double>& arg){
    const double c1 = 7.4;
    const double c2 = 8.3;
    const double c3 = 0.3125;
    const double c4 = 1.6;

    const double p = 46768052394588893.38251791464692105662;

    //----------------------------------------------------------------------------
    std::complex<double>  r[37];
    std::complex<double>  zh, s, t, v;
    std::complex<double>  hh, den;

    double xl = p;

    double x  = arg.real();
    double y  = arg.imag();
    double xa = fabs(x);
    double ya = fabs(y);

    if (ya < c1 && xa < c2) {
        zh = std::complex<double>(ya+c4, xa);

        for (int n=35; n>-1; n--) {
            t    = zh + std::conj(double(n+1)*r[n+1]);
            r[n] = 0.5 / (t.real()*t.real() + t.imag()*t.imag()) * t;
        }

        for (int n=32; n>-1; n--) {
            xl =  xl*c3;
            s  = (xl+s)*r[n];
        }

        v = m_2_SQRTPI * s;
    }
    else {
        zh = std::complex<double>(ya, xa);

        for (int n=8; n>-1; n--) {
            t    = zh + std::conj(double(n+1)*r[0]);
            r[0] = 0.5 / (t.real()*t.real() + t.imag()*t.imag()) * t;
        }

        v = m_2_SQRTPI * r[0];
    }

    if (ya == 0){
        v = std::complex<double>(exp(-xa*xa), v.imag());
    }

    if (y < 0) {
        hh = std::complex<double>(xa, ya);
        std::complex<double> hprod = -hh*hh;
        std::complex<double> anexp(exp(hprod.real())*cos(hprod.imag()),exp(hprod.real())*sin(hprod.imag()));
        v = 2.0 * anexp - v;
    }
    else if (x < 0){
        v = conj(v);
    }

    if(std::isnan(v.real())) spdlog::warn("wErrF: the result is NaN, do not trust anything!");

    return v;

}

std::complex<double> nwwerf(const std::complex<double> z) {
    std::complex<double>  zh,r[38],s,t,v;

    const double z1 = 1;
    const double hf = z1/2;
    const double z10 = 10;
    const double c1 = 74/z10;
    const double c2 = 83/z10;
    const double c3 = z10/32;
    const double c4 = 16/z10;
    const double p = pow(2.0*c4,33);

    double x=z.real();
    double y=z.imag();
    double xa=(x >= 0) ? x : -x;
    double ya=(y >= 0) ? y : -y;
    if(ya < c1 && xa < c2){
        zh = std::complex<double>(ya+c4,xa);
        r[37]= std::complex<double>(0,0);
        //       do 1 n = 36,1,-1
        for(int n = 36; n>0; n--){
            t=zh+double(n)*std::conj(r[n+1]);
            r[n]=hf*t/std::norm(t);
        }
        double xl=p;
        s=std::complex<double>(0,0);
        //    do 2 n = 33,1,-1
        for(int k=33; k>0; k--){
            xl=c3*xl;
            s=r[k]*(s+xl);
        }
        v=m_2_SQRTPI*s;
    }
    else{
        zh=std::complex<double>(ya,xa);
        r[1]=std::complex<double>(0,0);
        //       do 3 n = 9,1,-1
        for(int n=9;n>0;n--){
            t=zh+double(n)*std::conj(r[1]);
            r[1]=hf*t/std::norm(t);
        }
        v=m_2_SQRTPI*r[1];
    }
    if(ya == 0) v= std::complex<double>(exp(-xa*xa),v.imag());
    if(y < 0) {
        v=2.0*std::exp(std::complex<double>(-xa,-ya)*std::complex<double>(xa,ya))-v;
        if(x > 0) v=std::conj(v);
    }
    else{
        if(x < 0) v=std::conj(v);
    }
    return v;
}

void norm_convoluted_exp_sincos(double tau, double deltam, double& one, double& two){
    //double tau = 1.0/gamma();//0.5*(ctau_l() + ctau_h());//is this correct? 1/(0.5(gammal + gammah))
    double tau_2 = tau*tau;
    double deltam_2 = deltam*deltam;
    one = tau_2 * deltam / (tau_2 * deltam_2 + 1.0);//sin
    two = tau / (tau_2 * deltam_2 + 1.0);//cos
}

//calculates the gaussian convolution sin(angle) from 0 to pi
double convoluted_sin(double angle, double sigma){
    double b = -angle/sqrt(2.0)/sigma + TMath::Pi()/sqrt(2.0)/sigma;
    double c = angle/sqrt(2.0)/sigma;
    double a = sigma/sqrt(2.0);
    double sigma_2 = sigma*sigma;
    double result = 0.0;
    std::complex<double> arg1(b, a);
    std::complex<double> arg2(c, a);
    std::complex<double> erf1 = 1.0-cErrF_2(arg1);
    std::complex<double> erf2 = 1.0-cErrF_2(arg2);
    result = 0.5*sin(angle)*exp(-0.5*sigma_2)
            *(erf1.real()+erf2.real())
            +0.5*cos(angle)*exp(-0.5*sigma_2)
            *(-erf1.imag()+erf2.imag());
    spdlog::trace("convoluted sin:\t{0:f}", result);
    return result;
}

double convoluted_cos_sin(double angle, double sigma){
    return 0.5*convoluted_2_sin(angle, sigma);
}

//calculates the gaussian convolution sin(angle) from 0 to pi
double convoluted_cos(double angle, double sigma)
{
    double b = -angle/sqrt(2.0)/sigma + TMath::Pi()/sqrt(2.0)/sigma;
    double c = -angle/sqrt(2.0)/sigma;//also changed sign here
    double a = sigma/sqrt(2.0);
    double sigma_2 = sigma*sigma;
    double result = 0.0;
    std::complex<double> arg1(b, a);
    std::complex<double> arg2(c, a);
    std::complex<double> erf1 = 1.0-cErrF_2(arg1);
    std::complex<double> erf2 = 1.0-cErrF_2(arg2);
    result = 0.5*sin(angle)*exp(-0.5*sigma_2)
            *(erf1.imag()-erf2.imag())
            +0.5*cos(angle)*exp(-0.5*sigma_2)
            *(erf1.real()-erf2.real());
    return result;
}

//calculates the gaussian convolution cos(angle)^2 from 0 to pi
double convoluted_cos_2(double angle, double sigma){
    double b = -angle/sqrt(2.0)/sigma + TMath::Pi()/sqrt(2.0)/sigma;
    double c = -angle/sqrt(2.0)/sigma;//changed the sign here, more consistent
    double a = sigma*sqrt(2.0);//!!!
    double sigma_2 = sigma*sigma;
    double result = 0.0;
    std::complex<double> arg1(b, a);
    std::complex<double> arg2(c, a);
    std::complex<double> erf1 = 1.0-cErrF_2(arg1);
    std::complex<double> erf2 = 1.0-cErrF_2(arg2);
    result = 0.25*sin(2.0*angle)*exp(-2.0*sigma_2)
            *(erf1.imag()-erf2.imag())
            +0.25*cos(2.0*angle)*exp(-2.0*sigma_2)
            *(erf1.real()-erf2.real())
            -0.25*erf(-b)
            +0.25*erf(-c);
    return result;
}

//calculates the gaussian convolution cos(angle)^2*sin(angle) from 0 to pi
double convoluted_cos_2_sin(double angle, double sigma){
    double b = -angle/sqrt(2.0)/sigma + TMath::Pi()/sqrt(2.0)/sigma;
    double c = -angle/sqrt(2.0)/sigma;//changed the sign here, more consistent

    double a = 3.0*sigma/sqrt(2.0);
    double sigma_2 = sigma*sigma;
    double result = 0.0;
    std::complex<double> arg1(b, a);
    std::complex<double> arg2(c, a);
    std::complex<double> erf1 = 1.0-cErrF_2(arg1);
    std::complex<double> erf2 = 1.0-cErrF_2(arg2);

    result = 1.0/8.0*sin(3.0*angle)*exp(-4.5*sigma_2)
            *(erf1.real()-erf2.real())
            +1.0/8.0*cos(3.0*angle)*exp(-4.5*sigma_2)
            *(-erf1.imag()+erf2.imag());

    a = sigma/sqrt(2.0);
    arg1 = std::complex<double>(b, a);
    arg2 = std::complex<double>(c, a);

    erf1 = 1.0-cErrF_2(arg1);
    erf2 = 1.0-cErrF_2(arg2);

    result += 1.0/8.0*sin(angle)*exp(-0.5*sigma_2)
            *(erf1.real()-erf2.real())
            +1.0/8.0*cos(angle)*exp(-0.5*sigma_2)
            *(-erf1.imag()+erf2.imag());
    return result;
}

//gaussian convolution sin(angle)*sin(2*angle) from 0 to pi
double convoluted_2_sin_sin(double angle, double sigma){
    double b = -angle/(sqrt(2.0)*sigma) + MY_PI/(sqrt(2.0)*sigma);
    double c = -angle/(sqrt(2.0)*sigma);
    double a = 3.0*sigma/sqrt(2.0);
    double sigma_2 = sigma*sigma;
    double result = 0.0;
    std::complex<double> arg1(b, a);
    std::complex<double> arg2(c, a);
    std::complex<double> erf1 = 1.0-cErrF_2(arg1);
    std::complex<double> erf2 = 1.0-cErrF_2(arg2);
    result = 1.0/4.0*sin(3.0*angle)*exp(-4.5*sigma_2)
            *(-erf1.imag()+erf2.imag())
            +1.0/4.0*cos(3.0*angle)*exp(-4.5*sigma_2)
            *(-erf1.real()+erf2.real());

    a = sigma/sqrt(2.0);
    arg1 = std::complex<double>(b, a);
    arg2 = std::complex<double>(c, a);

    erf1 = 1.0-cErrF_2(arg1);
    erf2 = 1.0-cErrF_2(arg2);

    result += 1.0/4.0*sin(angle)*exp(-0.5*sigma_2)
            *(+erf1.imag()-erf2.imag())
            +1.0/4.0*cos(angle)*exp(-0.5*sigma_2)
            *(+erf1.real()-erf2.real());
    return result;//please recheck //I tried but gave up quickly as there is NOT A SINGLE COMMENT ANYHWERE
    //TODO
}

//gaussian convolution sin(angle)^3 from 0 to pi
double convoluted_sin_3(double angle, double sigma){
    double b = -angle/sqrt(2.0)/sigma + TMath::Pi()/sqrt(2.0)/sigma;
    double c = -angle/sqrt(2.0)/sigma;//changed the sign here, more consistent
    double a = 3.0*sigma/sqrt(2.0);
    double sigma_2 = sigma*sigma;
    double result = 0.0;
    std::complex<double> arg1(b, a);
    std::complex<double> arg2(c, a);

    std::complex<double> erf1 = 1.0-cErrF_2(arg1);
    std::complex<double> erf2 = 1.0-cErrF_2(arg2);

    result = 1.0/8.0*sin(3.0*angle)*exp(-4.5*sigma_2)
            *(-erf1.real()+erf2.real())
            +1.0/8.0*cos(3.0*angle)*exp(-4.5*sigma_2)
            *(erf1.imag()-erf2.imag());

    a = sigma/sqrt(2.0);
    arg1 = std::complex<double>(b, a);
    arg2 = std::complex<double>(c, a);

    erf1 = 1.0-cErrF_2(arg1);
    erf2 = 1.0-cErrF_2(arg2);

    result += 3.0/8.0*sin(angle)*exp(-0.5*sigma_2)
            *(erf1.real()-erf2.real())
            +3.0/8.0*cos(angle)*exp(-0.5*sigma_2)
            *(-erf1.imag()+erf2.imag());
    return result;//this should be correct now->recheck //TODO//I tried but gave up quickly as there is NOT A SINGLE COMMENT ANYHWERE
    //Also it is not used anywhere, so why bother

}

//calculates the gaussian convolution sin(angle)^2 from 0 to pi
double convoluted_sin_2(double angle, double sigma){
    double b = -angle/sqrt(2.0)/sigma + TMath::Pi()/sqrt(2.0)/sigma;
    double c = -angle/sqrt(2.0)/sigma;//changed the sign here, more consistent
    double a = sigma*sqrt(2.0);//!!!
    double sigma_2 = sigma*sigma;
    double result = 0.0;
    std::complex<double> arg1(b, a);
    std::complex<double> arg2(c, a);
    std::complex<double> erf1 = 1.0-cErrF_2(arg1);
    std::complex<double> erf2 = 1.0-cErrF_2(arg2);
    result = 0.25*sin(2.0*angle)*exp(-2.0*sigma_2)
            *(-erf1.imag()+erf2.imag())
            +0.25*cos(2.0*angle)*exp(-2.0*sigma_2)
            *(-erf1.real()+erf2.real())
            -0.25*erf(-b)
            +0.25*erf(-c);
    return result;
}

//calculates the gaussian convolution cos(2*angle) from 0 to pi
double convoluted_2_cos(double angle, double sigma){
    double b = -angle/sqrt(2.0)/sigma + TMath::Pi()/sqrt(2.0)/sigma;
    double c = -angle/sqrt(2.0)/sigma;//changed the sign here, more consistent
    double a = sigma*sqrt(2.0);//!!!
    double sigma_2 = sigma*sigma;
    double result = 0.0;
    std::complex<double> arg1(b, a);
    std::complex<double> arg2(c, a);
    std::complex<double> erf1 = 1.0-cErrF_2(arg1);
    std::complex<double> erf2 = 1.0-cErrF_2(arg2);
    result =
            0.5*sin(2.0*angle)*exp(-2.0*sigma_2)
            *(erf1.imag()-erf2.imag())
            +0.5*cos(2.0*angle)*exp(-2.0*sigma_2)
            *(erf1.real()-erf2.real());
    return result;
}

//calculates the gaussian convolution sin(2*angle) from 0 to pi
double convoluted_2_sin(double angle, double sigma){
    double b = -angle/sqrt(2.0)/sigma + TMath::Pi()/sqrt(2.0)/sigma;
    double c = -angle/sqrt(2.0)/sigma;//changed the sign here, more consistent
    double a = sigma*sqrt(2.0);//rechcek this for all other convs!
    double sigma_2 = sigma*sigma;
    double result = 0.0;
    std::complex<double> arg1(b, a);
    std::complex<double> arg2(c, a);
    std::complex<double> erf1 = 1.0-cErrF_2(arg1);
    std::complex<double> erf2 = 1.0-cErrF_2(arg2);
    result = 0.5*sin(2.0*angle)*exp(-2.0*sigma_2)
            *(erf1.real()-erf2.real())//ok
            +0.5*cos(2.0*angle)*exp(-2.0*sigma_2)
            *(-erf1.imag()+erf2.imag());

    return result;
}

//integrate(sqrt(alpha/%pi)*sin(d)*exp(-alpha*(t-d)^2),d,-inf,inf);
//integrate(1/sqrt(2*%pi)/sigma*sin(d)*exp(-(t-d)^2/2/sigma^2),d,-inf,inf); why does this not work?
double inf_convoluted_sin(double angle, double sigma){
    //%e^-(1/(4*alpha))*sin(t)
    return exp(-0.5*sigma*sigma)*sin(angle);
}

//integrate(sqrt(alpha/%pi)*cos(d)*exp(-alpha*(t-d)^2),d,-inf,inf);
double inf_convoluted_cos(double angle, double sigma){
    // %e^-(1/(4*alpha))*cos(t)
    return exp(-0.5*sigma*sigma)*cos(angle);
}

//integrate(sqrt(alpha/%pi)*cos(d)^2*exp(-alpha*(t-d)^2),d,-inf,inf);
double inf_convoluted_cos_2(double angle, double sigma){
    //%e^-(1/alpha)*cos(2*t)/2+1/2
    return 0.5*exp(-2.0*sigma*sigma)*cos(2.0*angle)+0.5;
}

//integrate(sqrt(alpha/%pi)*sin(d)^2*exp(-alpha*(t-d)^2),d,-inf,inf);
double inf_convoluted_sin_2(double angle, double sigma){
    // 1/2-%e^-(1/alpha)*cos(2*t)/2
    return 0.5-0.5*exp(-2.0*sigma*sigma)*cos(2.0*angle);
}

//integrate(sqrt(alpha/%pi)*sin(2*d)*exp(-alpha*(t-d)^2),d,-inf,inf);
double inf_convoluted_2_sin(double angle, double sigma){
    //%e^-(1/alpha)*sin(2*t)
    return exp(-2.0*sigma*sigma)*sin(2.0*angle);
}

//string(expand(integrate(integrate(sqrt(alpha/%pi)*exp(-alpha*(t-d)^2),d,-inf,inf),t)));
double int_inf_convoluted_const(double angle, double sigma){
    return angle;
}

//string(expand(integrate(integrate(sqrt(alpha/%pi)*sin(d)*exp(-alpha*(t-d)^2),d,-inf,inf),t)));
double int_inf_convoluted_sin(double angle, double sigma){
    //-%e^-(1/(4*alpha))*cos(t)
    return -exp(-0.5*sigma*sigma)*cos(angle);
}

//string(expand(integrate(integrate(sqrt(alpha/%pi)*cos(d)*exp(-alpha*(t-d)^2),d,-inf,inf),t)));
double int_inf_convoluted_cos(double angle, double sigma){
    //%e^-(1/(4*alpha))*sin(t)
    return exp(-0.5*sigma*sigma)*sin(angle);
}

//string(expand(integrate(integrate(sqrt(alpha/%pi)*cos(d)^2*exp(-alpha*(t-d)^2),d,-inf,inf),t)));
double int_inf_convoluted_cos_2(double angle, double sigma){
    //%e^-(1/alpha)*sin(2*t)/4+t/2
    return 0.25*exp(-2.0*sigma*sigma)*sin(2.0*angle)+0.5*angle;
}

//string(expand(integrate(integrate(sqrt(alpha/%pi)*sin(d)^2*exp(-alpha*(t-d)^2),d,-inf,inf),t)));
double int_inf_convoluted_sin_2(double angle, double sigma){
    //t/2-%e^-(1/alpha)*sin(2*t)/4
    return 0.5*angle-0.25*exp(-2.0*sigma*sigma)*sin(2.0*angle);
}

//string(expand(integrate(integrate(sqrt(alpha/%pi)*sin(2*d)*exp(-alpha*(t-d)^2),d,-inf,inf),t)));
double int_inf_convoluted_2_sin(double angle, double sigma){
    // -%e^-(1/alpha)*cos(2*t)/2
    return -0.5*exp(-2.0*sigma*sigma)*cos(2.0*angle);
}

//expand(integrate(integrate(sqrt(alpha/%pi)*sin(d)*exp(-alpha*(t-d)^2),d,0,%pi), t));
double int_convoluted_sin(double angle, double sigma){
    double alpha = 1.0/(2.0*sigma*sigma);
    double sa = sqrt(alpha);
    double sat = sa*angle;
    double sapi = sa*TMath::Pi();
    double overtwosa = 1.0/(2.0*sa);
    double expa = exp(-1.0/(4.0*alpha));
    double cost = cos(angle);
    double sint = sin(angle);
    double result = 0.0;
    result += 0.5*expa*cost*ErrF_2(std::complex<double>(sat - sapi, overtwosa)).real();
    result += -0.5*expa*sint*ErrF_2(std::complex<double>(sat - sapi, overtwosa)).imag();
    result += -0.5*expa*cost*ErrF_2(std::complex<double>(sat, overtwosa)).real();
    result += 0.5*expa*sint*ErrF_2(std::complex<double>(sat, overtwosa)).imag();
    result += 0.5*erf(sat - sapi);
    result += 0.5*erf(sat);
    return result;
}

double int_convoluted_cos_sin(double angle, double sigma){
    return 0.5*int_convoluted_2_sin(angle, sigma);
}

//expand(integrate(integrate(sqrt(alpha/%pi)*cos(d)^2*sin(d)*exp(-alpha*(t-d)^2),d,0,%pi), t));
double int_convoluted_cos_2_sin(double angle, double sigma){
    double alpha = 1.0/(2.0*sigma*sigma);
    double sa = sqrt(alpha);
    double sat = sa*angle;
    double sapi = sa*TMath::Pi();
    double overtwosa = 1.0/(2.0*sa);
    double threeovertwosa = 3.0/(2.0*sa);
    double expa = exp(-1.0/(4.0*alpha));
    double expninea = exp(-9.0/(4.0*alpha));
    double cost = cos(angle);
    double sint = sin(angle);
    double costhreet = cos(3.0*angle);
    double sinthreet = sin(3.0*angle);
    double result = 0.0;

    result += -(1.0/24.0)*expninea*sinthreet*ErrF_2(std::complex<double>(sat - sapi, threeovertwosa)).imag();
    result += (1.0/24.0)*expninea*costhreet*ErrF_2(std::complex<double>(sat - sapi, threeovertwosa)).real();

    result += -(1.0/8.0)*expa*sint*ErrF_2(std::complex<double>(sat - sapi, overtwosa)).imag();
    result += (1.0/8.0)*expa*cost*ErrF_2(std::complex<double>(sat - sapi, overtwosa)).real();

    result += (1.0/24.0)*expninea*sinthreet*ErrF_2(std::complex<double>(sat, threeovertwosa)).imag();
    result += -(1.0/24.0)*expninea*costhreet*ErrF_2(std::complex<double>(sat, threeovertwosa)).real();

    result += (1.0/8.0)*expa*sint*ErrF_2(std::complex<double>(sat, overtwosa)).imag();
    result += -(1.0/8.0)*expa*cost*ErrF_2(std::complex<double>(sat, overtwosa)).real();

    result += (1.0/6.0)*erf(sat-sapi);
    result += (1.0/6.0)*erf(sat);

    return result;

}

//expand(integrate(integrate(sqrt(alpha/%pi)*sin(d)*sin(2*d)*exp(-alpha*(t-d)^2),d,0,%pi), t));
double int_convoluted_2_sin_sin(double angle, double sigma){
    double alpha = 1.0/(2.0*sigma*sigma);
    double sa = sqrt(alpha);
    double sat = sa*angle;
    double sapi = sa*TMath::Pi();
    double overtwosa = 1.0/(2.0*sa);
    double threeovertwosa = 3.0/(2.0*sa);
    double expa = exp(-1.0/(4.0*alpha));
    double expninea = exp(-9.0/(4.0*alpha));
    double cost = cos(angle);
    double sint = sin(angle);
    double costhreet = cos(3.0*angle);
    double sinthreet = sin(3.0*angle);
    double result = 0.0;

    result += (1.0/12.0)*expninea*sinthreet*ErrF_2(std::complex<double>(sat - sapi, threeovertwosa)).real();
    result += (1.0/12.0)*expninea*costhreet*ErrF_2(std::complex<double>(sat - sapi, threeovertwosa)).imag();

    result += -(1.0/4.0)*expa*sint*ErrF_2(std::complex<double>(sat - sapi, overtwosa)).real();
    result += -(1.0/4.0)*expa*cost*ErrF_2(std::complex<double>(sat - sapi, overtwosa)).imag();

    result += -(1.0/12.0)*expninea*sinthreet*ErrF_2(std::complex<double>(sat, threeovertwosa)).real();
    result += -(1.0/12.0)*expninea*costhreet*ErrF_2(std::complex<double>(sat, threeovertwosa)).imag();

    result += (1.0/4.0)*expa*sint*ErrF_2(std::complex<double>(sat, overtwosa)).real();
    result += (1.0/4.0)*expa*cost*ErrF_2(std::complex<double>(sat, overtwosa)).imag();

    return result;
}

//expand(integrate(integrate(sqrt(alpha/%pi)*sin(d)^3*exp(-alpha*(t-d)^2),d,0,%pi), t));
double int_convoluted_sin_3(double angle, double sigma){
    double alpha = 1.0/(2.0*sigma*sigma);
    double sa = sqrt(alpha);
    double sat = sa*angle;
    double sapi = sa*TMath::Pi();
    double overtwosa = 1.0/(2.0*sa);
    double threeovertwosa = 3.0/(2.0*sa);
    double expa = exp(-1.0/(4.0*alpha));
    double expninea = exp(-9.0/(4.0*alpha));
    double cost = cos(angle);
    double sint = sin(angle);
    double costhreet = cos(3.0*angle);
    double sinthreet = sin(3.0*angle);
    double result = 0.0;
    result += (1.0/24.0)*expninea*sinthreet*ErrF_2(std::complex<double>(sat - sapi, threeovertwosa)).imag();
    result += -(1.0/24.0)*expninea*costhreet*ErrF_2(std::complex<double>(sat - sapi, threeovertwosa)).real();
    result += -(3.0/8.0)*expa*sint*ErrF_2(std::complex<double>(sat - sapi, overtwosa)).imag();
    result += (3.0/8.0)*expa*cost*ErrF_2(std::complex<double>(sat - sapi, overtwosa)).real();

    result += -(1.0/24.0)*expninea*sinthreet*ErrF_2(std::complex<double>(sat, threeovertwosa)).imag();
    result += (1.0/24.0)*expninea*costhreet*ErrF_2(std::complex<double>(sat, threeovertwosa)).real();
    result += +(3.0/8.0)*expa*sint*ErrF_2(std::complex<double>(sat, overtwosa)).imag();
    result += -(3.0/8.0)*expa*cost*ErrF_2(std::complex<double>(sat, overtwosa)).real();

    result += 1.0/3.0*erf(sat-sapi);
    result += 1.0/3.0*erf(sat);

    return result;
}

//expand(integrate(integrate(sqrt(alpha/%pi)*sin(d)^2*exp(-alpha*(t-d)^2),d,0,%pi), t));
double int_convoluted_sin_2(double angle, double sigma){//used by s wave

    double t = angle;
    double alpha = 1.0/(2.0*sigma*sigma);
    double sa = sqrt(alpha);
    double sat = sa*angle;
    double sapi = sa*TMath::Pi();
    double oversa = 1.0/(sa);
    double expa = exp(-1.0/(alpha));
    double costwot = cos(2.0*angle);
    double sintwot = sin(2.0*angle);
    double pi = TMath::Pi();
    double spi = sqrt(pi);
    double result = 0.0;
    result += (1.0/8.0)*expa*sintwot*ErrF_2(std::complex<double>(sat - sapi, oversa)).real();
    result += (1.0/8.0)*expa*costwot*ErrF_2(std::complex<double>(sat - sapi, oversa)).imag();

    result += -(1.0/8.0)*expa*sintwot*ErrF_2(std::complex<double>(sat, oversa)).real();
    result += -(1.0/8.0)*expa*costwot*ErrF_2(std::complex<double>(sat, oversa)).imag();

    result += -(t/4.0)*erf(sat-sapi);
    result += (pi/4.0)*erf(sat-sapi);
    result += (t/4.0)*erf(sat);

    result += -(1.0/(4.0*sa*spi))*exp(-alpha*t*t+2.0*pi*alpha*t-pi*pi*alpha);
    result += (1.0/(4.0*sa*spi))*exp(-alpha*t*t);

    return result;

}

//expand(integrate(integrate(sqrt(alpha/%pi)*sin(2*d)*exp(-alpha*(t-d)^2),d,0,%pi), t));
double int_convoluted_2_sin(double angle, double sigma){//used by s wave

    double alpha = 1.0/(2.0*sigma*sigma);
    double sa = sqrt(alpha);
    double sat = sa*angle;
    double sapi = sa*TMath::Pi();
    double oversa = 1.0/(sa);
    double expa = exp(-1.0/(alpha));
    double costwot = cos(2.0*angle);
    double sintwot = sin(2.0*angle);
    double result = 0.0;
    result += -(1.0/4.0)*expa*sintwot*ErrF_2(std::complex<double>(sat - sapi, oversa)).imag();
    result += (1.0/4.0)*expa*costwot*ErrF_2(std::complex<double>(sat - sapi, oversa)).real();

    result += (1.0/4.0)*expa*sintwot*ErrF_2(std::complex<double>(sat, oversa)).imag();
    result += -(1.0/4.0)*expa*costwot*ErrF_2(std::complex<double>(sat, oversa)).real();

    result += -(1.0/4.0)*erf(sat-sapi);
    result += (1.0/4.0)*erf(sat);

    return result;

}

//for n:0 thru 8 do print("n=",n,": ",expand(integrate(sqrt(alpha/%pi)*exp(-alpha*(t-d)^2)*cos(t)^n, t, -inf, inf)));
//for k:0 thru n/2 do print(expand(1/2^(n-1)*binomial(n,k)*exp(-1/4/alpha*(n-2*k)^2)*cos((n-2*k)*t)));
//for k:0 thru 8/2 do print(expand(1/(2^(8-1))*binomial(8,k)*exp(-1/4/alpha*(8-2*k)^2)*cos((8-2*k)*t)));
//for k:0 thru 8/2 do print(expand(1/(2^(8))*binomial(8,k)*exp(-1/4/alpha*(8-2*k)^2)*cos((8-2*k)*t)));

//for k:0 thru 8/2 do print(expand(1/2^(8)*binomial(8-1,k)*exp(-1/4/alpha*(8-2*k)^2)*cos((8-2*k)*t)));

//for n:0 thru 20 do print("n=",n,": ",string(expand(integrate(integrate(sqrt(alpha/%pi)*exp(-alpha*(t)^2)*cos(t+d)^n*sin(d), t, -inf, inf), d, 0, %pi))));
double int_convoluted_sin_cos_n(int n, double sigma){
    assert(n >= 0);
    assert(n <= 20);
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 1)
        return 0.0;
    switch (n) {
    case 0: return 2;
        break;
    case 2: return 1-exp(-1/alpha)/3;
        break;
    case 4: return -exp(-1/alpha)/3-exp(-4/alpha)/60+3.0/4.0;
        break;
    case 6: return -5*exp(-1/alpha)/16-exp(-4/alpha)/40-exp(-9/alpha)/560+5.0/8.0;
        break;
    case 8: return -7*exp(-1/alpha)/24-7*exp(-4/alpha)/240-exp(-9/alpha)/280-exp(-16/alpha)/4032+35.0/64.0;
        break;
    case 10: return -35*exp(-1/alpha)/128-exp(-4/alpha)/32-9*exp(-9/alpha)/1792-5*exp(-16/alpha)/8064-exp(-25/alpha)/25344+63.0/128.0;
        break;
    case 12: return -33*exp(-1/alpha)/128-33*exp(-4/alpha)/1024-11*exp(-9/alpha)/1792-11*exp(-16/alpha)/10752-exp(-25/alpha)/8448-exp(-36/alpha)/146432+231.0/512.0;
        break;
    case 14: return -1001*exp(-1/alpha)/4096-1001*exp(-4/alpha)/30720-143*exp(-9/alpha)/20480-13*exp(-16/alpha)/9216-91*exp(-25/alpha)/405504-7*exp(-36/alpha)/292864-exp(-49/alpha)/798720+429.0/1024.0;
        break;
    case 16: return -715*exp(-1/alpha)/3072-1001*exp(-4/alpha)/30720-39*exp(-9/alpha)/5120-65*exp(-16/alpha)/36864-35*exp(-25/alpha)/101376-15*exp(-36/alpha)/292864-exp(-49/alpha)/199680-exp(-64/alpha)/4177920+6435.0/16384.0;
        break;
    case 18: return -7293*exp(-1/alpha)/32768-663*exp(-4/alpha)/20480-663*exp(-9/alpha)/81920-17*exp(-16/alpha)/8192-85*exp(-25/alpha)/180224-51*exp(-36/alpha)/585728-51*exp(-49/alpha)/4259840-3*exp(-64/alpha)/2785280-exp(-81/alpha)/21168128+12155.0/32768.0;
        break;
    case 20: return -20995*exp(-1/alpha)/98304-4199*exp(-4/alpha)/131072-969*exp(-9/alpha)/114688-1615*exp(-16/alpha)/688128-323*exp(-25/alpha)/540672-4845*exp(-36/alpha)/37486592-19*exp(-49/alpha)/851968-19*exp(-64/alpha)/6684672-5*exp(-81/alpha)/21168128-exp(-100/alpha)/104595456+46189.0/131072.0;
        break;
    }
    return 0.0;
}

//for n:0 thru 20 do print("n=",n,": ",string(expand(integrate(integrate(sqrt(alpha/%pi)*exp(-alpha*(t)^2)*cos(t+d)^n*sin(d)^3, t, -inf, inf), d, 0, %pi))));
double int_convoluted_sin_3_cos_n(int n, double sigma){
    assert(n >= 0);
    assert(n <= 20);
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 1)
        return 0.0;
    switch (n) {
    case 0: return 4.0/3.0 ;
        break;
    case 2: return 2.0/3.0-2*exp(-1/alpha)/5 ;
        break;
    case 4: return -2*exp(-1/alpha)/5+exp(-4/alpha)/70+1.0/2.0 ;
        break;
    case 6: return -3*exp(-1/alpha)/8+3*exp(-4/alpha)/140+exp(-9/alpha)/2520+5.0/12.0 ;
        break;
    case 8: return -7*exp(-1/alpha)/20+exp(-4/alpha)/40+exp(-9/alpha)/1260+exp(-16/alpha)/36960+35.0/96.0 ;
        break;
    case 10: return -21*exp(-1/alpha)/64+3*exp(-4/alpha)/112+exp(-9/alpha)/896+exp(-16/alpha)/14784+exp(-25/alpha)/384384+21.0/64.0 ;
        break;
    case 12: return -99*exp(-1/alpha)/320+99*exp(-4/alpha)/3584+11*exp(-9/alpha)/8064+exp(-16/alpha)/8960+exp(-25/alpha)/128128+exp(-36/alpha)/3294720+77.0/256.0 ;
        break;
    case 14: return -3003*exp(-1/alpha)/10240+143*exp(-4/alpha)/5120+143*exp(-9/alpha)/92160+13*exp(-16/alpha)/84480+exp(-25/alpha)/67584+7*exp(-36/alpha)/6589440+exp(-49/alpha)/24893440+143.0/512.0 ;
        break;
    case 16: return -143*exp(-1/alpha)/512+143*exp(-4/alpha)/5120+13*exp(-9/alpha)/7680+13*exp(-16/alpha)/67584+5*exp(-25/alpha)/219648+exp(-36/alpha)/439296+exp(-49/alpha)/6223360+exp(-64/alpha)/171991040+2145.0/8192.0 ;
        break;
    case 18: return -21879*exp(-1/alpha)/81920+1989*exp(-4/alpha)/71680+221*exp(-9/alpha)/122880+51*exp(-16/alpha)/225280+255*exp(-25/alpha)/8200192+17*exp(-36/alpha)/4392960+9*exp(-49/alpha)/23429120+9*exp(-64/alpha)/343982080+exp(-81/alpha)/1111326720+12155.0/49152.0 ;
        break;
    case 20: return -4199*exp(-1/alpha)/16384+12597*exp(-4/alpha)/458752+323*exp(-9/alpha)/172032+323*exp(-16/alpha)/1261568+323*exp(-25/alpha)/8200192+323*exp(-36/alpha)/56229888+57*exp(-49/alpha)/79659008.0+exp(-64/alpha)/14483456.0+exp(-81/alpha)/222265344.0+exp(-100/alpha)/6816137216.0+46189.0/196608.0 ;
        break;
    }
    return 0.0;
}

//for n:0 thru 20 do print("n=",n,": ",string(expand(integrate(integrate(sqrt(alpha/%pi)*exp(-alpha*(t)^2)*cos(t+d)^n*sin(d)*cos(d)^2, t, -inf, inf), d, 0, %pi))));
double int_convoluted_cos_2_sin_cos_n(int n, double sigma){
    assert(n >= 0);
    assert(n <= 20);
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 1)
        return 0.0;
    switch (n) {
    case 0: return 2.0/3.0 ;
        break;
    case 2: return exp(-1/alpha)/15+1.0/3.0 ;
        break;
    case 4: return exp(-1/alpha)/15-13*exp(-4/alpha)/420+1.0/4.0 ;
        break;
    case 6: return exp(-1/alpha)/16-13*exp(-4/alpha)/280-11*exp(-9/alpha)/5040+5.0/24.0 ;
        break;
    case 8: return 7*exp(-1/alpha)/120-13*exp(-4/alpha)/240-11*exp(-9/alpha)/2520-61*exp(-16/alpha)/221760+35.0/192.0 ;
        break;
    case 10: return 7*exp(-1/alpha)/128-13*exp(-4/alpha)/224-11*exp(-9/alpha)/1792-61*exp(-16/alpha)/88704-97*exp(-25/alpha)/2306304+21.0/128.0 ;
        break;
    case 12: return 33*exp(-1/alpha)/640-429*exp(-4/alpha)/7168-121*exp(-9/alpha)/16128-61*exp(-16/alpha)/53760-97*exp(-25/alpha)/768768-47*exp(-36/alpha)/6589440+77.0/512.0 ;
        break;
    case 14: return 1001*exp(-1/alpha)/20480-1859*exp(-4/alpha)/30720-1573*exp(-9/alpha)/184320-793*exp(-16/alpha)/506880-97*exp(-25/alpha)/405504-329*exp(-36/alpha)/13178880-193*exp(-49/alpha)/149360640+143.0/1024.0 ;
        break;
    case 16: return 143*exp(-1/alpha)/3072-1859*exp(-4/alpha)/30720-143*exp(-9/alpha)/15360-793*exp(-16/alpha)/405504-485*exp(-25/alpha)/1317888-47*exp(-36/alpha)/878592-193*exp(-49/alpha)/37340160-253*exp(-64/alpha)/1031946240+2145.0/16384.0 ;
        break;
    case 18: return 7293*exp(-1/alpha)/163840-8619*exp(-4/alpha)/143360-2431*exp(-9/alpha)/245760-1037*exp(-16/alpha)/450560-8245*exp(-25/alpha)/16400384-799*exp(-36/alpha)/8785920.0-579*exp(-49/alpha)/46858240.0-759*exp(-64/alpha)/687964160.0-107*exp(-81/alpha)/2222653440.0+12155.0/98304.0 ;
        break;
    case 20: return 4199*exp(-1/alpha)/98304-54587*exp(-4/alpha)/917504-3553*exp(-9/alpha)/344064-19703*exp(-16/alpha)/7569408-31331*exp(-25/alpha)/49201152-15181*exp(-36/alpha)/112459776.0-3667*exp(-49/alpha)/159318016.0-253*exp(-64/alpha)/86900736.0-107*exp(-81/alpha)/444530688.0-397*exp(-100/alpha)/40896823296.0+46189.0/393216.0 ;
        break;
    }
    return 0.0;
}

//for n:0 thru 20 do print("n=",n,": ",string(expand(integrate(integrate(sqrt(alpha/%pi)*exp(-alpha*(t)^2)*cos(t+d)^n*sin(2*d)*sin(d), t, -inf, inf), d, 0, %pi))));
double int_convoluted_2_sin_sin_cos_n(int n, double sigma){
    assert(n >= 0);
    assert(n <= 20);
    double pi = TMath::Pi();
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 0)
        return 0.0;
    switch (n) {
    case 1: return pi*exp(-1/(4*alpha))/4 ;
        break;
    case 3: return 3*pi*exp(-1/(4*alpha))/16-pi*exp(-9/(4*alpha))/16 ;
        break;
    case 5: return 5*pi*exp(-1/(4*alpha))/32-5*pi*exp(-9/(4*alpha))/64 ;
        break;
    case 7: return 35*pi*exp(-1/(4*alpha))/256-21*pi*exp(-9/(4*alpha))/256 ;
        break;
    case 9: return 63*pi*exp(-1/(4*alpha))/512-21*pi*exp(-9/(4*alpha))/256 ;
        break;
    case 11: return 231*pi*exp(-1/(4*alpha))/2048-165*pi*exp(-9/(4*alpha))/2048 ;
        break;
    case 13: return 429*pi*exp(-1/(4*alpha))/4096-1287*pi*exp(-9/(4*alpha))/16384 ;
        break;
    case 15: return 6435*pi*exp(-1/(4*alpha))/65536-5005*pi*exp(-9/(4*alpha))/65536 ;
        break;
    case 17: return 12155*pi*exp(-1/(4*alpha))/131072-2431*pi*exp(-9/(4*alpha))/32768 ;
        break;
    case 19: return 46189*pi*exp(-1/(4*alpha))/524288-37791*pi*exp(-9/(4*alpha))/524288 ;
        break;
    }
    return 0.0;
}

//for n:0 thru 20 do print("n=",n,": ",string(expand(integrate(integrate(sqrt(alpha/%pi)*exp(-alpha*(t)^2)*cos(t+d)^n*sin(d)*cos(d), t, -inf, inf), d, 0, %pi))));
double int_convoluted_cos_sin_cos_n(int n, double sigma){
    assert(n >= 0);
    assert(n <= 20);
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 0)
        return 0.0;
    switch (n) {
    case 1: return 2*exp(-1/(4*alpha))/3 ;
        break;
    case 3: return exp(-1/(4*alpha))/2-exp(-9/(4*alpha))/10 ;
        break;
    case 5: return 5*exp(-1/(4*alpha))/12-exp(-9/(4*alpha))/8-exp(-25/(4*alpha))/168 ;
        break;
    case 7: return 35*exp(-1/(4*alpha))/96-21*exp(-9/(4*alpha))/160-exp(-25/(4*alpha))/96-exp(-49/(4*alpha))/1440 ;
        break;
    case 9: return 21*exp(-1/(4*alpha))/64-21*exp(-9/(4*alpha))/160-3*exp(-25/(4*alpha))/224-exp(-49/(4*alpha))/640-exp(-81/(4*alpha))/9856 ;
        break;
    case 11: return 77*exp(-1/(4*alpha))/256-33*exp(-9/(4*alpha))/256-55*exp(-25/(4*alpha))/3584-11*exp(-49/(4*alpha))/4608-exp(-81/(4*alpha))/3584-exp(-121/(4*alpha))/59904 ;
        break;
    case 13: return 143*exp(-1/(4*alpha))/512-1287*exp(-9/(4*alpha))/10240-715*exp(-25/(4*alpha))/43008-143*exp(-49/(4*alpha))/46080-39*exp(-81/(4*alpha))/78848-exp(-121/(4*alpha))/18432-exp(-169/(4*alpha))/337920 ;
        break;
    case 15: return 2145*exp(-1/(4*alpha))/8192-1001*exp(-9/(4*alpha))/8192-143*exp(-25/(4*alpha))/8192-91*exp(-49/(4*alpha))/24576-65*exp(-81/(4*alpha))/90112-35*exp(-121/(4*alpha))/319488-exp(-169/(4*alpha))/90112-exp(-225/(4*alpha))/1810432 ;
        break;
    case 17: return 12155*exp(-1/(4*alpha))/49152-2431*exp(-9/(4*alpha))/20480-221*exp(-25/(4*alpha))/12288-1547*exp(-49/(4*alpha))/368640-85*exp(-81/(4*alpha))/90112-85*exp(-121/(4*alpha))/479232-17*exp(-169/(4*alpha))/675840-exp(-225/(4*alpha))/425984-exp(-289/(4*alpha))/9338880 ;
        break;
    case 19: return 46189*exp(-1/(4*alpha))/196608-37791*exp(-9/(4*alpha))/327680-4199*exp(-25/(4*alpha))/229376-2261*exp(-49/(4*alpha))/491520-2907*exp(-81/(4*alpha))/2523136-323*exp(-121/(4*alpha))/1277952-323*exp(-169/(4*alpha))/7208960-171*exp(-225/(4*alpha))/28966912-exp(-289/(4*alpha))/1966080-exp(-361/(4*alpha))/46792704 ;
        break;
    }
    return 0.0;
}

//for n:0 thru 20 do print("n=",n,": ",string(expand(integrate(integrate(sqrt(alpha/%pi)*exp(-alpha*(t)^2)*cos(t+d)^n*sin(d)^2, t, -inf, inf), d, 0, %pi))));
double int_convoluted_sin_2_cos_n(int n, double sigma){
    assert(n >= 0);
    assert(n <= 20);
    double pi = TMath::Pi();
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 1)
        return 0.0;
    switch (n) {
    case 0: return pi/2 ;
        break;
    case 2: return pi/4-pi*exp(-1/alpha)/8 ;
        break;
    case 4: return 3*pi/16-pi*exp(-1/alpha)/8 ;
        break;
    case 6: return 5*pi/32-15*pi*exp(-1/alpha)/128 ;
        break;
    case 8: return 35*pi/256-7*pi*exp(-1/alpha)/64 ;
        break;
    case 10: return 63*pi/512-105*pi*exp(-1/alpha)/1024 ;
        break;
    case 12: return 231*pi/2048-99*pi*exp(-1/alpha)/1024 ;
        break;
    case 14: return 429*pi/4096-3003*pi*exp(-1/alpha)/32768 ;
        break;
    case 16: return 6435*pi/65536-715*pi*exp(-1/alpha)/8192 ;
        break;
    case 18: return 12155*pi/131072-21879*pi*exp(-1/alpha)/262144 ;
        break;
    case 20: return 46189*pi/524288-20995*pi*exp(-1/alpha)/262144 ;
        break;
    }
    return 0.0;
}

//for n:0 thru 20 do print("case ",n,": return ",string(float(expand(integrate(integrate(sqrt(alpha/%pi)*(phi+d)^n*exp(-alpha*phi^2), phi, -inf, inf), d, -%pi, %pi)))), "; break;");
double int_inf_convoluted_const_x_n(int n, double sigma){
    assert(n >= 0);
    assert(n <= 20);
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 1)
        return 0.0;
    switch (n) {
    case  0 : return  6.283185307179586 ; break;
    case  2 : return  3.141592653589793/alpha+20.67085112019988 ; break;
    case  4 : return  62.01255336059963/alpha+4.71238898038469/pow(alpha,2)+122.4078739141126 ; break;
    case  6 : return  918.0590543558442/alpha+232.5470751022486/pow(alpha,2)+11.78097245096172/pow(alpha,3)+862.9409222219404 ; break;
    case  8 : return  12081.17291110717/alpha+6426.41338049091/pow(alpha,2)+1085.219683810494/pow(alpha,3)+41.23340357836604/pow(alpha,4)+6624.244296321378 ; break;
    case  10 : return  149045.496667231/alpha+135913.1952499556/pow(alpha,2)+48198.10035368182/pow(alpha,3)+6104.360721434026/pow(alpha,4)+185.5503161026471/pow(alpha,5)+53491.63963161646 ; break;
    case  12 : return
                1765224.107843343/alpha+2459250.695009312/pow(alpha,2)+1495045.147749512/pow(alpha,3)+397634.327917875/pow(alpha,4)+40288.78076146457/pow(alpha,5)+1020.526738564559/pow(alpha,6)+446719.5800943512
                ; break;
    case  14 : return  2.0325740894292977E+7/alpha+4.0158848453436054E+7/pow(alpha,2)+3.7298635540974565E+7/pow(alpha,3)+1.7006138555650696E+7/pow(alpha,4)+3618472.384052663/pow(alpha,5)+305523.254107773/pow(alpha,6)+6633.423800669636/pow(alpha,7)+3821086.129251732 ; break;
    case  16 : return  2.2926516775510389E+8/alpha+6.0977222682878935E+8/pow(alpha,2)+8.0317696906872118E+8/pow(alpha,3)+5.5947953311461842E+8/pow(alpha,4)+2.0407366266780835E+8/pow(alpha,5)+3.6184723840526626E+7/pow(alpha,6)+2618770.749495197/pow(alpha,7)+49750.67850502227/pow(alpha,8)+3.3275831010180339E+7 ; break;
    case  18 : return  2.5456010722787962E+9/alpha+8.7693926666327229E+9/pow(alpha,2)+1.5549191784134127E+10/pow(alpha,3)+1.5360759533439291E+10/pow(alpha,4)+8.5600368566536617E+9/pow(alpha,5)+2.6019391990145564E+9/pow(alpha,6)+3.9544733911432672E+8/pow(alpha,7)+2.5041995292047825E+7/pow(alpha,8)+422880.7672926893/pow(alpha,9)+2.9384883679977304E+8 ; break;
    case  20 : return  2.7915639495978439E+10/alpha+1.2091605093324281E+11/pow(alpha,2)+2.7769743444336957E+11/pow(alpha,3)+3.6929330487318555E+11/pow(alpha,4)+2.9185443113534656E+11/pow(alpha,5)+1.3553391689701631E+11/pow(alpha,6)+3.5312031986626122E+10/pow(alpha,7)+4.6959371519826298E+9/pow(alpha,8)+2.6433217252717146E+8/pow(alpha,9)+4017367.289280548/pow(alpha,10)+2.623964937416502E+9
                ; break;
    }
    return 0.0;
}

//for n:0 thru 20 do print("case ",n,": return ",string(float(expand(integrate(integrate(sqrt(alpha/%pi)*(phi+d)^n*exp(-alpha*phi^2)*cos(d)^2, phi, -inf, inf), d, -%pi, %pi)))), "; break;");
double int_inf_convoluted_cos_2_x_n(int n, double sigma)
{
    assert(n >= 0);
    assert(n <= 20);
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 1)
        return 0.0;
    switch (n) {
    case  0 : return  3.141592653589793 ; break;
    case  2 : return  1.570796326794897/alpha+11.90622188689484 ; break;
    case  4 : return  35.71866566068451/alpha+2.356194490192345/pow(alpha,2)+87.4978246569714 ; break;
    case  6 : return  656.2336849272856/alpha+133.9449962275669/pow(alpha,2)+5.890486225480862/pow(alpha,3)+693.2958305395289 ; break;
    case  8 : return  9706.141627553405/alpha+4593.635794490999/pow(alpha,2)+625.0766490619789/pow(alpha,3)+20.61670178918302/pow(alpha,4)+5687.153431714451 ; break;
    case  10 : return  127960.9522135751/alpha+109194.0933099758/pow(alpha,2)+34452.26845868249/pow(alpha,3)+3516.056150973631/pow(alpha,4)+92.77515805132357/pow(alpha,5)+47830.36426946412 ; break;
    case  12 : return  1578402.020892315/alpha+2111355.71152399/pow(alpha,2)+1201135.026409734/pow(alpha,3)+284231.2147841306/pow(alpha,4)+23205.97059642596/pow(alpha,5)+510.2633692822797/pow(alpha,6)+410181.8769982042
                ; break;
    case  14 : return  1.8663275403418235E+7/alpha+3.5908645975300178E+7/pow(alpha,2)+3.2022228291447181E+7/pow(alpha,3)+1.3662910925410721E+7/pow(alpha,4)+2586504.054535588/pow(alpha,5)+175978.6103562302/pow(alpha,6)+3316.711900334818/pow(alpha,7)+3573008.555500609 ; break;
    case  16 : return  2.1438051333003667E+8/alpha+5.5989826210254693E+8/pow(alpha,2)+7.1817291950600362E+8/pow(alpha,3)+4.8033342437170768E+8/pow(alpha,4)+1.6395493110492867E+8/pow(alpha,5)+2.5865040545355879E+7/pow(alpha,6)+1508388.088767688/pow(alpha,7)+24875.33925251113/pow(alpha,8)+3.1522569930157386E+7 ; break;
    case  18 : return  2.4114765996571231E+9/alpha+8.2000546348738604E+9/pow(alpha,2)+1.4277405683614956E+10/pow(alpha,3)+1.3735057085552319E+10/pow(alpha,4)+7.3491013928871279E+9/pow(alpha,5)+2.0904253715878406E+9/pow(alpha,6)+2.8266794310281783E+8/pow(alpha,7)+1.4423961098841013E+7/pow(alpha,8)+211440.3836463447/pow(alpha,9)+2.810488910215596E+8 ; break;
    case  20 : return  2.6699644647049049E+10/alpha+1.1454513848371289E+11/pow(alpha,2)+2.5966839677100656E+11/pow(alpha,3)+3.390883849858551E+11/pow(alpha,4)+2.6096608462549408E+11/pow(alpha,5)+1.1636077205404617E+11/pow(alpha,6)+2.8370058614406403E+10/pow(alpha,7)+3.3566818243459616E+9/pow(alpha,8)+1.5225292270998847E+8/pow(alpha,9)+2008683.644640274/pow(alpha,10)+2.527977317637641E+9
                ; break;
    }
    return 0.0;
}

//for n:0 thru 20 do print("case ",n,": return ",string(float(expand(integrate(integrate(sqrt(alpha/%pi)*(phi+d)^n*exp(-alpha*phi^2)*sin(d)^2, phi, -inf, inf), d, -%pi, %pi)))), "; break;");
double int_inf_convoluted_sin_2_x_n(int n, double sigma){
    assert(n >= 0);
    assert(n <= 20);
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 1)
        return 0.0;
    switch (n) {
    case  0 : return  3.141592653589793 ; break;
    case  2 : return  1.570796326794897/alpha+8.764629233305042 ; break;
    case  4 : return  26.29388769991513/alpha+2.356194490192345/pow(alpha,2)+34.91004925714116 ; break;
    case  6 : return  261.8253694285586/alpha+98.60207887468172/pow(alpha,2)+5.890486225480862/pow(alpha,3)+169.6450916824115 ; break;
    case  8 : return  2375.031283553761/alpha+1832.777585999911/pow(alpha,2)+460.1430347485147/pow(alpha,3)+20.61670178918302/pow(alpha,4)+937.0908646069279 ; break;
    case  10 : return  21084.54445365589/alpha+26719.10193997981/pow(alpha,2)+13745.83189499933/pow(alpha,3)+2588.304570460395/pow(alpha,4)+92.77515805132357/pow(alpha,5)+5661.275362152337 ; break;
    case  12 : return
                186822.0869510286/alpha+347894.9834853215/pow(alpha,2)+293910.121339778/pow(alpha,3)+113403.1131337445/pow(alpha,4)+17082.81016503861/pow(alpha,5)+510.2633692822797/pow(alpha,6)+36537.70309614705
                ; break;
    case  14 : return  1662465.490874743/alpha+4250202.478135873/pow(alpha,2)+5276407.249527384/pow(alpha,3)+3343227.630239975/pow(alpha,4)+1031968.329517075/pow(alpha,5)+129544.6437515428/pow(alpha,6)+3316.711900334818/pow(alpha,7)+248077.5737511227 ; break;
    case  16 : return  1.4884654425067216E+7/alpha+4.9873964726242363E+7/pow(alpha,2)+8.5004049562717617E+7/pow(alpha,3)+7.9146108742910743E+7/pow(alpha,4)+4.0118731562879689E+7/pow(alpha,5)+1.0319683295170747E+7/pow(alpha,6)+1110382.66072751/pow(alpha,7)+24875.33925251113/pow(alpha,8)+1753261.080022954 ; break;
    case  18 : return  1.3412447262167311E+8/alpha+5.693380317588625E+8/pow(alpha,2)+1.2717861005191698E+9/pow(alpha,3)+1.6257024478869734E+9/pow(alpha,4)+1.2109354637665339E+9/pow(alpha,5)+5.1151382742671597E+8/pow(alpha,6)+1.1277939601150887E+8/pow(alpha,7)+1.0618034193206811E+7/pow(alpha,8)+211440.3836463447/pow(alpha,9)+1.2799945778213412E+7 ; break;
    case  20 : return  1.21599484892939E+9/alpha+6.3709124495299149E+9/pow(alpha,2)+1.8029037672363007E+10/pow(alpha,3)+3.0204919887330444E+10/pow(alpha,4)+3.0888346509852478E+10/pow(alpha,5)+1.9173144842970131E+10/pow(alpha,6)+6.941973372219717E+9/pow(alpha,7)+1.339255327636668E+9/pow(alpha,8)+1.1207924981718299E+8/pow(alpha,9)+2008683.644640274/pow(alpha,10)+9.5987619778861046E+7 ; break;
    }
    return 0.0;
}

//for n:0 thru 20 do print("case ",n,": return ",string(float(expand(integrate(integrate(sqrt(alpha/%pi)*(phi+d)^n*exp(-alpha*phi^2)*sin(d), phi, -inf, inf), d, -%pi, %pi)))), "; break;");
double int_inf_convoluted_sin_x_n(int n, double sigma){
    assert(n >= 0);
    assert(n <= 20);
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 0)
        return 0.0;
    switch (n) {
    case  1 : return  6.283185307179586 ; break;
    case  3 : return  9.424777960769379/alpha+24.31344151752212 ; break;
    case  5 : return  121.5672075876106/alpha+23.56194490192345/pow(alpha,2)+125.7705392201204 ; break;
    case  7 : return  1320.590661811266/alpha+638.2278398349555/pow(alpha,2)+82.46680715673207/pow(alpha,3)+758.2238083085194 ; break;
    case  9 : return  13648.02854955342/alpha+11885.31595630138/pow(alpha,2)+3829.367039009733/pow(alpha,3)+371.1006322052943/pow(alpha,4)+5026.084468678731 ; break;
    case  11 : return  138217.3228886668/alpha+187660.3925563593/pow(alpha,2)+108948.7295994294/pow(alpha,3)+26326.89839319191/pow(alpha,4)+2041.053477129119/pow(alpha,5)+35538.744393114 ; break;
    case  13 : return
                1386011.031330336/alpha+2695237.796329141/pow(alpha,2)+2439585.103232659/pow(alpha,3)+1062250.113594437/pow(alpha,4)+205349.807466897/pow(alpha,5)+13266.84760133927/pow(alpha,6)+263310.4159052214
                ; break;
    case  15 : return  1.3823796834663689E+7/alpha+3.6382789572466373E+7/pow(alpha,2)+4.7166661435754895E+7/pow(alpha,3)+3.201955447992897E+7/pow(alpha,4)+1.1153626192741588E+7/pow(alpha,5)+1796810.815335348/pow(alpha,6)+99501.35701004454/pow(alpha,7)+2021104.600121215 ; break;
    case  17 : return  1.3743511271761036E+8/alpha+4.7000909238989449E+8/pow(alpha,2)+8.2467656364292908E+8/pow(alpha,3)+8.018332444078331E+8/pow(alpha,4)+4.3546594092703295E+8/pow(alpha,5)+1.2640776351773798E+8/pow(alpha,6)+1.7454733634686239E+7/pow(alpha,7)+845761.5345853786/pow(alpha,8)+1.5948676302624345E+7 ; break;
    case  19 : return  1.363611830458992E+9/alpha+5.8753510678547668E+9/pow(alpha,2)+1.3395259133206604E+10/pow(alpha,3)+1.7627461547833862E+10/pow(alpha,4)+1.3711348479374664E+10/pow(alpha,5)+6.2053896582101974E+9/pow(alpha,6)+1.5439805401095133E+9/pow(alpha,7)+1.865474657207092E+8/pow(alpha,8)+8034734.578561097/pow(alpha,9)+1.2868057735972023E+8 ; break;
    }
    return 0.0;
}

//for n:0 thru 20 do print("case ",n,": return ",string(float(expand(integrate(integrate(sqrt(alpha/%pi)*(phi+d)^n*exp(-alpha*phi^2)*sin(2*d), phi, -inf, inf), d, -%pi, %pi)))), "; break;");
double int_inf_convoluted_2_sin_x_n(int n, double sigma){
    assert(n >= 0);
    assert(n <= 20);
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 0)
        return 0.0;
    switch (n) {
    case  1 : return  -3.141592653589793 ; break;
    case  3 : return  -4.71238898038469/alpha-26.29388769991513 ; break;
    case  5 : return  -131.4694384995756/alpha-11.78097245096172/pow(alpha,2)-174.5502462857058 ; break;
    case  7 : return  -1832.777585999911/alpha-690.214552122772/pow(alpha,2)-41.23340357836604/pow(alpha,3)-1187.515641776881 ; break;
    case  9 : return  -21375.28155198386/alpha-16494.99827399919/pow(alpha,2)-4141.287312736632/pow(alpha,3)-185.5503161026471/pow(alpha,4)-8433.817781462349 ; break;
    case  11 : return  -231929.9889902145/alpha-293910.121339778/pow(alpha,2)-151204.1508449926/pow(alpha,3)-28471.35027506435/pow(alpha,4)-1020.526738564559/pow(alpha,5)-62274.02898367599 ; break;
    case  13 : return  -2428687.130363355/alpha-4522634.785309188/pow(alpha,2)-3820831.577417114/pow(alpha,3)-1474240.470738678/pow(alpha,4)-222076.5321455019/pow(alpha,5)-6633.423800669636/pow(alpha,6)-474990.140249928
                ; break;
    case  15 : return  -2.4936982363121182E+7/alpha-6.3753037172038078E+7/pow(alpha,2)-7.9146108742910743E+7/pow(alpha,3)-5.0148414453599617E+7/pow(alpha,4)-1.5479524942756122E+7/pow(alpha,5)-1943169.656273142/pow(alpha,6)-49750.67850502227/pow(alpha,7)-3721163.606266804 ; break;
    case  17 : return  -2.5303912522615099E+8/alpha-8.4785740034611607E+8/pow(alpha,2)-1.4450688425662041E+9/pow(alpha,3)-1.3454838486294813E+9/pow(alpha,4)-6.8201843656895471E+8/pow(alpha,5)-1.7543461601790273E+8/pow(alpha,6)-1.8876505232367665E+7/pow(alpha,7)-422880.7672926893/pow(alpha,8)-2.9805438360381901E+7 ; break;
    case  19 : return  -2.5483649798125954E+9/alpha-1.0817422603417984E+10/pow(alpha,2)-2.4163935909864288E+10/pow(alpha,3)-3.0888346509852478E+10/pow(alpha,4)-2.3007773811564163E+10/pow(alpha,5)-9.7187627211076069E+9/pow(alpha,6)-2.1428085242186687E+9/pow(alpha,7)-2.017426496709294E+8/pow(alpha,8)-4017367.289280548/pow(alpha,9)-2.4319896978524876E+8 ; break;
    }
    return 0.0;
}

//for n:0 thru 20 do print("case ",n,": return ",string(float(expand(integrate(integrate(sqrt(alpha/%pi)*(phi+d)^n*exp(-alpha*phi^2)*cos(d), phi, -inf, inf), d, -%pi, %pi)))), "; break;");
double int_inf_convoluted_cos_x_n(int n, double sigma){
    assert(n >= 0);
    assert(n <= 20);
    double alpha = 1.0/(2.0*sigma*sigma);
    if (n % 2 == 1)
        return 0.0;
    switch (n) {
    case  0 : return  0.0 ; break;
    case  2 : return  -12.56637061435917 ; break;
    case  4 : return  -37.69911184307752/alpha-97.25376607008846 ; break;
    case  6 : return  -729.4032455256634/alpha-141.3716694115407/pow(alpha,2)-754.6232353207233 ; break;
    case  8 : return  -10564.72529449013/alpha-5105.822718679644/pow(alpha,2)-659.7344572538566/pow(alpha,3)-6065.790466468156 ; break;
    case  10 : return  -136480.2854955334/alpha-118853.1595630139/pow(alpha,2)-38293.67039009733/pow(alpha,3)-3711.006322052943/pow(alpha,4)-50260.8446867906 ; break;
    case  12 : return  -1658607.874663904/alpha-2251924.710676324/pow(alpha,2)-1307384.755193153/pow(alpha,3)-315922.780718303/pow(alpha,4)-24492.64172554942/pow(alpha,5)-426464.9327177554 ; break;
    case  14 : return  -1.9404154438628018E+7/alpha-3.7733329148607552E+7/pow(alpha,2)-3.4154191445257187E+7/pow(alpha,3)-1.4871501590322122E+7/pow(alpha,4)-2874897.304536558/pow(alpha,5)-185735.8664187498/pow(alpha,6)-3686345.822659835 ; break;
    case  16 : return  -2.2118074935461903E+8/alpha-5.8212463315946198E+8/pow(alpha,2)-7.5466658297207832E+8/pow(alpha,3)-5.1231287167886353E+8/pow(alpha,4)-1.784580190838654E+8/pow(alpha,5)-2.8748973045365565E+7/pow(alpha,6)-1592021.712160713/pow(alpha,7)-3.233767360193944E+7 ; break;
    case  18 : return  -2.4738320296894684E+9/alpha-8.4601636629215393E+9/pow(alpha,2)-1.4844178145556E+10/pow(alpha,3)-1.4432998399342041E+10/pow(alpha,4)-7.8383869366865463E+9/pow(alpha,5)-2.2753397433192844E+9/pow(alpha,6)-3.1418520542435235E+8/pow(alpha,7)-1.5223707622536814E+7/pow(alpha,8)-2.8707617035731125E+8 ; break;
    case  20 : return  -2.7272236397462036E+10/alpha-1.1750702138356006E+11/pow(alpha,2)-2.6790518266332812E+11/pow(alpha,3)-3.5254923095672754E+11/pow(alpha,4)-2.742269695875127E+11/pow(alpha,5)-1.2410779316420337E+11/pow(alpha,6)-3.0879610802190308E+10/pow(alpha,7)-3.7309493144141827E+9/pow(alpha,8)-1.6069469157122192E+8/pow(alpha,9)-2.5736123940656128E+9 ; break;
    }
    return 0.0;
}

//void fitter::ct_tagging_part(measurement* meas, double& one, double& two)
//this function calculates the convolution of the product exp(-t/tau) * sin(deltam*t) and exp(-t/tau) * cos(deltam*t)
void convoluted_exp_sincos(double ct, double sigma, double tau, double deltam, double& one, double& two){
    double sigma_2 = sigma*sigma;
    double tau_2 = tau*tau;
    double deltam_2 = deltam*deltam;
    const double sqrt2 = TMath::Sqrt(2.0);

    const bool cdf_way = false;
    if (!cdf_way){
        //safety for large negative decay times
        if (ct < - 6.0*sigma){
            one = 0.0;
            two = 0.0;
            return;
        }
        double c1 = 0.5;
        //      double exp1arg = 0.5*sigma_2/tau_2 - 0.5*deltam_2*sigma_2 - ct/tau;
        double exp1arg = 0.5*sigma_2*(1.0/tau_2 - deltam_2) - ct/tau;
        double exp1 = exp(exp1arg);
        double exp2arg = -deltam*(ct - sigma_2/tau);
        std::complex<double> exp2(cos(exp2arg), sin(exp2arg));
        std::complex<double> cerfarg(sigma/(tau*sqrt2) - ct/(sigma*sqrt2), +deltam*sigma/sqrt2);
        std::complex<double> cerf;
        if  (cerfarg.real() < -20.0) cerf = std::complex<double>(2.0,0.0);
        else                         cerf = cErrF_2(cerfarg);//best std::complex error function
        std::complex<double> c2(exp2*cerf);
        double im = -c2.imag();//exp*sin
        double re = +c2.real();//exp*cos
        one = c1*exp1*im ;
        two = c1*exp1*re ;
    }
    else{ //all three methods are analytically the same. There might be numerical differences, though
        if (ct < -6.0*sigma) {
            one = 0.0;
            two = 0.0;
            return;
        }

        //current cdf
        std::complex<double> z(deltam*sigma/sqrt2, (sigma/tau-ct/sigma)/sqrt2);
        if (ct<0) {//i do not quite get this part. the calculation should also be corerct for ct < 0, right?
            one= 2.0*nwwerf(z).real()/4.0*exp(-ct*ct/2.0/sigma/sigma);
            two= 2.0*nwwerf(z).imag()/4.0*exp(-ct*ct/2.0/sigma/sigma);
        }
        else {
            one= -2.0*nwwerf(std::conj(z)).real()/tau/4*exp(-ct*ct/2.0/sigma/sigma) +
                    exp(sigma*sigma/2 *(1/tau/tau - deltam*deltam) - ct/tau)*cos(deltam*ct - deltam/tau*sigma*sigma);
            two= +2.0*nwwerf(std::conj(z)).imag()/tau/4*exp(-ct*ct/2.0/sigma/sigma) +
                    exp(sigma*sigma/2 *(1/tau/tau - deltam*deltam) - ct/tau)*sin(deltam*ct - deltam/tau*sigma*sigma);
        }

    }

    if (std::isnan(one) || std::isnan(two)) {
        //void convoluted_exp_sincos(double ct, double sigma, double tau, double deltam, double& one, double& two)
        spdlog::error("ct= {0:f} sigma={1:f} tau={2:f} deltam={3:f} one={4:f} two={5:f}",
                ct, sigma, tau, deltam, one, two);
        spdlog::error("Error: NaN");
        assert(0);

    }
    return;
}


double daughtermom(double m, double m1, double m2){
    //E1+E2 = sqrt(m1*m1+p1*p1)+sqrt(m2*m2+p2*p2), where p1=-p2
    //=> p= sqrt(m*m-(m1+m2)*(m1+m2))*(m*m-(m1-m2)*(m1-m2))/(2.0*m)
    //E1 = (m*m+m1*m1-m2*m2)/(2.0*M)
    //E2 = (m*m-m1*m1+m2*m2)/(2.0*M)
    //So this returns daughter momentum and has nothing to do with the mother (=mom)

    double arg = (m*m-(m1+m2)*(m1+m2))*(m*m-(m1-m2)*(m1-m2));//original, reactivate

    if (arg < 0.0 || m1+m2>m)   return 0.0;

    return sqrt(arg)/(2.0*m);
    //all ok, confirmed that it gives breakup momentum
}


std::complex<double> mkpi_simple_bw_kstar_amp(double mkpi, double qsquared, double gammakstar, double mkstar){
    double gamma = gammakstar;
    std::complex<double> bwamp = 1.0/std::complex<double>(mkstar*mkstar-mkpi*mkpi,-mkstar*gamma);
    //combine everything, eq 28-30 in LHCb-ANA-2013-053
    return sqrt(mkstar*gamma/TMath::Pi()) * bwamp;
}

double mkpi_simple_bw_kstar_amp_squared(double mkpi, double qsquared, double gammakstar, double mkstar){
    std::complex<double> result = mkpi_simple_bw_kstar_amp(mkpi, qsquared, gammakstar, mkstar);
    return norm(result);
}

double gsl_mkpi_simple_bw_kstar_amp_squared(double x, void* par){
    double* v = (double*)par;
    return mkpi_simple_bw_kstar_amp_squared(sqrt(x), v[0], v[1], v[2]);
}

double int_mkpi_simple_bw_kstar_amp_squared(const double& mkpia, const double& mkpib, const double& qsquared, const double& gammakstar, const double& mkstar){
    double pars[3] = {qsquared, gammakstar, mkstar};
    ROOT::Math::GSLIntegrator gslint;
    gslint.SetFunction(gsl_mkpi_simple_bw_kstar_amp_squared, pars);
    return gslint.Integral(mkpia*mkpia, mkpib*mkpib);
}


std::complex<double> mkpi_simple_kstarzero_amp(double mkpi, double qsquared, double f800mag, double f800phase, double gammaf800, double mf800, double gammakstarzero, double mkstarzero){
    std::complex<double> gk(f800mag*cos(f800phase),f800mag*sin(f800phase));
    std::complex<double> denomkappa(mf800,-0.5*gammaf800);
    std::complex<double> denomkstar(mkstarzero,-0.5*gammakstarzero);

    std::complex<double> result = (-gk/(denomkappa*denomkappa-mkpi*mkpi)
                                   +1.0/(denomkstar*denomkstar-mkpi*mkpi));
    return result;
}

double mkpi_simple_kstarzero_amp_squared(double mkpi, double qsquared, double f800mag, double f800phase, double gammaf800, double mf800, double gammakstarzero, double mkstarzero)
{
    std::complex<double> result = mkpi_simple_kstarzero_amp(mkpi, qsquared, f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero);
    return norm(result);
}

double gsl_mkpi_simple_kstarzero_amp_squared(double x, void* par)
{
    double* v = (double*)par;
    return mkpi_simple_kstarzero_amp_squared(sqrt(x), v[0], v[1], v[2], v[3], v[4], v[5], v[6]);//this is for integrating over mkpi^2, only then is the norm 1
}

double gsl_mkpi_simple_re_bw_kstar_amp_kstarzero_amp_bar(double x, void* par)
{
    double* v = (double*)par;
    std::complex<double> result = mkpi_simple_bw_kstar_amp_kstarzero_amp_bar(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9]);
    return result.real();
}

double gsl_mkpi_simple_im_bw_kstar_amp_kstarzero_amp_bar(double x, void* par)
{
    double* v = (double*)par;
    std::complex<double> result = mkpi_simple_bw_kstar_amp_kstarzero_amp_bar(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9]);
    return result.imag();
}

double gsl_mkpi_simple_re_bw_kstar_amp_bar_kstarzero_amp(double x, void* par)
{
    double* v = (double*)par;
    std::complex<double> result = mkpi_simple_bw_kstar_amp_bar_kstarzero_amp(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9]);
    return result.real();
}

double gsl_mkpi_simple_im_bw_kstar_amp_bar_kstarzero_amp(double x, void* par)
{
    double* v = (double*)par;
    std::complex<double> result = mkpi_simple_bw_kstar_amp_bar_kstarzero_amp(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9]);
    return result.imag();
}



double int_mkpi_simple_kstarzero_amp_squared(const double& mkpia, const double& mkpib, const double& qsquared, const double& f800mag, const double& f800phase, const double& gammaf800, const double& mf800, const double& gammakstarzero, const double& mkstarzero){
    double pars[7] = {qsquared, f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero};
    ROOT::Math::GSLIntegrator gslint;//(1.0e-10, 1.0e-8, 10000);
    gslint.SetFunction(gsl_mkpi_simple_kstarzero_amp_squared, pars);
    return gslint.Integral(mkpia*mkpia, mkpib*mkpib);
}

std::complex<double> mkpi_simple_bw_kstar_amp_kstarzero_amp_bar(double mkpi, double qsquared,
                                                                double gammakstar, double mkstar, double asphase,
                                                                double f800mag, double f800phase, double gammaf800, double mf800,
                                                                double gammakstarzero, double mkstarzero){
    std::complex<double> A = mkpi_simple_bw_kstar_amp(mkpi, qsquared, gammakstar, mkstar);
    std::complex<double> B = std::complex<double>(cos(asphase),sin(asphase))
            *mkpi_simple_kstarzero_amp(mkpi, qsquared, f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero);
    return A*conj(B);
}

std::complex<double> mkpi_simple_bw_kstar_amp_bar_kstarzero_amp(double mkpi, double qsquared,
                                                                double gammakstar, double mkstar, double asphase,
                                                                double f800mag, double f800phase, double gammaf800, double mf800,
                                                                double gammakstarzero, double mkstarzero){
    std::complex<double> A = mkpi_simple_bw_kstar_amp(mkpi, qsquared, gammakstar, mkstar);
    std::complex<double> B = std::complex<double>(cos(asphase),sin(asphase))
            *mkpi_simple_kstarzero_amp(mkpi, qsquared, f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero);
    return conj(A)*B;
}

std::complex<double> int_mkpi_simple_bw_kstar_amp_kstarzero_amp_bar(const double& mkpia, const double& mkpib, const double& qsquared,
                                                                    const double& gammakstar, const double& mkstar, const double& asphase,
                                                                    const double& f800mag, const double& f800phase,
                                                                    const double& gammaf800, const double& mf800,
                                                                    const double& gammakstarzero, const double& mkstarzero){
    double pars[10] = {qsquared, gammakstar, mkstar, asphase,
                       f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero};
    ROOT::Math::GSLIntegrator gslintre;
    gslintre.SetFunction(gsl_mkpi_simple_re_bw_kstar_amp_kstarzero_amp_bar, pars);
    double result_re = gslintre.Integral(mkpia, mkpib);
    ROOT::Math::GSLIntegrator gslintim;
    gslintim.SetFunction(gsl_mkpi_simple_im_bw_kstar_amp_kstarzero_amp_bar, pars);
    double result_im = gslintim.Integral(mkpia, mkpib);
    return std::complex<double>(result_re, result_im);
}

std::complex<double> int_mkpi_simple_bw_kstar_amp_bar_kstarzero_amp(const double& mkpia, const double& mkpib, const double& qsquared,
                                                                    const double& gammakstar, const double& mkstar, const double& asphase,
                                                                    const double& f800mag, const double& f800phase,
                                                                    const double& gammaf800, const double& mf800,
                                                                    const double& gammakstarzero, const double& mkstarzero){
    double pars[10] = {qsquared, gammakstar, mkstar, asphase,
                       f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero};
    //double pars[8] = {qsquared, gammakstar, mkstar, asphase, a, r, gammakstarzero, mkstarzero};
    ROOT::Math::GSLIntegrator gslintre;
    gslintre.SetFunction(gsl_mkpi_simple_re_bw_kstar_amp_bar_kstarzero_amp, pars);
    double result_re = gslintre.Integral(mkpia, mkpib);
    ROOT::Math::GSLIntegrator gslintim;
    gslintim.SetFunction(gsl_mkpi_simple_im_bw_kstar_amp_bar_kstarzero_amp, pars);
    double result_im = gslintim.Integral(mkpia, mkpib);
    return std::complex<double>(result_re, result_im);
}

//the kstar amplitude
std::complex<double> mkpi_bw_kstar_amp(double mkpi, double qsquared, double gammakstar, double mkstar, double R){
    const double mk = PDGMASS_KST_KAON/1000.0;
    const double mpi = PDGMASS_KST_PION/1000.0;
    const double mb = PDGMASS_B/1000.0;
    double p     = daughtermom(mkpi, mk, mpi);
    double pzero = daughtermom(mkstar, mk, mpi);
    double k     = daughtermom(mb, mkpi, sqrt(qsquared));
    double phsp  = sqrt(k*p);
    //double gamma = gammakstar*pow(p/pzero, 3.0)*mkstar/mkpi*sqrt((1.0+R*R*pzero*pzero)/(1.0+R*R*p*p));//nominal
    double gamma = gammakstar*pow(p/pzero, 3.0)*mkstar/mkpi*((1.0+R*R*pzero*pzero)/(1.0+R*R*p*p));//squared??? !!!
    std::complex<double> bwamp = 1.0/std::complex<double>(mkstar*mkstar-mkpi*mkpi,-mkstar*gamma);
    //combine everything, eq 28-30 in LHCb-ANA-2013-053
    return phsp * bwamp
            * sqrt((1.0+R*R*pzero*pzero)/(1.0+R*R*p*p))
            * p/mkstar;
}

double mkpi_bw_kstar_amp_squared(double mkpi, double qsquared, double gammakstar, double mkstar, double R){
    std::complex<double> result = mkpi_bw_kstar_amp(mkpi, qsquared, gammakstar, mkstar, R);
    return norm(result);
}

double mkpi_kstarzero_lass_amp_squared(double mkpi, double qsquared, double asphase, double a, double r, double gammakstarzero, double mkstarzero, double R){
    std::complex<double> result = mkpi_kstarzero_lass_amp(mkpi, qsquared, asphase, a, r, gammakstarzero, mkstarzero, R);
    return norm(result);
}

//the combined s-wave amplitude
std::complex<double> mkpi_kstarzero_lass_amp(double mkpi, double qsquared, double asphase, double a, double r, double gammakstarzero, double mkstarzero, double R){
    const double mk = PDGMASS_KST_KAON/1000.0;
    const double mpi = PDGMASS_KST_PION/1000.0;
    const double mb = PDGMASS_B/1000.0;
    double p = daughtermom(mkpi, mk, mpi);
    double pzero = daughtermom(mkstarzero, mk, mpi);
    double k = daughtermom(mb, mkpi, sqrt(qsquared));
    double kzero = daughtermom(mb, mkstarzero, sqrt(qsquared));
    double phsp = sqrt(k*p);
    // const double R = 1.6;//[GeV^-1] from Belle, see Z paper LHCb-ANA-2013-053
    //as amplitude
    double asmag = 1.0;
    std::complex<double> as(asmag*cos(asphase),asmag*sin(asphase));
    //LASS shape
    double cotdeltaB = 1.0/(a*p) + 0.5*r*p;
    double deltaB = atan(1.0/cotdeltaB);
    double gamma = gammakstarzero*p/pzero*mkstarzero/mkpi;
    double cotdeltaR = (mkstarzero*mkstarzero - mkpi*mkpi)/(mkstarzero*gamma);
    std::complex<double> A = 1.0/std::complex<double>(cotdeltaB, -1.0);
    std::complex<double> B = 1.0/std::complex<double>(cotdeltaR, -1.0);
    std::complex<double> factor(cos(2.0*deltaB),sin(2.0*deltaB));
    std::complex<double> lass = A+factor*B;
    //combine everything, eq 52-53 in LHCb-ANA-2013-053
    return phsp * as * lass
            //* mkpi/p ///David: this is not derivable from given formulas, but somehow appears in eq 52 of LHCb-ANA-2013-053
            * k/mb
            * sqrt((1.0+R*R*kzero*kzero)/(1.0+R*R*k*k));
}


std::complex<double> mkpi_kstarzero_isobar_amp(double mkpi, double qsquared, double f800mag, double f800phase, double gammaf800, double mf800, double gammakstarzero, double mkstarzero, double R){
    const double mk = PDGMASS_KST_KAON/1000.0;
    const double mpi = PDGMASS_KST_PION/1000.0;
    const double mb = PDGMASS_B/1000.0;
    //const double R = 1.6;//[GeV^-1] from Belle, see Z paper LHCb-ANA-2013-053

    double p = daughtermom(mkpi, mk, mpi);
    double pzero = daughtermom(mf800, mk, mpi);
    double k = daughtermom(mb, mkpi, sqrt(qsquared));
    double kzero = daughtermom(mb, mf800, sqrt(qsquared));
    double phsp = sqrt(k*p);
    double gamma = gammaf800*pow(p/pzero, 1.0)*mf800/mkpi;
    std::complex<double> f800amp =
            phsp
            *k/mb*sqrt((1.0+R*R*kzero*kzero)/(1.0+R*R*k*k))
            *1.0/std::complex<double>(mf800*mf800-mkpi*mkpi,-mf800*gamma);

    p = daughtermom(mkpi, mk, mpi);
    pzero = daughtermom(mkstarzero, mk, mpi);
    k = daughtermom(mb, mkpi, sqrt(qsquared));
    kzero = daughtermom(mb, mkstarzero, sqrt(qsquared));
    phsp = sqrt(k*p);
    gamma = gammakstarzero*pow(p/pzero, 1.0)*mkstarzero/mkpi;
    std::complex<double> kstarzeroamp =
            phsp
            *k/mb*sqrt((1.0+R*R*kzero*kzero)/(1.0+R*R*k*k))
            *1.0/std::complex<double>(mkstarzero*mkstarzero-mkpi*mkpi,-mkstarzero*gamma);

    return f800mag*std::complex<double>(cos(f800phase),sin(f800phase))*f800amp + (1.0-f800mag)*kstarzeroamp;
}

double mkpi_kstarzero_isobar_amp_squared(double mkpi, double qsquared, double f800mag, double f800phase, double gammaf800, double mf800, double gammakstarzero, double mkstarzero, double R){
    std::complex<double> result = mkpi_kstarzero_isobar_amp(mkpi, qsquared, f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero, R);
    return norm(result);
}


std::complex<double> mkpi_bw_kstar_amp_kstarzero_isobar_amp_bar(double mkpi, double qsquared,
                                                                double gammakstar, double mkstar, double asphase,
                                                                double f800mag, double f800phase,
                                                                double gammaf800, double mf800,
                                                                double gammakstarzero, double mkstarzero, double R){
    std::complex<double> A = mkpi_bw_kstar_amp(mkpi, qsquared, gammakstar, mkstar, R);
    std::complex<double> B = std::complex<double>(cos(asphase),sin(asphase))
            *mkpi_kstarzero_isobar_amp(mkpi, qsquared, f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero, R);
    return A*conj(B);
}

std::complex<double> mkpi_bw_kstar_amp_bar_kstarzero_isobar_amp(double mkpi, double qsquared,
                                                                double gammakstar, double mkstar, double asphase,
                                                                double f800mag, double f800phase,
                                                                double gammaf800, double mf800,
                                                                double gammakstarzero, double mkstarzero, double R){
    std::complex<double> A = mkpi_bw_kstar_amp(mkpi, qsquared, gammakstar, mkstar, R);
    std::complex<double> B = std::complex<double>(cos(asphase),sin(asphase))
            *mkpi_kstarzero_isobar_amp(mkpi, qsquared, f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero, R);
    return conj(A)*B;
}

std::complex<double> mkpi_bw_kstar_amp_kstarzero_lass_amp_bar(double mkpi, double qsquared,
                                                              double gammakstar, double mkstar,
                                                              double asphase, double a, double r,
                                                              double gammakstarzero, double mkstarzero, double R){
    std::complex<double> A = mkpi_bw_kstar_amp(mkpi, qsquared, gammakstar, mkstar, R);
    std::complex<double> B = mkpi_kstarzero_lass_amp(mkpi, qsquared, asphase, a, r, gammakstarzero, mkstarzero, R);
    return A*conj(B);
}

std::complex<double> mkpi_bw_kstar_amp_bar_kstarzero_lass_amp(double mkpi, double qsquared,
                                                              double gammakstar, double mkstar,
                                                              double asphase, double a, double r,
                                                              double gammakstarzero, double mkstarzero, double R){
    std::complex<double> A = mkpi_bw_kstar_amp(mkpi, qsquared, gammakstar, mkstar, R);
    std::complex<double> B = mkpi_kstarzero_lass_amp(mkpi, qsquared, asphase, a, r, gammakstarzero, mkstarzero, R);
    return conj(A)*B;
}

double gsl_mkpi_bw_kstar_amp_squared(double x, void* par){
    double* v = (double*)par;
    return mkpi_bw_kstar_amp_squared(x, v[0], v[1], v[2], v[3]);
}

double gsl_mkpi_kstarzero_lass_amp_squared(double x, void* par){
    double* v = (double*)par;
    return mkpi_kstarzero_lass_amp_squared(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6]);
}

double gsl_mkpi_re_bw_kstar_amp_kstarzero_lass_amp_bar(double x, void* par){
    double* v = (double*)par;
    std::complex<double> result = mkpi_bw_kstar_amp_kstarzero_lass_amp_bar(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8]);
    return result.real();
}

double gsl_mkpi_im_bw_kstar_amp_kstarzero_lass_amp_bar(double x, void* par){
    double* v = (double*)par;
    std::complex<double> result = mkpi_bw_kstar_amp_kstarzero_lass_amp_bar(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8]);
    return result.imag();
}

double gsl_mkpi_re_bw_kstar_amp_bar_kstarzero_lass_amp(double x, void* par){
    double* v = (double*)par;
    std::complex<double> result = mkpi_bw_kstar_amp_bar_kstarzero_lass_amp(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8]);
    return result.real();
}

double gsl_mkpi_im_bw_kstar_amp_bar_kstarzero_lass_amp(double x, void* par){
    double* v = (double*)par;
    std::complex<double> result = mkpi_bw_kstar_amp_bar_kstarzero_lass_amp(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8]);
    return result.imag();
}

double gsl_mkpi_kstarzero_isobar_amp_squared(double x, void* par){
    double* v = (double*)par;
    return mkpi_kstarzero_isobar_amp_squared(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7]);
}

double gsl_mkpi_re_bw_kstar_amp_kstarzero_isobar_amp_bar(double x, void* par){
    double* v = (double*)par;
    std::complex<double> result = mkpi_bw_kstar_amp_kstarzero_isobar_amp_bar(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9], v[10]);
    return result.real();
}

double gsl_mkpi_im_bw_kstar_amp_kstarzero_isobar_amp_bar(double x, void* par){
    double* v = (double*)par;
    std::complex<double> result = mkpi_bw_kstar_amp_kstarzero_isobar_amp_bar(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9], v[10]);
    return result.imag();
}

double gsl_mkpi_re_bw_kstar_amp_bar_kstarzero_isobar_amp(double x, void* par){
    double* v = (double*)par;
    std::complex<double> result = mkpi_bw_kstar_amp_bar_kstarzero_isobar_amp(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9], v[10]);
    return result.real();
}

double gsl_mkpi_im_bw_kstar_amp_bar_kstarzero_isobar_amp(double x, void* par){
    double* v = (double*)par;
    std::complex<double> result = mkpi_bw_kstar_amp_bar_kstarzero_isobar_amp(x, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8], v[9], v[10]);
    return result.imag();
}

double int_mkpi_bw_kstar_amp_squared(double mkpia, double mkpib, double qsquared, double gammakstar,
                                     double mkstar, double R){
    double pars[4] = {qsquared, gammakstar, mkstar, R};
    ROOT::Math::GSLIntegrator gslint;
    gslint.SetFunction(gsl_mkpi_bw_kstar_amp_squared, pars);
    return gslint.Integral(mkpia, mkpib);
}

double int_mkpi_bw_kstar_amp_squared(double qsquared, double gammakstar, double mkstar, double R){
    double pars[4] = {qsquared, gammakstar, mkstar, R};
    ROOT::Math::GSLIntegrator gslint;//(1.0e-9, 1.0e-8, 100000);
    gslint.SetFunction(gsl_mkpi_bw_kstar_amp_squared, pars);
    return gslint.Integral();
}

double int_mkpi_kstarzero_lass_amp_squared(double mkpia, double mkpib, double qsquared, double asphase, double a, double r, double gammakstarzero, double mkstarzero, double R)
{
    double pars[7] = {qsquared, asphase, a, r, gammakstarzero, mkstarzero, R};
    ROOT::Math::GSLIntegrator gslint;
    gslint.SetFunction(gsl_mkpi_kstarzero_lass_amp_squared, pars);
    return gslint.Integral(mkpia, mkpib);
}

double int_mkpi_kstarzero_lass_amp_squared(double qsquared, double asphase, double a, double r, double gammakstarzero, double mkstarzero, double R)
{
    double pars[7] = {qsquared, asphase, a, r, gammakstarzero, mkstarzero, R};
    ROOT::Math::GSLIntegrator gslint;//(1.0e-9, 1.0e-8, 100000);
    gslint.SetFunction(gsl_mkpi_kstarzero_lass_amp_squared, pars);
    return gslint.Integral();
}

double int_mkpi_kstarzero_isobar_amp_squared(double mkpia, double mkpib, double qsquared, double f800mag, double f800phase, double gammaf800, double mf800, double gammakstarzero, double mkstarzero, double R)
{
    double pars[8] = {qsquared, f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero, R};
    ROOT::Math::GSLIntegrator gslint;
    gslint.SetFunction(gsl_mkpi_kstarzero_isobar_amp_squared, pars);
    return gslint.Integral(mkpia, mkpib);
}

double int_mkpi_kstarzero_isobar_amp_squared(double qsquared, double f800mag, double f800phase, double gammaf800, double mf800, double gammakstarzero, double mkstarzero, double R)
{
    double pars[8] = {qsquared, f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero, R};
    ROOT::Math::GSLIntegrator gslint;//(1.0e-9, 1.0e-8, 100000);
    gslint.SetFunction(gsl_mkpi_kstarzero_isobar_amp_squared, pars);
    return gslint.Integral();
}

std::complex<double> int_mkpi_bw_kstar_amp_kstarzero_lass_amp_bar(double mkpia, double mkpib, double qsquared, double gammakstar, double mkstar, double asphase, double a, double r, double gammakstarzero, double mkstarzero, double R)
{
    double pars[9] = {qsquared, gammakstar, mkstar, asphase, a, r, gammakstarzero, mkstarzero, R};
    ROOT::Math::GSLIntegrator gslintre;
    gslintre.SetFunction(gsl_mkpi_re_bw_kstar_amp_kstarzero_lass_amp_bar, pars);
    double result_re = gslintre.Integral(mkpia, mkpib);
    ROOT::Math::GSLIntegrator gslintim;
    gslintim.SetFunction(gsl_mkpi_im_bw_kstar_amp_kstarzero_lass_amp_bar, pars);
    double result_im = gslintim.Integral(mkpia, mkpib);
    return std::complex<double>(result_re, result_im);
}

std::complex<double> int_mkpi_bw_kstar_amp_bar_kstarzero_lass_amp(double mkpia, double mkpib, double qsquared, double gammakstar, double mkstar, double asphase, double a, double r, double gammakstarzero, double mkstarzero, double R)
{
    double pars[9] = {qsquared, gammakstar, mkstar, asphase, a, r, gammakstarzero, mkstarzero, R};
    ROOT::Math::GSLIntegrator gslintre;
    gslintre.SetFunction(gsl_mkpi_re_bw_kstar_amp_bar_kstarzero_lass_amp, pars);
    double result_re = gslintre.Integral(mkpia, mkpib);
    ROOT::Math::GSLIntegrator gslintim;
    gslintim.SetFunction(gsl_mkpi_im_bw_kstar_amp_bar_kstarzero_lass_amp, pars);
    double result_im = gslintim.Integral(mkpia, mkpib);
    return std::complex<double>(result_re, result_im);
}


std::complex<double> int_mkpi_bw_kstar_amp_kstarzero_lass_amp_bar(double qsquared, double gammakstar, double mkstar, double asphase, double a, double r, double gammakstarzero, double mkstarzero, double R)
{
    double pars[9] = {qsquared, gammakstar, mkstar, asphase, a, r, gammakstarzero, mkstarzero, R};
    ROOT::Math::GSLIntegrator gslintre;//(1.0e-9, 1.0e-8, 100000);
    gslintre.SetFunction(gsl_mkpi_re_bw_kstar_amp_kstarzero_lass_amp_bar, pars);
    double result_re = gslintre.Integral();//mkpisquareda, mkpisquaredb);
    ROOT::Math::GSLIntegrator gslintim;//(1.0e-9, 1.0e-8, 100000);
    gslintim.SetFunction(gsl_mkpi_im_bw_kstar_amp_kstarzero_lass_amp_bar, pars);
    double result_im = gslintim.Integral();//mkpisquareda, mkpisquaredb);
    return std::complex<double>(result_re, result_im);
}

std::complex<double> int_mkpi_bw_kstar_amp_bar_kstarzero_lass_amp(double qsquared, double gammakstar, double mkstar, double asphase, double a, double r, double gammakstarzero, double mkstarzero, double R)
{
    double pars[9] = {qsquared, gammakstar, mkstar, asphase, a, r, gammakstarzero, mkstarzero, R};
    ROOT::Math::GSLIntegrator gslintre;//(1.0e-9, 1.0e-8, 100000);
    gslintre.SetFunction(gsl_mkpi_re_bw_kstar_amp_bar_kstarzero_lass_amp, pars);
    double result_re = gslintre.Integral();//mkpisquareda, mkpisquaredb);
    ROOT::Math::GSLIntegrator gslintim;//(1.0e-9, 1.0e-8, 100000);
    gslintim.SetFunction(gsl_mkpi_im_bw_kstar_amp_bar_kstarzero_lass_amp, pars);
    double result_im = gslintim.Integral();//mkpisquareda, mkpisquaredb);
    return std::complex<double>(result_re, result_im);
}

///////////////
std::complex<double> int_mkpi_bw_kstar_amp_kstarzero_isobar_amp_bar(double mkpia, double mkpib, double qsquared,
                                                                    double gammakstar, double mkstar, double asphase,
                                                                    double f800mag, double f800phase, double gammaf800, double mf800, double gammakstarzero, double mkstarzero, double R)
{
    double pars[11] = {qsquared, gammakstar, mkstar, asphase,
                       f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero, R};
    ROOT::Math::GSLIntegrator gslintre;
    gslintre.SetFunction(gsl_mkpi_re_bw_kstar_amp_kstarzero_isobar_amp_bar, pars);
    double result_re = gslintre.Integral(mkpia, mkpib);
    ROOT::Math::GSLIntegrator gslintim;
    gslintim.SetFunction(gsl_mkpi_im_bw_kstar_amp_kstarzero_isobar_amp_bar, pars);
    double result_im = gslintim.Integral(mkpia, mkpib);
    return std::complex<double>(result_re, result_im);
}

std::complex<double> int_mkpi_bw_kstar_amp_bar_kstarzero_isobar_amp(double mkpia, double mkpib, double qsquared,
                                                                    double gammakstar, double mkstar, double asphase,
                                                                    double f800mag, double f800phase, double gammaf800, double mf800, double gammakstarzero, double mkstarzero, double R)
{
    double pars[11] = {qsquared, gammakstar, mkstar, asphase,
                       f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero, R};
    //double pars[8] = {qsquared, gammakstar, mkstar, asphase, a, r, gammakstarzero, mkstarzero};
    ROOT::Math::GSLIntegrator gslintre;
    gslintre.SetFunction(gsl_mkpi_re_bw_kstar_amp_bar_kstarzero_isobar_amp, pars);
    double result_re = gslintre.Integral(mkpia, mkpib);
    ROOT::Math::GSLIntegrator gslintim;
    gslintim.SetFunction(gsl_mkpi_im_bw_kstar_amp_bar_kstarzero_isobar_amp, pars);
    double result_im = gslintim.Integral(mkpia, mkpib);
    return std::complex<double>(result_re, result_im);
}


std::complex<double> int_mkpi_bw_kstar_amp_kstarzero_isobar_amp_bar(double qsquared,
                                                                    double gammakstar, double mkstar, double asphase,
                                                                    double f800mag, double f800phase, double gammaf800, double mf800, double gammakstarzero, double mkstarzero, double R)
{
    double pars[11] = {qsquared, gammakstar, mkstar, asphase,
                       f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero, R};
    //double pars[8] = {qsquared, gammakstar, mkstar, asphase, a, r, gammakstarzero, mkstarzero};
    ROOT::Math::GSLIntegrator gslintre;//(1.0e-9, 1.0e-8, 100000);
    gslintre.SetFunction(gsl_mkpi_re_bw_kstar_amp_kstarzero_isobar_amp_bar, pars);
    double result_re = gslintre.Integral();//mkpisquareda, mkpisquaredb);
    ROOT::Math::GSLIntegrator gslintim;//(1.0e-9, 1.0e-8, 100000);
    gslintim.SetFunction(gsl_mkpi_im_bw_kstar_amp_kstarzero_isobar_amp_bar, pars);
    double result_im = gslintim.Integral();//mkpisquareda, mkpisquaredb);
    return std::complex<double>(result_re, result_im);
}

std::complex<double> int_mkpi_bw_kstar_amp_bar_kstarzero_isobar_amp(double qsquared,
                                                                    double gammakstar, double mkstar, double asphase,
                                                                    double f800mag, double f800phase, double gammaf800, double mf800, double gammakstarzero, double mkstarzero, double R)
{
    double pars[11] = {qsquared, gammakstar, mkstar, asphase,
                       f800mag, f800phase, gammaf800, mf800, gammakstarzero, mkstarzero, R};
    ROOT::Math::GSLIntegrator gslintre;//(1.0e-9, 1.0e-8, 100000);
    gslintre.SetFunction(gsl_mkpi_re_bw_kstar_amp_bar_kstarzero_isobar_amp, pars);
    double result_re = gslintre.Integral();//mkpisquareda, mkpisquaredb);
    ROOT::Math::GSLIntegrator gslintim;//(1.0e-9, 1.0e-8, 100000);
    gslintim.SetFunction(gsl_mkpi_im_bw_kstar_amp_bar_kstarzero_isobar_amp, pars);
    double result_im = gslintim.Integral();//mkpisquareda, mkpisquaredb);
    return std::complex<double>(result_re, result_im);
}

void moments_to_observables(TRandom3* rnd,
                            std::vector<double> mom,
                            std::vector<double> momcov,
                            double fsext,
                            double dfsext,
                            std::vector<double>& obs,
                            std::vector<double>& obscov,
                            bool usefsext
                            ){

    assert(mom.size() == 18);
    assert(momcov.size() == 18*18);

    if(spdlog_debug()){
        std::cout << std::setw(10) << std::setprecision(6);
        for(UInt_t j = 0; j < 18; j++){
            for(UInt_t k = 0; k < 18; k++){
                std::cout << momcov.at(j*18+k);
            }
            std::cout << std::endl;
        }
    }

    //without s-waves, the matrix only comprises of 12 elements:
    bool swave = true;
    if(momcov.at(12) == 0.0 && momcov.at(13) == 0.0 && momcov.at(14) == 0.0 && momcov.at(15) == 0.0 && momcov.at(16) == 0.0 && momcov.at(17) == 0.0){
        swave = false;
    }
    spdlog::info("[MoM]\tSwave is "+std::string(swave ? "ON" : "OFF")+"!");
    const UInt_t N = swave ? 18 : 12;

    //convert 18x18 matrix to 12x12 matrix
    if(!swave){
        std::vector<double> momcov_tmp;
        for(UInt_t i = 0; i < 12; i++){
            for(UInt_t j = 0; j < 12; j++){
                momcov_tmp.push_back(momcov.at(i*18+j));
            }
        }
        momcov.clear();
        momcov = momcov_tmp;
    }

    assert(momcov.size() == N*N);

    //all off-diagonal matrix elements for S1c, S2s, S2c are forced to zero
    //all digaonal covariant elements are set equal to the one of S1s
    for (unsigned int i=0; i<N; i++){
        for (unsigned int j=0; j<N; j++){
            if (i==1 || i==2 || i==3 || j==1 || j==2 || j==3){
                if (i==j) momcov.at(i*N+j) = momcov.at(0);
                else momcov.at(i*N+j) = 0.0;
            }
        }
    }

    //cholesky decomposition of covariance matrix
    TMatrixDSym cov_sym(N, &momcov[0]);
    TDecompChol chol(cov_sym);
    bool success = chol.Decompose();
    assert(success);
    TMatrixD mU = chol.GetU();//gets upper triangular matrix
    //decomposition check
    
    TMatrixD mUT(mU);
    mUT.Transpose(mUT);
    TMatrixD check = mUT * mU;
    double distance = 0.0;
    for (unsigned int i=0; i<N; i++){
        for (unsigned int j=0; j<N; j++){
            distance += fabs(momcov.at(i+j*N) - check(i,j));
        }
    }
    spdlog::info("[MoM]\tDistance from Cholesky decomposition is {e:0}", distance);
    assert(fabs(distance) < 1.0e-8);

    //determine mean and rms (or rather confidence interval) from toys according to cov matrix
    unsigned int nruns = 10000;
    std::vector<std::vector< double > > randomizedobs;
    for (unsigned int l=0; l<nruns; l++){

        std::vector<double> rndobs(N, 0.0);

        //generate 10000 random moments according to covariance matrix
        //rndm gaussian(0,1) vector
        std::vector<double> rndvec(N, 0.0);
        for (unsigned int i=0;i<N;i++){
            //if (i!=1 && i!=2 && i!=3)
            rndvec.at(i) = rnd->Gaus(0.0, 1.0);
        }

        //random gaussian vector times upper triangular matrix
        std::vector<double> murndvec(N, 0.0);
        for (unsigned int i=0; i<N; i++){
            for (unsigned int j=0; j<N; j++){
                murndvec.at(i) += mU(j,i)*rndvec.at(j);
            }
        }

        //new random vector
        std::vector<double> result(N, 0.0);
        for (unsigned int i=0; i<N; i++){
            result.at(i) = mom.at(i) + murndvec.at(i);
        }

        //now get observables
        rndobs = std::vector<double>(N, 0.0);
        double j6c = 0.0;
        double pi = TMath::Pi();

        bool useLinearCombinationForS1s = false;
        bool useLinearCombinationForS6s = true;

        if(swave){

            double fl = (80.0*pi*result.at(0)+80.0*pi*result.at(12)-27.0)/(80.0*pi*result.at(0)+240.0*pi*result.at(12)-51.0);
            double fs = (80.0*pi*result.at(0)+240.0*pi*result.at(12)-45.0)/6.0;
            if (usefsext){
                fs = fsext+rnd->Gaus(0.0,dfsext);
            }
            rndobs.at(0) = 3.0/4.0*(1.0-fl); //S_1s
            rndobs.at(1) = 0.0; //S_1c
            rndobs.at(2) = 0.0; //S_2s
            rndobs.at(3) = 0.0; //S_2c
            rndobs.at(4) = -100.0*pi*result.at(4)/(9.0*fs-9.0); //S_3
            rndobs.at(5) = -100.0*pi*result.at(5)/(9.0*fs-9.0); //S_4
            rndobs.at(6) = -80*pi*result.at(6)/(9.0*fs-9.0); //S_5
            if(useLinearCombinationForS6s){
                rndobs.at(7) = -32.0*pi*(3.0*result.at(7) - 2.0*result.at(8))/(9.0*fs-9.0); //S_6s
                rndobs.at(8) = -32.0*pi*(8.0*result.at(8) - 2.0*result.at(7))/(9.0*fs-9.0); //S_6c
            }
            else{
                rndobs.at(7) = -(9.0*j6c*fs+320.0*pi*result.at(7)-9*j6c)/(36.0*fs-36.0); //S_6s
                double j6s = rndobs.at(7);
                rndobs.at(8) = -(18.0*j6s*fs+640.0*pi*result.at(8)-18.0*j6s)/(27.0*fs-27.0); //S_6c
                //no strange usage of S_6c:
                rndobs.at(7) = -80.0*pi*result.at(7)/(9.0*fs-9.0); //S_6s
                rndobs.at(8) = 0.0; //S_6c
            }
            rndobs.at(9) = -80.0*pi*result.at(9)/(9.0*fs-9.0);//S_7
            rndobs.at(10) = -100.0*pi*result.at(10)/(9.0*fs-9.0);//S_8
            rndobs.at(11) = -100.0*pi*result.at(11)/(9.0*fs-9.0);//S_9
            rndobs.at(12) = fs;//rnd->Gaus();
            rndobs.at(13) = 20.0*pi*result.at(13);//S_S1
            rndobs.at(14) = 20.0*pi*result.at(14);//S_S2
            rndobs.at(15) = 16.0*pi*result.at(15);//S_S3
            rndobs.at(16) = 16.0*pi*result.at(16);//S_S4
            rndobs.at(17) = 20.0*pi*result.at(17);//S_S5
        }
        else{
            if(useLinearCombinationForS1s){
                rndobs.at(0) = 32.0/9.0*pi*(21.0/16.0*result.at(0)-7.0/8.0*result.at(1)+1125.0/1184.0*result.at(2)-25.0/37.0*result.at(3)); //David recalculation
            }
            else{
                rndobs.at(0) = 15.0/8.0*(32.0/9.0*pi*result.at(0) - 2.0/5.0); //Thommis thesis
            }
            rndobs.at(1) = 0.0; //S_1c
            rndobs.at(2) = 0.0; //S_2s
            rndobs.at(3) = 0.0; //S_2c
            rndobs.at(4) = 100.0/9.0*pi*result.at(4); //S_3
            rndobs.at(5) = 100.0/9.0*pi*result.at(5); //S_4
            rndobs.at(6) = 80.0/9.0*pi*result.at(6); //S_5
            if(useLinearCombinationForS6s){
                rndobs.at(7) = 32.0/9.0*pi*(3.0*result.at(7) - 2.0*result.at(8)); //S_6s
                rndobs.at(8) = 32.0/9.0*pi*(8.0*result.at(8) - 2.0*result.at(7)); //S_6c
            }
            else{
                rndobs.at(7) = 80.0/9.0*pi*result.at(7);
                rndobs.at(8) = 0.0;
            }
            rndobs.at(9) = 80.0/9.0*pi*result.at(9); //S_7
            rndobs.at(10) = 100.0/9.0*pi*result.at(10); //S_8
            rndobs.at(11) = 100.0/9.0*pi*result.at(11); //S_9
        }
        randomizedobs.push_back(rndobs);
    }
    std::vector<double> res(18,0.0);
    std::vector<double> rescov(18*18,0.0);

    //determine mean value of moments M_i
    for (unsigned int i=0; i<nruns; i++){
        for (unsigned int j=0; j<N; j++){
            res.at(j) += randomizedobs.at(i).at(j)/nruns;
        }
    }

    //determine variance for each M_i
    for (unsigned int i=0; i<nruns; i++){
        for (unsigned int j=0; j<N; j++){
            for (unsigned int k=0; k<N; k++){
                rescov.at(j*18+k) += (randomizedobs.at(i).at(j) - res.at(j)) * (randomizedobs.at(i).at(k) - res.at(k)) / (nruns);
                if(std::isnan(rescov.at(j*18+k))){
                    spdlog::info("[MoM]\tCov element ({0:d},{1:d}) is nan at run={2:d}",j, k, i);
                    assert(0);
                }
            }
        }
    }

    obs = res;
    obscov = rescov;

    return;
}

///////////////

double threshold(double x, double thr, double n){
    return pow(x-thr,1.0/n);
}

double int_threshold(double a, double b, double thr, double n){
    return pow(b-thr,1.0/n+1.0)/(1.0/n+1.0)
            -pow(a-thr,1.0/n+1.0)/(1.0/n+1.0);
}

//int throw_multivariate_normal(unsigned int ndim, unsigned int ntoys, std::vector<double> xi, std::vector<double> cov, std::vector<std::vector<double> > results)
int throw_multivariate_normal(const unsigned int ndim, const unsigned int ntoys, const std::vector<double>& xi, const std::vector<double>& cov, std::vector<std::vector<double> >& results){
    assert(xi.size() == ndim);
    assert(cov.size() == ndim*ndim);
    results.clear();
    TRandom3* rnd = new TRandom3();
    
    TMatrixDSym cov_sym(ndim, &cov[0]);
    TDecompChol chol(cov_sym);
    bool success = chol.Decompose();

    assert(success);
    TMatrixD mU = chol.GetU();//gets upper triangular matrix
    
    for (unsigned int l=0;l<ntoys;l++){
        std::vector<double> result(ndim, 0.0);

        std::vector<double> rndvec(ndim, 0.0);
        for (unsigned int i=0;i<ndim;i++){
            rndvec.at(i) = rnd->Gaus(0.0, 1.0);
        }
        //random gaussian vector times upper triangular matrix
        std::vector<double> murndvec(ndim, 0.0);
        for (unsigned int i=0; i<ndim; i++){
            for (unsigned int j=0; j<ndim; j++){
                murndvec.at(i) += mU(j,i)*rndvec.at(j);
            }
        }

        //new random vector
        for (unsigned int i=0; i<ndim; i++){
            result.at(i) = xi.at(i) + murndvec.at(i);
        }
        results.push_back(result);
    }
    return 0;
}

int get_mean_cov(const std::vector<std::vector<double> >& meas, std::vector<double>& mean, std::vector<double>& cov){
    assert(meas.size() > 0);
    unsigned int ndim = meas.at(0).size();
    unsigned int ntoys = meas.size();

    std::vector<double> m(ndim, 0.0);
    
    for (unsigned int l=0; l<meas.size(); l++){
        for (unsigned int i=0; i<ndim; i++){
            m.at(i) += (meas.at(l).at(i))/(double(ntoys));
        }
    }

    std::vector<long double> c(ndim*ndim, 0.0);
    for (unsigned int l=0; l<meas.size(); l++){
        for (unsigned int i=0; i<ndim; i++){
            for (unsigned int j=0; j<ndim; j++){
                c.at(i*ndim+j) += (m.at(i)-meas.at(l).at(i))*(m.at(j)-meas.at(l).at(j))/double(ntoys);
            }
        }
    }
    mean.resize(m.size());
    cov.resize(c.size());
    for (unsigned int l=0; l<m.size(); l++){
        mean.at(l) = m.at(l);
    }
    for (unsigned int l=0; l<c.size(); l++){
        cov.at(l) = c.at(l);
    }

    return 0;
}

void sis_to_pis(TRandom3* rnd,
                std::vector<double> sis,
                std::vector<double> sicov,
                std::vector<double>& pis,
                std::vector<double>& picov){
    assert(sis.size()==8);
    assert(sicov.size()==8*8);
    
    //cholesky decomposition of covariance matrix
    TMatrixDSym cov_sym(8, &sicov[0]);
    TDecompChol chol(cov_sym);
    bool success = chol.Decompose();
    assert(success);
    TMatrixD mU = chol.GetU();//gets upper triangular matrix
    //decomposition check
    //determine mean and rms (or rather confidence interval) from toys according to cov matrix
    unsigned int nruns = 50000;
    std::vector<std::vector< double > > randomizedpis;
    for (unsigned int l=0; l<nruns; l++) {

        //generate random moments according to covariance matrix
        //rndm gaussian(0,1) vector
        std::vector<double> rndvec(8, 0.0);
        for (unsigned int i=0;i<8;i++){
            rndvec.at(i) = rnd->Gaus(0.0, 1.0);
        }
        //random gaussian vector times upper triangular matrix
        std::vector<double> murndvec(8, 0.0);
        for (unsigned int i=0; i<8; i++){
            for (unsigned int j=0; j<8; j++){
                murndvec.at(i) += mU(j,i)*rndvec.at(j);
            }
        }
        //new random vector
        std::vector<double> result(8, 0.0);
        for (unsigned int i=0; i<8; i++){
            result.at(i) = sis.at(i) + murndvec.at(i);
        }
        //now get observables
        //rndpi = std::vector<double>(8, 0.0);
        std::vector<double> rndpi(8, 0.0);

        double fl = result.at(0);
        double ft = 1.0-fl;
        double s3 = result.at(1);
        double s4 = result.at(2);
        double s5 = result.at(3);
        double afb = result.at(4);
        double s7 = result.at(5);
        double s8 = result.at(6);
        double s9 = result.at(7);

        rndpi.at(0) = fl;//fl
        rndpi.at(1) = 2.0*s3/(1.0-fl);//p1
        rndpi.at(2) = 0.5*4.0/3.0*afb/(1.0-fl);//p2
        rndpi.at(3) = -s9/(1-fl);//p3
        rndpi.at(4) = s4/sqrt(fl*ft);//p4
        rndpi.at(5) = s5/sqrt(fl*ft);//p5
        rndpi.at(6) = s7/sqrt(fl*ft);//p6
        rndpi.at(7) = s8/sqrt(fl*ft);//p8

        randomizedpis.push_back(rndpi);
    }
    std::vector<double> res(8,0.0);
    std::vector<double> rescov(8*8,0.0);
    for (unsigned int i=0; i<nruns; i++){
        for (unsigned int j=0; j<8; j++){
            res.at(j) += randomizedpis.at(i).at(j)/nruns;
        }
    }
    for (unsigned int i=0; i<nruns; i++){
        for (unsigned int j=0; j<8; j++){
            for (unsigned int k=0; k<8; k++){
                rescov.at(j*8+k) += (randomizedpis.at(i).at(j)-res.at(j))*(randomizedpis.at(i).at(k)-res.at(k))/(nruns);
            }
        }
    }
    pis = res;
    picov = rescov;
}

double rndgauss(TRandom3* rnd, double mean, double width, double min, double max){
    double r = 0.0;
    bool finished = false;
    while (!finished) {
        r = rnd->Gaus(mean, width);
        if (r > min && r < max) finished = true;
    }
    return r;
}

double rndgauss(TRandom3* rnd, double mean, double widthlo, double widthhi, double min, double max)
{
    double r = 0.0;
    bool finished = false;
    bool lo = rnd->Rndm() < 0.5;
    while (!finished){
        if (lo) r = mean - fabs(rnd->Gaus(0.0, widthlo));
        else        r = mean + fabs(rnd->Gaus(0.0, widthhi));

        if (r > min && r < max) finished = true;
    }
    return r;
}

std::vector<UInt_t> GetNOutOfM(UInt_t N, UInt_t M, Int_t Seed, bool DoubleCounts, bool PoissonFluctuate){

    TRandom3 * R = new TRandom3(Seed);

    if(PoissonFluctuate) N = TMath::Abs(R->Poisson(N));

    std::vector<UInt_t> v;

    if(DoubleCounts){
        while(v.size() < N){
            UInt_t I = R->Integer(M);
            v.push_back(I);
        }
    }
    else{
        assert(M > N);
        if(N <= M / 2){ //for small N get white list of Integers
            while(v.size() < N){
                UInt_t I = R->Integer(M);
                bool AlreadyChosen = false;
                for(UInt_t i = 0; i < v.size(); i++){
                    if(I == v.at(i)){
                        AlreadyChosen = true;
                        break;
                    }
                }
                if(!AlreadyChosen){
                    v.push_back(I);
                }
            }
        }
        else{ //large N: get black list and remove those numbers
            for(UInt_t i = 0; i < M; i++){
                v.push_back(i);
            }
            while(v.size() > N){
                UInt_t I = R->Integer(M);
                for(UInt_t i = 0; i < v.size(); i++){
                    if(I == v.at(i)){
                        v.erase(v.begin()+i);
                        break;
                    }
                }
            }
        }
    }

    assert(v.size() == N);

    return v;

}


}