FitFraction.h

Go to the documentation of this file.

#ifndef AMPGEN_FITFRACTION_H
#define AMPGEN_FITFRACTION_H
 
#include <memory>
#include <string>
#include <vector>
#include <algorithm>
#include "AmpGen/Types.h"
#include "AmpGen/ErrorPropagator.h"
#include "AmpGen/AmplitudeRules.h"
 
namespace AmpGen {
  class EventType;
  class Particle;
 
  class FitFraction {
  public:
    FitFraction(const std::string &line);
    FitFraction(const std::string &name, const double &frac, const double &err);
    FitFraction() = default;
 
    void setFracErr(const double &f, const double &e);
    double val() const;
    double err() const;
    std::string name() const;
 
  private:
    std::string m_name;
    double m_value;
    double m_error;
  };
 
  bool operator<(const FitFraction &lhs, const FitFraction &rhs);
  bool operator>(const FitFraction &lhs, const FitFraction &rhs);
  bool operator==(const FitFraction &lhs, const FitFraction &rhs);
  std::ostream &operator<<(std::ostream &os, const FitFraction &obj);
 
  template <class pdf_type> struct FitFractionCalculator {
    struct fcalc {
      std::string name;
      std::vector<size_t> i;
      std::vector<size_t> j;
      fcalc(const std::string &name, const std::vector<size_t> &i) : name(name), i(i), j(i) {};
      fcalc(const std::string &name, const std::vector<size_t> &i, const std::vector<size_t> &j) : name(name), i(i), j(j) {};
    };
    pdf_type *pdf;
    std::vector<fcalc> calculators;
    std::vector<size_t> normSet;
    bool recalculateIntegrals;
    template <class... ARGS> void emplace_back(ARGS &&...args) { calculators.emplace_back(args...); }
 
    FitFractionCalculator(pdf_type *pdf, const std::vector<size_t> &normSet, const bool &recalculateIntegrals = false)
        : pdf(pdf), normSet(normSet), recalculateIntegrals(recalculateIntegrals) {}
    std::vector<double> operator()() {
      if(recalculateIntegrals)
        pdf->prepare();
      else
        pdf->transferParameters();
      std::vector<double> rv;
      double sum = 0;
      for(size_t i = 0; i != calculators.size(); ++i) {
        auto v = getVal(i);
        rv.push_back(v);
        sum += v;
      }
      rv.push_back(sum);
      return rv;
    }
    real_t norm() const {
      complex_t sum = 0;
      for(auto &i : normSet) {
        for(auto &j : normSet) { sum += (*pdf)[i].coefficient * std::conj((*pdf)[j].coefficient) * (j >= i ? pdf->norm(i, j) : std::conj(pdf->norm(j, i))); }
      }
      return std::real(sum);
    }
    real_t getVal(const size_t &index, const bool &getImaginaryPart = false) const {
      complex_t sum = 0;
      for(auto &i : calculators[index].i) {
        for(auto &j : calculators[index].j) {
          sum += (*pdf)[i].coefficient * std::conj((*pdf)[j].coefficient) * (j >= i ? pdf->norm(i, j) : std::conj(pdf->norm(j, i)));
        }
      }
      return (getImaginaryPart ? std::imag(sum) : std::real(sum)) / norm();
    }
    std::vector<FitFraction> operator()(const std::string &name, const LinearErrorPropagator &linProp) {
      auto values = (*this)();
      auto errors = linProp.getVectorError(*this, calculators.size() + 1);
      std::vector<FitFraction> fractions;
      for(size_t i = 0; i < calculators.size(); ++i) fractions.emplace_back(calculators[i].name, values[i], errors[i]);
      std::sort(fractions.begin(), fractions.end());
      std::reverse(fractions.begin(), fractions.end());
      fractions.emplace_back("Sum_" + name, *values.rbegin(), *errors.rbegin());
      return fractions;
    }
  };
} // namespace AmpGen
 
#endif