CrossSectionAccumulator.hh

Go to the documentation of this file.

 
#pragma once
 
#include <iterator>
#include <map>
#include <ostream>
#include <vector>
 
#include "HEJ/event_types.hh"
 
namespace HEJ {
  class Event;
 
  template<typename T>
  struct XSWithError {
    T value = T{}; 
    T error = T{}; 
  };
 
  class CrossSectionAccumulator {
  public:
    void fill(Event const & ev);
    void fill(double wt, event_type::EventType type);
    void fill(double wt, double err, event_type::EventType type);
 
    void fill_correlated(std::vector<Event> const & evts);
    template<class ConstIt>
    void fill_correlated(ConstIt begin, ConstIt end);
    void fill_correlated(double sum, double sum2, event_type::EventType type);
 
    auto begin() const {
      return std::begin(xs_);
    }
    auto end() const {
      return std::end(xs_);
    }
    XSWithError<double> const & total() const {
      return total_;
    }
    XSWithError<double> const & at(event_type::EventType type) const {
      return xs_.at(type);
    }
    XSWithError<double> const & operator[](event_type::EventType type) const {
      return at(type);
    }
 
    // method to scale total_ and xs_ after they have been filled
    void scale(double xsscale);
 
  private:
    std::map<event_type::EventType, XSWithError<double>> xs_;
    XSWithError<double> total_;
  };
 
  std::ostream& operator<<(std::ostream& os, const CrossSectionAccumulator& xs);
 
// ------------ Implementation ------------
 
  template<class ConstIt>
  void CrossSectionAccumulator::fill_correlated(ConstIt begin, ConstIt end){
    if(std::distance(begin, end) < 2){ // only one event
      fill(*begin);
      return;
    }
    double sum = 0.;
    double sum2 = 0.;
    const auto type = begin->type();
    for(; begin != end; ++begin){
      double const wt = begin->central().weight;
      sum += wt;
      sum2 += wt*wt;
    }
    fill_correlated(sum, sum2, type);
  }
} // namespace HEJ